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

1、中科大李厚强图像分析大作业区域提取与分析 图像分析实验SA Hanxin Feng1、实验目的有一幅显微镜下获得的颗粒图，计算出每个粒子的面积、长轴、短轴。用列表的形式给出计算结果。实验目标很明确：1、对图像中的连通片进行提取2、统计出每个连通元的面积3、计算长轴和短轴语言C。2、实验方法1）图像二值化阈值与结果有很大的关系，图像中一些细小的点很多，而且明暗变化于0255之间，所以阈值过大时会排除掉很多点，程序中使用的折中的125。2）区域标记简单起见，考虑4-连通的情况。本图中，4-连通与8-连通有一点差别，如下一块，着色不同表示不连通。8连通时，左上角的小块跟大块是连通的。实际标记时利用了

2、两个模板，前者用于正向扫描，后者用于逆向扫描。0 A 0 0 1 AB 1 0 0 B 0实验中，正向扫描先考察A，A有标记（无标记，默认是-1），则当前位置表为A_mark，否则，如果B有标记，则标为B_mark，否则，用新的计数；逆向扫描类似。存储标记的节点如下：typedef struct PNode char pixel; /像素 int mark; /标记PNode;第一次正向扫描结果如下：（用不同的颜色表示不同的计数）部分放大：逆向扫描结果如下：还有六处没有正确标记的，其中一处如下：这是可以理解的，比如如下的情况：-1 2 2 2 2 2 2 -1-1 -1 1 1 1 -1 -1

3、 -1由于逆向扫描时，先考虑右边像素再考虑下方像素，则蓝色下划线的“2”处仍然不变。第二次正向扫描，得到正确的标定：（即，每个连通片只有一种颜色）3）整理标记，统计面积由于标记本身并没有完全按照计数的顺序排列（只是同一连通片当中的标记相等而已），因此进行一次整理，同时维护一个存储面积和整理过后的计数（依次排列）的链表，如下：typedef struct NUMNode int mark; /标记 int num; /像素个数（面积） NUMNode* next;NUMNode;将2）结构中的pix数组中的计数处理后存入链表，由于同一连通片拥有一样的标记，根据标记进行查找，统计出各连通片的面积（

4、包含的像素个数）。4）长轴计算重新定义一个结构体数组，每个结构体记录了一个连通片的边缘信息，结构如下：typedef struct EDGEENode int en; /边缘点个数（周长）， int x300; /边缘点的x坐标，数组长度只要比maxen大就可以了 int y300; /边缘点的y坐标 float lo; /连通片长轴 float sh; /连通片短轴EDGEENode;维护这一数组是为了方便理解和处理。数组的有效长度即连通片的个数（本图为30），首先根据计数值提取各个连通片的边缘像素的值以及个数（对于4-连通，边缘像素就是上下左右存在黑色像素），从而，可以将长轴与短轴的计算都

5、局限在一个小范围中。长轴的获取就是计算边缘点中相距最远的两点的距离，考虑到有很多连通片只有很少的几个像素组成，为了相互区别，因此使用的欧拉距离。5）短轴计算短轴的计算其实是一个很复杂的过程，按照定义，我的理解就是，在长轴的垂直方向上，距离最远的两个点的距离。垂直用斜率k1*k2=-1表示，对于k1或者k2不存在或者为0的情况，单独考虑。但是，关键问题在于，当长轴确定时，k1的值是精确的，但是计算得到的-1/k1与实际两点间的k2不可以能完全符合，或者说，连通片中甚至都不存在跟长轴完全垂直的两点，因此，短轴的选择是一个取舍问题。1、如果连通片很大，则这一问题不明显，我们还可以按照取mink1*k

6、2+1得到“最为”垂直的点，然后在这些点中寻找相距最远的两点。2、但是如果连通片很小，则会出现很多细节问题，比如短轴甚至不存在，或者说即使是mink1*k2+1对应的k2也显然不是垂直的关系。因此第二种方法是限定阈值，即fabs(k2+1/k1)=Threshold，理想情况自然是Threshold=0，但是问题是，差值选多少比较合适；|k1|大时，|k2|很小，反之亦然，所以按照百分比取阈值，实验中使用的(1/fabs(k1)*0.5;即50%。由于最终结果的好坏与否很难评定，因此实际上我也不知道这么选择是否合适，但是调整阈值大小可以看出变化，当阈值变小（比如用10%），则很多包含像素较少的

7、连通片的短轴都会变为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 PNode char p

8、ixel; int mark;PNode;PNode pix500500;typedef struct NUMNode int mark; int num; NUMNode* next;NUMNode;typedef struct EDGEENode int 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;

9、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+) p

10、ixji.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

11、=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

12、; 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)

13、=(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)

14、 /没有 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.xke

15、n=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;iim

16、ax;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) dma

17、x=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) /长轴为0 Eni.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; fo

18、r(j=0;jEni.en;j+) y=Eni.yj; x=Eni.xj; for(k=0;kdmax) dmax=dx; Eni.sh=dmax; continue; else k1=(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) dmax=distance; lx1=Eni.xj; lx2=Eni.xk; ly1=Eni.yj; ly2=Eni.yk; Eni.sh=(float)pow(double)dmax,0.5); np1=np; printf(MarktEdge_ntS_ntlongtshortn); while(np1!=NULL) printf(%dt%dt%dt%0.1ft%0.1fn,np1-mark,Ennp1-mark-1.en,np1-num,Ennp1-mark-1.lo,Ennp1-mark-1.sh); np1=np1-next; cvWaitKey(0);