像素交织——为何需要关注及处理方法
20682 评价本文：
未评分
像素交织——为何需要关注及处理方法
匿名 2015年6月19日 星期五
像素交织（pixel interleave）看起来只是栅格数据中一个不起眼的元数据项，但如果你不根据它来调整数据访问模式，它可能会对你的算法性能产生巨大影响。
我将首先快速回顾什么是像素交织，然后展示它如何影响性能，最后介绍如何更新你的代码以优化数据访问模式并最小化这种影响。
像素交织有三种形式——按像素波段交织（Band Interleaved by Pixel, BIP）、按行波段交织（Band Interleaved by Line, BIL）和波段序列（Band Sequential, BSQ）。ENVIRaster::GetData()的帮助文档中对此进行了讨论，因为你可以请求以特定的交织方式返回数据。如果你用不同的 INTERLEAVE 值调用 GetData()，你将得到相同的数据，但维度不同。BIP 交织最适合光谱处理，因为给定像素的所有波段值都是相邻的。BSQ 交织最适合空间处理，因为同一波段中的相邻像素是紧挨着的。BIL 交织是前两者之间的折衷方案，给定扫描行中的所有像素是相邻的，一次一个波段。
了解数据在磁盘上的存储方式，可以使请求相同像素集合的速度变得更快或更慢，这取决于你是否以错误的交织方式请求像素。让我们看看 ENVI 附带的 qb_boulder_msi 文件。它有 1024 列、1024 行和 4 个波段，并以 BSQ 交织方式存储。这意味着它在磁盘上存储为一个 [1024, 1024, 4] 数组，像素流存储为 [0, 0, 0], [1, 0, 0], [2, 0, 0], … , [1023, 0, 0], [0, 1, 0], [1, 1, 0], … , [1022, 1023, 0], [1023, 1023, 0], [0, 0, 1], [1, 0, 1], … , [1022, 1023, 3], [1023, 1023, 3]。如果你以 BSQ 交织方式请求数据，那么我们只需为每个波段读取一个子矩形区域，并逐片组装数据立方体。如果你以 BIL 或 BIP 交织方式请求数据，我们可以以类似的方式读取，但需要在复制到数据立方体之前对数据进行旋转。
当你想要使用 ENVI API 编写自己的算法时（最好封装在 ENVITaskFromProcedure 中），你必须对如何读取和处理数据做出一些决定。如果你幸运的话，可以将所有像素读入一个数组，然后随意迭代，但如果你想让代码更健壮，应该计划以瓦片（tile）方式访问数据，这样它可以处理任意维度的栅格。为了在这方面帮助你，ENVI API 包含了 ENVIRasterIterator 类，它可以通过 ENVIRaster::CreateIterator() 函数为你构建。这个函数有一个重要的关键字 MODE，用于控制生成的瓦片的形状以及它们遍历栅格的方式。
MODE 的默认值是 'spatial'，它只返回一个波段的像素。迭代这种类型的迭代器将以行主序返回跨越第一个波段的数据瓦片——瓦片从栅格/请求的 SUB_RECT 的左上角开始，从左向右移动到右边缘，然后向下移动一行并重复，直到覆盖第一个波段的所有像素，然后对每个后续波段遵循相同的模式。迭代器使用的默认瓦片大小是 1024x1024，但为了强调这一点，让我们考虑 qb_boulder_msi 上的一个 512x512 像素迭代器。第一个瓦片将覆盖波段 0 中从 [0,0] 到 [511,511] 的子矩形，然后是波段 0 中从 [512,0] 到 [1023,511]，接着是波段 0 中从 [0,512] 到 [511,1023]，然后是波段 0 中从 [512,512] 到 [1023, 1023]，接着是波段 1 中相同的 4 个子矩形序列，然后是波段 2，最后是波段 3。
MODE 的替代值是 'spectral'，它返回类似于 BIL 扫描行的像素。每个瓦片的列数与栅格相同，行数是栅格中的波段数。以这种方式迭代栅格会产生与栅格像素行数相同的瓦片数量。
我整理了以下简单过程来说明交织对性能的影响。它执行类似于光谱角制图（Spectral Angle Mapper, SAM）的计算，你传入一个 ENVIRaster 和一个参考光谱，并返回一个二维数组值。它做的第一件事是通过将每个元素除以数组与其自身点积的平方根来归一化输入光谱。然后我创建输出数组，其空间维度与输入栅格相同，以及一个名为 ones 的辅助数组，用于创建矩阵外积，这样我们就得到一个与瓦片大小相同的二维数组，其中每一行都相同。接着我创建迭代器，由于我正在做光谱点积，所以使用光谱模式，并用 foreach 进行迭代。正如 ENVIRasterIterator 文档所示，对此对象使用 foreach 将导致循环变量被设置为当前瓦片中的像素。我取出 nColumns * nBands 的瓦片并将其平方，然后对结果调用 Total() 以对每列求和，再取平方根。结果是一个一维数组，包含每个像素跨其波段的幅度。我需要将这个一维数组转换回二维数组，这是通过与该 ones 数组进行外积（# 运算符）来实现的，并用它除以原始瓦片。结果 normTile 现在是瓦片的归一化版本，其中每个像素的值在平方并求和后将归一化为 1。然后我取这个归一化瓦片与归一化光谱的点积，产生一个一维数组值，本质上是每个像素的光谱与参考光谱的点积。最后，我将这行像素复制到输出数组中的适当位置。
pro spectralDivergence, RASTER=oRaster, SPECTRUM=spectrumIn, RESULT=result
compile_opt idl2
refSpectrum = spectrumIn / Sqrt(Total(spectrumIn*spectrumIn))
result = DblArr(oRaster.nColumns, oRaster.nRows)
ones = Make_Array(oRaster.nBands, /DOUBLE, Value=1)
iter = oRaster.CreateTileIterator(MODE='spectral')
foreach tile, iter do begin
tileSq = Double(tile) * tile
norms = Sqrt(Total(tileSq, 2))
normTile = tile / (norms # ones)
tileDot = normTile # refSpectrum
result[*, iter.Current_Subrect[1]] = tileDot
endforeach
end
我使用“转换交织”（Convert Interleave）工具箱工具将 qb_boulder_msi 转换为 BIL 和 BIP 交织，然后对全部 3 个文件运行此代码。使用 tic/toc 的计时结果如下表：
| 交织 | 本地时间 (秒) | 网络时间 (秒) |
|---|---|---|
| BIP | 0.750 | 1.417 |
| BIL | 0.706 | 0.730 |
| BSQ | 0.738 | 1.671 |
我分别在本地硬盘和网络共享上运行了计时测试，后者更能代表大型数据仓库的情况。对于本地文件，结果看起来差异不大，大约在 5% 左右，但对于网络文件，差异达到 2 倍或更多。
让我们仔细看看数据在磁盘上的存储方式以及访问方式，以了解其原因。在 BIP 交织中，数据首先按波段存储，然后按每行的列存储，最后按行存储。因此，当我们为瓦片请求单行扫描像素时，我们可以从磁盘读取一个数据块。但这个数据块的维度是 nBands * nColumns，而我们想要的是 nColumns * nBands，所以数据需要转置。在上面详细说明的 BSQ 交织中，数据首先按列存储，然后按行存储，最后按波段存储。因此，当我们为瓦片请求单行扫描像素时，我们必须从每个波段读取一行扫描线，然后寻找到下一个波段并读取它，所以我们最终会有 nBands 次读取和寻道，而不是各一次。这就是为什么我们看到网络文件的性能下降如此显著，所有这些单独的寻道和读取累加起来，并且无法像本地文件那样被有效地缓存。
所以问题是如何调整你的代码以考虑像素交织，并减少数据访问中所有耗时的部分。
对于 BIP 交织，改动并不大，我们只是不能使用 ENVIRasterIterator，因为它会强制进行 BIL 请求。因此，我们手动跟踪当前像素行的 SUB_RECT 并将其传递给 ENVIRaster::GetData()。这将返回一个与原始版本转置的瓦片，因此我们将 # 矩阵运算符替换为 ## 运算符，以便在其他维度上执行计算。由于我们仍然为每行像素产生一个一维输出数组，所以我们永远不需要转置任何变量。
对于 BSQ 交织，我们需要做更多工作。由于我们试图避免每行扫描线进行 nBands 次寻道和读取，我们将一次性读取多行扫描线到一个三维数据数组中，然后从立方体中逐行处理扫描线，直到需要读取下一个立方体。那么最大的问题是一次读取多少行扫描线，这变成了时间和空间的权衡——你能负担多少内存用于这个立方体，以避免所有这些额外寻道和读取的开销？在我的例子中，我选择了一个我知道小到足以强制读取多个立方体的数字。我选择了一个要读取的最大像素数（五十万），然后除以列数和波段数以得到行数（对于 qb_boulder_msi 是 128）。我首先读取一个数据立方体，然后迭代所有扫描行。在每次迭代开始时，我检查当前行是否超过了最后一个 SUB_RECT，如果是，我将 SUB_RECT 向下平移，根据总行数进行钳位，并读取下一个数据立方体。然后我复制出一行扫描线的数据切片，它需要从一个 [1024, 1, 4] 数组重塑为一个 [1024, 4] 的二维数组，然后像 BIL 情况一样进行后续计算。
以下是更新后的代码：
pro spectralDivergence2, RASTER=oRaster, SPECTRUM=spectrumIn, RESULT=result
compile_opt idl2
refSpectrum = spectrumIn / Sqrt(Total(spectrumIn*spectrumIn))
result = DblArr(oRaster.nColumns, oRaster.nRows)
ones = Make_Array(oRaster.nBands, /DOUBLE, Value=1)
switch oRaster.Interleave of
'bip' : begin
for row = 0, oRaster.nRows-1 do begin
subRect = [0, row, oRaster.nColumns-1, row]
tile = oRaster.GetData(SUB_RECT=subRect)
tileSq = Double(tile) * tile
norms = Sqrt(Total(tileSq, 1))
normTile = tile / (norms ## ones)
tileDot = normTile ## refSpectrum
result[*, row] = Reform(tileDot)
endfor
break
end
'bil' : begin
iter = oRaster.CreateTileIterator(MODE='spectral')
foreach tile, iter do begin
tileSq = Double(tile) * tile
norms = Sqrt(Total(tileSq, 2))
normTile = tile / (norms # ones)
tileDot = normTile # refSpectrum
result[*, iter.Current_Subrect[1]] = tileDot
endforeach
break
end
'bsq' : begin
MaxTileSize = 512 * 1024L ; half a million pixels
numRows = Floor(MaxTileSize / Double(oRaster.nColumns * oRaster.nBands))
subRect = [0, 0, oRaster.nColumns-1, numRows-1]
bsqCube = oRaster.GetData(SUB_RECT=subRect)
for row = 0, oRaster.nRows-1 do begin
if (row gt subRect[-1]) then begin
subRect += [0, numRows, 0, numRows]
subRect[3] <= oRaster.nRows-1
bsqCube = oRaster.GetData(SUB_RECT=subRect)
endif
tile = Reform(bsqCube[*, row - subRect[1], *])
tileSq = Double(tile) * tile
norms = Sqrt(Total(tileSq, 2))
normTile = tile / (norms # ones)
tileDot = normTile # refSpectrum
result[*, row] = Reform(tileDot)
endfor
break
end
endswitch
end
更重要的是新的性能数据：
| 交织 | 本地时间 (秒) | 网络时间 (秒) |
|---|---|---|
| BIP | 0.648 | 1.277 |
| BIL | 0.720 | 0.740 |
| BSQ | 0.146 | 1.008 |
这些数字乍一看可能不太合理，但有几件事需要考虑。BIL 的数字差异不大，在我测量的标准偏差范围内，因为我为每个数据点多次运行了代码。BIP 的数字显示了一些改进，这很可能是由于手动更新 SUB_RECT 的逻辑比通过 ENVIRasterIterator 的 foreach 便利逻辑更简单。然而，BSQ 的数字最有趣。本地文件速度要快得多，快了 5 倍，但网络文件只快了 40%。如果我们分析整个算法的寻道和读取次数，我们会看到瓦片迭代器为每行扫描线执行 nBands 次寻道和读取，总共是 4 * 1024 = 4096 次寻道和读取，而缓存立方体读取版本为每个立方体执行 nBands 次寻道和读取，总共是 4 * 8 = 32 次寻道和读取。对于本地文件，由于可能发生的磁盘缓存，这带来了巨大的差异，而网络共享由于局域网中其他系统的文件请求，无法执行同样多的缓存。
你的实际效果将取决于所涉及的硬件类型及其负载，但从这里学到的重要一点是，通过一些预先考虑和一点点额外的代码，你可以仅仅因为改变了从存储介质获取数据到算法中的方式，就大大加速你的算法。