《图像处理实例》 之 中轴线提取

 

 

---

 

中轴线算法

这是http://www.imagepy.org/的作者原创，我只是对其理解之后改进和说明，欢迎大家使用这个小软件！

如果有需要C++版本的朋友，可以博文结尾留邮箱！

 

---

 

首先上效果图：

 

 

算法的流程：

　　　　

　　　　　　第一步：

　　　　　　　　　　距离变换

　　　　　　第二步：

　　　　　　　　　　把距离变换的图像进行像素值的排列（排列返回像素的位置信息）

　　　　　　第二步：

　　　　　　　　　　从小到大进行像素的查表操作

 

　　　　　　注释：这里为什么叫中轴线提取？因为提取的过程是绝对的按照对称来的，距离变换的结果就是前景到背景的距离，所以结果是绝对的集合中心。

 

 

 

 

Python代码：

 

import numpy as np
from skimage.data import horse, camera
import matplotlib.pyplot as plt
import scipy.ndimage as ndimg
from numba import jit
import cv2
from scipy.ndimage import label, generate_binary_structure

strc = np.ones((3, 3), dtype=np.bool)


# check whether this pixcel can be removed
def check(n):
    a = [(n >> i) & 1 for i in range(8)]
    a.insert(4, 0)  # make the 3x3 unit
    # if up, down, left, right all are 1, you cannot make a hole
    # if a[1] & a[3] & a[5] & a[7]:return False
    a = np.array(a).reshape((3, 3))
    # segments
    n = label(a, strc)[1]
    # if sum is 0, it is a isolate point, you cannot remove it.
    # if number of segments > 2, you cannot split them.
    return n < 2
    return a.sum() > 1 and n < 2
    if a.sum() == 1 or n > 2: return 2
    if a.sum() > 1 and n < 2: return 1
    return 0


lut = np.array([check(n) for n in range(256)])
lut = np.dot(lut.reshape((-1, 8)), [1, 2, 4, 8, 16, 32, 64, 128]).astype(np.uint8)
'''
lut = np.array([200, 206, 220, 204, 0, 207, 0, 204, 0, 207, 221, 51, 1, 207, 221, 51,
       0, 0, 221, 204, 0, 0, 0, 204, 1, 207, 221, 51, 1, 207, 221, 51], dtype=np.int8)
'''


@jit
def skel2dp(data, idx, lup):
    h, w = data.shape
    data = data.ravel()     
    for id in idx:          

        if data[id] == 0: continue  
        i2 = id - w
        i8 = id + w
        i1 = i2 - 1
        i3 = i2 + 1
        i4 = id - 1
        i6 = id + 1
        i7 = i8 - 1
        i9 = i8 + 1
        c = (data[i1] > 0) << 0 | (data[i2] > 0) << 1 \
            | (data[i3] > 0) << 2 | (data[i4] > 0) << 3 \
            | (data[i6] > 0) << 4 | (data[i7] > 0) << 5 \
            | (data[i8] > 0) << 6 | (data[i9] > 0) << 7
        if (lup[c // 8] >> c % 8) & 1: data[id] = 0
    return 0


def mid_axis(img):
    dis = ndimg.distance_transform_edt(img)
    idx = np.argsort(dis.flat).astype(np.int32)     
    skel2dp(dis, idx, lut)
    return dis


from time import time
img = ~horse()*255
#img = cv2.imread('123.jpg')
#img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
#ret2, img = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
dis = ndimg.distance_transform_edt(img)
plt.imshow(dis)
idx = np.argsort(dis.flat).astype(np.int32)
a = skel2dp(dis, idx, lut)
#mid_axis(img.copy())
t1 = time()
a = mid_axis(img)
t2 = time()
print(t2 - t1)
plt.imshow(a)
plt.show()

 

 

 

C++代码：

  

void center_axis(InputArray _src, Mat&  dst)
{
    typedef struct MyStruct
    {
        Point position;
        float num;
    }MyStruct;
    int wjy_array[] = { 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1,
        1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1,
        0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1,
        1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1,
        1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1,
        1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1,
        0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1,
        1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
        1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
        1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0,
        1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0 };
    uchar lut[] = { 200, 206, 220, 204, 0, 207, 0, 204, 0, 207, 221, 51, 1, 207, 221, 51,
        0, 0, 221, 204, 0, 0, 0, 204, 1, 207, 221, 51, 1, 207, 221, 51 };
    Mat src = _src.getMat();
    //Mat dst = _dst.getMat();
    distanceTransform(src, src, DIST_L2, DIST_MASK_3, 5);
    normalize(src, src, 0, 255, NORM_MINMAX);
    Mat img_row = src.reshape(0, 1);
    vector<MyStruct> my_vector;
    for (size_t j = 0; j < img_row.cols; j++)
    {
        if (img_row.at<float>(0, j) == 0) continue;
        MyStruct my_struct;
        my_struct.num = saturate_cast<float>(img_row.at<float>(0, j));
        my_struct.position = Point(saturate_cast<int>(j % src.cols), saturate_cast<int>(j / src.cols));
        my_vector.push_back(my_struct);
    }
    for (size_t i = 0; i < my_vector.size(); i++)
    {
        if (my_vector[i].num == 0) continue;
        for (size_t j = i; j < my_vector.size(); j++)
        {
            MyStruct temp;
            if (my_vector[i].num >= my_vector[j].num)
            {
                if (my_vector[j].num == 0) continue;
                temp = my_vector[j];
                my_vector[j] = my_vector[i];
                my_vector[i] = temp;
            }
        }
    }
    for (size_t i = 0; i < my_vector.size(); i++)
    {
        if (my_vector[i].position.y == 1
            || my_vector[i].position.x == 1
            || my_vector[i].position.y == src.rows - 1
            || my_vector[i].position.x == src.cols - 1
            || src.at<float>(my_vector[i].position.y, my_vector[i].position.x) == 0) continue;
        else
        {
            char num[] = { 1,1,1,1,1,1,1,1 };
            num[0] = src.at<float>(my_vector[i].position.y - 1, my_vector[i].position.x - 1)
    > 0 ? 1 : 0;
            num[1] = src.at<float>(my_vector[i].position.y - 1, my_vector[i].position.x)
        > 0 ? 1 : 0;
            num[2] = src.at<float>(my_vector[i].position.y - 1, my_vector[i].position.x + 1)
    > 0 ? 1 : 0;
            num[3] = src.at<float>(my_vector[i].position.y, my_vector[i].position.x - 1)
    > 0 ? 1 : 0;
            num[4] = src.at<float>(my_vector[i].position.y, my_vector[i].position.x + 1)
                        > 0 ? 1 : 0;
            num[5] = src.at<float>(my_vector[i].position.y + 1, my_vector[i].position.x - 1)
                        > 0 ? 1 : 0;
            num[6] = src.at<float>(my_vector[i].position.y + 1, my_vector[i].position.x)
                        > 0 ? 1 : 0;
            num[7] = src.at<float>(my_vector[i].position.y + 1, my_vector[i].position.x + 1)
                        > 0 ? 1 : 0;
            int sum = num[0] + num[1] * 2 + num[2] * 4 + num[3] * 8
                + num[4] * 16 + num[5] * 32 + num[6] * 64 + num[7] * 128;
            src.at<float>(my_vector[i].position.y, my_vector[i].position.x) = ((lut[uchar(sum / 8)] >> sum % 8) & 1) != 1 ? 255 : 0;
        }
    }
    dst = src.clone();
    dst.convertTo(dst, CV_8UC1);
}