立体重建¶
本节介绍如何将成对的立体图像转换为高程信息。
可用的地形重建的标准问题 OTB Applications 包含以下步骤：
- 极线几何变换的位移网格估计
- 使用这些网格对图像对进行核线重采样
- 密集视差图估计
- 视差在数字表面模型(DSM)上的投影
让我们进入第三维空间吧！
估计对极几何变换¶
这一步的目的是生成重采样网格，以将图像转换为核线几何。 Epipolar geometry 是立体视觉的几何学。立体校正的操作确定要应用于每个图像的变换，使得成对的共轭核线变为共线、平行于其中一个图像轴并且对齐。在此几何图形中，左侧图像的给定行上显示的对象也位于右侧图像的同一行上。
应用这种变换将高程(或立体声对应关系的确定)问题简化为一维问题。我们有两个传感器图像 image1 和 image2 在相同的区域(立体对)上，我们假设我们知道与每个图像相关联的定位函数(正向和反向)。
Forward功能允许从图像参考转到地理参考。对于第一个图像，将注意到此函数：
(长，横)=f_{1}(i，j，h)
哪里 h 是海拔假说， (i,j) 是图1中的像素坐标， (long,lat) 是地理坐标。正如您可以想象的那样，逆函数允许从地理坐标到图像几何图形。
对于第二幅图像，在这种情况下，反函数的表达式为：
(i，j)=f^{Inv}_{2}(经纬度，h)
联合使用来自图像对的正反函数，我们可以构造共局部化函数 g_{1 \rightarrow 2} 在第一个像素的位置和第二个像素的位置之间：
(i_{2}，j_{2})=g_{1right tarrow 2}(i_{1}，j_{1}，h)
此函数的表达式为：
g_{1 \rightarrow 2} (i_{1} , j_{1} , h) = f^{Inv}_{2} [ f_{1}(i_{1} , j_{1}, h) ]
表达式并不重要，您需要理解的是，如果我们能够确定图像1中给定像素的图像2中的相应像素，正如我们知道两个图像之间的共局部化函数的表达式一样，我们可以通过识别来确定关于海拔的信息(等式中的变量h)！
我们现在有了数学基础来理解如何通过检查两个2-D核线图像中对象的相对位置来提取3-D信息。
在VHR光学图像的情况下，两个极线网格的构造稍微复杂一些。这是因为大多数太空被动遥感使用推扫式传感器，该传感器对应于垂直于航天器飞行方向布置的一系列传感器。这种采集配置意味着立体声校正的策略略有不同 (see here )。
我们现在将解释如何使用 StereoRectificationGridGenerator 应用程序生成两个图像，这两个图像是 deformation grids 对极线几何中的两个图像进行重采样。
otbcli_StereoRectificationGridGenerator -io.inleft image1.tif
-io.inright image2.tif
-epi.elevation.default 50
-epi.step 10
-io.outleft grid_image1.tif
-io.outright grid_image2.tif
该应用程序估计要应用于两个输入图像中的每个像素的位移，以获得核线几何图形。该应用程序接受一个‘Step’参数来估计较粗网格上的位移。这里我们估计每10个像素的位移。这是因为在大多数情况下，使用一对VHR和两幅图像之间的小角度，这种网格非常平滑。此外，执行并不是 streamable 并且潜在地使用大量的内存。因此，通常以较粗的分辨率估计置换栅格是个好主意。
该应用程序以核线几何形式输出输出图像的大小。 Note these values ，我们将在下一步中使用它们对核极几何中的两个图像进行重新采样。
在我们的案例中，我们有：
Output parameters value:
epi.rectsizex: 4462
epi.rectsizey: 2951
epi.baseline: 0.2094
这个 epi.baseline 参数提供基线与传感器高度比(在文献中也称为B/H)的平均值(以每米像素为单位)。它可用于将差异近似转换为物理高程：
H=h_{ref}+frac{d}{B/H}
哪里 h_{REF} 是用于生成核极网的参考海拔高度(此处：50米)，以及 d 是图像1和图像2之间给定对象的视差值(以像素为单位)。
我们现在可以前进到核极几何中的重采样。
对极几何图形中的图像进行重采样¶
前一个应用程序生成两个位移网格。这个 GridBasedImageResampling 允许使用这些栅格对对极几何体中的两个输入图像进行重新采样。这些网格是中间结果，在大多数情况下，它们本身并不真正有用。这第二步 only 在于应用变换来对两个图像进行重新采样。这个应用程序显然可以在许多其他环境中使用。
生成极线图像的两个命令是：
otbcli_GridBasedImageResampling -io.in image1.tif
-io.out epi_image1.tif
-grid.in grid_image1.tif
-out.sizex 4462
-out.sizey 2951
otbcli_GridBasedImageResampling -io.in image2.tif
-io.out epi_image2.tif
-grid.in grid_image2.tif
-out.sizex 4462
-out.sizey 2951
如你所见，我们设置了 sizex 和 sizey 参数使用由 StereoRectificationGridGenerator 应用程序设置输出核线图像的大小。这两幅核线图像的大小应该相同。
图1：在契奥普斯金字塔上的极几何图形中重采样图像1和图像2的提取。©CNES 2012
我们在极线几何中获得了两幅图像，如 Figure 1 。请注意，该应用程序只允许使用 -out.ulx 和 -out.uly 参数。
视差估计：核线上的块匹配¶
最后，我们可以开始立体声通信查找过程！
事情正在变得有点复杂，但不要担心。首先，我们将描述 BlockMatching 申请。
极线几何中图像的重采样允许我们将搜索限制在一维直线上，而不是沿两个维度，但更重要的是，直线上的差异，即通过块匹配过程测量的直线上的偏移可以直接与局部高程联系起来
几乎完全的光谱 stereo correspondence algorithms 已经出版了，而且还在以显著的速度增加！这个 Orfeo ToolBox 为块匹配实施不同的本地策略：
- 平方距离和块匹配(SSD)
- 归一化互相关(NCC)
- Lp伪范数(LP)
另一个重要参数(在应用程序中是必需的！)就是差距的范围。从理论上讲，块匹配可以进行盲目探索和搜索立体对之间的无限范围的差异。我们现在需要评估将执行块匹配的一系列差异(在一般情况下从地球上最深的点， the Challenger Deep 。登上珠峰顶峰！)
我们故意夸大了这一点，但您可以想象，如果没有较小的范围，块匹配算法可能会花费大量时间。这就是为什么这些参数对于应用程序是必需的，因此我们需要手动估计它们。使用这两幅核线图像，这相当简单。
在我们的例子中，我们选择一个点 flat 区域。它的坐标是 [1525, 1970] 在极线图像1和 [1526, 1970] 在核线图像2中。然后我们选择更高区域上的第二个点(在我们的例子中是接近契奥普斯金字塔顶部的点！)。该像素的图像坐标为 [1661,1299] 在图1和 [1633,1300] 在图2中，我们检查了图1和图2中的列坐标之间的差异，以便推导出对水平勘探有用的视差间隔。在我们的例子中，这个间隔至少是 [-28,1] (视差范围的符号的约定是从图像1到图像2)。
请注意，使用中的外部DEM可以缩短该勘探间隔 StereoRectificationGridGenerator 申请。事实上，核线图像之间测量的差异是相对于计算核线网格时使用的参考高度(因此，定义了核线几何形状)。使用外部DEM应该会产生与参考高度偏差较小的极线图像，因此，差异接近于0。
关于垂直差异，我们在第一步中说过，我们将2D勘探问题归结为1D问题，但在一般情况下，这并不完全正确。由于视差误差，在垂直方向上可能存在小的差异(即，核线在垂直方向上表现出小的移位，大约1个像素)。事实上，沿着垂直方向的勘探通常比沿着水平方向的勘探要小。您还可以在核线对上估计它们(在本例中，我们使用的范围为 -1 至 1 )。
再说一次，注意最小和最大差异的符号(总是从Image1到Image2)。
对象的命令行 BlockMatching 应用程序为：
otbcli_BlockMatching -io.inleft epi_image1.tif
-io.inright epi_image2.tif
-io.out disparity_map_ncc.tif
-bm.minhd -45
-bm.maxhd 5
-bm.minvd -1
-bm.maxvd 1
-mask.inleft epi_mask_image1.tif
-mask.inright epi_mask_image2.tif
-io.outmetric 1
-bm.metric ncc
-bm.subpixel dichotomy
-bm.medianfilter.radius 5
-bm.medianfilter.incoherence 2.0
默认情况下，该应用程序创建一个两个波段的图像：水平和垂直差异。
这个 BlockMatching 应用程序提供了许多其他强大的功能，以提高输出视差图的质量。
以下是其中的一些功能：
io.outmetric ：如果激活最佳度量值图像，它将连接到输出图像(然后将有三个波段：水平视差、垂直视差和度量值)
bm.subpixel ：执行差异的亚像素估计
mask.inleft 和 mask.inright ：您可以指定一个无数据值，该值将丢弃具有该值的像素(例如，核线几何图形可以生成大部分具有黑色像素的图像)。此掩码可以使用 BandMath 应用程序：
otbcli_BandMath -il epi_image1.tif -out epi_mask_image1.tif -exp "im1b1<=0 ? 0 : 255"
otbcli_BandMath -il epi_image2.tif -out epi_mask_image2.tif -exp "im1b1<=0 ? 0 : 255"
mask.variancet ：块匹配算法很难在均匀区域上找到匹配。我们可以使用方差阈值来丢弃这些区域，从而加快计算时间。
bm.medianfilter.radius 和 bm.medianfilter.incoherence ：将中值过滤器应用于视差图。中值滤波属于非线性滤波器族。它用于平滑图像，而不会受到离群值或散粒噪声的影响。半径对应于计算中值的邻域。执行输入视差图和经中值滤波的视差图之间的不一致性的检测(绝对差大于其缺省值为1的阈值的情况)。必须在应用程序中定义这两个参数才能激活过滤器。
当然，所有这些参数都可以组合在一起来改善视差图。
图2：水平视差和最佳度量图
从视差到数字表面模型¶
使用前面的应用程序，我们评估了核线图像之间的差异。下一个(也是最后一个！)现在的步骤是将视差图转换为海拔信息，以生成海拔地图。它使用视差图(水平和垂直)作为输入，以产生具有规则采样的数字表面模型(DSM)。高程值是通过对每个匹配像素的“左-右”视线进行三角测量来计算的。当DSM单元上有多个可用高程时，将保留最高的高程。
首先，重要的一点是，修改由 BlockMatching 应用程序仅保留相关差异。为此，我们可以使用输出的最佳度量图像，并根据该值过滤视差。
例如，如果我们使用归一化互相关(NCC)，我们只能保留最佳度量值优于的差异 0.9 。低于此值的差异可被视为不准确，不会用于计算高程信息( -io.mask 参数可用于此目的)。
这种过滤可以很容易地用 OTB Applications 。
我们首先使用 BandMath 根据最佳度量值过滤差异的应用程序：
otbcli_BandMath -il disparity_map_ncc.tif
-out thres_disparity.tif uint8
-exp "im1b3>0.9 ? 255 : 0"
现在，我们可以使用 DisparityMapToElevationMap 应用程序从过滤后的视差图计算高程图。
otbcli_DisparityMapToElevationMap -io.in disparity_map_ncc.tif
-io.left image1.tif
-io.right image2.tif
-io.lgrid grid_image1.tif
-io.rgrid grid_image2.tif
-io.mask thres_disparity.tif
-io.out elevation_map.tif
-hmin 10
-hmax 400
-elev.default 50
它会生成WGS84中投影的高程地图(EPSG代码： 4326 )覆盖了立体声对覆盖的地面区域。像素值以米为单位表示。
图3：契奥普斯金字塔上方的高程图提取。
这个 Figure 3 显示了Cheops对的输出DEM。
一个应用程序在多立体声框架方案中统领它们¶
已使用All-in-One方法创建了融合一个或多个立体重建(S)的应用程序： StereoFramework 。它从一个或多个立体声对计算DSM。首先，用户必须选择他的输入数据，并使用 -input.co 字符串参数。每对图像由2个图像索引“a b”(从0开始)定义，中间用空格分隔。不同的对与Coma相连。例如，“0 1，0 2”将定义图像对“First with Second”和“First with Third”。如果不填，则图像按对处理(相当于使用“0 1，2 3，4 5”…)。除了通常的高程和投影参数外，主要参数还被分成组，详细说明如下：
output ：输出参数(DSM分辨率、NoData值、单元融合方法)
- 输出投影贴图选择。
- 输出DSM的空间采样距离，单位为米
- DSM空单元格填充浮点值(默认情况下为-32768)
- 在每个DSM细胞中选择融合策略(最大、最小、平均、最大)
- 输出DSM
- 输出DSM的范围
stereorect ：正、逆立体校正网格子采样参数
- 直接变形网格的步长(以像素为单位)
- 反极线网格的二次采样
bm ：块匹配参数。
- 块匹配指标选择(强大的固态硬盘、固态硬盘、NCC、LP标准)
- 匹配滤镜的块半径(以像素为单位， 2 默认情况下)
- 选定高程源以下的最小海拔高度(以米为单位，默认为-20.0)
- 选定高程源上方的最大海拔高度(以米为单位，默认为20.0)
postproc ：后处理参数
- 使用双射一致性。计算从右到左的相关性以验证从左到右的差异。如果找不到双射，则拒绝像素
- 使用中值差异过滤(默认情况下禁用)
- 使用块匹配度量输出丢弃低相关值的像素(默认情况下禁用，浮点值)
mask ：计算可选的中间掩码。
- 左侧输入图像的蒙版(所有配对必须具有相同的大小)
- 右输入图像的蒙版(所有配对必须具有相同的大小)
- 此参数允许丢弃其局部方差太小的像素。邻域的大小由半径参数给出。(默认情况下禁用)
立体重建的良好做法¶
在应用程序内部使用高度偏移参数来推导最小和最大水平视差探测，因此它们对计算时间有重要影响。建议选择距离要生成的DSM不太远的高程源(例如，SRTM高程模型)。因此，在极线几何图形中将考虑高程源的海拔高度，并且差异将显示高程偏移(如建筑物)。它允许您沿高程轴使用较小的勘探范围，导致沿水平方向的勘探差异较小，计算速度较快。
为了减少时间消耗，将所有传感器图像裁剪到相同的程度将是有用的。要做到这一点，最简单的方法是选择一幅图像作为参考，然后应用 ExtractROI 使用Fit模式选项在其他传感器图像上应用程序。
算法提纲¶
应用程序中使用了以下算法：对于每个传感器对
- 从立体对(正向和逆)计算极线变形网格
- 使用BCO插值器重采样为极几何图形
- 为每个核线图像创建遮罩：删除黑色边框并重新采样输入遮罩
- 使用块匹配算法计算水平视差
- 使用二分法将视差细化到亚像素精度
- 应用可选中值滤波
- 根据相关性分数(可选)和勘探范围过滤差异
- 平移传感器几何图形中的差异
- 将视差图转换为3D图
然后，将所有3D地图进行融合，以生成具有所需地理或地图投影和可参数化范围的DSM。