c++/matlab编程实现 canny 算子进行边缘检测。

C++代码实现：

#include <iostream>
#include <cmath>
#include <algorithm>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

// 高斯滤波
void gaussianBlur(Mat &image, int ksize, double sigma) {
    Mat kernel = getGaussianKernel(ksize, sigma, CV_64F);
    sepFilter2D(image, image, -1, kernel, kernel);
}

// sobel 算子计算水平方向和垂直方向梯度
void sobel(Mat &image, Mat &grad_x, Mat &grad_y) {
    Sobel(image, grad_x, CV_32F, 1, 0, 3);
    Sobel(image, grad_y, CV_32F, 0, 1, 3);
}

// 计算梯度幅值和方向
void gradient(Mat &grad_x, Mat &grad_y, Mat &grad_mag, Mat &grad_dir) {
    cartToPolar(grad_x, grad_y, grad_mag, grad_dir, true);
}

// 非最大抑制
void nonMaximumSuppression(Mat &grad_mag, Mat &grad_dir, Mat &grad_nms) {
    float pi = 3.14159265358979323846;
    grad_nms = Mat::zeros(grad_mag.size(), grad_mag.type());
    for (int i = 1; i < grad_mag.rows - 1; i++) {
        for (int j = 1; j < grad_mag.cols - 1; j++) {
            float angle = grad_dir.at<float>(i, j);
            float m1, m2;
            // 比较梯度方向上的两个像素
            if (angle < pi / 4 && angle >= -pi / 4) {
                m1 = grad_mag.at<float>(i, j - 1);
                m2 = grad_mag.at<float>(i, j + 1);
            } else if (angle < -pi / 4 && angle >= -3 * pi / 4) {
                m1 = grad_mag.at<float>(i - 1, j);
                m2 = grad_mag.at<float>(i + 1, j);
            } else if (angle < 3 * pi / 4 && angle >= pi / 4) {
                m1 = grad_mag.at<float>(i - 1, j - 1);
                m2 = grad_mag.at<float>(i + 1, j + 1);
            } else {
                m1 = grad_mag.at<float>(i - 1, j + 1);
                m2 = grad_mag.at<float>(i + 1, j - 1);
            }
            // 如果当前像素的梯度值不是最大值，就将其设为 0
            if (grad_mag.at<float>(i, j) < m1 || grad_mag.at<float>(i, j) < m2) {
                grad_nms.at<float>(i, j) = 0;
            } else {
                grad_nms.at<float>(i, j) = grad_mag.at<float>(i, j);
            }
        }
    }
}

// 双阈值处理
void doubleThreshold(Mat &grad_nms, Mat &grad_thres, float low_threshold, float high_threshold) {
    grad_thres = Mat::zeros(grad_nms.size(), grad_nms.type());
    for (int i = 0; i < grad_nms.rows; i++) {
        for (int j = 0; j < grad_nms.cols; j++) {
            float val = grad_nms.at<float>(i, j);
            if (val > high_threshold) {
                grad_thres.at<float>(i, j) = 255;
            } else if (val > low_threshold) {
                grad_thres.at<float>(i, j) = 127;
            }
        }
    }
}

// 连通域分析
void connectedComponents(Mat &grad_thres, Mat &edges, int min_size) {
    edges = Mat::zeros(grad_thres.size(), grad_thres.type());
    vector<vector<Point>> contours;
    findContours(grad_thres, contours, RETR_EXTERNAL, CHAIN_APPROX_SIMPLE);
    for (int i = 0; i < contours.size(); i++) {
        if (contours[i].size() < min_size) {
            continue;
        }
        for (int j = 0; j < contours[i].size(); j++) {
            Point p = contours[i][j];
            edges.at<float>(p) = 255;
        }
    }
}

int main() {
    Mat image = imread("test.jpg", IMREAD_GRAYSCALE);
    if (image.empty()) {
        cerr << "Failed to open image file!" << endl;
        return -1;
    }

    // 高斯滤波
    Mat blurred;
    gaussianBlur(image, 5, 1.4);

    // sobel 算子计算梯度
    Mat grad_x, grad_y;
    sobel(blurred, grad_x, grad_y);

    // 计算梯度幅值和方向
    Mat grad_mag, grad_dir;
    gradient(grad_x, grad_y, grad_mag, grad_dir);

    // 非最大抑制
    Mat grad_nms;
    nonMaximumSuppression(grad_mag, grad_dir, grad_nms);

    // 双阈值处理
    float low_threshold = 35, high_threshold = 70;
    Mat grad_thres;
    doubleThreshold(grad_nms, grad_thres, low_threshold, high_threshold);

    // 连通域分析
    int min_size = 10;
    Mat edges;
    connectedComponents(grad_thres, edges, min_size);

    // 显示结果
    namedWindow("Original", WINDOW_AUTOSIZE);
    imshow("Original", image);
    namedWindow("Edges", WINDOW_AUTOSIZE);
    imshow("Edges", edges);
    waitKey(0);

    return 0;
}

Matlab代码实现：

function edges = canny(image, low_threshold, high_threshold, min_size)
    % 高斯滤波
    blurred = imgaussfilt(image, 1.4);
    
    % sobel 算子计算梯度
    [grad_x, grad_y] = gradient(double(blurred));
    
    % 计算梯度幅值和方向
    grad_mag = hypot(grad_x, grad_y);
    grad_dir = atan2(grad_y, grad_x);
    
    % 非最大抑制
    grad_nms = non_maximum_suppression(grad_mag, grad_dir);
    
    % 双阈值处理
    grad_thres = double_threshold(grad_nms, low_threshold, high_threshold);
    
    % 连通域分析
    edges = connected_components(grad_thres, min_size);
end

function grad_nms = non_maximum_suppression(grad_mag, grad_dir)
    [rows, cols] = size(grad_mag);
    grad_nms = zeros(rows, cols);
    pi = 3.14159265358979323846;
    for i = 2 : rows - 1
        for j = 2 : cols - 1
            angle = grad_dir(i, j);
            if angle < pi / 4 && angle >= -pi / 4
                m1 = grad_mag(i, j - 1);
                m2 = grad_mag(i, j + 1);
            elseif angle < -pi / 4 && angle >= -3 * pi / 4
                m1 = grad_mag(i - 1, j);
                m2 = grad_mag(i + 1, j);
            elseif angle < 3 * pi / 4 && angle >= pi / 4
                m1 = grad_mag(i - 1, j - 1);
                m2 = grad_mag(i + 1, j + 1);
            else
                m1 = grad_mag(i - 1, j + 1);
                m2 = grad_mag(i + 1, j - 1);
            end
            if grad_mag(i, j) < m1 || grad_mag(i, j) < m2
                grad_nms(i, j) = 0;
            else
                grad_nms(i, j) = grad_mag(i, j);
            end
        end
    end
end

function grad_thres = double_threshold(grad_nms, low_threshold, high_threshold)
    [rows, cols] = size(grad_nms);
    grad_thres = zeros(rows, cols);
    for i = 1 : rows
        for j = 1 : cols
            val = grad_nms(i, j);
            if val > high_threshold
                grad_thres(i, j) = 255;
            elseif val > low_threshold
                grad_thres(i, j) = 127;
            end
        end
    end
end

function edges = connected_components(grad_thres, min_size)
    % 取整数类型
    grad_thres = uint8(grad_thres);
    % 连通域分析
    cc = bwconncomp(grad_thres);
    num_pixels = cellfun(@numel, cc.PixelIdxList);
    % 标记连通域是否符合条件
    labels = zeros(size(grad_thres));
    for i = 1 : cc.NumObjects
        if num_pixels(i) < min_size
            continue
        end
        labels(cc.PixelIdxList{i}) = 1;
    end
    % 输出边缘图像
    edges = uint8(labels) * 255;
end

以上代码均实现了高斯滤波、sobel算子计算梯度、计算梯度幅值和方向、非最大抑制、双阈值处理和连通域分析，可以实现 Canny 算子进行边缘检测。