%**************************************************************************
% 图像检索——形状特征提取
%利用HU的七个不变矩作为形状特征向量
%Image : 输入图像数据
%n : 返回七维形状特征行向量
%**************************************************************************
function n = CBIR_Shape(Image)
[M,N,O] = size(Image);
%--------------------------------------------------------------------------
%彩色图像灰度化
%--------------------------------------------------------------------------
% Gray = double(0.3*Image(:,:,1)+0.59*Image(:,:,2)+0.11*Image(:,:,3));
Gray=rgb2gray(Image);
%--------------------------------------------------------------------------
%用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);
%对图像二值化
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(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-3*u21)^2)+(u30-3*u21)*(u03+u21)*((u03+u21)^2-3*(u30+u12)^2);
% %--------------------------------------------------------------------------
% %内部归一化：
% %--------------------------------------------------------------------------
% en = mean(n);
% delta = sqrt(cov(n));
% n = abs(n-en)/(3*delta);
