深度剖析图像处理—边缘检测

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

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

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

先简单写一下反相的函数

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

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

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

//梯度算子加阈值
//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
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

看一下效果

画面确实变得更干净了

罗伯特算子

梯度算子的计算只涉及到了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
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

我们来看结果

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

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

我们来看下一个算子

索贝尔算子

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

因此，索贝尔算子(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
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

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

以下从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
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

我们看一下实现效果

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 对每个像素𝑆𝐽(𝑥,𝑦)进行过零点检测得到边缘点。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

我们来看程序

//沈俊算子
//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++;  // 尾列不做
    }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

当滤波系数为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;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

我们换一张图片试一下

总结

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

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

源码

IDP.h

#pragma once

#include 
#include 
#include 
#include 
#include 
 
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);//索贝尔＋沈俊算子
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

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;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931