3 SUBROUTINE DECOMPRESS (IBUF, OBUF, NIN, NOUT)GC**********************************************************************DC_TITLE DECOMPRESS - Decompresses image lines stored in a compressedC format.CC_ARGS INTEGER*4 NINBC Input - number of bytes in>C the compressed buffer. INTEGER*4 NOUTEC Input - known number of bytesCC in decompressed output line CHARACTER IBUF(1)?C Input - line of data toEC decompress. The first byte isCC the actual pixel value, theDC rest of the array is the bit>C stream of coded "first;C difference" values. CHARACTER OBUF(1)EC Output - decompressed line ofAC data. The routine assumesDC original image composed of 8@C bit pixels. This is the=C restored image line -=C compression and firstCC differences undone - of the<C line passed in ibuf.CBC_DESC This routine decompresses Huffman encoded compressed data.BC (Huffman compression encoding can be found in any standard<C data compression reference. One such book is: "Data>C Compression, Techniques and Applications, Hardware andDC Software Considerations" by Gilbert Held, John Wiley & Sons,DC publishers.) Huffman coding uses variable number of bits toFC encode different values of the original data. The compressionAC results because the most frequently occurring values have3C the smaller number of bits for their codes.CAC The routine, used in conjuction with DECMPINIT, provides<C a common data base from which to call the processingFC subroutines. DECMPINIT builds the Huffman tree from the first-?C difference histogram and is called only once per image.CC DECOMPRESS processes one compressed input line per call and'C returns it completely restored.C?C_CALLS This routine calls two subroutines DCMPRS and HUFFTREE.CFC DCMPRS is provided in two versions, one FORTRAN, the other VAXCC MACRO. The FORTRAN version is for compatability on non-VAX:C systems so that the algorithm might be more easilyEC understood, while the MACRO version is for speed when running;C on a VAX. If you intend to use the macro version of;C the DCMPRS routine, then remove the fortran routineDC of the same name from the DECOMP source file before buildingC the system.CAC An alternate version of the decompression software can beCC found in the DECOMPS.FOR file. DECOMPS.FOR has been adapted1C for the SUN fortran compiler version 3.4.C>C_LIMS The routines have been tested using version 4.8-276 ofDC the VAX/VMS fortran compiler and versions V3.31 and V4.01 ofEC the Microsoft fortran compilers. The hardware used in testingFC the software include a VAX-750, an IBM PC/XT and a DEC RAINBOW@C PC. Modifications of the code may be required for other>C computer systems. See the DECOMPS.FOR file if you have5C software which will run on a SUN workstation.C@C_HIST DEC87, DMcMacken, ISD, USGS, Flagstaff, Original versionCC_ENDHC***********************************************************************" INTEGER*2 TREE(5,1024)4C Huffman tree INTEGER*2 VALUE INTEGER*2 RIGHT INTEGER*2 LEFT INTEGER*2 PARENT INTEGER*2 BRANCHD PARAMETER (VALUE=1, RIGHT=2, LEFT=3, PARENT=4, BRANCH=5) COMMON /DECOM/ TREEHC***********************************************************************C Decompress an input lineHC***********************************************************************. CALL DCMPRS(IBUF, OBUF, NIN, NOUT, TREE) RETURN END SUBROUTINE DECMPINIT(HIST)HC***********************************************************************DC_TITLE DECMPINIT - initialize the decompression TREE from the first(C difference histogramC_ARGS INTEGER*4 HIST(512)CC Input, difference histogramCC table. The histogram of the?C whole image after firstBC difference has been run on?C each line. In a "firstHC difference" line the first pixelEC retains its actual value; allDC other pixels are obtained byGC subtracting the actual value ofGC the pixel from the actual valueBC of the preceding pixel and?C adding 256 to provide aBC positive number. The i-thAC element of the array HISTFC should be the number of pixelsBC in the image with value i.C@C Some computer hardware systems may need to swap the byte orderDC of the 32-bit words which make up the first-difference historgram.FC The histogram is configured on the CDROM in "least significant byteDC first" order. This is the order for integer values used by VAX and@C IBM PC computer systems. Users of other computer architectures@C (IBM Mainframes, MacIntosh, SUN, and Apollo) will need to swapAC byte pairs 1 and 4, and 2 and 3 before passing the histogram toC the DECMPINIT subroutine.CC@C_DESC This routine initializes the Huffman tree using the first?C difference histogram passed by the calling program. ThisAC subroutine must be called before the DECOMPRESS subroutineC can be utilized.C?C_HIST DEC87, DMcMacken, ISD, USGS, Flagstaff, Original VersionCC_ENDHC***********************************************************************" INTEGER*2 TREE(5,1024) COMMON /DECOM/ TREE DO 10 I = 1,1024 DO 20 J = 1,5 TREE(J,I) = 0 20 CONTINUE 10 CONTINUE CALL HUFFTREE (HIST, TREE) RETURN END% SUBROUTINE HUFFTREE(HIST, TREE)GC**********************************************************************AC_TITLE HUFFTREE creates a Huffman code tree from input histogramCC_ARGS INTEGER*4 HIST(512)<C In - image histogram" INTEGER*2 TREE(5,1024)?C Out - Huffman code treeCEC_DESC This routine creates a Huffman code tree from the input firstHC difference histogram of an image. It starts by making all validHC (non-zero) densities for the image leaf nodes. These are sortedGC in increasing order by histogram frequency. The two nodes withDC smallest frequencies are connected by branches to a new node>C which is given a frequency equal to the sum of the twoFC combining nodes. The new node takes the place of the two nodesDC and the frequency list is resorted. The process is repeatedAC until the frequency set is reduced to one node. The treeDC consists of a double dimensioned array whose first dimensionFC gives five parameters for each node. The first parameter givesDC a value to the node. This is a density value for leaf nodesFC and a -1 for all others. The second parameter is a pointer toCC the next node on the right branch from this node. The thirdAC points to the node on the left branch. The fourth pointsDC to the parent node from which this node branches. The fifthAC parameter notes whether this node is on the right or left"C branch of the parent node.C>C_CALLS This routine calls the specially provided sort routineC SORTFREQ.CAC_HIST Fall86, DMcMacken, ISD, USGS, Flagstaff, Original versionCC_ENDGC**********************************************************************# INTEGER*4 FREQLIST(512)AC Histogramm frequency list# INTEGER*2 NODELIST(512)6C Tree node list INTEGER*2 VALUE INTEGER*2 RIGHT INTEGER*2 LEFT INTEGER*2 PARENT INTEGER*2 BRANCH@ PARAMETER (VALUE=1,RIGHT=2,LEFT=3,PARENT=4,BRANCH=5) INTEGER*4 NUMNODES4C Node counter INTEGER*4 I5C Do loop index INTEGER*4 PTR<C Current node pointer INTEGER*4 NUMFREQ>C # non-zero frequencies INTEGER*4 INDEX>C Frequency list pointerGC**********************************************************************'C Create leaf nodes, and pointer tablesGC********************************************************************** NUMNODES = 0 DO 10 I = 1, 512 FREQLIST(I) = HIST(I) NUMNODES = NUMNODES + 1 PTR = NUMNODES NODELIST(I) = PTR TREE(VALUE, PTR) = I TREE(PARENT, PTR) = 0 10 CONTINUEGC**********************************************************************BC Make sure the last element is always 0. For the first difference*C histogram, there can only be 511 values.GC********************************************************************** FREQLIST(512) = 0GC**********************************************************************?C Sort freqlist in increasing order; skip over zero frequenciesGC********************************************************************** NUMFREQ = 5121 CALL SORTFREQ (FREQLIST, NODELIST, NUMFREQ) INDEX = 1& 20 IF (FREQLIST(INDEX) .EQ. 0) THEN INDEX = INDEX + 1 GOTO 20 ENDIF# NUMFREQ = NUMFREQ - INDEX + 1 30 IF (NUMFREQ .GT. 1) THEN NUMNODES = NUMNODES + 1 PTR = NUMNODES+ TREE(RIGHT, PTR) = NODELIST(INDEX), TREE(LEFT, PTR) = NODELIST(INDEX+1)- TREE(PARENT, TREE(RIGHT, PTR)) = PTR/ TREE(BRANCH, TREE(RIGHT, PTR)) = RIGHT, TREE(PARENT, TREE(LEFT, PTR)) = PTR- TREE(BRANCH, TREE(LEFT, PTR)) = LEFT TREE(VALUE, PTR) = -1. FREQLIST(INDEX + 1) = FREQLIST(INDEX) 1 + FREQLIST(INDEX + 1)" NODELIST(INDEX + 1) = PTR FREQLIST(INDEX) = 0 INDEX = INDEX + 1 NUMFREQ = NUMFREQ - 19 CALL SORTFREQ (FREQLIST(INDEX), NODELIST(INDEX), 1 NUMFREQ) GOTO 30 ENDIF" TREE(VALUE, 1024) = NUMNODES RETURN END6 SUBROUTINE SORTFREQ(FREQLIST, NODELIST, NUMFREQ)GC**********************************************************************CC_TITLE SORTFREQ sorts frequency and node lists in increasing freq.C order.CC_ARGS# INTEGER*4 FREQLIST(512)>C input - frequency list# INTEGER*2 NODELIST(512)9C input - node list INTEGER*4 NUMFREQEC input - # values in freq listCGC_DESC This routine uses an insertion sort to reorder a frequency listEC in order of increasing frequency. The corresponding elementsBC of the node list are reordered to maintain correspondence.CAC_HIST Fall86, DMcMacken, ISD, USGS, Flagstaff, Original versionCC_ENDHC*********************************************************************** INTEGER*4 I5C Do loop index INTEGER*4 J=C List position pointer INTEGER*4 TEMP1AC Temporary storage - freq. INTEGER*4 TEMP2@C Temporary storage - nodeHC***********************************************************************<C Save current element - starting with second - in temporary=C storage. Compare with all elements in first part of list -8C moving each up one element - until element is larger..C Insert current element at this point in listHC*********************************************************************** DO 20 I = 2, NUMFREQ TEMP1 = FREQLIST(I) TEMP2 = NODELIST(I) J = I+ 10 IF (FREQLIST(J-1) .GT. TEMP1) THEN' FREQLIST(J) = FREQLIST(J-1)' NODELIST(J) = NODELIST(J-1) J = J - 1! IF (J .GT. 1) GOTO 10 ENDIF FREQLIST(J) = TEMP1 NODELIST(J) = TEMP2 20 CONTINUE RETURN END4 SUBROUTINE DCMPRS(IBUF, OBUF, NIN, NOUT, TREE)HC***********************************************************************@C_TITLE DCMPRS decompresses Huffman coded compressed image lines?C (Remove this subroutine from the source code if you are?C planning to use the VAX/VMS DCMPRS.MAR macro routine inDC place of this fortran code. The DCMPRS macro version is more2C than twice as fast as the fortran version)CC_ARGS INTEGER*4 NINFC Input, # bytes in input buffer INTEGER*4 NOUTEC Input, # bytes in output line" INTEGER*2 TREE(5,1024)@C Input, Huffman code tree CHARACTER IBUF(1)EC Input, compressed data buffer CHARACTER OBUF(1)AC Output, decompressed lineCHC_DESC This subroutine follows a path from the root of the Huffman treeFC to one of its leaves. The choice at each branch is decided byGC the successive bits of the compressed input stream, left for 1,FC right for 0. Only leaf nodes have a value other than -1. TheDC routine traces a path through the tree until it finds a nodeEC with a value not equal to -1 (a leaf node). The value at theGC leaf node is subtracted from the preceding pixel value plus 256FC to restore the the umcompressed pixel. This algorithm is then:C repeated until the entire line has been processed.C7C_CALLS This routine uses the VAX/VMS implicit functionDC BTEST (VAR, IBIT) to test whether the bit number IBIT is set>C in the variable VAR. A fortran version of this routineC is available for BTEST.C@C_HIST DEC86, DMcMacken, ISD, USGS, Flagstaff, Original versionCC_ENDGC********************************************************************** INTEGER*4 PTRCC Pointer to position in tree INTEGER*4 K5C Do loop index INTEGER*4 BITFC Pointer to current bit in word INTEGER*4 NBYTEEC Counter, # bytes decompressed INTEGER*4 TEMP6C Working buffer CHARACTER TEMPB;C Byte working buffer INTEGER*2 VALUE;C Value index in tree INTEGER*2 RIGHTBC Right branch index in tree INTEGER*2 LEFTAC Left branch index in tree INTEGER*2 PARENTDC Parent pointer index in tree INTEGER*2 BRANCHCC Parent branch index in treeD PARAMETER (VALUE=1, RIGHT=2, LEFT=3, PARENT=4, BRANCH=5) INTEGER*2 DN=C Current density value INTEGER*2 DNL:C Last density value CHARACTER DNBDC Byte addressed density value LOGICAL END8C End of line flag LOGICAL BTESTC. EQUIVALENCE (DN, DNB), (TEMP, TEMPB)GC**********************************************************************C Start at root of treeGC********************************************************************** PTR = TREE(1,1024)GC***********************************************************************C Just copy first byte it is uncompressed.GC********************************************************************** OBUF(1) = IBUF(1)GC**********************************************************************C Count starts hereGC********************************************************************** NBYTE = 1GC**********************************************************************,C Initialize current and last density valuesGC********************************************************************** DN = 0 DNB = IBUF(1) DNL = DNGC**********************************************************************C Reset end of line flagGC********************************************************************** END = .FALSE.GC**********************************************************************>C K points to current input byte loop is done when k exceeds #C input bytes.GC********************************************************************** K = 2" IF (K .GT. NIN) END = .TRUE.GC**********************************************************************C This is the processing loopGC********************************************************************** 10 IF (.NOT. END) THEN TEMPB = IBUF(K) BIT = 7- 20 IF (BIT .GE. 0 .AND. .NOT. END) THEN& IF (BTEST(TEMP, BIT)) THEN$ PTR = TREE(LEFT, PTR) ELSE% PTR = TREE(RIGHT, PTR) ENDIFGC**********************************************************************)C We are at end leaf when value is not -1GC**********************************************************************. IF (TREE(VALUE, PTR) .NE. -1) THENGC**********************************************************************BC Compute value of pixel, reset or increment pointers and countersGC**********************************************************************$ DN = TREE(VALUE, PTR)" DN = DNL - DN + 256 DNL = DN NBYTE = NBYTE + 1 OBUF(NBYTE) = DNB" PTR = TREE(1, 1024)0 IF (NBYTE .EQ. NOUT) END = .TRUE. ENDIFGC**********************************************************************C Process next bit in byteGC********************************************************************** BIT = BIT - 1 GOTO 20 ENDIFGC**********************************************************************C Set index to next input byteGC********************************************************************** K = K + 1% IF (K .GT. NIN) END = .TRUE. GOTO 10 ENDIFGC***********************************************************************C Go back to caller with decompressed lineGC********************************************************************** RETURN END