图像重采样（CPU和GPU）
Posted 机器学习猪
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了图像重采样（CPU和GPU）相关的知识，希望对你有一定的参考价值。
1 前言
之前在写影像融合算法的时候，免不了要实现将多光谱影像重采样到全色大小。当时为了不影响融合算法整体开发进度，其中重采样功能用的是GDAL开源库中的Warp接口实现的。
后来发现GDAL Warp接口实现的多光谱到全色影像的重采样主要存在两个问题：1 与原有平台的已有功能不兼容，产生冲突；2 效率较低。因此，决定重新设计和开发一个这样的功能，方便后期软件系统的维护等。
2 图像重采样
图像处理从形式上来说主要包括两个方面：1 单像素或者邻域像素的处理，比如影像的相加或者滤波运算等；2 图像几何空间变换，如图像的重采样，配准等。
影像重采样的几何空间变换公式如下：
其中 为变换系数，常用的重采样算法主要包括以下三种：1 最近邻；2 双线性；3 三次卷积。
2.1 最近邻采样
最近邻采样的原理概况起来就是用采样点位置最近的一个像素值替代采样点位置的像素值。在这里插入一点：
通常图像空间变换有两种方法，直接法或者间接法。以图像重采样为例说明如下：直接法：从原始的图像行列初始值开始，根据变换公式，计算采样后的像素位置，并对位置赋值，但是这种方法会出现，原始图像的多个像素点对应到同一采样后的像素点，从而还要增加额外方法进行处理；间接法：是从重采样后图像的行列初始值开始，计算得到其在原始影像中的位置，并根据一定的算法进行计算，得到采样后的值。这种方法简单直接，本文就是采用这样的方法。
最近邻采样的实现算法如下：
首先遍历采样点，根据公式计算采样点在原始图像中的位置，假设位置为 。然后计算与 最近的点，并将其像素值赋为采样点的像素值。其公式如下：
2.2 双线性
双线性采样和最近邻赋值不同，是找到 最近的四个像素点，然后将距离作为权重分别插值 和 方向，从而插值到采样点位置。具体公式见代码部分。
2.3 三次卷积
三次卷积是利用 最近的16个像素点进行插值计算得到。同样也是分别插值 和 方向。具体公式见代码部分。
3 实验结果
主要实现了两个版本的差值结果：1 CPU版本；2 GPU版本（初学）。考虑到多光谱和全色影像范围大小不一致的事实，算法首先计算全色和多光谱的重叠区域。话不多说，先看看详细的代码实现过程。
CPU版本：
1 #ifndef RESAMPLECPU_H 2 #define RESAMPLECPU_H 3 4 #include <gdal_alg_priv.h> 5 #include <gdal_priv.h> 6 #include <assert.h> 7 8 9 template<typename T> 10 void ReSampleCPUKernel(const float *mssData, 11 T *resampleData, 12 int mssWidth, 13 int mssHeight, 14 int mssBandCount, 15 int mssOffsetX, 16 int mssOffsetY, 17 int panWidth, 18 int panHeight, 19 float radioX, 20 float radioY, 21 double dfDstNoDataValue, 22 int MethodType) 23 { 24 assert(mssData != NULL); 25 float eps = 0.00001f; 26 for(int idx = 0;idx < panHeight;idx++){ 27 for(int idy = 0;idy<panWidth;idy++){ 28 // 找到对应的MSS像素位置 29 float curX = (float)idx * radioX; 30 float curY = (float)idy * radioY; 31 if(mssData[int(curX)*mssWidth*mssBandCount + int(curY)] == dfDstNoDataValue) 32 { 33 int i = 0; 34 while(i < mssBandCount){ 35 resampleData[idx*panWidth*mssBandCount+i*panWidth + idy] = T(dfDstNoDataValue); 36 i++; 37 } 38 continue; 39 } 40 if(MethodType == 0){ // 最近邻 41 int nearX = (int)(curX + 0.5)>(int)curX?(int)(curX + 1):(int)curX; 42 int nearY = (int)(curY + 0.5)>(int)curY?(int)(curY + 1):(int)curY; 43 if(nearX >= mssHeight - 1){ 44 nearX = mssHeight - 1; 45 } 46 if(nearY >= mssWidth - 1){ 47 nearY = mssWidth - 1; 48 } 49 if(nearX < mssHeight && nearY < mssWidth){ 50 int i = 0; 51 while(i < mssBandCount){ 52 float tmp = mssData[nearX*mssWidth*mssBandCount + i*mssWidth + nearY]; 53 resampleData[idx*panWidth*mssBandCount+i*panWidth + idy] = T(tmp); 54 i++; 55 } 56 } 57 } 58 59 if(MethodType == 1){ // 双线性 60 float dataX = curX - (int)curX; 61 float dataY = curY - (int)curY; 62 int preX = (int)curX; 63 int preY = (int)curY; 64 int postX = (int)curX + 1; 65 int postY = (int)curY + 1; 66 if(postX >= mssHeight - 1){ 67 postX = mssHeight - 1; 68 } 69 if(postY >= mssWidth - 1){ 70 postY = mssWidth - 1; 71 } 72 73 float Wx1 = 1 - dataX; 74 float Wx2 = dataX; 75 float Wy1 = 1 - dataY; 76 float Wy2 = dataY; 77 // 双线性差值核心代码 78 int i = 0; 79 while(i < mssBandCount){ 80 float pMssValue[4] = {0,0,0,0}; 81 pMssValue[0] = mssData[preX*mssWidth*mssBandCount + i*mssWidth + preY]; 82 pMssValue[1] = mssData[preX*mssWidth*mssBandCount + i*mssWidth + postY]; 83 pMssValue[2] = mssData[postX*mssWidth*mssBandCount + i*mssWidth + preY]; 84 pMssValue[3] = mssData[postX*mssWidth*mssBandCount + i*mssWidth + postY]; 85 float tmp = Wy1*(Wx1*pMssValue[0] + Wx2*pMssValue[2]) + Wy2*(Wx1*pMssValue[1] + Wx2*pMssValue[3]); 86 resampleData[idx*panWidth*mssBandCount+i*panWidth + idy] = T(tmp); 87 i++; 88 } 89 } 90 91 92 if(MethodType == 2){ // 双三次卷积 93 float dataX = curX - (int)curX; 94 float dataY = curY - (int)curY; 95 int preX1 = (int)curX - 1; 96 int preX2 = (int)curX; 97 int postX1 = (int)curX + 1; 98 int postX2 = (int)curX + 2; 99 int preY1 = (int)curY - 1; 100 int preY2 = (int)curY; 101 int postY1 = (int)curY + 1; 102 int postY2 = (int)curY + 2; 103 if(preX1 < 0) preX1 = 0; 104 if(preY1 < 0) preY1 = 0; 105 if(postX1 > mssHeight - 1) postX1 = mssHeight - 1; 106 if(postX2 > mssHeight - 1) postX2 = mssHeight - 1; 107 if(postY1 > mssWidth - 1) postY1 = mssWidth - 1; 108 if(postY2 > mssWidth - 1) postY2 = mssWidth - 1; 109 110 float Wx1 = -1.0f*dataX + 2*dataX*dataX - dataX*dataX*dataX; 111 float Wx2 = 1 - 2*dataX*dataX + dataX*dataX*dataX; 112 float Wx3 = dataX + dataX*dataX - dataX*dataX*dataX; 113 float Wx4 = -1.0f*dataX*dataX + dataX*dataX*dataX; 114 float Wy1 = -1.0f*dataY + 2*dataY*dataY - dataY*dataY*dataY; 115 float Wy2 = 1 - 2*dataY*dataY + dataY*dataY*dataY; 116 float Wy3 = dataY + dataY*dataY - dataY*dataY*dataY; 117 float Wy4 = -1.0f*dataY*dataY + dataY*dataY*dataY; 118 119 int i = 0; 120 while(i < mssBandCount){ 121 float *pMssValue = new float[16]; 122 memset(pMssValue,0,sizeof(float)*16); 123 pMssValue[0] = mssData[preX1*mssWidth*mssBandCount + i*mssWidth + preY1]; 124 pMssValue[1] = mssData[preX1*mssWidth*mssBandCount + i*mssWidth + preY2]; 125 pMssValue[2] = mssData[preX1*mssWidth*mssBandCount + i*mssWidth + postY1]; 126 pMssValue[3] = mssData[preX1*mssWidth*mssBandCount + i*mssWidth + postY2]; 127 128 pMssValue[4] = mssData[preX2*mssWidth*mssBandCount + i*mssWidth + preY1]; 129 pMssValue[5] = mssData[preX2*mssWidth*mssBandCount + i*mssWidth + preY2]; 130 pMssValue[6] = mssData[preX2*mssWidth*mssBandCount + i*mssWidth + postY1]; 131 pMssValue[7] = mssData[preX2*mssWidth*mssBandCount + i*mssWidth + postY2]; 132 133 pMssValue[8] = mssData[postX1*mssWidth*mssBandCount + i*mssWidth + preY1]; 134 pMssValue[9] = mssData[postX1*mssWidth*mssBandCount + i*mssWidth + preY2]; 135 pMssValue[10] = mssData[postX1*mssWidth*mssBandCount + i*mssWidth + postY1]; 136 pMssValue[11] = mssData[postX1*mssWidth*mssBandCount + i*mssWidth + postY2]; 137 138 pMssValue[12] = mssData[postX2*mssWidth*mssBandCount + i*mssWidth + preY1]; 139 pMssValue[13] = mssData[postX2*mssWidth*mssBandCount + i*mssWidth + preY2]; 140 pMssValue[14] = mssData[postX2*mssWidth*mssBandCount + i*mssWidth + postY1]; 141 pMssValue[15] = mssData[postX2*mssWidth*mssBandCount + i*mssWidth + postY2]; 142 143 float tmp = Wy1*(Wx1*pMssValue[0] + Wx2*pMssValue[4] + Wx3*pMssValue[8] + Wx4*pMssValue[12])+ 144 Wy2*(Wx1*pMssValue[1] + Wx2*pMssValue[5] + Wx3*pMssValue[9] + Wx4*pMssValue[13])+ 145 Wy3*(Wx1*pMssValue[2] + Wx2*pMssValue[6] + Wx3*pMssValue[10] + Wx4*pMssValue[14])+ 146 Wy4*(Wx1*pMssValue[3] + Wx2*pMssValue[7] + Wx3*pMssValue[11] + Wx4*pMssValue[15]); 147 resampleData[idx*panWidth*mssBandCount+i*panWidth + idy] = T(tmp); 148 delete []pMssValue;pMssValue = NULL; 149 i++; 150 } 151 } 152 } 153 } 154 } 155 156 int ReSampleCPUApp(const char *mssfileName, 157 const char *panfileName, 158 const char *resamplefileName, 159 int MethodType = 1) 160 { 161 GDALAllRegister(); 162 GDALDataset *poPANDS = (GDALDataset*)GDALOpen(panfileName,GA_ReadOnly); 163 GDALDataset *poMSSDS = (GDALDataset*)GDALOpen(mssfileName,GA_ReadOnly); 164 if(!poPANDS || !poMSSDS) 165 return -1; 166 167 //MSS info 168 int mssBandCount = poMSSDS->GetRasterCount(); 169 int mssWidth = poMSSDS->GetRasterXSize(); 170 int mssHeight = poMSSDS->GetRasterYSize(); 171 double adfMssGeoTransform[6] = {0}; 172 poMSSDS->GetGeoTransform(adfMssGeoTransform); 173 GDALDataType mssDT = poMSSDS->GetRasterBand(1)->GetRasterDataType(); 174 175 int bSrcHasNoData; 176 double dfSrcNoDataValue = 0; 177 dfSrcNoDataValue = GDALGetRasterNoDataValue(poMSSDS->GetRasterBand(1),&bSrcHasNoData); 178 if(!bSrcHasNoData) dfSrcNoDataValue = 0.0; 179 //DT = mssDT; 180 181 int *pBandMap= new int[mssBandCount]; 182 for(int i = 0;i<mssBandCount;i++){ 183 pBandMap[i] = i+1; 184 } 185 186 // PAN Info 187 int panBandCount = poPANDS->GetRasterCount(); 188 int panWidth = poPANDS->GetRasterXSize(); 189 int panHeidht = poPANDS->GetRasterYSize(); 190 double adfPanGeoTransform[6] = {0}; 191 poPANDS->GetGeoTransform(adfPanGeoTransform); 192 GDALDataType panDT = poPANDS->GetRasterBand(1)->GetRasterDataType(); 193 194 // 创建新数据集=======投影信息 195 double adfResampleGeoTransform[6] = {0}; 196 adfResampleGeoTransform[1] = adfPanGeoTransform[1]; 197 adfResampleGeoTransform[5] = adfPanGeoTransform[5]; 198 adfResampleGeoTransform[2] = adfPanGeoTransform[2]; 199 adfResampleGeoTransform[4] = adfPanGeoTransform[4]; 200 if(adfMssGeoTransform[0] >= adfPanGeoTransform[0]){ 201 adfResampleGeoTransform[0] = adfMssGeoTransform[0]; 202 }else{ 203 adfResampleGeoTransform[0] = adfPanGeoTransform[0]; 204 } 205 if(adfMssGeoTransform[3] > adfPanGeoTransform[3]){ 206 adfResampleGeoTransform[3] = adfPanGeoTransform[3]; 207 }else{ 208 adfResampleGeoTransform[3] = adfMssGeoTransform[3]; 209 } 210 211 // 创建新数据集=======影像大小 212 double panEndX = adfPanGeoTransform[0] + panWidth*adfPanGeoTransform[1] + 213 panHeidht*adfPanGeoTransform[2]; 214 double panEndY = adfPanGeoTransform[3] + panHeidht*adfPanGeoTransform[4] + 215 panHeidht*adfPanGeoTransform[5]; 216 217 double mssEndX = adfMssGeoTransform[0] +mssWidth*adfMssGeoTransform[1] + 218 mssHeight*adfMssGeoTransform[2]; 219 double mssEndY = adfMssGeoTransform[3] + mssWidth*adfMssGeoTransform[4] + 220 mssHeight*adfMssGeoTransform[5]; 221 double resampleEndXY[2] = {0}; 222 if(panEndX > mssEndX) 223 resampleEndXY[0] = mssEndX; 224 else 225 resampleEndXY[0] = panEndX; 226 if(panEndY >= mssEndY) 227 resampleEndXY[1] = panEndY; 228 else 229 resampleEndXY[1] = mssEndY; 230 231 // 创建新数据集=======MSS AND PAN 有效长宽 232 int resampleWidth = static_cast<int>((resampleEndXY[0] - adfResampleGeoTransform[0])/adfResampleGeoTransform[1] + 0.5); 233 int resampleHeight = static_cast<int>((resampleEndXY[1] - adfResampleGeoTransform[3])/adfResampleGeoTransform[5] + 0.5); 234 int mssEffectiveWidth = static_cast<int>((resampleEndXY[0] - adfResampleGeoTransform[0])/adfMssGeoTransform[1] + 0.5); 235 int mssEffectiveHeight = static_cast<int>((resampleEndXY[1] - adfResampleGeoTransform[3])/adfMssGeoTransform[5] + 0.5); 236 int panEffectiveWidth = resampleWidth; 237 int panEffectiveHeight = resampleHeight; 238 239 // 创建新数据集=======位置增益大小 240 int mssGainX = static_cast<int>((adfResampleGeoTransform[0] - adfMssGeoTransform[0])/adfMssGeoTransform[1] + 0.5); 241 int mssGainY = static_cast<int>((adfResampleGeoTransform[3] - adfMssGeoTransform[3])/adfMssGeoTransform[5] + 0.5); 242 int panGainX = static_cast<int>((adfResampleGeoTransform[0] - adfPanGeoTransform[0])/adfPanGeoTransform[1] + 0.5); 243 int panGainY = static_cast<int>((adfResampleGeoTransform[3] - adfPanGeoTransform[3])/adfPanGeoTransform[5] + 0.5); 244 245 246 // 创建新数据集=======创建文件 247 GDALDriver *poOutDriver = (GDALDriver*)GDALGetDriverByName("GTIFF"); 248 if(!poOutDriver){ 249 return -1; 250 } 251 GDALDataset *poOutDS = poOutDriver->Create(resamplefileName,resampleWidth, 252 resampleHeight,mssBandCount,mssDT,NULL); 253 poOutDS->SetGeoTransform(adfResampleGeoTransform); 254 poOutDS->SetProjection(poPANDS->GetProjectionRef()); 255 256 257 // 重采样核心代码============图像分块 258 int iNumRow = 256; 259 if(iNumRow > resampleHeight){ 260 iNumRow = 1; 261 } 262 int loopNum = (resampleHeight + iNumRow - 1)/iNumRow; //分块数 263 int nLineSpace,nPixSpace,nBandSpace; 264 nLineSpace = sizeof(float)*mssEffectiveWidth*mssBandCount; 265 nPixSpace = 0; 266 nBandSpace = sizeof(float)*mssEffectiveWidth; 267 268 // 重采样采样比例 269 float radioX = float(adfPanGeoTransform[1]/adfMssGeoTransform[1]); 270 float radioY = float(adfPanGeoTransform[5]/adfMssGeoTransform[5]); 271 272 int mssCurPosX = mssGainX; 273 int mssCurPosY = mssGainY; 274 int mssCurWidth = 0; 275 int mssCurHeight = 0; 276 277 // 重采样核心代码============ 278 for(int i = 0;i<loopNum;i++){ 279 int tmpRowNum = iNumRow; 280 int startR = i*iNumRow; 281 int endR = startR + iNumRow - 1; 282 if(endR>resampleHeight -1){ 283 tmpRowNum = resampleHeight - startR; 284 //endR = startR + tmpRowNum - 1; 285 } 286 //计算读取的MSS影像区域大小 287 int mssCurWidth = mssEffectiveWidth; 288 int mssCurHeight = 0; 289 if(MethodType == 0) 290 mssCurHeight = int(tmpRowNum*radioY); 291 else if(MethodType == 1) 292 mssCurHeight = int(tmpRowNum*radioY)+1; 293 else if(MethodType == 2) 294 mssCurHeight = int(tmpRowNum*radioY)+2; 295 296 if(mssCurHeight + mssCurPosY > mssHeight - 1){ 297 mssCurHeight = mssHeight - mssCurPosY; 298 } 299 300 //创建数据 301 /*float *resampleBuf = (float *)malloc(sizeof(cl_float)*tmpRowNum*resampleWidth*trueBandCount);*/ 在 iPhone 上处理 GPU（金属）和 CPU（OpenCV）上的摄像头馈送数据