Obtener esqueleto de imagen binaria

Se le llama esqueleto de una imagen a un conjunto de curvas centradas que surge de reducir la forma original. El cálculo del esqueleto es una herramienta de análisis no escalar de formas, que conserva las propiedades topológicas de la forma original así como las propiedades geométricas, según el método utilizado. El esqueleto también es conocido como eje medio y tiene diversas aplicaciones, por ejemplo: reconocimiento de letras, números o símbolos, identificación de huellas dactilares, y muchas otras.

En color verde podemos ver el esqueleto de una estrella.

Para esta tarea usaremos los algoritmos de adelgazamiento propuestos por Zhang Suen y Guo Hall ambos son similares con pequeñas variaciones, las implementaciones de dichos algoritmos las he tomado de el módulo ximgproc por lo que si tienes este módulo compilado en tu versión OpenCV puedes usarlo directamente.

En esta web puedes encontrar una explicación del algoritmo de Zhang Suen: http://www.tecnohobby.net/ppal/index.php/vision-computacional/topicos-generales/2-algoritmo-de-zhang-suen

El código C++ para ambos algoritmos es el siguiente:

// Applies a thinning iteration to a binary image
static void thinningIteration(Mat img, int iter, int thinningType) {
 Mat marker = Mat::zeros(img.size(), CV_8UC1);

 if (thinningType == THINNING_ZHANGSUEN) {
  for (int i = 1; i < img.rows - 1; i++)
  {
   for (int j = 1; j < img.cols - 1; j++)
   {
    uchar p2 = img.at<uchar>(i - 1, j);
    uchar p3 = img.at<uchar>(i - 1, j + 1);
    uchar p4 = img.at<uchar>(i, j + 1);
    uchar p5 = img.at<uchar>(i + 1, j + 1);
    uchar p6 = img.at<uchar>(i + 1, j);
    uchar p7 = img.at<uchar>(i + 1, j - 1);
    uchar p8 = img.at<uchar>(i, j - 1);
    uchar p9 = img.at<uchar>(i - 1, j - 1);

    int A = (p2 == 0 && p3 == 1) + (p3 == 0 && p4 == 1) +
     (p4 == 0 && p5 == 1) + (p5 == 0 && p6 == 1) +
     (p6 == 0 && p7 == 1) + (p7 == 0 && p8 == 1) +
     (p8 == 0 && p9 == 1) + (p9 == 0 && p2 == 1);
    int B = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9;
    int m1 = iter == 0 ? (p2 * p4 * p6) : (p2 * p4 * p8);
    int m2 = iter == 0 ? (p4 * p6 * p8) : (p2 * p6 * p8);

    if (A == 1 && (B >= 2 && B <= 6) && m1 == 0 && m2 == 0)
     marker.at<uchar>(i, j) = 1;
   }
  }
 }
 if (thinningType == THINNING_GUOHALL) {
  for (int i = 1; i < img.rows - 1; i++)
  {
   for (int j = 1; j < img.cols - 1; j++)
   {
    uchar p2 = img.at<uchar>(i - 1, j);
    uchar p3 = img.at<uchar>(i - 1, j + 1);
    uchar p4 = img.at<uchar>(i, j + 1);
    uchar p5 = img.at<uchar>(i + 1, j + 1);
    uchar p6 = img.at<uchar>(i + 1, j);
    uchar p7 = img.at<uchar>(i + 1, j - 1);
    uchar p8 = img.at<uchar>(i, j - 1);
    uchar p9 = img.at<uchar>(i - 1, j - 1);

    int C = ((!p2) & (p3 | p4)) + ((!p4) & (p5 | p6)) +
     ((!p6) & (p7 | p8)) + ((!p8) & (p9 | p2));
    int N1 = (p9 | p2) + (p3 | p4) + (p5 | p6) + (p7 | p8);
    int N2 = (p2 | p3) + (p4 | p5) + (p6 | p7) + (p8 | p9);
    int N = N1 < N2 ? N1 : N2;
    int m = iter == 0 ? ((p6 | p7 | (!p9)) & p8) : ((p2 | p3 | (!p5)) & p4);

    if ((C == 1) && ((N >= 2) && ((N <= 3)) & (m == 0)))
     marker.at<uchar>(i, j) = 1;
   }
  }
 }

 img &= ~marker;
}

// Apply the thinning procedure to a given image
void thinning(InputArray input, OutputArray output, int thinningType) {
 Mat processed = input.getMat().clone();
 // Enforce the range of the input image to be in between 0 - 255
 processed /= 255;

 Mat prev = Mat::zeros(processed.size(), CV_8UC1);
 Mat diff;

 do {
  thinningIteration(processed, 0, thinningType);
  thinningIteration(processed, 1, thinningType);
  absdiff(processed, prev, diff);
  processed.copyTo(prev);
 } while (countNonZero(diff) > 0);

 processed *= 255;

 output.assign(processed);
}

Como podemos ver y como se explica en la web antes mencionada este es un algoritmo iterativo que busca adelgazar la figura hasta obtener el respectivo esqueleto de grosor 1 pixel, para especificar el método a utilizar usamos la thinningType con valor igual a 1 usaremos Zhang Suen y con valor igual a 2 Guo Hall.

int main(int argc, char** argv)
{
 Mat image = imread("../star.jpg", CV_LOAD_IMAGE_GRAYSCALE);

 // clonar imagen original
 Mat src = image.clone();

 // obtener el esqueleto
 thinning(src, src, 2);

 // 1. dibujar el esqueleto a color verde sobre la imagen original
 cvtColor(image, image, CV_GRAY2BGR);

 for (int i = 0; i < image.cols; i++)
 {
  for (int j = 0; j < image.rows; j++)
  {
   Scalar intensity = src.at<uchar>(j, i);
   if (intensity.val[0] == 255) {
    image.at<Vec3b>(j, i) = Vec3b(0, 255, 0);
   }
  }
 }
 // end 1

 imshow("OpenCV Skeleton Final", src);
 imshow("Original Skeleton Final", image);
 waitKey(0);

 return 0;
}

El proceso de adelgazamiento para llegar al eje medio se puede apreciar en la siguiente animación:

Otro método utilizado para calcular el esqueleto de una imagen es mediante operaciones morfológicas de erosión y dilatación, puedes ver el siguiente enlace: http://felix.abecassis.me/2011/09/opencv-morphological-skeleton/ 

• Descargar código fuente: opencv-esqueleto.zip
• Ver en GitHub: opencv-esqueleto-src