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