Matlab 亚像素边缘提取代码详解及实现

这段代码实现了在 Matlab 中进行亚像素边缘提取的过程。具体步骤如下：

1. 读取图像并转换为灰度图像。
2. 使用 Canny 边缘检测算法提取图像的边缘。
3. 显示边缘图像。
4. 获取边缘像素的坐标。
5. 计算边缘参数，其中使用了 Zernike 矩计算函数。
6. Zernike 矩计算函数中，首先将图像转换为双精度类型并获取图像尺寸。
7. 然后计算图像的中心。
8. 接着计算归一化的 Zernike 矩，其中使用了 Zernike 多项式计算函数。
9. Zernike 多项式计算函数中，根据给定的 n、m 和 rho 计算 Zernike 多项式。
10. 最后将 Zernike 矩存储起来。

这段代码的输出是显示亚像素边缘提取后的图像。

% 读取图像
img = imread('lena.jpg');
% 将图像转换为灰度图像
grayImg = rgb2gray(img);
% Canny边缘检测
edges = edge(grayImg, 'canny');
% 显示边缘图像
figure;
imshow(edges);
title('Canny边缘图像');
% 获取边缘像素的坐标
[row, col] = find(edges);
% 计算边缘参数
edgeParams = zeros(length(row), 25); % 修改数组大小为1×25
for i = 1:length(row)
    % 获取当前像素的邻域
    neighborhood = grayImg(row(i)-1:row(i)+1, col(i)-1:col(i)+1);
    % 计算Zernike矩
    zernikeMoments = ZernikeMomentsnew(neighborhood);
    % 存储边缘参数
    edgeParams(i, :) = zernikeMoments;
end
% Zernike矩计算函数
function zernikeMoments = ZernikeMomentsnew(img)
    % 将图像转换为双精度类型
    img = im2double(img);
    % 获取图像尺寸
    [M, N] = size(img);
    % 计算图像的中心
    centerX = (N+1)/2;
    centerY = (M+1)/2;
    % 计算归一化的Zernike矩
    zernikeMoments = zeros(1, 25); % 修改数组大小为1×25
    for n = 0:4
        for m = -n:2:n
            % 初始化Zernike矩
            moment = 0;
            % 计算每个像素的Zernike矩
            for i = 1:M
                for j = 1:N
                    % 计算像素的极坐标
                    rho = sqrt((i-centerY)^2 + (j-centerX)^2) / centerY;
                    theta = atan2(j-centerX, i-centerY);
                    % 计算Zernike多项式
                    R = ZernikePolynomialnew(n, m, rho);
                    % 计算Zernike矩
                    moment = moment + img(i, j) * R * exp(-1i * m * theta);
                end
            end
            % 归一化Zernike矩
            moment = moment * (n+1) / pi;
            % 存储Zernike矩
            zernikeMoments(n^2+n+m+1) = moment;
        end
    end
end
% Zernike多项式计算函数
function R = ZernikePolynomialnew(n, m, rho)
    R = 0;
    for s = 0:(n-abs(m))/2
        R = R + (-1)^s * factorial(n-s) / (factorial(s) * factorial((n+abs(m))/2-s) * factorial((n-abs(m))/2-s)) * rho^(n-2*s);
    end
    R = R * sqrt((n+1) / pi);
end