Skip to content

使用 Motion JPEG2000 压缩高光谱数据

原文链接: https://www.nv5geospatialsoftware.com/Learn/Blogs/Blog-Details/compressing-hyperspectral-data-using-motion-jpeg2000

11302 评价本文:

尚未评分

使用 Motion JPEG2000 压缩高光谱数据

匿名作者 2016年8月12日，星期五

正如我之前在博客中提到的，你必须关注栅格数据的交错方式（interleave），以便高效地遍历数据。那篇博文是特定于 ENVI 的，而这次我要介绍一种使用纯 IDL 存储数据立方体（data cube）以支持高效 BIL 模式访问的方法。

随着越来越多的高光谱（HSI）传感器被制造出来，具有数百甚至数千个波段的数据立方体将变得更加普遍。这些传感器通常是推扫式扫描仪（pushbroom scanners），具有固定的样本数和波段数，但采集的行数变化很大。这种采集模式天然适合 BIP 和 BIL 存储交错方式，因此每个像素扫描行都可以写入磁盘，而无需使用大量内存来保存多行数据以进行 BSQ 或分块输出。如果你使用 HSI 数据，你可能也对你分析中的大部分（如果不是全部）波段感兴趣，因此 BIP 和 BIL 交错方式对于这些数据访问模式来说效率更高。

HSI 数据的挑战在于，由于波段数量多，数据量可能很大。某种形式的数据压缩会很有用，前提是它是无损的——如果你打算通过压缩或插值伪影丢失细节，那么花钱收集所有这些波段还有什么意义呢？

在 IDL 中压缩数据立方体有几种选择：

zip 选项的缺点在于你必须解压缩整个文件才能使用它，因此没有任何空间节省，实际上比使用原始的 ENVI 文件更糟糕。IDL 对 TIFF 格式的支持提供了 4 种不同的压缩模式：

  • 无
  • LZW
  • PackBits
  • JPEG

“无”选项毫无意义，因为它不执行任何压缩。JPEG 选项是有损的，会引入伪影，但首先它不能用于此类数据。第一，它只支持 8 位值，因此任何其他数据类型都会失败。第二，如果你向该方法发送超过 3 个波段的数据，你会收到一个奇怪的错误信息“Bogus input colorspace”。此消息实际上来自 libTIFF 库，试图告诉你它无法对超过 3 个波段的数据进行 JPEG 压缩。

IDLffMJPEG2000 对象是 IDL 对 Motion JPEG2000 压缩格式的接口。Motion JPEG2000 标准是 JPEG2000 的扩展，JPEG2000 是基于小波的、备受尊敬的 JPEG 格式的后继者。JPEG2000 格式同时支持有损和无损压缩，我们在这里将使用后者，以避免向数据引入任何伪影。与制作视频格式不同，Motion JPEG2000 只执行帧内压缩，不进行帧间压缩。虽然这可能导致比其他格式更低的压缩比，但也意味着解压缩动画的单个帧更快，因为它不需要为帧间插值加载其他帧。创建 Motion JPEG2000 文件时有几个可调节的参数——BIT_RATE、N_LAYERS、N_LEVELS、PROGRESSION 和 TILE_DIMENSIONS。我没有尝试改变这些值如何影响压缩比或解压缩所需的时间，所以你自己尝试调整这些值可能是有价值的。

首先，我将展示用于执行不同压缩格式创建的代码，然后展示几个不同数据立方体的压缩比和压缩时间。

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

IDLffMJPEG2000 类要求你一次添加一帧数据（或一帧的子集），因此需要循环遍历栅格的行数。我必须在其中加入一个 Reform() 调用，以便将 3D 切片降维为 2D，因为最后一个维度始终是 1。这个对象实际上是线程化的，所以 SetData() 调用不会阻塞，它们将数据排队以压缩到输出文件中，因此 Commit() 函数有一个可选的超时参数，即等待多少毫秒后返回。在我的测试中，提交几乎不需要什么时间，处理器能跟上压缩请求。

现在来看有趣的部分，压缩性能：

传感器 行数 原始大小 GZIP ZIP TIFF LZW TIFF PackBits MJP2
Hyperion
254 列
242 波段		 1310 161046160 96668737
60.0%
7.2 秒		 96668847
60.0%
10.6 秒		 118720382
73.7 %
4.8 秒		 133119274
82.7 %
2.8 秒		 64289348
39.9 %
13.2 秒
AVIRIS
614 列
224 波段		 1024 281674956 213266080
75.7 %
15.2 秒		 213266228
75.7 %
24.0 秒		 271876484
96.5 %
7.8 秒		 277500090
98.5 %
3.8 %		 129922225
60.9 %
19.5 秒
AVIRIS
952 列
224 波段		 3830 1633479680 784610303
48.0 %
60.5 秒		 784610437
48.0 %
92.3 秒		 1015098818
62.1 %
44.8 秒		 1076409756
65.9 %
25.4 秒		 486091874
29.8 %
109.1 秒

GZIP 和 ZIP 结果几乎相同并不奇怪，尽管 GZIP 速度要快得多。TIFF 的 LZW 模式比 PackBits 模式更高效但更慢也不奇怪。有点令人惊讶的是，LZW 压缩的效果比 GZIP 和 ZIP 差那么多，而 Motion JPEG2000 比 GZIP 和 ZIP 稍微好一些。

现在，假设你想读取其中一个 Motion JPEG2000 压缩文件，但只需要其中的空间和/或光谱子集。这里有一个函数可以做到这一点：

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

我使用了与 ENVIRaster::GetData() 相同的关键字解释：

  • SUB_RECT = [x1, y1, x2, y2]
  • BANDS – 任何不重复的波段索引数组

使用 ENVI ROI 进行基本特征提取 最小面积边界框