【opencv】經典的細化提取骨架理論及原始碼
本章我們學習一下Hilditch演算法的基本原理，從網上找資料的時候，竟然發現兩個有很大差別的演算法描述，而且都叫Hilditch演算法。不知道那一個才是正宗的，兩個演算法實現的效果接近，第一種演算法更好一些。
第一種演算法描述參考paper和程式碼：
Linear Skeletons from Square Cupboards
Speedup Method for Real-Time Thinning Algorithm
第二種演算法描述參考資料：
下面我們分別看一下這兩種演算法描述：
一、第一種演算法描述
假設當前被處理的畫素為p0，我們使用下圖所示的8鄰域表示方式。
我們處理的為二值影象，背景為黑色，值為0，要細化的前景物體畫素值為255。
對於Hilditch演算法來說，它並不是一個完全的並行演算法，而是序列並行相結合。當前畫素是否是能夠刪除的骨架點，不僅是由它周圍的8鄰域決定，而且和前面畫素的判定結果有關。一個畫素判定為可以刪除，我們並不直接刪除它，而是在目地影象中設定畫素值為GRAY=128,這個資訊可能會影響之後其它畫素的判定。
當影象一次掃描迭代完成後，我們把所有置為GRAY的畫素設定為0，從而刪除它。
演算法的描述如下。
迭代掃描當前影象
對於當前畫素點，掃描它的8鄰域，如果鄰域的畫素值為255，則b[i]=1(i=0…8)，畫素值為128(GRAY，表示該畫素點在前面的迴圈中被標記為刪除
下面會根據b[i]的值進行6個條件判斷，如果條件滿足，則會標記該畫素值為GRAY(128)。
1. b[0]=1，即當前畫素必須為前景點。
2. 1-abs(b1) + 1 – abs(b3) + 1 – abs(b5) + 1 – abs(b7) >= 1，該條件表示當前畫素為邊界點，即東西南北四個點至少有一個b[i]=0。
3. abs(b1)+…+abs(b8)>=2， 該條件表示不能刪除端點，即p0點周圍只有一個點為1或-1的情況。
4. 統計b1到b8等於1的數量，該數量值必須大於1，該條件表示不能刪除端點。
5. 連通性檢測，使用下面的公式：首先根據當前畫素周圍3*3域的值，記錄d[9]陣列，如果b[i]等於0，則d[i]=0， 否則d[i]=1，最後計算 d1-d1*d2*d3+d3-d3*d4*d5+d5-d5*d6*d7+d7-d7*d8*d1是否為1，為1則滿足連通性，可以刪除。
6.最後一個條件保證當輪廓是2個畫素寬時，只刪除一邊。統計sum的值，當值為8時候，可以刪除。
sum = 0; for (i = 1; i <= 8; i++) { if (b[i] != -1) { sum++; } else { copy = b[i]; b[i] = 0; if (func_nc8(b) == 1) sum++; b[i] = copy; }
當這6個條件都滿足時候，標記當前畫素值為GRAY(128)，然後在判斷別的畫素。當所有畫素都掃描一遍後，完成一次迭代。
此時我們會把剛才標記為GARY的畫素，都設定為0，真正的刪除它，如果上一次迴圈已經沒有標記刪除的畫素，則退出迭代，否則進行下一次迭代。
演算法程式碼：
int gThin::func_nc8(int *b)
//端點的連通性檢測
{
int n_odd[4] = { 1, 3, 5, 7 }; //四鄰域
int i, j, sum, d[10];
for (i = 0; i <= 9; i++) {
j = i;
if (i == 9) j = 1;
if (abs(*(b + j)) == 1)
{
d[i] = 1;
}
else
{
d[i] = 0;
}
}
sum = 0;
for (i = 0; i < 4; i++)
{
j = n_odd[i];
sum = sum + d[j] - d[j] * d[j + 1] * d[j + 2];
}
return (sum);
}
void gThin::cvHilditchThin(cv::Mat& src, cv::Mat& dst)
{
if(src.type()!=CV_8UC1)
{
printf("只能處理二值或灰度影象\n");
return;
}
//非原地操作時候，copy src到dst
if(dst.data!=src.data)
{
src.copyTo(dst);
}
//8鄰域的偏移量
int offset[9][2] = {{0,0},{1,0},{1,-1},{0,-1},{-1,-1},
{-1,0},{-1,1},{0,1},{1,1} };
//四鄰域的偏移量
int n_odd[4] = { 1, 3, 5, 7 };
int px, py;
int b[9]; //3*3格子的灰度資訊
int condition[6]; //1-6個條件是否滿足
int counter; //移去畫素的數量
int i, x, y, copy, sum;
uchar* img;
int width, height;
width = dst.cols;
height = dst.rows;
img = dst.data;
int step = dst.step ;
do
{
counter = 0;
for (y = 0; y < height; y++)
{
for (x = 0; x < width; x++)
{
//前面標記為刪除的畫素，我們置其相應鄰域值為-1
for (i = 0; i < 9; i++)
{
b[i] = 0;
px = x + offset[i][0];
py = y + offset[i][1];
if (px >= 0 && px < width && py >= 0 && py <height)
{
// printf("%d\n", img[py*step+px]);
if (img[py*step+px] == WHITE)
{
b[i] = 1;
}
else if (img[py*step+px] == GRAY)
{
b[i] = -1;
}
}
}
for (i = 0; i < 6; i++)
{
condition[i] = 0;
}
//條件1，是前景點
if (b[0] == 1) condition[0] = 1;
//條件2，是邊界點
sum = 0;
for (i = 0; i < 4; i++)
{
sum = sum + 1 - abs(b[n_odd[i]]);
}
if (sum >= 1) condition[1] = 1;
//條件3， 端點不能刪除
sum = 0;
for (i = 1; i <= 8; i++)
{
sum = sum + abs(b[i]);
}
if (sum >= 2) condition[2] = 1;
//條件4， 孤立點不能刪除
sum = 0;
for (i = 1; i <= 8; i++)
{
if (b[i] == 1) sum++;
}
if (sum >= 1) condition[3] = 1;
//條件5， 連通性檢測
if (func_nc8(b) == 1) condition[4] = 1;
//條件6，寬度為2的骨架只能刪除1邊
sum = 0;
for (i = 1; i <= 8; i++)
{
if (b[i] != -1)
{
sum++;
} else
{
copy = b[i];
b[i] = 0;
if (func_nc8(b) == 1) sum++;
b[i] = copy;
}
}
if (sum == 8) condition[5] = 1;
if (condition[0] && condition[1] && condition[2] &&condition[3] && condition[4] && condition[5])
{
img[y*step+x] = GRAY; //可以刪除，置位GRAY，GRAY是刪除標記，但該資訊對後面畫素的判斷有用
counter++;
//printf("----------------------------------------------\n");
//PrintMat(dst);
}
}
}
if (counter != 0)
{
for (y = 0; y < height; y++)
{
for (x = 0; x < width; x++)
{
if (img[y*step+x] == GRAY)
img[y*step+x] = BLACK;
}
}
}
}while (counter != 0);
}二、第二種演算法描述
第二種演算法描述和Zhang的並行演算法很相似，特別是前2個條件一模一樣，不同的是3，4兩個條件，還有就是該描述演算法並沒有像zhang演算法那樣，把一次迭代分成2個階段。
此時我們使用的8鄰域標記為：
下面看下它的演算法描述：
複製目地影象到臨時影象，對臨時影象進行一次掃描，對於不為0的點，如果滿足以下四個條件，則在目地影象中刪除該點(就是設定該畫素為0)
a. 2<= p2+p3+p4+p5+p6+p7+p8+p9<=6
大於等於2會保證p1點不是端點或孤立點，因為刪除端點和孤立點是不合理的，小於等於6保證p1點是一個邊界點，而不是一個內部點。等於0時候，周圍沒有等於1的畫素，所以p1為孤立點，等於1的時候，周圍只有1個灰度等於1的畫素，所以是端點（注：端點是周圍有且只能有1個值為1的畫素)。
b. p2->p9的排列順序中，01模式的數量為1，比如下面的圖中，有p2p3 => 01, p6p7=>01,所以該畫素01模式的數量為2。
之所以要01模式數量為1，是要保證刪除當前畫素點後的連通性。比如下面的圖中，01模式數量大於1，如果刪除當前點p1，則連通性不能保證。
c. p2.p4.p8 = 0 or A(p2)!=1，A(p2)表示p2周圍8鄰域的01模式和。這個條件保證2個畫素寬的垂直條不完全被腐蝕掉。
d.p2.p4.p6 = 0 or A(p4)!=1，A(p4)表示p4周圍8鄰域的01模式和。這個條件保證2個畫素寬的水平條不完全被腐蝕掉。
演算法程式碼：
void gThin::cvHilditchThin1(cv::Mat& src, cv::Mat& dst)
{
//http://cgm.cs.mcgill.ca/~godfried/teaching/projects97/azar/skeleton.html#algorithm
//演算法有問題，得不到想要的效果
if(src.type()!=CV_8UC1)
{
printf("只能處理二值或灰度影象\n");
return;
}
//非原地操作時候，copy src到dst
if(dst.data!=src.data)
{
src.copyTo(dst);
}
int i, j;
int width, height;
//之所以減2，是方便處理8鄰域，防止越界
width = src.cols -2;
height = src.rows -2;
int step = src.step;
int p2,p3,p4,p5,p6,p7,p8,p9;
uchar* img;
bool ifEnd;
int A1;
cv::Mat tmpimg;
while(1)
{
dst.copyTo(tmpimg);
ifEnd = false;
img = tmpimg.data+step;
for(i = 2; i < height; i++)
{
img += step;
for(j =2; j<width; j++)
{
uchar* p = img + j;
A1 = 0;
if( p[0] > 0)
{
if(p[-step]==0&&p[-step+1]>0) //p2,p3 01模式
{
A1++;
}
if(p[-step+1]==0&&p[1]>0) //p3,p4 01模式
{
A1++;
}
if(p[1]==0&&p[step+1]>0) //p4,p5 01模式
{
A1++;
}
if(p[step+1]==0&&p[step]>0) //p5,p6 01模式
{
A1++;
}
if(p[step]==0&&p[step-1]>0) //p6,p7 01模式
{
A1++;
}
if(p[step-1]==0&&p[-1]>0) //p7,p8 01模式
{
A1++;
}
if(p[-1]==0&&p[-step-1]>0) //p8,p9 01模式
{
A1++;
}
if(p[-step-1]==0&&p[-step]>0) //p9,p2 01模式
{
A1++;
}
p2 = p[-step]>0?1:0;
p3 = p[-step+1]>0?1:0;
p4 = p[1]>0?1:0;
p5 = p[step+1]>0?1:0;
p6 = p[step]>0?1:0;
p7 = p[step-1]>0?1:0;
p8 = p[-1]>0?1:0;
p9 = p[-step-1]>0?1:0;
//計算AP2,AP4
int A2, A4;
A2 = 0;
//if(p[-step]>0)
{
if(p[-2*step]==0&&p[-2*step+1]>0) A2++;
if(p[-2*step+1]==0&&p[-step+1]>0) A2++;
if(p[-step+1]==0&&p[1]>0) A2++;
if(p[1]==0&&p[0]>0) A2++;
if(p[0]==0&&p[-1]>0) A2++;
if(p[-1]==0&&p[-step-1]>0) A2++;
if(p[-step-1]==0&&p[-2*step-1]>0) A2++;
if(p[-2*step-1]==0&&p[-2*step]>0) A2++;
}
A4 = 0;
//if(p[1]>0)
{
if(p[-step+1]==0&&p[-step+2]>0) A4++;
if(p[-step+2]==0&&p[2]>0) A4++;
if(p[2]==0&&p[step+2]>0) A4++;
if(p[step+2]==0&&p[step+1]>0) A4++;
if(p[step+1]==0&&p[step]>0) A4++;
if(p[step]==0&&p[0]>0) A4++;
if(p[0]==0&&p[-step]>0) A4++;
if(p[-step]==0&&p[-step+1]>0) A4++;
}
//printf("p2=%d p3=%d p4=%d p5=%d p6=%d p7=%d p8=%d p9=%d\n", p2, p3, p4, p5, p6,p7, p8, p9);
//printf("A1=%d A2=%d A4=%d\n", A1, A2, A4);
if((p2+p3+p4+p5+p6+p7+p8+p9)>1 && (p2+p3+p4+p5+p6+p7+p8+p9)<7 && A1==1)
{
if(((p2==0||p4==0||p8==0)||A2!=1)&&((p2==0||p4==0||p6==0)||A4!=1))
{
dst.at<uchar>(i,j) = 0; //滿足刪除條件，設定當前畫素為0
ifEnd = true;
//printf("\n");
//PrintMat(dst);
}
}
}
}
}
//printf("\n");
//PrintMat(dst);
//PrintMat(dst);
//已經沒有可以細化的畫素了，則退出迭代
if(!ifEnd) break;
}
}