影像處理（三）映像分割(1)Random Walks分割_影像處理

基於隨機遊走的映像分割演算法

基於隨機遊走的映像分割演算法是屬於圖論分割方法中的一種，這個演算法比較偏，網上的paper比較少，剛開始學習找個資料都不容易，其實這個演算法的原理就是通過求解一個鄰接矩陣方程組，跟三維空間三角網格曲面的調和場求解有點類似。

1、演算法開始前，先簡單描述一下隨機遊走模型

一維隨機遊走問題：設一個質點（隨機遊走者）沿著一條直線運動，單位時間內只能運動一個單位長度，且只能停留在該直線上的整數點，假設在時刻t，該質點位於直線上的點i，那麼在時刻t +1，該質點的位置有三種可能：①以p 的機率跳到整數點i-1，②或以q的機率跳到點i+1，③或以r=1-p-q的機率繼續停留在點i ，由於每一步的結果都是獨立的，且每種情況發生的機率之和都為1，則該過程服從伯努利分布，稱為貝努利隨機遊走過程。當 p=q=0.5時，即質點在下一時刻到達其相鄰點的機率是相等的，稱為簡單的隨機遊走。

例子1：如下圖所示，假設某一時刻一質點位於刻度2的位置，質點左右遊走的機率各為0.5，那麼下一時刻該質點既有可能往左走，也有可能往右走，當質點運動到位置0、5位置時，質點停止運動，求質點到最後運動到位置5的機率。該問題便是隨機遊走問題。

 

對於一維的簡單隨機遊走問題，滿足：

   

,    

其中，x為當前的位置點，x-1、x+1為位置x的左右鄰接頂點。根據該公式，我們可以列出由n個未知數組成的n個方程組，可以發現該方程組的係數矩陣即為拉普拉斯鄰接矩陣。拉普拉斯矩陣是非滿秩矩陣，需要添加邊界約束條件，方程組才有唯一解。

如例子1的問題，設添加邊界約束條件：

則最後可以列出如下方程組，求出各點到位置5的機率。

2、基於隨機遊走的映像分割演算法

①參考文獻：《Random Walks for Image Segmentation》

②文獻概述：隨機遊走演算法是一種基於圖論的分割演算法，屬於一種互動映像分割。它的分割思想是，以映像的像素為圖的頂點，相鄰像素之間的四鄰域或八鄰域關係為圖的邊，並根據像素屬性及相鄰像素之間特徵的相似性定義圖中各邊的權值，以此構建網狀圖，然後由通過使用者手工指定前景和背景標記，即前景物體和背景物體的種子像素，以邊上的權重為轉移機率，未標記像素節點為初始點，計算每個未標記節點首次到達各種子像素的機率，根據機率大小，劃分未標記節點，得到最終分割結果。

例子2：如圖下所示，圖中的小圓圈代表映像上的每個像素點。L1，L2，L3三個種子點分別由使用者互動輸入，作為標記的種子點。現要把映像分割成對應的三部分。

 

③演算法流程：

A.計算圖中任意一點vi與其各個鄰接頂點串連邊的權重:

其中，表示個像素點的灰階值、或紋理資訊等參數；

B.對於圖中任意一點vi的機率，其滿足隨機遊走機率公式：

其中，Ni為vi點的鄰接頂點（可為四鄰接頂點或八鄰接頂點），根據上式，可構建圖的拉普拉斯矩陣，然而拉普拉斯是非滿秩矩陣，需要添加邊界約束條件，才可根據方程組解出個各未知點的機率。也就是將映像分割問題轉換為Dirichlet問題進行求解。

C.添加邊界約束條件：以已標記的K類頂點作為邊界約束條件，求解未知點到各個類的機率。如下圖所示：求解各未知點遊走到L1的機率，則以，作為約束條件，可求得個未知點的機率，如下圖所示:

 

到達L1的機率

 

到達L2的機率

 

到達L3的機率

(5) 每一個未標記點，根據獲得的對 K 類標記的隸屬度值進行判斷，若未標記點到達第k類的機率最大，則將未標記節點vi判別為屬於類別k，完成分割。

最後貼一下自己寫的部分重要函數代碼：

[cpp]  view plain  copy //根據鄰接關係，構造拉普拉斯矩陣   void CRandomWalk::ComputeCoff()   {       int height=m_image->Height;       int width=m_image->Width;       m_A.resize(height*width,width*height);       int vn=height*width;       typedef Eigen::Triplet<double> Tri;       std::vector<Tri> tripletList;       for (int i=0;i<height;i++)       {           for (int j=0;j<width;j++)           {               int idex=i*width+j;               Eigen::Vector2i nei[4]={Eigen::Vector2i(i-1,j),Eigen::Vector2i(i,j-1),Eigen::Vector2i(i+1,j),Eigen::Vector2i(i,j+1)};               BYTE *data=GetpData(Eigen::Vector2i(i,j),m_image);               float sumw=0;               for (int k=0;k<4;k++)               {                   if (nei[k][0]>=0&&nei[k][0]<height&&nei[k][1]>=0&&nei[k][1]<width)                   {                       int idexnei=nei[k][0]*width+nei[k][1];                       BYTE *neidata=GetpData(nei[k],m_image);                       float w=-GetGrad(data,neidata);                       w=exp(w/(2*50*50));                       sumw+=w;                       tripletList.push_back(Tri(idex,idexnei,w));                   }                  }               //計算A               tripletList.push_back(Tri(idex,idex,-sumw));              }       }       m_A.setFromTriplets(tripletList.begin(),tripletList.end());       m_B.resize(height*width);       m_B.setZero();             }   void CRandomWalk::AddConstrain()   {          for (int i=0;i<m_front.size();i++)       {           int indexf=m_image->Width*m_front[i].y+m_front[i].x;           float a=m_A.coeff(indexf,indexf) +1;           m_A.coeffRef(indexf,indexf)=a;       }          for (int j=0;j<m_back.size();j++)       {           int indexb=m_image->Width*m_back[j].y+m_back[j].x;           float b=m_A.coeff(indexb,indexb) +1;           m_A.coeffRef(indexb,indexb)=b;       }          m_MatricesCholesky=new Eigen::SparseLU<Eigen::SparseMatrix<double>>(m_A);      }   void CRandomWalk::Solver()   {       ComputeCoff();       AddConstrain();             Eigen::VectorXd b=m_B;       for (int i=0;i<m_front.size();i++)       {           int indexf=m_image->Width*m_front[i].y+m_front[i].x;           b(indexf)+=1;       }       Eigen::VectorXd x=m_MatricesCholesky->solve(b);             Eigen::VectorXd bb=m_B;       for (int j=0;j<m_back.size();j++)       {           int indexfb=m_image->Width*m_back[j].y+m_back[j].x;           bb(indexfb)+=1;       }       Eigen::VectorXd y=m_MatricesCholesky->solve(bb);          //比較機率大小       for (int i=0;i<m_image->Height*m_image->Width;i++)       {           if (x(i)>y(i))           {               BYTE *data=(BYTE*)m_image->Scan0+i*4;               for (int k=0;k<3;k++)               {                   data[k]=255;               }           }       }       //int indexb=m_image->Width*m_front.y+m_front.x;      }   BYTE* CRandomWalk::GetpData(Eigen::Vector2i pt,BitmapData*image)   {       return (BYTE*)image->Scan0+image->Width*4*pt[0]+4*pt[1];      }   float CRandomWalk::GetGrad(BYTE*data1,BYTE*data2)   {       float sum0=0;       for (int i=0;i<3;i++)       {           sum0+=(data1[i]-data2[i])*(data1[i]-data2[i]);          }       return sum0;      }  

本文地址：http://blog.csdn.net/hjimce/article/details/45201263     作者：hjimce     聯絡qq：1393852684
更多資源請關注我的部落格：http://blog.csdn.net/hjimce                  原創文章，著作權，轉載請註明出處 。 參考文獻： 1、Random Walks for Image Segmentation 2、《隨機遊走映像分割演算法的研究》 這個演算法的效果感覺不是很好，所以分割效果就不貼了，具體可以看一下原版的文獻