OpenCV Thinning - eiichiromomma/CVMLAB GitHub Wiki

(OpenCV) OpenCVによる細線化(Thinning)

白黒画像(背景が黒)の細線化。OpenCVに関数が無い(?)ので作ってみた(かなり適当)

線抽出はKasvandの反復型線検出オペレータも参照。

OpenMPを使うと若干速くなる。

アルゴリズム

画像処理ハンドブックやWebに載っているので省略。

OpenCVらしいサンプルソース

2.x

//
//  main.cpp
//  Thinning
//
//  Created by Eiichiro Momma on 2014/08/11.
//  Copyright (c) 2014 Eiichiro Momma. All rights reserved.
//

#include <iostream>
#include <stdio.h>
#include <opencv2/opencv.hpp>

int main(int argc, const char * argv[])
{
    cv::Mat *kpb = new cv::Mat[8];
    cv::Mat *kpw = new cv::Mat[8];
    kpb[0]=(cv::Mat_<float>(3,3) << 1,1,0,1,0,0,0,0,0);
    kpb[1]=(cv::Mat_<float>(3,3) << 1,1,1,0,0,0,0,0,0);
    kpb[2]=(cv::Mat_<float>(3,3) << 0,1,1,0,0,1,0,0,0);
    kpb[3]=(cv::Mat_<float>(3,3) << 0,0,1,0,0,1,0,0,1);
    kpb[4]=(cv::Mat_<float>(3,3) << 0,0,0,0,0,1,0,1,1);
    kpb[5]=(cv::Mat_<float>(3,3) << 0,0,0,0,0,0,1,1,1);
    kpb[6]=(cv::Mat_<float>(3,3) << 0,0,0,1,0,0,1,1,0);
    kpb[7]=(cv::Mat_<float>(3,3) << 1,0,0,1,0,0,1,0,0);

    kpw[0]=(cv::Mat_<float>(3,3) << 0,0,0,0,1,1,0,1,0);
    kpw[1]=(cv::Mat_<float>(3,3) << 0,0,0,0,1,0,1,1,0);
    kpw[2]=(cv::Mat_<float>(3,3) << 0,0,0,1,1,0,0,1,0);
    kpw[3]=(cv::Mat_<float>(3,3) << 1,0,0,1,1,0,0,0,0);
    kpw[4]=(cv::Mat_<float>(3,3) << 0,1,0,1,1,0,0,0,0);
    kpw[5]=(cv::Mat_<float>(3,3) << 0,1,1,0,1,0,0,0,0);
    kpw[6]=(cv::Mat_<float>(3,3) << 0,1,0,0,1,1,0,0,0);
    kpw[7]=(cv::Mat_<float>(3,3) << 0,0,0,0,1,1,0,0,1);

    cv::Mat src = cv::imread("thinning_test.png",cv::IMREAD_GRAYSCALE);
    cv::Mat src_w(src.rows,src.cols, CV_32FC1);
    cv::Mat src_b(src.rows,src.cols, CV_32FC1);
    cv::Mat src_f(src.rows,src.cols, CV_32FC1);
    src.convertTo(src_f, CV_32FC1);
    src_f.mul(1./255.);
    cv::threshold(src_f, src_f, 0.5, 1.0, CV_THRESH_BINARY);
    cv::threshold(src_f, src_w, 0.5, 1.0, CV_THRESH_BINARY);
    cv::threshold(src_f, src_b, 0.5, 1.0, CV_THRESH_BINARY_INV);
    
    double sum=1;
    while (sum>0) {
        sum=0;
        for (int i=0; i<8; i++) {
            cv::filter2D(src_w, src_w, CV_32FC1, kpw[i]);
            cv::filter2D(src_b, src_b, CV_32FC1, kpb[i]);
            cv::threshold(src_w, src_w, 2.99, 1.0, CV_THRESH_BINARY);
            cv::threshold(src_b, src_b, 2.99, 1.0, CV_THRESH_BINARY);
            cv::bitwise_and(src_w, src_b, src_w);
            sum += cv::sum(src_w).val[0];
            cv::bitwise_xor(src_f, src_w, src_f);
            src_f.copyTo(src_w);
            cv::threshold(src_f, src_b, 0.5, 1.0, CV_THRESH_BINARY_INV);
        }
    }
    cv::imshow("Result", src_f);
    cv::waitKey(0);
    return 0;
    
}

Python (cv2)

# encoding: utf-8
# thinning.py
# Hilditch Thinning Algorithms
# 2012-3-6
# Eiichiro Momma
__author__ = 'momma'
import cv2
import numpy as np

def Init(kpw, kpb):
    kpw.append(np.array([[0.,0.,0.],[0.,1.,1.],[0.,1.,0.]]))
    kpw.append(np.array([[0.,0.,0.],[0.,1.,0.],[1.,1.,0.]]))
    kpw.append(np.array([[0.,0.,0.],[1.,1.,0.],[0.,1.,0.]]))
    kpw.append(np.array([[1.,0.,0.],[1.,1.,0.],[0.,0.,0.]]))
    kpw.append(np.array([[0.,1.,0.],[1.,1.,0.],[0.,0.,0.]]))
    kpw.append(np.array([[0.,1.,1.],[0.,1.,0.],[0.,0.,0.]]))
    kpw.append(np.array([[0.,1.,0.],[0.,1.,1.],[0.,0.,0.]]))
    kpw.append(np.array([[0.,0.,0.],[0.,1.,1.],[0.,0.,1.]]))
    kpb.append(np.array([[1.,1.,0.],[1.,0.,0.],[0.,0.,0.]]))
    kpb.append(np.array([[1.,1.,1.],[0.,0.,0.],[0.,0.,0.]]))
    kpb.append(np.array([[0.,1.,1.],[0.,0.,1.],[0.,0.,0.]]))
    kpb.append(np.array([[0.,0.,1.],[0.,0.,1.],[0.,0.,1.]]))
    kpb.append(np.array([[0.,0.,0.],[0.,0.,1.],[0.,1.,1.]]))
    kpb.append(np.array([[0.,0.,0.],[0.,0.,0.],[1.,1.,1.]]))
    kpb.append(np.array([[0.,0.,0.],[1.,0.,0.],[1.,1.,0.]]))
    kpb.append(np.array([[1.,0.,0.],[1.,0.,0.],[1.,0.,0.]]))

if __name__ == '__main__':
    kpw = []
    kpb = []
    Init(kpw, kpb)
    src = cv2.imread('thinning_test.jpg)',0)
    src_w = np.array(src, dtype=np.float32)/255.
    thresh, src_b = cv2.threshold(src_w, 0.5, 1.0, cv2.THRESH_BINARY_INV)
    thresh, src_f = cv2.threshold(src_w, 0.5, 1.0, cv2.THRESH_BINARY)
    thresh, src_w = cv2.threshold(src_w, 0.5, 1.0, cv2.THRESH_BINARY)
    th = 1.
    while th > 0:
        th = 0.
        for i in range(8):
            src_w = cv2.filter2D(src_w, cv2.CV_32F, kpw[i])
            src_b = cv2.filter2D(src_b, cv2.CV_32F, kpb[i])
            thresh, src_w = cv2.threshold(src_w, 2.99, 1, cv2.THRESH_BINARY)
            thresh, src_b = cv2.threshold(src_b, 2.99, 1, cv2.THRESH_BINARY)
            src_w = np.array(np.logical_and(src_w,src_b), dtype=np.float32)
            th += np.sum(src_w)
            src_f = np.array(np.logical_xor(src_f, src_w), dtype=np.float32)
            src_w = src_f.copy()
            thresh, src_b = cv2.threshold(src_f, 0.5, 1.0, cv2.THRESH_BINARY_INV)
            cv2.imshow('result', src_w)
            cv2.waitKey(10)
    cv2.imshow('result', src_f)
    cv2.waitKey(0)

APIを使ってそれらしい処理に。

同じサンプルでだいたい1秒。

  #include <cv.h>
  #include <highgui.h>

  void myThinningInit(CvMat** kpw, CvMat** kpb)
  {
    //cvFilter2D用のカーネル
    //アルゴリズムでは白、黒のマッチングとなっているのをkpwカーネルと二値画像、
    //kpbカーネルと反転した二値画像の2組に分けて畳み込み、その後でANDをとる
    for (int i=0; i<8; i++){
      *(kpw+i) = cvCreateMat(3, 3, CV_8UC1);
      *(kpb+i) = cvCreateMat(3, 3, CV_8UC1);
      cvSet(*(kpw+i), cvRealScalar(0), NULL);
      cvSet(*(kpb+i), cvRealScalar(0), NULL);
    }
    //cvSet2Dはy,x(row,column)の順となっている点に注意
    //kernel1
    cvSet2D(*(kpb+0), 0, 0, cvRealScalar(1));
    cvSet2D(*(kpb+0), 0, 1, cvRealScalar(1));
    cvSet2D(*(kpb+0), 1, 0, cvRealScalar(1));
    cvSet2D(*(kpw+0), 1, 1, cvRealScalar(1));
    cvSet2D(*(kpw+0), 1, 2, cvRealScalar(1));
    cvSet2D(*(kpw+0), 2, 1, cvRealScalar(1));
    //kernel2
    cvSet2D(*(kpb+1), 0, 0, cvRealScalar(1));
    cvSet2D(*(kpb+1), 0, 1, cvRealScalar(1));
    cvSet2D(*(kpb+1), 0, 2, cvRealScalar(1));
    cvSet2D(*(kpw+1), 1, 1, cvRealScalar(1));
    cvSet2D(*(kpw+1), 2, 0, cvRealScalar(1));
    cvSet2D(*(kpw+1), 2, 1, cvRealScalar(1));
    //kernel3
    cvSet2D(*(kpb+2), 0, 1, cvRealScalar(1));
    cvSet2D(*(kpb+2), 0, 2, cvRealScalar(1));
    cvSet2D(*(kpb+2), 1, 2, cvRealScalar(1));
    cvSet2D(*(kpw+2), 1, 0, cvRealScalar(1));
    cvSet2D(*(kpw+2), 1, 1, cvRealScalar(1));
    cvSet2D(*(kpw+2), 2, 1, cvRealScalar(1));
    //kernel4
    cvSet2D(*(kpb+3), 0, 2, cvRealScalar(1));
    cvSet2D(*(kpb+3), 1, 2, cvRealScalar(1));
    cvSet2D(*(kpb+3), 2, 2, cvRealScalar(1));
    cvSet2D(*(kpw+3), 0, 0, cvRealScalar(1));
    cvSet2D(*(kpw+3), 1, 0, cvRealScalar(1));
    cvSet2D(*(kpw+3), 1, 1, cvRealScalar(1));
    //kernel5
    cvSet2D(*(kpb+4), 1, 2, cvRealScalar(1));
    cvSet2D(*(kpb+4), 2, 2, cvRealScalar(1));
    cvSet2D(*(kpb+4), 2, 1, cvRealScalar(1));
    cvSet2D(*(kpw+4), 0, 1, cvRealScalar(1));
    cvSet2D(*(kpw+4), 1, 1, cvRealScalar(1));
    cvSet2D(*(kpw+4), 1, 0, cvRealScalar(1));
    //kernel6
    cvSet2D(*(kpb+5), 2, 0, cvRealScalar(1));
    cvSet2D(*(kpb+5), 2, 1, cvRealScalar(1));
    cvSet2D(*(kpb+5), 2, 2, cvRealScalar(1));
    cvSet2D(*(kpw+5), 0, 2, cvRealScalar(1));
    cvSet2D(*(kpw+5), 0, 1, cvRealScalar(1));
    cvSet2D(*(kpw+5), 1, 1, cvRealScalar(1));
    //kernel7
    cvSet2D(*(kpb+6), 1, 0, cvRealScalar(1));
    cvSet2D(*(kpb+6), 2, 0, cvRealScalar(1));
    cvSet2D(*(kpb+6), 2, 1, cvRealScalar(1));
    cvSet2D(*(kpw+6), 0, 1, cvRealScalar(1));
    cvSet2D(*(kpw+6), 1, 1, cvRealScalar(1));
    cvSet2D(*(kpw+6), 1, 2, cvRealScalar(1));
    //kernel8
    cvSet2D(*(kpb+7), 0, 0, cvRealScalar(1));
    cvSet2D(*(kpb+7), 1, 0, cvRealScalar(1));
    cvSet2D(*(kpb+7), 2, 0, cvRealScalar(1));
    cvSet2D(*(kpw+7), 1, 1, cvRealScalar(1));
    cvSet2D(*(kpw+7), 1, 2, cvRealScalar(1));
    cvSet2D(*(kpw+7), 2, 2, cvRealScalar(1));
  }

  void main(void)
  {
    //白黒それぞれ8個のカーネルの入れ物
    CvMat** kpb = new CvMat *[8];
    CvMat** kpw = new CvMat *[8];
    myThinningInit(kpw, kpb);
    IplImage* src = cvLoadImage("thinning_test.jpg",0);
    IplImage* dst = cvCloneImage(src);
    //32Fの方が都合が良い
    IplImage* src_w = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, 1);
    IplImage* src_b = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, 1);
    IplImage* src_f = cvCreateImage(cvGetSize(src), IPL_DEPTH_32F, 1);
    cvScale(src, src_f, 1/255.0, 0);
    //原画像を2値化(しきい値は用途に合わせて考える)
    //src_f:2値化した画像(32F)
    //src_w:作業バッファ
    //src_b:作業バッファ(反転)
    cvThreshold(src_f,src_f,0.5,1.0,CV_THRESH_BINARY);
    cvThreshold(src_f,src_w,0.5,1.0,CV_THRESH_BINARY);
    cvThreshold(src_f,src_b,0.5,1.0,CV_THRESH_BINARY_INV);
    //デバッグ用
    //cvNamedWindow("src",1);
    //cvShowImage("src",src);
    //1ターンでマッチしてなければ終了
    double sum=1;
    while(sum>0){
      sum=0;
      for (int i=0; i<8; i++){
        cvFilter2D(src_w, src_w, *(kpw+i));
        cvFilter2D(src_b, src_b, *(kpb+i));
        //各カーネルで注目するのは3画素ずつなので、マッチした注目画素の濃度は3となる
        //カーネルの値を1/9にしておけば、しきい値は0.99で良い
        cvThreshold(src_w,src_w,2.99,1,CV_THRESH_BINARY); //2.5->2.99に修正
        cvThreshold(src_b,src_b,2.99,1,CV_THRESH_BINARY); //2.5->2.99
        cvAnd(src_w, src_b, src_w);
        //この時点でのsrc_wが消去候補点となり、全カーネルで候補点が0となった時に処理が終わる
        sum += cvSum(src_w).val[0];
        //原画像から候補点を消去(二値画像なのでXor)
        cvXor(src_f, src_w, src_f);
        //作業バッファを更新
        cvCopyImage(src_f, src_w);
        cvThreshold(src_f,src_b,0.5,1,CV_THRESH_BINARY_INV);
      }
    }
    //8Uの画像に戻して表示
    cvConvertScaleAbs(src_f, dst, 255, 0);
    cvNamedWindow("dst",1);
    cvShowImage("dst", dst);
    cvWaitKey(0);
  }

処理の結果

処理の過程

処理前の画像

処理後の画像

旧サンプルソース

アルゴリズムに忠実に作った場合。 やっている事はひたすら比較なので、かなり泥臭い。

  #include <cv.h>
  #include <highgui.h>
  #include <stdio.h>

  //白黒パターン埋め込み泥臭い関数
  //部分入れ替えは頭が混乱するので全取り替え
  void setPattern(CvPoint* w3p, CvPoint* b3p, int *pCount)
  {
    if (*pCount == 8){
      *pCount=1;
    }else{
      (*pCount)++;
    }
    switch(*pCount){
      case 1:
        w3p[0]=cvPoint(1,1);w3p[1]=cvPoint(2,1);w3p[2]=cvPoint(1,2);
        b3p[0]=cvPoint(0,0);b3p[1]=cvPoint(1,0);b3p[2]=cvPoint(0,1);
        break;
      case 2:
        b3p[0]=cvPoint(0,0);b3p[1]=cvPoint(1,0);b3p[2]=cvPoint(2,0);
        w3p[0]=cvPoint(1,1);w3p[1]=cvPoint(0,2);w3p[2]=cvPoint(1,2);
        break;
      case 3:
        b3p[0]=cvPoint(1,0);b3p[1]=cvPoint(2,0);b3p[2]=cvPoint(2,1);
        w3p[0]=cvPoint(0,1);w3p[1]=cvPoint(1,1);w3p[2]=cvPoint(1,2);
        break;
      case 4:
        b3p[0]=cvPoint(2,0);b3p[1]=cvPoint(2,1);b3p[2]=cvPoint(2,2);
        w3p[0]=cvPoint(0,0);w3p[1]=cvPoint(0,1);w3p[2]=cvPoint(1,1);
        break;
      case 5:
        b3p[0]=cvPoint(2,1);b3p[1]=cvPoint(2,2);b3p[2]=cvPoint(1,2);
        w3p[0]=cvPoint(1,0);w3p[1]=cvPoint(1,1);w3p[2]=cvPoint(0,1);
        break;
      case 6:
        b3p[0]=cvPoint(0,2);b3p[1]=cvPoint(1,2);b3p[2]=cvPoint(2,2);
        w3p[0]=cvPoint(2,0);w3p[1]=cvPoint(1,0);w3p[2]=cvPoint(1,1);
        break;
      case 7:
        b3p[0]=cvPoint(0,1);b3p[1]=cvPoint(0,2);b3p[2]=cvPoint(1,2);
        w3p[0]=cvPoint(1,0);w3p[1]=cvPoint(1,1);w3p[2]=cvPoint(2,1);
        break;
      case 8:
        b3p[0]=cvPoint(0,0);b3p[1]=cvPoint(0,1);b3p[2]=cvPoint(0,2);
        w3p[0]=cvPoint(1,1);w3p[1]=cvPoint(2,1);w3p[2]=cvPoint(2,2);
        break;
    }
  }

  //ROIとして取得した3x3の画像と白黒テーブルを比較、一致すれば1を返す
  int myMatching(IplImage *win,CvPoint *w3p,CvPoint *b3p)
  {
    if (((unsigned char*)(win->imageData+w3p[0].y*win->widthStep))[w3p[0].x] ==255 &&
      ((unsigned char*)(win->imageData+w3p[1].y*win->widthStep))[w3p[1].x] ==255 &&
      ((unsigned char*)(win->imageData+w3p[2].y*win->widthStep))[w3p[2].x] ==255 &&
      ((unsigned char*)(win->imageData+b3p[0].y*win->widthStep))[b3p[0].x] ==0 &&
      ((unsigned char*)(win->imageData+b3p[1].y*win->widthStep))[b3p[1].x] ==0 &&
      ((unsigned char*)(win->imageData+b3p[2].y*win->widthStep))[b3p[2].x] ==0 ){
        return 1;
    }
    return 0;
  }

  int main(void)
  {
    //白黒テーブル
    CvPoint *white3Points;
    CvPoint *black3Points;
    //8パターンを1ターンとするカウント
    int patternCount=0;
    int turnCount=1;
    IplImage *src;
    IplImage *dst;
    IplImage *ROIimg;
    IplImage *bufferimg;
    char fname[1024];
    long isChanged=1;
    int x,y;

    dst=cvLoadImage("thinning_test.png",0);
    src=cvCloneImage(dst);
    //作業用画像
    cvThreshold(dst,src,128,255,CV_THRESH_BINARY);
    cvCopy(src,dst,NULL);
    bufferimg=cvCloneImage(src);
    //3x3のROIコピー用
    ROIimg=cvCreateImage(cvSize(3,3),IPL_DEPTH_8U,1);

    white3Points = (CvPoint*)cvAlloc(sizeof(CvPoint)*3);
    black3Points = (CvPoint*)cvAlloc(sizeof(CvPoint)*3);
    //終了条件は1ターン(8パターン)完了時で変更が無かった場合
    while(patternCount !=8 || isChanged>0){
      //1ターンでカウントをリセット
      if (patternCount==8){
        turnCount++;
        isChanged=0;
      }
      //パターン変更
      setPattern(white3Points,black3Points,&patternCount);
      for (y=0; y<dst->height-3; y++){
        for (x=0; x<dst->width-3; x++){
          //ROIとして3x3を切り出し
          cvSetImageROI(dst,cvRect(x,y,3,3));
          cvCopy(dst,ROIimg,NULL);
          //マッチング関数に放り込み、一致してれば3x3の中心に相当するbufferimgの画素を0に
          if (myMatching(ROIimg,white3Points,black3Points)){
            isChanged++;
            ((unsigned char*)(bufferimg->imageData+(y+1)*bufferimg->widthStep))[x+1] = 0;
          }
          cvResetImageROI(dst);
        }
      }
      //デバッグ用画像
      //各パターンで削られていく様子が分かる
      sprintf(fname,"thinning_turn%02d_pattern%d.jpg",turnCount,patternCount);
      printf("turn%02d, pattern%d\n",turnCount,patternCount);
      cvSaveImage(fname,bufferimg);
      //各パターンでのスキャン終了時にbufferimg->dstを行なう
      cvCopy(bufferimg,dst,NULL);
    }
    cvCopy(bufferimg,dst,NULL);
    cvNamedWindow("src",1);
    cvNamedWindow("thinning",1);
    cvShowImage("src",src);
    cvShowImage("thinning",dst);
    cvWaitKey(0);
    cvSaveImage("thinning_result.jpg",dst);
    cvReleaseImage(&src);
    cvReleaseImage(&dst);
    cvReleaseImage(&ROIimg);
    cvDestroyWindow("src");
    cvDestroyWindow("thinning");
    cvFree(&white3Points);
    cvFree(&black3Points);
    return 0;
  }
⚠️ **GitHub.com Fallback** ⚠️