中科大李厚强图像分析大作业.doc

区域提取与分析 图像分析实验sa11010902 hanxin feng1、 实验目的有一幅显微镜下获得的颗粒图，计算出每个粒子的面积、长轴、短轴。用列表的形式给出计算结果。实验目标很明确：1、 对图像中的连通片进行提取2、 统计出每个连通元的面积3、 计算长轴和短轴语言c。2、 实验方法1） 图像二值化阈值与结果有很大的关系，图像中一些细小的点很多，而且明暗变化于0255之间，所以阈值过大时会排除掉很多点，程序中使用的折中的125。2） 区域标记简单起见，考虑4-连通的情况。本图中，4-连通与8-连通有一点差别，如下一块，着色不同表示不连通。8连通时，左上角的小块跟大块是连通的。实际标记时利用了两个模板，前者用于正向扫描，后者用于逆向扫描。0 a 0 0 1 ab 1 0 0 b 0实验中，正向扫描先考察a，a有标记（无标记，默认是-1），则当前位置表为a_mark，否则，如果b有标记，则标为b_mark，否则，用新的计数；逆向扫描类似。存储标记的节点如下：typedef struct pnodechar pixel;/像素int mark;/标记pnode;第一次正向扫描结果如下：（用不同的颜色表示不同的计数）部分放大：逆向扫描结果如下：还有六处没有正确标记的，其中一处如下：这是可以理解的，比如如下的情况：-1222222-1-1-1111-1-1-1由于逆向扫描时，先考虑右边像素再考虑下方像素，则蓝色下划线的“2”处仍然不变。第二次正向扫描，得到正确的标定：（即，每个连通片只有一种颜色）3） 整理标记，统计面积由于标记本身并没有完全按照计数的顺序排列（只是同一连通片当中的标记相等而已），因此进行一次整理，同时维护一个存储面积和整理过后的计数（依次排列）的链表，如下：typedef struct numnodeint mark;/标记int num;/像素个数（面积）numnode* next;numnode;将2）结构中的pix数组中的计数处理后存入链表，由于同一连通片拥有一样的标记，根据标记进行查找，统计出各连通片的面积（包含的像素个数）。4） 长轴计算重新定义一个结构体数组，每个结构体记录了一个连通片的边缘信息，结构如下：typedef struct edgeenodeint en;/边缘点个数（周长），int x300;/边缘点的x坐标，数组长度只要比maxen大就可以了int y300;/边缘点的y坐标float lo;/连通片长轴float sh;/连通片短轴edgeenode;维护这一数组是为了方便理解和处理。数组的有效长度即连通片的个数（本图为30），首先根据计数值提取各个连通片的边缘像素的值以及个数（对于4-连通，边缘像素就是上下左右存在黑色像素），从而，可以将长轴与短轴的计算都局限在一个小范围中。长轴的获取就是计算边缘点中相距最远的两点的距离，考虑到有很多连通片只有很少的几个像素组成，为了相互区别，因此使用的欧拉距离。5） 短轴计算短轴的计算其实是一个很复杂的过程，按照定义，我的理解就是，在长轴的垂直方向上，距离最远的两个点的距离。垂直用斜率k1*k2=-1表示，对于k1或者k2不存在或者为0的情况，单独考虑。但是，关键问题在于，当长轴确定时，k1的值是精确的，但是计算得到的-1/k1与实际两点间的k2不可以能完全符合，或者说，连通片中甚至都不存在跟长轴完全垂直的两点，因此，短轴的选择是一个取舍问题。1、 如果连通片很大，则这一问题不明显，我们还可以按照取mink1*k2+1得到“最为”垂直的点，然后在这些点中寻找相距最远的两点。2、 但是如果连通片很小，则会出现很多细节问题，比如短轴甚至不存在，或者说即使是mink1*k2+1对应的k2也显然不是垂直的关系。因此第二种方法是限定阈值，即fabs(k2+1/k1)=threshold，理想情况自然是threshold=0，但是问题是，差值选多少比较合适；|k1|大时，|k2|很小，反之亦然，所以按照百分比取阈值，实验中使用的(1/fabs(k1)*0.5;即50%。由于最终结果的好坏与否很难评定，因此实际上我也不知道这么选择是否合适，但是调整阈值大小可以看出变化，当阈值变小（比如用10%），则很多包含像素较少的连通片的短轴都会变为0。3、 结果分析及结论程序运行前，请将图像放在d:目录下，并且命名为1.png。运行后在d:目录下形成2.png,3.png,4.png，即上面贴出的几张图。实验结果如下：mark是标记计数，edge_n为周长，s_n为面积，long为长轴，short为短轴。由于只是实现功能，所以代码部分很简陋，也没有定义函数，结构体定义也在sd1sd2.cpp中，只有一个主函数，具体实现都在主函数中。基于vc6，用到了少量opencv的存储读取等功能函数。4、 参考资料课程ppt5、 c代码#include#include#includetypedef struct pnodechar pixel;int mark;pnode;pnode pix500500;typedef struct numnodeint mark;int num;numnode* next;numnode;typedef struct edgeenodeint en;int x300;int y300;float lo;float sh;edgeenode;edgeenode en100;void main()int i,j,k,imax,imin,jmax,jmin;iplimage* img=cvloadimage(d:1.png);int dn=img-nchannels;int height=img-height;int width=img-width;int width0=img-widthstep;char* p=img-imagedata;imin=0;jmin=0;imax=width;jmax=height;int flag=0;int i0=0;int i1=0;cvthreshold(img, img, 125, 255, cv_thresh_binary); cvnamedwindow(img1,cv_window_autosize);for(j=jmin;jjmax;j+)for(i=imin;iimax;i+)pixji.pixel=*(p+width0*j+i*dn);pixji.mark=-1;/正向扫描for(j=jmin;jjmax;j+)for(i=imin;i=jmin&pixj-1i.pixel=-1)pixji.mark=pixj-1i.mark;else if(i-1=imin&pixji-1.pixel=-1)pixji.mark=pixji-1.mark;else/都黑i0+;pixji.mark=i0;for(j=jmin;jjmax;j+)for(i=imin;ijmin;j-)for(i=imax;iimin;i-)if(pixji.pixel=0) /黑点continue;/如果是白点，则考虑ab位置是否有白点，有，则编号if(i+1=imax&pixji+1.pixel=-1)pixji.mark=pixji+1.mark;else if(j+1=jmax&pixj+1i.pixel=-1)pixji.mark=pixj+1i.mark;for(j=jmin;jjmax;j+)for(i=imin;iimax;i+)if(pixji.pixel=-1)*(p+width0*j+i*dn)=(pixji.mark*50)%255;*(p+width0*j+i*dn+1)=(pixji.mark*70)%255;cvsaveimage(d:3.png,img);/第二次正向扫描for(j=jmin;jjmax;j+)for(i=imin;i=jmin&pixj-1i.pixel=-1)pixji.mark=pixj-1i.mark;else if(i-1=imin&pixji-1.pixel=-1)pixji.mark=pixji-1.mark;for(j=jmin;jjmax;j+)for(i=imin;iimax;i+)if(pixji.pixel=-1)*(p+width0*j+i*dn)=(pixji.mark*50)%255;*(p+width0*j+i*dn+1)=(pixji.mark*70)%255;cvsaveimage(d:4.png,img);/统计标志和面积numnode* np=null;numnode* np1,*np2;flag=0;int ken;for(j=jmin;jjmax;j+)for(i=imin;imark=pixji.mark;np-num=1;np-next=null;continue;while(np1!=null)if(np1-mark=pixji.mark)np1-num+;flag=1; /已经有了break;np2=np1;np1=np1-next;if(flag=0) /没有np2-next=new numnode;np2-next-mark=pixji.mark;np2-next-num=1;np2-next-next=null;/统计en数组，存储各连通元的边缘坐标与周长 for(j=0;j100;j+)enj.en=0;for(j=jmin;jjmax;j+)for(i=imin;imark=pixji.mark)if(j-1jmin&pixj-1i.pixel=0)|(i-1imin&pixji-1.pixel=0)|(j+1jmax&pixj+1i.pixel=0)|(i+1mark-1.en;ennp1-mark-1.xken=i;ennp1-mark-1.yken=j;ennp1-mark-1.en+;break;np2=np1;np1=np1-next;np1=np;i=1; while(np1!=null)eni-1.en=ennp1-mark-1.en;for(j=0;jmark-1.en;j+)eni-1.xj=ennp1-mark-1.xj;eni-1.yj=ennp1-mark-1.yj;np1-mark=i;np1=np1-next;i+;i0=i-1;/i-1即为连通元个数，也是en数组的长度for(j=jmin;jjmax;j+)for(i=imin;iimax;i+)if(pixji.pixel=-1)*(p+width0*j+i*dn)=(pixji.mark*50)%255;*(p+width0*j+i*dn+1)=(pixji.mark*70)%255;cvshowimage(img1,img);printf(白像素个数：%dn,i1);int distance,dmax=0,dx,dy;int lx1,lx2,ly2,ly1,y,x;float k1,k2,dk2;/利用en数组统计各连通元长轴for(i=0;ii0;i+)dmax=0;for(j=0;jeni.en;j+)for(k=0;kdmax)dmax=distance;lx1=eni.xj;lx2=eni.xk;ly1=eni.yj;ly2=eni.yk;eni.lo=(float)pow(double)dmax,0.5);dmax=0;if(eni.lo=0)/长轴为0eni.sh=0;continue;else if(lx2=lx1)/长轴为y轴方向dmax=0;for(j=0;jeni.en;j+)y=eni.yj;x=eni.xj;for(k=0;kdmax)dmax=dx;eni.sh=dmax;continue; else if(ly2=ly1)/长轴为x轴方向dmax=0;for(j=0;jeni.en;j+)y=eni.yj;x=eni.xj;for(k=0;kdmax)dmax=dx;eni.sh=dmax;continue;elsek1=(float)(ly2-ly1)/(lx2-lx1);dmax=0;k2=0;dk2=(1/fabs(k1)*0.5;for(j=0;jeni.en;j+)for(k=0;keni.en;k+)if(k=j)continue;k2=(float)(eni.yk-eni.yj)/(eni.xk-eni.xj);/printf(%0.1f ,k1*k2);if(fabs(k2+1/k1)dmax)d