深度剖析图像处理

什么是边缘检测

边缘检测(Edge Detection)就是提取图像中的边缘点(Edge Point)。边缘点是与周围像素相比灰度值有阶跃变化或屋顶状变化的像素。边缘常存在于目标与背景之间、目标与目标之间、目标与其影子之间。

在图像处理和图像分析中，经常要用到边缘(Edge)、边界(Boundary)、轮廓(Contour)等术语。一般来说，边缘指的是边缘点，它不能被称为边缘线。边界指的是图像中不同区域之间的分界线，比如不同灰度、不同颜色的区域之间的分界线，它是线而不是点，可以被称为边界线。轮廓一般是指目标的轮廓，目标就是语义明确的区域，轮廓一般是在二值图像中围绕白色区域的闭合曲线。

边缘检测算法和边界线检测算法一般作用于灰度图像，对于二值图像进行边缘检测是没有意义的。轮廓一定是闭合的，但边界线不一定闭合，比如道路区域与道边植被的边界线；边缘点最多是断断续续的线段，不保证连续，更不保证闭合。掌握边缘、边界、轮廓的准确术语是非常必要的。

边缘类型

边缘检测是一种邻域运算，即一个像素是否是边缘点是由其所处的邻域决定的。在一定大小的邻域内，边缘分为阶跃边缘(step edge)和屋顶状边缘(roof edge)两种类型。下面以一维信号为例，分析这两种不同类型的边缘的导数特征

求导与差分

在边缘检测中，导数的计算通常采用两种方法：

将邻域从离散空间变换到连续空间，得到解析描述，然后进行求导操作**。**

具体做法是，先将**邻域按照一定的数学模型(曲线拟合、曲面拟合)得到其在连续空间中的解析描述，然后对此解析描述进行求导，得到边缘点。解析描述求导得到的导数位置是有小数位的，比如在位置4.17处取得导数最大值，即边缘点的位置是在4.17而不是像素的整数坐标4。这样得到的边缘点位置精度能够小于1个像素，因此又将此方法称之为亚像素(sub pixel)边缘检测。在已知数学模型的指导下，目前工业界做到的边缘检测最高精度为1/50个像素。**亚像素边缘检测能够在大大节省硬件成本的同时，得到高的边缘检测精度，是图像测量中的常用方法。
直接用差分(difference)代替求导。导数的公式见下式，如果令dx=1，即得到是差分描述，式（4-2）是差分描述的x方向偏导数，式(4-3)是差分描述的y方向偏导数。

边缘强度与边缘方向

导数是有大小也有方向的，因此边缘也有强弱与方向，分别叫做边缘强度(edge intensity)和边缘方向(edge direction)，边缘强度即边缘的幅值(magnitude)。用M(x,y)代表边缘的强度，θ(x,y)代表边缘的方向，有：

实例

我们该如何提取这张图片的边缘呢?

test

首先当然是需要我们写一个函数来把24位彩色图像转化为8位灰度值图像

//24位彩色图像转8位灰度值
//rgbImage原始图像
//grayImage输出灰度图像
//width,height图片的宽和高
void convertToGray(uint8_t* rgbImage, uint8_t* grayImage, int width, int height)
{
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            // 获取当前像素的 RGB 分量
            uint8_t r = rgbImage[3 * (y * width + x) + 0];
            uint8_t g = rgbImage[3 * (y * width + x) + 1];
            uint8_t b = rgbImage[3 * (y * width + x) + 2];

            // 计算灰度值（常用的加权平均法）
            // 这里使用的加权系数是常见的：R: 0.299, G: 0.587, B: 0.114
            uint8_t gray = (uint8_t)(0.299 * r + 0.587 * g + 0.114 * b);

            // 将灰度值写入灰度图像数组
            grayImage[y * width + x] = gray;
        }
    }
}

我们来看看效果

那接下来就需要用边缘检测来提取边缘了

一阶微分算子

根据边缘类型及其导数特征，可以设计不同的检测算法。下面讲述几种常用的边缘检测算法，习惯上称为边缘检测算子(Operator)。当使用差分时，一般写成模板的表示形式。

对于阶跃边缘而言，边缘点处的导数特征是“一阶导数取极值”。若边缘点处的一阶导数为正值，则其为最大值；反之，则为最小值，即在边缘点处的导数绝对值最大。

基于一阶导数的边缘检测算子称为一阶微分算子，常用的一阶微分算子有梯度算子、罗伯特算子、索贝尔算子、Prewitt算子，Robinson算子、Kirsch算子等。

梯度算子

梯度的本意是一个向量（矢量），表示某一函数在该点处的方向导数沿着该方向取得最大值，即函数在该点处沿着该方向（此梯度的方向）变化最快，变化率最大（为该梯度的模）。梯度的含义和边缘点是一致的，因此产生了边缘检测的梯度算子(Gradient Operator)

写成模板的情况就是如下图所示：

但是上面我们只是给出了像素(x,y)的边缘强度，称为梯度值；但是它是不是边缘点，还需要一定的约束条件，比如，设定当Gradient(x,y)≥threshold时，像素(x,y)才是边缘点，threshold称为阈值。

这就是用梯度算子计算的结果，我们上代码看看

void RmwGradientGryImg(uint8_t* pGryImg, int width, int height, uint8_t* pGrdImg)
{
    uint8_t* pGry, * pGrd;
    int dx, dy;
    int x, y;

    for (y = 0, pGry = pGryImg, pGrd = pGrdImg; y < height - 1; y++)
    {
        for (x = 0; x < width - 1; x++, pGry++)
        {
            dx = *pGry - *(pGry + 1);
            dy = *pGry - *(pGry + width);
            int gradient = (int)(sqrt(dx * dx * 1.0 + dy * dy));
            *pGrd++ = (gradient > 255) ? 255 : gradient;
        }
        *pGrd++ = 0; //尾列不做,边缘强度赋0
        pGry++;
    }
    memset(pGrd, 0, width); //尾行不做,边缘强度赋0
}

这里简单说一下这个函数

函数参数：
pGryImg：输入的灰度图像数据指针。
width：图像的宽度。
height：图像的高度。
pGrdImg：输出的梯度图像数据指针。

双层循环：
外层循环 for (y = 0, pGry = pGryImg, pGrd = pGrdImg; y<height-1; y++) 遍历图像的每一行，pGry 和 pGrd 分别指向当前行的灰度图像数据和梯度图像数据。
内层循环 for (x = 0; x<width-1; x++, pGry++) 遍历当前行的每个像素，pGry 指向当前像素的灰度值。

计算梯度：
梯度的计算采用的是简单的基于像素差值的方法，分别计算水平方向和垂直方向的梯度。
dx = *pGry-*(pGry+1)：水平方向的梯度，计算当前像素和右侧像素的灰度差值。
dy = *pGry-*(pGry+width)：垂直方向的梯度，计算当前像素和下方像素的灰度差值。
sqrt(dx*dx*1.0+dy*dy)：使用欧式距离公式计算梯度幅值。
min(255, ...)：确保梯度幅值不超过255，限定在0到255之间。

内存清零：
memset(pGrd, 0, width)：清零输出图像的最后一行，因为最后一行的梯度值未计算。

返回：
函数返回，处理完成。

看看实际运行的效果是什么样子的

效果并不是很好，我们把原始图像做一下图像增强（灰度值均衡化）试试

换一种图像增强的方法，加一下反相试试

先简单写一下反相的函数

void invertImage(uint8_t *image, int width, int height) {
    for (int i = 0; i < width * height; i++) {
        image[i] = 255 - image[i];
    }
}

这次使用了线性拉伸并简单处理了一下参数

效果还算不错，那如果加一下阈值呢？就像上面说的，我们重新写一下这个函数

//梯度算子加阈值
//pGryImg：输入的灰度图像数据指针。
//width：图像的宽度。
//height：图像的高度。
//pGrdImg：输出的梯度图像数据指针
void RmwGradientGryImgPlus(uint8_t* pGryImg, int width, int height, uint8_t* pGrdImg, int threshold)
{
    uint8_t* pGry, * pGrd;
    int dx, dy;
    int x, y;

    for (y = 0, pGry = pGryImg, pGrd = pGrdImg; y < height - 1; y++)
    {
        for (x = 0; x < width - 1; x++, pGry++)
        {
            dx = *pGry - *(pGry + 1);
            dy = *pGry - *(pGry + width);
            int gradient = (int)(sqrt(dx * dx * 1.0 + dy * dy));
            *(pGrd++) = (gradient > threshold) ? min(255, gradient) : 0;
        }
        *(pGrd++) = 0; //尾列不做,边缘强度赋0
        pGry++;
    }
    memset(pGrd, 0, width); //尾行不做,边缘强度赋0
}

看一下效果

画面确实变得更干净了

罗伯特算子

梯度算子的计算只涉及到了3个像素，只在水平和垂直方向做差分。罗伯特算子（Roberts Operator)给出了一个4个像素之间进行运算的算子，分别在两个对角线方向做差分。

由于对角线上2个像素之间的距离为√2，所以罗伯特算子的∆x和∆y采用对角线差分后，不再采用√(∆_x^2+∆_y2 )，其描述见下式。罗伯特算子取∆x绝对值与∆y绝对值中的最大值。

我们来看演示

那么罗伯特算子好在哪呢？

罗伯特算子去掉了梯度算子的开方运算，计算复杂度也降低了不少。

void RmwRobertsGryImg(uint8_t *pGryImg, int width, int height, uint8_t *pRbtImg)
{
    uint8_t *pGry, *pRbt;
    int dx, dy;
    int x, y;

    for (y = 0, pGry = pGryImg, pRbt = pRbtImg; y < height - 1; y++)
    {
        for (x = 0; x < width - 1; x++, pGry++)
        {
            dx = *pGry - *(pGry + width + 1);
            dy = *(pGry + 1) - *(pGry + width);
            *pRbt++ = (uint8_t)(dx > dy ? dx : dy); // 使用三目运算符选择较大的值
        }
        *pRbt++ = 0; // 尾列不做, 边缘强度赋0
        pGry++;
    }
    memset(pRbt, 0, width); // 尾行不做, 边缘强度赋0
}

我们来看结果

为什么看起来的效果好像还没有梯度算子的结果好呢

是因为图像太复杂没有滤波，那么在边缘计算前能不能先滤波？

我们来看下一个算子

索贝尔算子

在一幅噪声较大的图像中，如果不进行图像平滑就进行边缘检测，必然会在边缘图像中产生噪声干扰。

因此，索贝尔算子(Sobel Operator)中，在求∆x和∆y前，先进行滤波。在求∆x前，先执行如下图的所示的高斯均值滤波；在求∆y前，先执行如下图的所示的高斯均值滤波。

另外，索贝尔算子进一步拉大进行差分的2个像素之间的距离，∆x和∆y采用如下的模板形式

我们来看示意图

在图像1中，红色方块代表当前像素(x,y)，先执行图4-8(a)所示的高斯滤波，用D ̅、E ̅代表滤波后的值，则得到：

执行∆x模板，则有∆x=D ̅-E ̅=(A+2D+F)-(C+2E+H)

另外，索贝尔算子在对∆x和∆y的使用上，采用了它们的绝对值相加的形式

看效果：

对原始灰度图像执行索贝尔算子得到的结果，图中虚线框所示的边缘变成了双线宽。

我们来写代码

//索贝尔算子
void RmwSobelGryImg(uint8_t* pGryImg, int width, int height, uint8_t* pSbImg)
{
    uint8_t* pGry, * pSb;
    int dx, dy;
    int x, y;

    memset(pSbImg, 0, width); // 首行不做, 边缘强度赋0
    for (y = 1, pGry = pGryImg + width, pSb = pSbImg + width; y < height - 1; y++)
    {
        *pSb++ = 0; // 首列不做, 边缘强度赋0
        pGry++;
        for (x = 1; x < width - 1; x++, pGry++)
        {
            // 求dx
            dx = *(pGry - 1 - width) + (*(pGry - 1) * 2) + *(pGry - 1 + width);
            dx -= *(pGry + 1 - width) + (*(pGry + 1) * 2) + *(pGry + 1 + width);
            // 求dy
            dy = *(pGry - width - 1) + (*(pGry - width) * 2) + *(pGry - width + 1);
            dy -= *(pGry + width - 1) + (*(pGry + width) * 2) + *(pGry + width + 1);
            // 结果
            *pSb++ = (uint8_t)min(255, abs(dx) + abs(dy));
        }
        *pSb++ = 0; // 尾列不做, 边缘强度赋0
        pGry++;
    }
    memset(pSb, 0, width); // 尾行不做, 边缘强度赋0
}

梯度算子、罗伯特算子、索贝尔算子的比较

以下从4个方面对梯度算子、罗伯特算子、索贝尔算子进行比较。

•1. 偏导数∆x和∆y的求取

梯度算子在3个像素之间进行运算，只在水平和垂直方向做差分，做差分的2个像素之间的距离为1。

罗伯特算子在4个像素之间进行运算，分别在两个对角线方向做差分，做差分的2个像素之间的距离为√2。

索贝尔算子在8个像素之间进行运算，只在水平和垂直方向做差分，做差分的2个像素之间的距离为2。

•2. 是否“先平滑后求导”

索贝尔算子在差分之前，进行了加权均值滤波对图像进行平滑（加权函数采用了高斯函数），因此索贝尔具有滤除噪声的效果。梯度算子和罗伯特算子都没有进行平滑。“先平滑后求导”是边缘检测的通用策略，一般在执行梯度算子和罗伯特算子前是需要使用另外的步骤做图像平滑，索贝尔算子则是把平滑写到了算子中。

•3. 边缘强度的大小

按照式(4-4)边缘强度的定义，梯度算子是严格遵守的，罗伯特算子是取∆x绝对值和∆y绝对值中的最大值，索贝尔算子是取∆x与∆y的绝对值之和。而且，索贝尔算子在高斯滤波后没有除以4，所以又相当于∆x、∆y放大了4倍。对于边缘强度，罗伯特算子、梯度算子、索贝尔算子之间的数值关系大致如下：

•4. 邻域与边缘宽度

梯度算子、罗伯特算子的计算只涉及到了2行2列，所以它们得到的边缘宽度是1个像素；索贝尔算子涉及到了3行3列，所以它得到的边缘宽度是2个像素，边缘变成了双线宽。

方向模板

若是能根据边缘的具体走向求偏导数，则边缘强度值应该会更准确。

因此在实际应用中，是先假定了有限的几个边缘方向，再对这些假定的每个边缘方向设置一个特定的模板，计算每个模板的边缘强度，从中选择最大的边缘强度作为边缘强度的结果，而且该最大边缘强度对应模板的方向就认为是边缘的方向。

常用基于方向模板的边缘检测算子有：Prewitt算子，Robinson算子、Kirsch算子。Prewitt算子使用4个方向模板，Robinson算子和Kirsch算子都使用8个方向模板。这些算子都是先采用了均值滤波，然后进行差分计算。

Prewitt算子

Prewitt算子设定了0°、45°、90°和135°，共计4种边缘方向；根据这4种边缘方向，分别设计了4个模板

对每个像素分别计算这个4个模板的值，取绝对值最大者作为该像素的边缘强度Prewitt(x,y)。同时该最大值对应模板的方向作为该像素的边缘方向(与边缘的走向相差90°，因为显然边缘走向的法线方向上的导数最大)。若把这些模板中为“0”点(空白处)的连成一条直线，可以发现这些模板强调了水平线、135°斜线、竖直线、45°斜线的检测。Prewitt算子强调对直线的检测，对于上述走向的直线，总有一个模板的输出值最大。

void RmwPrewittGryImg(uint8_t *pGryImg, int width, int height, uint8_t *pPRTImg)
{
    uint8_t *pGry, *pPRT;
    int dx, dy, d45, d135, v1, v2;
    int x, y;

    memset(pPRTImg, 0, width); // 首行不做, 边缘强度赋0
    for (y = 1, pGry = pGryImg + width, pPRT = pPRTImg + width; y < height - 1; y++)
    {
        *pPRT++ = 0; // 首列不做, 边缘强度赋0
        pGry++;
        for (x = 1; x < width - 1; x++, pGry++)
        {
            // 求dx
            dx = *(pGry - 1 - width) + *(pGry - 1) + *(pGry - 1 + width);
            dx -= *(pGry + 1 - width) + *(pGry + 1) + *(pGry + 1 + width);
            // 求dy
            dy = *(pGry - width - 1) + *(pGry - width) + *(pGry - width + 1);
            dy -= *(pGry + width - 1) + *(pGry + width) + *(pGry + width + 1);
            // 求45度
            d45 = *(pGry - width - 1) + *(pGry - width) + *(pGry - 1);
            d45 -= *(pGry + width + 1) + *(pGry + width) + *(pGry + 1);
            // 求135度
            d135 = *(pGry - width) + *(pGry - width + 1) + *(pGry + 1);
            d135 -= *(pGry + width - 1) + *(pGry + width) + *(pGry - 1);
            // 结果
            v1 = abs(dx) > abs(dy) ? abs(dx) : abs(dy);
            v2 = abs(d45) > abs(d135) ? abs(d45) : abs(d135);
            *pPRT++ = (uint8_t)((v1 > v2) ? ((v1 > 255) ? 255 : v1) : ((v2 > 255) ? 255 : v2));
        }
        *pPRT++ = 0; // 尾列不做, 边缘强度赋0
        pGry++;
    }
    memset(pPRT, 0, width); // 尾行不做, 边缘强度赋0
}

我们看一下实现效果

Robinson算子

Robinson算子除(a)外的7个模板，都是由其上个模板顺时针旋转1个像素得到的。若是把模板中的负数值合并成一个区域，可以看出该算子强调了对角点的检测；对于各种形状的角点，总有一个模板的输出值最大。

有兴趣大家可以查一下啊，这里就不再过多解释

二阶微分算子

一阶微分算子能够得到边缘强度，但是需要再加上一定的条件约束，比如设置个阈值，才能判定一个像素是不是边缘点。

通过对边缘类型及其导数的分析可知，阶跃边缘的导数特征除了“一阶导数取极值”外，还有“二阶导数过零点”。因此可以采用二阶导数，利用过零点得到边缘点，这样就不需要其他的条件了。

拉普拉斯算子

拉普拉斯算子(Laplacian Operator)是近似给出二阶导数的流行方法，其使用3×3的邻域，给出了4-邻接和8-邻接的邻域的2种模板

对如图所示的原始灰度图像执行4邻域拉普拉斯算子，得到的结果如图所示，图中虚线框所示的位置上发生了过零点（导数由负数变到了正数），此处即边缘。

在拉普拉斯算子的结果图像中，可以发现过零点位置刚好就是边缘的位置。由于过零点是在像素之间，不在整数坐标上，所以在提取边缘点时，往往采取下面的策略：

当一个像素的二阶导数大于0，其邻域内有像素的二阶导数小于0或等于0，则该像素被标记为边缘点。

沈俊算子

唯一一个以国人命名的算子

沈俊教授同样提出了先滤波后求导的边缘检测方法（J. Shen and S. Castan, An optimal linear operator for step edge detection, CVGIP: Graphical Models and Image Processing, Vol. 54 No.2, Mar. 1992, pp.112 – 133），即沈俊算子(ShenJun Edge Operator)。

沈教授在阶跃边缘和可加白噪声模型下，就信噪比最大准则，证明了图像平滑的最佳滤波器是对称的指数函数，形式如下：

显然，当a_0越大时，c2就越小，T(j,i)就越陡越窄，相当于滤波邻域就越小，压制噪声的能力就弱，图像模糊程度就越小，边缘定位的精度就越高。

在算子实现上，沈教授对图像分别按行、按列各进行两次先正方向再反方向的递推滤波实现(|j|、|i|的优点)，等价于用上述指数函数进行图像滤波；证明了滤波结果减去原始灰度值得到的差值乘以2c1〖ln〗^c2，约等于其二阶导数的值。沈俊算子的实现过程如下：

沈俊算子的实现过程
step.1 对每行从左向右进行：
𝑔1(0,𝑦)=𝑔(0,𝑦)，
𝑔1(𝑥,𝑦)=𝑔1(𝑥−1,𝑦)+𝑎_0×(𝑔(𝑥,𝑦)−𝑔1(𝑥−1,𝑦))，𝑥=1,2,⋯𝑤𝑖𝑑𝑡ℎ−1。
step.2 对每行从右向左进行：
𝑔2(𝑤𝑖𝑑𝑡ℎ−1,𝑦)=𝑔1(𝑤𝑖𝑑𝑡ℎ−1,𝑦)，
𝑔2(𝑥,𝑦)=𝑔2(𝑥+1,𝑦)+𝑎_0×(𝑔1(𝑥,𝑦)−𝑔2(𝑥+1,𝑦))，𝑥=𝑤𝑖𝑑𝑡ℎ−2, ⋯, 1,0。
step.3 对每列从上向下进行：
𝑔3(𝑥,0)=𝑔2(𝑥,0)，
𝑔3(𝑥,𝑦)=𝑔3(𝑥,𝑦−1)+𝑎_0×(𝑔2(𝑥,𝑦)−𝑔3(𝑥,𝑦−1))，𝑦=1,2,⋯ℎ𝑒𝑖𝑔ℎ𝑡−1。
step.4 对每列从下向上进行：
𝑔4(𝑥,ℎ𝑖𝑔ℎ𝑡−1)=𝑔3(𝑥,ℎ𝑖𝑔ℎ𝑡−1)，
𝑔4(𝑥,𝑦)=𝑔4(𝑥,𝑦+1)+𝑎_0×(𝑔3(𝑥,𝑦)−𝑔4(𝑥,𝑦+1))，𝑦=ℎ𝑒𝑖𝑔ℎ𝑡−2, ⋯, 1,0。
step.5 对每个像素(𝑥,𝑦)执行𝑆𝐽(𝑥,𝑦)=𝑔4(𝑥,𝑦)−𝑔(𝑥,𝑦)，得到二阶导数𝑆𝐽(𝑥,𝑦)。
step.6 对每个像素𝑆𝐽(𝑥,𝑦)进行过零点检测得到边缘点。

我们来看程序

//沈俊算子
//pGryImg 和 pTmpImg 是指向 uint8_t 类型的指针，它们分别指向原始灰度图像数据和辅助图像数据。
//width 和 height 是整型参数，表示图像的宽度和高度。
//a0 是双精度浮点型参数，表示滤波系数。
//pSJImg 是指向 uint8_t 类型的指针，它指向了输出的图像数据。
void RmwShenJunGryImg(uint8_t* pGryImg,uint8_t* pTmpImg, int width, int height, double a0, uint8_t* pSJImg)
{
    uint8_t* pGry, * pCur, * pSJ, * pEnd;
    int LUT[512], * ALUT; // a0查找表
    int x, y, pre, dif;

    // Step 1: 初始化查找表
    a0 = (a0 < 0.01) ? 0.01 : ((a0 > 0.99) ? 0.99 : a0); // 安全性检查
    // a0查找表, 进行了四舍五入
    ALUT = LUT + 256;
    for (ALUT[0] = 0, dif = 1; dif < 256; dif++)
    {
        ALUT[dif] = (int)(dif * a0 + 0.5);
        ALUT[-dif] = (int)(-dif * a0 - 0.5);
    }

    // Step 2: 递推实现指数滤波
    // 按行滤波
    for (y = 0, pGry = pGryImg, pCur = pTmpImg; y < height; y++)
    {
        // 从左向右: p1(y,x) = p1(y,x-1) + a * [p(y,x) - p1(y,x-1)]
        *(pCur++) = pre = *(pGry++);
        for (x = 1; x < width; x++, pGry++)
            *(pCur++) = pre = pre + ALUT[*pGry - pre];
        pCur--; // 回到行尾
        // 从右向左: p2(y,x) = p2(y,x+1) - a * [p1(y,x) - p2(y,x+1)]
        for (x = width - 2, pCur = pCur - 1; x >= 0; x--)
            *(pCur--) = pre = pre + ALUT[*pCur - pre];
        pCur += (width + 1); // 回到下一行的开始
    }
    // 按列滤波
    for (x = 0, pCur = pTmpImg; x < width; x++, pCur = pTmpImg + x)
    {
        // 从上向下: p3(y,x) = p3(y-1,x) + a * [p2(y,x) - p3(y-1,x)]
        pre = *pCur;
        for (y = 1, pCur += width; y < height; y++, pCur += width)
            *pCur = pre = pre + ALUT[*pCur - pre];
        pCur -= width; // 回到列尾
        // 从下向上: p4(i,j) = p4(i+1,j) + a * [p3(i,j) - p4(i+1,j)]
        for (y = height - 2, pCur -= width; y >= 0; y--, pCur -= width)
            *pCur = pre = pre + ALUT[*pCur - pre];
    }

    // Step 3: 正导数=1，负导数为0，0必须也是0
    pEnd = pTmpImg + width * height;
    for (pCur = pTmpImg, pGry = pGryImg; pCur < pEnd; pGry++)
    {
        *(pCur++) = (*pCur > *pGry);
    }

    // Step 4: 过零点检测
    memset(pSJImg, 0, width * height); // 边缘强度赋0
    pSJ = pSJImg + width;
    pCur = pTmpImg + width; // 首行不做 
    for (y = 1; y < height - 1; y++)
    {
        pSJ++; pCur++;  // 首列不做
        for (x = 1; x < width - 1; x++, pGry++, pCur++, pSJ++)
        {
            if (*pCur) // 正导数
            {
                // 下面使用4邻域, 边缘为8连通, 不能保证4连通; 使用8邻域才能保证边缘4连通
                if ((!*(pCur - 1)) || // 左, 必须<=0, 不能<0
                    (!*(pCur + 1)) || // 右, 必须<=0, 不能<0
                    (!*(pCur - width)) || // 上, 必须<=0, 不能<0
                    (!*(pCur + width)))   // 下, 必须<=0, 不能<0
                {
                    *pSJ = 255; // 周围有导数小于等于0
                }
            }
        }
        pSJ++; pCur++;  // 尾列不做
    }
}

当滤波系数为0.1时，我们来看效果

沈俊算子能够得到闭合的边缘。一种边缘检测方法能够得到一般要用图像分割才能得到的目标轮廓，会具有很高的实用价值。

沈俊算子只需要一个参数a_0，且a_0语义明确，而且沈俊算子的代码量非常小，所以沈俊算子使用起来非常方便。

测试

那么该如何更好的得到我们原始图像的边缘呢

可能一个边缘处理都没有办法得到很好的效果，我们可以结合多个处理办法

我们可以看到效果比较好的是索贝尔算子，那么我们结合一下

//沈俊算子加索贝尔算子
//pGryImg：指向原始灰度图像数据的指针
//pTmpImg：指向辅助图像数据的指针
//width：图像的宽度
//height：图像的高度
//a0：这是沈俊算子的参数，用于控制边缘检测的灵敏度。
//grdThre：这是Sobel算子的梯度阈值
//pEdgeImg：最终边缘图像数据的指针
void RmwExtractRiceEdge(uint8_t* pGryImg,uint8_t* pTmpImg,int width,int height,double a0, int grdThre, uint8_t* pEdgeImg)
{
    // step.1------------沈俊算子-----------------------//
    RmwShenJunGryImg(pGryImg, pTmpImg, width, height, a0, pEdgeImg);
    // step.2------------Sobel算子----------------------//
    RmwSobelGryImg(pGryImg, width, height, pTmpImg);
    // step.3------------二者融合-----------------------//
    for (int i = 0; i < width * height; i++)
    {
        *(pEdgeImg + i) = (pEdgeImg[i] && (pTmpImg[i] > grdThre)) * 255;
    }
    // step.4------------结束---------------------------//
    return;
}

我们换一张图片试一下

总结

那么到目前为止，我们已经可以成功的进行边缘检测

接下来就是边缘增强，边缘分割等相关内容，会尽快更新

源码

IDP.h

#pragma once

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include <nmmintrin.h>
 
uint8_t* readGrayScaleBMP(const char* filename, int* width, int* height);//读取8位灰度图片
void saveGrayScaleBMP(const char* filename, const uint8_t* imageData, int width, int height);// 将8位灰度图像数据保存为BMP文件
uint8_t* readColorBMP(const char* filename, int* width, int* height);//读取24位彩色图像的BMP文件
void saveColorBMP(const char* filename, const uint8_t* imageData, int width, int height);//将24位彩色图像数据保存为BMP文件
void convertToGray(uint8_t* rgbImage, uint8_t* grayImage, int width, int height);//24位彩色图像转8位灰度值
void LinearStretchDemo(uint8_t* pGryImg, int width, int height, double k, double b);//灰度线性拉伸
void GetHistogram(uint8_t* pImg, int width, int height, int* histogram);//统计图像灰度值
void GetBrightContrast(int* histogram, double* bright, double* contrast);//亮度和对比度
void RmwHistogramEqualize(uint8_t* pGryImg, int width, int height);//直方图均衡化
void RmwLogTransform(uint8_t* pGryImg, int width, int height);//对数变换
void RmwAvrFilterBySumCol(uint8_t* pGryImg, int width, int height, int M, int N, uint8_t* pResImg);//基于列积分的快速均值滤波
void RmwDoSumGryImg(uint8_t* pGryImg, int width, int height, int* pSumImg);//基于列积分的积分图实现
void RmwDoSumGryImg_SSE(uint8_t* pGryImg, int width, int height, int* pSumImg);//基于SSE的积分图实现
void RmwAvrFilterBySumImg(int* pSumImg, int width, int height, int M, int N, uint8_t* pResImg);//基于积分图的快速均值滤波  
void GetMedianGry(int* histogram, int N, int* medGry);//求灰度值中值
double RmwMedianFilter(uint8_t* pGryImg, int width, int height, int M, int N, uint8_t* pResImg);//中值滤波
void RmwBinImgFilter(uint8_t* pBinImg, int width, int height, int M, int N, double threshold, uint8_t* pResImg);//二值滤波
void RmwGradientGryImg(uint8_t* pGryImg, int width, int height, uint8_t* pGrdImg);//梯度算子
void RmwGradientGryImgPlus(uint8_t* pGryImg, int width, int height, uint8_t* pGrdImg, int threshold);//梯度算子加阈值
void invertImage(uint8_t* image, int width, int height);//反相
void RmwRobertsGryImg(uint8_t* pGryImg, int width, int height, uint8_t* pRbtImg);//罗伯特算子
void RmwSobelGryImg(uint8_t* pGryImg, int width, int height, uint8_t* pSbImg);//索贝尔算子
void RmwPrewittGryImg(uint8_t* pGryImg, int width, int height, uint8_t* pPRTImg); //Prewitt算子
void RmwShenJunGryImg(uint8_t* pGryImg, uint8_t* pTmpImg, int width, int height, double a0, uint8_t* pSJImg);//沈俊算子
void RmwExtractRiceEdge(uint8_t* pGryImg, uint8_t* pTmpImg, int width, int height, double a0, int grdThre, uint8_t* pEdgeImg);//索贝尔＋沈俊算子

IDP.c

#define _CRT_SECURE_NO_WARNINGS 1

#include "IDP.h"

//读取8位灰度图片
//filename：字符数组的指针，用于指定要保存的图像文件的名称或路径。
//imageData：无符号 8 位整型数据的指针，代表要保存的图像数据。
//width：图像的宽度。
//height：图像的高度。
uint8_t* readGrayScaleBMP(const char* filename, int* width, int* height) 
{
    FILE* file = fopen(filename, "rb");
    if (!file) {
        fprintf(stderr, "Error opening file %s\n", filename);
        return NULL;
    }

    // 读取BMP文件头部信息
    uint8_t bmpHeader[54];
    fread(bmpHeader, 1, 54, file);

    // 从文件头部提取图像宽度和高度信息
    *width = *(int*)&bmpHeader[18];
    *height = *(int*)&bmpHeader[22];

    // 分配存储图像数据的内存
    uint8_t* imageData = (uint8_t*)malloc(*width * *height);
    if (!imageData) {
        fprintf(stderr, "内存分配失败\n");
        fclose(file);
        return NULL;
    }

    // 计算调色板的大小
    int paletteSize = *(int*)&bmpHeader[46];
    if (paletteSize == 0)
        paletteSize = 256;

    // 读取调色板数据
    uint8_t palette[1024];
    fread(palette, 1, paletteSize * 4, file);

    // 读取图像数据
    fseek(file, *(int*)&bmpHeader[10], SEEK_SET);
    fread(imageData, 1, *width * *height, file);

    fclose(file);

    return imageData;
}

// 将8位灰度图像数据保存为BMP文件
//filename：字符数组的指针，用于指定要保存的图像文件的名称或路径。
//imageData：无符号 8 位整型数据的指针，代表要保存的图像数据。
//width：图像的宽度。
//height：图像的高度。
void saveGrayScaleBMP(const char* filename, const uint8_t* imageData, int width, int height) 
{
    FILE* file = fopen(filename, "wb");
    if (!file) {
        fprintf(stderr, "Error creating file %s\n", filename);
        return;
    }

    // BMP文件头部信息
    uint8_t bmpHeader[54] = {
        0x42, 0x4D,             // 文件类型标识 "BM"
        0x36, 0x00, 0x0C, 0x00, // 文件大小（以字节为单位，此处假设图像数据大小不超过4GB）
        0x00, 0x00,             // 保留字段
        0x00, 0x00,             // 保留字段
        0x36, 0x00, 0x00, 0x00, // 位图数据偏移（以字节为单位）
        0x28, 0x00, 0x00, 0x00, // 位图信息头大小（40字节）
        0x00, 0x00, 0x00, 0x00, // 图像宽度
        0x00, 0x00, 0x00, 0x00, // 图像高度
        0x01, 0x00,             // 目标设备的级别（此处为1，不压缩）
        0x08, 0x00,             // 每个像素的位数（8位）
        0x00, 0x00, 0x00, 0x00, // 压缩类型（此处为不压缩）
        0x00, 0x00, 0x00, 0x00, // 图像数据大小（以字节为单位，此处为0，表示不压缩）
        0x00, 0x00, 0x00, 0x00, // 水平分辨率（像素/米，此处为0，表示未知）
        0x00, 0x00, 0x00, 0x00, // 垂直分辨率（像素/米，此处为0，表示未知）
        0x00, 0x00, 0x00, 0x00, // 使用的颜色索引数（0表示使用所有调色板项）
        0x00, 0x00, 0x00, 0x00  // 重要的颜色索引数（0表示所有颜色都重要）
    };

    // 更新BMP文件头部信息中的宽度和高度
    *(int*)&bmpHeader[18] = width;
    *(int*)&bmpHeader[22] = height;

    // 写入BMP文件头部信息
    fwrite(bmpHeader, 1, 54, file);

    // 写入调色板数据
    for (int i = 0; i < 256; i++) {
        fputc(i, file);  // 蓝色分量
        fputc(i, file);  // 绿色分量
        fputc(i, file);  // 红色分量
        fputc(0, file);  // 保留字节
    }

    // 写入图像数据
    fwrite(imageData, 1, width * height, file);

    fclose(file);
}

// 读取24位彩色图像的BMP文件
//filename：字符数组的指针，用于指定要读取的 BMP 格式图像文件的名称或路径。
//width：整型变量的指针，用于存储读取的图像的宽度。
//height：整型变量的指针，用于存储读取的图像的高度。
uint8_t* readColorBMP(const char* filename, int* width, int* height) 
{
    FILE* file = fopen(filename, "rb");
    if (!file) {
        fprintf(stderr, "Error opening file %s\n", filename);
        return NULL;
    }

    // 读取BMP文件头部信息
    uint8_t bmpHeader[54];
    fread(bmpHeader, 1, 54, file);

    // 从文件头部提取图像宽度和高度信息
    *width = *(int*)&bmpHeader[18];
    *height = *(int*)&bmpHeader[22];

    // 分配存储图像数据的内存
    uint8_t* imageData = (uint8_t*)malloc(*width * *height * 3);
    if (!imageData) {
        fprintf(stderr, "Memory allocation failed\n");
        fclose(file);
        return NULL;
    }

    // 读取图像数据
    fseek(file, *(int*)&bmpHeader[10], SEEK_SET);
    fread(imageData, 1, *width * *height * 3, file);

    fclose(file);

    return imageData;
}

//将24位彩色图像数据保存为BMP文件
//filename：字符数组的指针，用于指定要保存的图像文件的名称或路径。
//imageData：无符号 8 位整型数据的指针，代表要保存的图像数据。
//width：图像的宽度。
//height：图像的高度。
void saveColorBMP(const char* filename, const uint8_t* imageData, int width, int height) 
{
    FILE* file = fopen(filename, "wb");
    if (!file) {
        fprintf(stderr, "Error creating file %s\n", filename);
        return;
    }

    // BMP文件头部信息
    uint8_t bmpHeader[54] = {
        0x42, 0x4D,             // 文件类型标识 "BM"
        0x00, 0x00, 0x00, 0x00, // 文件大小（占位，稍后计算）
        0x00, 0x00,             // 保留字段
        0x00, 0x00,             // 保留字段
        0x36, 0x00, 0x00, 0x00, // 位图数据偏移（以字节为单位）
        0x28, 0x00, 0x00, 0x00, // 位图信息头大小（40字节）
        0x00, 0x00, 0x00, 0x00, // 图像宽度
        0x00, 0x00, 0x00, 0x00, // 图像高度
        0x01, 0x00,             // 目标设备的级别（此处为1，不压缩）
        0x18, 0x00,             // 每个像素的位数（24位）
        0x00, 0x00, 0x00, 0x00, // 压缩类型（此处为不压缩）
        0x00, 0x00, 0x00, 0x00, // 图像数据大小（占位，稍后计算）
        0x00, 0x00, 0x00, 0x00, // 水平分辨率（像素/米，此处为0，表示未知）
        0x00, 0x00, 0x00, 0x00, // 垂直分辨率（像素/米，此处为0，表示未知）
        0x00, 0x00, 0x00, 0x00, // 使用的颜色索引数（0表示使用所有调色板项）
        0x00, 0x00, 0x00, 0x00  // 重要的颜色索引数（0表示所有颜色都重要）
    };

    // 更新BMP文件头部信息中的宽度和高度
    *(int*)&bmpHeader[18] = width;
    *(int*)&bmpHeader[22] = height;

    // 计算图像数据大小
    uint32_t imageDataSize = width * height * 3 + 54; // 加上文件头部大小
    bmpHeader[2] = (uint8_t)(imageDataSize & 0xFF);
    bmpHeader[3] = (uint8_t)((imageDataSize >> 8) & 0xFF);
    bmpHeader[4] = (uint8_t)((imageDataSize >> 16) & 0xFF);
    bmpHeader[5] = (uint8_t)((imageDataSize >> 24) & 0xFF);

    // 写入BMP文件头部信息
    fwrite(bmpHeader, 1, 54, file);

    // 写入图像数据
    fwrite(imageData, width * height * 3, 1, file);

    fclose(file);
}

//24位彩色图像转8位灰度值
//rgbImage原始图像
//grayImage输出灰度图像
//width,height图片的宽和高
void convertToGray(uint8_t* rgbImage, uint8_t* grayImage, int width, int height)
{
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            // 获取当前像素的 RGB 分量
            uint8_t r = rgbImage[3 * (y * width + x) + 0];
            uint8_t g = rgbImage[3 * (y * width + x) + 1];
            uint8_t b = rgbImage[3 * (y * width + x) + 2];

            // 计算灰度值（常用的加权平均法）
            // 这里使用的加权系数是常见的：R: 0.299, G: 0.587, B: 0.114
            uint8_t gray = (uint8_t)(0.299 * r + 0.587 * g + 0.114 * b);

            // 将灰度值写入灰度图像数组
            grayImage[y * width + x] = gray;
        }
    }
}

//灰度线性拉伸
//pGryImg：灰度图像数据的指针。
//width：图像的宽度。
//height：图像的高度。
//k：线性拉伸的斜率。它控制着拉伸的速率或程度。当(k) 大于 1 时，图像的对比度增加；当(k) 小于 1 时，对比度降低。
//b：线性拉伸的偏移。它控制着拉伸后灰度值的起始位置。当(b) 大于 0 时，图像的整体亮度增加；当(b) 小于 0 时，整体亮度减小。
void LinearStretchDemo(uint8_t* pGryImg, int width, int height, double k, double b)
{
    uint8_t* pCur, * pEnd;
    int LUT[256];    //因为只有[0,255]共256个灰度值

    //step1. 生成查找表
    for (int g = 0; g < 256; g++)
    {
        LUT[g] = max(0, min(255, k * g + b));
    }

    //step2. 进行变换
    for (pCur = pGryImg, pEnd = pGryImg + width * height; pCur < pEnd; pCur++)
    {
        *pCur = LUT[*pCur];
    }
    //step3. 结束
    return;
}

//统计图像灰度值
//pImg：灰度图像数据的指针。
//width：图像的宽度。
//height：图像的高度。
//* histogram：数组首元素地址，需要一个能储存256个变量的整型数组
void GetHistogram(uint8_t* pImg, int width, int height, int* histogram)
{
    uint8_t* pCur;
    uint8_t* pEnd = pImg + width * height;

    // 初始化直方图数组
    memset(histogram, 0, sizeof(int) * 256);

    // 直方图统计
    for (pCur = pImg; pCur < pEnd;)
    {
        histogram[*pCur]++;
        pCur++;
    }

    // 函数结束
    return;
}

//亮度和对比度
//储存histogram灰度直方图的指针
//接收亮度的变量地址
//接收对比度的变量地址
void GetBrightContrast(int* histogram, double* bright, double* contrast)
{
    int g;
    double sum, num; //书上说图像很亮时，int有可能会溢出，所以我这里直接用double
    double fsum;

    //step.1 求亮度
    for (sum = num = 0, g = 0; g < 256; g++)
    {
        sum += histogram[g] * g;
        num += histogram[g];
    }
    *bright = sum * 1.0 / num;

    //step.2 求对比度
    for (fsum = 0.0, g = 0; g < 256; g++)
    {
        fsum += histogram[g] * (g - *bright) * (g - *bright);
    }
    *contrast = sqrt(fsum / (num - 1)); //即Std Dev

    //step.3 结束
    return;
}

//pGryImg：灰度图像数据的指针。
//width：图像的宽度。
//height：图像的高度。
void RmwHistogramEqualize(uint8_t* pGryImg, int width, int height)
{
    uint8_t* pCur, * pEnd = pGryImg + width * height; // 指针变量，指向当前像素和图像末尾
    int histogram[256], LUT[256], A, g; // 直方图数组、查找表数组、累积直方图、灰度级

    // step.1-------------求直方图--------------------------//
    memset(histogram, 0, sizeof(int) * 256); // 初始化直方图数组为0
    for (pCur = pGryImg; pCur < pEnd;)
        histogram[*(pCur++)]++; // 统计每个灰度级出现的频率

    // step.2-------------求LUT[g]-------------------------//
    A = histogram[0]; // 初始化累积直方图的值为第一个灰度级的频率
    LUT[0] = 255 * A / (width * height); // 计算第一个灰度级对应的均衡化后的灰度值
    for (g = 1; g < 256; g++) {
        A += histogram[g]; // 更新累积直方图的值
        LUT[g] = 255 * A / (width * height); // 计算当前灰度级对应的均衡化后的灰度值
    }

    // step.3-------------查表------------------------------//
    for (pCur = pGryImg; pCur < pEnd;)
        *(pCur++) = LUT[*pCur]; // 使用查找表对每个像素进行灰度映射

    // step.4-------------结束------------------------------//
    return;
}

//对数变换
//pGryImg：灰度图像数据的指针。
//width：图像的宽度。
//height：图像的高度。
void RmwLogTransform(uint8_t* pGryImg, int width, int height)
{
    uint8_t* pCur, * pEnd = pGryImg + width * height; // 指向灰度图像数据的当前指针和结束指针
    int histogram[256], LUT[256], gmax, g; // 声明直方图数组、查找表数组、最大灰度值、当前灰度值
    double c; // 声明常数c

    // step.1-------------求直方图--------------------------//
    memset(histogram, 0, sizeof(int) * 256); // 初始化直方图数组为0
    for (pCur = pGryImg; pCur < pEnd;)
        histogram[*(pCur++)]++; // 遍历图像数据，统计每个灰度级的像素数量

    // step.2-------------最大值---------------------------//
    for (gmax = 255; gmax >= 0; gmax++)
        if (histogram[gmax]) break; // 从最大灰度级开始向低灰度级搜索，找到第一个非零灰度级，即最大灰度值

    // step.3-------------求LUT[g]-------------------------//
    c = 255.0 / log(1 + gmax); // 计算常数c
    for (g = 0; g < 256; g++)
    {
        LUT[g] = (int)(c * log(1 + g)); // 根据对数变换公式计算查找表中每个灰度级的映射值
    }

    // step.4-------------查表------------------------------//
    for (pCur = pGryImg; pCur < pEnd;)
        *(pCur++) = LUT[*pCur]; // 使用查找表将图像数据进行对数变换

    // step.5-------------结束------------------------------//
    return; // 函数结束
}

//基于列积分的快速均值滤波
//原始灰度图像
//图像的宽度和高度
//滤波邻域：M列N行
//结果图像
void RmwAvrFilterBySumCol(uint8_t* pGryImg,int width, int height,int M, int N,uint8_t* pResImg) 
{
    uint8_t* pAdd, * pDel, * pRes;
    int halfx, halfy;
    int x, y;
    int sum, c;
    int sumCol[4096]; // 约定图像宽度不大于4096

    // step.1------------初始化--------------------------//
    M = M / 2 * 2 + 1; // 奇数化
    N = N / 2 * 2 + 1; // 奇数化
    halfx = M / 2; // 滤波器的半径x
    halfy = N / 2; // 滤波器的半径y
    c = (1 << 23) / (M * N); // 乘法因子
    memset(sumCol, 0, sizeof(int) * width);
    for (y = 0, pAdd = pGryImg; y < N; y++) {
        for (x = 0; x < width; x++) sumCol[x] += *(pAdd++);
    }
    // step.2------------滤波----------------------------//
    for (y = halfy, pRes = pResImg + y * width, pDel = pGryImg; y < height - halfy; y++) {
        // 初值
        for (sum = 0, x = 0; x < M; x++) sum += sumCol[x];
        // 滤波
        pRes += halfx; // 跳过左侧
        for (x = halfx; x < width - halfx; x++) {
            // 求灰度均值
            // *(pRes++)=sum/(N*M);
            *(pRes++) = (sum * c) >> 23; // 用整数乘法和移位代替除法
            // 换列,更新灰度和
            sum -= sumCol[x - halfx]; // 减左边列
            sum += sumCol[x + halfx + 1]; // 加右边列
        }
        pRes += halfx; // 跳过右侧
        // 换行,更新sumCol
        for (x = 0; x < width; x++) {
            sumCol[x] -= *(pDel++); // 减上一行
            sumCol[x] += *(pAdd++); // 加下一行
        }
    }
    // step.3------------返回----------------------------//
    return;
}

//基于列积分的积分图实现
//pGryImg, // 原始灰度图像
//width,       // 图像的宽度 
//height,      // 图像的高度
//pSumImg     // 计算得到的积分图
void RmwDoSumGryImg(uint8_t* pGryImg,int width,int height, int* pSumImg)
{
    uint8_t* pGry;
    int* pRes;
    int x, y;
    int sumCol[4096]; // 约定图像宽度不大于4096

    memset(sumCol, 0, sizeof(int) * width);
    for (y = 0, pGry = pGryImg, pRes = pSumImg; y < height; y++)
    {
        // 最左侧像素的特别处理
        sumCol[0] += *(pGry++);
        *(pRes++) = sumCol[0];
        // 正常处理
        for (x = 1; x < width; x++)
        {
            sumCol[x] += *(pGry++);       // 更新列积分
            int temp = *(pRes - 1);
            *(pRes++) = temp + sumCol[x];
        }
    }
    return;
}

//基于SSE的积分图实现
//pGryImg原始灰度图像
//width图像的宽度，必须是4的倍数
//height图像的高度
//pSumImg计算得到的积分图
void RmwDoSumGryImg_SSE(uint8_t* pGryImg,int width,int height,int* pSumImg)
{
    int sumCol[4096]; //约定图像宽度不大于4096
    __m128i* pSumSSE, A;
    uint8_t* pGry;
    int* pRes;
    int x, y;

    memset(sumCol, 0, sizeof(int) * width);
    for (y = 0, pGry = pGryImg, pRes = pSumImg; y < height; y++)
    {
        // 0:需要特别处理
        sumCol[0] += *(pGry++);
        *(pRes++) = sumCol[0];
        // 1
        sumCol[1] += *(pGry++);
        *(pRes++) = *(pRes - 1) + sumCol[1];
        // 2
        sumCol[2] += *(pGry++);
        *(pRes++) = *(pRes - 1) + sumCol[2];
        // 3
        sumCol[3] += *(pGry++);
        *(pRes++) = *(pRes - 1) + sumCol[3];
        // [4...width-1]
        for (x = 4, pSumSSE = (__m128i*)(sumCol + 4); x < width; x += 4, pGry += 4)
        {
            // 把变量的低32位(有4个8位整数组成)转换成32位的整数
            A = _mm_cvtepu8_epi32(_mm_loadl_epi64((__m128i*)pGry));
            // 4个32位的整数相加
            *(pSumSSE++) = _mm_add_epi32(*pSumSSE, A);
            // 递推
            *(pRes++) = *(pRes - 1) + sumCol[x + 0];
            *(pRes++) = *(pRes - 1) + sumCol[x + 1];
            *(pRes++) = *(pRes - 1) + sumCol[x + 2];
            *(pRes++) = *(pRes - 1) + sumCol[x + 3];
        }
    }
    return;
}

//基于积分图的快速均值滤波
//pSumImg计算得到的积分图
//width,height,图像的宽度和高度
//M, N,滤波邻域：M列N行
//pResImg 结果图像
void RmwAvrFilterBySumImg(int* pSumImg,int width, int height,int M, int N,uint8_t* pResImg)
{
    // 没有对边界上邻域不完整的像素进行处理，可以采用变窗口的策略
    int* pY1, * pY2;
    uint8_t* pRes;
    int halfx, halfy;
    int y, x1, x2;
    int sum, c;

    // step.1------------初始化--------------------------//
    M = M / 2 * 2 + 1; // 奇数化
    N = N / 2 * 2 + 1; // 奇数化
    halfx = M / 2;      // 滤波器的半径x
    halfy = N / 2;      // 滤波器的半径y
    c = (1 << 23) / (M * N); // 乘法因子
    // step.2------------滤波----------------------------//
    for (y = halfy + 1, pRes = pResImg + y * width, pY1 = pSumImg, pY2 = pSumImg + N * width;
        y < height - halfy;
        y++, pY1 += width, pY2 += width)
    {
        pRes += halfx + 1; // 跳过左侧
        for (x1 = 0, x2 = M; x2 < width; x1++, x2++) // 可以简化如此，但不太容易读
        {
            sum = *(pY2 + x2) - *(pY2 + x1) - *(pY1 + x2) + *(pY1 + x1);
            *(pRes++) = (uint8_t)((sum * c) >> 23); // 用整数乘法和移位代替除法
        }
        pRes += halfx; // 跳过右侧
    }
    // step.3------------返回----------------------------//
    return;
}

void GetMedianGry(int* histogram, int N, int* medGry)
{
    int g;
    int num;

    // step.1-------------求灰度中值------------------------//
    num = 0;
    for (g = 0; g < 256; g++)
    {
        num += histogram[g];
        if (2 * num > N) break;  //num>N/2
    }
    *medGry = g;
    // step.2-------------结束------------------------------//
    return;
}

//中值滤波
//pGryImg：指向待处理灰度图像数据的指针。
//width、height：表示图像的宽度和高度。
//M、N：分别表示中值滤波器的水平和垂直邻域大小（以像素为单位）。
//pResImg：指向存储结果图像数据的指针。
double RmwMedianFilter(uint8_t* pGryImg, int width, int height, int M, int N, uint8_t* pResImg) 
{
    uint8_t* pCur, * pRes;
    int halfx, halfy, x, y, i, j, y1, y2;
    int histogram[256];
    int wSize, j1, j2;
    int num, med, v;
    int dbgCmpTimes = 0; // 搜索中值所需比较次数的调试

    M = M / 2 * 2 + 1; // 奇数化
    N = N / 2 * 2 + 1; // 奇数化
    halfx = M / 2;      // x半径
    halfy = N / 2;      // y半径
    wSize = (halfx * 2 + 1) * (halfy * 2 + 1); // 邻域内像素总个数

    for (y = halfy, pRes = pResImg + y * width; y < height - halfy; y++) {
        // step.1----初始化直方图
        y1 = y - halfy;
        y2 = y + halfy;
        memset(histogram, 0, sizeof(int) * 256);

        for (i = y1, pCur = pGryImg + i * width; i <= y2; i++, pCur += width) {
            for (j = 0; j < halfx * 2 + 1; j++) {
                histogram[*(pCur + j)]++;
            }
        }

        // step.2-----初始化中值
        num = 0; // 记录着灰度值从0到中值的个数
        for (i = 0; i < 256; i++) {
            num += histogram[i];
            if (num * 2 > wSize) {
                med = i;
                break;
            }
        }

        // 滤波
        pRes += halfx; // 没有处理图像左边界侧的像素
        for (x = halfx; x < width - halfx; x++) {
            // 赋值
            *(pRes++) = med;

            // step.3-----直方图递推: 减去当前邻域最左边的一列,添加邻域右侧的一个新列
            j1 = x - halfx;     // 最左边列
            j2 = x + halfx + 1; // 右边的新列

            for (i = y1, pCur = pGryImg + i * width; i <= y2; i++, pCur += width) {
                // 减去最左边列
                v = *(pCur + j1);
                histogram[v]--;  // 更新直方图
                if (v <= med) num--; // 更新num

                // 添加右边的新列
                v = *(pCur + j2);
                histogram[v]++; // 更新直方图
                if (v <= med) num++; // 更新num
            }

            // step.4-----更新中值
            if (num * 2 < wSize) { // 到上次中值med的个数不够了,则med要变大
                for (med = med + 1; med < 256; med++) {
                    dbgCmpTimes += 2; // 总的比较次数,调试用
                    num += histogram[med];
                    if (num * 2 > wSize) break;
                }
                dbgCmpTimes += 1; // 总的比较次数,调试用
            }
            else { // 到上次中值med的个数多了,则med要变小
                while ((num - histogram[med]) * 2 > wSize) { // 若减去后,仍变小
                    dbgCmpTimes++; // 总的比较次数,调试用
                    num -= histogram[med];
                    med--;
                }
                dbgCmpTimes += 2; // 总的比较次数,调试用
            }
        }
        pRes += halfx; // 没有处理图像右边界侧的像素
    }
    // 返回搜索中值需要的平均比较次数
    return dbgCmpTimes * 1.0 / ((width - halfx * 2) * (height - halfy * 2));
}

//二值滤波
//pBinImg,  原始二值图像
// width, height,图像的宽度和高度
// M, N, 滤波邻域：M列N行
// threshold, 灰度阈值,大于等于该值时结果赋255
// pResImg 结果图像
void RmwBinImgFilter(uint8_t* pBinImg,int width, int height,int M, int N,double threshold,uint8_t* pResImg )
{
    // 没有对边界上邻域不完整的像素进行处理，可以采用变窗口的策略
    uint8_t* pAdd, * pDel, * pRes;
    int halfx, halfy;
    int x, y, sum, sumThreshold;
    int sumCol[4096]; //约定图像宽度不大于4096

    // step.1------------初始化--------------------------//
    M = M / 2 * 2 + 1; //奇数化
    N = N / 2 * 2 + 1; //奇数化
    halfx = M / 2; //滤波器的x半径
    halfy = N / 2; //滤波器的y半径
    sumThreshold = max(1, (int)(threshold * M * N)); //转换成邻域内灰度值之和的阈值
    memset(sumCol, 0, sizeof(int) * width);
    for (y = 0, pAdd = pBinImg; y < N; y++)
    {
        for (x = 0; x < width; x++)
            sumCol[x] += *(pAdd++);
    }
    // step.2------------滤波----------------------------//
    for (y = halfy, pRes = pResImg + y * width, pDel = pBinImg; y < height - halfy; y++)
    {
        //初值
        for (sum = 0, x = 0; x < M; x++)
            sum += sumCol[x];
        //滤波
        pRes += halfx; //跳过左侧
        for (x = halfx; x < width - halfx; x++)
        {
            //求灰度均值
            /*if (sum>=sumThreshold)
            {
                *(pRes++) = 255;
            }
            else  *(pRes++) = 0;*/
            *(pRes++) = (sum >= sumThreshold) * 255; //请理解这个表达式的含义
            //换列,更新灰度和
            sum -= sumCol[x - halfx];     //减左边列
            sum += sumCol[x + halfx + 1]; //加右边列
        }
        pRes += halfx; //跳过右侧
        //换行,更新sumCol
        for (x = 0; x < width; x++)
        {
            sumCol[x] -= *(pDel++); //减上一行
            sumCol[x] += *(pAdd++); //加下一行
        }
    }
    // step.3------------返回----------------------------//
    return;
}

//梯度算子加阈值
//pGryImg：输入的灰度图像数据指针。
//width：图像的宽度。
//height：图像的高度。
//pGrdImg：输出的梯度图像数据指针
void RmwGradientGryImg(uint8_t* pGryImg, int width, int height, uint8_t* pGrdImg)
{
    uint8_t* pGry, * pGrd;
    int dx, dy;
    int x, y;

    for (y = 0, pGry = pGryImg, pGrd = pGrdImg; y < height - 1; y++)
    {
        for (x = 0; x < width - 1; x++, pGry++)
        {
            dx = *pGry - *(pGry + 1);
            dy = *pGry - *(pGry + width);
            int gradient = (int)(sqrt(dx * dx * 1.0 + dy * dy));
            *pGrd++ = (gradient > 255) ? 255 : gradient;
        }
        *pGrd++ = 0; //尾列不做,边缘强度赋0
        pGry++;
    }
    memset(pGrd, 0, width); //尾行不做,边缘强度赋0
}

//梯度算子加阈值
//pGryImg：输入的灰度图像数据指针。
//width：图像的宽度。
//height：图像的高度。
//pGrdImg：输出的梯度图像数据指针
void RmwGradientGryImgPlus(uint8_t* pGryImg, int width, int height, uint8_t* pGrdImg, int threshold)
{
    uint8_t* pGry, * pGrd;
    int dx, dy;
    int x, y;

    for (y = 0, pGry = pGryImg, pGrd = pGrdImg; y < height - 1; y++)
    {
        for (x = 0; x < width - 1; x++, pGry++)
        {
            dx = *pGry - *(pGry + 1);
            dy = *pGry - *(pGry + width);
            int gradient = (int)(sqrt(dx * dx * 1.0 + dy * dy));
            *(pGrd++) = (gradient > threshold) ? min(255, gradient) : 0;
        }
        *(pGrd++) = 0; //尾列不做,边缘强度赋0
        pGry++;
    }
    memset(pGrd, 0, width); //尾行不做,边缘强度赋0
}

//反相
void invertImage(uint8_t* image, int width, int height) {
    for (int i = 0; i < width * height; i++) {
        image[i] = 255 - image[i];
    }
}

//罗伯特算子
void RmwRobertsGryImg(uint8_t* pGryImg, int width, int height, uint8_t* pRbtImg)
{
    uint8_t* pGry, * pRbt;
    int dx, dy;
    int x, y;

    for (y = 0, pGry = pGryImg, pRbt = pRbtImg; y < height - 1; y++)
    {
        for (x = 0; x < width - 1; x++, pGry++)
        {
            dx = *pGry - *(pGry + width + 1);
            dy = *(pGry + 1) - *(pGry + width);
            *pRbt++ = (uint8_t)(dx > dy ? dx : dy); // 使用三目运算符选择较大的值
        }
        *pRbt++ = 0; // 尾列不做, 边缘强度赋0
        pGry++;
    }
    memset(pRbt, 0, width); // 尾行不做, 边缘强度赋0
}

//索贝尔算子
void RmwSobelGryImg(uint8_t* pGryImg, int width, int height, uint8_t* pSbImg)
{
    uint8_t* pGry, * pSb;
    int dx, dy;
    int x, y;

    memset(pSbImg, 0, width); // 首行不做, 边缘强度赋0
    for (y = 1, pGry = pGryImg + width, pSb = pSbImg + width; y < height - 1; y++)
    {
        *pSb++ = 0; // 首列不做, 边缘强度赋0
        pGry++;
        for (x = 1; x < width - 1; x++, pGry++)
        {
            // 求dx
            dx = *(pGry - 1 - width) + (*(pGry - 1) * 2) + *(pGry - 1 + width);
            dx -= *(pGry + 1 - width) + (*(pGry + 1) * 2) + *(pGry + 1 + width);
            // 求dy
            dy = *(pGry - width - 1) + (*(pGry - width) * 2) + *(pGry - width + 1);
            dy -= *(pGry + width - 1) + (*(pGry + width) * 2) + *(pGry + width + 1);
            // 结果
            *pSb++ = (uint8_t)min(255, abs(dx) + abs(dy));
        }
        *pSb++ = 0; // 尾列不做, 边缘强度赋0
        pGry++;
    }
    memset(pSb, 0, width); // 尾行不做, 边缘强度赋0
}

//Prewitt算子
void RmwPrewittGryImg(uint8_t* pGryImg, int width, int height, uint8_t* pPRTImg)
{
    uint8_t* pGry, * pPRT;
    int dx, dy, d45, d135, v1, v2;
    int x, y;

    memset(pPRTImg, 0, width); // 首行不做, 边缘强度赋0
    for (y = 1, pGry = pGryImg + width, pPRT = pPRTImg + width; y < height - 1; y++)
    {
        *pPRT++ = 0; // 首列不做, 边缘强度赋0
        pGry++;
        for (x = 1; x < width - 1; x++, pGry++)
        {
            // 求dx
            dx = *(pGry - 1 - width) + *(pGry - 1) + *(pGry - 1 + width);
            dx -= *(pGry + 1 - width) + *(pGry + 1) + *(pGry + 1 + width);
            // 求dy
            dy = *(pGry - width - 1) + *(pGry - width) + *(pGry - width + 1);
            dy -= *(pGry + width - 1) + *(pGry + width) + *(pGry + width + 1);
            // 求45度
            d45 = *(pGry - width - 1) + *(pGry - width) + *(pGry - 1);
            d45 -= *(pGry + width + 1) + *(pGry + width) + *(pGry + 1);
            // 求135度
            d135 = *(pGry - width) + *(pGry - width + 1) + *(pGry + 1);
            d135 -= *(pGry + width - 1) + *(pGry + width) + *(pGry - 1);
            // 结果
            v1 = abs(dx) > abs(dy) ? abs(dx) : abs(dy);
            v2 = abs(d45) > abs(d135) ? abs(d45) : abs(d135);
            *pPRT++ = (uint8_t)((v1 > v2) ? ((v1 > 255) ? 255 : v1) : ((v2 > 255) ? 255 : v2));
        }
        *pPRT++ = 0; // 尾列不做, 边缘强度赋0
        pGry++;
    }
    memset(pPRT, 0, width); // 尾行不做, 边缘强度赋0
}

//沈俊算子
//pGryImg 和 pTmpImg 是指向 uint8_t 类型的指针，它们分别指向原始灰度图像数据和辅助图像数据。
//width 和 height 是整型参数，表示图像的宽度和高度。
//a0 是双精度浮点型参数，表示滤波系数。
//pSJImg 是指向 uint8_t 类型的指针，它指向了输出的图像数据。
void RmwShenJunGryImg(uint8_t* pGryImg,uint8_t* pTmpImg, int width, int height, double a0, uint8_t* pSJImg)
{
    uint8_t* pGry, * pCur, * pSJ, * pEnd;
    int LUT[512], * ALUT; // a0查找表
    int x, y, pre, dif;

    // Step 1: 初始化查找表
    a0 = (a0 < 0.01) ? 0.01 : ((a0 > 0.99) ? 0.99 : a0); // 安全性检查
    // a0查找表, 进行了四舍五入
    ALUT = LUT + 256;
    for (ALUT[0] = 0, dif = 1; dif < 256; dif++)
    {
        ALUT[dif] = (int)(dif * a0 + 0.5);
        ALUT[-dif] = (int)(-dif * a0 - 0.5);
    }

    // Step 2: 递推实现指数滤波
    // 按行滤波
    for (y = 0, pGry = pGryImg, pCur = pTmpImg; y < height; y++)
    {
        // 从左向右: p1(y,x) = p1(y,x-1) + a * [p(y,x) - p1(y,x-1)]
        *(pCur++) = pre = *(pGry++);
        for (x = 1; x < width; x++, pGry++)
            *(pCur++) = pre = pre + ALUT[*pGry - pre];
        pCur--; // 回到行尾
        // 从右向左: p2(y,x) = p2(y,x+1) - a * [p1(y,x) - p2(y,x+1)]
        for (x = width - 2, pCur = pCur - 1; x >= 0; x--)
            *(pCur--) = pre = pre + ALUT[*pCur - pre];
        pCur += (width + 1); // 回到下一行的开始
    }
    // 按列滤波
    for (x = 0, pCur = pTmpImg; x < width; x++, pCur = pTmpImg + x)
    {
        // 从上向下: p3(y,x) = p3(y-1,x) + a * [p2(y,x) - p3(y-1,x)]
        pre = *pCur;
        for (y = 1, pCur += width; y < height; y++, pCur += width)
            *pCur = pre = pre + ALUT[*pCur - pre];
        pCur -= width; // 回到列尾
        // 从下向上: p4(i,j) = p4(i+1,j) + a * [p3(i,j) - p4(i+1,j)]
        for (y = height - 2, pCur -= width; y >= 0; y--, pCur -= width)
            *pCur = pre = pre + ALUT[*pCur - pre];
    }

    // Step 3: 正导数=1，负导数为0，0必须也是0
    pEnd = pTmpImg + width * height;
    for (pCur = pTmpImg, pGry = pGryImg; pCur < pEnd; pGry++)
    {
        *(pCur++) = (*pCur > *pGry);
    }

    // Step 4: 过零点检测
    memset(pSJImg, 0, width * height); // 边缘强度赋0
    pSJ = pSJImg + width;
    pCur = pTmpImg + width; // 首行不做 
    for (y = 1; y < height - 1; y++)
    {
        pSJ++; pCur++;  // 首列不做
        for (x = 1; x < width - 1; x++, pGry++, pCur++, pSJ++)
        {
            if (*pCur) // 正导数
            {
                // 下面使用4邻域, 边缘为8连通, 不能保证4连通; 使用8邻域才能保证边缘4连通
                if ((!*(pCur - 1)) || // 左, 必须<=0, 不能<0
                    (!*(pCur + 1)) || // 右, 必须<=0, 不能<0
                    (!*(pCur - width)) || // 上, 必须<=0, 不能<0
                    (!*(pCur + width)))   // 下, 必须<=0, 不能<0
                {
                    *pSJ = 255; // 周围有导数小于等于0
                }
            }
        }
        pSJ++; pCur++;  // 尾列不做
    }
}

//沈俊算子加索贝尔算子
//pGryImg：指向原始灰度图像数据的指针
//pTmpImg：指向辅助图像数据的指针
//width：图像的宽度
//height：图像的高度
//a0：这是沈俊算子的参数，用于控制边缘检测的灵敏度。
//grdThre：这是Sobel算子的梯度阈值
//pEdgeImg：最终边缘图像数据的指针
void RmwExtractRiceEdge(uint8_t* pGryImg,uint8_t* pTmpImg,int width,int height,double a0, int grdThre, uint8_t* pEdgeImg)
{
    // step.1------------沈俊算子-----------------------//
    RmwShenJunGryImg(pGryImg, pTmpImg, width, height, a0, pEdgeImg);
    // step.2------------Sobel算子----------------------//
    RmwSobelGryImg(pGryImg, width, height, pTmpImg);
    // step.3------------二者融合-----------------------//
    for (int i = 0; i < width * height; i++)
    {
        *(pEdgeImg + i) = (pEdgeImg[i] && (pTmpImg[i] > grdThre)) * 255;
    }
    // step.4------------结束---------------------------//
    return;
}