文章目录
图像处理特征分析
一、颜色特征描述
三原色RGB
色调饱和度亮度HSV
%拆分一个HSV图像的图像阵列
RGB=reshape(ones(64,1)*reshape(jet(64),1,192),[64,64,3]); %调整颜色条尺寸为正方形
HSV=rgb2hsv(RGB); %将RGB图像转换为HSV图像
H=HSV(:,:,1); %提取H矩阵
S=HSV(:,:,2); %提取S矩阵
V=HSV(:,:,3); %提取V矩阵
figure(1)
subplot(121), imshow(H) ,title('H'); %显示H图像
subplot(122), imshow(S) ,title('S'); %显示S图像
figure(2)
subplot(121), imshow(V) ,title('V'); %显示V图像
subplot(122), imshow(RGB),title('RGB'); %显示RGB图像
颜色矩
一阶矩 定义每个颜色分量平均强度 mean2()
二阶矩 反应待测区域的颜色方差 即不均匀性 std()
三阶矩 定义颜色分量的偏斜度 即颜色的不对称性
颜色矩不能区分颜色区域的空间分布位置
彩色图像的颜色矩一共9个分量
%【例11-1】颜色矩求法
close all; %关闭当前所有图形窗口，清空工作空间变量，清除工作空间所有变量
clear all;
clc;
I=imread('hua.jpg'); %I为花的彩色图像，以下是求花的图像的RGB分量均值
R=I(:,:,1); %红色分量
G=I(:,:,2); %绿色分量
B=I(:,:,3); %蓝色分量
R=double(R); G=double(G); B=double(B); %利用double()函数将变量类型转为double型
Ravg1=mean2(R) %红色分量均值
Gavg1=mean2(G) %绿色分量均值
Bavg1=mean2(B) %蓝色分量均值
Rstd1=std(std(R)) %红色分量的方差
Gstd1= std(std(G)) %绿色分量的方差
Bstd1=std(std(B)) %蓝色分量的方差
J=imread('yezi.jpg'); %J为叶子的彩色图像以下是求叶子的图像的RGB分量均值
R=J(:,:,1); %红色分量
G=J(:,:,2); %绿色分量
B=J(:,:,3); %蓝色分量
R=double(R); G=double(G); B=double(B); %利用double()函数将变量类型转为double型
fprintf('yezi');fprintf('红色分量均值');
Ravg2=mean2(R) %红色分量均值
fprintf('绿色分量均值'); Gavg2=mean2(G) %绿色分量均值
fprintf('蓝色分量均值');Bavg2=mean2(B) %蓝色分量均值
fprintf('红色分量方差');Rstd2=std(std(R)) %红色分量的方差
fprintf('绿色分量方差');Gstd2= std(std(G)) %绿色分量的方差
fprintf('蓝色分量均值');Bstd2=std(std(B)) %蓝色分量的方差
set(0,'defaultFigurePosition',[100,100,1000,500]); %修改图形图像位置的默认设置
set(0,'defaultFigureColor',[1 1 1])
K=imread('flower1.jpg');figure;subplot(131),imshow(K); %显示原图像
subplot(132),imshow(I); %显示花的图像
subplot(133),imshow(J); %显示叶子的图像
%基于颜色特征识别花朵和叶子
close all; %关闭当前所有图形窗口，清空工作空间变量，清除工作空间所有变量
clear all;
clc;
I=imread('hua.jpg'); %I为花的彩色图像，以下是求花的图像的RGB分量均值
R=I(:,:,1); %红色分量
G=I(:,:,2); %绿色分量
B=I(:,:,3); %蓝色分量
R=double(R); G=double(G); B=double(B); %利用double()函数将变量类型转为double型
Ravg1=mean2(R); %红色分量均值
Gavg1=mean2(G); %绿色分量均值
Bavg1=mean2(B); %蓝色分量均值
Rstd1=std(std(R)); %红色分量的方差
Gstd1= std(std(G)); %绿色分量的方差
Bstd1=std(std(B)); %蓝色分量的方差
J=imread('yezi.jpg'); %J为叶子的彩色图像以下是求叶子的图像的RGB分量均值
R=J(:,:,1); %红色分量
G=J(:,:,2); %绿色分量
B=J(:,:,3); %蓝色分量
R=double(R); G=double(G); B=double(B); %利用double()函数将变量类型转为double型
Ravg2=mean2(R); %红色分量均值
Gavg2=mean2(G); %绿色分量均值
Bavg2=mean2(B); %蓝色分量均值
Rstd2=std(std(R)); %红色分量的方差
Gstd2= std(std(G)); %绿色分量的方差
Bstd2=std(std(B)); %蓝色分量的方差
set(0,'defaultFigurePosition',[100,100,1000,500]); %修改图形图像位置的默认设置
set(0,'defaultFigureColor',[1 1 1])
K=imread('flower1.jpg');figure;subplot(131),imshow(K); %显示原图像
subplot(132),imshow(I); %显示花的图像
subplot(133),imshow(J); %显示叶子的图像
颜色直方图
描述不同色彩在整幅图像中的占比，不关心每种色彩所处的空间位置。适用于描述那些难以自动分割的图像和不需要考虑物体空间位置的图像
%绘制彩色图像R,G,B分量的灰度直方图
I=imread('huangguahua.jpg'); %读入要处理的图像，并赋值给I
R=I(:,:,1); %图像的R分量
G=I(:,:,2); %图像的G分量
B=I(:,:,3); %图像的B分量
set(0,'defaultFigurePosition',[100,100,1000,500]); %修改图形图像位置的默认设置
set(0,'defaultFigureColor',[1 1 1])
figure;subplot(121);imshow(I),title('原始图像'); %显示彩色图像
subplot(122);imshow(R),title('红色分量灰度图'); %R分量灰度图
figure;subplot(121);imshow(G),title('绿色分量灰度图'); %G分量灰度图
subplot(122);imshow(B),title('蓝色分量灰度图'); %B分量灰度图
figure;subplot(131);
imhist(I(:,:,1)) ,title('红色分量直方图'); %显示红色分辨率下的直方图
subplot(132);imhist(I(:,:,2)) ,title('绿色分量直方图'); %显示绿色分辨率下的直方图
subplot(133);imhist(I(:,:,3)) ,title('蓝色分量直方图'); %显示蓝色分辨率下的直方图
二、纹理特征描述
对应物体表面性质，图像纹理特征常常具有周期性，反应物体质地，如粗糙度、光滑度、颗粒度、随机性和规范性
1.灰度差分
图像中一点和与它只有微小距离的点的灰度差值
相关纹理特征 均值 对比度 熵
%例【11-5】灰度差分统计特征计算，分析纹理图像
J=imread('stone.jpg'); %读入纹理图像，分别输入wall.jpg和stone.jpg两幅图进行对比
A=double(J);
[m,n]=size(A); %求A矩阵的大小，赋值给m n
B=A;
C=zeros(m,n); %新建全零矩阵C，以下求解归一化的灰度直方图
for i=1:m-1
for j=1:n-1
B(i,j)=A(i+1,j+1);
C(i,j)=abs(round(A(i,j)-B(i,j)));
end
end
h=imhist(mat2gray(C))/(m*n);
mean=0;con=0;ent=0; % 均值mean、对比度con和熵ent初始值赋零
for i=1:256 %循环求解均值mean、对比度con和熵ent
mean=mean+(i*h(i))/256;
con=con+i*i*h(i);
if(h(i)>0)
ent=ent-h(i)*log2(h(i));
end
end
imshow(J)
从灰度均值特征来看，墙面纹理图的灰度均值和熵都大于大理石纹理图说明墙面纹理较为粗糙。对比度差异大可以用来区分墙面和大理石面
2.自相关函数
纹理粗糙性的大小与局部结构的空间重复周期有关，周期大的纹理粗，周期小的纹理细，通常采用自相关函数作为纹理测度
%【例11-6】
%步骤1：定义自相关函数zxcor()，建立zxcor.m文件
function [epsilon,eta,C]=zxcor(f,D,m,n)
%自相关函数zxcor(),f为读入的图像数据，D为偏移距离，【m，n】是图像的尺寸数据，返回图像相关函数C的
%值，epsilon和eta是自相关函数C的偏移变量。
for epsilon=1:D %循环求解图像f(x,y)与偏离值为D的像素之间的相关值
for eta=1:D
temp=0;
fp=0;
for x=1:m
for y=1:n
if(x+ epsilon -1)>m|(y+ eta -1)>n
f1=0;
else
f1=f(x,y)*f(x+ epsilon -1,y+ eta -1);
end
temp=f1+temp;
fp=f(x,y)*f(x,y)+fp;
end
end
f2(epsilon, eta)=temp;
f3(epsilon, eta)=fp;
C(epsilon, eta)= f2(epsilon, eta)/ f3(epsilon, eta); %相关值C
end
end
epsilon =0:(D-1); % 方向的取值范围
eta =0:(D-1);
%步骤2：调用zxcor()函数，分析不同图像的纹理特征。
f11=imread('wall.jpg'); %读入砖墙面图像，图像数据赋值给f
f1=rgb2gray(f11); %彩色图像转换成灰度图像
f1=double(f1); %图像数据变为double类型
[m,n]=size(f1); %图像大小赋值为[m,n]
D=20; %偏移量为20
[epsilon1,eta1,C1]=zxcor(f1,D,m,n); %调用自相关函数
f22=imread('stone.jpg'); %读入大理石图像，图像数据赋值给f
f2=rgb2gray(f22);
f2=double(f2);
[m,n]=size(f2);
[epsilon2,eta2,C2]=zxcor(f2,20,m,n); %调用自相关函数
set(0,'defaultFigurePosition',[100,100,1000,500]); %修改图形图像位置的默认设置
set(0,'defaultFigureColor',[1 1 1]);
figure;subplot(121);imshow(f11),title('墙面纹理图');
subplot(122);imshow(f22),title('大理石面纹理图');
figure;subplot(121);mesh(epsilon1,eta1,C1); %显示自相关函数与x，y的三维图像
xlabel(' epsilon ');ylabel(' eta '),title('砖墙面图像纹理的自相关函数的三维图'); %标示坐标轴变量
subplot(122);mesh(epsilon2,eta2,C2);
xlabel(' epsilon ');ylabel(' eta '),title('大理石图像纹理的自相关函数的三维图');
3.灰度共生矩
由于纹理是由灰度分布在空间位置上反复出现而形成的，因而在图像空间中相隔某距离的两像素之间会存在一定的灰度关系，即图中灰度的空间相关特性，灰度共生矩就是一种通过研究灰度的空间相关特性来描述纹理的常用方法
灰度共生矩是对图像上保持某距离的两像素分别具有某灰度的状况进行统计得到的
一幅图像的灰度共生矩阵能反应出图像灰度关于方向，相邻间隔和变化幅度的综合信息，他是分析图像的局部模式和他们排列规则的基础
一般来说，如果图像是由具有相似灰度值的像素块构成，则灰度共生矩对角元素会有比较大的值，如果图像像素灰度值在局部有变化，那么偏离对角线的元素会有比较大的值
能量：元素值的平方和，反应灰度均匀程度和纹理粗细程度 若共生矩阵的所有制=值相等则能量值小，相反，若其中一些值大而其它值小，能量值大。当共生矩阵中元素集中分布时，能量值最大 能量值大表明一种较均一和规则变化的纹理模式
对比度：反映了图像的清晰度和纹理沟深浅程度 纹理钩纹越深，对比度越大 灰度共生矩种远离对角线元素值越大，对比度越大
相关：它度量空间灰度共生矩阵元素在行或列方向上的相似程度。当矩阵元素值均匀相等时，相关值就大
熵：不确定性 当共生矩阵种所有元素有最大的随机性，空间共生矩阵几乎所有值都相等时，共生矩阵中元素散列分布时，熵较大它表示了图像中纹理的非均匀程度或者复杂度
均匀度：翻译图像纹理的粗糙度 粗纹理均匀度大 细纹理均匀的小
%实现图像灰度共生矩阵
I = imread('circuit.tif'); %读入图像circuit.tif
glcm = graycomatrix(I,'Offset',[0 2]); %图像I的灰度共生矩阵，2表示当前像素与邻居的距离为2，offset为[0 2]表示角度为0为水平方向
set(0,'defaultFigurePosition',[100,100,1000,500]); %修改图形图像位置的默认设置
set(0,'defaultFigureColor',[1 1 1])
imshow(I);
glcm
%【例11-8】遥感图像基于灰度共生矩的纹理特征统计
I=imread('hill.jpg');
HSV=rgb2hsv(I);
Hgray=rgb2gray(HSV);
%计算64位灰度共生矩阵
glcms1=graycomatrix(Hgray,'numlevels',64,'offset',[0 1;-1 1;-1 0;-1 -1]);
%纹理特征统计值(包括对比度、相关性、熵、平稳度、二阶矩也叫能量)
stats=graycoprops(glcms1,{'contrast','correlation','energy','homogeneity'});
ga1=glcms1(:,:,1);%0度
ga2=glcms1(:,:,2);%45度
ga3=glcms1(:,:,3);%90度
ga4=glcms1(:,:,4);%135度
energya1=0;energya2=0;energya3=0;energya4=0;
for i=1:64
for j=1:64
energya1=energya1+sum(ga1(i,j)^2);
energya2=energya2+sum(ga2(i,j)^2);
energya3=energya3+sum(ga3(i,j)^2);
energya4=energya4+sum(ga4(i,j)^2);
j=j+1;
end
i=i+1;
end
s1=0;s2=0;s3=0;s4=0;s5=0;
for m=1:4
s1=stats.Contrast(1,m)+s1;%对比度
m=m+1;
end
for m=1:4
s2=stats.Correlation(1,m)+s2;%相关性
m=m+1;
end
for m=1:4
s3=stats.Energy(1,m)+s3;%熵
m=m+1;
end
for m=1:4
s4=stats.Homogeneity(1,m)+s4;%平稳度
m=m+1;
end
s5=0.000001*(energya1+energya2+energya3+energya4);%二阶矩 （能量)
I=imread('hill.jpg');
J=imread('sea.jpg');
K=imread('house.jpg');
set(0,'defaultFigurePosition',[100,100,1000,500]); %修改图形图像位置的默认设置
set(0,'defaultFigureColor',[1 1 1])
figure;subplot(131);imshow(I);
subplot(132);imshow(J);
subplot(133);imshow(K);
a=rgb2gray(K);
glcm=graycomatrix(a);%实现图像灰度共生矩
4 频谱分析法
将空间的纹理图像变换到频率域中
基于傅里叶变换分析纹理图像
%【例11-9】基于傅里叶变换分析纹理图像
I = imread('wall.jpg'); %读入图像
I=rgb2gray(I); %图像变为灰度图像
wall=fft2(I); %对图像做快速傅里叶变换
s=fftshift(wall); %将变换后的图象频谱中心从矩阵的原点移到矩阵的中心
s=abs(s);
[nc,nr]=size(s);
x0=floor(nc/2+1);
y0=floor(nr/2+1);
rmax=floor(min(nc,nr)/2-1);
srad=zeros(1,rmax);
srad(1)=s(x0,y0);
thetha=91:270; %thetha取值91-270
for r=2:rmax %循环求解纹理频谱能量
[x,y]=pol2cart(thetha,r);
x=round(x)'+x0;
y=round(y)'+y0;
for j=1:length(x)
srad(r)=sum(s(sub2ind(size(s),x,y)));
end
end
[x,y]=pol2cart(thetha,rmax);
x=round(x)'+x0;
y=round(y)'+y0;
sang=zeros(1,length(x));%循环求解纹理频谱能量 sr
for th=1:length(x)
vx=abs(x(th)-x0);
vy=abs(y(th)-y0);
if((vx==0)&(vy==0))
xr=x0;
yr=y0;
else
m=(y(th)-y0)/(x(th)-x0)
xr=(x0:x(th)).'
yr=round(y0+m*(xr-x0));
end
for j=1:length(xr)
sang(th)=sum(s(sub2ind(size(s),xr,yr)));
end
end
set(0,'defaultFigurePosition',[100,100,1000,500]); %修改图形图像位置的默认设置
set(0,'defaultFigureColor',[1 1 1])
figure;subplot(121);
imshow('wall.jpg');subplot(122); %显示原图
imshow(log(abs(wall)),[]); %显示频谱图
figure;subplot(121);plot(srad); %显示
subplot(122);plot(sang); %显示
Gabor变换 Gabor函数可以在频域不同尺度 不同方向上提取相关特征
%【例11-10】
close all; %关闭当前所有图形窗口，清空工作空间变量，清除工作空间所有变量
clear all;
clc;
I = imread('wall.jpg'); %读取图像，并赋值给I
I=rgb2gray(I); %彩色图像变为灰度图像
[G,gabout] = gaborfilter(I,2,4,16,pi/10); %调用garborfilter()函数对图像做小波变换
J=fft2(gabout); %对滤波后的图像做FFT变换，变换到频域
A=double(J);
[m,n]=size(A);
B=A;
C=zeros(m,n);
for i=1:m-1
for j=1:n-1
B(i,j)=A(i+1,j+1);
C(i,j)=abs(round(A(i,j)-B(i,j)));
end
end
h=imhist(mat2gray(C))/(m*n); %对矩阵C归一化处理后求其灰度直方图，得到归一化的直方图
mean=0;con=0;ent=0;
for i=1:256 %求图像的均值、对比度和熵
mean=mean+(i*h(i))/256;
con=con+i*i*h(i);
if(h(i)>0)
ent=ent-h(i)*log2(h(i));
end
end
set(0,'defaultFigurePosition',[100,100,1000,500]);%修改图形图像位置的默认设置
set(0,'defaultFigureColor',[1 1 1])
figure;subplot(121);imshow(I);
subplot(122);imshow(uint8(gabout));
mean,con,ent
小波变换
三、形状特征描述
形状特征两类表示方法：
轮廓特征：针对物体外边界
区域特征：关系到整个形状区域
边界表示方式
链码 归一化 一阶差分
%【例11-11】
I=[1 1 1 1;1 1 0 1;0 1 0 1;0 1 1 1]; %图像数据赋值给I，I为4 4大小的矩阵
%跟踪目标的边界，返回值为一个p 1的数组单元，p为目标的个数，其中每一个单元又是一个Q 2的矩阵，即
%Q个点的x,y坐标。
g=boundaries(I,4); %追踪4连接的目标边界
c=fchcode(g{:},4); %求4方向freeman链码
c.x0y0 %显示代码开始处的坐标（1 2）
c.fcc %Freeman链码表示 大小为n 2
c.diff %代码c.fcc的一阶差分（1 n）
c.mm %最小幅度的整数
c.diffmm %代码c.mm的一阶差分（1 n）
形状数描述
多边形近似
基于最小周长多边形法
基于聚合的最小均方误差线段逼近法
基于分裂的最小均方误差线段逼近法
%【例11-14】基于最小周长多边形法描述图像边界
I=imread('leaf1.bmp'); %读入图像数据赋值给I
I=rgb2gray(I); %将彩色图像变为灰度图像
bwI=im2bw(I,graythresh(I)); %对图像进行二值化处理得到二值化图像赋值给bwI
bwIsl=~bwI; %对二值图像取反
h=fspecial('average'); %选择中值滤波
bwIfilt=imfilter(bwIsl,h); %对图像进行中值滤波
bwIfiltfh=imfill(bwIfilt,'holes'); %填充二值图像的空洞区域
bdI=boundaries(bwIfiltfh,4,'cw'); %追踪4连接目标边界
d=cellfun('length',bdI); %求bdI中每一个目标边界的长度，返回值d是一个向量
[dmax,k]=max(d); %返回向量d中最大的值，存在max_d中，k为其索引
B4=bdI{k(1)}; %若最大边界不止一条，则取出其中的一条即可。B4是一个坐标数组
[m,n]=size(bwIfiltfh); %求二值图像的大小
xmin=min(B4(:,1));
ymin=min(B4(:,2));
%生成一幅二值图像,大小为m n，xmin,ymin是B4中最小的x和y轴坐标
bim=bound2im(B4,m,n,xmin,ymin);
[x,y]=minperpoly(bwIfiltfh,2); %使用大小为2的方形单元
b2=connectpoly(x,y); %按照坐标(X,Y)顺时针或者逆时针连接成多边形
B2=bound2im(b2,m,n,xmin,ymin);
set(0,'defaultFigurePosition',[100,100,1000,500]);%修改图形图像位置的默认设置
set(0,'defaultFigureColor',[1 1 1])
figure,subplot(121);imshow(bim); %显示原图像边界
subplot(122),imshow(B2); %显示按大小为2的正方形单元近似的边界
标记图
边界特征描述
简单的边界特征描述:边界长度 边界直径 长轴 短轴和离心率
形状数 c=fchode(g{:},cn) 形状数定义为最小数量级的差分码
傅里叶形状描述子
%【例11-13】利用edge()函数提取图像轮廓，绘制出对象边界和提取边界坐标信息
I= imread('leaf1.bmp'); %读入图像
c= im2bw(I, graythresh(I)); %I转换为二值图像
set(0,'defaultFigurePosition',[100,100,1000,500]); %修改图形图像位置的默认设置
set(0,'defaultFigureColor',[1 1 1])
figure;subplot(131);imshow(I),title('原图像'); %显示原图
c=flipud(c); %实现矩阵c上下翻转
b=edge(c,'canny'); %基于canny算子进行轮廓提取
[u,v]=find(b); %返回边界矩阵b中非零元素的位置
xp=v; %行值v赋给xp
yp=u; %列值u赋给yp
x0=mean([min(xp),max(xp)]); %x0为行值的均值
y0=mean([min(yp),max(yp)]); %y0为列值的均值
xp1=xp-x0;
yp1=yp-y0;
[cita,r]=cart2pol(xp1,yp1); %直角坐标转换成极坐标
q=sortrows([cita,r]); %从r列开始比较数值并按升序排序
cita=q(:,1); %赋角度值
r=q(:,2); %赋半径模值
subplot(132);polar(cita,r),title('极坐标下的轮廓图'); %画出极坐标下的轮廓图
[x,y]=pol2cart(cita,r);
x=x+x0;
y=y+y0;
subplot(133);plot(x,y),title('直角坐标下的轮廓图'); %画出直角坐标下的轮廓图
Matlab 中提供了bwboundaries()和edge()函数提取目标轮廓
区域特征描述
简单区域描述
区域面积 区域重心
%【例11-14】利用regionprops()求区域的面积和重心
I= imread('leaf1.bmp'); %读入图像
I= im2bw(I); %转换为二值图像
C=bwlabel(I,4); %对二值图像进行4连通的标记
Ar=regionprops(C,'Area'); %求C的面积
Ce=regionprops(C,'Centroid'); %求C的重心
四叉树描述
%【例11-15】利用函数qtdecomp()对图像进行四叉树分解 此外qtsetblk也可用于获取四叉树分解块值
close all; %关闭当前所有图形窗口，清空工作空间变量，清除工作空间所有变量
clear all;
clc;
I = imread('liftingbody.png'); %读入图像
S = qtdecomp(I,.27); %四叉树分解，阈值为0.27
blocks = repmat(uint8(0),size(S)); %矩阵扩充为S的大小
for dim = [512 256 128 64 32 16 8 4 2 1];
numblocks = length(find(S==dim));
if (numblocks > 0)
values = repmat(uint8(1),[dim dim numblocks]);%左上角元素为1
values(2:dim,2:dim,:) = 0;%其他地方元素为0
blocks = qtsetblk(blocks,S,dim,values);
end
end
blocks(end,1:end) = 1; blocks(1:end,end) = 1;
set(0,'defaultFigurePosition',[100,100,1000,500]); %修改图形图像位置的默认设置
set(0,'defaultFigureColor',[1 1 1])
figure;subplot(121);imshow(I);
subplot(122), imshow(blocks,[]) %显示四叉树分解的图像
拓扑描述
区域的拓扑描述用于描述物体平面区结构形状的整体特性
只要图形不撕裂或折叠 拓扑描述的性质将不受图形变形影响
孔：H 欧拉数： EUL=C-H 图形中孔洞的数量 H 和连通分量的数量 C
A 欧拉数 0 B 欧拉数-1
欧拉数常常被用来识别数字
%【例11-16】利用函数bweuler()计算二值图形欧拉数
I=imread('number.jpg'); %读入图像
K=im2bw(I); %I转换为二值图像
J=~K; %图像取反
EUL=bweuler(J) %求图像的欧拉数
set(0,'defaultFigurePosition',[100,100,1000,500]);
set(0,'defaultFigureColor',[1 1 1])
figure;subplot(131);imshow(I),title('原图'); %绘出原图
subplot(132);imshow(K),title('二值图'); %二值图
subplot(133);imshow(J),title('取反后的图'); %取反后的图
几何形状描述
曲线长度
曲线偏心率
矩描述
若用矩描述区域目标 具有平移 旋转和缩放不变性
不变矩常用于识别二维物体 但并不能区别所有的形状 而且对噪声十分敏感
%【例11-17】求同一幅图像旋转镜像和尺寸缩小后的图像7阶矩 证明矩的不变性
close all; %关闭当前所有图形窗口，清空工作空间变量，清除工作空间所有变量
clear all;
clc
I=imread('cameraman.tif'); %读入要处理的图像，并赋值给I
%set(0,'defaultFigurePosition',[100,100,1000,500]); %修改图形图像位置的默认设置
%set(0,'defaultFigureColor',[1 1 1])
Image=I; %图像I数据赋给Image
figure;subplot(121);imshow(Image);
Image1=imrotate(I,10,'bilinear'); %图像顺时针旋转10度--旋转变换
subplot(122);imshow(Image1);
Image2=fliplr(I); %对图像做镜像变换---镜像变换
figure;subplot(121);imshow(Image2);
Image3=imresize(I,0.3,'bilinear'); %图像缩小1/3---尺寸变换
subplot(122);imshow(Image3);
%调用自定义函数Moment_Seven()求解图像七阶矩%
display('原图像');
Moment_Seven(Image);
display('旋转变化后的图像');
Moment_Seven(Image1);
display('镜像变化后的图像');
Moment_Seven(Image2);
display('尺度变化后的图像');
Moment_Seven(Image3);
%求7阶矩函数Moment_Seven()的函数清单：
function Moment_Seven(J) %J为要求解的图像
A=double(J); %将图像数据转换为double类型
[m,n]=size(A); %求矩阵A的大小
[x,y]=meshgrid(1:n,1:m); %生成网格采样点的数据，x,y的行数等于m，列数等于n
x=x(:); %矩阵赋值
y=y(:);
A=A(:);
m00=sum(A); %求矩阵A中每列的和，得到m00是行向量
if m00==0 %如果m00=0，则赋值m00=eps，即m00=0
m00=eps;
end
m10=sum(x.*A); %以下为7阶矩求解过程，参见7阶矩的公式
m01=sum(y.*A);
xmean=m10/m00;
ymean=m01/m00;
cm00=m00;
cm02=(sum((y-ymean).^2.*A))/(m00^2);
cm03=(sum((y-ymean).^3.*A))/(m00^2.5);
cm11=(sum((x-xmean).*(y-ymean).*A))/(m00^2);
cm12=(sum((x-xmean).*(y-ymean).^2.*A))/(m00^2.5);
cm20=(sum((x-xmean).^2.*A))/(m00^2);
cm21=(sum((x-xmean).^2.*(y-ymean).*A))/(m00^2.5);
cm30=(sum((x-xmean).^3.*A))/(m00^2.5);
Mon(1)=cm20+cm02; %1阶矩Mon(1)
Mon(2)=(cm20-cm02)^2+4*cm11^2; %2阶矩Mon(2)
Mon(3)=(cm30-3*cm12)^2+(3*cm21-cm03)^2; %3阶矩Mon(3)
Mon(4)=(cm30+cm12)^2+(cm21+cm03)^2; %4阶矩Mon(4)
Mon(5)=(cm30-3*cm12)*(cm30+cm12)*((cm30+cm12)^2-3*(cm21+cm03)^2)+(3*(cm30+cm12)^2-(cm21+cm03)^2); %5阶矩Mon(5)
Mon(6)=(cm20-cm02)*((cm30+cm12)^2-(cm21+cm03)^2)+4*cm11*(cm30+cm12)*(cm21+cm03); %6阶矩Mon(6)
Mon(7)=(3*cm21-cm03)*(cm30+cm12)*((cm30+cm12)^2-3*(cm21+cm03)^2)+(3*cm12-cm30)*(cm21+cm03)*(3*(cm30+cm12)^2-(cm21+cm03)^2); %7阶矩Mon(7)
Moment=abs(log(Mon)) %采用log函数缩小不变矩的动态范围值
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
