分水岭算法代码C.doc

1、 弄了快一周了，分水岭算法代码终于出来了，这个算法太有难度了，通过自己敲一遍后，终于明白了为什么这么多人可以不厌其烦的该进这样经典的算法，他们对算 法其实根本就没有改良和创新，只是对图像进展一些预处理和改变一下微分算子，其实最核心的还是分水岭算法本身，到至今为止我还没有看到有那篇论文能对分水 岭算法本身有所改良，即便是效率上的进步也没有，甚至还没有发现有那个论文能把分水岭算法本身流程给描绘清楚也许没找到，通过这段时间的读论文，我发 现很多论文只有开头背景和结论，中间过程就写得很随意，我想一个算法结论固然重要，但没有过程拿来的结论，那只能认为这是拍脑壳拍出来的。下面是我这 几天用c+实现的分水岭

2、算法，该算法基于的是LUV颜色空间，假设有人要用，需要把RGB转化为LUV。/ ImageSegWatershed.cpp: implementation of the CImageSegWatershed class./#include ImageSegWatershed.h#include math.h#include qstring.h#include qqueue.husing namespace std;#define MAXV 256static const double RGB33= double(3.2405),double(-1.5371),double(-0.4985),

3、double(-0.9693),double(1.8760),double(0.0416), double(0.0556),double(-2040),double(1.0573);static const int XYZ33=4125,3576,1804,2125,7154,721,193,1192,9502;static const double Xn=double(0.9505);static const double Yn=double(1.0);static const double Zn=double(1.0888);static const double Un_prime=dou

4、ble(0.1978);static const double Vn_prime=double(0.4683);static const double Lt=double(0.008856);/ Construction/Destruction/CImageSegWatershed:CImageSegWatershed(int ImageHeight,int ImageWidth,int *ImageData)PicHeight=ImageHeight;PicWidth=ImageWidth;PicData=ImageData;luvData=new MyluvImageHeight*Imag

5、eWidth;RgbtoLuvPcm(ImageData,ImageWidth,ImageHeight,luvData);CImageSegWatershed:CImageSegWatershed()void CImageSegWatershed:WaterShedArith()long Piclen=PicHeight*PicWidth;int i,j,temp;/梯度模数组double *deltar=new doublePiclen;/梯度角度数组double *deltasita=new doublePiclen;/个点标识数组int *flag=new intPiclen;int *

6、gradientfre=new int256; /图像中各点梯度值的频数int *gradientadd=new int257; /各梯度起终位置memset(gradientfre,0,sizeof(int)*256);memset(gradientadd,0,sizeof(int)*257);/得到各点的梯度值GetGradient(deltar,deltasita);/以下统计各梯度频数；MyImageGraPt* graposarr = new MyImageGraPtPiclen;long xstart, imagepos, deltapos;xstart = imagepos =

7、deltapos = 0;for (int y=0; yPicHeight; y+)xstart = y*PicWidth;for (int x=0; x255)deltardeltapos = 255; /范围在0-255之间int tempi = (int)(deltardeltapos);gradientfretempi +; /灰度值频数；/统计各梯度的累加频数int added=0;gradientadd0=0; /第一个初始位置为0for(i=1;i256;i+)added+=gradientfrei-1;gradientaddi=added;gradientadd256=Picl

8、en; /最后累加和为总的像素点的个数memset(gradientfre,0,256*sizeof(int);MyImageGraPt *graposarr=new MyImageGraPtPiclen;/自左上至右下sorting. (注意是对所求得的各个像素的累加梯度梯度值进展排序)for (y=0; yPicHeight; y+)xstart = y*PicWidth;for (int x=0; xPicWidth; x+)deltapos = xstart + x;int tempi = (int)(deltardeltapos); /当前点的梯度值，由于前面的步骤，最大只能为255

9、0-255)；/根据梯度值决定在排序数组中的位置；int tempos = gradientaddtempi + gradientfretempi;gradientfretempi +;/梯度内指针后移(统计各个像素点的总的个数；graposarrtempos.gradient = tempi;/根据当前点的梯度将该点信息放后排序后数组中的适宜位置中去；graposarrtempos.x = x;graposarrtempos.y = y;int rgnumber=0;/记录分割后的区域数/下面是分水岭算法的核心步骤/开场吞没泛洪(泛洪函数各个形参的意义)/第一个参数是按一定规那么排好序的个

10、像素的梯度值和该像素在图片中的位置， /第二个参数是各个像素点的梯度累加和，倒数第二个参数是记录在泛洪过程中各个像素点的信息属于哪个区域，/最后一个参数是记录整个泛洪过程产生多少个区域FloodVincent(graposarr,gradientadd,0,255,flag,rgnumber); /区域增长步骤/以下准备计算各个区域的LUV的均值(也可以用其他的颜色空间)MyRgnInfo *rginfoarrr=new MyRgnInforgnumber+1;/分割后各个区的一些统计信息，第一个元素不用，图像各点所属区域的信息存放在flag数组中 /清空该数组(对各区域要统计的数据进展初始化

11、)for(int i=0;i=rgnumber;i+)rginfoarrri.isflag=FALSE;rginfoarrri.l=0;rginfoarrri.u=0;rginfoarrri.v=0;rginfoarrri.ptcount=0;long xstart;for(int y=0;yPicHeight;y+)xstart=y*PicWidth;for(int x=0;xPicWidth;x+)int pos=xstart+x;int rgid=flagpos;/当前位置点所属区域在区统计信息数组中的位置。/以下将该点的信息加到其所属区域信息中区rginfoarrrrgid.ptcou

12、nt+;rginfoarrrrgid.l+=luvDatapos.l;rginfoarrrrgid.u+=luvDatapos.u;rginfoarrrrgid.v+=luvDatapos.v;/求出各个区域的LUV均值for(i=0;i=rgnumber;i+)rginfoarrri.l=double(rginfoarrri.l)/rginfoarrri.ptcount);rginfoarrri.u=double(rginfoarrri.u)/rginfoarrri.ptcount);rginfoarrri.v=double(rginfoarrri.v)/rginfoarrri.ptcoun

13、t);/根据各个区域灰度均值和各区之间的邻接关系(用flag计算)进展区域交融(主要是消除不必要的过分割)int *mergearr=new int rgnumber+1;memset(mergearr,-1,sizeof(int)*(rgnumber+1);int mergednum=0;/第一个参数是各个区域的信息，第二个参数是图像所分的区域数，第三个参数是各个像素点的标志/倒数第二个参数和倒数第一个参数用于区域交融MergeRgs(rginfoarrr,rgnumber,flag,PicWidth,PicHeight,mergearr,mergednum);/确定合并后各像素点所属区域f

14、or(y=0;yPicHeight;y+)xstart=y*PicWidth;for(int x=0;xPicWidth;x+)int pos=xstart+x;int rgid=flagpos;/改点所属区域flagpos=FindMergedRgn(rgid,mergearr);delete mergearr;mergearr=NULL;/用各区均值来代替原像素点值double *luvbuff=NULL;luvbuff=new doublePicWidth*PicHeight*3;for(y=1;yPicHeight-1;y+)xstart=y*PicWidth;for(int x=1;

15、xPicWidth-1;x+)int pos=xstart+x;int rgid=flagpos;/改点所属区域luvDatapos.l=rginfoarrrrgid.l; /luvData用于全局保存luv值luvDatapos.u=rginfoarrrrgid.u;luvDatapos.v=rginfoarrrrgid.v;luvbuffpos*3=rginfoarrrrgid.l; /luvbuff用于专递参数给LuvToRgb();luvbuffpos*3+1=rginfoarrrrgid.u;luvbuffpos*3+2=rginfoarrrrgid.v;/求各像素的梯度值void

16、CImageSegWatershed:GetGradient(double *deltar,double *deltasita)/下面计算各像素在程度和垂直方向上的梯度,边缘点梯度计为0；int* deltaxarr;int* deltayarr;int grawidth = PicWidth;int graheight = PicHeight;int deltacount = grawidth * graheight;deltaxarr = new intdeltacount;deltayarr = new intdeltacount; /暂不计算边缘点；for (int y=1; ygra

17、height-1; y+)for (int x=1; xgrawidth-1; x+)int inarrpos = (y)*PicWidth + (x)*3 + 1;/在输入块中的位置；int deltaarrpos = y*grawidth + x;/在梯度数组中的位置；/卷积计算；deltaxarrdeltaarrpos = (int) ( (PicData(y-1)*PicWidth + (x+1)*3 + 1 /右上+ PicData(y)*PicWidth + (x+1)*3 + 1 /右+ PicData(y+1)*PicWidth + (x+1)*3 + 1 /右下- PicDa

18、ta(y-1)*PicWidth + (x-1)*3 + 1 /左上- PicData(y)*PicWidth + (x-1)*3 + 1 /左- PicData(y+1)*PicWidth + (x-1)*3 + 1 ) / 3 );/左下deltayarrdeltaarrpos = (int) ( (PicData(y-1)*PicWidth + (x+1)*3 + 1 /右上+ PicData(y-1)*PicWidth + (x)*3 + 1 /上+ PicData(y-1)*PicWidth + (x-1)*3 + 1 /左上- PicData(y+1)*PicWidth + (x-

19、1)*3 + 1 /左下- PicData(y+1)*PicWidth + (x)*3 + 1 /下- PicData(y+1)*PicWidth + (x+1)*3 + 1) / 3 );/右下/边缘赋为其内侧点的值；for (y=0; ygraheight; y+)int x1 = 0;int pos1 = y*grawidth + x1;deltaxarrpos1 = deltaxarrpos1+1;deltayarrpos1 = deltayarrpos1+1;int x2 = grawidth-1;int pos2 = y*grawidth + x2;deltaxarrpos2 =

20、deltaxarrpos2-1;deltayarrpos2 = deltayarrpos2-1;for (int x=0; xgrawidth; x+)int y1 = 0;int pos1 = x;int inner = x + grawidth;/下一行；deltaxarrpos1 = deltaxarrinner;deltayarrpos1 = deltayarrinner;int y2 = graheight-1;int pos2 = y2*grawidth + x;inner = pos2 - grawidth;/上一行；deltaxarrpos2 = deltaxarrinner;

21、deltayarrpos2 = deltayarrinner;for (y=0; ygraheight; y+)for (x=0; xgrawidth; x+)int temppos = y*grawidth + x;if ( (deltaxarrtemppos)=0 )if (deltayarrtemppos!=0)deltasitatemppos = 0;/程度方向;deltartemppos = (double) abs(deltayarrtemppos);elsedeltasitatemppos = -1;/无确定方向;deltartemppos = (double) abs(delt

22、ayarrtemppos);continue;deltasitatemppos = (double) ( atan(double)deltayarrtemppos/ (double)deltaxarrtemppos ) + PI/2. );deltartemppos = (double) sqrt(double)( deltayarrtemppos*deltayarrtemppos+ deltaxarrtemppos*deltaxarrtemppos ) );delete deltaxarr;deltaxarr = NULL; /删除程度和垂直梯度数组；delete deltayarr;del

23、tayarr = NULL;void CImageSegWatershed:FloodVincent(MyImageGraPt *imiarr,int *graddarr,int minh,int maxh,int *flagarr,int &outrgnumber) const int INIT = -2;const int MASK = -1;const int WATERSHED = 0;int h = 0;int imagelen = PicWidth * PicHeight;for (int i=0; iimagelen; i+)flagarri = INIT;int* imd =

24、new intimagelen;/间隔 数组，直接存取；for (i=0; iimagelen; i+)imdi = 0;/保存分水岭的队列std:queue myqueue;int curlabel = 0;/各盆地标记； /下面是从在imiarr记录中的像素从开场泛起 for (h=minh; h=maxh; h+)int stpos = graddarrh;int edpos = graddarrh+1;for (int ini=stpos; ini=0)if (flagarrleft=0)imdipos = 1;myqueue.push(ipos);/点位置压入fifo;continu

25、e;int right = ipos + 1;if (x+1=0)imdipos = 1;myqueue.push(ipos);/点位置压入fifo;continue;int up = ipos - PicWidth;if (y-1=0)if (flagarrup=0)imdipos = 1;myqueue.push(ipos);/点位置压入fifo;continue;int down = ipos + PicWidth;if (y+1=0)imdipos = 1;myqueue.push(ipos);/点位置压入fifo;continue; /以下根据先进先出队列扩展现有盆地；int cur

26、dist = 1;myqueue.push(-99);/特殊标记；while (TRUE)int p = myqueue.front();myqueue.pop();if (p = -99)if ( myqueue.empty() )break;elsemyqueue.push(-99);curdist = curdist + 1;p = myqueue.front();myqueue.pop(); A00j=a00*j;A01j=a01*j;A02j=a02*j; A10j=a10*j;A11j=a11*j;A12j=a12*j; A20j=a20*j;A21j=a21*j;A22j=a22

27、j; double *my_pow=new doubleMAXV; for(j=0;jLt) l_star=my_pow(int)(tval*255+0.5); else l_star=double(903.3*tval); my_temp=x+15*y+3*z; if(my_temp) u_prime=double(x0)int pos = tempstr.Find( );ASSERT(pos=0);CString idstr = tempstr.Left(pos);tempstr.Delete(0, pos+1);int idint = (int) strtol(idstr, NULL,

28、 10);/判断该区是否已被合并，假设是，那么一直找到该区当前的区号；idint = FindMergedRgn(idint, mergearr);if (idint=curid)continue;/这个邻区已被合并到当前区，跳过；double tl, tu, tv;tl = rginfoarridint.l;/当前处理的邻区的LUV值;tu = rginfoarridint.u;tv = rginfoarridint.v;double tempdis = pow(tl-cl, 2)+ pow(tu-cu, 2) + pow(tv-cv, 2);if (tempdisNearMeasureBias)MergeTwoRgn(curid, idint, neiarr, rginfoarr, mergearr);cl = rginfoarrcurid.l;/当前区的LUV值刷新；cu = rginfoarrcurid.u;cv = rginfoarrcurid.v;loopmerged = TRUE;void CImageSegWatershed:LuvToRgb(double *LuvData,int PicWidth,int PicHeight,int *PicData)int i,j;int r,g,b;12 / 12