SUBROUTINE DECOMPRESS (IBUF, OBUF, NIN, NOUT)
C**********************************************************************
C_TITLE DECOMPRESS - Decompresses image lines stored in a compressed
C                    format. (alternate fortran souce code which will
C                    run on a SUN workstation)
C
C_ARGS
      INTEGER*4       NIN
C                                       Input - number of bytes in
C                                       the compressed buffer.
      INTEGER*4       NOUT
C                                       Input - known number of bytes
C                                       in decompressed output line
      CHARACTER       IBUF(1)
C                                       Input - line of data to
C                                       decompress. The first byte is
C                                       the actual pixel value, the
C                                       rest of the array is the bit
C                                       stream of coded "first
C                                       difference" values.
      CHARACTER       OBUF(1)
C                                       Output - decompressed line of
C                                       data. The routine assumes
C                                       original image composed of 8
C                                       bit pixels.  This is the
C                                       restored image line -
C                                       compression and first
C                                       differences undone - of the
C                                       line passed in ibuf.
C
C_DESC  This routine decompresses Huffman encoded compressed data.
C       (Huffman compression encoding can be found in any standard
C       data compression reference.  One such book is: "Data
C       Compression, Techniques and Applications, Hardware and
C       Software Considerations" by Gilbert Held, John Wiley & Sons,
C       publishers.)  Huffman coding uses variable number of bits to
C       encode different values of the original data.  The compression
C       results because the most frequently occurring values have
C       the smaller number of bits for their codes.
C
C       The routine, used in conjuction with DECMPINIT,  provides
C       a common data base from which to call the processing
C       subroutines. DECMPINIT builds the Huffman tree from the first-
C       difference histogram and is called only once per image.
C       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.
C
C
C_LIMS  The fortran code found in the DECOMPS.FOR file was adapted
C       to run on a SUN work station using fortran version 3.4
C
C_HIST  DEC87, DMcMacken, ISD, USGS, Flagstaff, Original version
C
C_END
C***********************************************************************
      INTEGER*2       TREE(5,1024)
C                                       Huffman tree
      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)
      COMMON /DECOM/ TREE
C***********************************************************************
C Decompress an input line
C***********************************************************************
      CALL DCMPRS(IBUF, OBUF, NIN, NOUT, TREE)
      RETURN
      END
      SUBROUTINE DECMPINIT(HIST)
C***********************************************************************
C_TITLE DECMPINIT - initialize the decompression TREE from the first
C                   difference histogram
C_ARGS
      INTEGER*4       HIST(512)
C                                       Input, difference histogram
C                                       table. The histogram of the
C                                       whole image after first
C                                       difference has been run on
C                                       each line.  In a "first
C                                       difference" line the first pixel
C                                       retains its actual value; all
C                                       other pixels are obtained by
C                                       subtracting the actual value of
C                                       the pixel from the actual value
C                                       of the preceding pixel and
C                                       adding 256 to provide a
C                                       positive number.  The i-th
C                                       element of the array HIST
C                                       should be the number of pixels
C                                       in the image with value i.
C
C Some computer hardware systems may need to swap the byte order
C of the 32-bit words which make up the first-difference historgram.
C The histogram  is configured on the CDROM in "least significant byte
C 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 swap
C byte pairs 1 and 4, and 2 and 3 before passing the histogram to
C the DECMPINIT subroutine.
C
C
C_DESC This routine initializes the Huffman tree using the first
C      difference histogram passed by the calling program. This
C      subroutine must be called before the DECOMPRESS subroutine
C      can be utilized.
C
C_HIST DEC87, DMcMacken, ISD, USGS, Flagstaff, Original Version
C
C_END
C***********************************************************************
      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)
C**********************************************************************
C_TITLE HUFFTREE creates a Huffman code tree from input histogram
C
C_ARGS
      INTEGER*4       HIST(512)
C                                       In - image histogram
      INTEGER*2       TREE(5,1024)
C                                       Out - Huffman code tree
C
C_DESC  This routine creates a Huffman code tree from the input first
C       difference histogram of an image.  It starts by making all valid
C       (non-zero) densities for the image leaf nodes.  These are sorted
C       in increasing order by histogram frequency.  The two nodes with
C       smallest frequencies are connected by branches to a new node
C       which is given a frequency equal to the sum of the two
C       combining nodes. The new node takes the place of the two nodes
C       and the frequency list is resorted.  The process is repeated
C       until the frequency set is reduced to one node.  The tree
C       consists of a double dimensioned array whose first dimension
C       gives five parameters for each node. The first parameter gives
C       a value to the node.  This is a density value for leaf nodes
C       and a -1 for all others.  The second parameter is a pointer to
C       the next node on the right branch from this node. The third
C       points to the node on the left branch.  The fourth points
C       to the parent node from which this node branches.  The fifth
C       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 routine
C               SORTFREQ.
C
C_HIST  Fall86, DMcMacken, ISD, USGS, Flagstaff, Original version
C
C_END
C**********************************************************************
      INTEGER*4       FREQLIST(512)
C                                       Histogramm frequency list
      INTEGER*2       NODELIST(512)
C                                       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       NUMNODES
C                                       Node counter
      INTEGER*4       I
C                                       Do loop index
      INTEGER*4       PTR
C                                       Current node pointer
      INTEGER*4       NUMFREQ
C                                       # non-zero frequencies
      INTEGER*4       INDEX
C                                       Frequency list pointer
C**********************************************************************
C Create leaf nodes, and pointer tables
C**********************************************************************
      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 CONTINUE
C**********************************************************************
C Make sure the last element is always 0. For the first difference
C histogram, there can only be 511 values.
C**********************************************************************
      FREQLIST(512) = 0
C**********************************************************************
C Sort freqlist in increasing order; skip over zero frequencies
C**********************************************************************
      NUMFREQ = 512
      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 - 1
         CALL SORTFREQ (FREQLIST(INDEX), NODELIST(INDEX),
     1   NUMFREQ)
         GOTO 30
      ENDIF
      TREE(VALUE, 1024) = NUMNODES
      RETURN
      END
      SUBROUTINE SORTFREQ(FREQLIST, NODELIST, NUMFREQ)
C**********************************************************************
C_TITLE SORTFREQ sorts frequency and node lists in increasing freq.
C       order.
C
C_ARGS
      INTEGER*4       FREQLIST(512)
C                                       input - frequency list
      INTEGER*2       NODELIST(512)
C                                       input - node list
      INTEGER*4       NUMFREQ
C                                       input - # values in freq list
C
C_DESC  This routine uses an insertion sort to reorder a frequency list
C       in order of increasing frequency.  The corresponding elements
C       of the node list are reordered to maintain correspondence.
C
C_HIST  Fall86, DMcMacken, ISD, USGS, Flagstaff, Original version
C
C_END
C***********************************************************************
      INTEGER*4       I
C                                       Do loop index
      INTEGER*4       J
C                                       List position pointer
      INTEGER*4       TEMP1
C                                       Temporary storage - freq.
      INTEGER*4       TEMP2
C                                       Temporary storage - node
C***********************************************************************
C Save current element - starting with second - in temporary
C storage.  Compare with all elements in first part of list -
C moving each up one element - until  element is larger.
C Insert current element at this point in list
C***********************************************************************
      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
      END
      SUBROUTINE DCMPRS(IBUF, OBUF, NIN, NOUT, TREE)
C***********************************************************************
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 in
C       place of this fortran code. The DCMPRS macro version is more
C       than twice as fast as the fortran version)
C
C_ARGS
      INTEGER*4       NIN
C                                       Input, # bytes in input buffer
      INTEGER*4       NOUT
C                                       Input, # bytes in output line
      INTEGER*2       TREE(5,1024)
C                                       Input, Huffman code tree
      CHARACTER       IBUF(1)
C                                       Input, compressed data buffer
      CHARACTER       OBUF(1)
C                                       Output, decompressed line
C
C_DESC  This subroutine follows a path from the root of the Huffman tree
C       to one of its leaves.  The choice at each branch is decided by
C       the successive bits of the compressed input stream, left for 1,
C       right for 0.  Only leaf nodes have a value other than -1.  The
C       routine traces a path through the tree until it finds a node
C       with a value not equal to -1 (a leaf node).  The value at the
C       leaf node is subtracted from the preceding pixel value plus 256
C       to restore the the umcompressed pixel.  This algorithm is then
C       repeated until the entire line has been processed.
C
C_CALLS This routine uses the VAX/VMS implicit function
C       BTEST (VAR, IBIT) to test whether the bit number IBIT is set
C       in the variable VAR. A fortran version of this routine
C       is available for BTEST.
C
C_HIST  DEC86, DMcMacken, ISD, USGS, Flagstaff, Original version
C
C_END
C**********************************************************************
      INTEGER*4       PTR
C                                       Pointer to position in tree
      INTEGER*4       K
C                                       Do loop index
      INTEGER*4       BIT
C                                       Pointer to current bit in word
      INTEGER*4       NBYTE
C                                       Counter, # bytes decompressed
      INTEGER*4       TEMP
C                                       Working variable
      INTEGER*2       VALUE
C                                       Value index in tree
      INTEGER*2       RIGHT
C                                       Right branch index in tree
      INTEGER*2       LEFT
C                                       Left branch index in tree
      INTEGER*2       PARENT
C                                       Parent pointer index in tree
      INTEGER*2       BRANCH
C                                       Parent branch index in tree
      PARAMETER       (VALUE=1, RIGHT=2, LEFT=3, PARENT=4, BRANCH=5)
C
      INTEGER*2       DN
C                                       Current density value
      INTEGER*2       DNL
C                                       Last density value
      LOGICAL       END
C                                       End of line flag
      LOGICAL       BTEST
C
C**********************************************************************
C Start at root of tree
C**********************************************************************
      PTR = TREE(1,1024)
C**********************************************************************
C Just copy first byte it is uncompressed.
C**********************************************************************
      OBUF(1) = IBUF(1)
C**********************************************************************
C Count starts here
C**********************************************************************
      NBYTE = 1
C**********************************************************************
C Initialize current and last density values
C**********************************************************************
      DN = ICHAR(IBUF(1))
      IF (DN.LT.0) DN = 256 + DN
      DNL = DN
C**********************************************************************
C Reset end of line flag
C**********************************************************************
      END = .FALSE.
C**********************************************************************
C K points to current input byte loop is done when k exceeds #
C input bytes.
C**********************************************************************
      K = 2
      IF (K .GT. NIN) END = .TRUE.
C**********************************************************************
C This is the processing loop
C**********************************************************************
   10 IF (.NOT. END) THEN
         TEMP = ICHAR(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)
            ENDIF
C**********************************************************************
C We are at end leaf when value is not -1
C**********************************************************************
            IF (TREE(VALUE, PTR) .NE. -1) THEN
C**********************************************************************
C Compute value of pixel, reset or increment pointers and counters
C**********************************************************************
               DN = TREE(VALUE, PTR)
               DN = DNL - DN + 256
               DNL = DN
               NBYTE = NBYTE + 1
               OBUF(NBYTE) = CHAR(DN)
               PTR = TREE(1, 1024)
               IF (NBYTE .EQ. NOUT) END = .TRUE.
            ENDIF
C**********************************************************************
C Process next bit in byte
C**********************************************************************
            BIT = BIT - 1
            GOTO 20
         ENDIF
C**********************************************************************
C Set index to next input byte
C**********************************************************************
         K = K + 1
         IF (K .GT. NIN) END = .TRUE.
         GOTO 10
      ENDIF
C**********************************************************************
C Go back to caller with decompressed line
C**********************************************************************
      RETURN
      END