%**************************************************************************
% 图像检索——形状特征提取
% 利用HU的七个不变矩作为形状特征向量
% Image : 输入图像数据
% n: 返回七维形状特征行向量
%**************************************************************************
function n = HuSquare(Image)
% Image = imread('E:\\1\\1.jpg');
% M = 256;
% N = 256;
% 如果图像不是彩色图像的处理,变为3通道
if length(size(Image))==2
Image1=zeros(size(Image,1),size(Image,2),3);
Image1(:,:,1)=Image;
Image1(:,:,2)=Image;
Image1(:,:,3)=Image;
Image=Image1;
end
[M,N,O] = size(Image);
%--------------------------------------------------------------------------
%彩色图像灰度化
%--------------------------------------------------------------------------
Gray = double(0.3*Image(:,:,1)+0.59*Image(:,:,2)+0.11*Image(:,:,3));
%--------------------------------------------------------------------------
%用Canny边缘检测提取边缘,保留边缘灰度图像
%--------------------------------------------------------------------------
% BW = uint8(edge(Gray,'canny'));
Egray = uint8(edge(Gray,'canny'));
for i = 1:M
for j = 1:N
if Egray(i,j)==0
Gray(i,j)=0;
end
end
end
%--------------------------------------------------------------------------
%Otsu提出的类判别分析法自动为每一幅廓图像选定阈值,然后用该阈值对图像二值化
%--------------------------------------------------------------------------
% 计算灰度级归一化直方图
for i = 0:255
h(i+1) = size(find(Gray==i),1);
end
p = h/sum(h);
% 计算灰度均值
ut = 0;
for i = 0:255
ut = i*p(i+1)+ut;
end
% 计算直方图的零阶累积矩和一阶累积矩:
for k = 0:254
w(k+1) = sum(p(1:k+1));
u(k+1) = sum((0:k).*p(1:k+1));
end
% 计算类分离指标
deltaB = zeros(1,255);
for k = 0:254
if w(k+1)~=0&w(k+1)~=1
deltaB(k+1) = (ut*w(k+1)-u(k+1))^2/(w(k+1)*(1-w(k+1)));
end
end
[value,thresh] = max(deltaB);
% deltaB = zeros(1,255);
% delta1 = zeros(1,255);
% delta2 = zeros(1,255);
% deltaW = zeros(1,255);
% for k = 0:254
% if w(k+1)~=0&w(k+1)~=1
% deltaB(k+1) = (ut*w(k+1)-u(k+1))^2/(w(k+1)*(1-w(k+1)));
% delta1(k+1) = 0;
% delta2(k+1) = 0;
% for i = 0:k
% delta1(k+1) = (i-u(k+1)/w(k+1))^2*p(i+1)+delta1(k+1);
% end
% for i = k+1:255
% delta2(k+1) = (i-(ut-u(k+1))/(1-w(k+1)))^2*p(k+1)+delta2(k+1);
% end
% deltaW(k+1) = delta1(k+1)+delta2(k+1);
% end
% end
% for i = 1:255
% if deltaB==0
% yita=0;
% else
% yita(i) = 1/(1+deltaW(i)./deltaB(i));
% end
% end
% % D的最大值作为最佳阈值
% [value,thresh] = max(yita);
%对图像二值化
for i = 1:M
for j = 1:N
if Gray(i,j)>=thresh
BW(i,j) = 1;
else
BW(i,j) = 0;
end
end
end
%--------------------------------------------------------------------------
%计算图像质心:(I,J)
%--------------------------------------------------------------------------
m00 = sum(sum(BW)); %零阶矩
m01 = 0; %一阶矩
m10 = 0; %一阶矩
for i = 1:M
for j = 1:N
m01 = BW(i,j)*j+m01;
m10 = BW(i,j)*i+m10;
end
end
I = (m10)/(m00);
J = m01/m00;
%--------------------------------------------------------------------------
%中心矩:
%--------------------------------------------------------------------------
u11 = 0;
u20 = 0; u02 = 0;
u30 = 0; u03 = 0;
u12 = 0; u21 = 0;
for i = 1:M
for j = 1:N
u20 = BW(i,j)*(i-I)^2+u20;
u02 = BW(i,j)*(j-J)^2+u02;
u11 = BW(i,j)*(i-I)*(j-J)+u11;
u30 = BW(i,j)*(i-I)^3+u30;
u03 = BW(i,j)*(j-J)^3+u03;
u12 = BW(i,j)*(i-I)*(j-J)^2+u12;
u21 = BW(i,j)*(i-I)^2*(j-J)+u21;
end
end
u20 = u20/m00^2;
u02 = u02/m00^2;
u11 = u11/m00^2;
u30 = u30/m00^(5/2);
u03 = u03/m00^(5/2);
u12 = u12/m00^(5/2);
u21 = u21/m00^(5/2);
%--------------------------------------------------------------------------
%7个Hu不变矩:
%--------------------------------------------------------------------------
n(1) = u20+u02;
n(2) = (u20-u02)^2+4*u11^2;
n(3) = (u30-3*u12)^2+(u03-3*u21)^2;
n(4) = (u30+u12)^2+(u03+u21)^2;
% n(5) = (u30-3*u12)*(u30+u12)*((u30+u12)^2-3*(u03-3*u21)^2)+(u03-3*u21)*(u03+u21)*((u03+u21)^2-3*(u30+u12)^2);
% n(5) = (u30-3*u12)*(u30+u12)*((u30+u12)^2-3*(u03+u21)^2)+(u03-3*u21)*(u03+u21)*((u03+u21)^2-3*(u30+u12)^2);
n(5) = (u03-3*u12)*(u30+u12)*((u30+3*u12)^2-3*(u03+u21)^2)+(3*u21-u03)*(u03+u21)*((u03+u21)^2-3*(u30+u12)^2);
n(6) = (u20-u02)*((u30+u12)^2-(u03+u21)^2)+4*u11*(u30+u12)*(u03+u21);
n(7) = (3*u21-u03)*(u30+u12)*((u30+u12)^2-3*(u03+u21)^2)+(3*u21-u30)*(u03+u21)*(3*(u30+u12)^2-(u03+u21)^2);
% n(7) = (3*u21-u03)*(u30+u12)*((u30+u12)^2-3*(u03+u21)^2)+(u30-3*u12)*(u03+u21)*((u03+u21)^2-3*(u30+u12)^2);
% %--------------------------------------------------------------------------
% %内部归一化:
% %--------------------------------------------------------------------------
% en = mean(n);
% delta = sqrt(cov(n));
% n = abs(n-en)/(3*delta);
end
Image = imread('../images/2.jpg');
n = HuSquare(Image)