在MATLAB中，写canny算法对图像进行边缘检测算法，不要使用MATLAB自带的canny算子

可以按照以下步骤编写Canny算法的MATLAB代码：

1. 读取图像并将其转换为灰度图像。

img = imread('image.jpg');
gray = rgb2gray(img);

2. 对灰度图像进行高斯滤波以平滑噪声。

sigma = 1.4;
kernelSize = 2 * ceil(3 * sigma) + 1; % 根据sigma计算高斯核大小
gaussianFilter = fspecial('gaussian', kernelSize, sigma); % 创建高斯滤波器
smoothed = imfilter(gray, gaussianFilter);

3. 计算图像的梯度幅值和方向。

sobelX = [-1 0 1; -2 0 2; -1 0 1]; % Sobel算子
sobelY = sobelX';
gradientX = imfilter(smoothed, sobelX);
gradientY = imfilter(smoothed, sobelY);
gradientMag = sqrt(gradientX.^2 + gradientY.^2); % 梯度幅值
gradientDir = atan2(gradientY, gradientX); % 梯度方向

4. 对梯度幅值进行非极大值抑制，以保留图像中的细线条。

gradientDir(gradientDir < 0) = gradientDir(gradientDir < 0) + pi; % 调整梯度方向到[0, pi)
gradientDir(gradientDir >= pi/2) = gradientDir(gradientDir >= pi/2) - pi/2; % 将梯度方向离散为0°、45°、90°、135°
quantizedDir = zeros(size(gradientDir));
quantizedDir(gradientDir < pi/4) = 0; % 0°
quantizedDir(gradientDir >= pi/4 & gradientDir < pi/2) = 45; % 45°
quantizedDir(gradientDir >= pi/2 & gradientDir < 3*pi/4) = 90; % 90°
quantizedDir(gradientDir >= 3*pi/4) = 135; % 135°
nms = zeros(size(gradientMag));
for i = 2:size(gradientMag, 1)-1
    for j = 2:size(gradientMag, 2)-1
        switch quantizedDir(i, j)
            case 0 % 梯度方向为0°
                if gradientMag(i, j) > gradientMag(i, j-1) && gradientMag(i, j) > gradientMag(i, j+1)
                    nms(i, j) = gradientMag(i, j);
                end
            case 45 % 梯度方向为45°
                if gradientMag(i, j) > gradientMag(i-1, j-1) && gradientMag(i, j) > gradientMag(i+1, j+1)
                    nms(i, j) = gradientMag(i, j);
                end
            case 90 % 梯度方向为90°
                if gradientMag(i, j) > gradientMag(i-1, j) && gradientMag(i, j) > gradientMag(i+1, j)
                    nms(i, j) = gradientMag(i, j);
                end
            case 135 % 梯度方向为135°
                if gradientMag(i, j) > gradientMag(i-1, j+1) && gradientMag(i, j) > gradientMag(i+1, j-1)
                    nms(i, j) = gradientMag(i, j);
                end
        end
    end
end

5. 进行双阈值处理以得到强边缘和弱边缘。

thresholdLow = 0.1;
thresholdHigh = 0.3;
maxValue = max(nms(:));
thresholdLow = thresholdLow * maxValue;
thresholdHigh = thresholdHigh * maxValue;
strongEdges = nms > thresholdHigh;
weakEdges = nms > thresholdLow & nms <= thresholdHigh;

6. 对弱边缘进行连接，以找到与强边缘相邻的弱边缘，并将其标记为强边缘。

mask = strongEdges;
while true
    [r, c] = find(weakEdges);
    if isempty(r) || isempty(c)
        break;
    end
    idx = sub2ind(size(weakEdges), r, c);
    [maxVal, maxIdx] = max(nms(idx));
    if maxVal > thresholdHigh
        mask(r(maxIdx), c(maxIdx)) = true;
    end
    weakEdges(r(maxIdx), c(maxIdx)) = false;
end

7. 返回检测到的边缘图像。

result = mask;

完整代码如下：

img = imread('image.jpg');
gray = rgb2gray(img);

sigma = 1.4;
kernelSize = 2 * ceil(3 * sigma) + 1;
gaussianFilter = fspecial('gaussian', kernelSize, sigma);
smoothed = imfilter(gray, gaussianFilter);

sobelX = [-1 0 1; -2 0 2; -1 0 1];
sobelY = sobelX';
gradientX = imfilter(smoothed, sobelX);
gradientY = imfilter(smoothed, sobelY);
gradientMag = sqrt(gradientX.^2 + gradientY.^2);
gradientDir = atan2(gradientY, gradientX);

gradientDir(gradientDir < 0) = gradientDir(gradientDir < 0) + pi;
gradientDir(gradientDir >= pi/2) = gradientDir(gradientDir >= pi/2) - pi/2;
quantizedDir = zeros(size(gradientDir));
quantizedDir(gradientDir < pi/4) = 0;
quantizedDir(gradientDir >= pi/4 & gradientDir < pi/2) = 45;
quantizedDir(gradientDir >= pi/2 & gradientDir < 3*pi/4) = 90;
quantizedDir(gradientDir >= 3*pi/4) = 135;
nms = zeros(size(gradientMag));
for i = 2:size(gradientMag, 1)-1
    for j = 2:size(gradientMag, 2)-1
        switch quantizedDir(i, j)
            case 0
                if gradientMag(i, j) > gradientMag(i, j-1) && gradientMag(i, j) > gradientMag(i, j+1)
                    nms(i, j) = gradientMag(i, j);
                end
            case 45
                if gradientMag(i, j) > gradientMag(i-1, j-1) && gradientMag(i, j) > gradientMag(i+1, j+1)
                    nms(i, j) = gradientMag(i, j);
                end
            case 90
                if gradientMag(i, j) > gradientMag(i-1, j) && gradientMag(i, j) > gradientMag(i+1, j)
                    nms(i, j) = gradientMag(i, j);
                end
            case 135
                if gradientMag(i, j) > gradientMag(i-1, j+1) && gradientMag(i, j) > gradientMag(i+1, j-1)
                    nms(i, j) = gradientMag(i, j);
                end
        end
    end
end

thresholdLow = 0.1;
thresholdHigh = 0.3;
maxValue = max(nms(:));
thresholdLow = thresholdLow * maxValue;
thresholdHigh = thresholdHigh * maxValue;
strongEdges = nms > thresholdHigh;
weakEdges = nms > thresholdLow & nms <= thresholdHigh;

mask = strongEdges;
while true
    [r, c] = find(weakEdges);
    if isempty(r) || isempty(c)
        break;
    end
    idx = sub2ind(size(weakEdges), r, c);
    [maxVal, maxIdx] = max(nms(idx));
    if maxVal > thresholdHigh
        mask(r(maxIdx), c(maxIdx)) = true;
    end
    weakEdges(r(maxIdx), c(maxIdx)) = false;
end

result = mask;