转自http://blog.csdn.net/tt2com/article/details/6321610
一直没看分水岭，vc图像处理的书对该算法的介绍也很少，cximage也找到相关的算法。就看了matlab何opencv。matlab对分水岭的算法是封装的也就没看源码，而opencv的写的实在太高升了，水平差在加没有过opencv，没办法看懂。网上找了几个算法的源码，一个老外用template写的，一个是国人写的。老外写的是NB，看不动，国人写的看了一些也是云里雾里的，最后么办法自己想办法写了个。
首先声明，这个只是简单的实现，仅仅针对一些简单的情况。效果还可以但有不足，一些地方的分割有可能将原本的单一前景给分割掉。和我一样的菜鸟，可以看看。高人路过可以指导。
先介绍下什么叫分水岭，
定义：
所谓图像分割是指将图像中具有特殊涵义的不同区域区分开来，这些区域是互相不交叉的，每一个区域都满足特定区域的一致性
本代码，没有完全的实现，只能满足现行对工作的需要。
分水岭的原理网上有很多，我仅仅说明下我对这个算法的简单理解：该算法可以解决针对，你可以用分割成前景，但有些干扰或者其他原因，分割出来的前景可以连成一片，无法区分那个前景是那个。这个时候，就可以用这个分水岭算法来帮忙。
（先看几张分水岭原理的图片，这个打字太麻烦，直接看图片吧）
算法的理解我花了半天多的时间才看明白。我的理解（可能有错误或者不足）就是M为灰度值最小（或者最大相对背景）为最低的盆地，其实也就是各个前景最为远离背景的地方（这个肯定有不足的理解，希望各位看客指导一下）。根据这个盆地我们不断的加水，使得盆地的水位不断的提升（也就是不断的靠近背景），在个提升地盘的时候，可能会出现和其他盆地相冲突的地方。这个时候我们需要本着和谐的态度来防止这种事情的发生，就是在冲突的地方建设堤坝。这个堤坝其实就是我们需要寻找的东西。
提升盆地水位的时候有三种可能
1：盆地上方实力，也就说这个区域的灰度和该盆地相邻，和其他势力没关系。这种情况，肯定是被吞并的。
2：盆地相邻实力，也就说这个区域的灰度不仅仅和该盆地相邻，也和其他的盆地有关系。这种情况，就是我们需要寻找的堤坝。
3：当水位提到一定高度的时候出现的新盆地，该地区和原始的盆地均不关。该水位点为该区域最低水位。这个时候就有新的盆地生成，也就是新的恶势力出现了。
针对上面的理解，我写了个依据图像距离变化的分水岭，适合一些小部分区域相交叉的情况。
首先我们需要对图像进行二值化，分割出感兴趣的位置区域。二值化，大家随便，只要能分割进行。
第二进行距离变换，代码来自一个vc图像算法原理和实现，大概这个书名吧
- //计算二值化图像的距离值
- bool processing::distance(BYTE *&pdata)
- {
- // 存储象素值的数组
- if (pdata!=NULL)
- {
- return false;
- }
- int * pnBinary = new int[m_size];
- int * pnStore = new int[m_size];
- pdata = new BYTE[m_size];
- // 将图象二值化
- for (int j = 0; j < m_height; ++j)
- {
- for(int i = 0; i < m_biwidth; ++i)
- {
- // 白色象素为背景，存成0
- if(m_pbinary[i+j*m_biwidth])
- {
- pnBinary[i+j*m_biwidth] = 0;
- pnStore [i+j*m_biwidth] = 0;
- }
- // 黑象素存成1
- else
- {
- pnBinary[i+j*m_biwidth] = 1;
- pnStore [i+j*m_biwidth] = 1;
- }
- }
- }
- int s = 1;
- while(s == 1)
- {
- s = 0;
- // 进行距离的累加
- for (int j = 1; j < m_height - 1; ++j)
- {
- for(int i = 1; i < m_biwidth - 1; ++i)
- {
- // 如果是背景，进行下一个循环
- if(pnBinary[i+j*m_biwidth]== 0)
- continue;
- // 如果当前象素值和四周的值一样，象素值增加一
- if(pnBinary[i+j*m_biwidth]==pnBinary[m_biwidth* (j-1) + i] && pnBinary[i+j*m_biwidth]==pnBinary[m_biwidth * (j+1) + i])
- if(pnBinary[i+j*m_biwidth]==pnBinary[m_biwidth * j + i-1] && pnBinary[i+j*m_biwidth]==pnBinary[m_biwidth* j + i+1])
- {
- pnStore[m_biwidth * j + i]++;
- s=1;
- }
- }
- }
- // 在进行下一轮循环前将当前的结果储存
- for (int j = 0; j < m_height; j++)
- for(int i = 1; i < m_biwidth; i++)
- pnBinary[m_biwidth * j + i] = pnStore[m_biwidth * j + i];
- }
- for (int j = 0; j < m_height; ++j)
- {
- for(int i = 0; i < m_biwidth; ++i)
- {
- //pdata[i+j*m_biwidth] =max(0,((25 - pnStore[j * m_biwidth + i]) * 10 + 5)); //
- pdata[i+j*m_biwidth]=max(0,255 - min(255,pnBinary[i+j*m_biwidth]));
- }
- }
- delete pnStore;
- delete pnBinary;
- return true;
- }
这个函数是写在一个类里面的所以，一些图像的基本参数为类的成员。对这个函数我就不再说明了。
- //分水岭处理二值化距离图像
- bool processing::Watershed(BYTE* Data,int * hist)
- {
- //修改成try机制
- if (Data==NULL||hist==NULL)
- {
- return false;
- }
- queue<POINT > *statistic = new queue<POINT>[255];//记录同级盆地队列，背景点不处理
- int * wdata = new int[m_size];
- memset(wdata,255,m_size*sizeof(int));
- bool cue = true;//变化标记
- int start_queue=0;//标记起始队列
- int lab_cout=256;//盆地标号
- POINT pt;
- const int WSHED = 255;//水坝
- //所有数据压入队列//边间点不进行处理
- for (int i=1;i<m_biwidth-1;++i)
- {
- for (int j=1;j<m_height-1;++j)
- {
- wdata[i+j*m_biwidth]=Data[i+j*m_biwidth];
- if (Data[i+j*m_biwidth]<255)
- {
- pt.x=i;
- pt.y=j;
- statistic[Data[i+j*m_biwidth]].push(pt);
- }
- }
- }
- //完成寻找第一低位盆地起始水位 代码说明1
- for(;start_queue<256&&hist[start_queue]==0;++start_queue);
- //完成寻找所有起始水位盆地 代码说明2
- while (hist[start_queue]>0)
- {
- pt=statistic[start_queue].front();
- statistic[start_queue].pop();
- wdata[pt.x+pt.y*m_biwidth] = lab_cout;
- --hist[start_queue];
- cue = true;
- while(cue)
- {
- cue = false;
- int size = (int)statistic[start_queue].size();
- for (int i=0;i<size;++i)
- {
- int t;
- pt=statistic[start_queue].front();
- statistic[start_queue].pop();
- //8个方向寻找lab
- t = wdata[pt.x+(pt.y-1)*m_biwidth];
- if (t>255)
- {
- wdata[pt.x+pt.y*m_biwidth]=lab_cout;
- --hist[start_queue];
- cue = true;
- continue;
- }
- t = wdata[pt.x-1+(pt.y-1)*m_biwidth];
- if (t>255)
- {
- wdata[pt.x+pt.y*m_biwidth]=lab_cout;
- --hist[start_queue];
- cue = true;
- continue;
- }
- t = wdata[pt.x-1+pt.y*m_biwidth];
- if (t>255)
- {
- wdata[pt.x+pt.y*m_biwidth]=lab_cout;
- --hist[start_queue];
- cue = true;
- continue;
- }
- t = wdata[pt.x-1+(pt.y+1)*m_biwidth];
- if (t>255)
- {
- wdata[pt.x+pt.y*m_biwidth]=lab_cout;
- --hist[start_queue];
- cue = true;
- continue;
- }
- t = wdata[pt.x+(pt.y+1)*m_biwidth];
- if (t>255)
- {
- wdata[pt.x+pt.y*m_biwidth]=lab_cout;
- --hist[start_queue];
- cue = true;
- continue;
- }
- t = wdata[pt.x+1+(pt.y+1)*m_biwidth];
- if (t>255)
- {
- wdata[pt.x+pt.y*m_biwidth]=lab_cout;
- --hist[start_queue];
- cue = true;
- continue;
- }
- t = wdata[pt.x+1+(pt.y)*m_biwidth];
- if (t>255)
- {
- wdata[pt.x+pt.y*m_biwidth]=lab_cout;
- --hist[start_queue];
- cue = true;
- continue;
- }
- t = wdata[pt.x+1+(pt.y-1)*m_biwidth];
- if (t>255)
- {
- wdata[pt.x+pt.y*m_biwidth]=lab_cout;
- --hist[start_queue];
- cue = true;
- continue;
- }
- statistic[start_queue].push(pt);//没有处理到继续压入队列
- }
- }
- ++lab_cout;
- }
- //完成寻找盆地和涨水过程//背景不处理
- for (int k=start_queue+1;k<255;++k)
- {
- while (hist[k]>0)
- {
- cue = true;
- while(cue)
- {
- cue = false;
- int size = (int)statistic[k].size();
- for (int i=0;i<size;++i)
- {
- int t;
- int lab=0;
- pt=statistic[k].front();
- statistic[k].pop();
- //8个方向寻找lab
- t = wdata[pt.x+(pt.y-1)*m_biwidth];
- if (t>255)
- {
- if (lab == 0) lab = t;
- else if( t != lab ) lab = WSHED;
- }
- t = wdata[pt.x-1+(pt.y-1)*m_biwidth];
- if (t>255)
- {
- if (lab == 0) lab = t;
- else if( t != lab ) lab = WSHED;
- }
- t = wdata[pt.x-1+pt.y*m_biwidth];
- if (t>255)
- {
- if (lab == 0) lab = t;
- else if( t != lab ) lab = WSHED;
- }
- t = wdata[pt.x-1+(pt.y+1)*m_biwidth];
- if (t>255)
- {
- if (lab == 0) lab = t;
- else if( t != lab ) lab = WSHED;
- }
- t = wdata[pt.x+(pt.y+1)*m_biwidth];
- if (t>255)
- {
- if (lab == 0) lab = t;
- else if( t != lab ) lab = WSHED;
- }
- t = wdata[pt.x+1+(pt.y+1)*m_biwidth];
- if (t>255)
- {
- if (lab == 0) lab = t;
- else if( t != lab ) lab = WSHED;
- }
- t = wdata[pt.x+1+(pt.y)*m_biwidth];
- if (t>255)
- {
- if (lab == 0) lab = t;
- else if( t != lab ) lab = WSHED;
- }
- t = wdata[pt.x+1+(pt.y-1)*m_biwidth];
- if (t>255)
- {
- if (lab == 0) lab = t;
- else if( t != lab ) lab = WSHED;
- }
- if (lab == WSHED)
- {
- --hist[k];
- wdata[pt.x+pt.y*m_biwidth]=WSHED;
- continue;
- }
- else
- {
- if (lab == 0)
- {
- statistic[k].push(pt);//没有处理到继续压入队列
- }
- if (lab>255)
- {
- wdata[pt.x+pt.y*m_biwidth]=lab;
- --hist[k];
- cue = true;
- }
- }
- }
- }
- if (!statistic[k].size())
- {
- pt=statistic[k].front();
- statistic[k].pop();
- wdata[pt.x+pt.y*m_biwidth] = ++lab_cout;;
- --hist[k];
- }
- }
- }
- for (int i=0;i<m_biwidth;++i)
- {
- for (int j=0;j<m_height;++j)
- {
- if (wdata[i+j*m_biwidth]>255)
- {
- Data[i+j*m_biwidth] = 0;
- }
- else
- Data[i+j*m_biwidth] = wdata[i+j*m_biwidth];
- }
- }
- delete []statistic;
- return true;
- }
这个分水岭有二个参数，第一是有二值得到的距离数据，第二个是这个距离数据的直方图统计。
对直方图的统计来决定循环的次数，当一个灰度被处理，如被水位吞没或者成为堤坝，均认为被处理。
说明1：首先找到直方图统计中最低的灰度做为起点水位
说明2：根据这个水位来寻找所以的盆地。此时的盆地就2中情况。一独立，二和相关的盆地相连接。由于是根据直方图统计来做的所以不知道是否连通。所以会出现如下情况，第一在8个方向能找到起始的盆地，那么我们将这个点淹没做为盆地，并将记录标记设置为真。二8方向没有盆地标记，那么其有可能是此刻盆地的势力范围，但暂时没有发展过来，或者和此刻的盆地无关。那么当标记为真的情况均处理完毕的时候，说明低水位的第一个盆地已经找完了。剩下的是这个水位下的其他盆地。此时我们增加一个盆地标记，开始寻找第二个盆地。如何循环直到所有的初始水位盆地均被找出来。
...（太晚了明天再接着写）
7696
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
