GDAL 遥感图像处理锐化（Laplace算子、Sobel算子）

GDAL 图像锐化

简介
拉普拉斯（Laplace）算子
- 部分代码：
索贝尔（Sobel）算子
- 部分代码：
处理效果
- 原图
- （Laplace）
- （Sobel）
结尾
参考文章

简介

图像锐化(image sharpening)是补偿图像的轮廓，增强图像的边缘及灰度跳变的部分，使图像变得清晰，分为空间域处理和频域处理两类。图像锐化是为了突出图像上地物的边缘、轮廓，或某些线性目标要素的特征。这种滤波方法提高了地物边缘与周围像元之间的反差，因此也被称为边缘增强。

各种算子是图像锐化的核心，GDAL中没有提供相应的方法，所以我参考了网上C++的图像处理算法和GDAL进行了融合：

拉普拉斯（Laplace）算子

拉普拉斯算子是一个是n维欧几里德空间中的一个二阶微分算子，它的定义如下：

在x方向上

在y方向上

合起来就是

拉普拉斯强调的是图像中灰度的突变，并不强调图像的灰度缓变（灰度缓变由一阶微分，也就是梯度，图像应用是sobel算子，具体下面介绍）

根据上边的表达式，可以确定拉普拉斯算子的模板：

部分代码：

可根据我的另一篇博客（图像平滑），通过改动算法核心写出锐化函数。

    int laplace4[3][3] = { 0, -1, 0, -1, 5, -1, 0, -1, 0 };//laplace锐化模板,4邻域int laplace8[3][3] = { -1, -1, -1, -1, 9, -1, -1, -1, -1 };//laplace锐化模板,8邻域int m,n;//平滑算子也是通过模板进行处理的，所以可以把平滑处理和锐化处理通过一个函数实现int biBitCount = 8;//红色波段处理//分配新像素组的空间int lineByte = (nImgSizeX * biBitCount / 8 + 3) / 4 * 4;rbandsmooth = new unsigned char[lineByte * nImgSizeY];//进行模板操作for (int i = 0; i < nImgSizeY; i++){for (int j = 0; j < nImgSizeX; j++){if (i == 0 || j == 0 || i == nImgSizeY - 1 || j == nImgSizeX - 1)rbandsmooth[i*nImgSizeX + j] = rband[i*nImgSizeX + j];else{int sum = 0;for (m = i - 1; m < i + 2; m++)for (n = j - 1; n < j + 2; n++){sum += (*(rband + m*nImgSizeX  + n  ))*laplace8[n - j + 1][m - i + 1] / 1;}//8位BMP中一个像素值，对应调色板中索引号为该像素值的项所存放的RGB色彩//所以像素值范围为0~255,像素值小于0就取0，大于255就取255sum = (sum > 0) ? sum : 0;sum = (sum > 255) ? 255 : sum;rbandsmooth[i*nImgSizeX + j] = sum;}}}//绿色波段处理//分配新像素素组的空间gbandsmooth = new unsigned char[lineByte * nImgSizeY];//进行模板操作for (int i = 0; i < nImgSizeY; i++){for (int j = 0; j < nImgSizeX; j++){if (i == 0 || j == 0 || i == nImgSizeY - 1 || j == nImgSizeX - 1)gbandsmooth[i*nImgSizeX  + j] = gband[i*nImgSizeX  + j];else{int sum = 0;for (m = i - 1; m < i + 2; m++)for (n = j - 1; n < j + 2; n++){sum += (*(gband + m*nImgSizeX  + n  ))*laplace8[n - j + 1][m - i + 1] / 1;}//8位BMP中一个像素值，对应调色板中索引号为该像素值的项所存放的RGB色彩//所以像素值范围为0~255,像素值小于0就取0，大于255就取255sum = (sum > 0) ? sum : 0;sum = (sum > 255) ? 255 : sum;gbandsmooth[i*nImgSizeX + j]  = sum;}}}//蓝色波段处理//分配新像素素组的空间bbandsmooth = new unsigned char[lineByte * nImgSizeY];//进行模板操作for (int i = 0; i < nImgSizeY; i++){for (int j = 0; j < nImgSizeX; j++){if (i == 0 || j == 0 || i == nImgSizeY - 1 || j == nImgSizeX - 1)bbandsmooth[i*nImgSizeX + j] = bband[i*nImgSizeX + j];else{int sum = 0;for (m = i - 1; m < i + 2; m++)for (n = j - 1; n < j + 2; n++){sum += (*(bband + m*nImgSizeX  + n  ))*laplace8[n - j + 1][m - i + 1] / 1;}//8位BMP中一个像素值，对应调色板中索引号为该像素值的项所存放的RGB色彩//所以像素值范围为0~255,像素值小于0就取0，大于255就取255sum = (sum > 0) ? sum : 0;sum = (sum > 255) ? 255 : sum;bbandsmooth[i*nImgSizeX + j] = sum;}}}

索贝尔（Sobel）算子

Sobel算子是像素图像边缘检测中最重要的算子之一，在机器学习、数字媒体、计算机视觉等信息科技领域起着举足轻重的作用。在技术上，它是一个离散的一阶差分算子，用来计算图像亮度函数的一阶梯度之近似值。在图像的任何一点使用此算子，将会产生该点对应的梯度矢量或是其法矢量。
索贝尔（Sobel）算子定义：

该算子包含两组3x3的矩阵，分别为横向及纵向，将之与图像作平面卷积，即可分别得出横向及纵向的亮度差分近似值。如果以A代表原始图像，Gx及Gy分别代表经横向及纵向边缘检测的图像，其公式如下:

图像的每一个像素的横向及纵向梯度近似值可用以下的公式结合，来计算梯度的大小。

可用以下公式计算梯度方向。

在以上例子中，如果以上的角度Θ等于零，即代表图像该处拥有纵向边缘，左方较右方暗。

部分代码：

与拉普拉斯（Laplace）算子相比较，思路大同小异，主要是Sobel算法有两个算子，需要分别进行计算。

    int sobelx[3][3] = { 1,2,1,0,0,0,-1,-2,-1};//sobelxint sobely[3][3] = { 1,0,-1,2,0,-2,1,0,-1};//sobelyint m,n;//平滑算子也是通过模板进行处理的，所以可以把平滑处理和锐化处理通过一个函数实现unsigned char *imagedatasobelx;unsigned char *imagedatasobely;//分配新像素组的空间int lineByte = (nImgSizeX * biBitCount / 8 + 3) / 4 * 4;rbandsmooth = new unsigned char[lineByte * nImgSizeY];imagedatasobelx = new unsigned char[lineByte * nImgSizeY];imagedatasobely = new unsigned char[lineByte * nImgSizeY];//进行模板操作for (int i = 0; i < nImgSizeY; i++){for (int j = 0; j < nImgSizeX; j++){if (i == 0 || j == 0 || i == nImgSizeY - 1 || j == nImgSizeX - 1){*(imagedatasobelx + i*nImgSizeX + j) = rband[i*nImgSizeX + j];*(imagedatasobely + i*nImgSizeX + j) = rband[i*nImgSizeX + j];}else{int sumx = 0;int sumy = 0;for (int m = i - 1; m < i + 2; m++)for (int n = j - 1; n < j + 2; n++){sumx += (*(rband + m*nImgSizeX + n))*sobelx[n - j + 1][m - i + 1] / 1;sumy += (*(rband + m*nImgSizeX + n))*sobely[n - j + 1][m - i + 1] / 1;}//8位BMP中一个像素值，对应调色板中索引号为该像素值的项所存放的RGB色彩//所以像素值范围为0~255,像素值小于0就取0，大于255就取255sumx = (sumx > 0) ? sumx : 0;sumx = (sumx >255) ? 255 : sumx;imagedatasobelx[i*nImgSizeX + j] = sumx;//sobelxsumy = (sumy > 0) ? sumy : 0;sumy = (sumy >255) ? 255 : sumy;imagedatasobely[i*nImgSizeX + j] = sumy;//sobelyrbandsmooth[i*nImgSizeX + j] =imagedatasobelx[i*nImgSizeX + j] + imagedatasobely[i*nImgSizeX + j];}}}//gbandgbandsmooth = new unsigned char[lineByte * nImgSizeY];imagedatasobelx = new unsigned char[lineByte * nImgSizeY];imagedatasobely = new unsigned char[lineByte * nImgSizeY];//进行模板操作for (int i = 0; i < nImgSizeY; i++){for (int j = 0; j < nImgSizeX; j++){if (i == 0 || j == 0 || i == nImgSizeY - 1 || j == nImgSizeX - 1){*(imagedatasobelx + i*nImgSizeX + j) = gband[i*nImgSizeX + j];*(imagedatasobely + i*nImgSizeX + j) = gband[i*nImgSizeX + j];}else{int sumx = 0;int sumy = 0;for (int m = i - 1; m < i + 2; m++)for (int n = j - 1; n < j + 2; n++){sumx += (*(gband + m*nImgSizeX + n))*sobelx[n - j + 1][m - i + 1] / 1;sumy += (*(gband + m*nImgSizeX + n))*sobely[n - j + 1][m - i + 1] / 1;}//8位BMP中一个像素值，对应调色板中索引号为该像素值的项所存放的RGB色彩//所以像素值范围为0~255,像素值小于0就取0，大于255就取255sumx = (sumx > 0) ? sumx : 0;sumx = (sumx >255) ? 255 : sumx;imagedatasobelx[i*nImgSizeX + j] = sumx;//sobelxsumy = (sumy > 0) ? sumy : 0;sumy = (sumy >255) ? 255 : sumy;imagedatasobely[i*nImgSizeX + j] = sumy;//sobelygbandsmooth[i*nImgSizeX + j] =imagedatasobelx[i*nImgSizeX + j] + imagedatasobely[i*nImgSizeX + j];}}}//bbandbbandsmooth = new unsigned char[lineByte * nImgSizeY];imagedatasobelx = new unsigned char[lineByte * nImgSizeY];imagedatasobely = new unsigned char[lineByte * nImgSizeY];//进行模板操作for (int i = 0; i < nImgSizeY; i++){for (int j = 0; j < nImgSizeX; j++){if (i == 0 || j == 0 || i == nImgSizeY - 1 || j == nImgSizeX - 1){*(imagedatasobelx + i*nImgSizeX + j) = bband[i*nImgSizeX + j];*(imagedatasobely + i*nImgSizeX + j) = bband[i*nImgSizeX + j];}else{int sumx = 0;int sumy = 0;for (int m = i - 1; m < i + 2; m++)for (int n = j - 1; n < j + 2; n++){sumx += (*(bband + m*nImgSizeX + n))*sobelx[n - j + 1][m - i + 1] / 1;sumy += (*(bband + m*nImgSizeX + n))*sobely[n - j + 1][m - i + 1] / 1;}//8位BMP中一个像素值，对应调色板中索引号为该像素值的项所存放的RGB色彩//所以像素值范围为0~255,像素值小于0就取0，大于255就取255sumx = (sumx > 0) ? sumx : 0;sumx = (sumx >255) ? 255 : sumx;imagedatasobelx[i*nImgSizeX + j] = sumx;//sobelxsumy = (sumy > 0) ? sumy : 0;sumy = (sumy >255) ? 255 : sumy;imagedatasobely[i*nImgSizeX + j] = sumy;//sobelybbandsmooth[i*nImgSizeX + j] =imagedatasobelx[i*nImgSizeX + j] + imagedatasobely[i*nImgSizeX + j];}}}

处理效果

原图

（Laplace）

这里有个问题，在拉普拉斯算子中间是，4和8，但是用原版算子的处理结果是这样的：

看起来是做了边缘提取；
后来经过修改，将4邻域和8邻域改成了5和9，就出了这个结果：

看起来比原来清楚许多，实际应用效果会更明显。

（Sobel）

索贝尔算子还是有问题，可能是因为从灰度图像的锐化移植为真彩色图像时的问题，还需要进行修改，结果也像是做了边缘提取：

结尾

因为之前没有太涉及到算法问题，而且GDAL的锐化方面还没有太多的过来人，所以经验不多，还需要人们探索，如有不足请多指教。

参考文章

https://www.cnblogs.com/fydeblog/p/6748411.html
https://blog.csdn.net/wanty_chen/article/details/80336986
https://blog.csdn.net/White_Idiot/article/details/51794364