基于OpenCL的数字地形分析之坡度坡向提取
2024-05-09 22:38:50
转自：http://blog.csdn.net/zhouxuguang236/article/details/23568875
又有一段时间没有发表博客了，可能最近工作有点忙。今天就把最近的学习和研究成果和大家分享一下。对于GIS稍微有点了解的人都知道地形分析中的坡度和坡向，这是数字地形分析中最基本的分析了，对于数字地形分析中很多计算都是邻域分析，所以非常适合数据并行。
一、相关概念和公式
坡度严格地讲，是地表任意一点过该点的切平面与水平面之间的夹角。坡度表示了地表的倾斜程度。坡度的计算公式如下：
其中fx和fy分别代表x和y方向上的高程变化率。坡度有两种表示方法，一种是坡度角，另外一种是坡度百分比，即高程增量与水平增量的百分比。
两种表示方法如下：
图：坡度的两种表示方式
而坡向的定义是地面任意一点切平面法线在水平面上的投影与正北方向的方向角。坡向的计算公式如下：
由上面两个公式可以看出，都需要求fx和fy。对于这两个分量的计算方法有很多种。具体的计算方法有二阶差分、三阶差分等，更多的计算方法见下图，图中Slope-we对应fx，Slope-sn对应fy。
各个商业软件基本上都是选取上述的算法进行计算，ArcGIS采用的算法二，ERDAS采用的是算法4。本文采用算法作为计算方法。
二、CPU上的实现
坡度和坡向算法只是最后的计算方法不一样，前面的计算过程都差不多。所以本文就主要讲解坡度的提取。根据前面的讲解，计算要保证高程的单位和水平面上的单位要保持一样。所以，有两种方法解决这个问题，一个是让用户在界面去输入Z方向和高程和水平面方向单位的换算系数；另外一种方法是程序自动换算。我这里就采用程序自动换算，数据的IO操作采用GDAL，其具体的方法如下：
- //投影
- const char pszWkt = poSrcDS->GetProjectionRef();
- double dbNres = 0;
- double dbEres = 0;
- double dfGeotrans[6];
- if (poSrcDS->GetGeoTransform(dfGeotrans) != CE_None)
- {
- dbNres = fabs(dfGeotrans[5]);
- dbEres = fabs(dfGeotrans[1]);
- }
- else
- {
- OGRSpatialReference pSrs = (OGRSpatialReference*)OSRNewSpatialReference(pszWkt);
- if (pSrs != NULL)
- {
- if (pSrs->IsGeographic())
- {
- dbNres = fabs(dfGeotrans[5]);
- dbNres = 110000;
- dbEres = fabs(dfGeotrans[1]);
- dbEres = 110000;
- }
- else if (pSrs->IsGeocentric())
- {
- dbNres = fabs(dfGeotrans[5]);
- dbEres = fabs(dfGeotrans[1]);
- }
- else if (pSrs->IsProjected())
- {
- dbNres = fabs(dfGeotrans[5]);
- dbEres = fabs(dfGeotrans[1]);
- }
- }
- OSRDestroySpatialReference((OGRSpatialReferenceH)pSrs);
- }
- SlopeOption slopeOption;
- slopeOption.dbEwres = dbEres;
- slopeOption.dbNsres = dbNres;
- slopeOption.slopeType = eSlopeType;
//投影 const char* pszWkt = poSrcDS->GetProjectionRef();double dbNres = 0;
double dbEres = 0;double dfGeotrans[6];
if (poSrcDS->GetGeoTransform(dfGeotrans) != CE_None)
{dbNres = fabs(dfGeotrans[5]);dbEres = fabs(dfGeotrans[1]);
}else
{OGRSpatialReference* pSrs = (OGRSpatialReference*)OSRNewSpatialReference(pszWkt);if (pSrs != NULL){if (pSrs->IsGeographic()){dbNres = fabs(dfGeotrans[5]);dbNres *= 110000;dbEres = fabs(dfGeotrans[1]);dbEres *= 110000;}else if (pSrs->IsGeocentric()){dbNres = fabs(dfGeotrans[5]);dbEres = fabs(dfGeotrans[1]);}else if (pSrs->IsProjected()){dbNres = fabs(dfGeotrans[5]);dbEres = fabs(dfGeotrans[1]);}}OSRDestroySpatialReference((OGRSpatialReferenceH)pSrs);
}SlopeOption slopeOption;
slopeOption.dbEwres = dbEres;
slopeOption.dbNsres = dbNres;
slopeOption.slopeType = eSlopeType;</pre><br><br><br><br>此外，还需要包装两个结构用于传参数；<br>//坡度的表达方式<br>typedef enum <br>{<br><span style="white-space:pre;"> </span>DEGREE_SLOPE,<span style="white-space:pre;"> </span>//度数方式<br><span style="white-space:pre;"> </span>PERCENT_SLOPE<span style="white-space:pre;"> </span>//百分比方式<br>}SLOPE_TYPE;<br><br><br>//坡度算法结构体<br>struct SlopeOption<br>{<br><span style="white-space:pre;"> </span>double dbNsres;<span style="white-space:pre;"> </span>//南北方向分辨率<br><span style="white-space:pre;"> </span>double dbEwres;<span style="white-space:pre;"> </span>//东西方向分辨率<br><span style="white-space:pre;"> </span>SLOPE_TYPE slopeType;<span style="white-space:pre;"> </span>//坡度方式<br>} ;<br>最终，我设计的计算坡度的函数如下：<br>bool ExtractSlope(const char* pszDEMfile,const char* pszOutSlpoeFile,const char* pszFormat, SLOPE_TYPE eSlopeType ,double dbScale)<br>pszDEMfile为输入DEM数据；pszOutSlpoeFile为输出坡度数据；pszFormat为输出影像格式；<br>eSlopeType是坡度的类型，即坡度百分比还是坡度角度；dbScale可选，一般情况下为1.0。<br><br><br> 在GIS和遥感领域，一般情况下，栅格数据都是非常大的，不可能将所有像素全部读进来一次性处理，所以就必然涉及到影像分块。关于影像分块，对于不同的算法可以采用不同的策略，其主要策略如下：<br> 对于单点运算，即运算过程之和当前计算的像素相关的算法，分块可以按照按照行来分块，也可以矩形分块，类似于地图切片。一般来说，影像是按照行优先的顺序类存储的，所以按照行分块可以减少文件指针移动的次数，提高速度。<br> 对于邻域相关的运算，如本文提到的坡度、卷积等运算，也可以按照行分块或者矩形分块。但是要注意一个问题，这样简单的分块会导致块与块之间的结果不连续。针对这种情况，为了缝合块边缘的结果，可以再读取数据的时候边缘的数据重复读取，这样就保证了最后合成的结果具有连续性，也就是说保证有一定的重叠度，一般重叠度为邻域窗口大小的一半。<br> 还有一种情况是全局相关的运算，如求图像的最值，这个也可以采用上述两种方法都可以，即分好块之后，最后求得各个块的最值。<br> 总结一下，分块的思想体现了算法设计中的分治法思想，即所谓的分而治之。<br> 对于本文，我采用的分块读取、分块处理和分块写入数据的思路如下：<br>1、按照行来分块，窗口大小为3，所以块之间重叠一个像素；<br>2、对于分块后只有一个块的情况，不需要做特别处理<br>3、对于分块后有多个块的情况，也根据块的索引做不同处理：为了说明方便，我先定义几个变量<br><span style="white-space:pre;"> </span>//实际的块的大小<br><span style="white-space:pre;"> </span>int nRealWidth = nXsize;<br><span style="white-space:pre;"> </span>int nRealHeight = nSubHeight;<br><span style="white-space:pre;"> </span>//读取数据块的大小<br><span style="white-space:pre;"> </span>int nReadWidth = nXsize;<br><span style="white-space:pre;"> </span>int nReadHeight = nSubHeight;<br><span style="white-space:pre;"> </span>int nYOffset = i*nSubHeight;<span style="white-space:pre;"> </span>//某一块读取数据的Y方向上的偏移量<br>nSubHeight为分块的高<br> a、对于第0个块，块的实际读取数据的高度nReadHeight为nReadHeight += 1;，即向下要多读取一行数据；写入数据的时候，实际写入的高度为nRealHeight，Y方向上的偏移量为0。<br> b、对于介于0和最后一个块之间的块，块的最顶部要和上一块重叠一行像素，块的最下部要和相邻块重叠一行像素，即多读取两行数据，nReadHeight += 2;nYOffset -= 1; 写入数据的时候，实际写入的高度为nRealHeight= nReadHeight - 2;，Y方向上的偏移量为nYOffset = i*nSubHeight;<br> c、对于最后一个块，块的实际读取数据的高度nReadHeight为nReadHeight += 1，即向上要多读取一行的数据，nReadHeight = nRealHeight + 1;nYOffset -= 1; 实际写入的高度为nRealHeight= nReadHeight - 1;，Y方向上的偏移量为nYOffset = i*nSubHeight;<br><br><br>这样，就保证了分块之间结果缝合在一起了。<br><p>为了方便直观表达其意思，其代码如下：</p><p></p><div class="dp-highlighter bg_cpp"><div class="bar"><div class="tools"><b>[cpp]</b> <a title="view plain" class="ViewSource" onclick="dp.sh.Toolbar.Command('ViewSource',this);return false;" href="#">view plain</a><span class="tracking-ad" data-mod="popu_168"> <a title="copy" class="CopyToClipboard" onclick="dp.sh.Toolbar.Command('CopyToClipboard',this);return false;" href="#">copy</a><div style="left: 243px; top: 4693px; width: 16px; height: 16px; position: absolute; z-index: 99;"><embed name="ZeroClipboardMovie_2" width="16" height="16" align="middle" id="ZeroClipboardMovie_2" pluginspage="http://www.macromedia.com/go/getflashplayer" src="http://csdnimg.cn/public/highlighter/ZeroClipboard.swf" type="application/x-shockwave-flash" wmode="transparent" flashvars="id=2&width=16&height=16" allowfullscreen="false" allowscriptaccess="always" bgcolor="#ffffff" quality="best" menu="false" loop="false"></div></span><span class="tracking-ad" data-mod="popu_169"> <a title="print" class="PrintSource" onclick="dp.sh.Toolbar.Command('PrintSource',this);return false;" href="#">print</a></span><a title="?" class="About" onclick="dp.sh.Toolbar.Command('About',this);return false;" href="#">?</a></div></div><ol class="dp-cpp"><li class="alt"><span><span class="comment">//分块处理</span><span> </span></span></li><li><span> </span></li><li class="alt"><span> <span class="datatypes">int</span><span> nSubHeight = 2000; </span></span></li><li><span> <span class="datatypes">int</span><span> nYBlockCount = (nYsize+nSubHeight-1)/nSubHeight; </span><span class="comment">//计算分块的数量</span><span> </span></span></li><li class="alt"><span> <span class="keyword">for</span><span> (</span><span class="datatypes">int</span><span> i = 0; i < nYBlockCount; i ++) </span></span></li><li><span> { </span></li><li class="alt"><span> <span class="comment">//实际的块的大小</span><span> </span></span></li><li><span> <span class="datatypes">int</span><span> nRealWidth = nXsize; </span></span></li><li class="alt"><span> <span class="datatypes">int</span><span> nRealHeight = nSubHeight; </span></span></li><li><span> <span class="comment">//读取数据块的大小</span><span> </span></span></li><li class="alt"><span> <span class="datatypes">int</span><span> nReadWidth = nXsize; </span></span></li><li><span> <span class="datatypes">int</span><span> nReadHeight = nSubHeight; </span></span></li><li class="alt"><span> <span class="datatypes">int</span><span> nYOffset = i*nSubHeight; </span></span></li><li><span> <span class="keyword">if</span><span> (1 == nYBlockCount) </span></span></li><li class="alt"><span> { </span></li><li><span> nRealHeight = nYsize; </span></li><li class="alt"><span> nReadHeight = nRealHeight; </span></li><li><span> } </span></li><li class="alt"><span> <span class="keyword">else</span><span> </span></span></li><li><span> { </span></li><li class="alt"><span> <span class="keyword">if</span><span> (i == 0) </span></span></li><li><span> { </span></li><li class="alt"><span> nReadHeight += 1; </span></li><li><span> } </span></li><li class="alt"><span> <span class="keyword">else</span><span> </span><span class="keyword">if</span><span> (i > 0 && i < nYBlockCount-1) </span></span></li><li><span> { </span></li><li class="alt"><span> nReadHeight += 2; </span></li><li><span> nYOffset -= 1; </span></li><li class="alt"><span> } </span></li><li><span> <span class="keyword">else</span><span> </span><span class="keyword">if</span><span>(i == nYBlockCount-1) </span></span></li><li class="alt"><span> { </span></li><li><span> nRealHeight = nYsize-nSubHeight*(nYBlockCount-1); </span></li><li class="alt"><span> nReadHeight = nRealHeight + 1; </span></li><li><span> nYOffset -= 1; </span></li><li class="alt"><span> } </span></li><li><span> } </span></li><li class="alt"><span> </span></li><li><span> <span class="comment">//读取数据</span><span> </span></span></li><li class="alt"><span> <span class="datatypes">float</span><span>* poData = </span><span class="keyword">new</span><span> </span><span class="datatypes">float</span><span>[nReadWidth*nReadHeight]; </span></span></li><li><span> <span class="datatypes">float</span><span>* poOutData = </span><span class="keyword">new</span><span> </span><span class="datatypes">float</span><span>[nReadWidth*nReadHeight]; </span></span></li><li class="alt"><span> poBand->RasterIO(GF_Read,0,nYOffset,nReadWidth,nReadHeight,poData,nReadWidth,nReadHeight,GDT_Float32,0,0); </span></li><li><span> </span></li><li class="alt"><span> <span class="comment">//中间处理过程</span><span> </span></span></li><li><span> </span></li><li class="alt"><span> <span class="comment">//写入数据</span><span> </span></span></li><li><span> <span class="datatypes">int</span><span> pBandList[] = {1}; </span></span></li><li class="alt"><span> <span class="keyword">if</span><span> (1 == nYBlockCount) </span></span></li><li><span> { </span></li><li class="alt"><span> poDstDS->RasterIO(GF_Write,0,nYOffset,nRealWidth,nRealHeight,poOutData,nRealWidth,nRealHeight, </span></li><li><span> GDT_Float32,1,pBandList,0,0,0); </span></li><li class="alt"><span> } </span></li><li><span> </span></li><li class="alt"><span> <span class="keyword">else</span><span> </span></span></li><li><span> { </span></li><li class="alt"><span> <span class="keyword">if</span><span> (i == 0) </span></span></li><li><span> { </span></li><li class="alt"><span> poDstDS->RasterIO(GF_Write,0,nYOffset,nRealWidth,nRealHeight,poOutData,nRealWidth,nRealHeight, </span></li><li><span> GDT_Float32,1,pBandList,0,0,0); </span></li><li class="alt"><span> } </span></li><li><span> <span class="keyword">else</span><span> </span><span class="keyword">if</span><span> (i > 0 && i < nYBlockCount-1) </span></span></li><li class="alt"><span> { </span></li><li><span> poDstDS->RasterIO(GF_Write,0,nYOffset+1,nRealWidth,nRealHeight,poOutData+nRealWidth,nRealWidth,nRealHeight, </span></li><li class="alt"><span> GDT_Float32,1,pBandList,0,0,0); </span></li><li><span> } </span></li><li class="alt"><span> <span class="keyword">else</span><span> </span><span class="keyword">if</span><span>(i == nYBlockCount-1) </span></span></li><li><span> { </span></li><li class="alt"><span> poDstDS->RasterIO(GF_Write,0,nYOffset+1,nRealWidth,nRealHeight,poOutData+nRealWidth,nRealWidth,nRealHeight, </span></li><li><span> GDT_Float32,1,pBandList,0,0,0); </span></li><li class="alt"><span> } </span></li><li><span> } </span></li><li class="alt"><span> </span></li><li><span> <span class="keyword">if</span><span> (poData != NULL) </span></span></li><li class="alt"><span> { </span></li><li><span> <span class="keyword">delete</span><span> []poData; </span></span></li><li class="alt"><span> poData = NULL; </span></li><li><span> } </span></li><li class="alt"><span> </span></li><li><span> <span class="keyword">if</span><span> (poOutData != NULL) </span></span></li><li class="alt"><span> { </span></li><li><span> <span class="keyword">delete</span><span> []poOutData; </span></span></li><li class="alt"><span> poOutData = NULL; </span></li><li><span> } </span></li><li class="alt"><span> </span></li><li><span> } </span></li></ol></div><pre class="cpp" style="display: none;" name="code">//分块处理int nSubHeight = 2000;
int nYBlockCount = (nYsize+nSubHeight-1)/nSubHeight; //计算分块的数量
for (int i = 0; i < nYBlockCount; i ++)
{//实际的块的大小int nRealWidth = nXsize;int nRealHeight = nSubHeight;//读取数据块的大小int nReadWidth = nXsize;int nReadHeight = nSubHeight;int nYOffset = i*nSubHeight;if (1 == nYBlockCount){nRealHeight = nYsize;nReadHeight = nRealHeight;}else{if (i == 0){nReadHeight += 1;}else if (i > 0 && i < nYBlockCount-1){nReadHeight += 2;nYOffset -= 1;}else if(i == nYBlockCount-1){nRealHeight = nYsize-nSubHeight*(nYBlockCount-1);nReadHeight = nRealHeight + 1;nYOffset -= 1;}}//读取数据float* poData = new float[nReadWidth*nReadHeight];float* poOutData = new float[nReadWidth*nReadHeight];poBand->RasterIO(GF_Read,0,nYOffset,nReadWidth,nReadHeight,poData,nReadWidth,nReadHeight,GDT_Float32,0,0);//中间处理过程//写入数据int pBandList[] = {1};if (1 == nYBlockCount){poDstDS->RasterIO(GF_Write,0,nYOffset,nRealWidth,nRealHeight,poOutData,nRealWidth,nRealHeight,GDT_Float32,1,pBandList,0,0,0);}else{if (i == 0){poDstDS->RasterIO(GF_Write,0,nYOffset,nRealWidth,nRealHeight,poOutData,nRealWidth,nRealHeight,GDT_Float32,1,pBandList,0,0,0);}else if (i > 0 && i < nYBlockCount-1){poDstDS->RasterIO(GF_Write,0,nYOffset+1,nRealWidth,nRealHeight,poOutData+nRealWidth,nRealWidth,nRealHeight,GDT_Float32,1,pBandList,0,0,0);}else if(i == nYBlockCount-1){poDstDS->RasterIO(GF_Write,0,nYOffset+1,nRealWidth,nRealHeight,poOutData+nRealWidth,nRealWidth,nRealHeight,GDT_Float32,1,pBandList,0,0,0);}}if (poData != NULL){delete []poData;poData = NULL;}if (poOutData != NULL){delete []poOutData;poOutData = NULL;}}</pre><br><br><p></p>讲到了这里，现在就应该讲怎么计算坡度了吧，采用ARCGIS的计算方法，根据公式来计算，代码如下：<br><div class="dp-highlighter bg_cpp"><div class="bar"><div class="tools"><b>[cpp]</b> <a title="view plain" class="ViewSource" onclick="dp.sh.Toolbar.Command('ViewSource',this);return false;" href="#">view plain</a><span class="tracking-ad" data-mod="popu_168"> <a title="copy" class="CopyToClipboard" onclick="dp.sh.Toolbar.Command('CopyToClipboard',this);return false;" href="#">copy</a><div style="left: 243px; top: 6546px; width: 16px; height: 16px; position: absolute; z-index: 99;"><embed name="ZeroClipboardMovie_3" width="16" height="16" align="middle" id="ZeroClipboardMovie_3" pluginspage="http://www.macromedia.com/go/getflashplayer" src="http://csdnimg.cn/public/highlighter/ZeroClipboard.swf" type="application/x-shockwave-flash" wmode="transparent" flashvars="id=3&width=16&height=16" allowfullscreen="false" allowscriptaccess="always" bgcolor="#ffffff" quality="best" menu="false" loop="false"></div></span><span class="tracking-ad" data-mod="popu_169"> <a title="print" class="PrintSource" onclick="dp.sh.Toolbar.Command('PrintSource',this);return false;" href="#">print</a></span><a title="?" class="About" onclick="dp.sh.Toolbar.Command('About',this);return false;" href="#">?</a></div></div><ol class="dp-cpp"><li class="alt"><span><span class="datatypes">float</span><span> SlopeCal (</span><span class="datatypes">float</span><span>* afRectData, </span><span class="datatypes">float</span><span> fDstNoDataValue,</span><span class="keyword">void</span><span>* pData) </span></span></li><li><span>{ </span></li><li class="alt"><span> <span class="keyword">const</span><span> </span><span class="datatypes">double</span><span> radiansToDegrees = 180.0 / M_PI; </span></span></li><li><span> SlopeOption *psData = (SlopeOption*)pData; </span></li><li class="alt"><span> </span></li><li><span> <span class="datatypes">double</span><span> dx =((afRectData[0] + afRectData[3]*2 + afRectData[6]) - </span></span></li><li class="alt"><span> (afRectData[2]+ afRectData[5]*2 + afRectData[8])) / (psData->dbEwres*8); </span></li><li><span> </span></li><li class="alt"><span> <span class="datatypes">double</span><span> dy =((afRectData[6] + afRectData[7]*2 + afRectData[8]) - </span></span></li><li><span> (afRectData[0]+ afRectData[1]*2 + afRectData[2])) / (psData->dbNsres*8); </span></li><li class="alt"><span> </span></li><li><span> <span class="datatypes">double</span><span> key = (dx *dx + dy * dy); </span></span></li><li class="alt"><span> </span></li><li><span> <span class="keyword">if</span><span>(psData->slopeType == DEGREE_SLOPE) </span></span></li><li class="alt"><span> { </span></li><li><span> <span class="keyword">return</span><span> (</span><span class="datatypes">float</span><span>)(atan(sqrt(key) ) * radiansToDegrees); </span></span></li><li class="alt"><span> <span class="comment">//return key;</span><span> </span></span></li><li><span> } </span></li><li class="alt"><span> <span class="keyword">else</span><span> </span><span class="keyword">if</span><span> (psData->slopeType == PERCENT_SLOPE) </span></span></li><li><span> <span class="keyword">return</span><span> (</span><span class="datatypes">float</span><span>)(100*(sqrt(key) )); </span></span></li><li class="alt"><span> </span></li><li><span> <span class="keyword">return</span><span> 0; </span></span></li><li class="alt"><span> </span></li><li><span>} </span></li></ol></div><pre class="cpp" style="display: none;" name="code">float SlopeCal (float* afRectData, float fDstNoDataValue,void* pData)
{ const double radiansToDegrees = 180.0 / M_PI; SlopeOption psData = (SlopeOption)pData;
double dx =((afRectData[0] + afRectData[3]*2 + afRectData[6]) -(afRectData[2]+ afRectData[5]*2 + afRectData[8])) / (psData->dbEwres*8);double dy =((afRectData[6] + afRectData[7]*2 + afRectData[8]) -(afRectData[0]+ afRectData[1]*2 + afRectData[2])) / (psData->dbNsres*8);double key = (dx *dx + dy * dy);if(psData->slopeType == DEGREE_SLOPE)
{return (float)(atan(sqrt(key) ) * radiansToDegrees);//return key;
}
else if (psData->slopeType == PERCENT_SLOPE)return (float)(100*(sqrt(key) ));return 0;
}
最后以90米分辨率的DEM数据作为测试数据，其得到的结果在ARCGIS中分级渲染如下：
在arcgis中也基于同样的数据生成坡度数据，发现其边缘像素的值不一样，这是因为边缘处理的策略不一样导致的。
三、基于opencl的GPU实现
好了，上一段我们讲解了如何在CPU上实现坡度提取。我们可以知道，图像其实很适合做并行计算，这是因为数据量大和各个像素的计算过程都一样。本文采用OpenCL异构计算开放标准来实现，所谓异构计算，其实就是可以发挥你机器所有的计算资源，包括CPU、GPU、APU等设备。这儿我基于本机上的显卡计算。关于opencl的入门程序，我这边就不多描述了。
这里DEM只有一个波段，所以就不用opencl中图像对象，而是直接采用缓冲区对象，这样所也是为了节约内存以及显存，因为图像对象中每个像素需要存储四个值，所以这里没必要。
opencl的计算函数声明如下：
void SlopeCal_OpenCL(float* poDataIn,float poDataOut,int nWidth,int nHeight,const SlopeOption pSlopeType)
函数中poDataIn接收前面分块的输入数据，poDataOut接收分块的输出数据，nWidth是分块的宽度，nHeight是分块数据的高，pSlopeType算法参数结构体。
这里最主要的工作就是需要把poDataIn，nWidth，nHeight，pSlopeType这几个参数传到内核函数中去，至于poDataOut可以不用传输到显存中去，因为是输出参数。poDataIn就必须作为缓冲区对象，缓冲区对象时opencl内核中可用的一块连续的内存区；其他几个参数是普通参数，SlopeOption结构体要传输到内核函数中的话，就必须在.cl文件中声明和主机端一样的结构体。由于结构体中有double型的参数，opencl是默认禁用掉了double类型，所以需要编译器打开，即在cl文件中要声明
#pragma OPENCL EXTENSION cl_khr_fp64: enable
所以算法的内核函数可以声明如下：
__kernel void slope_kernel( __global const float pSrcData,
__global float *pDestData,const int nWidth,const int nHeight, struct SlopeOption slopeType)
这样，可以用nWidth nHeight个线程并行地计算各个像素的值，所以主机端设置为一个二维的全局N-D Range空间。
其内核函数的实现如下：
- __kernel void slope_kernel( __global const float pSrcData,
- __global float *pDestData,const int nWidth,const int nHeight
- , struct SlopeOption slopeType)
- {
- int j = (int)get_global_id(0);
- int i = (int)get_global_id(1);
- if (j >= nWidth || i >= nHeight)
- return;
- int nTopTmp = i-1;
- int nBottomTmp = i+1;
- int nLeftTep = j-1;
- int nRightTmp = j+1;
- //处理边界情况
- if (0 == i)
- {
- nTopTmp = i;
- }
- if (0 == j)
- {
- nLeftTep = j;
- }
- if (i == nHeight-1)
- {
- nBottomTmp = i;
- }
- if (j == nWidth-1)
- {
- nRightTmp = j;
- }
- float dbRectData[9];
- dbRectData[0] = pSrcData[nTopTmp*nWidth+nLeftTep];
- dbRectData[1] = pSrcData[nTopTmp*nWidth+j];
- dbRectData[2] = pSrcData[nTopTmp*nWidth+nRightTmp];
- dbRectData[3] = pSrcData[i*nWidth+nLeftTep];
- dbRectData[4] = pSrcData[i*nWidth+j];
- dbRectData[5] = pSrcData[i*nWidth+nRightTmp];
- dbRectData[6] = pSrcData[nBottomTmp*nWidth+nLeftTep];
- dbRectData[7] = pSrcData[nBottomTmp*nWidth+j];
- dbRectData[8] = pSrcData[nBottomTmp*nWidth+nRightTmp];
- double dx = ((dbRectData[0] + dbRectData[3]*2 + dbRectData[6]) -
- (dbRectData[2]+ dbRectData[5]*2 + dbRectData[8])) / (slopeType.dbEwres*8);
- double dy =((dbRectData[6] + dbRectData[7]*2 + dbRectData[8]) -
- (dbRectData[0]+ dbRectData[1]*2 + dbRectData[2])) / (slopeType.dbNsres*8);
- double fTmp = (dx *dx + dy dy);
- //计算坡度
- double radiansToDegrees = 180.0/M_PI;
- double fValue = 0;
- if(slopeType.slopeType == DEGREE_SLOPE)
- {
- fValue = atan(sqrt(fTmp) ) * radiansToDegrees;
- }
- else if (slopeType.slopeType == PERCENT_SLOPE)
- fValue = 100*sqrt(fTmp);
- pDestData[i*nWidth+j] = fValue;
- }
__kernel void slope_kernel( __global const float *pSrcData, __global float *pDestData,const int nWidth,const int nHeight , struct SlopeOption slopeType)
{ int j = (int)get_global_id(0); int i = (int)get_global_id(1);if (j >= nWidth || i >= nHeight)return;int nTopTmp = i-1;
int nBottomTmp = i+1;
int nLeftTep = j-1;
int nRightTmp = j+1;//处理边界情况
if (0 == i)
{nTopTmp = i;
}
if (0 == j)
{nLeftTep = j;
}
if (i == nHeight-1)
{nBottomTmp = i;
}
if (j == nWidth-1)
{nRightTmp = j;
}
float dbRectData[9];
dbRectData[0] = pSrcData[nTopTmp*nWidth+nLeftTep];
dbRectData[1] = pSrcData[nTopTmp*nWidth+j];
dbRectData[2] = pSrcData[nTopTmp*nWidth+nRightTmp];dbRectData[3] = pSrcData[i*nWidth+nLeftTep];
dbRectData[4] = pSrcData[i*nWidth+j];
dbRectData[5] = pSrcData[i*nWidth+nRightTmp];dbRectData[6] = pSrcData[nBottomTmp*nWidth+nLeftTep];
dbRectData[7] = pSrcData[nBottomTmp*nWidth+j];
dbRectData[8] = pSrcData[nBottomTmp*nWidth+nRightTmp];double dx = ((dbRectData[0] + dbRectData[3]*2 + dbRectData[6]) -(dbRectData[2]+ dbRectData[5]*2 + dbRectData[8])) / (slopeType.dbEwres*8);double dy =((dbRectData[6] + dbRectData[7]*2 + dbRectData[8]) -(dbRectData[0]+ dbRectData[1]*2 + dbRectData[2])) / (slopeType.dbNsres*8);double fTmp = (dx *dx + dy * dy);//计算坡度
double radiansToDegrees = 180.0/M_PI;
double fValue = 0;
if(slopeType.slopeType == DEGREE_SLOPE)
{fValue = atan(sqrt(fTmp) ) * radiansToDegrees;
}
else if (slopeType.slopeType == PERCENT_SLOPE)fValue = 100*sqrt(fTmp);pDestData[i*nWidth+j] = fValue;
}
而在主机端，需要将数据拷贝到GPU设备端，以及设置设备端内核函数的参数，其主机端主要代码如下：
- //opencl平台搭建
- cl_int status = 0; //状态号码
- static cl_context cxGPUContext = NULL; // OpenCL context
- static cl_command_queue cqCommandQueue = NULL;// OpenCL command que
- static cl_platform_id cpPlatform = NULL; // OpenCL platform
- static cl_device_id cdDevice = NULL; // OpenCL device
- static cl_program cpProgram = NULL; // OpenCL program
- static cl_kernel ckKernel = NULL; // OpenCL kernel
- static bool bInit = 0; //是否初始化了
- if (!bInit)
- {
- OpenCLInit(&cpPlatform,&cdDevice,&cxGPUContext);
- BuildKernel(cpPlatform,cdDevice,cxGPUContext,&cpProgram,&cqCommandQueue);
- ckKernel = clCreateKernel(cpProgram,”slope_kernel”,&status);
- bInit = 1;
- }
- cl_int errNum;
- cl_mem bufIn = clCreateBuffer(cxGPUContext,CL_MEM_READ_WRITE|CL_MEM_COPY_HOST_PTR,
- sizeof(float)*nWidth*nHeight,poDataIn,&errNum);
- cl_mem bufOut = clCreateBuffer(cxGPUContext,CL_MEM_WRITE_ONLY,
- sizeof(float)*nWidth*nHeight,NULL,&errNum);
- //设置参数
- status = clSetKernelArg(ckKernel,0,sizeof(cl_mem),&bufIn);
- status = clSetKernelArg(ckKernel,1,sizeof(cl_mem),&bufOut);
- status = clSetKernelArg(ckKernel,2,sizeof(cl_int),&nWidth);
- status = clSetKernelArg(ckKernel,3,sizeof(cl_int),&nHeight);
- SlopeOption slopeOpt;
- memcpy(&slopeOpt,pSlopeType,sizeof(SlopeOption));
- status = clSetKernelArg(ckKernel,4,sizeof(SlopeOption),&slopeOpt);
- //执行核函数
- size_t globalThreads[] = {nWidth,nHeight};
- status = clEnqueueNDRangeKernel(cqCommandQueue,ckKernel,2,
- NULL,globalThreads,NULL,0,NULL,NULL);
- status = clFinish(cqCommandQueue);
- status = clEnqueueReadBuffer(cqCommandQueue,bufOut,CL_TRUE,0,sizeof(float)*nWidth*nHeight,poDataOut,0,NULL,NULL);
//opencl平台搭建 cl_int status = 0; //状态号码 static cl_context cxGPUContext = NULL; // OpenCL context static cl_command_queue cqCommandQueue = NULL;// OpenCL command que static cl_platform_id cpPlatform = NULL; // OpenCL platform static cl_device_id cdDevice = NULL; // OpenCL device static cl_program cpProgram = NULL; // OpenCL program static cl_kernel ckKernel = NULL; // OpenCL kernelstatic bool bInit = 0; //是否初始化了
if (!bInit)
{OpenCLInit(&cpPlatform,&cdDevice,&cxGPUContext);BuildKernel(cpPlatform,cdDevice,cxGPUContext,&cpProgram,&cqCommandQueue);ckKernel = clCreateKernel(cpProgram,"slope_kernel",&status);bInit = 1;
}cl_int errNum;
cl_mem bufIn = clCreateBuffer(cxGPUContext,CL_MEM_READ_WRITE|CL_MEM_COPY_HOST_PTR,sizeof(float)*nWidth*nHeight,poDataIn,&errNum);
cl_mem bufOut = clCreateBuffer(cxGPUContext,CL_MEM_WRITE_ONLY,sizeof(float)*nWidth*nHeight,NULL,&errNum);//设置参数
status = clSetKernelArg(ckKernel,0,sizeof(cl_mem),&bufIn);
status = clSetKernelArg(ckKernel,1,sizeof(cl_mem),&bufOut);
status = clSetKernelArg(ckKernel,2,sizeof(cl_int),&nWidth);
status = clSetKernelArg(ckKernel,3,sizeof(cl_int),&nHeight);
SlopeOption slopeOpt;
memcpy(&slopeOpt,pSlopeType,sizeof(SlopeOption));
status = clSetKernelArg(ckKernel,4,sizeof(SlopeOption),&slopeOpt);//执行核函数
size_t globalThreads[] = {nWidth,nHeight};
status = clEnqueueNDRangeKernel(cqCommandQueue,ckKernel,2,NULL,globalThreads,NULL,0,NULL,NULL);
status = clFinish(cqCommandQueue);status = clEnqueueReadBuffer(cqCommandQueue,bufOut,CL_TRUE,0,sizeof(float)*nWidth*nHeight,poDataOut,0,NULL,NULL);</pre><p><br></p><p>当然，最后也别忘了释放之前在GPU设备上申请的内存。</p><br><h1><a name="t3"></a>四、性能测试</h1> 通过对几组数据进行测试，分别用一副90米分辨率、宽和高为6001的DEM数据以及12001*12001的数据进行测试。测试环境为N卡GT750M。分别对这两组不同数据进行大量测试，并且取得其平均值，其测试结果见下表。<br><p><br></p><p>数据<span style="white-space:pre;"> </span>CPU时间（毫秒）<span style="white-space:pre;"> </span>GPU时间（毫秒）<span style="white-space:pre;"> </span>加速比</p>6001*6001<span style="white-space:pre;"> </span>5846<span style="white-space:pre;"> </span>1385<span style="white-space:pre;"> </span>4.221<br>12001*12001<span style="white-space:pre;"> </span>23620<span style="white-space:pre;"> </span>5895<span style="white-space:pre;"> </span>4.007<br><br><br>从上表可以看出，其加速比大致维持在4-5之间。不过这个时间是包括了IO时间，如果撇开IO时间，那么统计时间会更短。我觉得这种加速效果不是特别明显，看到有些论文有提到可以提高至少10倍，几十倍，不知道他们是怎么做到的，我觉得程序性能还有提升的空间，还有待挖掘。关于本文的代码已经上传，下载地址为：<a href="http://download.csdn.net/detail/zhouxuguang236/7184841" target="_blank">http://download.csdn.net/detail/zhouxuguang236/7184841</a><br>其中opencl和GDAL环境得自己去下载配置了。<br><br><br><h1><a name="t4"></a>参考文献</h1>1、地理信息系统算法基础、张宏等<br>2、数字高程模型、李志林、朱庆<br>3、OpenCL编程指南、苏金国等翻译<br><br><br><h1><a name="t5"></a>后记</h1> 其实在GIS和遥感这个领域，有很多算法可以进行并行化改造，从而提高我们数据处理的速度，关于OpenCL这个开放标准目前也在发展中，虽说没有CUDA发展得好，但是这个事开放标准，有很好的前途，希望opencl的明天会更好。对于GIS中算法的GPU并行化还有待去探索。 </div></div>
</article>
基于OpenCL的数字地形分析之坡度坡向提取相关推荐
- ArcGIS数字地形分析
数字地形分析 转自:kaitexue 文章目录 数字地形分析原理 基于ArcGIS的数字地形分析操作 DEM的建立 1 栅格表面的创建 2 TIN的创建 3 等高线的创建 4 Terrain(地形数据 ...
- 【ArcGIS|空间分析】数字地形分析
文章目录数字地形分析原理基于ArcGIS的数字地形分析操作DEM的建立1 栅格表面的创建2 TIN的创建3 等高线的创建4 Terrain(地形数据集)的建立基本因子分析1 坡度(栅格表面与TIN表面 ...
- [GIS原理] 9 数字地形分析DTA、数字地形模型DTM、数字高程模型DEM、数字地表模型DSM、不规则三角网TIN
在知识传播途中,向涉及到的相关著作权人谨致谢意! 文章目录 1 数字地形分析(DTA) 1.1 数字地形模型(DTM) 1.1.1 DSM与DEM对比 1.2 数字地形分析研究与应用进展 1.2.1 ...
- gis实验——数字地形分析
跟着汤国安老师学GIS-<GIS实验> http://www.icourse163.org/course/NJNU-1206774803 地势图 分级设色--结果清晰美观(某地区5米分辨率 ...
- GIS—DEM与数字地形分析
DEM:在数学中定义为二维空间的连续函数:H=f(x,y).对点的高程取样然后用矩阵表示.等间距连续的采样可能会错失很多重要的地形特征点. 优点:数据结构简单:便于在GIS环境下处理与显示. 缺点:数 ...
- 【QGIS入门实战精品教程】10.1：QGIS基于DEM数据的地形分析案例教程
本文讲解QGIS中基于DEM数据的地形分析方法,包括:坡度分析.坡向分析.山体阴影.地貌分析.强度指数(地形复杂性). 文章目录 一.加载DEM 二.坡度分析 三.坡向分析 四.山体阴影 五.地貌分析 ...
- GIS空间分析 数字地形分析2 基本地形因子的提取
目录 一.实验名称 二.实验目的 三.实验背景 四.实验准备 1.数据 2.软件 五.实验步骤: 本文数据免费下载 其他GIS空间分析文章 一.实验名称 数字地形分析之基本地形因子的提取 二.实验目的 ...
- GIS空间分析 数字地形分析3 可视性分析
目录 一.实验名称 二.实验目的 三.实验背景 1.可视性分析 2.通视分析 3.视域分析 四.实验准备 1.数据准备 2.软件准备 五.实验步骤 本文数据免费下载 其他GIS空间分析文章 一.实验名 ...
- GIS空间分析 数字地形分析1 地势图的制作
目录 一.实验名称 二.实验目的 三.实验准备 1.数据 2.软件 四.实验步骤 本文数据免费下载 其他GIS空间分析文章 一.实验名称 数字地形分析之地势图的制作 二.实验目的 0 掌握数字高程模型 ...
最新文章
- 【C#】关闭 Window 之后，无法设置 Visibility，也无法调用 Show、ShowDialogor 或 WindowInteropHelper.EnsureHandle...
- C# IPGlobalStatistics获取本机网络流量信息
- PostgreSQL的实践一：初识
- 怎么解决表字段变化引起的MBG 文件变化的问题？
- CentOS 6.5 安装配置Tomcat7服务器
- 跳槽季：跳和不跳之外的第三选择
- 【C语言深入】[001] static 关键字：
- 如何去掉ArrayList重复的id
- Python excel转图片保存
- java各种优秀开源库收集
- verilog语法记录（一）
- g120宏的说明书_西门子变频器G120操作说明书
- 【STM32f401学习之路-01】GPIO实战—点灯、检测按键
- 用IMAP4访问Exchange邮箱
- 软件工程项目之随心所欲—食堂点菜系统
- 浪里淘沙，沙里淘金——我是面试官
- 如何在iPhone手机里创建网页快捷方式图标（Web Clip）
- linux查看gc日志,将GC日志输出到文件
- 2017-8-24股票复盘
- 使用vtk显示标准格式网格文件stl