一种基于GPU‑CPU协同的无人机影像快速镶嵌方法与流程
本发明涉及到无人机影像的快速镶嵌技术领域，具体地说，是一种基于GPU-CPU协同的无人机影像快速镶嵌方法。
背景技术：
无人机等低空遥感系统以其机动性强、获取影像分辨率高、硬件价格低廉及维护成本低等优势，在应急测绘等领域中扮演着重要的角色。无人机影像具有分辨率高、框幅小、重叠度高等特点，在利用大序列无人机影像来制作测区数字影像地图(Digital Orthophoto Map，DOM)时，传统方法是先对测区内所有影像进行正射纠正，然后对所有纠正后影像进行边缘提取、生成镶嵌线网络，最后进行镶嵌。传统模式中单片的全幅正射纠正需要对每个像素进行纠正耗时巨大、基于边缘提取生成镶嵌线网络的方法复杂度高，效率低下，同时先全幅纠正再镶嵌的作业模式会造成巨大数据冗余、时间和资源开销，很难满足应急测绘的需求。
在航空影像拍摄过程中，由于受到相机内部和外部因素的影响，如光学透镜成像的不均匀性、曝光时间、拍摄角度、大气条件、光照差异、地物属性等，会造成影像内部或影像之间色彩的不一致，如影像的亮度、色调均衡等。经过影像间匀光后相邻影像间的色彩不一致的问题得到了很大的改善，但是依旧存在一定的偏差。在两幅影像经过镶嵌线镶嵌之后，在镶嵌线的两侧依然存在着色彩过渡不自然的现象。传统减弱镶嵌线两侧色彩过渡不自然的方法有直接平均融合法和线性加权融合法。但是整个测区所有镶嵌线两侧缓冲区域的色彩融合计算任务是非常大的，如果仅用CPU对每条镶嵌线两侧的缓冲区每个像素灰度进行加权内插将非常的耗时。
2007年，NVIDIA推出CUDA(Compute Unified Device Architecture，统一计算设备结构)以来，越来越多的学者将GPU强大的并行能力成功地运用于遥感数据的快速处理中，逐渐形成了GPU-CPU协同处理技术。所谓GPU-CPU协同就是将CPU作为主机端，而GPU作为设备端，CPU主要负责逻辑性强的串行计算以及对GPU端的调度，GPU则主要负责高度并行化的计算，并将计算结果反馈给CPU端，两者各司其职，协同计算。基于GPU-CPU协同的航空影像的正射纠正和卫星影像的正射纠正较传统方法在效率上分别提升了60倍与110倍。
然而，影像镶嵌比正射纠正的算法和过程都更加复杂繁琐，如何设计优化影像镶嵌的算法流程，将GPU-CPU协同处理技术应用其中，提高大序列无人机影像镶嵌效率，对于应急测绘具有重大意义。
技术实现要素：
针对现有技术的不足，本发明的目的是提供一种基于GPU-CPU协同的无人机影像快速镶嵌方法，该方法从匀光后影像出发，先自动生成测区Voronoi图镶嵌线网络，将正射纠正嵌入影像镶嵌过程中，且只对每张影像有效区域进行纠正直接生成DOM全景图，能够解决大序列、高重叠度无人机影像纠正镶嵌效率慢，数据冗余资源开销大的问题。
为达到上述目的，本发明采用的技术方案如下：
一种基于GPU-CPU协同的无人机影像快速镶嵌方法，其关键在于按照以下步骤进行：
步骤1：数据预处理，CPU端计算测区所有影像像主点对应像底点的位置，以测区所有像底点自动生成Voronoi图镶嵌线网络；
步骤2：根据Voronoi图镶嵌线网络，CPU端计算每张影像的有效区域对应的局部DEM和原始影像局部纹理；
步骤3：CPU端计算镶嵌后DOM全景影像图的大小，并在硬盘上创建DOM全景影像存储指针；
步骤4：根据步骤2获得的局部DEM和原始影像局部纹理，以及影像对应的外参数，GPU端依次对测区每张影像的有效区域进行正射纠正，并将纠正结果拷出内存，按照步骤3创建的DOM全景影像存储指针写到对应的位置；
步骤5：对于正射纠正后的影像，GPU端依次对每条镶嵌线两侧的影像进行线性加权融合处理，并将融合处理后的纹理拷出内存，按照步骤3创建的DOM全景影像存储指针写到对应的位置，对原始镶嵌结果进行更新。
进一步的描述是，步骤1中所有影像像底点位置的计算方程式为：
其中，(X，Y)为影像像底点的位置，(x,y)为影像的像主点，Xs,Ys,Zs为影像外方位元素中的三个角元素，a1,b1,c1,a2,b2,c2,a3,b3,c3为影像定向结果中的旋转矩阵，f为相机焦距，Z为点(Xs,Ys)在数字高程模型DEM上内插的高程值。
进一步的描述是，步骤2中所述局部DEM和原始影像局部纹理的计算步骤为：
步骤2-1：根据步骤1得到的Voronoi图镶嵌线网络，得到当前影像的有效区域范围，即有效区域多边形PolygonS；
步骤2-2：遍历多边形PolygonS的各个角点从中选出坐标的四个极值Xmin,Xmax,Ymin,Ymax，得到多边形PolygonS的最小外接矩形ABCD，提取矩形ABCD所覆盖的局部DEM数据；
步骤2-3：将矩形ABCD四个角点按照公式反算到原始影像像方坐标系中得到四边形abcd,计算四边形abcd的最小外接矩形A'B'C'D'，得到原始影像上矩形A'B'C'D'所覆盖的影像纹理数据；其中，(XA,YA)为矩形ABCD的角点A，(xA,yA)为四边形abcd的角点a，Xs,Ys,Zs为影像外方位元素中的三个角元素，a1,b1,c1,a2,b2,c2,a3,b3,c3为影像定向结果中的旋转矩阵，f为相机焦距，ZA为点(XA,YA)在数字高程模型DEM 上内插的高程值。
进一步的描述是，步骤3中计算DOM全景影像图大小的具体步骤如下：
步骤3-1：将所有原始影像的四个角点投影到地面，得到每张影像纠正后四个方向的极值Ximin，Ximax，Yimin，Yimax，按照公式求得镶嵌后DOM影像在四个方向上的极值Xmin,Xmax,Ymin,Ymax，其中，min()为该一维数组中的最小值，max()为该一维数组中的最大值，i＝1～n，n为影像张数；
步骤3-2：按照计算镶嵌后DOM全景图的宽度width和高度height，M为纠正后正射影像的地面分辨率。
进一步的描述是，步骤4中对影像有效区域的正射纠正按照如下步骤进行：
步骤4-1：获取当前影像有效区域对应的局部DEM、原始影像局部纹理及当前影像对应的外参数；
步骤4-2：根据当前影像有效区域多边形PolygonS的最小外接矩形大小和纠正后正射影像的地面分辨率，确定当前任务GPU线程格网中线程块的大小和个数；
步骤4-3：GPU多线程并行，利用反解法数字微分纠正对该影像块进行正射纠正，利用线程索引计算该像素对应的地面点的坐标(X,Y)；
步骤4-4：判断点(X,Y)是否在当前影像有效区域多边形PolygonS内部，如果在内部则根据点(X,Y)在局部DEM块中双线性内插得到该点的高程Z，再计算该像素在原始影像块上的像点坐标p(x,y)，内插灰度值，否则不进行处理；
步骤4-5：线程按索引号进行灰度赋值；
步骤4-6：将纠正后影像从显存拷贝到内存并写入硬盘上创建好的DOM全景影像图的对应位置。
进一步的描述是，步骤5中对镶嵌线两侧影像线性加权融合处理的步骤为：
步骤5-1：计算当前镶嵌线对应的局部DEM以及该镶嵌线两侧影像上局部的原始纹理数据，并将其拷入显存；
步骤5-2：设置线程格网中线程块的大小和数量；
步骤5-3：所有线程并行运行，利用线程索引计算该当前像素对应的地面点的坐标(X,Y)与高程Z；
步骤5-4：线程按索引号计算点(X,Y)在两原始影像上像点灰度值g1、g2；
步骤5-5：线程按索引号对灰度值g1、g2进行加权内插得到最终灰度值G；
步骤5-6：线程按索引进行灰度赋值；
步骤5-7：将融合后的局部影像从显存拷贝到内存，并写入硬盘上创建好的DOM全景影像图的对应位置，覆盖原来影像对当前镶嵌线两侧的影像纹理进行更新。
本发明的显著效果是：
(1)在整个DOM全景影像图生产模式上，本方法突破了先单张全幅正射纠正，然后进行镶嵌线网络选择，再进行镶嵌的传统模式，从匀光后原始影像出发先计算镶嵌线网络及每张影像的镶嵌有效区域，再对每张影像的有效区域进行正射纠正，直接镶嵌输出测区DOM全景图。本方法没有中间数据产生，并且只对每张影像的有效区域进行正射纠正，大大减少了数据冗余以及时间和资源的开销；
(2)充分利用了正射纠正和镶嵌线两侧融合处理的高度并行化特点，成功将GPU-CPU协同处理技术应用于无人机影像的快速纠正并镶嵌，极大提高了从原始影像到DOM全景图的制作效率，较传统方法提高了十多倍；
(3)根据Voronoi图镶嵌线网络镶嵌后的DOM影像每个像素理论上投影变形最小，有效减少了中心投影造成的影像边缘投影变形大的影响。
附图说明
图1是本发明的方法流程图；
图2是基于投影光线法计算影像像底点位置的示意图；
图3是测区所有影像像底点分布示意图；
图4是Voronoi图镶嵌线网络示意图；
图5是影像有效区域正射纠正流程图；
图6是镶嵌线两侧线性加权融合示意图；
图7是单条镶嵌线两侧影像融合处理流程图；
图8是镶嵌后带镶嵌线的DOM全景图；
图9是镶嵌后不带镶嵌线的DOM全景图；
图10是经过镶嵌线两侧影像纹理线性加权融合前后对比图。
具体实施方式
下面结合附图对本发明的具体实施方式以及工作原理作进一步详细说明。
本实施例结合某地230张分辨率为0.1m、航向重叠约为70％，旁向重叠约为45％、框幅为3744*5616的无人机影像为例，对本发明方法进行详细说明，镶嵌后DOM全景图的框幅为76104*66850，大小为14.5GB。
如图1所示，一种基于GPU-CPU协同的无人机影像快速镶嵌方法，具体步骤如下：
步骤1：数据预处理，包括局部DEM、内外方位元素等参数读入；然后CPU端计算测区所有影像像主点对应像底点的位置，以测区所有像底点自动生成Voronoi图镶嵌线网络；
具体实施时，首先进行所有影像像底点位置的计算：
根据测区影像定向结果，利用投影光线法将影像的像主点(x,y)投影到地面，即得到该影像像底点的位置(X，Y)，所述投影光线法的投影方式如图2所示。所述影像像底点位置的计算方程式为：
其中，(X，Y)为影像像底点的位置，(x,y)为影像的像主点，Xs,Ys,Zs为影像外方位元素中的三个角元素，a1,b1,c1,a2,b2,c2,a3,b3,c3为影像定向结果中的旋转矩阵，f为相机焦距，Z为点(Xs,Ys)在数字高程模型DEM上内插的高程值。
将测区内的所有影像像主点投影到地面，即得到本测区所有影像像底点的位置，如附图3所示。
然后，以测区所有像底点为点集，生成Voronoi图镶嵌线网络：
以测区所有像底点为散点集，先根据逐点内插法进行Delaunay三角剖分，再根据三角剖分结果生成Voronoi图镶嵌线网络，每张影像对应的泰森多边形即是该影像在镶嵌时有效区域范围，以下称为影像有效区域，对应的泰森多边形则称为影像有效区域多边形；同时得到每条镶嵌线与两侧影像之间的拓扑关系，即为本测区的Voronoi图镶嵌线网络，如附图4所示。
步骤2：根据Voronoi图镶嵌线网络，CPU端计算每张影像的有效区域对应的局部DEM和原始影像局部纹理，具体实施时，所述局部DEM和原始影像局部纹理的计算步骤如下：
进入步骤2-1：根据步骤1得到的Voronoi图镶嵌线网络，从第i张影像对应的有效多边形开始，得到当前影像的有效区域范围，即有效区域多边形PolygonS；
步骤2-2：遍历多边形PolygonS的各个角点从中选出坐标的四个极值Xmin,Xmax,Ymin,Ymax，得到多边形PolygonS的最小外接矩形ABCD，提取矩形ABCD所覆盖的局部DEM数据；
步骤2-3：将矩形ABCD四个角点按照公式反算到原始影像像方坐标系中得到四边形abcd,计算四边形abcd的最小外接矩形A'B'C'D'，得到原始影像上矩形A'B'C'D'所覆盖的局部影像纹理数据；其中，(XA，YA)为矩形ABCD的角点A，(xA，yA)为四边形abcd的角点a，Xs,Ys,Zs为影像外方位元素中的三个角元素，a1,b1,c1,a2,b2,c2,a3,b3,c3为影像定向结果中的旋转矩阵，f为相机焦距，ZA为点(XA，YA)在数字高程模型DEM上内插的高程值。
之后，重复以上步骤，得到测区所有影像有效区域对应的局部DEM范围和原始影像局部纹理范围。
步骤3：CPU端计算镶嵌后DOM全景影像图的大小，并在硬盘上创建镶嵌后DOM全景影像存储指针，即创建影像拼接结果指针，其中DOM全景影像图大小的计算步骤如下：
步骤3-1：将所有原始影像的四个角点投影到地面，得到每张影像纠正后四个方向的极值Ximin，Ximax，Yimin，Yimax，按照公式求得镶嵌后DOM影像在四个方向上的极值Xmin,Xmax,Ymin,Ymax，其中，min()为该一维数组中的最小值，max()为该一维数组中的最大值，i＝1～n，n为影像张数；
步骤3-2：按照计算镶嵌后DOM全景图的宽度width和高度height，M为纠正后正射影像的地面分辨率；
步骤3-3：在计算后将每张影像对应的有效纹理和局部DEM读入显存。
步骤4：根据步骤2获得的局部DEM和原始影像局部纹理，以及影像对应的外参数，GPU端依次对测区每张影像的有效区域进行正射纠正并写入DOM全景影像结果中，如图5所示，所述的正射纠正的具体步骤如下：
步骤4-1：先将当前的第i影像有效区域对应的局部DEM、原始影像局部纹理及当前影像对应的外参数从内存拷入显存；
步骤4-2：根据当前影像有效区域多边形PolygonS的四个极值Xmin,Xmax,Ymin,Ymax和纠正后正射影像的地面分辨率M，按照公式按照求得该影像有效多边形PolygonS最小外接矩形ABCD的宽度width和高度height，width*height则是该块需要正射纠正像素的个数；
接着，即可根据计算任务确定GPU线程块的大小和数目，本例中，假定线程格网设置为二维，线程块的大小设为dimBlock(N,N)，即每个线程块中有N*N个线程，格网dimGrid的大小可以根据公式计算得到；
步骤4-3：GPU多线程并行，即所有线程同时执行，利用反解法数字微分纠正对该影像块进行正射纠正，按照计算公式利用线程索引计算该像素对应的地面点的坐标(X,Y)，其中M为纠正后正射影像地面几何分辨率，(X0,Y0)为矩形块A'B'C'D'左上角的坐标，Height为矩形块A'B'C'D'的行数，(row,col)为当前像素在矩形块A'B'C'D'中的行列位置；
步骤4-4：判断点(X,Y)是否在当前影像有效区域多边形PolygonS内部，如果在内部则进行正射纠正，如果在外部则不进行处理。因为多边形PolygonS是凸多边形，因此可以通过判断点(X,Y)在多边形PolygonS每条边的左右两侧的方法来判断点是否在其内部。具体方法如下：设多边形PolygonS有n个顶点P1,P2,P3,···,Pn，n条边分别为P1P2,P2P3,P3P4,···,PnP1，依次判断点(X,Y)在边P1P2,P2P3,P3P4,···,PnP1的左右侧，如果全部在左侧或是右侧则点(X,Y)在多边形PolygonS内部，如果有一条边与其他边不一样，则点在多边形PolygonS外部，如果点(X,Y)在多边形PolygonS的任意边上，则视为其在内部。
如果在有效多边形内部，根据点坐标(X,Y)在局部DEM块中双线性内插得到该点的高程Z，再根据公式计算该像素在原始影像块上的像点坐标p(x,y)；
步骤4-5：线程按索引号进行灰度赋值，即将双线性内插三个波段上灰度值作为纠正后影像上像素(col,row)的灰度值；
步骤4-6：所有线程执行完毕后，将所计算得到的数据从显存拷贝到内存，释放开辟的所有显存；CPU将获得到的纠正后影像数据从内存写到硬盘上创建好的DOM全景图的对应位置，释放所有已开辟内存。
在步骤4的进行过程中，重复上述的步骤4-1至步骤4-6，直至完成测区所有影像有效区域多边形的正射纠正和写出。
步骤5：对于正射纠正后的影像，GPU端依次对每条镶嵌线两侧的影像进行线性加权融合处理，并将融合处理后的纹理拷出内存，按照步骤3创建的DOM全景影像存储指针写到对应的位置，对原始镶嵌结果进行更新。
根据步骤1中计算的镶嵌线网络可知，每条镶嵌线对应两张影像，镶嵌线两侧影像的融合过程只与镶嵌线两侧一定范围内的局部影像有关，因此为了提高GPU-CPU协同镶嵌线两侧融合计算效率，不能把两侧影像的所有纹理数据都读入显存，只需将有用的局部纹理读入即可。如附图6所示，L0L1为影像L和R之间的镶嵌线，在L0L1两侧分别选择宽度为L的缓冲区域，计算包含该区域的最小外接矩形RECline，然后按照步骤2的方法计算得到当前镶嵌线对应的局部DEM和两侧影像上局部的原始纹理数据，将其与对应的外方位元素等参数从内存拷入显存。
而镶嵌线两侧影像线性加权融合，是指镶嵌影像在缓冲区左侧的像素取左影像L的像素灰度值，在缓冲区右侧的像素取右影像R上的像素灰度值，缓冲区内部的像素则取两影像对应像素灰度的距离加权值。
如图7所示，对第j条镶嵌线两侧影像线性加权融合处理的步骤具体为：
步骤5-1：计算第j条镶嵌线对应的局部DEM以及该镶嵌线两侧影像上局部的原始纹理数据，并将其从内存拷入显存，其中局部DEM的计算方法如步骤2所述；
步骤5-2：根据矩形RECline的大小来确定设置线程格网中线程块的大小和数量，设置方法与步骤4相似，在此不做赘述；
步骤5-3：每一个线程对应一个像素，线程格网中所有线程并行计算，对于每一个线程，根据线程索引计算像素(row,col)对应的物方坐标(X,Y)与高程Z；
步骤5-4：通过计算点(X,Y)与镶嵌线L0L1的距离大小来判断该点是否在镶嵌线两侧的缓冲区域内部，如果是则根据公式分别计算出该点在两侧影像中位置，并线程按索引号内插计算出点(X,Y)在两原始影像上像点灰度值g1、g2，；
步骤5-5：线程按索引号对灰度值g1、g2进行加权内插得到最终灰度值G，即按照公式进行距离加权得到像素(row,col)的最终灰度值G，式中L为镶嵌线一侧缓冲区大小，d为当前像素与镶嵌线的距离，g1为该像素在影像L上的灰度值，g2为当前像素在影像R上的灰度值；
步骤5-6：线程按索引进行灰度赋值，即将双线性内插的灰度值作为纠正后影像上像素(col,row)的灰度值；
步骤5-7：CPU将融合后的局部影像覆盖原来影像，对当前镶嵌线两侧的影像纹理进行更新，经过镶嵌线两侧影像纹理线性加权融合前后对比如图10所示，即：待所有线程计算完毕后，将线性加权融合后的矩形RECline的影像纹理从显存拷出到内存，再由内存写到镶嵌好的DOM全景图对应位置，更新当前镶嵌线两侧的影像纹理，并释放所开辟的显存和内存。
在步骤5的运行过程中，重复上述的步骤5-1至步骤5-7，直至完成所有镶嵌线两侧影像纹理的线性加权融合处理，镶嵌后的DOM全景影像图如图8与图9所示，其中图9不含镶嵌线。
本方案通过计算测区所有影像像底点位置并自动生成Voronoi图镶嵌线网络；计算每张影像的有效区域多边形对应的局部DEM和局部原始影像纹理范围；计算镶嵌后测区DOM全景影像图的大小并创建结果影像指针；完成每张影像有效区域多边形的正射纠正并写入结果影像指针；完成每条镶嵌线两侧影像的快速融合后将计算结果拷出并更新原DOM全景影像图等技术手段，突破了先单张全幅正射纠正，然后进行镶嵌线网络选择，再进行镶嵌的传统模式，从匀光后原始影像出发先计算镶嵌线网络及每张影像的镶嵌有效区域，再对每张影像的有效区域进行正射纠正，直接镶嵌输出测区DOM全景图。因此，本方法没有中间数据产生，并且只对每张影像的有效区域进行正射纠正，大大减少了数据冗余以及时间和资源的开销。
- 还没有人留言评论。精彩留言会获得点赞！