Hu不变矩，用matlab实现

实现Hu不变矩前需要先进行图像的预处理，包括灰度化、二值化和图片剪裁等。接下来是Hu不变矩的具体实现代码：

1. 首先导入图片，将其转换为灰度图像

I = imread('test.jpg'); % 导入图片
I = rgb2gray(I); % 转换为灰度图像

2. 对图像进行二值化处理

BW = imbinarize(I, graythresh(I)); % 对图像进行二值化
imshow(BW); % 显示二值化处理后的图像

3. 图片剪裁

x = 100; % 剪裁起点的x坐标
y = 100; % 剪裁起点的y坐标
w = 100; % 剪裁的宽度
h = 100; % 剪裁的高度
cropImage = imcrop(BW, [x y w h]); % 进行剪裁操作
imshow(cropImage); % 显示剪裁后的图像

4. 计算Hu不变矩

% 计算Hu不变矩
moment = hu_moments(cropImage); % hu_moments是自定义函数

5. 自定义函数hu_moments的实现

function [moment] = hu_moments(I)
    moment = zeros(1,7);

    % 计算几何中心
    [m, n] = size(I);
    x_sum = 0;
    y_sum = 0;
    for i = 1:m
        for j = 1:n
            x_sum = x_sum + i*I(i,j);
            y_sum = y_sum + j*I(i,j);
        end
    end
    x_center = x_sum/sum(I(:));
    y_center = y_sum/sum(I(:));
    
    % 矩阵归一化
    mu = zeros(3,3);
    for i = 0:2
        for j = 0:2
            mu(i+1,j+1) = moment_calc(I,i,j,x_center,y_center);
        end
    end
    
    % 计算二阶中心距
    m20 = mu(2,0) - mu(1,1)^2/mu(0,0);
    m02 = mu(0,2) - mu(1,1)^2/mu(0,0);
    m11 = (mu(1,1) - x_center*y_center)/mu(0,0);
    m30 = moment_calc(I,3,0,x_center,y_center) - 3*m11*m20 + 2*m11^3/mu(0,0);
    m03 = moment_calc(I,0,3,x_center,y_center) - 3*m11*m02 + 2*m11^3/mu(0,0);
    m21 = moment_calc(I,2,1,x_center,y_center) - 2*m11*m20 - m02*m11/mu(0,0) + 2*m11^3/mu(0,0)^2;
    m12 = moment_calc(I,1,2,x_center,y_center) - 2*m11*m02 - m20*m11/mu(0,0) + 2*m11^3/mu(0,0)^2;
    
    % 计算Hu不变矩
    moment(1) = m20 + m02;
    moment(2) = (m20 - m02)^2 + 4*m11^2;
    moment(3) = (m30 - 3*m12)^2 + (3*m21 - m03)^2;
    moment(4) = (m30 + m12)^2 + (m21 + m03)^2;
    moment(5) = (m30 - 3*m12)*(m30 + m12)*((m30 + m12)^2 - 3*(m21 + m03)^2)...
                + (3*m21 - m03)*(m21 + m03)*(3*(m30 + m12)^2 - (m21 + m03)^2);
    moment(6) = (m20 - m02)*((m30 + m12)^2 - (m21 + m03)^2) + 4*m11*(m30 + m12)*(m21 + m03);
    moment(7) = (3*m21 - m03)*(m30 + m12)*((m30 + m12)^2 - 3*(m21 + m03)^2)...
                - (m30 - 3*m12)*(m21 + m03)*(3*(m30 + m12)^2 - (m21 + m03)^2);
end

function [Iij] = moment_calc(I,i,j,x,y)
    [m, n] = size(I);
    Iij = 0;
    for p = 1:m
        for q = 1:n
            Iij = Iij + (p-x)^i*(q-y)^j*I(p,q);
        end
    end
end