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:
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
end
pro compress_to_tiff, raster, file
compile_opt idl2
dat = raster.GetData()
WRITE_TIFF, file, dat, /SIGNED, /SHORT,
COMPRESSION=1
end
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
endelse
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]))
endfor
print, mj2k.Commit(30000)
end
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
|
Hyperion
254 cols
242 bands
|
1310
|
161046160
|
96668737
60.0%
7.2 sec
|
96668847
60.0%
10.6 sec
|
118720382
73.7 %
4.8 sec
|
133119274
82.7 %
2.8 sec
|
64289348
39.9 %
13.2 sec
|
AVIRIS
614 cols
224 bands
|
1024
|
281674956
|
213266080
75.7 %
15.2 sec
|
213266228
75.7 %
24.0 sec
|
271876484
96.5 %
7.8 sec
|
277500090
98.5 %
3.8 %
|
129922225
60.9 %
Sec
|
AVIRIS
952 cols
224 bands
|
3830
|
1633479680
|
784610303
48.0 %
60.5 sec
|
784610437
48.0 %
92.3 sec
|
1015098818
62.1 %
44.8 sec
|
1076409756
65.9 %
25.4 sec
|
486091874
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]
endif
if (~ISA(bands))
then begin
bands = INDGEN(dims[1])
endif
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
endelse
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]
endfor
return, outData
end
I used the same keyword interpretation as
ENVIRaster::GetData():
- SUB_RECT = [x1, y1, x2, y2]
- BANDS – any array of nonrepeating band indices