watersed标志分水岭算法



/*====================================================================
函数名：            Watershed
功能：              用标记-分水岭算法对输入图像进行分割
算法实现：          无
输入参数说明：      OriginalImage --输入图像（灰度图，0~255）
SeedImage       --标记图像（二值图，0-非标记，1-标记）
LabelImage      --输出图像（1-第一个分割区域，2-第二个分割区域，...）
row             --图像行数
col             --图像列数
返回值说明：        无
====================================================================*/
#include"stdafx.h"
#include<vector>
#include<queue>
#include "Windows.h"
using namespace std;

void Watershed(unsigned char **OriginalImage, char** SeedImage, int **LabelImage, int row, int col)
{
 // using namespace std;

 //标记区域标识号，从1开始
 int Num = 0;
 int i, j;

 //保存每个队列种子个数的数组
 vector<int*> SeedCounts;
 //临时种子队列
 queue<POINT> quetem;
 //保存所有标记区域种子队列的数组，里面放的是种子队列的指针
 vector<queue<POINT>*> vque;

 int* array;
 //指向种子队列的指针
 queue<POINT> *pque;
 POINT temp;

 for (i = 0; i<row; i++)
 {
  for (j = 0; j<col; j++)
   LabelImage[i][j] = 0;
 }

 int m, n, k = 0;
 BOOL up, down, right, left, upleft, upright, downleft, downright;//8 directions...

 //预处理,提取区分每个标记区域，并初始化每个标记的种子队列
 //种子是指标记区域边缘的点，他们可以在水位上升时向外淹没（或者说生长）
 //pan's words：我的理解是梯度值较小的象素点，或者是极小灰度值的点。
 for (i = 0; i<row; i++)
 {
  for (j = 0; j<col; j++)
  {
   //如果找到一个标记区域
   if (SeedImage[i][j] == 1)
   {
    //区域的标识号加一
    Num++;
    //分配数组并初始化为零，表示可有256个灰阶
    array = new int[256];
    ZeroMemory(array, 256 * sizeof(int));
    //种子个数数组进vector，每次扫描则生成一个数组，并用区域标识号来做第一维。灰度级做第二维。
    //表示某个盆地区域中某灰阶所对应的点的数目。
    SeedCounts.push_back(array);
    //分配本标记号的优先队列，256个种子队列，
    //表示对应一个灰阶有一个队列，并且每个队列可以存储一个集合的点信息
    pque = new queue<POINT>[256];
    //加入到队列数组中，对应的是本标记号Num的
    vque.push_back(pque);
    //当前点放入本标记区域的临时种子队列中
    temp.x = i;
    temp.y = j;
    quetem.push(temp);
    //当前点标记为已处理
    LabelImage[i][j] = Num;
    SeedImage[i][j] = 127;//表示已经处理过

    //让临时种子队列中的种子进行生长直到所有的种子都生长完毕
    //生长完毕后的队列信息保存在vque中，包括区域号和灰阶，对应点数存储在seedcounts中
    while (!quetem.empty())
    {
     up = down = right = left = FALSE;
     upleft = upright = downleft = downright = FALSE;
     //队列中取出一个种子
     temp = quetem.front();
     m = temp.x;
     n = temp.y;
     quetem.pop();
     //注意到127对扫描过程的影响，影响下面的比较，但是不影响while语句中的扫描

     if (m>0)
     {
      //上方若为可生长点则加为新种子
      if (SeedImage[m - 1][n] == 1)
      {
       temp.x = m - 1;
       temp.y = n;
       quetem.push(temp);//如果这样的话，那么这些标记过的区域将再次在while循环中被扫描到，不会，因为值是127
       //新种子点标记为已淹没区域，而且是当前区域，并记录区域号到labelImage
       LabelImage[m - 1][n] = Num;
       SeedImage[m - 1][n] = 127;
      }
      else//否则上方为不可生长
      {
       up = TRUE;
      }
     }
     if (m>0 && n>0)
     {
      if (SeedImage[m - 1][n - 1] == 1)//左上方若为可生长点则加为新种子
      {
       temp.x = m - 1;
       temp.y = n - 1;
       quetem.push(temp);
       //新种子点标记为已淹没区域，即下一个循环中以127来标识不再扫描，而且是当前区域
       LabelImage[m - 1][n - 1] = Num;
       SeedImage[m - 1][n - 1] = 127;
      }
      else//否则左上方为不可生长
      {
       upleft = TRUE;
      }
     }

     if (m<row - 1)
     {
      if (SeedImage[m + 1][n] == 1)//下方若为可生长点则加为新种子
      {
       temp.x = m + 1;
       temp.y = n;
       quetem.push(temp);
       //新种子点标记为已淹没区域，而且是当前区域
       LabelImage[m + 1][n] = Num;
       SeedImage[m + 1][n] = 127;
      }
      else//否则下方为不可生长
      {
       down = TRUE;
      }
     }
     if (m<(row - 1) && n<(col - 1))
     {
      if (SeedImage[m + 1][n + 1] == 1)//下方若为可生长点则加为新种子
      {
       temp.x = m + 1;
       temp.y = n + 1;
       quetem.push(temp);
       //新种子点标记为已淹没区域，而且是当前区域
       LabelImage[m + 1][n + 1] = Num;
       SeedImage[m + 1][n + 1] = 127;
      }
      else//否则下方为不可生长
      {
       downright = TRUE;
      }
     }

     if (n<col - 1)
     {
      if (SeedImage[m][n + 1] == 1)//右方若为可生长点则加为新种子
      {
       temp.x = m;
       temp.y = n + 1;
       quetem.push(temp);
       //新种子点标记为已淹没区域，而且是当前区域
       LabelImage[m][n + 1] = Num;
       SeedImage[m][n + 1] = 127;
      }
      else//否则右方为不可生长
      {
       right = TRUE;
      }
     }
     if (m>0 && n<(col - 1))
     {
      if (SeedImage[m - 1][n + 1] == 1)//右上方若为可生长点则加为新种子
      {
       temp.x = m - 1;
       temp.y = n + 1;
       quetem.push(temp);
       //新种子点标记为已淹没区域，而且是当前区域
       LabelImage[m - 1][n + 1] = Num;
       SeedImage[m - 1][n + 1] = 127;
      }
      else//否则右上方为不可生长
      {
       upright = TRUE;
      }
     }

     if (n>0)
     {
      if (SeedImage[m][n - 1] == 1)//左方若为可生长点则加为新种子
      {
       temp.x = m;
       temp.y = n - 1;
       quetem.push(temp);
       //新种子点标记为已淹没区域
       LabelImage[m][n - 1] = Num;
       SeedImage[m][n - 1] = 127;
      }
      else//否则左方为不可生长
      {
       left = TRUE;
      }
     }
     if (m<(row - 1) && n>0)
     {
      if (SeedImage[m + 1][n - 1] == 1)//左下方若为可生长点则加为新种子
      {
       temp.x = m + 1;
       temp.y = n - 1;
       quetem.push(temp);
       //新种子点标记为已淹没区域
       LabelImage[m + 1][n - 1] = Num;
       SeedImage[m + 1][n - 1] = 127;
      }
      else//否则左方为不可生长
      {
       downleft = TRUE;
      }
     }

     //上下左右只要有一点不可生长，那么本点为初始种子队列中的一个
     //这里可否生长是由seedimage中的值来决定的。
     if (up || down || right || left ||
      upleft || downleft || upright || downright)
     {
      temp.x = m;
      temp.y = n;
      //下面这个矢量数组：第一维是标记号；第二维是该图像点的灰度级
      //m,n点对应的是while循环里面扫描的像素点。
      //Num是当前的区域号
      //这样这个二维信息就表示了，某个区域中对应某个灰度级对应的成员点的集合与个数
      //分别由下面两个量来表达
      vque[Num - 1][OriginalImage[m][n]].push(temp);//这两句中我把Num-1改成了Num...pan's codes...
      SeedCounts[Num - 1][OriginalImage[m][n]]++;
     }
    }//while结束，扫描到quetem为空而止。也就是对应所有的节点都得到不可生长为止（或者是周围的点要么不可生长，要么已生长）
   }//if结束
   // if(Num==5)
   // return;
  }
 }
 //在上述过程中，如果标记的点为0则表示，没有扫描到的点，或者表明不是输入的种子点
 //这里相当于是找seedimage传过来的初始区域的分水岭界线的所有的点；并且用标号记录每个区域，同时集水盆的边缘点进入队列。
 //上面是找集水盆的程序。同时也是连通区域。
 /*************************************/
 //test 这里测试一下剩下的非水盆地的点数。
 int seednum;
 for (i = 0; i<row; i++)
 {
  for (j = 0; j<col; j++)
  {
   if (SeedImage[i][j] == 0)
    seednum++;
  }
 }

 CString str;
 str.Format("pre region num:%d", Num);
 AfxMessageBox(str);
 /*************************************/
 bool actives;//在某一水位处，所有标记的种子生长完的标志
 int WaterLevel;

 //淹没过程开始，水位从零开始上升，水位对应灰度级，采用四连通法
 for (WaterLevel = 0; WaterLevel<180; WaterLevel++)//第二维。。。
 {
  actives = true;
  while (actives)
  {
   actives = false;
   //依次处理每个标记号所对应的区域，且这个标记号对应的区域的点的个数在SeedCounts里面
   for (i = 0; i<Num; i++)//第一维。。。
   {
    if (!vque[i][WaterLevel].empty())//对应的分水岭不为空集,i表示区域号，waterlevel表示灰阶
    {
     actives = true;
     while (SeedCounts[i][WaterLevel]>0)
     {
      SeedCounts[i][WaterLevel]--;//取出一个点，个数少一
      temp = vque[i][WaterLevel].front();//取出该区域的一个点，清空这个边缘点，表示当前
      //灰度级该像素已经处理掉了。
      vque[i][WaterLevel].pop();
      m = temp.x;
      n = temp.y;//当前种子的坐标
      if (m>0)
      {
       if (!LabelImage[m - 1][n])//上方若未处理，表示没有标号，应该在输入前已经作过初始化为0
        //本函数中在开头也作过初始化
       {
        temp.x = m - 1;
        temp.y = n;
        LabelImage[m - 1][n] = i + 1;//上方点标记为已淹没区域
        //注意到这个标记是与扫描点的区域号相同，一定在这个标号所属的区域吗？是的
        //这样在下一轮至少会扫描到这个点，确保不遗漏，但是下一轮的处理会使它合理
        //归类吗？问题还有这样标记并没有一定将它加入到种子队列。也就是说它
        //只是被淹没而不能向上淹没。只有满足下述可生长条件才行。
        if (OriginalImage[m - 1][n] <= WaterLevel)//上方若为可生长点则加入当前队列，当前高度的队列
        {
         vque[i][WaterLevel].push(temp);
        }
        else//否则加入OriginalImage[m-1][n]对应的灰度级的队列，为什么？
        {
         vque[i][OriginalImage[m - 1][n]].push(temp);
         SeedCounts[i][OriginalImage[m - 1][n]]++;
        }
       }
      }

      if (m<row - 1)
      {
       if (!LabelImage[m + 1][n])//下方若未处理
       {
        temp.x = m + 1;
        temp.y = n;
        LabelImage[m + 1][n] = i + 1;//下方点标记为已淹没区域

        if (OriginalImage[m + 1][n] <= WaterLevel)//下方若为可生长点则加入当前队列
        {
         vque[i][WaterLevel].push(temp);
        }
        else//否则加入OriginalImage[m+1][n]级队列
        {
         vque[i][OriginalImage[m + 1][n]].push(temp);
         SeedCounts[i][OriginalImage[m + 1][n]]++;
        }
       }
      }

      if (n<col - 1)
      {
       if (!LabelImage[m][n + 1])//右边若未处理
       {
        temp.x = m;
        temp.y = n + 1;
        LabelImage[m][n + 1] = i + 1;//右边点标记为已淹没区域

        if (OriginalImage[m][n + 1] <= WaterLevel)//右边若为可生长点则加入当前队列
        {
         vque[i][WaterLevel].push(temp);
        }
        else//否则加入OriginalImage[m][n+1]级队列
        {
         vque[i][OriginalImage[m][n + 1]].push(temp);
         SeedCounts[i][OriginalImage[m][n + 1]]++;
        }
       }
      }

      if (n>0)
      {
       if (!LabelImage[m][n - 1])//左边若未处理
       {
        temp.x = m;
        temp.y = n - 1;
        LabelImage[m][n - 1] = i + 1;//左边点标记为已淹没区域

        if (OriginalImage[m][n - 1] <= WaterLevel)//左边若为可生长点则加入当前队列
        {
         vque[i][WaterLevel].push(temp);
        }
        else//否则加入OriginalImage[m][n-1]级队列
        {
         vque[i][OriginalImage[m][n - 1]].push(temp);
         SeedCounts[i][OriginalImage[m][n - 1]]++;
        }
       }
      }
     }//while循环结束
     SeedCounts[i][WaterLevel] = vque[i][WaterLevel].size();
    }//if结束
   }//for循环结束
  }//while循环结束
 }//for循环结束

 /**********************************
 //test whether the origional segmentation num is changed...
 int nonzeronum=Num;
 for(m=0;m<Num;m++)
 {
 for(i=0;i<row;i++)
 {
 for(j=0;j<col;j++)
 {
 if(LabelImage[i][j]==m)
 break;
 }
 if(LabelImage[i][j]==m)
 break;
 }
 if(j==col-1 && i==row-1 && LabelImage[i][j]!=m)
 nonzeronum--;
 }
 str.Format("new region num:%d",nonzeronum);//这表示的是可扩张的点
 AfxMessageBox(str);
 /*********************************/
 while (!vque.empty())
 {
  pque = vque.back();
  delete[] pque;
  vque.pop_back();
 }
 while (!SeedCounts.empty())
 {
  array = SeedCounts.back();
  delete[] array;
  SeedCounts.pop_back();
 }
}

//在主程序中的调用代码如下(用于处理8位图)
/**************************************/
//get one image...
CAPCountMDoc *pDoc = GetDocument();
HDIB m_hDIBOptFlo = OnReadimage();//HDIB 是Dib图像句柄
LPSTR lpDIB2 = (LPSTR) ::GlobalLock((HGLOBAL)m_hDIBOptFlo);
int m_imageH = viewdib.DIBHeight(lpDIB2);
int m_imageW = viewdib.DIBWidth(lpDIB2);
// viewdib.MedianFilter(lpDIB2,m_imageW,m_imageH,3,3,1,1);//filter
LPSTR lpDIBBits = viewdib.FindDIBBits(lpDIB2);
DWORD DWidth = WIDTHBYTES(m_imageW * 8);
//watershed...
LPSTR lpDIBG = viewdib.MakeGray8Bmp(lpDIB2);
lpDIBG = viewdib.GetGrad(lpDIBG);
LPSTR lpDIBBitsG = viewdib.FindDIBBits(lpDIBG);
int i, k, num = 0;
int m = 0, n = 0;
//get oritional image
unsigned char** ori_image;
ori_image = alloc_uszchar_space(m_imageH, m_imageW);
for (i = 0; i<m_imageH; i++)
{
 for (k = 0; k<m_imageW; k++)
  ori_image[i][k] = *(lpDIBBitsG + i*DWidth + k);
}
//get seedimage and initialize it
char** seed_image;
seed_image = alloc_char_space(m_imageH, m_imageW);
for (i = 0; i<m_imageH; i++)
{
 for (k = 0; k<m_imageW; k++)
  seed_image[i][k] = 0;
}
unsigned char uszThre = 3;
unsigned char uszTem = 0;
for (i = 2, num = 0; i<m_imageH - 2; i++)
{
 for (k = 3; k<m_imageW - 3; k++)
 {
  uszTem = (BYTE)*(lpDIBBitsG + i*DWidth + k);
  if (uszTem<uszThre)
  {
   seed_image[i][k] = 1;
   num++;
  }
 }
}
//get labelimage, and initialize it
int** labelimage = alloc_int_space(m_imageH, m_imageW);
for (i = 0; i<m_imageH; i++)
{
 for (k = 0; k<m_imageW; k++)
  labelimage[i][k] = 0;
}
//call watershed..
viewdib.Watershed(ori_image, seed_image, labelimage, m_imageH, m_imageW);
//get the maximum number in labelimage[][]
int nregion = 0;
for (i = 0; i<m_imageH; i++)
{
 for (k = 0; k<m_imageW; k++)
 {
  if (nregion<labelimage[i][k])
   nregion = labelimage[i][k];
 }
}

/*看看算法，算法思想是这样的：
首先准备好山和初始的水。这山就是我们的初始图像了，比如用自然获取的图像的梯度来表征山地的每一点的高度吧；而初始的水就是在阈值记为Thre底下，所有的低于这个高度的整个山地都加水，直到这个阈值Thre高度。从而有三个初始量：unsigned char** Ori_image、char** Seed_image和int** Label_image。最后一个是为最终的结果做准备的。当然要做好初始化，比如，Ora_image赋值为原图像（256色灰度图）的梯度值，Seed_image则是初始状态下有水的置位，无水的复位，而Label_image则全初始化为0，最终得到的是各点对应的区域号。
接下来是考虑将已加的水进行记录，记录成连通的区域，也就是看看有多少个互不相关的集水盆，有五个，那么我们就涨出五个湖，而且尽可能的高，只要大家想到不溢出。在算法上，有多少个连通的区域就记录成多少个数据结构，工夫就在于如何将这些连通的区域连接成一块，并由一个数据结构来表达了。很好，我们准备用一个向量容器来实现初始保存，保存所有标记区域种子队列的数组，里面放的是种子队列的指针vector<queue<POINT>*> vque; ，而且这个队列是由一系列属于同一个区域的图像点组成，我们来自一个集水盆：）；其保存方式是这样的：queue<point> *pque = new queue<point>[256]；vque.push_back(pque)，这样便将一个成员放进到这个区域来了，即容器－－集水盆的保管都，容器中的每个指针，都指向一个集水盆，也就是我们要的连通区域；所以我们可以方便地由这个容器数据结构直接读值的方便性进行操作，一个脚标就可以得到一个区域（队列指针）的指针；而每个队列还不简单，并不是一列整形数那么易搞，所以说啊，这个算法，真头痛，这个队列的一个成员是一个点；而注意到vque里存放的一256个队列的的起始指针，真够残忍的。也就是说vque[m][n]就表达了一个队列，这个队列里可以存储操作一系列的点；显然容量取256是因为所有的初始或者是最终的区域中可能有0 - 256之间的不同的灰阶的点，那么我一个区域分用256个队列来记录这些成员点啦，很有可能，这里就只有一个集水盆，那么，256个灰阶的点都存在一个区域就有可能了
统计初始连通区域的方法是，八连通邻域法（还有其他方法：四连通域，六连通域），即从逐一扫描输入的Seed_image的每个像素点，将所有的标记了的初始集水盆一一纳入各自的区域，这是整修图像的扫描，形成外循环。先创建一个临时队列quetem，用来处理当前初始集水盆的连通连接，将逐一扫描到的属于一个特定的初始集水盆区域的可生长点暂存，并形成一个内循环。对当前扫描点的处理是，首先判断该点是否为某个初始集水盆的点，如果不是跳过；接下来是，如果是初始集水盆的点，那么它的八连通域中是否存在不可生长的点（这里的不可生长是指Seed_image中没有标记的点），扫描的八连通邻域中的点是可生长的，即有标记的，则将之加入到临时队列中quetem；如果扫描到的连通邻域中有不可生长的至少一个点存在，那么加入到种子队列，记当前区域号为Num，当前扫描点为(m，n)，从而当前的灰阶为Ori_image[m][n]，将当前点添加到种子队列：vque[Num－1][Ori_image[m][n]].push(POINT(m，n))。这里有两个循环，一个是quetem，另一个是Seed_image。直到两个循环完整结束，那么就得到了各个连通初始集水盆的记录，保存标记是区域号Num；而我们同时得到了初始的分水岭，那就放在了保存地点vque，这里面标识了它们对应的区域号，和区域里面对应的点的灰阶，即是特定区域特定灰阶对应的点的集合；我们可以获取这些分水岭的点所在的区域号，可以得到该区域的所有的灰阶的点信息。一句话，统计连通区域的功能有两个，一是标记初始区域，二是找分水岭
初始的区域标记好了，分岭也找到了，那么可以开始“水漫梁山”了。这就是淹没过程。淹没过程由也是由一个内嵌循环的循环来实现的：外循环是做水位上升（这里循环次数一定要256以内），waterlevel的上升，原来是已经做过了初始的水位分割，那么现在可以从Thre开始了，让水位慢慢上升，让它原本的湖慢慢扩张，尽量利用其应有的空间，而又不至于淹没到其它的邻居湖泊。内循环是扫描每个初始区域（当前Num，从而有Num个循环）的分水岭的点（在vque[][]中），按照给定的水位进行扩张。扩张过程是这样的：扫描到的分水岭的当前点，对其进行四连通邻域进行逐一检查，如果四连通域中有点没有标记的（那这一定是高度较高的点，较低的前面一定扫描过），那么先对该点以本区域号做标记Num（注意是当前的Num）；再判断它在当前水位下是否可生长（灰阶是否小于等于waterlevel），如果可生长那么加入到vque[Num][waterlevel]种子队列中，将会再次进入内循环，否则如果在当前水位下不可生长，则加入到这个邻域点的分水岭集合中vque[Num][Ori_image[邻域点]]队列中。如此往复循环，直到对应区域完成，一个水位，扫描所有的区域的分水岭，这样各自同时在一个水位下扩张，保证了不出现跳跃的情况出现（就是一个水位一个区域全局扩张）。
最终，所有的区域在每个水位都扩张完毕，得到了分割图，我们的大湖泊形成了
这是分水岭算法的一种实现方式。仔细考察不难发现这种实现方式不能产生新的集水盆，也就是说，由于初始集水盆的局限性，很可能会浪费大部分没有发掘出来的起始点较高的集水盆地。这样，我们就要对算法进行修改了。实现方式是：在淹没的过程中，我们是由阈值Thre的水位开始淹没的，那么我们可以对初始区域之外的没有标记的点（从Seed_image中考察），对之进行标记，条件是先把这一轮的内循环做好，然后在剩下的没标记区域中发掘新的集水盆，并加入到我们的种子队列中，下一个水位开始，就又多了一个新成员了，让它们不断膨胀，成长，拥有自己的小天的成员就会逐一的被分割出来。不过话说回来，我们这里是采用梯度图像，一般情况下，阈值初始分割能够满足我们的要求，把灰阶变化平滑的先截取出来，梯度信息已然足够强大；而如果采用了新盆地扩张，则比较适用于原始图像。*/
分水岭算法主要的分割目的在于找到图像的连通区域。利用梯度信息作为输入图像，会有一个矛盾点，如果对原始图像进行梯度计算时不作滤波平滑处理，很容易将物体分割成多个物体，那是因为噪声的影响；而如果进行滤波处理，又容易造成将某些原本几个的物体合成一个物体。当然这里的物体主要还是指图像变化不大或者说是灰度值相近的目标区域。