# An operator for edge reconnection and skeletonization

Bernard Allouche

► To cite this version:

Bernard Allouche. An operator for edge reconnection and skeletonization. Acta Stereologica, 1992, 11 (2), pp.175-184. hal-05063140

HAL Id: hal-05063140

https://hal.science/hal-05063140v1

Submitted on 12 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# AN OPERATOR FOR EDGE RECONNECTION AND SKELETONIZATION

Bernard Allouche

Ecole d'Architecture de Lyon
3 rue Maurice Audin, Boîte Postale 170
69512 Vaulx en Velin, France

## ABSTRACT

On the basis of an analogy between classical image processing problems and a simple problem of strength of materials, an operator called $\Delta^{-1}$ is introduced. The methodology is outlined and its advantages presented, starting by an elliptic linear partial differential equation in the form $\Delta A = B$, where $\Delta$ represents the Laplacian. The operator's properties will not be described with mathematical formalism, but will be applied to a series of image processing problems. They will also be compared with those of other morphological operators, in order to prove their advantages.

**Keywords** image analysis, Laplacian, skeletonization, edge reconnection

## INTRODUCTION

Edge extractions achieved in a grey tone image with a classical operator, such as gradient, Laplacian or Sobel filters, are generally discontinuous and thick. This process leads to many difficulties, particularly for pattern recognition or edge reconnection and vectorization applications. An example is given in Babu and Nevatia (1980).

The numerous methods described in the literature for edge thinning and reconnection can be grouped into two families. The first one needs morphological operators, whereas the second one is based on the nearest edges in a small neighbourhood (the notion of proximity is here linked with the Euclidean distance but also concerns the edge itself, the alignment of pixels and the degree of parallelism in their alignments ... ). Of course, both families have advantages and drawbacks. The methods belonging to the "morphologic" family are rather global, and relatively systematic and quick. On the other hand, some configurations cause mistakes as shown in figure 1.
<figure><img src="image_1.png" /><figcaption>Figure 1 - a reconnection based on mathematical morphology.</figcaption></figure>

Both rectangles are subjected to mathematical morphology operations but proper results are obtained only in the first case. A process using the second family of methods could take into account edge directions, alignments, or other properties, so that case 2 of figure 1 could be easily solved.

A set of common inadequacies can be listed for the two families

* the difficulty in determining the appropriate distance between two pixels to avoid their reconnection,

* the high sensitivity to noise,

* the inappropriateness of grey tone images.

The purpose of this article is to present a natural way to soften these drawbacks. Some vocabulary is informally introduced leading to analogies that will be investigated more thoroughly afterwards. Some pixels have to be added during edge reconnection to "knit" the segments together. Consequently, a pixel "joint" is defined as a pixel having "the highest potential" of belonging to the examined edge. Figure 2 gives an illustration of such "edges" pixels (dotted) and "joints" pixels (hatched).

<figure><img src="image_2.png" /><figcaption>Figure 2 - the pixels "joint".</figcaption></figure>

The "potential" mentioned above, should take into account the relative proximity of the segments, their alignment and orientation, and so on... Next, the hatched pixels in figure 2 indicate the line along which the potential's value is the highest from one segment to another. In the following part, this potential will be presented with more details.
ANALOGIES WITH SOME PROBLEMS OF PHYSICS

The pixels corresponding to edges of binary images can be considered as a set of points with masses µ(x,y,z) distributed on a rectangular plane. There results a gravitational potential G, given by the equation $ΔG(x,y,z) = -4π g . μ(x,y,z)$, where g is a constant.

An interpretation is also possible in terms of electrostatics. The pixels are conceived as electrical particles $q(x,y,z)$. This case is governed by $ΔV(x,y,z) = k.q(x,y,z)$, where V is the electrostatic potential and k is a constant.

Numerous other problems of physics are also based on the two- or three-dimensional version of Poisson's equation $ΔA = B$. In the given examples, $μ, q, G$ and $V$ are functions of three real coordinates. They do not correspond to the pixel's potential, whose nature is two-dimensional. A stationary problem of mechanics governed by Poisson's equation is therefore selected: the bending of a loaded thin plate.

Structures have received major interest in engineering for a long time. Numerous models have been derived from experiments to meet the needs of the specialists in strength of materials. The selected problem, mentioned previously, concerns a rectangular, thin elastic plate, which is supposed to have also a constant thickness and isotropic properties. This plate is then subjected to a load $q(x,y,z)$ and undergoes small deflections, as shown in figure 3.

A mathematical justification of the classical approach is proposed in Ciarlet and Destuynder (1979), which gives the following fourth order differential equation:

$$ \Delta\Delta F(x,y) = \frac{-h^3 Eq(x,y)}{12(1-v^2)} \quad [1] $$

where v is the Poisson's coefficient, E is Young's modulus of the plate's material, h the thickness, $q(x,y)$ is the load at each point, and $F(x,y)$ the elastic deflection of the plate.

Using appropriate coordinates and units for E, v and h (see figure 3), this equation can be split into the following system:

$$ \begin{cases} \Delta F(x,y) = M(x,y) \\ \Delta M(x,y) = q(x,y) \end{cases} \quad [2] $$

where $M(x,y)$ is an unknown quantity having the dimensions of a bending moment.

A condition must be added concerning the boundary of the plate. To that end, it is assumed that $M(x,y)=0$ holds on the boundary of the plate. Actually, $M(x,y)$ can be associated with a flexural couple from which bending stresses can be derived.

If the thin plate cannot resist the load, it starts to crack where the stresses resulting from $M(x,y)$ are highest. Gradually, the fissures will propagate in the direction of the highest breaking "potential", thereby following the lines of highest
bending moment. Even if the reality is a bit more complicated, this description of fissuration enables the establishment of the parallel between technical mechanics and image processing. Indeed, looking at figure 2, one can imagine the dotted pixels to be fissures. Upon propagating, the fissures will reach hatched pixels.

<figure><img src="image_3.png" /><figcaption>Figure 3 - appropriate axis and units for a rectangular plate.</figcaption></figure>

## APPLICATION TO IMAGE PROCESSING

The parallelism between thin plate cracking and edge reconnection provides a model for a pixel's potential. A grey tone function $g(x,y)$ exists for every image, so that a potential function $p(x,y)$ can be defined for each pixel by:

$$ \Delta p(x,y) = g(x,y) \quad [3] $$

with $\{x,y\}$ in $[0,1]$ and $p(x,y)=0$ on the boundary.

In this equation, $g(x,y)$ takes the place of the load $q(x,y)$ in equation [2], while $p(x,y)$ replaces $M(x,y)$. Next, equation [3] will be solved.

The unicity of the solution of [3] (given $g(x,y)$) is equivalent to the following definition of the operator $\Delta^{-1}$:

$$ p(x,y) = \Delta^{-1} g(x,y) \quad [4] $$

Generally, the solution of equation [3] is approximated by a method of finite differences. This method uses Taylor's formula, which gives for a size of discretization equal to $(N+1)^2$ (number of lines and rows of the considered image):

$$ (4p_{ij} - p_{i,j-1} - p_{i,j+1} - p_{i-1,j} - p_{i+1,j}) \cdot N = g_{ij} \quad [5] $$

for $(i,j)$ from 1 to $N-1$.

In this set of equations, the discretization step is $1/N$, $g_{ij}$ is the value of pixel's grey level and $p_{ij}$ is its potential. This can be transformed into unit steps by changing the units. As a result, the following system is obtained:

$$ 4p_{ij} - p_{i,j-1} - p_{i,j+1} - p_{i-1,j} - p_{i+1,j} = g_{ij} \quad [6] $$

for every $(i,j)$ from 1 to $N-1$.
(N-1)² equations with (N-1)² unknown quantities exist, which are solved following Jacobi's iterative method:

$$4p_{ij}^{(k+1)} = p_{ij-1}^{(k)} + p_{i-1j}^{(k)} + p_{i+1j}^{(k)} + p_{ij+1}^{(k)} + g_{ij} \quad [7]$$

where $p_{ij}^{(k)}$ is the kth iterated value of $p_{ij}$.

Ciarlet (1988) describes this method in more detail. It is assumed that the iterations converge regardless of the value of $p_{ij}^{(0)}$. The limit of convergence will be denoted $\Delta^{-1}_N(g_{ij})$. For N infinite, this limit equals the value of $\Delta^{-1}g(x,y)$.

$\Delta^{-1}_N(g_{ij})$ is an approximation of $\Delta^{-1}(g_{ij})$. So, $\Delta^{-1}$ and $\Delta^{-1}_N$ will represent the same operator in what follows because N is a constant depending only on the stage of numerization. $\Delta^{-1}$ is the operation by which the values of $p_{ij} = p_{ij}(\text{infinite})$ - the solution of equation [6] - are calculated. When the initial image is represented by $[g_{ij}]$ and if $[p_{ij}^{(k)}]$ corresponds to the potentials at the iteration k, then:

$$[p_{ij}^{(k+1)}] = [p_{ij}^{(k)}] * \begin{array}{c|ccc} 0 & 1/4 & 0 \\ \hline 1/4 & 0 & 1/4 \\ \hline 0 & 1/4 & 0 \end{array} + 1/4 [g_{ij}] \quad [8]$$

Consequently, the calculation of $\Delta^{-1}$ implies repeating a 3x3 convolution and adding a constant until it converges.

<table><caption>Figure 4 - the values of $g_{ij}$, $p_{ij}$ and the extraction of the crest line.</caption><tbody><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>255,00</td><td>255,00</td><td>255,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>255,00</td><td>255,00</td><td>255,00</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>12,96</td><td>23,09</td><td>29,20</td><td>31,84</td><td>32,50</td><td>31,84</td><td>29,20</td><td>23,09</td><td>12,95</td><td>0,00</td></tr><tr><td>28,73</td><td>50,20</td><td>61,88</td><td>65,67</td><td>66,30</td><td>65,67</td><td>61,88</td><td>50,19</td><td>28,73</td><td>0,00</td></tr><tr><td>51,76</td><td>87,09</td><td>102,44</td><td>102,65</td><td>101,36</td><td>102,64</td><td>102,44</td><td>87,08</td><td>51,76</td><td>0,00</td></tr><tr><td>91,22</td><td>143,94</td><td>158,15</td><td>141,11</td><td>133,85</td><td>141,11</td><td>158,14</td><td>143,93</td><td>91,22</td><td>0,00</td></tr><tr><td>169,18</td><td>239,29</td><td>245,09</td><td>169,78</td><td>151,82</td><td>169,78</td><td>245,09</td><td>239,28</td><td>169,18</td><td>0,00</td></tr><tr><td>91,22</td><td>143,93</td><td>158,14</td><td>141,11</td><td>133,85</td><td>141,10</td><td>158,14</td><td>143,93</td><td>91,22</td><td>0,00</td></tr><tr><td>51,76</td><td>87,08</td><td>102,44</td><td>102,64</td><td>101,36</td><td>102,64</td><td>102,43</td><td>87,08</td><td>51,75</td><td>0,00</td></tr><tr><td>28,73</td><td>50,19</td><td>61,87</td><td>65,66</td><td>66,29</td><td>65,66</td><td>61,87</td><td>50,19</td><td>28,72</td><td>0,00</td></tr><tr><td>12,95</td><td>23,09</td><td>29,20</td><td>31,84</td><td>32,49</td><td>31,84</td><td>29,20</td><td>23,09</td><td>12,95</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>25,25</td><td>81,14</td><td>86,95</td><td>11,64</td><td>10,71</td><td>11,64</td><td>86,95</td><td>81,14</td><td>25,25</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr><tr><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td><td>0,00</td></tr></tbody></table>
Next, the result of this operator applied to a small size image (9x9 pixels) is considered. This image reveals an edge having a discontinuity of three pixels represented by the first array of values in figure 4. This array gives the values of the $g_{ij}$. The second array shows the values of the potential $p_{ij}$. Note that basically a smoothing operation is performed. To get access to the pixels "joints", the line of highest potential is considered. Such lines are simply found by application of the Top Hat transformation, which is the difference between an image and the opening of its grey tone function (Serra, 1986). The result of the application of $\Delta^{-1}$ is represented by the third array of values. Figures 5.b and 5.c illustrate the "potential" in the case of a binary discontinuous edge with 512 x 512 pixels after 128 iterations and 25600 iterations, respectively.

<figure><img src="image_4.png" /><figcaption>Figure 5 a - b - c - visualization of the evolution of the $p_{ij}$.</figcaption></figure>

## EDGE THINNING AND SKELETONIZATION OF GREY TONE IMAGES

The previous notion of line of highest potential enables a treatment of the problem of discontinuous binarized edges. But this method also allows to handle grey tone images of edges, especially those obtained by Laplacian or Sobel filters. The thinning stage can be avoided as it is implicit in the notion of line of highest potential. So, $\Delta^{-1}$ can be seen as a "one-pass" thinning operator.

<figure><img src="image_5.png" /><figcaption>Figure 6 - edge thinning.</figcaption></figure>
Figure 6 demonstrates the effects on a binary image. The image on the right can be obtained after the calculation of $Δ^{-1}$, followed by the detection of lines of highest potential. Next, binarization is performed, whereupon masking is achieved with the image on the left.

Edge thinning is a special case of skeletonization. Figure 7 shows the skeleton obtained by following the process outlined, for simple binary images.

<figure><img src="image_6.png" /><figcaption>Figure 7 - skeletonization.</figcaption></figure>

The method of skeletonization described can be applied to a grey tone image. The skeleton is produced in grey levels, whereby its form will depend on grey levels of all pixels of the image. Figure 8 presents such a skeleton after binarization on the left, and superimposition onto the original image on the right.

<figure><img src="image_7.png" /><figcaption>Figure 8 - skeletonization of a grey tone image.</figcaption></figure>

## EDGE RECONNECTION

It will be demonstrated that some analogy exists with basic operators of mathematical morphology (erosion, dilation, ...). To that end, the result of $Δ^{-1}$ for an image made of a single point P can be observed more closely. Figure 9 explains this result in 3D. A surface is displayed with a vertical axis of symmetry (invariant
by rotation and translation of the Laplacian). The intersection with a horizontal plane is a circle with a radius between zero and infinity. For a plane $z=h$, ($0 < h < H$, $H$ being the vertical distance of the point of turning tangent to a plane of which the intersection circle has an infinitely large radius), parts of the surface placed above this horizontal plane are projected on a disk with a radius $R$. This disk can be considered the dilation of the point $P$ by a ball with a radius $R$. For accuracy reasons, the maximum radius of dilation is actually $R^\circ$ (which is a function of the grey tone of $P$, and of the number of iterations). Pixels beyond $R^\circ$ are not influenced by $P$. This influence zone is bounded by a black curve on figure 5 and 9.

<figure><img src="image_8.png" /><figcaption>Figure 9 - a grid representation for the $p_{ij}$.</figcaption></figure>

In doing so, $\Delta^{-1}$ remains the dilation operator of mathematical morphology, in spite of a completely different mathematical formulation. Nevertheless, seeing the analogy makes it worth continuing the comparison with the methods of image processing derived from mathematical morphology. A possible application of $\Delta^{-1}$, seen as a dilation operator, is the reconnection of edges shown in figure 10. It is to be noted that the larger the distance between segments (a: $3>2>1$), the narrower the bridge formed by the pixels "joints" (c: $3<2<1$).

## CONCLUSION

In spite of some inconveniencies (the uncertainty in the choice of the number of iterations and the "barbed" aspect of the result of a thinning down operation), the operator $\Delta^{-1}$ seems to be rich in possible applications. Noteworthy is the long computational time about 1 second per iteration on a grey tone image (512x512 pixels with 256 grey levels) with a RISC machine (15 MIPS). However, the very simple algorithm is similar to a convolution and should be treated in parallel (skeletonization in real time).

All characteristics seen in 2D are preserved in 3D. This generalization allows to realize a thinning down operation and a skeletonization of 3D grey tone images made of voxels. Obviously, calculation times will be much more significant in such cases.
<figure><img src="image_9.png" /><figcaption>a) original</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>b) dilated with $\Delta^{-1}$ (120 iterations)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>c) eroded with $\Delta^{-1}$ (30 iterations)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>d) thinning down with $\Delta^{-1}$ (60 iterations)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>e) classical dilation and erosion</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>f) morphological skeleton</figcaption></figure>

Figure 10 - a comparison between classical mathematical morphology and $\Delta^{-1}$.
## ACKNOWLEDGEMENT

The assistances of Professor Michel Jourlin and Professor Jean-Marie Becker are
deeply acknowledged.

## REFERENCES

- Babu R, Nevatia R. Linear feature extraction and description. Comp Graphics and Image Processing, 1980; 13: 257-69.

- Ciarlet P G. Introduction à l'analyse numérique matricielle et à l'optimisation. Paris: Masson, 1988.

- Ciarlet P G, Destuynder P. A justification of the two-dimensional linear plate model. J Mécanique 1979; 18: 315-44.

- Jain A K. Partial differential equations and finite difference method in image processing - Part I: Image representation. J Optimiz Theory and Appl 1977; 23: 65-91.

- Jain A K, Jain J R. Partial differential equations and finite difference method in image processing - Part II: Image restoration. IEEE Transaction on Automatic Control 1978; 23: 817-33.

- Serra J. Introduction to mathematical morphology. Computer graphics and image processing 1986; 35: 283-305.

Received: 1991-09-12

Accepted: 1992-07-30