Skeleton extraction is a representation that reduces the binary image to 1 pixel wide. This is useful for feature extraction and / or representing the topology of objects.
# pip install scikit-image
from skimage.morphology import skeletonize
from skimage import data
import matplotlib.pyplot as plt
from skimage.util import invert
# Invert the horse image
image = invert(data.horse())
# perform skeletonization
skeleton = skeletonize(image)
# display results
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 4),
sharex=True, sharey=True)
ax = axes.ravel()
ax[0].imshow(image, cmap=plt.cm.gray)
ax[0].axis('off')
ax[0].set_title('original', fontsize=20)
ax[1].imshow(skeleton, cmap=plt.cm.gray)
ax[1].axis('off')
ax[1].set_title('skeleton', fontsize=20)
fig.tight_layout()
plt.show()
1. Zhang's method and Li's method
skeletonize works by continuously traversing the image and removing pixels on the object boundary. This continues until no more pixels can be removed. The image is associated with a mask that assigns a number in the range of [0... 255] to each pixel, corresponding to each possible mode of its 8 adjacent pixels. The lookup table is then used to assign values 0, 1, 2, or 3 to the pixels, which are selectively deleted during the iteration.
skeletonize(..., method='lee ') uses an octree data structure to check the 3x3x3 neighborhood of pixels. The algorithm scans the image through iteration and deletes pixels in each iteration until the image stops changing. Each iteration consists of two steps: first, assemble the candidate deletion list; Then recheck the order of pixels from this list to better maintain the connectivity of the image.
Note that Lee's method is intended for 3-D images and is automatically selected for these images. For illustration purposes, we apply this algorithm to two-dimensional images.
import matplotlib.pyplot as plt
from skimage.morphology import skeletonize
blobs = data.binary_blobs(200, blob_size_fraction=.2,
volume_fraction=.35, seed=1)
skeleton = skeletonize(blobs)
skeleton_lee = skeletonize(blobs, method='lee')
fig, axes = plt.subplots(1, 3, figsize=(8, 4), sharex=True, sharey=True)
ax = axes.ravel()
ax[0].imshow(blobs, cmap=plt.cm.gray)
ax[0].set_title('original')
ax[0].axis('off')
ax[1].imshow(skeleton, cmap=plt.cm.gray)
ax[1].set_title('skeletonize')
ax[1].axis('off')
ax[2].imshow(skeleton_lee, cmap=plt.cm.gray)
ax[2].set_title('skeletonize (Lee 94)')
ax[2].axis('off')
fig.tight_layout()
plt.show()
2. Skeletonization of central shaft
The central axis of an object is a collection of all points with more than one nearest point on the object boundary. It is often called a topological skeleton because it is the 1-pixel wide skeleton of the object and has the same connectivity as the original object. Here, we use the central axis transformation to calculate the width of the foreground object. Because the function media_ Axis returns the distance transformation (using the keyword parameter return_distance=True) in addition to the central axis, so you can use this function to calculate the distance from all points on the central axis to the background. This gives an estimate of the local width of the object. skeletonize should be preferred for skeletons with fewer branches.
from skimage.morphology import medial_axis, skeletonize
# Generate the data
blobs = data.binary_blobs(200, blob_size_fraction=.2,
volume_fraction=.35, seed=1)
# Compute the medial axis (skeleton) and the distance transform
skel, distance = medial_axis(blobs, return_distance=True)
# Compare with other skeletonization algorithms
skeleton = skeletonize(blobs)
skeleton_lee = skeletonize(blobs, method='lee')
# Distance to the background for pixels of the skeleton
dist_on_skel = distance * skel
fig, axes = plt.subplots(2, 2, figsize=(8, 8), sharex=True, sharey=True)
ax = axes.ravel()
ax[0].imshow(blobs, cmap=plt.cm.gray)
ax[0].set_title('original')
ax[0].axis('off')
ax[1].imshow(dist_on_skel, cmap='magma')
ax[1].contour(blobs, [0.5], colors='w')
ax[1].set_title('medial_axis')
ax[1].axis('off')
ax[2].imshow(skeleton, cmap=plt.cm.gray)
ax[2].set_title('skeletonize')
ax[2].axis('off')
ax[3].imshow(skeleton_lee, cmap=plt.cm.gray)
ax[3].set_title("skeletonize (Lee 94)")
ax[3].axis('off')
fig.tight_layout()
plt.show()
3. Morphological refinement
The morphological refinement implemented in the thin function works the same as skeletonization: pixels are removed from the boundary at each iteration until no pixels can be removed without changing the connectivity. Different removal rules can accelerate skeletonization and lead to different final skeletons.
The thin function also takes the optional max_num_iter keyword parameter to limit the number of refinement iterations, resulting in a relatively thick skeleton.
from skimage.morphology import skeletonize, thin
skeleton = skeletonize(image)
thinned = thin(image)
thinned_partial = thin(image, max_num_iter=25)
fig, axes = plt.subplots(2, 2, figsize=(8, 8), sharex=True, sharey=True)
ax = axes.ravel()
ax[0].imshow(image, cmap=plt.cm.gray)
ax[0].set_title('original')
ax[0].axis('off')
ax[1].imshow(skeleton, cmap=plt.cm.gray)
ax[1].set_title('skeleton')
ax[1].axis('off')
ax[2].imshow(thinned, cmap=plt.cm.gray)
ax[2].set_title('thinned')
ax[2].axis('off')
ax[3].imshow(thinned_partial, cmap=plt.cm.gray)
ax[3].set_title('partially thinned')
ax[3].axis('off')
fig.tight_layout()
plt.show()
4. The principle of zhang-suen image skeleton extraction algorithm and OPENCV implementation
Each iteration step of the algorithm is to corrode the target pixels that meet the specific conditions, and the effect is that the target becomes thinner and thinner. Continuously iterate until the target after the last corrosion is corroded. In this round of operation, no new pixels are corroded, and the algorithm ends.
#include<opencv2/opencv.hpp>
#include<iostream>
using namespace std;
using namespace cv;
void thinImage(Mat & srcImg) {
vector<Point> deleteList;
int neighbourhood[9];
int nl = srcImg.rows;
int nc = srcImg.cols;
bool inOddIterations = true;
while (true) {
for (int j = 1; j < (nl - 1); j++) {
uchar* data_last = srcImg.ptr<uchar>(j - 1);
uchar* data = srcImg.ptr<uchar>(j);
uchar* data_next = srcImg.ptr<uchar>(j + 1);
for (int i = 1; i < (nc - 1); i++) {
if (data[i] == 255) {
int whitePointCount = 0;
neighbourhood[0] = 1;
if (data_last[i] == 255) neighbourhood[1] = 1;
else neighbourhood[1] = 0;
if (data_last[i + 1] == 255) neighbourhood[2] = 1;
else neighbourhood[2] = 0;
if (data[i + 1] == 255) neighbourhood[3] = 1;
else neighbourhood[3] = 0;
if (data_next[i + 1] == 255) neighbourhood[4] = 1;
else neighbourhood[4] = 0;
if (data_next[i] == 255) neighbourhood[5] = 1;
else neighbourhood[5] = 0;
if (data_next[i - 1] == 255) neighbourhood[6] = 1;
else neighbourhood[6] = 0;
if (data[i - 1] == 255) neighbourhood[7] = 1;
else neighbourhood[7] = 0;
if (data_last[i - 1] == 255) neighbourhood[8] = 1;
else neighbourhood[8] = 0;
for (int k = 1; k < 9; k++) {
whitePointCount += neighbourhood[k];
}
if ((whitePointCount >= 2) && (whitePointCount <= 6)) {
int ap = 0;
if ((neighbourhood[1] == 0) && (neighbourhood[2] == 1)) ap++;
if ((neighbourhood[2] == 0) && (neighbourhood[3] == 1)) ap++;
if ((neighbourhood[3] == 0) && (neighbourhood[4] == 1)) ap++;
if ((neighbourhood[4] == 0) && (neighbourhood[5] == 1)) ap++;
if ((neighbourhood[5] == 0) && (neighbourhood[6] == 1)) ap++;
if ((neighbourhood[6] == 0) && (neighbourhood[7] == 1)) ap++;
if ((neighbourhood[7] == 0) && (neighbourhood[8] == 1)) ap++;
if ((neighbourhood[8] == 0) && (neighbourhood[1] == 1)) ap++;
if (ap == 1) {
if (inOddIterations && (neighbourhood[3] * neighbourhood[5] * neighbourhood[7] == 0)
&& (neighbourhood[1] * neighbourhood[3] * neighbourhood[5] == 0)) {
deleteList.push_back(Point(i, j));
}
else if (!inOddIterations && (neighbourhood[1] * neighbourhood[5] * neighbourhood[7] == 0)
&& (neighbourhood[1] * neighbourhood[3] * neighbourhood[7] == 0)) {
deleteList.push_back(Point(i, j));
}
}
}
}
}
}
if (deleteList.size() == 0)
break;
for (size_t i = 0; i < deleteList.size(); i++) {
Point tem;
tem = deleteList[i];
uchar* data = srcImg.ptr<uchar>(tem.y);
data[tem.x] = 0;
}
deleteList.clear();
inOddIterations = !inOddIterations;
}
}
int main()
{
string srcImg = "opencv.png";
Mat img = imread(srcImg);
cvtColor(img, img, COLOR_BGR2GRAY);
threshold(img, img, 127, 255, THRESH_OTSU);
namedWindow("src", 0);
imshow("src", img);
thinImage(img);
namedWindow("dst", 0);
imshow("dst", img);
waitKey(0);
}
Original drawing:
After refinement:
5. Skeleton extraction based on ITK
Template: template < typename tinputimage, typename toutputimage > class binarythinningimagefilter: public ITK:: imagetoimagefilter < tinputimage, toutputimage > this function calculates a pixel wide edge of the input image. This function is parameterized by the type of input image and the type of output image.
Suppose the input is a binary image. If the foreground pixel value of the input image is not 1, they are internally rescaled to 1 to simplify the calculation.
Function generates the skeleton of the object. The output background value is 0 and the foreground value is 1.
This function is a continuous thinning algorithm, and the calculation time depends on the size of the image.
- python implementation
#!/usr/bin/env python
import itk
import argparse
parser = argparse.ArgumentParser(description="Thin Image.")
parser.add_argument("input_image", nargs="?")
args = parser.parse_args()
PixelType = itk.UC
Dimension = 2
ImageType = itk.Image[PixelType, Dimension]
if args.input_image:
image = itk.imread(args.input_image)
else:
# Create an image
start = itk.Index[Dimension]()
start.Fill(0)
size = itk.Size[Dimension]()
size.Fill(100)
region = itk.ImageRegion[Dimension]()
region.SetIndex(start)
region.SetSize(size)
image = ImageType.New(Regions=region)
image.Allocate()
image.FillBuffer(0)
# Draw a 5 pixel wide line
image[50:55, 20:80] = 255
# Write Image
itk.imwrite(image, "input.png")
image = itk.binary_thinning_image_filter(image)
# Rescale the image so that it can be seen (the output is 0 and 1, we want 0 and 255)
image = itk.rescale_intensity_image_filter(image, output_minimum=0, output_maximum=255)
# Write Image
itk.imwrite(image, "outputPython.png")
- C + + implementation
#include"itkBinaryThinningImageFilter.h"
#include"itkImage.h"
#include"itkImageFileReader.h"
#include"itkImageFileWriter.h"
#include"itkRescaleIntensityImageFilter.h"
// PNG correspondence
#include <itkPNGImageIOFactory.h>
// BMP correspondence
#include <itkBMPImageIOFactory.h>
// JPG correspondence
#include <itkJPEGImageIOFactory.h>
// mhd format image
// #include <itkMetaImageIOFactory.h>
using ImageType=itk::Image<unsigned char,2>;
static void WriteImage(const ImageType::Pointer image,const std::string&fileName);
static void CreateImage(ImageType::Pointer image);
int main(int argc, char*argv[])
{
ImageType::Pointer image = ImageType::New();
if (argc == 2) {
itk::ObjectFactoryBase::RegisterFactory(itk::PNGImageIOFactory::New());
// Or itk::BMPImageIOFactory::New() / / or itk::JPGImageIOFactory::New()
using ImageReader = itk::ImageFileReader<ImageType>;
ImageReader::Pointer reader = ImageReader::New();
std::string fileName = argv[1];
reader->SetFileName(fileName);
reader->Update();
image = reader->GetOutput();
}
else {
itk::ObjectFactoryBase::RegisterFactory(itk::PNGImageIOFactory::New());
// Or itk::BMPImageIOFactory::New() / / or itk::JPGImageIOFactory::New()
CreateImage(image);
WriteImage(image, "input.png");
}
using BinaryThinningImageFilterType = itk::BinaryThinningImageFilter<ImageType, ImageType>;
BinaryThinningImageFilterType::Pointer binaryThinningImageFilter = BinaryThinningImageFilterType::New();
binaryThinningImageFilter->SetInput(image);
binaryThinningImageFilter->Update();
//Rescale the image so that it can be seen(the output is 0 and 1,we want 0 and 255)
using RescaleType = itk::RescaleIntensityImageFilter<ImageType, ImageType>;
RescaleType::Pointer rescaler = RescaleType::New();
rescaler->SetInput(binaryThinningImageFilter->GetOutput());
rescaler->SetOutputMinimum(0);
rescaler->SetOutputMaximum(255);
rescaler->Update();
WriteImage(rescaler->GetOutput(), "output.png");
return EXIT_SUCCESS;
}
void CreateImage(ImageType::Pointer image)
{//Create an image
ImageType::IndexType start;
start.Fill(0);
ImageType::SizeType size;
size.Fill(100);
ImageType::RegionType region(start,size);
image->SetRegions(region);
image->Allocate();
image->FillBuffer(0);
//Draw a 5pixel wide line
for(unsigned int i=20;i<80;++i){
for(unsigned int j=50;j<55;++j){
itk::Index<2>index;
index[0]=i;
index[1]=j;
image->SetPixel(index,255);
}
}
}
void WriteImage(const ImageType::Pointer image, const std::string&fileName) {
using WriterType = itk::ImageFileWriter<ImageType>;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(fileName);
writer->SetInput(image);
writer->Update();
}
Reference catalogue
https://www.freesion.com/article/8054209684/
https://scikit-image.org/docs/dev/auto_examples/edges/plot_skeleton.html