**[研究內容]** 二值圖像距離變換 **[正文]** 二值圖像距離變換的概念由Rosenfeld和Pfaltz于1966年在論文[1]中提出，目前廣泛應用于計算機圖形學，目標識別及GIS空間分析等領域，其主要思想是通過表識空間點(目標點與背景點)距離的過程，最終將二值圖像轉換為灰度圖像。 距離變換按照距離的類型可以分為歐式距離變換(Eudlidean Distance Transfrom)和非歐式距離變換兩種，其中，非歐式距離變換又包括棋盤距離變換(Chessboard Distance Transform)，城市街區距離變換(Cityblock Distance Transform)，倒角距離變換(Chamfer Distance Transform)等； 距離變換的主要過程： 假設一幅二值圖像I，包含一個連通區域S，其中有目標集O和背景集B，距離圖為D，則距離變換的定義如公式1-(1)： ??????????????????????????????????????? 1-(1) 具體步驟如下： 1，將圖像中的目標像素點分類，分為內部點，外部點和孤立點。 以中心像素的四鄰域為例，如果中心像素為目標像素(值為1)且四鄰域都為目標像素(值為1)，則該點為內部點。如果該中心像素為目標像素，四鄰域為背景像素(值為0)，則該中心點為孤立點，如下圖Fig.1所示。除了內部點和孤立點之外的目標區域點為邊界點。  2，計算圖像中所有的內部點和非內部點，點集分別為S1，S2。 3，對于S1中的每一個內部點(x,y)，使用距離公式disf()計算騎在S2中的最小距離，這些最小距離構成集合S3。 4，計算S3中的最大最小值Max,Min。 5，對于每一個內部點，轉換后的灰度值G計算如下公式1-(2)所示： ??? 1-(2) 其中，S3(x,y)表示S1中的點(x,y)在S2中的最短距離 6，對于孤立點保持不變。 在以上距離變換的過程中，距離函數disf()的選取如果是歐式距離，則該距離變換稱為歐式距離變換，依次類推。對于距離的求取，目前主要的距離公式如下： 歐式距離： ???? 1-(3) 棋盤距離： ?????????????? 1-(4) 城市街區距離： ? 1-(5) 對于歐式距離變換，由于其結果準確，而計算相比非歐式距離變換較為復雜，因此，出現了較多的快速歐式距離變換算法，這里筆者介紹一種基于3*3模板的快速歐式距離變換算法(文獻2)，具體過程如下：? 1，按照從上到下，從左到右的順序，使用模板如圖Fig.2，依次循環遍歷圖像I，此過程稱為前向循環。  對于p對應的像素(x,y)，我們計算五個距離：d0,d1,d2,d3,d4： ????d0=p(x,y) ????d1=p(x-1,y)+disf((x-1,y),(x,y)) ????d2=p(x-1,y-1)+disf((x-1,y-1),(x,y)) ????d3=p(x,y-1)+disf((x,y-1),(x,y)) ????d4=p(x+1,y-1)+disf((x-1,y-1),(x,y)) ????則p(x,y)變換后的像素值為： ????p(x,y)=Min(d0,d1,d2,d3,d4); ????使用上述算法得到圖像I'。 ????2，按照從下到上，從右到左的順序，使用Fig.2所示模板依次循環遍歷圖像I’，此過程稱為后向循環。 對于p對應的像素(x,y)，我們計算五個距離：d0,d5,d6,d7,d8： ????d0=p(x,y) ????d5=p(x+1,y)+disf((x+1,y),(x,y)) ????d6=p(x+1,y+1)+disf((x+1,y+1),(x,y)) ????d7=p(x,y+1)+disf((x,y+1),(x,y)) ????d8=p(x-1,y+1)+disf((x-1,y+1),(x,y)) ????則p(x,y)后向變換后的像素值為： ????p(x,y)=Min(d0,d5,d6,d7,d8); ????使用上述算法得到的圖像即為距離變換得到的灰度圖像。 ????以上過程即文獻2所示快速歐式距離變換算法。如果我們需要非歐氏距離變換的快速算法，只需要修改文獻2算法中的歐式距離公式disf()為非歐式距離公式，如棋盤距離，城市街區距離等，過程依次類推。 對于歐式距離變換算法，相關學者研究了速度更快的倒角距離變換算法，來近似歐式距離變換的效果。具體過程如下： ????1，使用前向模板如圖Fig.3中左邊3*3模板，對圖像從上到下，從左到右進行掃描，模板中心0點對應的像素值如果為0則跳過，如果為1則計算模板中每個元素與其對應的像素值的和，分別為Sum1,Sum2,Sum3,Sum4，Sum5，而中心像素值為這五個和值中的最小值。 ????2，使用后向模板如圖Fig.3中右邊的3*3模板，對圖像從下到上，從右到左進行掃描，方法同上。 ????3，一般我們使用的倒角距離變換模板為3*3和5*5，分別如下圖所示：  [實驗結果] ????實驗采用512*512大小的圖像進行測試，測試PC為64位，Intel(R)?Core(TM)?i5-3230?CPU,?2.6GHz,?4G?RAM，運行環境為VS2008，C#。 ????實驗結果如下：  ???????????????????????????????????????????????????????????????????? (a)原圖  ??????????????????????????????????????????? (b)Euclidean Distance Transfrom  ?????????????????????????????????????????????????? (c) Cityblock??Distance Transfrom  ????????????????????????????????????????????? (d) Chessboard Distance Transform  ????????????????????????????????????????????? (e)?Chamfer?Distance?Transform 對于以上歐式距離變換與非歐式距離變換，我們做了時間分析，結果如下：  ?對于Table?1的數據，是通過計算50張512*512大小的圖像得到的平均結果，代碼未曾優化，距離變換結果均做了均衡化處理，對于不同配置，不同程序語言可能存在一定差異，總體而言，基于3*3模板的倒角距離變換速度最快，大概是歐氏距離快速算法的一半。 [參考文獻] [1]?Rosenfeld?A,PfaltzJ.L,?Sequential?operations?in?digital?pic?ture?processing.?Journal?of?ACM,1966,?13(4):471-494. [2]?Frank?Y.Shih,Yi-Ta?Wu,?Fast?Euclidean?distance?transformation?in?two?scans?using?a?3*3?neighborhood.?Journal?of?Computer?Vision?and?Image?Understanding?2004,195–205. 附錄 本人使用C#編寫的代碼如下： 本人的完整demo，包含參考文獻，測試圖像，地址：[http://download.csdn.net/detail/trent1985/6841125](http://download.csdn.net/detail/trent1985/6841125) ~~~ /// <summary> /// Distance transform for binary image. /// </summary> /// <param name="src">The source image.</param> /// <param name="distanceMode">One parameter to choose the mode of distance transform from 1 to 3, 1 means Euclidean Distance, 2 means CityBlock Distance, 3 means ChessBoard Distance.</param> /// <returns>The result image.</returns> public Bitmap DistanceTransformer(Bitmap src, int distanceMode) { int w = src.Width; int h = src.Height; double p0 = 0, p1 = 0, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0, p7 = 0, p8 = 0, min = 0; int mMax = 0, mMin = 255; System.Drawing.Imaging.BitmapData srcData = src.LockBits(new Rectangle(0, 0, w, h), System.Drawing.Imaging.ImageLockMode.ReadOnly, System.Drawing.Imaging.PixelFormat.Format24bppRgb); unsafe { byte* p = (byte*)srcData.Scan0.ToPointer(); int stride = srcData.Stride; byte* pTemp; for (int y = 0; y < h; y++) { for (int x = 0; x < w; x++) { if (x > 0 && x < w - 1 && y > 0 && y < h - 1) { p0 = (p + x * 3 + y * stride)[0]; if (p0 != 0) { pTemp = p + (x - 1) * 3 + (y - 1) * stride; p2 = pTemp[0] + GetDistance(x, y, x - 1, y - 1, distanceMode); pTemp = p + x * 3 + (y - 1) * stride; p3 = pTemp[0] + GetDistance(x, y, x, y - 1, distanceMode); pTemp = p + (x + 1) * 3 + (y - 1) * stride; p4 = pTemp[0] + GetDistance(x, y, x + 1, y - 1, distanceMode); pTemp = p + (x - 1) * 3 + y * stride; p1 = pTemp[0] + GetDistance(x, y, x - 1, y, distanceMode); min = GetMin(p0, p1, p2, p3, p4); pTemp = p + x * 3 + y * stride; pTemp[0] = (byte)Math.Min(min, 255); pTemp[1] = (byte)Math.Min(min, 255); pTemp[2] = (byte)Math.Min(min, 255); } } else { pTemp = p + x * 3 + y * stride; pTemp[0] = 0; pTemp[1] = 0; pTemp[2] = 0; } } } for (int y = h - 1; y > 0; y--) { for (int x = w - 1; x > 0; x--) { if (x > 0 && x < w - 1 && y > 0 && y < h - 1) { p0 = (p + x * 3 + y * stride)[0]; if (p0 != 0) { pTemp = p + (x + 1) * 3 + y * stride; p5 = pTemp[0] + GetDistance(x, y, x + 1, y, distanceMode); pTemp = p + (x + 1) * 3 + (y + 1) * stride; p6 = pTemp[0] + GetDistance(x, y, x + 1, y + 1, distanceMode); pTemp = p + x * 3 + (y + 1) * stride; p7 = pTemp[0] + GetDistance(x, y, x, y + 1, distanceMode); pTemp = p + (x - 1) * 3 + (y + 1) * stride; p8 = pTemp[0] + GetDistance(x, y, x - 1, y + 1, distanceMode); min = GetMin(p0, p5, p6, p7, p8); pTemp = p + x * 3 + y * stride; pTemp[0] = (byte)Math.Min(min, 255); pTemp[1] = (byte)Math.Min(min, 255); pTemp[2] = (byte)Math.Min(min, 255); mMax = (int)Math.Max(min, mMax); } } else { pTemp = p + x * 3 + y * stride; pTemp[0] = 0; pTemp[1] = 0; pTemp[2] = 0; } } } mMin = 0; for (int y = 0; y < h; y++) { for (int x = 0; x < w; x++) { pTemp = p + x * 3 + y * stride; if (pTemp[0] != 0) { int temp = pTemp[0]; pTemp[0] = (byte)((temp - mMin) * 255 / (mMax - mMin)); pTemp[1] = (byte)((temp - mMin) * 255 / (mMax - mMin)); pTemp[2] = (byte)((temp - mMin) * 255 / (mMax - mMin)); } } } src.UnlockBits(srcData); return src; } } /// <summary> /// Chamfer distance transform(using two 3*3 windows:forward window434 300 000;backward window 000 003 434). /// </summary> /// <param name="src">The source image.</param> /// <returns>The result image.</returns> public Bitmap ChamferDistancetransfrom(Bitmap src) { int w = src.Width; int h = src.Height; double p0 = 0, p1 = 0, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0, p7 = 0, p8 = 0, min = 0; int mMax = 0, mMin = 255; System.Drawing.Imaging.BitmapData srcData = src.LockBits(new Rectangle(0, 0, w, h), System.Drawing.Imaging.ImageLockMode.ReadOnly, System.Drawing.Imaging.PixelFormat.Format24bppRgb); unsafe { byte* p = (byte*)srcData.Scan0.ToPointer(); int stride = srcData.Stride; byte* pTemp; for (int y = 0; y < h; y++) { for (int x = 0; x < w; x++) { if (x > 0 && x < w - 1 && y > 0 && y < h - 1) { p0 = (p + x * 3 + y * stride)[0]; if (p0 != 0) { pTemp = p + (x - 1) * 3 + (y - 1) * stride; p2 = pTemp[0] + 4; pTemp = p + x * 3 + (y - 1) * stride; p3 = pTemp[0] + 3; pTemp = p + (x + 1) * 3 + (y - 1) * stride; p4 = pTemp[0] + 4; pTemp = p + (x - 1) * 3 + y * stride; p1 = pTemp[0] + 3; min = GetMin(p0, p1, p2, p3, p4); pTemp = p + x * 3 + y * stride; pTemp[0] = (byte)Math.Min(min, 255); pTemp[1] = (byte)Math.Min(min, 255); pTemp[2] = (byte)Math.Min(min, 255); } } else { pTemp = p + x * 3 + y * stride; pTemp[0] = 0; pTemp[1] = 0; pTemp[2] = 0; } } } for (int y = h - 1; y > 0; y--) { for (int x = w - 1; x > 0; x--) { if (x > 0 && x < w - 1 && y > 0 && y < h - 1) { p0 = (p + x * 3 + y * stride)[0]; if (p0 != 0) { pTemp = p + (x + 1) * 3 + y * stride; p5 = pTemp[0] + 3; pTemp = p + (x + 1) * 3 + (y + 1) * stride; p6 = pTemp[0] + 4; pTemp = p + x * 3 + (y + 1) * stride; p7 = pTemp[0] + 3; pTemp = p + (x - 1) * 3 + (y + 1) * stride; p8 = pTemp[0] + 4; min = GetMin(p0, p5, p6, p7, p8); pTemp = p + x * 3 + y * stride; pTemp[0] = (byte)Math.Min(min, 255); pTemp[1] = (byte)Math.Min(min, 255); pTemp[2] = (byte)Math.Min(min, 255); mMax = (int)Math.Max(min, mMax); } } else { pTemp = p + x * 3 + y * stride; pTemp[0] = 0; pTemp[1] = 0; pTemp[2] = 0; } } } mMin = 0; for (int y = 0; y < h; y++) { for (int x = 0; x < w; x++) { pTemp = p + x * 3 + y * stride; if (pTemp[0] != 0) { int temp = pTemp[0]; pTemp[0] = (byte)((temp - mMin) * 255 / (mMax - mMin)); pTemp[1] = (byte)((temp - mMin) * 255 / (mMax - mMin)); pTemp[2] = (byte)((temp - mMin) * 255 / (mMax - mMin)); } } } src.UnlockBits(srcData); return src; } } private double GetDistance(int x, int y, int dx, int dy, int mode) { double result = 0; switch (mode) { case 1: result = EuclideanDistance(x, y, dx, dy); break; case 2: result = CityblockDistance(x, y, dx, dy); break; case 3: result = ChessboardDistance(x, y, dx, dy); break; } return result; } private double EuclideanDistance(int x, int y, int dx, int dy) { return Math.Sqrt((x - dx) * (x - dx) + (y - dy) * (y - dy)); } private double CityblockDistance(int x, int y, int dx, int dy) { return Math.Abs(x - dx) + Math.Abs(y - dy); } private double ChessboardDistance(int x, int y, int dx, int dy) { return Math.Max(Math.Abs(x - dx), Math.Abs(y - dy)); } private double GetMin(double a, double b, double c, double d, double e) { return Math.Min(Math.Min(Math.Min(a, b), Math.Min(c, d)), e); } ~~~ 最后，分享一個專業的圖像處理網站（微像素），里面有很多源代碼下載： [http://www.zealpixel.com/portal.php](http://www.zealpixel.com/portal.php)