|
|
Posted on 2007-04-30 00:01 cacar 阅读(448) 评论(0) 编辑 收藏 引用 所属分类: 图像
看别人的代码真的很辛苦，e文的资料，作者功底一般都比较nb，好在关键的地方总有些注释的。国人更nb，看过得很多源码根本就没有注释，只好自己边理解边写注释，费解死。
Mean Shift做图像分割的资料能找到的少之又少，看Rutgers上的论文，对我这个小本科来说还是有难度的。昨天上pudn淘了篇均值漂移分割图像的代码，只能做灰度图的，完全没有注释，想改一下却有最后分割的几句看不太懂，还没注释。。。
说明：下文注释不一定对哦，我的极其有限的理解，而已。
1 //MeanShift 分割滤波算法
2//pImg 图像数据指针
3//width height 宽度 高度
4//bUseFilter true - filter false - image segmentation
5//w - window size
6// h1 space variance
7// h2 range variance
8void MeanShiftFilterOrSegmentation(bool bUseFilter,BYTE *pImg,int width,int height,int w,int h1,double h2,BYTE *pResImg)
9   {
10 BYTE *pCur,*CovG;
11bool find,bOutOfImg;
12int count,fwin = 2*w+1,top=-1,size=width*height;
13int *G,*CovX,*CovY,x,y,i,j,new_x,new_y,new_g,old_x,old_y,old_g;
14double sum,grayDif,*PosWeight,*pWCur,numerator_x,numerator_y,numerator_g,denominator,g;
15//memory alloc
16G = new int[width*height];
17CovX = new int[width*height];
18CovY = new int[width*height];
19CovG = new BYTE[width*height];
20PosWeight = new double[fwin*fwin];
21//distance weight
22
23// 根据到当前点的距离设置窗口邻域内各点的权重
24// weight = exp(-0.5*(distence^2)/(h1^2));
25for(pWCur=PosWeight,sum=0.0,i=0;i<fwin;i++)
26  {
27for(j=0;j<fwin;j++,pWCur++)
28{
29*pWCur = exp(-0.5 * ( (i-w)*(i-w)+(j-w)*(j-w) )/(h1*h1) );
30 }
31}
32//Iteration Until Convergence
33// 逐行扫描,y为行号,pCur指向图像数据
34for(pCur=pImg,y=0;y<height;y++)
35{
36// 逐列扫描,x为列号
37for(x=0;x<width;x++,pCur++)
38{
39// old_x纪录当前像元列号
40old_x = x;
41// old_y记录当前像元行号
42old_y = y;
43// old_g指向当前像元灰度值
44old_g = *pCur;
45// 初始化计数器
46count=0;
47// 扫描图像越界标志，超出图像范围时为true
48bOutOfImg = false;
49while(count<1000)
50{
51count++;
52//
53numerator_x = numerator_y = numerator_g = denominator = 0.0;
54// 用i=old_y-w,i<old_y+w+1限定窗口行值范围，pWCur指向窗口权重数组
55for(pWCur=PosWeight,i=old_y-w;i<old_y+w+1;i++)
56{
57// 用j=old_x-w,j<old_x+w+1限定窗口列值范围,当j++时扫描当前行的
58// 下一个像元，同时pWCur指向权重数组中下一个值.
59// 当扫描完窗口内当前行后，由于列值j不再满足j<old_x+w+1条件，结束
60// 本层循环，由i++将扫描像元推向窗口的下一行
61for(j=old_x-w;j<old_x+w+1;j++,pWCur++)
62{
63//
64if(i<0 || i>height-1 || j<0 || j>width-1 ) continue;
65// *(pImg+i*width+j)是窗口内当前像元的灰度值，
66// grayDif为窗口内当前像元灰度值与pCur指向的作为窗口中心的
67// 当前点灰度值
68grayDif = (*(pImg+i*width+j) - old_g);
69// grayDif = grayDif^2;
70grayDif = grayDif*grayDif;
71// 用g表示当前窗口点与窗口中心样本点的差异（带权重）
72// g为平移向量 ？？？？？？
73// g为差异？？？？
74// exp(-0.5*grayDif/(h2*h2)) 高斯核？？？
75// g = *pWCur*exp(-0.5*grayDif/(h2*h2));
76// *pWCur为当前窗口位置的灰度值
77g = *pWCur*exp(-0.5*grayDif/(h2*h2));
78// numerator_x是列值平移向量的和
79numerator_x += j*g;
80// numerator_y是行值平移向量的和
81numerator_y += i*g;
82//
83numerator_g += *(pImg+i*width+j)*g;
84// denominator为平移向量（距离，权重~）的和
85denominator += g;
86}
87}
88new_x = (int)(numerator_x/denominator+0.5);
89new_y = (int)(numerator_y/denominator+0.5);
90new_g = (int)(numerator_g/denominator+0.5);
91if(new_y<0 || new_y>height-1 || new_x<0 || new_x>width-1 )
92{//超出图像范围
93new_x = old_x;
94new_y = old_y;
95new_g = old_g;
96bOutOfImg = true;
97break;
98}
99// 如果窗体移动距离<1.1 ???
100// 则结束迭代
101if( abs(new_x-old_x) + abs(new_y-old_y) < 1.1) break;
102// 变换窗口中心，继续迭代
103old_x = new_x;
104old_y = new_y;
105old_g = new_g;
106}
107if(count>999 ) printf("No Converge!\n");
108if(bOutOfImg) printf("Out Of Image!\n");
109// 将收敛点的列号(new_x),行号(new_y),灰度(new_g)赋给ConX,ConY,ConV中
110// 的相应点
111*(CovX + y*width + x) = new_x;
112*(CovY + y*width + x) = new_y;
113*(CovG + y*width + x) = new_g;
114}
115}
116//仅用于图像分割
117if(!bUseFilter)
118{
119for(i=0;i<size;i++)
120 {
121 for(find=false,j=0;j<=top;j++)
122 {
123 if( abs(*(CovG+i) - *(G+j))<h2 )
124 {
125 find = true;
126 *(CovG+i) = *(G+j);
127 break;
128 }
129 }
130 if(!find)
131 {
132 top++;
133 *(G+top) = *(CovG+i);
134 }
135 }
136}
137memcpy(pResImg,CovG,width*height);
138delete G;
139delete CovX;
140delete CovY;
141delete CovG;
142delete PosWeight;
143 }
144
145
EDISON的分割流程有个大概的印象了，Filter里面的一大块比较费力，希望可以尽快搞得定。
|