9241 Rate this article:
No rating

Compressing hyperspectral data using Motion JPEG2000


As I’ve blogged about before, you have to pay attention to a raster’s interleave so that you can iterate over the data efficiently.  That post was ENVI specific, this time I’m here to present a way to store a data cube for efficient BIL mode access that is pure IDL.

As more hyperspectral (HSI) sensors are produced, data cubes with hundreds or even thousands of bands will become more common.  These sensors are usually pushbroom scanners, with a fixed number of samples and bands and a highly variable number of lines collected.  The collection modality lends itself to the BIP and BIL storage interleaves, so each scanline of pixels can be written to disk without having to use large amounts of memory to hold many lines for BSQ or tiled output.  If you are using HSI data, you are probably also interested in most if not all of the bands for your analysis, so BIP and BIL interleaves are more efficient for those data access patterns.

The challenge of HSI data is that it can be large, due to the high number of bands.  Some form of data compression would be useful, as long as it is lossless – why spend the money collecting all those bands if you are going to lose the details with compression or interpolation artifacts?

There a few options for compressing a datacube in IDL:

The zip option is not attractive in that you have to unzip the entire file to use it, so there isn’t any space savings, in fact it’s worse than using the original ENVI file.  The IDL support for the TIFF format allows for 4 different compression modes:

  • None
  • LZW
  • PackBits
  • JPEG

The None option is pointless, since it doesn’t perform any compression.  The JPEG option is lossy and would introduce artifacts, but it can’t be used for this data in the first place.  Firstly, it only supports 8-bit values, so any other datatypes will fail.  Secondly, if you send more than 3 bands of data into the method, you get back the curious error message “Bogus input colorspace”.  This message actually comes from the libTIFF library, trying to tell you it can’t JPEG compress more than 3 band data.

The IDLffMJPEG2000 object is IDL’s interface to the Motion JPEG2000 compression format.  The Motion JPEG2000 standard is an extension to JPEG2000, the wavelet based successor to the venerable JPEG format.  The JPEG2000 format supports both lossy and lossless compression, and we’ll use the latter here to avoid introducing any artifacts to the data.  Unlike make video formats, Motion JPEG2000 only performs intraframe compression, no interframe compression.  While this likely results in lower compression ratios than those other formats, it also means that decompressing a single frame of the animation is faster because it doesn’t need to load other frames for the interframe interpolation.  There are a few degrees of freedom available when creating Motion JPEG2000 files – BIT_RATE, N_LAYERS, N_LEVELS, PROGRESSION, and TILE_DIMENSIONS.  I did not experiment with how changing these values effects the compression ratio or time needed to decompress, so playing with these one your own may be worthwhile.

First I’ll present the code I used to perform the creation of the different compressed formats , and then show the compression ratios and time to compress for a couple different datacubes.

pro compress_to_zip, raster, file   compile_opt idl2
  FILE_GZIP, raster.URI, file
pro compress_to_tiff, raster, file
  compile_opt idl2
  dat = raster.GetData()
pro compress_to_mj2k, raster, file
  compile_opt idl2
  dat = raster.GetData()
  signed = !True
  if (raster.Data_Type eq 'byte') then begin
    depth = 8
  endif else if (raster.Data_Type eq 'int') then begin
    depth = 16
  endif else if (raster.Data_Type eq 'long') then begin
    depth = 24
  endif else if (raster.Data_Type eq 'uint') then begin
    depth = 16
    signed = !False
  endif else if (raster.Data_Type eq 'ulong') then begin
    depth = 24
    signed = !False
  endif else begin
    Message, 'Unsupported pixel datatype : ' + raster.Data_Type
  mj2k = IDLffMJPEG2000(file, /WRITE, /REVERSIBLE, $
                        BIT_DEPTH=depth, SIGNED=signed, $
                        DIMENSIONS=[raster.nColumns, raster.nBands])
  for row = 0, raster.nRows-1 do begin
    !null = mj2k.SetData(Reform(dat[*, *, row]))
  print, mj2k.Commit(30000)


The IDLffMJPEG2000 class requires you to add data one frame at a time (or subset of a frame), hence the for loop over the number of rows in the raster.  I did have to stick a Reform() call in there, to collapse the 3D slice down to 2D, since the last dimension is always 1.  This object is actually threaded, so the SetData() calls don’t block, they queue up the data to compressed into the output file, so the Commit() function has an optional timeout argument, which is how many milliseconds to wait before returning.  In my testing it doesn’t take much time at all to do the commit, the processor can keep up with the compression requests.

Now the interesting part, compression performance:


Sensor Rows Size GZIP ZIP TIFF LZW TIFF PackBits MJP2
254 cols
242 bands
1310 161046160 96668737
7.2 sec
10.6 sec
73.7 %
4.8 sec
82.7 %
2.8 sec
39.9 %
13.2 sec
614 cols
224 bands
1024 281674956 213266080
75.7 %
15.2 sec
75.7 %
24.0 sec
96.5 %
7.8 sec
98.5 %
3.8 %
60.9 %
952 cols
224 bands
3830 1633479680 784610303
48.0 %
60.5 sec
48.0 %
92.3 sec
62.1 %
44.8 sec
65.9 %
25.4 sec
29.8 %
109.1 sec


It isn’t surprising that the GZIP and ZIP results are almost identical, though GZIP is way faster.  It also isn’t surprising that the LZW mode for TIFF is more efficient than the PackBits mode, though slower.  It is somewhat surprising that the LZW compression was so much work than GZIP and ZIP, and that Motion JPEG2000 is a bit better than GZIP and ZIP.

Now let’s say you wanted to read one of the Motion JPEG2000 compressed files back, but only need a spatial and/or spectral subset of it.  Here is a function that will do just that:

function read_mj2k, file, SUB_RECT=subRect, BANDS=bands
  compile_opt idl2
  mj2k = IDLffMJPEG2000(file)
  mj2k.GetProperty, DIMENSIONS=dims, N_FRAMES=nFrames, $
                    BIT_DEPTH=depth, SIGNED=signed
  if (~ISA(subRect)) then begin
    subRect = [0, 0, dims[0]-1, nFrames-1]
  if (~ISA(bands)) then begin
    bands = INDGEN(dims[1])
  startFrame = subRect[1]
  endFrame = subRect[3]
  leftCol = subRect[0]
  rightCol = subRect[2]
  topRow = MIN(bands)
  bottomRow = MAX(bands)
  if (depth le 8) then begin
    type = 1 ; Byte
  endif else if (depth le 16) then begin
    type = signed ? 2 : 12 ; Int or UInt
  endif else begin
    type = signed ? 3 : 13 ; Long or ULong
  outData = MAKE_ARRAY(rightCol-leftCol+1, N_ELEMENTS(bands), $
                       endFrame-startFrame+1, TYPE=type)
  region = [leftCol, 0, rightCol-leftCol+1, dims[1]]
  for frame = 0, endFrame-startFrame do begin
    dat = mj2k.GetData(frame+startFrame, REGION=region)
    outData[*, *, frame] = dat[*, bands]
  return, outData


I used the same keyword interpretation as ENVIRaster::GetData():

  • SUB_RECT = [x1, y1, x2, y2]
  • BANDS – any array of nonrepeating band indices