Mathematica's morphological thinning algorithm

Question

Could anyone cite some papers that describe Mathematica's morphological thinning algorithm?

Answer 1

Prologue

Some five years ago I have asked exactly this question to the Wolfram support people. Below I have taken their respective answers (one sentence for each Method) but have added a lot of further reading. Finally, in an Add-On I demonstrate my own implementation of Rosenfeld's 1971 variant of a 2D thinning algorithm, in order to let you compare a few methods.

Guo/Hall:

Method -> Morphological implements the algorithm described in: Lam L, Lee SW, and Suen CY, "Thinning Methodologies - A Comprehensive Survey," IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol 14, No. 9, September 1992, page 879, bottom of first column through top of second column

However, as this paper is a survey, it might be better to use the original reference:

Guo Z and Hall RW: "Parallel Thinning with Two-subiteration Algorithms." Communications of the ACM, Vol. 32, No. 3, March 1989, see section 3 starting on page 360 entitled "A TWO-SUBITERATION THINNING ALGORITHM" (referred to as Algorithm 1 throughout the paper)

Besides, there are further roots cited by Guo and Hall:

Zhang TY and Suen CY: "A Fast Parallel Algorithm for Thinning Digital Patterns." Communications of the ACM, Vol 27, No. 3, March 1984

Lü, HE and Wang PSP: "A Comment on 'A Fast Parallel Algorithm for Thinning Digital Patterns'" Communications of the ACM, Vol 29, No. 3, March 1986

Gonzalez/Woods:

Method -> MedialAxis is described in Gonzalez/Woods "Digital image processing", in a section named like "morphological medial axis by thinning".

In fact, in the first edition from 1992 for the medial axis computation the description is given in 8.1.5 "The Skeleton of a Region", in the second edition from 2001 you find it in under 11.1.5 "Skeletons" (sorry, the first two paragraphs are missing there), in the third edition as of 2010 it is given in 11.1.7 "Skeletons". Let's see where the authors will include it in the upcoming fourth edition :-)

As the copy referenced above was somehow cropped, here is the neighborhood as is defined by Gonzalez/Woods:

{{p9, p2, p3}, {p8, p1, p4}, {p7, p6, p5}}

Anyway, Rafael Gonzalez and Richard Woods cite Harry Blum who originally has proposed the concept of medial axis transform back in 1967 (naming it as "medial axis function - MAF"):

Harry Blum: "A Transformation for Extracting New Descriptors of Shape", in W. Wathen-Dunn (ed.), Proc. Models for the Perception of Speech and Visual Form, pp. 362-380, MIT Press, Cambridge, MA, November 1967.

Since Harry Blum is rather giving an abstract description of the general properties of the MAF and is not proposing a certain computational algorithm to obtain the MAF, and as Gonzalez and Woods do not cite anything else here, we need to assume that the given medial axis algorithm is their original work.

Add-On: Stefanelli/Rosenfeld

I was just thinking if I could please some of you by adding my own implementation of the good old Rosenfeld thinning. It was published decades ago:

Stefanelli R and Rosenfeld A: "Some Parallel Thinning Algorithms for Digital Pictures", Journal of the ACM, Vol. 18, No. 2, April 1971

Clear[Rosenfeld];
Rosenfeld[image_Image, nh_: 8, borders_: "NSEW"] :=
  Module[{imagedata, table, table4, table8, coords, deletions,
    deletionindices, northshift, southshift, eastshift, westshift, 
    directionalsequence, check}, 
   imagedata = ArrayPad[ImageData[image, "Bit"], 1];
   northshift = {-1, 0};
   southshift = {1, 0};
   eastshift = {0, 1};
   westshift = {0, -1};
   directionalsequence = {northshift, southshift, westshift, eastshift};
   directionalsequence = 
    directionalsequence[[Flatten[StringPosition["NSEW", #] & /@ 
         StringCases[borders, RegularExpression["[NSEW]"]], 1][[All, 1]]]];
   table4 = 
    IntegerDigits[
     53069725914310846919685767315128399372436785243025614538690558020\
      5726755357445973110133053045854705883358494050138831185007206823686093\
      00909165535, 2, 512];
   table8 = 
    IntegerDigits[
     15543186043038566186790018003123690632026825017137796030838268923\
      1709659598404465396940268644806626670319011047096356514674276473698052\
      70924042222591, 2, 512];
   Switch[nh, 4, table = table4, 8, table = table8, _, 
    Print["The neighborhood " <> ToString[nh] <> 
       " is not permitted, use either 4 or 8 (8 is preset)."];
   ];
   check = {{1}, {1}, {1}, {1}};
   While[Flatten[check] != {}, 
    Map[Function[ds,  coords = Position[imagedata - RotateLeft[imagedata, ds], 1];
       deletions = 
        Map[table[[FromDigits[Flatten[imagedata[[{#[[1]] - 1, #[[1]], #[[1]] + 1},
                     {#[[2]] - 1, #[[2]], #[[2]] + 1}]]], 2] + 1]] &, coords]; 
       deletionindices = Flatten[Position[deletions, 1]];
       imagedata = ReplacePart[imagedata, coords[[deletionindices]] -> 0];
       check[[1]] = deletionindices;
       check = RotateRight[check, 1];
      ],
      directionalsequence];
   ];
   Image[ArrayPad[imagedata, -1, "Bit"]]
  ];

If you try to copy this code above, please make sure that no \[Times] symbols are being inserted in the huge numbers when pasting into your Mathematica front end.

To be honest, it is not parallelized, to do so, one would have to replace the inner of the two nested Map inside the While statement with ParallelMap, however, the overhead data traffic is too high for typical images. If someone has an advice for a better parallelization, let me know. Also, please note that a single channel binary image is being expected. As you see I am using two different look-up tables to find the simple points in either 4- or 8-neighborhoods.

Some Examples

Let's see how it works:

HighlightImage[#, Rosenfeld@#] &@
 Binarize@
  ColorNegate@
   ImageResize[
    ImagePad[
     Import["http://www.carstylisten.de/img/Leipzig_Skyline_einzeln.jpg"], -1],
    Scaled[0.5]]

Here the variant with a 4-neighborhood (ensuring the skeleton's pixels are always edge-connected):

HighlightImage[#, Rosenfeld[#, 4]] &@
 Binarize@
  ColorNegate@
   ImageResize[
    ImagePad[
     Import["http://www.carstylisten.de/img/Leipzig_Skyline_einzeln.jpg"], -1],
    Scaled[0.5]]

Please compare with Thinning, first using the Guo/Hall variant:

HighlightImage[#, Thinning[#, Method -> "Morphological"]] &@
 Binarize@
  ColorNegate@
   ImageResize[
    ImagePad[
     Import["http://www.carstylisten.de/img/Leipzig_Skyline_einzeln.jpg"], -1],
    Scaled[0.5]]

Finally, here is the Gonzalez/Woods variant

HighlightImage[#, Thinning[#, Method -> "MedialAxis"]] &@
 Binarize@
  ColorNegate@
   ImageResize[
    ImagePad[
     Import["http://www.carstylisten.de/img/Leipzig_Skyline_einzeln.jpg"], -1],
    Scaled[0.5]]

Now it is up to you to decide what do you like best...

Related Question

Has anybody some 3D thinning algorithm implemented in Mathematica?