Frédérique Robert ^(1,3) Guy Lefebvre ⁽²⁾

Guy Lefebvre received his Ph.D degree in Geotechnical Engineering from the University of Laval, Canada in 1970. He has been Professor at the Civil Engineering Department, Faculty of Applied Sciences, University of Sherbrooke since 1972 and actually he is involved in teaching, research and consulting on Geotechnical and Geoenvironmental Engineering. His main topics of interest are embankment dams, slope stability, wave velocity testings, soil decontamination, site remediation and valorisation of industrial wastes.

1. Introduction

Some thresholding methods have been implemented in order to obtain a particle map, in other words, a binary image on which white points have a high probability of belonging to particles and black points to the voids between them. But most automatic thresholding algorithms, Cocquerez and Philipp (1995), for example, do not give good results because they rely on the histogram shape and some points of medium grey-levels belong either to the particles or to the voids. That is why a unique threshold cannot discriminate the particles from the background. For example, thresholding by entropy maximization gives particle maps on which there are too many black points and these can belong either to voids or to particles (though the grey-leveled image has been inverted as void surface is clearly less important than the particle area). Another way of investigation is to use gradient computation in order to detect particle boundaries. There again, the results are not as good as expected. Actually, a lot of noisy points are detected in addition to the searched boundaries. These noisy points represent either the asperities on the rock surface or an edge separating two sides of a stone (see Figure 1).

The aim of this paper is to explain how distance maps can be used for particle disconnection or reconstruction on a grey-level image without prior knowledge of the particle shape. The following method of particle disconnection consists in using information from both results of thresholding and gradient computation. The resulting particle map is a combination of thresholding information (pointing out the particle inside) and gradient information (describing the boundaries). And then, a process which takes into account only points of the gradient image representing the particle boundaries has to be set up. Though this process using distance maps enables the selection of most points belonging to particle boundaries, it does not allow to find out every boundary point. That is why some of the particles remain connected, as a few points of their boundaries have not been detected. It is then necessary to rebuild particles. This rebuilding is also based on two distance maps, one giving a marker for each particle and the second its print. The distance map giving markers is unique although a different map is used for each particle print. The whole process, based on distance mapping has been implemented for an application in civil engineering. The purpose is to obtain the granulometry of a riprap (set of big stones covering an embankment) in order to determine the stability of the structure.

2. Segmentation of Riprap Images

2.1 Automatic Threshold Computation

Several methods of automatic thresholding have been used to find out the first particle map (Figure 5b), Labouré and Zeboudj (1987), Cocquerez and Philipp (1995), Sahoo et al (1988). But the one giving the best results is based on the computation of void area. As there are a lot of different grey-levels on considered particles, the grey-level is not a good criterion to make the difference between particles and voids. And then, for a subset of grey-levels neighbouring the threshold, points of same value belong either to particles or to voids. By computing the void area, the threshold is underestimated and that makes sure that black points belong to voids on the particle map, though white points have a greater probability of belonging to particles. The area ratio has been fixed to 1/8. This ratio is an empirical value, established by experiments from the density of voids in an elementary volume.

2.2 Watershed

A watershed algorithm has been implemented to find the boundaries (Vincent and Soille (1991)). But it failed to give the expected particle boundaries. Most of these boundaries are not really visible and this method does not enable them to be pointed out. Figure 2 shows the result of the watershed transformation and the thresholded image containing the found boundaries.

Figure 2. a. Watershed Transformation On A Riprap Image

Figure 2 b. Thresholded Watershed Image

2.3 Other Methods

Some other methods such as a morphological opening can be used but they do not give results good enough to be compared.

3. Distance Map

A certain binary image represents a set of objects on a background. Let us assume that objects are white and the background black. A distance map is computed from this binary image. It gives, for each point of the background, its Euclidean distance to the object set, in other words, its Euclidean distance to the nearest object (Figure 1).

A reverse distance map can also be computed, giving for each white point its distance to the black set (Figure 3). A point is assimilated to a square or to the center of this square; for distance computation, the square centers are only considered. Then, the distance values are given in Table 1.

The main problem to solve when computing a distance map is to minimize the operation number because it can be very long to evaluate distances. The primary algorithm is to compute for a point M of the background its distance to each point of the objects and then, to find out the one giving the minimal value. This value is the distance value at point M. This method requires a lot of operations : two square evaluations, an addition and a square root computation for each point of the objects, per point of the background. It can be improved by taking into account only points belonging to object boundaries but it remains very long.

Some authors, Borgefors (1986), Danielsson (1980), have proposed other solutions enabling the computation of distance map by incremental algorithms. But some of these algorithms do not give an accurate evaluation of the Euclidean distance. The algorithm proposed by Danielsson (1980) gives good values in most cases. Tables 2 and 3 give X and Y coordinates for the example of Figure 4. Then, if a point M has (x, y) as relative coordinates, its distance value d on the distance map is given by the classical formula :

This algorithm requires a few square root computations per point of the image, whatever the number of the object points is. It will be used to compute the following distance maps.

	4				1	1	1					1	1	1
	3	3	2	1	0	0	0	1	2	2	1	0	0	0	1	2	3
	2	2	2	1	0	0	0	1	2	2	1	0	0	0	1	2	3
	1	1	1		1	1	1			2	1	0	0	0	1	2	3
1	0	0	0	1	2		1	1			1	0	0	0	1	2	3
	1	1	1		2	1	0	0	1	1	0	0	0	0	1	2	3
	2	2	2		2	1	0	0	1	1	0	0	1	1
	3	3	3	3	2	1	0	0	1		1	1		2
	4	4	4				1	1			2	2

1	0	3	2	1	0	0	0	1	2	2	1	0	0	0	1	2	3
1	0	0	2	1	0	0	0	1	2	2	1	0	0	0	1	2	3
1	0	0	0	1	0	0	0	1	2	2	1	0	0	0	1	2	3
1	0	0	0	1	0	0	0	1	2	2	1	0	0	0	1	2	3
1	0	0	0	1	2	1	0	0	1	1	1	0	0	0	1	2	3
1	0	0	0	1	2	1	0	0	1	1	0	0	0	0	1	2	3
1	0	0	0	1	2	1	0	0	1	1	0	0	1	1	1	2	3
1	0	0	0	3	2	1	0	0	1	1	0	0	1	0	1	2	3
1	0	0	0	3	2	1	0	0	1	2	0	0	1	2	1	2	3

4	4	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1
3	3	3	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
2	2	2	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0
1	1	1	1	1	1	1	1	1	1	0	0	0	0	0	0	0	0
0	0	0	0	0	0	1	1	1	1	1	0	0	0	0	0	0	0
1	1	1	1	1	0	1	0	0	0	0	0	0	0	0	0	0	0
2	2	2	2	2	2	1	0	0	0	0	0	0	1	1	1	1	1
3	3	3	3	0	0	0	0	0	0	1	1	1	1	2	2	2	2
4	4	4	4	1	1	1	1	1	1	1	2	2	2	2	3	3	3

4. Particle Disconnection

The main steps of the particle disconnection process are summarized in Figure 5. First, two different particle maps are computed, one by thresholding and the second one by combining thresholding and gradient computation.

A particle map is a binary image on which white points are supposed to belong to objects.

The thresholding process yields a particle map which has no noisy point but most particles remain connected. This first particle map becomes a good basis for the process. Nevertheless, such information is also included in the gradient image because it describes every particle boundary but some points represent edges or asperities of particles and they must not be taken into account.

Figure 5 shows the two particle maps (b and d) extracted from the initial image (a). The gradient image (c) is filtered by a high-pass filter and the remaining points (having a high contrast on the initial image) are taken as points of voids. Edges and boundaries detected by gradient computation are associated with the thresholded image. The superposition of the two images allows to get a particle map : the mixed image (d).

This is a binary image, white points representing the inside of particles and black points the voids between them.

Then, distance maps are computed from the binary images and their combination gives the final particle map (e).

4.1 Gradient Computation

As boundaries lie in high contrast variation regions of the image, computing a gradient image is a good solution. As well, the interesting values on such an image are the non-zero values. But particle surfaces are not very smooth, and this implies that small shadows will be detected as contrast variation, because of the presence of small holes or peaks. Furthermore, particle grey-level is not homogeneous over the whole surface and this also creates non-zero values on the gradient image.

That is why the gradient image needs to be filtered. Obviously, it can be high-pass filtered in order to retain only the greatest values. But, because of contrast grey-level or non-contrast grey-level, some of the particles may remain connected.

For example, if two particle surfaces belong approximately to the same plane and if both grey-levels are very similar on the interface, the gradient values, between the two particles, will be low enough not to be kept. Thus the common boundary will be lost, and the particles still connected.

The gradient image is computed by convolutions. Table 4 shows the masks used.

Table 4 a and b. Mask A and Mask B

1	-1	-1	-1	-1		1	1	1	1	1
1	4	-9	-4	1		-1	4	9	4	1
1	9	16	-9	1		-1	-9	16	9	1
1	4	9	4	1		-1	-4	-9	4	1
1	1	1	1	1		-1	-1	-1	-1	1

Each mask is convoluted to the initial image so that the values of the two convoluted images are considered as gradient coordinates.

The module of each gradient vector is then computed and the boundaries appear on the new image. But there is a lot of noise due to particle defaults or grey-level inhomogeneity.

This gradient image is very interesting because it contains all the boundaries. So the required information is computed this way. The next step of the algorithm is to separate edges from noise.

Furthermore, edges are highlighted and completely detected. So, noise elimination is the true issue, to get a good boundary image. Mask values are chosen so that they detect all grey- level variations in the initial image.

A conventional method of thresholding does not yield a convenient distinction between good values (representing parts of boundaries) and bad values. Furthermore, particle edges may be detected too, if a particle has two visible sides, and the particle is then cut into two or more parts.

Applying most of the usual methods of thresholding (histogram thresholding, entropy maximisation, and so on) shows that it is impossible to compute a global threshold, which can determine over the whole image whether a gradient value must be taken into account or not. The gradient image is then high-pass filtered before being filtered more accurately.

4.2 Distance Maps for Gradient Filtering

The contrast gradient computed in this way is very noisy (asperities on particle surfaces). Noise can be eliminated by thresholding the gradient image using windows of decreasing sizes, Robert and Lefebvre (1995). This method presents two drawbacks : firstly, there is some residual noise, and secondly, some true boundary parts may disappear. However, it is efficient enough to partially disconnect the particles. The solution we consider, using distance mapping is based on two different binary images. The first is the thresholded image on that black points have a probability , whereas for white points, is lower. Actually, we assume that boundary points belong to the particles in order to get disconnected particles on the final map (see Figure 5e), in other words we are looking for the inner boundaries. The second one is the mixed image on which for a white point , whereas for a black point is lower.

Let BP be the set of black points of the mixed image where is lowest and let WP be the set of white points of the thresholded image where is lowest.

Let BD be the point set of particle boundaries on the thresholded image.

Let GP be the set of points remaining on the gradient image after high-pass filtering, then :

GP = BP U WP U BD

BP (resp.WP) is the point set where gradient values correspond to asperities or edges on the particle surfaces (resp. to particle boundaries)(Figure 6).

The distance maps (see Figures 7 and 8) are computed, using Danielsson’s.

Let DP be the set of points where the value is equal to 1 on the gradient distance map. Considering the mixed image and B the ball of radius 1, V is the point set representing voids and :

DP = (V Å B) \ V

where :

is the Minkowski addition (Serra, 1982)

The filtering method consists in sharing DP in three disjoined subsets DP1, DP2 and DP3 (see Figure 9) (corresponding to BD, BP and WP), such that :

- For each point of DP1, both values on the initial distance map and the gradient distance map are equal to 1. DP1 is the set of points neighbouring the boundaries marking off particles on the thresholded image.

- For each point of DP2 and DP3, the value on the gradient distance map is higher than 1, but points of DP2 can be related to a point of DP1, whereas points of DP3 cannot.

- A point P2 of DP2 is related to a point P1 of DP1 if P1 itself or another point P’2 of DP2 belongs to the ball centered at DP2.

- Points belonging to DP2 are neighbouring searched boundaries. Points of DP3 are near asperities or edges.

- DP2 represents prolongations of the boundaries.

Note that particle disconnection is equivalent to gradient filtering, because the points belonging to the particles boundaries are particular points of the gradient image because of their values. That is why distance maps are used to take into account the particle topology, and then, the gradient image is filtered. This process can use any kind of discrete distance computation but Danielsson’s algorithm gives better results because particles are not regular polygons.

5. Particle Rebuilding

When particles are connected, it is impossible to compute their Féret diameter (maximum diameter) without rebuilding them. The first step of the rebuilding process is to find out some markers giving a good approximation to the particle shape and then, to compute a dilation of this approximation in order to obtain a shadow (projection of the particle in the acquisition plane) of the particle. We will do it with distance maps.

The distance map giving particle markers is unique. It is computed from the final particle map (Figure 5e). Markers are then extracted from this map by thresholding. Each particle has a specific thresholding . First, a pixel (i,j) where the absolute maximum value Max is reached on the distance map. Max is the radius of the inscribed disk in the largest particle and it gives an estimation of the sieving diameter as particles are almost convex. Then a square window, centered in (i,j) and whose side is equal to 2.Max, is considered. Inside this window, the distance map is thresholded according to a threshold value equal to . The remaining points (whose distance values are larger than ) describe the marker M of the biggest particle. The marker is written in a new binary image in which all points belong to the background but points belonging to the marker (object). Then, a new distance map (Figure 10) is computed, out of which the dilated set will be determined.

This dilation is defined as follows :

where : P is the shadow of the particle (the opening of the particle by B)

M is the marker (the eroded particle by B)

B is the disk of radius

is the Minkowski addition (Serra (1982))

The value allows to obtain a good marker which does not take into account the remaining connected shapes. The marker has an inscribed disk of radius value  ; dilating it by B rebuilds a particle with an inscribed disk of radius value Max (as assumed) (Figure 11). An opening of the biggest particle is done in the same way by disk B. Figure 12 shows the biggest rebuilt particle.

Figure 12. Particle Rebuilding

Then, the same process can be used for other particles. The maximum value research on the distance map is done by decreasing progressively this value. But it is necessary to check if the new point N (center of the new marker) belongs or not to a particle already considered. First, the distances NNp, between N and the previous considered points Np, are compared to the sums d(N)+d(Np) (d(M) is the value at the point M on the distance map). If NNp<d(N)+d(Np) then N is rejected else it can be the center of a new marker. That is the distance condition of acceptance (Figure 13).

Second, the assumption that particles are almost convex is used : lines are drawn between N and the Nps and, if one of these is totally included in the particle set, then N is rejected. That is the convexity condition of acceptance (Figure 14). If N can be considered, the new marker M(N) is found by considering the window centered at N and whose width and height is 2.d(N), to threshold the distance map. The threshold value is . Then, a new binary image is considered, only containing the new marker as an object. The associated distance map is computed to determine the shadow of the particle.

This process is applied to all particles. Its results are good, although it is time consuming. Some problems may remain on very elongated particles because a good deal of information about thin endings is lost by opening. Figure 15 shows the result of the reconstruction process on the whole image.

Figure 15. Reconstruction of Every Particle

6. Application

An image processing project has been launched by the University of Sherbrooke and funded by Hydro Quebec Company, (Lefebvre et al (1989), Robert (1994)). The aim was to obtain a granulometry of the riprap of existing embankment dams by processing video recordings. These video recordings are transformed into an image set which is representative of the riprap surface. Each image is then processed in order to obtain a « good » particle map, where rock diameters can now be estimated. A distance map is computed, giving for each white point (belonging to the rocks) its distance to the voids using the previous algorithm. Sieving diameters correspond to local maxima on this map assuming that rocks have almost no convex irregularities, and that the inner diameter (measured from the distance map) is approximately equal to the outer diameter (experimentally measured by a siever).

As some particle pairs may remain connected, these local maxima may not be unique for some connected sets. But, assuming that particles have almost no convex irregularities, we can check whether this property holds or not. This is done in the same way as the convexity condition rebuilding was checked. Every point M (center of a marker) is then a center of a particle and d(M) is its sieving diameter value. That gives a first set of diameter values.

In order to make the granulometry evaluation more accurate, a second set of diameter values is computed by determining the Féret diameter of each particle, after rebuilding. In order to do that, several different directions are chosen, given by their angle q. For each of them, the coordinate system is rotated, in order to find the diameter in the given direction just by computing the minimum xmin(q) and the maximum xmax(q) of the X-coordinate. The Féret diameter is the maximum value of the differences xmax(q)-xmin(q).

When doing hand experiments, three diameters are measured on each stone : the maximum diameter d₁, the sieving diameter d₂ and the minimum diameter d₃. Experiments show that the sieving and the minimal diameter are almost identical whereas the maximum diameter is very different. The mass of a stone is evaluated according to the following formula :

where :

t is the ratio between the average volume of stones and the parallelepiped having edge lengths d₁, d₂, d₃

m is the specific mass of the stone

According to this formula, the average diameter is given by :

By image processing, only two diameters are computed. And the average diameter can be estimated by :

where : d_F is the Féret diameter

d_S is the sieving diameter

The mass is then estimated by :

Granulometry curves (experimental and image processing ones) can be drawn, giving for each diameter value d the percentage of rock mass, representing stones of average diameter smaller than d.

7. Conclusion

The existing methods for particles disconnection are not efficient enough when particles have a non homogeneous surface or do not belong to the same plane (that induces shadows and grey-levels variations). A new method based on distance mapping has been studied and implemented. This process, less time-consuming, gives better results than the previous ones. It is sufficient to assume that a good approximation of the particle boundaries has been obtained by thresholding. Then, a detection of prolongations of these boundaries from the gradient image is computed in order to complete the disconnection. Such a method can be applied to every set of particles (sharing a common boundary), provided that there is sufficient contrast between them. That implies that particles must not perfectly fit into each other. But they can cover each other partially. An application in civil engineering has been presented . The considered images represent an embankment made of a stone piling. The stone surfaces are not homogeneous and have asperities and colour differences. This method for disconnection is well-adapted to such images and yields good results. The rock particles are enough disconnected to enable sieving diameters computation. This process has the advantage that it does not depend on the particle shape. It is not restricted to spherical or ellipsoidal particles. And it can be used for any kind of shape provided that particles have almost no convex irregularities. Furthermore, particle disconnection is equivalent to gradient filtering.

But, as this disconnection has not been totally completed, it is necessary to compute a particle rebuilding to enable the maximum Féret diameter computation. Using markers extracted from the distance map associated with the final particle map yields a good approximation of the particle shapes. Furthermore, the dilation obtained from distance maps, when rebuilding, is more accurate than dilations with a disk of unit radius because of the errors of approximation for disk drawing on a grid, and it is faster than dilation with a large discrete radius disk. Using distance maps reduce these errors, particularly for dilation with large disks.

References

Borgefors, G., Distance Transformations in Digital Images, CVGIP. 34, No.3, 1986, pp. 344-371.

Cocquerez, J.P. and Philipp, S., Analyse d’images: filtrage et segmentation, Éditions Masson, France, 1995.

Danielsson, P. E., Euclidean Distance Mapping, Computer Graphics and Image Processing, 14, 1980, pp. 227-248.

Labouré, M. J. and Zeboudj, R., Seuillage automatique. Publication No.11, Department of Mathematics, University of Saint-Etienne, France, 1987.

Lefebvre, G., Rohan, K., Ben Belfadhel, M. and Dascal, O., Field Performance and Analysis Steep Riprap, ASCE, Journal of Geotechnical Engineering, Vol. 118, No. 9, 1989.

Riley, J. Sedim. Petrol. V 11, 2, 1941.

Robert, F. and Lefebvre, G., Distance Mapping for Image Filtering, Proceedings of International Conference for Stereology, Copenhagen, 1995.

Sahoo et al, A Survey of Thresholding Techniques, CVGIP, Vol. 41, 1988.

Serra, J., Image Analysis and Mathematical Morphology, AC. Press, New York, 1982.

Vincent, L. and Soille, P., Watersheds in Digital Spaces: An Efficient Algorithm Based On Immersion Simulations, IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 13, No. 6, 1991, pp. 583-598.