f20-4

896 Chapter20. Less-NumericalAlgorithmsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).20.4 Huffman Codingand Compression of Data A lossless data compression algorithm takes a string of symbols (typically ASCII charactersorbytes)andtranslatesit reversibly intoanotherstring,onethatis ontheaverage ofshorterlength. Thewords“ontheaverage”arecrucial;itisobvious that no reversible algorithm can make all strings shorter — there just aren’t enough short strings to be in one-to-one correspondence with longer strings. Compression algorithms are possible only when, on the input side, some strings, or some inputsymbols, are more common than others. These can then be encoded in fewer bits than rarer input strings or symbols, giving a net average gain. There exist many, quite different, compression techniques, corresponding to differentwaysofdetectingandusingdeparturesfromequiprobabilityininputstrings. Inthissectionandthenextweshallconsideronly variablelengthcodes withdefined wordinputs. In these, the input is sliced into fixed units, for example ASCII characters, while the corresponding output comes in chunks of variable size. The simplest such method is Huffman coding [1], discussed in this section. Another example, arithmetic compression , is discussed in §20.5. At the opposite extremefrom defined-word,variable lengthcodes are schemes thatdivideupthe inputintounitsofvariablelength(wordsorphrasesofEnglishtext, forexample)andthentransmitthese,oftenwithafixed-lengthoutputcode. Themost widely used code of this type is the Ziv-Lempel code [2]. References [3-6]give the flavorofsome othercompressiontechniques,with referencesto thelargeliterature. The idea behind Huffman coding is simply to use shorter bit patterns for more commoncharacters. We can make this idea quantitative by consideringthe conceptofentropy. Suppose the input alphabet has N chcharacters, and that these occur in the input string with respective probabilities pi,i=1 ,...,N ch, so that/summationtextpi=1. Then the fundamental theorem of informationtheory says that strings consisting ofindependentlyrandomsequencesofthese characters(a conservative,but not always realistic assumption) require, on the average, at least H=−/summationdisplay p ilog2pi (20.4.1 ) bits per character. Here His the entropy of the probability distribution. Moreover, coding schemes exist which approach the bound arbitrarily closely. For the case of equiprobable characters, with all pi=1 /N ch, one easily sees that H=l o g2Nch, which is the case of no compression at all. Any other set of pi’s gives a smaller entropy, allowing some useful compression. Noticethattheboundof(20.4.1)wouldbeachievedifwecouldencodecharacter iwith a code of length Li=−log2pibits: Equation (20.4.1) would then be the average/summationtextpiLi. The trouble with such a scheme is that −log2piis not generally an integer. How can we encode the letter “Q” in 5.32 bits? Huffman coding makesa stab at this by, in effect, approximating all the probabilities p iby integer powers of 1/2, so that all the Li’s are integral. If all the pi’s are in fact of this form, then a Huffman code does achieve the entropy bound H. The construction of a Huffman code is best illustrated by example. Imagine a language, Vowellish, with the Nch=5character alphabet A, E, I, O, and U, occurringwiththerespectiveprobabilities0.12,0.42,0.09,0.30,and0.07. Thenthe constructionofaHuffmancodeforVowellishisaccomplishedinthefollowingtable: 20.4HuffmanCodingandCompressionofData 897Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).NodeStage: 1234 5 1A: 0.12 0.12 2E: 0.42 0.42 0.42 0.42 3I: 0.09 4O: 0.30 0.30 0.30 5U: 0.07 6 UI: 0.16 7 AUI: 0.28 8 AUIO: 0.58 9 EAUIO: 1.00 Here is how it works, proceedingin sequence through N chstages, represented by the columns of the table. The first stage starts with Nchnodes, one for each letterofthealphabet,containingtheirrespectiverelativefrequencies. At eachstage, the two smallest probabilities are found, summed to make a new node, and thendroppedfrom the list of active nodes. (A “block” denotes the stage where a node is dropped.) All active nodes (including the new composite) are then carried over to the next stage (column). In the table, the names assigned to new nodes (e.g., AUI)are inconsequential. In the example shown, it happens that (after stage 1) the two smallest nodes are always an original node and a composite one; this need not be trueingeneral: Thetwo smallest probabilitiesmightbebothoriginalnodes,orboth composites,oroneofeach. At thelast stage,all nodeswill havebeencollectedinto one grand composite of total probability 1. Now, to see the code, you redraw the data in the above table as a tree (Figure 20.4.1). As shown, each node of the tree corresponds to a node (row) in the table, indicated by the integer to its left and probabilityvalue to its right. Terminal nodes,so called, are shown as circles; these are single alphabetic characters. The branches ofthetreearelabeled0and1. Thecodeforacharacteris thesequenceofzerosand onesthatleadtoit,fromthetopdown. Forexample,Eis simply0,whileUis 1010. Any string of zeros and ones can now be decoded into an alphabetic sequence. Consider, for example, the string 1011111010. Starting at the top of the tree wedescend through 1011 to I, the first character. Since we have reached a terminal node, we reset to the top of the tree, next descending through11 to O. Finally 1010 gives U. The string thus decodes to IOU. These ideas are embodied in the following routines. Input to the first routine hufmakis an integer vector of the frequency of occurrence of the nchin ≡N ch alphabetic characters, i.e., a set of integers proportional to the pi’s.hufmak, along withhufapp,whichitcalls,performstheconstructionoftheabovetable,andalsothe tree of Figure 20.4.1. The routine utilizes a heap structure (see §8.3) for efficiency; for a detailed description, see Sedgewick [7]. 898 Chapter20. Less-NumericalAlgorithmsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).EEAUIO A UAUIAUIO UI IO1.00 0.58 0.28 0.30 0.09 0.07 3 50.16 60.120.4229 8 7 4 1 1 01 01 01 0 Figure 20.4.1. Huffman code for the fictitious language Vowellish, in tree form. A letter (A, E, I, O, or U) is encoded or decoded by traversing the tree from the top down; the code is the sequence of 0’s and 1’s on the branches. The value to the right of each node is its probability; to the left, its node number in the accompanying table. SUBROUTINE hufmak(nfreq,nchin,ilong,nlong) INTEGER ilong,nchin,nlong,nfreq(nchin),MC,MQPARAMETER (MC=512,MQ=2*MC-1) C USES hufapp Given the frequency of occurrence table nfreq(1:nchin) ofnchin characters, construct in the common block /hufcom/ the Huffman code. Returned values ilong andnlong are the character number that produced the longest code symbol, and the length of that symbol. You should check that nlong is not larger than your machine’s word length. INTEGER ibit,j,k,n,nch,node,nodemx,nused,ibset,index(MQ), * iup(MQ),icod(MQ),left(MQ),iright(MQ),ncod(MQ),nprob(MQ) COMMON /hufcom/ icod,ncod,nprob,left,iright,nch,nodemx SAVE /hufcom/ nch=nchin Initialization. nused=0do 11j=1,nch nprob(j)=nfreq(j) icod(j)=0ncod(j)=0 if(nfreq(j).ne.0)then nused=nused+1index(nused)=j endif enddo 11 do12j=nused,1,-1 Sort nprob into a heap structure in index . call hufapp(index,nprob,nused,j) enddo 12 k=nch 1 if(nused.gt.1)then Combine heap nodes, remaking the heap at each stage. node=index(1) index(1)=index(nused) nused=nused-1call hufapp(index,nprob,nused,1) k=k+1 20.4HuffmanCodingandCompressionofData 899Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).nprob(k)=nprob(index(1))+nprob(node) left(k)=node Store left and right children of a node. iright(k)=index(1) iup(index(1)) = -k Indicate whether a node is a left or right child of its parent. iup(node)=k index(1)=k call hufapp(index,nprob,nused,1) goto 1endif nodemx=k iup(nodemx)=0do 13j=1,nch Make the Huffman code from the tree. if(nprob(j).ne.0)then n=0 ibit=0node=iup(j) 2 if(node.ne.0)then if(node.lt.0)then n=ibset(n,ibit)node = -node endif node=iup(node)ibit=ibit+1 goto 2 endificod(j)=nncod(j)=ibit endif enddo 13 nlong=0 do14j=1,nch if(ncod(j).gt.nlong)then nlong=ncod(j)ilong=j-1 endif enddo 14 returnEND SUBROUTINE hufapp(index,nprob,m,l) INTEGER m,l,MC,MQPARAMETER (MC=512,MQ=2*MC-1) INTEGER index(MQ),nprob(MQ) Used by hufmak to maintain a heap structure in the array index(1:l) . INTEGER i,j,k,n n=m i=lk=index(i) 2 if(i.le.n/2)then j=i+i if (j.lt.n.and.nprob(index(j)).gt.nprob(index(j+1))) j=j+1if (nprob(k).le.nprob(index(j))) goto 3 index(i)=index(j) i=j goto 2endif 3 index(i)=k returnEND 900 Chapter20. Less-NumericalAlgorithmsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Once the code is constructed, one encodes a string of characters by repeated calls to hufenc, which simply does a table lookup of the code and appends it to the output message. SUBROUTINE hufenc(ich,code,lcode,nb) INTEGER ich,lcode,nb,MC,MQPARAMETER (MC=512,MQ=2*MC-1) Huffman encode the single character ich(in the range 0..nch-1 ), write the result to the character array code(1:lcode) starting at bit nb(whose smallest valid value is zero), and increment nbappropriately. This routine is called repeatedly to encode consecutive characters in a message, but must be preceded by a single initializing call to hufmak . INTEGER k,l,n,nc,nch,nodemx,ntmp,ibset INTEGER icod(MQ),left(MQ),iright(MQ),ncod(MQ),nprob(MQ)LOGICAL btestCHARACTER*1 code(*) COMMON /hufcom/ icod,ncod,nprob,left,iright,nch,nodemx SAVE /hufcom/k=ich+1 Convert character range 0..nch-1 to array index range 1..nch . if(k.gt.nch.or.k.lt.1)pause ’ich out of range in hufenc.’ do 11n=ncod(k),1,-1 Loop over the bits in the stored Huffman code for ich. nc=nb/8+1if (nc.gt.lcode) pause ’lcode too small in hufenc.’ l=mod(nb,8) if (l.eq.0) code(nc)=char(0)if(btest(icod(k),n-1))then Set appropriate bits in code. ntmp=ibset(ichar(code(nc)),l) code(nc)=char(ntmp) endifnb=nb+1 enddo 11 return END Decoding a Huffman-encoded message is slightly more complicated. The codingtreemustbetraversedfromthetopdown,usingupavariablenumberofbits: SUBROUTINE hufdec(ich,code,lcode,nb) INTEGER ich,lcode,nb,MC,MQPARAMETER (MC=512,MQ=2*MC-1) Starting at bit number nbin the character array code(1:lcode) , use the Huffman code stored in common block /hufcom/ to decode a single character (returned as ichin the range 0..nch-1 ) and increment nbappropriately. Repeated calls, starting with nb =0 will return successive characters in a compressed message. The returned value ich=nch indicates end-of-message. This routine must be preceded by a single initializing call to hufmak . Parameters: MCis the maximum value of nch, the input alphabet size. INTEGER l,nc,nch,node,nodemx INTEGER icod(MQ),left(MQ),iright(MQ),ncod(MQ),nprob(MQ)LOGICAL btestCHARACTER*1 code(lcode) COMMON /hufcom/ icod,ncod,nprob,left,iright,nch,nodemx SAVE /hufcom/node=nodemx Setnode to the top of the decoding tree. 1 continue Loop until a valid character is obtained. nc=nb/8+1if (nc.gt.lcode)then Ran out of input; with ich=nch indicating end of message. ich=nch return endifl=mod(nb,8) Now decoding this bit. nb=nb+1 20.4HuffmanCodingandCompressionofData 901Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).if(btest(ichar(code(nc)),l))then Branch left or right in tree, depending on its value. node=iright(node) else node=left(node) endif if(node.le.nch)then If we reach a terminal node, we have a complete character and can return. ich=node-1 return endif goto 1 END For simplicity, hufdecquits when it runs out of code bytes; if your coded message is not an integralnumberof bytes,and if N chis less than 256, hufdeccan return a spurious final character or two, decoded from the spurious trailing bits in yourlast codebyte. Ifyouhaveindependentknowledgeofthenumberofcharacters sent,youcanreadilydiscardthese. Otherwise,youcanfixthisbehaviorbyprovidinga bit, not byte, count, and modifyingthe routine accordingly. (When N chis 256 or larger, hufdecwill normally run out of code in the middle of a spurious character, and it will be discarded.) Run-Length Encoding For the compression of highly correlated bit-streams (for examplethe black or white values along a facsimile scan line), Huffman compression is often combined withrun-lengthencoding : Insteadofsendingeachbit,theinputstreamisconvertedto aseriesofintegersindicatinghowmanyconsecutivebitshavethesamevalue. Theseintegers are then Huffman-compressed. The Group 3 CCITT facsimile standard functions in this manner, with a fixed, immutable, Huffman code, optimized for a set of eight standard documents [8,9]. CITED REFERENCES AND FURTHER READING: Gallager,R.G. 1968, Information Theory and ReliableCommunication (New York: Wiley). Hamming, R.W. 1980, Codingand Information Theory (Englewood Cliffs, NJ: Prentice-Hall). Storer, J.A. 1988, Data Compression: Methods and Theory (Rockville, MD: Computer Science Press). Nelson, M. 1991, The Data Compression Book (Redwood City, CA: M&T Books). Huffman,D.A.1952, ProceedingsoftheInstituteofRadioEngineers ,vol.40,pp.1098–1101.[1] Ziv,J., andLempel, A. 1978, IEEE Transactionson Information Theory ,vol. IT-24,pp. 530–536. [2] Cleary, J.G., and Witten, I.H. 1984, IEEE Transactions on Communications , vol. COM-32, pp. 396–402. [3] Welch, T.A. 1984, Computer , vol. 17, no. 6, pp. 8–19. [4] Bentley, J.L., Sleator, D.D., Tarjan, R.E., and Wei, V.K. 1986, Communications of the ACM , vol. 29, pp. 320–330. [5] Jones, D.W. 1988, Communications of the ACM , vol. 31, pp. 996–1007. [6] Sedgewick, R. 1988, Algorithms , 2nd ed. (Reading, MA: Addison-Wesley), Chapter 22. [7] Hunter, R., and Robinson, A.H. 1980, Proceedings of the IEEE , vol. 68, pp. 854–867. [8] Marking, M.P. 1990, The C Users’ Journal , vol. 8, no. 6, pp. 45–54. [9] 902 Chapter20. Less-NumericalAlgorithmsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).20.5 Arithmetic Coding Wesawintheprevioussectionthataperfect(entropy-bounded)codingscheme would use Li=−log2pibits to encode character i(in the range 1≤i≤Nch), ifpiis its probability of occurrence. Huffman coding gives a way of rounding the Li’s to close integer values and constructing a code with those lengths. Arithmetic coding[1], which we now discuss, actually does manage to encode characters using noninteger numbers of bits! It also provides a convenient way to output the resultnot as a stream of bits, but as a stream of symbols in any desired radix. This latter property is particularly useful if you want, e.g., to convert data from bytes (radix 256) to printable ASCII characters (radix 94), or to case-independentalphanumericsequences containing only A-Z and 0-9 (radix 36). In arithmetic coding, an input message of any length is represented as a real number Rin the range 0≤R< 1. The longer the message, the more precision requiredof R. This is bestillustratedbyanexample,so letus returnto thefictitious language,Vowellish,oftheprevioussection. RecallthatVowellishhasa5characteralphabet (A, E, I, O, U), with occurrence probabilities 0.12, 0.42, 0.09, 0.30, and 0.07,respectively. Figure20.5.1showshowamessagebeginning“IOU”isencoded: The interval [0,1)is divided into segments corresponding to the 5 alphabetical characters; the length of a segment is the corresponding p i. We see that the first messagecharacter,“I”,narrowstherangeof Rto0.37≤R< 0.46. This intervalis nowsubdividedintofivesubintervals,againwithlengthsproportionaltothe pi’s. The second message character, “O”, narrows the range of Rto0.3763≤R< 0.4033. The“U” characterfurthernarrowsthe rangeto 0.37630 ≤R< 0.37819.Anyvalue ofRin this range can be sent as encoding “IOU”. In particular, the binary fraction .011000001 is in this range, so “IOU” can be sent in 9 bits. (Huffman coding took 10 bits for this example, see §20.4.) Ofcoursethereistheproblemofknowingwhentostopdecoding. Thefraction .011000001 representsnotsimply“IOU,”but“IOU ...,”wheretheellipsesrepresent an infinite string of successor characters. To resolve this ambiguity, arithmetic coding generally assumes the existence of a special Nch+1th character, EOM (end of message), which occurs only once at the end of the input. Since EOMhas a low probability of occurrence, it gets allocated only a very tiny piece of the number line. In the above example, we gave Ras a binary fraction. We could just as well have output it in any other radix, e.g., base 94 or base 36, whatever is convenient for the anticipated storage or communication channel. You might wonder how one deals with the seemingly incredible precision requiredof Rforalongmessage. Theansweristhat Risneveractuallyrepresented all at once. At any give stage we have upper and lower bounds for Rrepresented as a finite number of digits in the output radix. As digits of the upper and lower boundsbecomeidentical,we can left-shift them away andbringin new digits at the low-significance end. The routines below have a parameter NWKfor the number of working digits to keep around. This must be large enough to make the chance of an accidental degeneracy vanishingly small. (The routines signal if a degeneracy ever occurs.) Since the process of discarding old digits and bringingin new ones is performedidentically onencodingand decoding,everythingstays synchronized.