Compressing hyperspectral data using Motion JPEG2000
Anonym
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