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

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

《中科大李厚强图像分析大作业.docx》由会员分享，可在线阅读，更多相关《中科大李厚强图像分析大作业.docx（16页珍藏版）》请在冰豆网上搜索。

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

区域提取与分析

——图像分析实验

SAHanxinFeng

1、实验目的

有一幅显微镜下获得的颗粒图，计算出每个粒子的面积、长轴、短轴。

用列表的形式给出计算结果。

实验目标很明确：

1、对图像中的连通片进行提取

2、统计出每个连通元的面积

3、计算长轴和短轴

语言C。

2、实验方法

1）图像二值化

阈值与结果有很大的关系，图像中一些细小的点很多，而且明暗变化于0~255之间，所以阈值过大时会排除掉很多点，程序中使用的折中的125。

2）区域标记

简单起见，考虑4-连通的情况。

本图中，4-连通与8-连通有一点差别，如下一块，着色不同表示不连通。

8连通时，左上角的小块跟大块是连通的。

实际标记时利用了两个模板，前者用于正向扫描，后者用于逆向扫描。

0A001A

B100B0

实验中，正向扫描先考察A，A有标记（无标记，默认是-1），则当前位置表为A_mark，否则，如果B有标记，则标为B_mark，否则，用新的计数；逆向扫描类似。

存储标记的节点如下：

typedefstructPNode{

charpixel;//像素

intmark;//标记

}PNode;

第一次正向扫描结果如下：

（用不同的颜色表示不同的计数）

部分放大：

逆向扫描结果如下：

还有六处没有正确标记的，其中一处如下：

这是可以理解的，比如如下的情况：

-1222222-1

-1-1111-1-1-1

由于逆向扫描时，先考虑右边像素再考虑下方像素，则蓝色下划线的“2”处仍然不变。

第二次正向扫描，得到正确的标定：

（即，每个连通片只有一种颜色）

3）整理标记，统计面积

由于标记本身并没有完全按照计数的顺序排列（只是同一连通片当中的标记相等而已），因此进行一次整理，同时维护一个存储面积和整理过后的计数（依次排列）的链表，如下：

typedefstructNUMNode{

intmark;//标记

intnum;//像素个数（面积）

NUMNode*next;

}NUMNode;

将2）结构中的pix数组中的计数处理后存入链表，由于同一连通片拥有一样的标记，根据标记进行查找，统计出各连通片的面积（包含的像素个数）。

4）长轴计算

重新定义一个结构体数组，每个结构体记录了一个连通片的边缘信息，结构如下：

typedefstructEDGEENode{

inten;//边缘点个数（周长），

intx[300];//边缘点的x坐标，数组长度只要比max{en}大就可以了

inty[300];//边缘点的y坐标

floatlo;//连通片长轴

floatsh;//连通片短轴

}EDGEENode;

维护这一数组是为了方便理解和处理。

数组的有效长度即连通片的个数（本图为30），首先根据计数值提取各个连通片的边缘像素的值以及个数（对于4-连通，边缘像素就是上下左右存在黑色像素），从而，可以将长轴与短轴的计算都局限在一个小范围中。

长轴的获取就是计算边缘点中相距最远的两点的距离，考虑到有很多连通片只有很少的几个像素组成，为了相互区别，因此使用的欧拉距离。

5）短轴计算

短轴的计算其实是一个很复杂的过程，按照定义，我的理解就是，在长轴的垂直方向上，距离最远的两个点的距离。

垂直用斜率k1*k2=-1表示，对于k1或者k2不存在或者为0的情况，单独考虑。

但是，关键问题在于，当长轴确定时，k1的值是精确的，但是计算得到的-1/k1与实际两点间的k2不可以能完全符合，或者说，连通片中甚至都不存在跟长轴完全垂直的两点，因此，短轴的选择是一个取舍问题。

1、如果连通片很大，则这一问题不明显，我们还可以按照取min{k1*k2+1}得到“最为”垂直的点，然后在这些点中寻找相距最远的两点。

2、但是如果连通片很小，则会出现很多细节问题，比如短轴甚至不存在，或者说即使是min{k1*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、参考资料

课程PPT

5、C代码

#include

#include

#include

typedefstructPNode{

charpixel;

intmark;

}PNode;

PNodepix[500][500];

typedefstructNUMNode{

intmark;

intnum;

NUMNode*next;

}NUMNode;

typedefstructEDGEENode{

inten;

intx[300];

inty[300];

floatlo;

floatsh;

}EDGEENode;

EDGEENodeEn[100];

voidmain（）{

inti,j,k,imax,imin,jmax,jmin;

IplImage*img=cvLoadImage（"d:

\\1.png"）;

intdn=img->nChannels;

intheight=img->height;

intwidth=img->width;

intwidth0=img->widthStep;

char*p=img->imageData;

imin=0;

jmin=0;

imax=width;

jmax=height;

intflag=0;

inti0=0;

inti1=0;

cvThreshold（img,img,125,255,CV_THRESH_BINARY）;

cvNamedWindow（"img1",CV_WINDOW_AUTOSIZE）;

for（j=jmin;j

for（i=imin;i

pix[j][i].pixel=*（p+width0*j+i*dn）;

pix[j][i].mark=-1;

}

}

//正向扫描

for（j=jmin;j

for（i=imin;i

if（pix[j][i].pixel==0）//黑点

continue;

i1++;

//如果是白点，则考虑AB位置是否有白点，有，则编号

if（j-1>=jmin&&pix[j-1][i].pixel==-1）{

pix[j][i].mark=pix[j-1][i].mark;

}

elseif（i-1>=imin&&pix[j][i-1].pixel==-1）{

pix[j][i].mark=pix[j][i-1].mark;

}

else{

//都黑

i0++;

pix[j][i].mark=i0;

}

}

}

for（j=jmin;j

for（i=imin;i

if（pix[j][i].pixel==-1）{

*（p+width0*j+i*dn）=（pix[j][i].mark*50）%255;

*（p+width0*j+i*dn+1）=（pix[j][i].mark*70）%255;

}

}

}

cvSaveImage（"d:

\\2.png",img）;

//逆向扫描

for（j=jmax;j>jmin;j--）{

for（i=imax;i>imin;i--）{

if（pix[j][i].pixel==0）//黑点

continue;

//如果是白点，则考虑AB位置是否有白点，有，则编号

if（i+1<=imax&&pix[j][i+1].pixel==-1）{

pix[j][i].mark=pix[j][i+1].mark;

}

elseif（j+1<=jmax&&pix[j+1][i].pixel==-1）{

pix[j][i].mark=pix[j+1][i].mark;

}

}

}

for（j=jmin;j

for（i=imin;i

if（pix[j][i].pixel==-1）{

*（p+width0*j+i*dn）=（pix[j][i].mark*50）%255;

*（p+width0*j+i*dn+1）=（pix[j][i].mark*70）%255;

}

}

}

cvSaveImage（"d:

\\3.png",img）;

//第二次正向扫描

for（j=jmin;j

for（i=imin;i

if（pix[j][i].pixel==0）//黑点

continue;

//如果是白点，则考虑AB位置是否有白点，有，则编号

if（j-1>=jmin&&pix[j-1][i].pixel==-1）{

pix[j][i].mark=pix[j-1][i].mark;

}

elseif（i-1>=imin&&pix[j][i-1].pixel==-1）{

pix[j][i].mark=pix[j][i-1].mark;

}

}

}

for（j=jmin;j

for（i=imin;i

if（pix[j][i].pixel==-1）{

*（p+width0*j+i*dn）=（pix[j][i].mark*50）%255;

*（p+width0*j+i*dn+1）=（pix[j][i].mark*70）%255;

}

}

}

cvSaveImage（"d:

\\4.png",img）;

//统计标志和面积

NUMNode*np=NULL;

NUMNode*np1,*np2;

flag=0;

intken;

for（j=jmin;j

for（i=imin;i

if（pix[j][i].pixel==0）//黑点

continue;

else{

np1=np;

np2=np;

flag=0;

if（np==NULL）{

np=newNUMNode;

np->mark=pix[j][i].mark;

np->num=1;

np->next=NULL;

continue;

}

while（np1!

=NULL）{

if（np1->mark==pix[j][i].mark）{

np1->num++;

flag=1;//已经有了

break;

}

np2=np1;

np1=np1->next;

}

if（flag==0）{//没有

np2->next=newNUMNode;

np2->next->mark=pix[j][i].mark;

np2->next->num=1;

np2->next->next=NULL;

}

}

}

}

//统计En数组，存储各连通元的边缘坐标与周长

for（j=0;j<100;j++）{

En[j].en=0;

}

for（j=jmin;j

for（i=imin;i

if（pix[j][i].pixel==0）//黑点

continue;

else{

np1=np;

np2=np;

flag=0;

while（np1!

=NULL）{

if（np1->mark==pix[j][i].mark）{

if（（j-1>jmin&&pix[j-1][i].pixel==0）||（i-1>imin&&pix[j][i-1].pixel==0）||（j+1

//周围有黑点，则说明是边界点

ken=En[np1->mark-1].en;

En[np1->mark-1].x[ken]=i;

En[np1->mark-1].y[ken]=j;

En[np1->mark-1].en++;

}

break;

}

np2=np1;

np1=np1->next;

}

}

}

}

np1=np;

i=1;

while（np1!

=NULL）{

En[i-1].en=En[np1->mark-1].en;

for（j=0;jmark-1].en;j++）{

En[i-1].x[j]=En[np1->mark-1].x[j];

En[i-1].y[j]=En[np1->mark-1].y[j];

}

np1->mark=i;

np1=np1->next;

i++;

}

i0=i-1;

//i-1即为连通元个数，也是En数组的长度

for（j=jmin;j

for（i=imin;i

if（pix[j][i].pixel==-1）{

*（p+width0*j+i*dn）=（pix[j][i].mark*50）%255;

*（p+width0*j+i*dn+1）=（pix[j][i].mark*70）%255;

}

}

}

cvShowImage（"img1",img）;

printf（"白像素个数：

%d\n",i1）;

intdistance,dmax=0,dx,dy;

intlx1,lx2,ly2,ly1,y,x;

floatk1,k2,dk2;

//利用En数组统计各连通元长轴

for（i=0;i

dmax=0;

for（j=0;j

for（k=0;k

dx=En[i].x[j]-En[i].x[k];

dy=En[i].y[j]-En[i].y[k];

if（（distance=dx*dx+dy*dy）>dmax）{

dmax=distance;

lx1=En[i].x[j];

lx2=En[i].x[k];

ly1=En[i].y[j];

ly2=En[i].y[k];

}

}

}

En[i].lo=（float）pow（（double）dmax,0.5）;

dmax=0;

if（En[i].lo==0）{

//长轴为0

En[i].sh=0;

continue;

}

elseif（lx2==lx1）{

//长轴为y轴方向

dmax=0;

for（j=0;j

y=En[i].y[j];

x=En[i].x[j];

for（k=0;k

if（En[i].y[k]==y）{

if（（dx=abs（En[i].x[k]-x））>dmax）

dmax=dx;

}

}

}

En[i].sh=dmax;

continue;

}

elseif（ly2==ly1）{

//长轴为x轴方向

dmax=0;

for（j=0;j

y=En[i].y[j];

x=En[i].x[j];

for（k=0;k

if（En[i].x[k]==x）{

if（（dx=abs（En[i].y[k]-y））>dmax）

dmax=dx;

}

}

}

En[i].sh=dmax;

continue;

}

else

k1=（float）（ly2-ly1）/（lx2-lx1）;

dmax=0;

k2=0;

dk2=（1/fabs（k1））*0.5;

for（j=0;j

for（k=0;k

if（k==j）

continue;

k2=（float）（En[i].y[k]-En[i].y[j]）/（En[i].x[k]-En[i].x[j]）;

//printf（"%0.1f",k1*k2）;

if（fabs（k2+1/k1）<=dk2）{

//k1*k2=-1，则垂直

dx=En[i].x[j]-En[i].x[k];

dy=En[i].y[j]-En[i].y[k];

if（（distance=dx*dx+dy*dy）>dmax）{

dmax=distance;

lx1=En[i].x[j];

lx2=En[i].x[k];

ly1=En[i].y[j];

ly2=En[i].y[k];

}

}

}

}

En[i].sh=（float）pow（（double）dmax,0.5）;

}

np1=np;

printf（"Mark\tEdge_n\tS_n\tlong\tshort\n"）;

while（np1!

=NULL）{

printf（"%d\t%d\t%d\t%0.1f\t%0.1f\n",np1->mark,En[np1->mark-1].en,np1->num,En[np1->mark-1].lo,En[np1->mark-1].sh）;

np1=np1->next;

}

cvWaitKey（0）;

}