基于matlab的工业摄影作业文库.docx
- 文档编号:9205276
- 上传时间:2023-02-03
- 格式:DOCX
- 页数:24
- 大小:1.70MB
基于matlab的工业摄影作业文库.docx
《基于matlab的工业摄影作业文库.docx》由会员分享,可在线阅读,更多相关《基于matlab的工业摄影作业文库.docx(24页珍藏版)》请在冰豆网上搜索。
基于matlab的工业摄影作业文库
工业摄影测量技术
基于Matlab图像处理软件的设计
姓名:
学号:
杨潇
M201370503
学科专业:
指导教师:
精密仪器及机械
唐立新教授
提交日期:
2013.11.07
一.摘要
基于Matlab的数字图像处理
Matlab对于技术计算来说是一种高性能的语言。
它以易于应用的环境集成了计算、可视化和编程,在该环境下,问题及其解以我们熟悉的数学表示法来表示。
本软件主要通过matlab实现以下几个功能:
1.Log算子边缘提取
2.Canny算子边缘提取
3.几种滤波器的实现
4.数出图像中米粒个数
5.将全方位图像展开
6.分割图像中环形针并标出其质心及方向。
本次简易的图像处理软件是利用Matlab的GUI界面来实现。
关键字:
matlab,Log,Canny,滤波器,数米粒,图像展开,方位辨识
2.基本图像处理
本软件界面如图2.1所示:
图2.1软件界面
界面设计是通过任务要求来设计,包含Log算子边缘提取、Canny算子边缘提取、滤波、数数目、图像展开及方位分析。
左右两个显示屏分别显示图像处理前后的图像。
1.Log算子边缘提取的实现
Log算子可有效的滤除噪声信号并提取出清晰的边缘轮廓。
具体分为以下4步:
(1)对图像进行高斯平滑滤波,去除噪声信号
关键代码:
n1=floor((n+1)/2);%?
?
?
?
?
?
?
?
fori=1:
n
forj=1:
n
b(i,j)=exp(-((i-n1)^2+(j-n1)^2)/(2*sigma));%?
?
?
?
?
?
end
end
A=conv2(I,b,'same');%?
?
?
?
其中,n为高斯模板大小,sigma为方差
(2)利用二维拉普拉斯函数进行二阶求导争强边缘信号
关键代码:
fori=2:
w-1
forj=2:
h-1
f=A(i+1,j+1)+A(i+1,j-1)+A(i-1,j-1)+A(i-1,j+1)+A(i+1,j)+A(i-1,j)+A(i,j+1)+A(i,j-1)-8*A(i,j);%?
?
?
?
?
?
?
?
?
其中,w,h分别表示图像的长宽。
通过图像A与拉普拉斯算子
进行卷积,从而突出了边缘信号。
(3)对边缘进行过零检测
(4)使用线性内插方法在亚像素分辨率水平上估计边缘的位置
由于对于后两步算法掌握不是很好,顾只对其进行阈值比较二值化
关键代码:
e=zeros(w,h);%?
?
?
?
e(find(y>s))=1;%?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
Y为图像,S为二值化阈值,经测试,选取S=120为较好效果。
下面是log算子边缘提取效果:
图2.2Log算子效果图
2.Canny算子边缘提取的实现
图象边缘检测必须满足两个条件:
一能有效地抑制噪声;二必须尽量精确确定边缘的位置。
根据对信噪比与定位乘积进行测度,得到最优化逼近算子。
这就是Canny边缘检测算子。
Canny算子边缘提取可分为4步
(1)用高斯滤波器平滑图象,去除噪声信号
关键代码:
k1=filter2(fspecial('gaussian'),I);%?
?
?
?
?
3*3?
?
?
?
?
?
?
其中,I为输入的图像
(2)用一阶偏导的有限差分来计算梯度的幅值和方向
关键代码:
%?
?
x,y?
?
?
?
?
?
?
?
?
?
fori=1:
m-1
forj=1:
n-1
Gx(i,j)=(k1(i,j+1)-k1(i,j)+k1(i+1,j+1)-k1(i+1,j))/2;%?
?
X?
?
?
Gy(i,j)=(k1(i,j)-k1(i+1,j)+k1(i,j+1)-k1(i+1,j+1))/2;%?
?
Y?
?
?
M(i,j)=sqrt(Gx(i,j)^2+Gy(i,j)^2);%?
?
?
?
?
a(i,j)=atan2(Gy(i,j),Gx(i,j));%?
?
?
?
?
?
end
End
(3)对梯度幅值进行非极大值抑制
仅仅得到全局的梯度并不足以确定边缘,因此为确定边缘,必须保留局部梯度最大的点,而抑制非极大值。
(non-maximasuppression,NMS)
解决方法:
利用梯度的方向。
图2.3非极大值抑制梯度方向示意图
四个扇区的标号为0到3,对应3*3邻域的四种可能组合。
在每一点上,邻域的中心象素M与沿着梯度线的两个象素相比。
如果M的梯度值不比沿梯度线的两个相邻象素梯度值大,则令M=0。
即:
关键代码:
%?
?
?
?
?
?
?
?
ifa(i,j)>=-1/8*pi&&a(i,j)<1/8*pi
a(i,j)=0;
......
elseif(a(i,j)>=5/8*pi&&a(i,j)<7/8*pi)||(a(i,j)>=-3/8*pi&&a(i,j)<-1/8*pi)
a(i,j)=3;
else
a(i,j)=0;
......
ifa(i,j)==0%?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
ifM(i,j) M1(i,j)=0; end end ...... ifa(i,j)==3 ifM(i,j) M1(i,j)=0; end end 经过非极大值抑制后效果如图2.4 图2.4非极大值抑制效果图 (4)用双阈值算法检测和连接边缘 双阈值算法。 双阈值算法对非极大值抑制图象作用两个阈值τ1和τ2,且2τ1≈τ2,从而可以得到两个阈值边缘图象N1[i,j]和N2[i,j]。 由于N2[i,j]使用高阈值得到,因而含有很少的假边缘,但有间断(不闭合)。 双阈值法要在N2[i,j]中把边缘连接成轮廓,当到达轮廓的端点时,该算法就在N1[i,j]的8邻点位置寻找可以连接到轮廓上的边缘,这样,算法不断地在N1[i,j]中收集边缘,直到将N2[i,j]连接起来为止。 设置高低阈值分别为40,20. 关键代码: %step4? ? ? ? ? ? M2=M1;%? ? ? ? ? ? ? ? ? ? M3=M1;%? ? ? ? ? ? ? ? ? ? x1=str2num(get(handles.lowgate,'string'));%? ? ? ? ? x2=str2num(get(handles.highgate,'string'));%? ? 高? ? fori=2: m-1 forj=2: n-1 ifM2(i,j) M2(i,j)=0; end ifM3(i,j) M3(i,j)=0; end end end 至此获得了2个高低阈值处理后的图像,如图2.5 图2.5双阈值图像 连接图像关键代码: CANNY=M3;%? ? ? ? ? ? ? fori=2: (m-2) forj=2: (n-2) ifisequal(CANNY(i,j),0) if~isequal(M2(i,j),0) if~isequal(sum(sum((CANNY(i-1: i+1,j-1: j+1)))),CANNY(i,j)) CANNY(i,j)=M2(i,j); end end end end end 连接后效果如图2.6 图2.6Canny算子效果图 3.几种滤波器的实现 (1)理想低通滤波器 一个理想的低通滤波器的传递函数由下式确定: H(u,v)=1D(u,v)≤D0 0D(u,v)>D0 式中 是从点(u,v)到频率平面的原点的距离。 D0为滤波器的截止频率。 关键代码: DD=sqrt((i-a0)^2+(j-b0)^2); ifDD<=D0%? ? ? ? H=1; else H=0; end f(i,j)=H*f(i,j); 设置截止频率为40,效果如图2.7 图2.7理想低通滤波效果图 (2)理想高通滤波器 一个二维理想高通滤波器的传递函数为 H(u,v)=0D(u,v)≤D0 1D(u,v)>D0 式中D(u,v),D0定义同上 关键代码: DD=sqrt((i-a0)^2+(j-b0)^2); ifDD>=D0%? ? ? ? H=1; else H=0; end f(i,j)=H*f(i,j); 设置截止频率为40,效果如图2.8 图2.8理想高通滤波效果图 (3)Buttonworth低通滤波器 一个N阶Butterworth低通滤波器的传递函数为: 式中n为阶数,其余同上 关键代码: DD=sqrt((i-a0)^2+(j-b0)^2); H=1/(1+(DD/D0)^(2*N)); f(i,j)=H*f(i,j); 设置截止频率为40,阶数为2,效果如图2.9 图2.9Butterworth低通滤波效果图 (4)Buttonworth高通滤波器 一个N阶Butterworth高通滤波器的传递函数为: 式中n为阶数,其余同上 关键代码: DD=sqrt((i-a0)^2+(j-b0)^2); H=1/(1+(D0/DD)^(2*N)); f(i,j)=H*f(i,j); 设置截止频率为40,阶数为2,效果如图2.10 图2.10Butterworth高通滤波效果图 (5)Guass低通滤波器 关键代码: DD=sqrt((i-a0)^2+(j-b0)^2); H=exp(-DD^2/(2*D0^2)); f(i,j)=H*f(i,j); 设置截止频率为40,效果如图2.11 图2.11Guass低通滤波效果图 (6)Guass高通滤波器 关键代码: DD=sqrt((i-a0)^2+(j-b0)^2); H=exp(-DD^2/(2*D0^2)); f(i,j)=(1-H)*f(i,j); 设置截止频率为40,效果如图2.12 图2.12Guass高通滤波效果图 4.数米粒的实现 数米粒算法思想分为以下几步: (1)将米粒在图像中提取出来: 先用Canny算子提取边缘,然后对米粒进行填充。 此时与边缘相交的米粒是不会填充的。 而后进行开启运算,平滑图像轮廓,消弱狭窄部分,去除与边界相交米粒。 关键代码: edI=edge(I,'canny'); fhI=imfill(edI,'hole');%? ? ? se=strel('disk',4);%? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? opI=imopen(fhI,se);%? 行? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? %? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 效果如图2.13 图2.13去除边缘米粒效果图 (2)对图像进行逐点扫描,若发现白点,这找出其所在连通域,并标示为已辨识过,米粒数+1.这样,下次扫描到这个连通域的另外的白点时,不会增加米粒数。 关键代码: ricearr=zeros(row*col,1); visited=zeros(size(opI)); fori=1: row forj=1: col ifopI(i,j)&&visited(i,j)==0%? ? ? ? ? ? ? ? ? ? +1,? ? ? ? ? ? ? ? ? ? if? ? ? ? count=count+1;%? ? ? +1 pt=[i;j]; [c,domain]=finddomain(opI,pt);%? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ricearr(count)=c;%? ? ? ? ? ? ? ? ? ? ? ? form=1: c visited(domain(1,m),domain(2,m))=1;%? ? ? ? ? ? ? ? ? ? ? ? ? ? ? end end end end 以上程序调用了finddomain函数,此函数作用是输入一幅图像,及图像中的一个点,返回图像中该点所在的连通域和连通域中像素的个数并输出。 通过这个函数,可以找出白点所在的连通域。 此函数算法为网上查找,因个人能力有限,对其算法不然完全理解。 具体代码见邮件中finddomain.m程序。 %? ? ? ? ? ? ? ? ? ? ? ? function[c,domain]=finddomain(input,pt) %? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? %? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? [rowcol]=size(input); iplus=[-1-101110-1]; jplus=[01110-1-1-1]; visited=zeros(size(input)); stack=zeros(2,row*col); domain=zeros(2,row*col); stack(: 1)=pt; domain(: 1)=pt; visited(pt (1),pt (2))=1;%? ? (pt (1),pt (2))? ? ? count=1; eos=1; sos=1; whilesos~=eos+1 m=stack(1,sos); n=stack(2,sos); fori=1: 8%? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ifm+iplus(i)>0&&m+iplus(i)<=row&&n+jplus(i)>0&&n+jplus(i)<=col ifinput(m+iplus(i),n+jplus(i))==1&&visited(m+iplus(i),n+jplus(i))==0 eos=eos+1; stack(: eos)=[m+iplus(i);n+jplus(i)]; visited(m+iplus(i),n+jplus(i))=1; count=count+1; domain(: count)=[m+iplus(i);n+jplus(i)]; end end end sos=sos+1;%? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1. end c=count; 数出米粒数为 5.图像展开的实现 图像展开的关键是要找到环形图的圆心及展开半径。 此段程序的思想是对原图进行二值化处理后标记连通域,从而找出中心点所在的连通域。 通过连通域属性结构体的Centroid变量获取重心坐标,从而确定了图像展开的原点。 利用原点即可求出展开半径。 关键代码如下: globalI; [iheight,iwidth,icolors]=size(I);%? ? ? ? ? ? ? ? ? ? Ibw=im2bw(I,0.43);%? ? ? ? ? IL=bwlabel(Ibw,8);%? ? ? ? ? ? ? ? 8? ? ? cSect_ID=IL(round(iheight/2),round(iwidth/2));%? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? cSect_stats=regionprops(IL);%? ? ? ? ? ? ? ? ? ? ? [ix,iy]=size(cSect_stats); ifcSect_ID==0%? ? ? ? ? ? ? ? ? icenterY=round(cSect_stats(round((ix+1)/2)).Centroid (1)); icenterX=round(cSect_stats(round((iy+1)/2)).Centroid (2)); else icenterY=round(cSect_stats(cSect_ID).Centroid (1));%? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? icenterX=round(cSect_stats(cSect_ID).Centroid (2)); end Rmax=min(min(icenterX,icenterY),min((iheight-icenterX),(iwidth-icenterY)));%? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 找到原点及展开半径后,对其进行极坐标展开,然后用线性加权内插法对个点颜色进行差值,恢复为彩色图像。 关键代码: Cout=360;%? ? ? ? ? ? ? ? 360? forrr=1: 1: Rmax%? ? ? ? ? ? ? ? ? ? ? ? ? ? forangl=1: 1: Cout Ixx=icenterX+(rr-1)*sind(360*angl/Cout); Iyy=icenterY+(rr-1)*cosd(360*angl/Cout); IPlr((Rmax-rr+1),(Cout-angl+1),1: 3)=round(...%? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (fix(Ixx)+1-Ixx)*(fix(Iyy)+1-Iyy)*I(fix(Ixx),fix(Iyy),1: 3)+... (fix(Ixx)+1-Ixx)*(Iyy-fix(Iyy))*I(fix(Ixx),fix(Iyy)+1,1: 3)+... (Ixx-fix(Ixx))*(fix(Iyy)+1-Iyy)*I(fix(Ixx)+1,fix(Iyy),1: 3)+... (Ixx-fix(Ixx))*(Iyy-fix(Iyy))*I(fix(Ixx)+1,fix(Iyy)+1,1: 3)); end end 效果如图2.14 图2.14图像展开效果图 可以看到展开效果并不是太好,因为圆心特征还不是特别明显,顾圆心连通域没找好。 要突出圆心特征,可在二值化后对图像进行多次收缩及扩张操作,然后再标记连通域,这样可以得到较准确的原点坐标。 然而由于时间原因没有继续完善程序。 6.分析图像方位的实现 (1)标出质心位置 算法思想: 首先用Log算子对回形针进行边缘提取 效果如图2.15 图2.15回形针Log边缘提取图 由于提取的是双边缘,对后续处理有些许繁琐。 之后进行填充操作,一个回形针包含3个连通域,全部都填充后表现为最外圈以内全部填充。 最后用标记连通域的函数bwlabel获取连通域及其个数,并依次在图上标出其质心。 关键代码: globalI; A=edge(I,'log'); %? ? ? ? M=imfill(A,'hole'); figure(4),imshow(M); title('? ? ? '); [IL,number]=bwlabel(M,8);%? ? ? ? ? ? ? ? 8? ? ? ? ? ? cSect_stats=regionprops(IL);%? ? ? ? ? ? ? ? ? ? Q=zeros(2,number); fori=1: number Q(1,i)=cSect_stats(i).Centroid (1);%? ? ? ? ? ? ? ? X? ? Q(2,i)=cSect_stats(i).Centroid (2);%? ? ? ? ? ? ? ? Y? ? end holdon; fori=1: number plot(Q(1,i),Q(2,i),'*r');%? ? ? ? ? ? ? ? ? ? ? ? ? end holdoff; 效果如图2.16 图2.16回形针质心图 可以看到中间的回形针外边缘提取不是很好,导致外边缘未填充,因而找出来3个连通域的中心。 左边的一个中心标到了外面去,不知道具体原因。 (2)标出角度信息 与质心算法相似,先用log算子提取出边缘,如何直接进行连通域标记,并标出每个连通域的中心。 关键代码: [IL,number]=bwlabel(A,8);%? ? ? ? ? ? ? ? 8? ? ? ? ? ? cSect_stats=regionprops(IL);%? ? ? ? ? ? ? ? ? ? Q=zeros(2,number); fori=1: number Q(1,i)=cSect_stats(i).Centroid (1);%? ? ? ? ? ? ? ? X? ? Q(2,i)=cSect_stats(i).Centroid (2);%? ? ? ? ? ? ? ? Y? ? end holdon; fori=1: number plot(Q(1,i),Q(2,i),'*y');%? ? ? ? ? ? ? ? ? ? ? ? ? end holdoff; 效果如图2.17 图2.16回形针角度图 可以看到一个回形针里包含了3个连通域,故每个回形针标出了3个中心,至于左边的单独点不知其出现原因。 后续处理思想: 把任意两个中心逐个连接,若连线经过大于2次白线,证明这2个中心不在同一回形针内;反之在同一回形针内,则计算出这条线的斜率即为回形针方向。 缺点是这样无法判断回形针头的方向。 3.学习心得和建议 通过学习工业摄影测量,让我感觉受益匪浅,这门课程的学习让我了解了图像处理的基础,是一门非常实用的课程。 并且我觉得唐老师这种独特的课程考核方式非常好,我们的学习模式不再是考前对书本知识的突击,而是在项目中为实现各种处理的功能,主动去了解这种处理的方式和处理后的效果,可以我们自己主动学习。 授课的过程中,唐老师还讲述有关这门课程的相关实际项目,加深我们对课程印象的同时,还让提高我们学习的积极性。 不过图像处理的算法对于我们大部分初学者来讲还是有些晦涩难懂,而且由于自己导师额外会安排任务,基本课下没有太多时间回去复习。 我就是作业开始得较晚,虽然加班加点,但几个程序都编的不完善,最后一问还没做完。 在此向唐老师致歉。 最后希望唐老师能降低下作业难度,祝唐老师工作顺利。 参考文献 [1]余成波.数字图像处理及MATLAB实现[M].重庆大学出版社,2002
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 基于 matlab 工业 摄影 作业 文库