LAB EXERCISE

什么是骨干提取（教案）

代码总体框架

  • 骨干提取Thinning：
    • 骨干提取主函数：SequentialThinning。
    • 自研骨干提取算法：solve，配合三个检查函数“CheckEdge,CheckConnectivity和CheckEndPoint”。
    • Zhang-suen骨干提取算法：zhangsuen_solve，配合两个检查函数“UpDownCheck和LeftRightCheck”。
  • 特征点提取Feature Extraction：
    • FeatureExtraction：配合四个检查函数。新增Check_3BranchPoint和Check_4BranchPoint。
  • 提取后图像优化Optimization：
    • removeLinesUseFeatureExtraction: 利用特征点进行优化，配合ColorPixelDetect。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
class LWCThinning{
    public:
        Mat SequentialThinning(Mat& Image, int num_iterations);
        Mat FeatureExtraction(Mat& Image);
        //Optimization
        Mat removeLinesUseFeatureExtraction(Mat& Image);
    private:
        //VSP-slice based thinning
        Mat solve(Mat& Image, int iter_time);        
        int CheckEdge(Mat& image3x3);
        int CheckConnectivity(Mat& image3x3);
        int CheckEndPoint(Mat& image3x3);
        //zhang-suen algorithm
        Mat zhangsuen_solve(Mat& Image, int iter_time);
        int UpDownCheck(Mat& image3x3);
        int LeftRightCheck(Mat& image3x3);
        //Optimization functions
        int ColorPixelDetect(Mat& image3x3);
        //Pixel Analyze
        int Check_3BranchPoint(Mat& image3x3);
        int Check_4BranchPoint(Mat& image3x3);
};

时间复杂度上的考虑

本代码中除了对原始图像的遍历使用了双重for循环，在其他方法中对3x3的卷积核进行处理时，均没有使用for循环。这是出于时间复杂度上的考量。3x3的卷积核很小，一共也才9个元素，如果再用一个常规的双重for循环进行处理，时间复杂度将会增加O(n^2)，与外部遍历结合起来就是O(n^4),处理大图片将会非常缓慢；但相对地，如果将其展开书写，将时间复杂度降低为O(1)的情况下代码也不会变得过于冗杂。

三次判断

总体思路：三个check方法返回0或1。在solve()中调用这些check，三次判断全1才会执行像素的删除。

  • return 0: 该项check不通过，不可以删除。
  • return 1: 该项check通过，可以删除。

连续性判断 CheckConnectivity

先说结论：判断中心注目像素周围一圈的八个像素，从一个节点开始顺时针旋转遍历。如果像素值“从黑转变为白”的次数超过1次，那么就是不连续的。（另外，前提是周围一圈有两个及以上的黑色像素。如果只有一个黑像素，那么这个中心像素就是端点end-point）

我们可以通过画图来探索这个规律：

  • ①很好地展示了斜线连续性被破坏的情况。
  • ③很好地展示了直线连续性被破坏的情况。
  • ⑤展示了两条直线交点处连续性被破坏的情况。
  • ②中即便删除了中心像素，其连续性没有被破坏。

观察上述所有破快连续性的情况，我们发现他们周围一圈像素从黑变白的转变次数都大于一。(图中红色箭头和数字)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
int LWCThinning::CheckConnectivity(Mat& image3x3){
    //======== iff time which around 8-pixel transfer from B to W is only 1, then it will have connectivity. ========
    int B = 0;
    if(image3x3.at<uchar>(0,0) == 0)  B++;
    if(image3x3.at<uchar>(0,1) == 0)  B++;
    if(image3x3.at<uchar>(0,2) == 0)  B++;
    if(image3x3.at<uchar>(1,0) == 0)  B++;
    if(image3x3.at<uchar>(1,2) == 0)  B++;
    if(image3x3.at<uchar>(2,0) == 0)  B++;
    if(image3x3.at<uchar>(2,1) == 0)  B++;
    if(image3x3.at<uchar>(2,2) == 0)  B++;
    int T = 0;   //transfer from B to W 
    if(image3x3.at<uchar>(0,2) == 0 && image3x3.at<uchar>(0,1) == 255)    T++;    
    if(image3x3.at<uchar>(0,1) == 0 && image3x3.at<uchar>(0,0) == 255)    T++;
    if(image3x3.at<uchar>(0,0) == 0 && image3x3.at<uchar>(1,0) == 255)    T++;
    if(image3x3.at<uchar>(1,0) == 0 && image3x3.at<uchar>(2,0) == 255)    T++;
    if(image3x3.at<uchar>(2,0) == 0 && image3x3.at<uchar>(2,1) == 255)    T++;
    if(image3x3.at<uchar>(2,1) == 0 && image3x3.at<uchar>(2,2) == 255)    T++;
    if(image3x3.at<uchar>(2,2) == 0 && image3x3.at<uchar>(1,2) == 255)    T++;
    if(image3x3.at<uchar>(1,2) == 0 && image3x3.at<uchar>(0,2) == 255)    T++;

    if(B >= 2 && T == 1)
        return 1;   //preserve connectivity
    else
        return 0;
}

端点判断 CheckEndPoint

不难发现，如果中心注目像素的周围只存在一个黑色像素，那么它就是端点end-point。我们不能删除端点像素。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
int LWCThinning::CheckEndPoint(Mat& image3x3){
    //======== iff around 8-pixel only have one B, then this center pixel is endpoint.========
    int B = 0;
    if(image3x3.at<uchar>(0,0) == 0)  B++;
    if(image3x3.at<uchar>(0,1) == 0)  B++;
    if(image3x3.at<uchar>(0,2) == 0)  B++;
    if(image3x3.at<uchar>(1,0) == 0)  B++;
    if(image3x3.at<uchar>(1,2) == 0)  B++;
    if(image3x3.at<uchar>(2,0) == 0)  B++;
    if(image3x3.at<uchar>(2,1) == 0)  B++;
    if(image3x3.at<uchar>(2,2) == 0)  B++;

    if(B == 1)
        return 0;
    else
        return 1;  //not end point
}

孤立点判断 CheckEdge

这个是最简单的，只需要查看周围一圈是否有黑色像素即可。它的深层含义是：画面上孤立的黑色像素点不能删去，因为他完整地表征着一部分的画面信息。

1
2
3
4
5
6
7
8
int LWCThinning::CheckEdge(Mat& image3x3){
    if(image3x3.at<uchar>(0,0) == 0 || image3x3.at<uchar>(0,1) == 0 || image3x3.at<uchar>(0,2) == 0 ||   
       image3x3.at<uchar>(1,0) == 0 ||                                 image3x3.at<uchar>(1,2) == 0 ||
       image3x3.at<uchar>(2,0) == 0 || image3x3.at<uchar>(2,1) == 0 || image3x3.at<uchar>(2,2) == 0 )
        return 1;  //boundary have B-pixel
    else
        return 0;
}

slove()主方法 和 SequentialThinning()接口方法

  • 使用3x3的卷积核遍历原始图像。具体而言，当检测到3x3的中心是黑色像素时，从原图中抽取这个3x3的正方形像素区域构成一个新图像kernel_nowscanning，然后进行三次check。如果check通过，则将中心的黑色像素删除（即变为白色）。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
Mat LWCThinning::solve(Mat& Image, int iter_time){

    int img_width = Image.size().width;
    int img_height = Image.size().height;
    
    for(int y = 1; y < img_height-1; y++){
        for(int x = 1; x < img_width-1; x++){
            Mat kernel_nowscanning = Mat::zeros(3,3,CV_8UC1);   
            kernel_nowscanning.at<uchar>(0,0) = Image.at<uchar>(y-1,x-1);
            kernel_nowscanning.at<uchar>(0,1) = Image.at<uchar>(y-1,x);
            kernel_nowscanning.at<uchar>(0,2) = Image.at<uchar>(y-1,x+1);
            kernel_nowscanning.at<uchar>(1,0) = Image.at<uchar>(y,x-1);
            kernel_nowscanning.at<uchar>(1,1) = Image.at<uchar>(y,x);     //center-pixel
            kernel_nowscanning.at<uchar>(1,2) = Image.at<uchar>(y,x+1);
            kernel_nowscanning.at<uchar>(2,0) = Image.at<uchar>(y+1,x-1);
            kernel_nowscanning.at<uchar>(2,1) = Image.at<uchar>(y+1,x);
            kernel_nowscanning.at<uchar>(2,2) = Image.at<uchar>(y+1,x+1);
            //if center-pixel is black, go to 'mittsu no check', find out whether it can be change to white.
            int check1 = 0;
            int check2 = 0;
            int check3 = 0;
            if(kernel_nowscanning.at<uchar>(1,1) == 0){
                check1 = CheckEdge(kernel_nowscanning);
                check2 = CheckConnectivity(kernel_nowscanning);
                check3 = CheckEndPoint(kernel_nowscanning);
                //change center-pixel from B to W (in original-image)
                if(check1 == 1 && check2 == 1 && check3 == 1){
                    //kernel_nowscanning.at<uchar>(1,1) = 255;
                    Image.at<uchar>(y,x) = 255;
                }
            }
        }
    }

    string file_dic = "result/process_time" + to_string(iter_time) + ".png";
    imwrite(file_dic, Image);
    cout << "Image saved to " << file_dic << endl;

    return Image;
}
  • SequentialThinning()接口方法可以被用户(主函数)调用，对给定的图样进行N次骨干提取，将结果保存在result路径下。
1
2
3
4
5
6
7
8
9
10
void LWCThinning::SequentialThinning(Mat& Image, int num_iterations){
    Mat processedImage = Image.clone();
    for(int i=0; i < num_iterations; i++){
        processedImage = solve(processedImage, i);
        cout << "||||" << i << "||||";
    }
    imshow("Final Processed Image", processedImage);
    waitKey(0);
    cout << "Processing Completed." << endl;
}

初步测试结果

ABCD英文字符

  • 字母A：
  • 字母C：
  • 字母E：

复杂的彩色图样

首先将彩色图片读取为灰度图，然后进行二值化处理，最后再调用SequentialThinning。

1
2
3
4
5
6
7
8
9
10
int main(){
    //Complex image
    Mat image_ena = imread("source/ENANA.png", IMREAD_GRAYSCALE);
    Mat binary_image;
    double threshold_value = 128;
    threshold(image_ena, binary_image, threshold_value, 255, THRESH_BINARY);

    LWCThinning handle;
    handle.SequentialThinning(binary_image,5);
}

测试原图与处理后图像如下：

  • 原图：
  • 转化为二值化图像：
  • 我们把处理后的图像放大观察，可以看到骨干提取的一些细节。比如衣服的领带和眼睛等等：

问题：对于复杂图像，处理后的图像中存在很多垂直水平面的直线（噪声），如何去除？

优化一：“Zhang-Suen算法”

问题就出在我们对原始图像的处理顺序上。我们总是从上往下、从左往右地进行逐行扫描，这种顺序性可能导致水平或垂直方向的像素被优先移除或保留，从而产生水平或垂直的直线。这直接导致了很多倾斜的骨干被强制地转化为了垂直的直线。

从上方测试图“字母A”可以看出这种状况，并且在彩色测试图中也很明显能发现许多垂直直线。

Zhang-Suen算法

与先前自己实现的算法相比，其实就是check部分的逻辑不同。Zhang-Suen算法会先后进行两次遍历，一次是UpDown一次是LeftRight。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
//zhangsuen -----------------------------------------------------
int LWCThinning::UpDownCheck(Mat& image3x3){
    /* ---- */
    if((B >= 2 && B <= 6) && T == 1 && ((P2*P4*P6 == 0) || (P4*P6*P8 == 0)))     
        return 1;
    else
        return 0;
}

int LWCThinning::LeftRightCheck(Mat& image3x3){
    /* ---- */
    if((B >= 2 && B <= 6) && T == 1 && ((P2*P4*P8 == 0) || (P2*P6*P8 == 0)))     
        return 1;
    else
        return 0;
}

Mat LWCThinning::zhangsuen_solve(Mat& Image, int iter_time){
    int img_width = Image.size().width;
    int img_height = Image.size().height;
    
    for(int y = 1; y < img_height-1; y++){
        for(int x = 1; x < img_width-1; x++){
            Mat kernel_nowscanning = Mat::zeros(3,3,CV_8UC1);   
            kernel_nowscanning.at<uchar>(0,0) = Image.at<uchar>(y-1,x-1);
            kernel_nowscanning.at<uchar>(0,1) = Image.at<uchar>(y-1,x);
            kernel_nowscanning.at<uchar>(0,2) = Image.at<uchar>(y-1,x+1);
            kernel_nowscanning.at<uchar>(1,0) = Image.at<uchar>(y,x-1);
            kernel_nowscanning.at<uchar>(1,1) = Image.at<uchar>(y,x);     //center-pixel
            kernel_nowscanning.at<uchar>(1,2) = Image.at<uchar>(y,x+1);
            kernel_nowscanning.at<uchar>(2,0) = Image.at<uchar>(y+1,x-1);
            kernel_nowscanning.at<uchar>(2,1) = Image.at<uchar>(y+1,x);
            kernel_nowscanning.at<uchar>(2,2) = Image.at<uchar>(y+1,x+1);
            if(kernel_nowscanning.at<uchar>(1,1) == 0){
                int firstcheck = UpDownCheck(kernel_nowscanning);
                if(firstcheck == 1){
                    Image.at<uchar>(y,x) = 255;
                }
            }
        }
    }

    for(int y = 1; y < img_height-1; y++){
        for(int x = 1; x < img_width-1; x++){
            Mat kernel_nowscanning = Mat::zeros(3,3,CV_8UC1);   
            kernel_nowscanning.at<uchar>(0,0) = Image.at<uchar>(y-1,x-1);
            kernel_nowscanning.at<uchar>(0,1) = Image.at<uchar>(y-1,x);
            kernel_nowscanning.at<uchar>(0,2) = Image.at<uchar>(y-1,x+1);
            kernel_nowscanning.at<uchar>(1,0) = Image.at<uchar>(y,x-1);
            kernel_nowscanning.at<uchar>(1,1) = Image.at<uchar>(y,x);     //center-pixel
            kernel_nowscanning.at<uchar>(1,2) = Image.at<uchar>(y,x+1);
            kernel_nowscanning.at<uchar>(2,0) = Image.at<uchar>(y+1,x-1);
            kernel_nowscanning.at<uchar>(2,1) = Image.at<uchar>(y+1,x);
            kernel_nowscanning.at<uchar>(2,2) = Image.at<uchar>(y+1,x+1);
            if(kernel_nowscanning.at<uchar>(1,1) == 0){
                int secondcheck = LeftRightCheck(kernel_nowscanning);
                if(secondcheck == 1){
                    Image.at<uchar>(y,x) = 255;
                }
            }
        }
    }

    string file_dic = "result/process_time" + to_string(iter_time+1) + ".png";
    imwrite(file_dic, Image);
    cout << "Image saved to " << file_dic << endl;

    return Image;
}

优化结果

貌似没有解决垂直直线较多的问题…

优化二：利用“特性点提取”消除直线

刚刚在做“特征点提取”的工作时，突然来了灵感：我发现这些多余的直线都是由非特征点组成的，亦或是原理特征点的。换言之，我们只要提取了骨架图的特征点，然后去除远离特征点的点即可！

特征点提取

这是特征点提取后的图片，不同的特征点标上了不同的颜色：

  • 孤立点 isolated point：红色。
  • 端点 end point：绿色。
  • 3分歧点：蓝色。图片中的阶梯状斜线属于这一类。
  • 4分歧点：黄色。十字路口。

我们可以放大图片观察一下细节：（由于4分歧点的黄色标注看的不明显，我在边上放了一双眼睛。）

特征点提取的核心代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
int LWCThinning::Check_3BranchPoint(Mat& image3x3){
    /* ----- */
    if(B >= 3 && (T == 2 || T == 3)){
        return 1;
    }
    else{
        return 0;
    }
}

int LWCThinning::Check_4BranchPoint(Mat& image3x3){
    /* ----- */
    if(T == 4 || 
      (image3x3.at<uchar>(0,1) == 0 && image3x3.at<uchar>(1,0) == 0 && image3x3.at<uchar>(1,2) == 0 && image3x3.at<uchar>(2,1) == 0)){
        return 1;
    }
    else{
        return 0;
    }
}

Mat LWCThinning::FeatureExtraction(Mat& Image){

    Mat rgbImage;   //After PixelAnalyzing
    cvtColor(Image, rgbImage, COLOR_GRAY2BGR);

    int img_width = Image.size().width;
    int img_height = Image.size().height;
    
    for(int y = 1; y < img_height-1; y++){
        for(int x = 1; x < img_width-1; x++){
            Mat kernel_nowscanning = Mat::zeros(3,3,CV_8UC1);   
            kernel_nowscanning.at<uchar>(0,0) = Image.at<uchar>(y-1,x-1);
            kernel_nowscanning.at<uchar>(0,1) = Image.at<uchar>(y-1,x);
            kernel_nowscanning.at<uchar>(0,2) = Image.at<uchar>(y-1,x+1);
            kernel_nowscanning.at<uchar>(1,0) = Image.at<uchar>(y,x-1);
            kernel_nowscanning.at<uchar>(1,1) = Image.at<uchar>(y,x);     //center-pixel
            kernel_nowscanning.at<uchar>(1,2) = Image.at<uchar>(y,x+1);
            kernel_nowscanning.at<uchar>(2,0) = Image.at<uchar>(y+1,x-1);
            kernel_nowscanning.at<uchar>(2,1) = Image.at<uchar>(y+1,x);
            kernel_nowscanning.at<uchar>(2,2) = Image.at<uchar>(y+1,x+1);

            int check_iso = 0;       //isolated point
            int check_endpoint = 0;  //end point
            int check_3branchpoint = 0;  //3branchpoint
            int check_4branchpoint = 0;  //4branchpoint
            if(kernel_nowscanning.at<uchar>(1,1) == 0){
                check_iso = CheckEdge(kernel_nowscanning);
                check_endpoint = CheckEndPoint(kernel_nowscanning);
                check_3branchpoint = Check_3BranchPoint(kernel_nowscanning);
                check_4branchpoint = Check_4BranchPoint(kernel_nowscanning);
                if(check_iso == 0){
                    rgbImage.at<Vec3b>(y, x)[0] = 0;    
                    rgbImage.at<Vec3b>(y, x)[1] = 0;  
                    rgbImage.at<Vec3b>(y, x)[2] = 255;  //Red
                }
                if(check_endpoint == 0){
                    rgbImage.at<Vec3b>(y, x)[0] = 0;    
                    rgbImage.at<Vec3b>(y, x)[1] = 255;  //Green
                    rgbImage.at<Vec3b>(y, x)[2] = 0;  
                }
                if(check_3branchpoint == 1){
                    rgbImage.at<Vec3b>(y, x)[0] = 255;  //Blue
                    rgbImage.at<Vec3b>(y, x)[1] = 0;  
                    rgbImage.at<Vec3b>(y, x)[2] = 0;  
                }
                if(check_4branchpoint == 1){
                    rgbImage.at<Vec3b>(y, x)[0] = 0;  
                    rgbImage.at<Vec3b>(y, x)[1] = 255;  
                    rgbImage.at<Vec3b>(y, x)[2] = 255;  
                }
            }
        }
    }

    return rgbImage;
}

利用“特性点提取”消除直线

仍然使用3x3的扫描核kernel（中心像素为注目像素），如果kernel内不存在彩色像素，即不存在特征点，那么就删去这个中心像素。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
Mat LWCThinning::removeLinesUseFeatureExtraction(Mat& Image){

    int img_width = Image.size().width;
    int img_height = Image.size().height;
    
    for(int y = 1; y < img_height-1; y++){
        for(int x = 1; x < img_width-1; x++){
            Mat kernel_nowscanning = Mat::zeros(3,3,CV_8UC3);   
            kernel_nowscanning.at<Vec3b>(0,0) = Image.at<Vec3b>(y-1,x-1);
            kernel_nowscanning.at<Vec3b>(0,1) = Image.at<Vec3b>(y-1,x);
            kernel_nowscanning.at<Vec3b>(0,2) = Image.at<Vec3b>(y-1,x+1);
            kernel_nowscanning.at<Vec3b>(1,0) = Image.at<Vec3b>(y,x-1);
            kernel_nowscanning.at<Vec3b>(1,1) = Image.at<Vec3b>(y,x);     //center-pixel
            kernel_nowscanning.at<Vec3b>(1,2) = Image.at<Vec3b>(y,x+1);
            kernel_nowscanning.at<Vec3b>(2,0) = Image.at<Vec3b>(y+1,x-1);
            kernel_nowscanning.at<Vec3b>(2,1) = Image.at<Vec3b>(y+1,x);
            kernel_nowscanning.at<Vec3b>(2,2) = Image.at<Vec3b>(y+1,x+1);

            if(kernel_nowscanning.at<Vec3b>(1,1)[0] == 0 && kernel_nowscanning.at<Vec3b>(1,1)[1] == 0 && kernel_nowscanning.at<Vec3b>(1,1)[2] == 0){
                int is_coloured = ColorPixelDetect(kernel_nowscanning);
                if(!is_coloured){
                    //remove(change to W)
                    Image.at<Vec3b>(y,x)[0] = 255;
                    Image.at<Vec3b>(y,x)[1] = 255;
                    Image.at<Vec3b>(y,x)[2] = 255;
                }
            }
            
        }
    }
    return Image;
}
int LWCThinning::ColorPixelDetect(Mat& image3x3){
    for (int i = 0; i < 3; i++) {  
        for (int j = 0; j < 3; j++) {  
            Vec3b pixel = image3x3.at<Vec3b>(i, j);
            int B = pixel[0];
            int G = pixel[1];
            int R = pixel[2];
            //have colour
            if (!((B == 255 && G == 255 && R == 255) || (B == 0 && G == 0 && R == 0))) {  
                return 1;
            }
        }
    }
    return 0;
}

优化结果

动漫人物图片：ShinonomeEnana

  • 原图：
  • 初步骨架提取：
  • 特征点提取：
  • 利用特征点去直线：

环境景色图片：ひびきの南

  • 原图：
  • 初步骨架提取：
  • 特征点提取：
  • 利用特征点去直线：

包含大量文字的图片：日本年号对照表

  • 原图：
  • 初步骨架提取：
  • 特征点提取：

我们发现对文字进行处理后效果已经很好了，因此无需再继续优化。

backtop