home *** CD-ROM | disk | FTP | other *** search
/ Voyagers to the Outer Planets 2: Uranus / VoyagestotheOuterPlanetsVol2.cdr / software / decomps.for < prev    next >
Text File  |  1988-09-14  |  20KB  |  427 lines
  1.       SUBROUTINE DECOMPRESS (IBUF, OBUF, NIN, NOUT)
  2. C**********************************************************************
  3. C_TITLE DECOMPRESS - Decompresses image lines stored in a compressed 
  4. C                    format. (alternate fortran souce code which will
  5. C                    run on a SUN workstation)
  6. C
  7. C_ARGS
  8.       INTEGER*4       NIN             
  9. C                                       Input - number of bytes in
  10. C                                       the compressed buffer.
  11.       INTEGER*4       NOUT            
  12. C                                       Input - known number of bytes
  13. C                                       in decompressed output line
  14.       CHARACTER       IBUF(1)         
  15. C                                       Input - line of data to 
  16. C                                       decompress. The first byte is 
  17. C                                       the actual pixel value, the 
  18. C                                       rest of the array is the bit 
  19. C                                       stream of coded "first
  20. C                                       difference" values.
  21.       CHARACTER       OBUF(1)         
  22. C                                       Output - decompressed line of 
  23. C                                       data. The routine assumes 
  24. C                                       original image composed of 8 
  25. C                                       bit pixels.  This is the 
  26. C                                       restored image line - 
  27. C                                       compression and first 
  28. C                                       differences undone - of the 
  29. C                                       line passed in ibuf.
  30. C
  31. C_DESC  This routine decompresses Huffman encoded compressed data.  
  32. C       (Huffman compression encoding can be found in any standard 
  33. C       data compression reference.  One such book is: "Data 
  34. C       Compression, Techniques and Applications, Hardware and 
  35. C       Software Considerations" by Gilbert Held, John Wiley & Sons, 
  36. C       publishers.)  Huffman coding uses variable number of bits to 
  37. C       encode different values of the original data.  The compression 
  38. C       results because the most frequently occurring values have
  39. C       the smaller number of bits for their codes.
  40. C
  41. C       The routine, used in conjuction with DECMPINIT,  provides 
  42. C       a common data base from which to call the processing
  43. C       subroutines. DECMPINIT builds the Huffman tree from the first-
  44. C       difference histogram and is called only once per image.  
  45. C       DECOMPRESS processes one compressed input line per call and 
  46. C       returns it completely restored.
  47. C
  48. C_CALLS This routine calls two subroutines DCMPRS and HUFFTREE.
  49. C
  50. C
  51. C_LIMS  The fortran code found in the DECOMPS.FOR file was adapted
  52. C       to run on a SUN work station using fortran version 3.4
  53. C
  54. C_HIST  DEC87, DMcMacken, ISD, USGS, Flagstaff, Original version
  55. C
  56. C_END
  57. C***********************************************************************
  58.       INTEGER*2       TREE(5,1024)    
  59. C                                       Huffman tree
  60.       INTEGER*2       VALUE
  61.       INTEGER*2       RIGHT
  62.       INTEGER*2       LEFT
  63.       INTEGER*2       PARENT
  64.       INTEGER*2       BRANCH
  65.       PARAMETER       (VALUE=1, RIGHT=2, LEFT=3, PARENT=4, BRANCH=5)
  66.       COMMON /DECOM/ TREE
  67. C***********************************************************************
  68. C Decompress an input line
  69. C***********************************************************************
  70.       CALL DCMPRS(IBUF, OBUF, NIN, NOUT, TREE)
  71.       RETURN
  72.       END
  73.       SUBROUTINE DECMPINIT(HIST)
  74. C***********************************************************************
  75. C_TITLE DECMPINIT - initialize the decompression TREE from the first
  76. C                   difference histogram
  77. C_ARGS
  78.       INTEGER*4       HIST(512)       
  79. C                                       Input, difference histogram 
  80. C                                       table. The histogram of the 
  81. C                                       whole image after first 
  82. C                                       difference has been run on 
  83. C                                       each line.  In a "first 
  84. C                                       difference" line the first pixel
  85. C                                       retains its actual value; all 
  86. C                                       other pixels are obtained by 
  87. C                                       subtracting the actual value of 
  88. C                                       the pixel from the actual value 
  89. C                                       of the preceding pixel and 
  90. C                                       adding 256 to provide a 
  91. C                                       positive number.  The i-th 
  92. C                                       element of the array HIST 
  93. C                                       should be the number of pixels 
  94. C                                       in the image with value i.
  95. C
  96. C Some computer hardware systems may need to swap the byte order
  97. C of the 32-bit words which make up the first-difference historgram.  
  98. C The histogram  is configured on the CDROM in "least significant byte 
  99. C first" order. This is the order for integer values used by VAX and 
  100. C IBM PC computer systems. Users of other computer architectures
  101. C (IBM Mainframes, MacIntosh, SUN, and Apollo) will need to swap
  102. C byte pairs 1 and 4, and 2 and 3 before passing the histogram to
  103. C the DECMPINIT subroutine.
  104. C                                     
  105. C
  106. C_DESC This routine initializes the Huffman tree using the first
  107. C      difference histogram passed by the calling program. This 
  108. C      subroutine must be called before the DECOMPRESS subroutine
  109. C      can be utilized.
  110. C
  111. C_HIST DEC87, DMcMacken, ISD, USGS, Flagstaff, Original Version
  112. C
  113. C_END
  114. C***********************************************************************
  115.       INTEGER*2       TREE(5,1024)    
  116.       COMMON /DECOM/ TREE
  117.  
  118.       DO 10 I = 1,1024
  119.       DO 20 J = 1,5
  120.       TREE(J,I) = 0
  121.    20 CONTINUE
  122.    10 CONTINUE
  123.       CALL HUFFTREE (HIST, TREE)
  124.       RETURN
  125.       END
  126.       SUBROUTINE HUFFTREE(HIST, TREE)
  127. C**********************************************************************
  128. C_TITLE HUFFTREE creates a Huffman code tree from input histogram
  129. C
  130. C_ARGS
  131.       INTEGER*4       HIST(512)               
  132. C                                       In - image histogram
  133.       INTEGER*2       TREE(5,1024)            
  134. C                                       Out - Huffman code tree
  135. C
  136. C_DESC  This routine creates a Huffman code tree from the input first
  137. C       difference histogram of an image.  It starts by making all valid
  138. C       (non-zero) densities for the image leaf nodes.  These are sorted
  139. C       in increasing order by histogram frequency.  The two nodes with
  140. C       smallest frequencies are connected by branches to a new node 
  141. C       which is given a frequency equal to the sum of the two 
  142. C       combining nodes. The new node takes the place of the two nodes 
  143. C       and the frequency list is resorted.  The process is repeated 
  144. C       until the frequency set is reduced to one node.  The tree 
  145. C       consists of a double dimensioned array whose first dimension 
  146. C       gives five parameters for each node. The first parameter gives 
  147. C       a value to the node.  This is a density value for leaf nodes 
  148. C       and a -1 for all others.  The second parameter is a pointer to 
  149. C       the next node on the right branch from this node. The third 
  150. C       points to the node on the left branch.  The fourth points
  151. C       to the parent node from which this node branches.  The fifth
  152. C       parameter notes whether this node is on the right or left 
  153. C       branch of the parent node.
  154. C
  155. C_CALLS This routine calls the specially provided sort routine
  156. C               SORTFREQ.
  157. C
  158. C_HIST  Fall86, DMcMacken, ISD, USGS, Flagstaff, Original version
  159. C
  160. C_END
  161. C**********************************************************************
  162.       INTEGER*4       FREQLIST(512)          
  163. C                                       Histogramm frequency list
  164.       INTEGER*2       NODELIST(512)          
  165. C                                       Tree node list
  166.       INTEGER*2       VALUE
  167.       INTEGER*2       RIGHT
  168.       INTEGER*2       LEFT
  169.       INTEGER*2       PARENT
  170.       INTEGER*2       BRANCH
  171.       PARAMETER       (VALUE=1,RIGHT=2,LEFT=3,PARENT=4,BRANCH=5)
  172.       INTEGER*4       NUMNODES               
  173. C                                       Node counter
  174.       INTEGER*4       I                       
  175. C                                       Do loop index
  176.       INTEGER*4       PTR                     
  177. C                                       Current node pointer
  178.       INTEGER*4       NUMFREQ                
  179. C                                       # non-zero frequencies
  180.       INTEGER*4       INDEX                   
  181. C                                       Frequency list pointer
  182. C**********************************************************************
  183. C Create leaf nodes, and pointer tables
  184. C**********************************************************************
  185.       NUMNODES = 0
  186.       DO 10 I = 1, 512
  187.          FREQLIST(I) = HIST(I)
  188.          NUMNODES = NUMNODES + 1
  189.          PTR = NUMNODES
  190.          NODELIST(I) = PTR
  191.          TREE(VALUE, PTR) = I
  192.          TREE(PARENT, PTR) = 0
  193.    10 CONTINUE
  194. C**********************************************************************
  195. C Make sure the last element is always 0. For the first difference
  196. C histogram, there can only be 511 values.
  197. C**********************************************************************
  198.       FREQLIST(512) = 0
  199. C**********************************************************************
  200. C Sort freqlist in increasing order; skip over zero frequencies
  201. C**********************************************************************
  202.       NUMFREQ = 512
  203.       CALL SORTFREQ (FREQLIST, NODELIST, NUMFREQ)
  204.       INDEX = 1
  205.    20 IF (FREQLIST(INDEX) .EQ. 0) THEN
  206.          INDEX = INDEX + 1
  207.          GOTO 20
  208.       ENDIF
  209.  
  210.       NUMFREQ = NUMFREQ - INDEX + 1
  211.    30 IF (NUMFREQ .GT. 1) THEN
  212.          NUMNODES = NUMNODES + 1
  213.          PTR = NUMNODES
  214.          TREE(RIGHT, PTR) = NODELIST(INDEX)
  215.          TREE(LEFT, PTR) = NODELIST(INDEX+1)
  216.          TREE(PARENT, TREE(RIGHT, PTR)) = PTR
  217.          TREE(BRANCH, TREE(RIGHT, PTR)) = RIGHT
  218.          TREE(PARENT, TREE(LEFT, PTR)) = PTR
  219.          TREE(BRANCH, TREE(LEFT, PTR)) = LEFT
  220.          TREE(VALUE, PTR) = -1
  221.  
  222.          FREQLIST(INDEX + 1) = FREQLIST(INDEX)
  223.      1   + FREQLIST(INDEX + 1)
  224.          NODELIST(INDEX + 1) = PTR
  225.          FREQLIST(INDEX) = 0
  226.          INDEX = INDEX + 1
  227.          NUMFREQ = NUMFREQ - 1
  228.          CALL SORTFREQ (FREQLIST(INDEX), NODELIST(INDEX),
  229.      1   NUMFREQ)
  230.          GOTO 30
  231.       ENDIF
  232.       TREE(VALUE, 1024) = NUMNODES
  233.       RETURN
  234.       END
  235.       SUBROUTINE SORTFREQ(FREQLIST, NODELIST, NUMFREQ)
  236. C**********************************************************************
  237. C_TITLE SORTFREQ sorts frequency and node lists in increasing freq. 
  238. C       order.
  239. C
  240. C_ARGS
  241.       INTEGER*4       FREQLIST(512)          
  242. C                                       input - frequency list
  243.       INTEGER*2       NODELIST(512)          
  244. C                                       input - node list
  245.       INTEGER*4       NUMFREQ                
  246. C                                       input - # values in freq list
  247. C
  248. C_DESC  This routine uses an insertion sort to reorder a frequency list
  249. C       in order of increasing frequency.  The corresponding elements
  250. C       of the node list are reordered to maintain correspondence.
  251. C
  252. C_HIST  Fall86, DMcMacken, ISD, USGS, Flagstaff, Original version
  253. C
  254. C_END
  255. C***********************************************************************
  256.       INTEGER*4       I                       
  257. C                                       Do loop index
  258.       INTEGER*4       J                       
  259. C                                       List position pointer
  260.       INTEGER*4       TEMP1                   
  261. C                                       Temporary storage - freq.
  262.       INTEGER*4       TEMP2                   
  263. C                                       Temporary storage - node
  264. C***********************************************************************
  265. C Save current element - starting with second - in temporary
  266. C storage.  Compare with all elements in first part of list -
  267. C moving each up one element - until  element is larger.
  268. C Insert current element at this point in list
  269. C***********************************************************************
  270.       DO 20 I = 2, NUMFREQ
  271.          TEMP1 = FREQLIST(I)
  272.          TEMP2 = NODELIST(I)
  273.          J = I
  274.    10    IF (FREQLIST(J-1) .GT. TEMP1) THEN
  275.             FREQLIST(J) = FREQLIST(J-1)
  276.             NODELIST(J) = NODELIST(J-1)
  277.             J = J - 1
  278.             IF (J .GT. 1) GOTO 10
  279.          ENDIF
  280.          FREQLIST(J) = TEMP1
  281.          NODELIST(J) = TEMP2
  282.    20 CONTINUE
  283.       RETURN
  284.       END
  285.       SUBROUTINE DCMPRS(IBUF, OBUF, NIN, NOUT, TREE)
  286. C***********************************************************************
  287. C_TITLE DCMPRS decompresses Huffman coded compressed image lines
  288. C       (Remove this subroutine from the source code if you are
  289. C       planning to use the VAX/VMS DCMPRS.MAR macro routine in
  290. C       place of this fortran code. The DCMPRS macro version is more 
  291. C       than twice as fast as the fortran version)
  292. C
  293. C_ARGS
  294.       INTEGER*4       NIN             
  295. C                                       Input, # bytes in input buffer
  296.       INTEGER*4       NOUT            
  297. C                                       Input, # bytes in output line
  298.       INTEGER*2       TREE(5,1024)    
  299. C                                       Input, Huffman code tree
  300.       CHARACTER       IBUF(1)         
  301. C                                       Input, compressed data buffer
  302.       CHARACTER       OBUF(1)         
  303. C                                       Output, decompressed line
  304. C
  305. C_DESC  This subroutine follows a path from the root of the Huffman tree
  306. C       to one of its leaves.  The choice at each branch is decided by
  307. C       the successive bits of the compressed input stream, left for 1,
  308. C       right for 0.  Only leaf nodes have a value other than -1.  The
  309. C       routine traces a path through the tree until it finds a node 
  310. C       with a value not equal to -1 (a leaf node).  The value at the 
  311. C       leaf node is subtracted from the preceding pixel value plus 256 
  312. C       to restore the the umcompressed pixel.  This algorithm is then 
  313. C       repeated until the entire line has been processed.
  314. C
  315. C_CALLS This routine uses the VAX/VMS implicit function
  316. C       BTEST (VAR, IBIT) to test whether the bit number IBIT is set 
  317. C       in the variable VAR. A fortran version of this routine
  318. C       is available for BTEST.
  319. C
  320. C_HIST  DEC86, DMcMacken, ISD, USGS, Flagstaff, Original version
  321. C
  322. C_END
  323. C**********************************************************************
  324.       INTEGER*4       PTR             
  325. C                                       Pointer to position in tree
  326.       INTEGER*4       K               
  327. C                                       Do loop index
  328.       INTEGER*4       BIT             
  329. C                                       Pointer to current bit in word
  330.       INTEGER*4       NBYTE           
  331. C                                       Counter, # bytes decompressed
  332.       INTEGER*4       TEMP            
  333. C                                       Working variable
  334.       INTEGER*2       VALUE           
  335. C                                       Value index in tree
  336.       INTEGER*2       RIGHT           
  337. C                                       Right branch index in tree
  338.       INTEGER*2       LEFT            
  339. C                                       Left branch index in tree
  340.       INTEGER*2       PARENT          
  341. C                                       Parent pointer index in tree
  342.       INTEGER*2       BRANCH          
  343. C                                       Parent branch index in tree
  344.       PARAMETER       (VALUE=1, RIGHT=2, LEFT=3, PARENT=4, BRANCH=5)
  345. C
  346.       INTEGER*2       DN              
  347. C                                       Current density value
  348.       INTEGER*2       DNL             
  349. C                                       Last density value
  350.       LOGICAL       END             
  351. C                                       End of line flag
  352.       LOGICAL       BTEST
  353. C
  354. C**********************************************************************
  355. C Start at root of tree
  356. C**********************************************************************
  357.       PTR = TREE(1,1024)
  358. C**********************************************************************
  359. C Just copy first byte it is uncompressed.
  360. C**********************************************************************
  361.       OBUF(1) = IBUF(1)
  362. C**********************************************************************
  363. C Count starts here
  364. C**********************************************************************
  365.       NBYTE = 1
  366. C**********************************************************************
  367. C Initialize current and last density values
  368. C**********************************************************************
  369.       DN = ICHAR(IBUF(1))
  370.       IF (DN.LT.0) DN = 256 + DN
  371.       DNL = DN
  372. C**********************************************************************
  373. C Reset end of line flag
  374. C**********************************************************************
  375.       END = .FALSE.
  376. C**********************************************************************
  377. C K points to current input byte loop is done when k exceeds # 
  378. C input bytes.
  379. C**********************************************************************
  380.       K = 2
  381.       IF (K .GT. NIN) END = .TRUE.
  382. C**********************************************************************
  383. C This is the processing loop
  384. C**********************************************************************
  385.    10 IF (.NOT. END) THEN
  386.          TEMP = ICHAR(IBUF(K))
  387.          BIT = 7                         
  388.    20    IF (BIT .GE. 0 .AND. .NOT. END) THEN 
  389.             IF (BTEST(TEMP, BIT)) THEN
  390.                PTR = TREE(LEFT, PTR)
  391.             ELSE
  392.                PTR = TREE(RIGHT, PTR)
  393.             ENDIF
  394. C**********************************************************************
  395. C We are at end leaf when value is not -1
  396. C**********************************************************************
  397.             IF (TREE(VALUE, PTR) .NE. -1) THEN
  398. C**********************************************************************
  399. C Compute value of pixel, reset or increment pointers and counters
  400. C**********************************************************************
  401.                DN = TREE(VALUE, PTR)
  402.                DN = DNL - DN + 256
  403.                DNL = DN
  404.                NBYTE = NBYTE + 1
  405.                OBUF(NBYTE) = CHAR(DN)
  406.                PTR = TREE(1, 1024)
  407.                IF (NBYTE .EQ. NOUT) END = .TRUE.
  408.             ENDIF
  409. C**********************************************************************
  410. C Process next bit in byte
  411. C**********************************************************************
  412.             BIT = BIT - 1
  413.             GOTO 20
  414.          ENDIF
  415. C**********************************************************************
  416. C Set index to next input byte
  417. C**********************************************************************
  418.          K = K + 1
  419.          IF (K .GT. NIN) END = .TRUE.    
  420.          GOTO 10
  421.       ENDIF
  422. C**********************************************************************
  423. C Go back to caller with decompressed line
  424. C**********************************************************************
  425.       RETURN
  426.       END
  427.