提取图像的骨架(Skeleton)算法

图像的骨架似乎没有严格的数学定义，可认为是图像细化（Thinning）的产物（中轴可以看作一种骨架，其有严格的数学定义）。目前已经有许多细化算法，这些算法得到的骨架可能略有差异。本文实现了Khalid Sheed 的 K3M算法。该算法属于迭代腐蚀边界的一类算法，该类算法的思想是，假定从二值图像中物体的边界处同时开始燃烧，物体就会被逐步细化，但在燃烧过程中要保证满足一定条件的点被保留或者被“烧掉”，以确定燃烧结束后，剩下的最后一像素宽的图像为图像的骨架。这些条件的确定没有统一的标准，各个算法采取了不同的方案。一般来讲，为了满足计算的速度要求和算法的准确，迭代中算法会对图像边界上某点的3*3邻域内进行检查，判断是否满足要求。

K3M算法在每次迭代中需要进行六次检查

Phase 0. 标记出图像的边界，

Phase 1. 如果该点的邻域中有3个点（非零，以下皆如此）相邻，删除该点

Phase 2. 如果该点的邻域中有3或4个点相邻，删除该点。

Phase 3. 如果该点的邻域中有3，4，5个点相邻，删除该点。

Phase 4. 如果该点的邻域中有3，4，5，6个点相邻，删除该点。

Phase 5. 如果该点的邻域中有3，4，5，6，7个点相邻，删除该点。

Phase 6. 剩余的边界点取消标记，如果Phase 5中没有点被修改，停止迭代，否则返回Phase 0。

具体的步骤可以阅读论文：http://matwbn.icm.edu.pl/ksiazki/amc/amc20/amc2029.pdf。论文中算法实现的一个小技巧是，对邻域中的8个点的值看作二进制，即二进制编码，这样不同的值就对应邻域中不同的状态。迭代中通过计算即可判断该点是否满足条件，是否可以删除。具体细节请移步阅读论文。

算法的测试结果如下：

参考：

http://homepages.inf.ed.ac.uk/rbf/HIPR2/skeleton.htm

http://home.agh.edu.pl/~saeed/arts/2001%20CAIP.pdf

http://www.cs.sunysb.edu/~algorith/files/thinning.shtml

http://matwbn.icm.edu.pl/ksiazki/amc/amc20/amc2029.pdf

#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/objdetect/objdetect.hpp>
#include <opencv2/highgui/highgui.hpp>
#include<vector>
#include<iostream>
#include<algorithm>

using std::vector;

vector<int> GetFlags(int a[],int length)
{
	vector<int> vec;
	int neighbour[]={1,2,4,8,16,32,64,128,1,2,4,8,16,32,64};
	for(int i=0;i<length;i++)
	{
		for(int j=0;j<8;j++)
		{
			int sum=0;
			for(int k=j;k<=j+a[i];k++)
				sum+=neighbour[k];
			vec.push_back(sum);
			std::cout<<sum<<" ";
		}
	}
	std::cout<<std::endl;
	return vec;
}
void skeleton(cv::Mat &Input) //Input-binary image
{
	int a0[]={1,2,3,4,5,6};
	int a1[]={2};
	int a2[]={2,3};
	int a3[]={2,3,4};
	int a4[]={2,3,4,5};
	int a5[]={2,3,4,5,6};
	vector<int> A0=GetFlags(a0,6);

	vector<int> A1=GetFlags(a1,1);

	vector<int> A2=GetFlags(a2,2);
	vector<int> A3=GetFlags(a3,3);
	vector<int> A4=GetFlags(a4,4);
	vector<int> A5=GetFlags(a5,5);
	vector<cv::Point2i> border;
	bool modify=true;
	int neighbour[3][3]={
		{128,1,2},
		{64,0,4},
		{32,16,8}
	};
	int row=Input.rows;
	int col=Input.cols;
	while(modify)
	{
		modify=false;
		// flag the border Pharse 0
		for(int m=1;m<row-1;++m)
		{
			for(int n=1;n<col-1;++n)
			{
				int weight=0;
				for(int j=-1;j<=1;++j)
				{
					for(int k=-1;k<=1;k++)
					{
						weight+=neighbour[j+1][k+1]*Input.at<uchar>(m+j,n+k);
					}
				}
				if(std::find(A0.begin(),A0.end(),weight)!=A0.end())
					border.push_back(cv::Point2i(m,n));
			}
		}
		//Pharse 1
		vector<cv::Point2i>::iterator first=border.begin();
		while(first!=border.end())
		{
			int weight=0;
			for(int j=-1;j<=1;++j)
			{
				for(int k=-1;k<=1;k++)
				{
					weight+=neighbour[j+1][k+1]*Input.at<uchar>((*first).x+j,(*first).y+k);
				}
			}
			if(std::find(A1.begin(),A1.end(),weight)!=A1.end())
			{
				Input.at<uchar>((*first).x,(*first).y)=0;
				first=border.erase(first);
			}
			else
			++first;
		}
		//Pharse2
		first=border.begin();
		while(first!=border.end())
		{
			int weight=0;
			for(int j=-1;j<=1;++j)
			{
				for(int k=-1;k<=1;k++)
				{
					weight+=neighbour[j+1][k+1]*Input.at<uchar>((*first).x+j,(*first).y+k);
				}
			}
			if(std::find(A2.begin(),A2.end(),weight)!=A2.end())
			{
				Input.at<uchar>((*first).x,(*first).y)=0;
				first=border.erase(first);
			}
			else
			++first;
		}
		//Pharse3
		first=border.begin();
		while(first!=border.end())
		{
			int weight=0;
			for(int j=-1;j<=1;++j)
			{
				for(int k=-1;k<=1;k++)
				{
					weight+=neighbour[j+1][k+1]*Input.at<uchar>((*first).x+j,(*first).y+k);
				}
			}
			if(std::find(A3.begin(),A3.end(),weight)!=A3.end())
			{
				Input.at<uchar>((*first).x,(*first).y)=0;
				first=border.erase(first);
			}
			else
			++first;
		}
		//Pharse4
		first=border.begin();
		while(first!=border.end())
		{
			int weight=0;
			for(int j=-1;j<=1;++j)
			{
				for(int k=-1;k<=1;k++)
				{
					weight+=neighbour[j+1][k+1]*Input.at<uchar>((*first).x+j,(*first).y+k);
				}
			}
			if(std::find(A4.begin(),A4.end(),weight)!=A4.end())
			{
				Input.at<uchar>((*first).x,(*first).y)=0;
				first=border.erase(first);
			}
			else
			++first;
		}
		//Pharse5
		first=border.begin();
		while(first!=border.end())
		{
			int weight=0;
			for(int j=-1;j<=1;++j)
			{
				for(int k=-1;k<=1;k++)
				{
					weight+=neighbour[j+1][k+1]*Input.at<uchar>((*first).x+j,(*first).y+k);
				}
			}
			if(std::find(A5.begin(),A5.end(),weight)!=A5.end())
			{
				Input.at<uchar>((*first).x,(*first).y)=0;
				first=border.erase(first);
				modify=true;
			}
			else
			++first;
		}
		//Pharse6
		border.clear();
	}
	for(int m=1;m<row-1;++m)
	{
		for(int n=1;n<col-1;++n)
		{
			int weight=0;
			for(int j=-1;j<=1;++j)
			{
				for(int k=-1;k<=1;k++)
				{
					weight+=neighbour[j+1][k+1]*Input.at<uchar>(m+j,n+k);
				}
			}
			if(std::find(A0.begin(),A0.end(),weight)!=A0.end())
					Input.at<uchar>(m,n)=0;;
		}
	}

}
int main()
{
	cv::Mat raw=cv::imread("test.bmp");
	cv::Mat image(raw.rows,raw.cols,CV_8UC1);
	cv::cvtColor(raw,image,CV_RGB2GRAY);
	cv::Mat binaryImage(image.rows,image.cols,CV_8UC1);
	cv::threshold(image,binaryImage,150,1,CV_THRESH_BINARY_INV);

	cv::imshow("input",image);
	skeleton(binaryImage);
		for(int p=0;p<binaryImage.rows;p++)
		{
			for(int q=0;q<binaryImage.cols;q++)
			{
				if(binaryImage.at<uchar>(p,q)==1)
					binaryImage.at<uchar>(p,q)=0;
				else
					binaryImage.at<uchar>(p,q)=255;
			}
		}
	cv::imshow("output",binaryImage);

	cv::waitKey(0);
	return 0;
}