Abstract

This paper deals with reconstructions of series of partially focussed images. Some of these methods are known in case of images which were acquired in small field of view (by confocal microscope or CCD camera, e.g.). In this case, recorded images do not differ in any geometrical transformation from each other. In case of larger samples (oversized for microscope or CCD camera), it is necessary to use wider viewing field (standard cameras, e.g.), and taken images primarily differ in scaling but may also differ in shifting and rotation too. These images cannot be used for reconstruction directly; they must be registered; that is, we must determine all transformations which the images differ and eliminate their effects. There are several ways to do this. This paper deals with the registration based on phase correlation. After this registration, it is necessary to identify the sharp parts and to compose a and model. Present methods are very sensitive to noise and their results are not satisfactory in many cases. We introduce a new method for reconstruction which is significantly better.

1. Introduction

The three-dimensional reconstruction of general surfaces plays an important role in many branches; for example, the morphological analysis of fracture surfaces reveals information on mechanical properties of natural or construction materials.

There are more techniques capable of producing digital three-dimensional (3D) replicas of solid surfaces. Mechanical engineers can use contacting electronic profilometers to determine digital two-dimensional (2D) profiles that can be combined into 3D surface profiles; see [1, 2], for example. The contacting mode of atomic force microscopes actually belongs to this mechanical category [3]. Besides the mechanical tools, there exist different optical devices [2], light section microscopy [4, 5], coherence scanning interferometry [6], speckle metrology [7], stereo projection [8], photogrammetry [9], and various kinds of light profilometry [10], to mention some of them.

3D laser scanning techniques are the next possibilities of how to obtain the data. These techniques have also been tested in some rock engineering projects, such as 3D digital fracture mapping [1113].

However, these devices are not of universal use. Each of them has technical limits [14, 15]. For example, very rough surfaces can hardly be measured by atomic force microscopes, which work in the nanoregions. On the other hand, plane surfaces with microscopically small irregularities may be measured, for example, by the microscopic sectional technique within the so-called confocal microscopes [1619]. However, confocal microscope is often not suitable for technical purposes due to the small size of the visual field (maximal visual field is approximately 2 cm [5, 2022]).

In this paper, how to perform the 3D reconstruction of larger surfaces using a standard camera will be shown.

2. Materials and Methods

2.1. Equipment

In technical practice the confocal microscope is used as a standard tool for imaging microscopic three-dimensional surfaces. The depth of the optical field of the microscope is very small and its advanced hardware is capable of removing nonsharp points from the images. The points of the object close to the focal plane are visible as sharp points. The parts lying above or beneath the focal plane are invisible and represented as black regions. To create a 2D and 3D reconstruction, it is necessary to obtain a series of images of the same object, each of them with different focusing, and each point of the object focussed in one of the images (in the ideal case). The sharp parts are identified and composed of a 2D and 3D model.

However, the confocal microscope is often not suitable for technical purposes due to the small size of the visual field (maximal visual field is approximately 2 cm [5, 2022]). Nevertheless, the same principle may be used even in the case of a CCD camera, classical microscope, or camera. For nondestructive scanning, the camera must be mounted on a stand which enables a movement in the direction approximately orthogonal to the surface with a controlled step.

The difference between confocal microscope and standard camera concerns nonsharp regions that are displayed by classic camera, whereas they are missing in the case of a confocal microscope. However, the sharp and blurred areas can be detected by software. The next dissimilarity lies in a central projection which caused different scaling of partial images in the image series (see Figure 1). Different image scaling (including possibly shift and rotation) may be also quantified and corrected. This problem has not been solved in available sources.

2.2. Present Methods

Probably the first attempts to carry out nonconfocal reconstructions come from the seventies and eighties of the last century [2328]. Blurred areas detectors (or so-called focusing criteria) are based on various principles. Present methods work in three steps.

The First Step. Three-dimensional matrix is stated. Its element determines statistical range, variance, or standard Fourier transform of certain neighborhood of the pixel in th image in the series of images [29].

The Second Step. Maxima in the columns are found. Height is assigned to the pixel if and only if the maximum is detected on the th image. In this way, a stair-approximation is obtained.

The Third Step. The stair-approximation from the second step is obviously refined by various interpolation of detector values adjacent with these maxima, that is, values ; : inverse proportionality [29, 30], parabolic fits [3133], or Gaussian fits [34].

In [35], it is stated: “We have verified that there are only small differences between the three-point Gaussian and the three-point parabolic approximations.” However, this statement does not indicate an accuracy of the reconstruction. As shown below, this fact can only mean that both reconstructions are roughly equally bad.

A noise is the fundamental problem of the second and the third step. This problem was solved by several ways: varying size of computational pixel windows [33, 36], averaging multiple snapshots taken at each vertical position [26], or input data averaging [37]. However, these methods are able to decrease a noise but the noise cannot be zeroed. Therefore, the maxima computed in the second step need not indicate the height correctly; moreover, any interpolation (the third step) is not useable for noise data from point of view of numerical mathematics.

Except for these inaccuracies, methods cited above must assume the same scale of all images in the series because they are not able to detect image transformations. In the following sections, we will introduce the reconstruction method which is able to process a series of noisy images acquired using central projection with variable center, that is, in various scaling (various shift and rotation eventually).

2.3. The Fourier Transform and Phase Correlation
2.3.1. Standard Fourier Transform and Inverse Transform

Standard (continuous) Fourier transform of function is function (if this integral exists and is finite).

Standard (continuous) Fourier transform of function is function (if this integral exists and is finite).

Inverse (continuous) Fourier transform of function is function (if this integral exists and is finite).

Inverse (continuous) Fourier transform of function is function (if this integral exists and is finite).

Standard discrete Fourier transform of function is function Function is also called the Fourier spectrum of function . It is possible to obtain the function from its Fourier spectrum using inverse discrete Fourier transform (see [38] for proof). Function is called amplitude spectrum of .

2.3.2. δ-Distribution

One-dimensional -distribution is a limit of a sequence of functions ; such that

Two-dimensional -distribution is a limit of a sequence of functions ; such that

Example 1. One of the well-known examples (which we demonstrate in for simplicity) is the series of expanding rectangular signals with constant unitary intensity on ; and zeroed elsewhere. Inverse Fourier transform givesBecause for each (see [36] for proof), condition (a) is fulfilled. Because for each , there is It means that also condition (b) holds and limit is the (one-dimensional) -distribution. Expanded unitary signal and its Fourier transform are illustrated in Figure 2.

2.3.3. Phase Correlation

For image processing, it is necessary to transform the images so that the studied structures are at the same position in all the images. This is the task of image registration, to find the transformation. In some applications we assume that images were shifted only; in others we allow shift, rotation, and scale change (i.e., similarity), general linear transformation, or even general transformations. The methods used for registration depend on the expected transformation and on the structures in the image. Some methods use corresponding structures or points in the images and then find a global transformation using the measurements of positions of the structures or points [3941]. These methods require these structures to be clearly visible. Other methods are based on correlation and work with the image as a whole. The phase correlation proved to be a powerful tool (not only) for registration of particular focussed images. For functions ; , it is defined as and its modification aswhere strip means complex conjugation and is a bounded real function such that and are arbitrary constants. It can be proved that for real functions ; the phase-correlation function is real [42]. This is of great value, since it enables us to search for extremes of the phase-correlation function.

2.4. Image Transformation

Identical Images. Let be the infinity periodic expansion of an image. Denote as its value in the point , . It is evident that the value of the phase correlation of the with itself is According to Example 1 in Section 2.3.2, one hasIt means that the inverse Fourier transform of the convolution of two identical images is the two-dimensional -distribution .

Shifted Images. If two functions are shifted in arguments, that is, , their Fourier transforms are shifted in phase; that is,and their phase-correlation function is the -distribution shifted in arguments by the opposite shift vectorThis is the main idea of phase correlation. The task to find a shift between two images is converted by the phase correlation to the task of finding the only nonzero point in a matrix. If the images are not identical (up to a shift), that is, if the images are not ideal, the phase-correlation function is more complicated, but it still has a global maximum at the coordinates corresponding to the shift vector.

Rotated Images. The phase-correlation function can be also used for estimation of image rotation and rescale. Let be function rotated and shifted in arguments; that is,Their Fourier spectra ; and amplitude spectra ; are related as follows:The shift results in a phase shift and the spectra are rotated in the same way as the original functions. A crucial step here is transformation of the amplitude spectra into the polar coordinate system to obtain functions such that . The rotation around an unknown center of rotation was transformed to a shift. This shift is estimated with the standard phase correlation (see the previous paragraph) after rotating back by the measured angle; the shift is measured with another computation of the phase correlation.

Scaled Images. Let be function rotated, shifted, and scaled in arguments; that is,Their Fourier spectra and amplitude spectra are related as follows:The shift results in a phase shift, and the spectra are rotated in the same way as the originalfunctions and scaled with a reciprocal factor. A crucial step here is transformation of theamplitude spectra into the logarithmic-polar coordinate systemto obtain such that .

Both rotation and scale change were transformed to a shift. The unknown angle and unknown factor can be estimated by means of the phase correlation applied on the amplitude spectra in the logarithmic-polar coordinate system ; . After rotating function back by the estimated angle and scaling by factor , the shift vector .

2.5. Proposed Method of Surface Approximation
2.5.1. Preprocessing: Image Registration

Before profile calculation, it is necessary to identify all geometric transformation in the image series and to eliminate them. The images are analyzed as geometrically similar. Generally, similarity is a combination of rotation, scale change, shift, and axial symmetry. Axial symmetry is not possible in our case.

The discrete Fourier transform described above is used for the image registration. It either works with periodic functions or makes them periodic. In general case, an image has not the same values on the edges and by periodizing an image, we obtain a function with jumps at the edges of the original image. These jumps are often the most contrast structures in the function and may lead to incorrect registration. Therefore, it is necessary to remove such edges from the image used for the shift estimation, to smooth them out. This is done by multiplying the image by a suitable function, a so-called window function. Such function must be zero or almost zero at the image edges and one on a large part of the image. We can primarily use the Gaussian window function and the Hanning window function. Let be a given number and setLet be the distance of point from set ; that is,Functionsare called rectangular or circular Gaussian window function. Functionsare called rectangular or circular Hanning window function.

Let be the image series to be registered; image was acquired by means of the biggest angle of view. This image will be not transformed or (formally) it will be transformed by identity mapping to the image . Now we must find transform to obtain image which differs from in focussed and blurred parts only. In the same way, transforms must be found.

After the multiplying both images ; by the chosen window function, we determine the rotation angle and scale factor by means of the method described in Section 2.4. Then, image is rotated by angle and scaled by factor to compensate the rotation and scale change found by the phase correlation, creating image . Only shift and different focussed and blurred parts remain between image and image . Now we can apply phase correlation to find the shift and image is shifted by vector to compensate the shift, creating image which differs from in focussed and blurred parts only.

2.5.2. The First Step: Matrix of Sharpness Detectors

As was said in Section 2.2, statistical range, variance, or standard Fourier transform of certain neighborhood of the pixel may be used as the sharpness detectors (focusing criteria). The neighborhood of the pixel is a square of pixels where is obviously ten to twenty. In case of the standard Fourier transform, algorithm FFT is used. It requires the square with the side ; that is, or is used in this case. However, standard Fourier transform suffers from jumps at the edges of the square and it is necessary to use a suitable window function like in Section 2.5.1. Therefore, cosine Fourier transform is preferable. This transform is obtained from the standard Fourier transform which is applied to even extension of the neighborhood to be processed. This extension is illustrated in Figure 3. Due to this extension, there are no jumps at the edges. Sine frequencies in (5) are zeroed and application of any window function is not necessary.

However, low frequencies in amplitude spectrum detect blurred parts of image and very height frequencies are given by noise. Therefore, suitable weight must be assigned to each frequency in sharpness detector calculation. In our software, the following detectors may be used:where is the cosine spectrum of the pixel neighborhood and is volume of the annulus with the center . Elements ; ; which are summed in (26); (27); (28) are illustrated as pixel values in Figure 4. Size of detectors really used in software is or pixels. In Figure 4, higher resolution is used for better illustration.

2.5.3. The Second Step: Profile Heights Calculation

The main imperfection of current methods is the hypothesis that the profile height in given pixel is exactly determined by values of chosen sharpness detector. This hypothesis implies that these values can be interpolated. However, this conclusion is quite false. We have available the series ; for assessment of the height of pixel . This series is not deterministic but it is a random variable. It cannot be interpolated, and it must be processed by statistical method. One of the possibilities is a regression analysis but it would be very complicated. Direct calculation of the expected value is much easier.

For each pixel , we can construct theoretically infinitely many probability distribution functions using different exponents applied to series members :Expected values of random variables given by these probability functions estimate the height of surface in its pixel :

3. Results and Discussion

3.1. The Preprocessing: Image Registration

In Figure 5, we can see the photos from Figure 1 but they are displayed in so-called supplementary pseudocolors. If these images were completely identical, then arithmetic mean of the blue-green image on Figure 5(a) and orange image on Figure 5(b) would be “perfectly grey.” Arithmetic mean of these images is constructed in Figure 6(a). It is clear that components of this mean are very different, values of the orange image are bigger in yellow parts of the mean, and values of the blue-green image are bigger in blue-violet parts of the mean.

In Figure 6(b), the same construction after the registration is realized. Very low color saturation of arithmetic mean confirms very good conformity. Of course, arithmetic mean cannot be “perfectly grey” in our case because components of the mean differs in blurred and sharp parts.

In Table 1, indicated and applied transforms in separated images of the blue marble are summarized. All transformations have been detected with subpixel precision; they are listed with a precision of one thousandths pixels. It is obvious that the scaling plays most important role; however, nor shifts are negligible. Rotation angle between the first and the last image is over five arcminutes; it means approximately one pixel on the image periphery (used data resolution is pixels). This transformation is marginal in our case.

3.2. The First Step: Sharpness Detectors

As is clear from (26), (27), and (28) and Figure 4, detector increased high frequencies (right bottom corner of squares in Figure 4) too much; it means that it is very noise sensitive. The disadvantages of the detector are too sharp cuts which remove low and height frequencies.

Differences between the sharpness detectors stated in Section 2.5.2 are best demonstrated by their maxima. The graphical representation of detector (i.e., corresponding stair-approximation of the surface) is shown in Figure 7.

Differences between representations ; ; and are visually quite miniscule (therefore, and are not illustrated); however, differences exist. For their quantification, Root Mean Square Error is as follows:Average Deviation is as follows:Pearson Correlation Coefficient is as follows:and Difference of surface Information Entropy is as follows:where ; are width and height of surface domain in pixels; ; are height of first (second) surface in the pixel ; ; are the average height of the first (second) surface; andis the Information Entropy of the surface .

Values of these characteristics are summarized in Table 2. Note that it is ; for a pair of identical surfaces.

Similar summary is made in Table 3 for the present reconstruction methods: parabolic interpolation of the stair-approximations , , and . In Table 4, there are summarized , , , and for parabolic and Gaussian interpolation of the same stair-approximation and proposed method of the expected values of columns ; .

As was already stated in [35], only small differences are between parabolic and Gaussian interpolation. Really, values of , , and for parabolic and Gaussian interpolation are approximately thousand times smaller than for interpolation and proposed method. Correlation between interpolations is stated even as one in Table 4. Therefore, it might seem that interpolations are much better than expected value estimation but the opposite is true. Small differences between these methods mean only that these reconstructions are roughly equally bad as shown in Figures 8 and 9. Results illustrated in Figures 8 and 9 have been obtained by combination of present and proposed method (they are not realized by present methods only due to different scaling of the input data). Visualization of the parabolic interpolation of is shown in Figure 8; Gaussian interpolation is illustrated in Figure 9. Practically the same and large noise is evident in these reconstructions.

Low-pass filters are commonly used to reduce the noise, but these filters are not ableto differentiate whether high-frequency information is caused by noise or by small details of useful information. Therefore, the loss of useful information is inevitable as shown in Figure 10.

In Figure 11, there is illustrated significantly better result of proposed method applied to the same matrix . There is no visible noise at the output, and high-frequency useful information (small surface details) is retained.

4. Conclusion

Reconstructions of three-dimensional object surfaces are an important task in many branches of research. Even though the standard method of imaging object surfaces is the use of confocal microscopes or laser scanners, there exist sophisticated mathematical methods that are able to process images acquired by classic cameras and to construct 3D reconstruction similar to these from confocal microscopes or laser scanners. This enables us to obtain similar results with substantially less expensive equipment.

Conflicts of Interest

The author declares that they have no conflicts of interest.

Acknowledgments

This work was supported by the Ministry of Education of the Czech Republic, Project LO1202 with financial subsidy, under the NPU I Programme. Thanks are due to Pavel Štarha for providing data.