3D Reconstruction of the Surface Using a Standard Camera

Copyright © 2017 Dalibor Martišek. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. This paper deals with 3D 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 2D and 3Dmodel. Present methods are very sensitive to noise and their results are not satisfactory in many cases. We introduce a new method for 3D reconstruction which is significantly better.


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 3D data.These techniques have also been tested in some rock engineering projects, such as 3D digital fracture mapping [11][12][13].
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 [16][17][18][19].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,[20][21][22]).
In this paper, how to perform the 3D reconstruction of larger surfaces using a standard camera will be shown.

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,[20][21][22]).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.

Present Methods.
Probably the first attempts to carry out nonconfocal reconstructions come from the seventies and eighties of the last century [23][24][25][26][27][28].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 max  {  } in the columns { 1 ;  2 ; . . .;   } 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.
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 3D 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).

The Fourier Transform and Phase Correlation
(if this integral exists and is finite).Standard (continuous) Fourier transform of function (; ) : R 2 → C is function (; )  −(+) d d (2) (if this integral exists and is finite).

𝛿-Distribution.
One-dimensional -distribution () is a limit of a sequence of functions   ();  ∈ N such that Two-dimensional -distribution (; ) is a limit of a sequence of functions   (; );  ∈ N such that  Because It means that also condition (b) holds and limit () = lim →∞   () is the (one-dimensional) -distribution.Expanded unitary signal and its Fourier transform are illustrated in Figure 2.

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 [39][40][41].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  1 ;  2 , it is defined as and its modification as where strip means complex conjugation and (; ) is a bounded real function such that (; ) = (−; −) and ;  > 0 are arbitrary constants.It can be proved that for real functions  1 ;  2 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.

Image Transformation
Identical Images.Let  be the infinity periodic expansion of an image.Denote  +  as its value (; ) in the point (; ),  +  ̸ = 0.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 has It 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,  2 (; ) =  1 ( −  0 ;  −  0 ), their Fourier transforms are shifted in phase; that is, and their phase-correlation function is the -distribution shifted in arguments by the opposite shift vector This 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  2 be function  1 rotated and shifted in arguments; that is, Their Fourier spectra  1 ;  2 and amplitude spectra  1 ;  2 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 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 system exp  = √  2 +  2 ;  = exp  cos ;  = exp  sin  (21) to obtain 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  1 1 ;  1 2 .After rotating function  2 back by the estimated angle  and scaling by factor , the shift vector ( 0 ;  0 ).

Proposed Method of Surface Approximation
2.5.1.Preprocessing: Image Registration.Before 3D 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  ∈ R + be a given number and set Let (; S) be the distance of point  = (; ) from set S; that is, Functions  GR (; ) =  − 2 (;R)/ 2 or  GC (; ) =  − 2 (;C)/ 2 (24) are called rectangular or circular Gaussian window function.Functions are called rectangular or circular Hanning window function.
Let { 1 ;  2 ; . . .;   } be the image series to be registered; image  1 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  * 1 .Now we must find transform  2 →  * 1 to obtain image  * 2 which differs from  * 1 in focussed and blurred parts only.In the same way, transforms  3 →  * 2 ; . . .;   →  * −1 ; . . .;   →  * −1 must be found.
After the multiplying both images   ;  * −1 by the chosen window function, we determine the rotation angle   and

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  = 2  ; that is,  = 8 or  = 16 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 (A  ) is volume of the annulus A  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 16×16 or 32×32 pixels.In Figure 4, higher resolution is used for better illustration.

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 {  };  = 1; 2; . . .;  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 (; ):

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 bluegreen image on Figure 5(a) and orange image on Figure 5(b) would be "perfectly grey."Arithmetic mean of these images 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 1600 × 1200 pixels).This transformation is marginal in our case.

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.
Similar summary is made in Table 3 for the present reconstruction methods: parabolic interpolation of the stairapproximations  ,  , and  .In Table 4, there are summarized RMSE, AD, PCC, and DIE for parabolic and Gaussian interpolation of the same stair-approximation   and proposed method of the expected values of columns {    };  = 1; 2; . . .; .
As was already stated in [35], only small differences are between parabolic and Gaussian interpolation.Really, values of RMSE, AD, and DIE 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 highfrequency 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.

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 Figure 9: 3D reconstruction of the data from Figure 1: preprocessing according to Section 2.5.1 (new method), the first step according to Section 2.5.2 (new method), and the second and the third steps according to Section 2.2 (present method).Expression (28) with  = 16 was used in calculation of the matrix of sharpness detectors, its columns {    };  = 1; 2; . . .;  was interpolated by the Gaussian interpolation (present method).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.1: preprocessing according to Section 2.5.1, the first step according to Section 2.5.2-expression(28) with  = 16 was usedand the second step according to Section 2.5.3, the expected values of its columns {    };  = 1; 2; . . .;  was calculated according to expression (29) where  = 5 was used.

Conflicts of Interest
The author declares that they have no conflicts of interest.

Figure 1 :
Figure 1: Different scaling and different sharp and nonsharp regions in images acquired by classic camera placed in different distances from the 3D relief, the first and the fifteenth images from the series of fifteen images of blue marble.Frame size 10 × 7.5 cm.Locality Nedvědice, Czech Republic, photo Pavel Štarha.

Figure 2 :
Figure 2: Fourier transforms of the fifth member of the series of expanding rectangular signals.The series converges to distribution.
) → R + 0 such that   1 (; ) =   2 (;  + ).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 ( 0 ;  0 ) is measured with another computation of the phase correlation.Scaled Images.Let  2 be function  1 rotated, shifted, and scaled in arguments; that is,  2 (; ) =  1 ( ( cos  −  sin ) −  0 ;  ( sin  +  cos ) −  0 ) .
scale factor   by means of the method described in Section 2.4.Then, image   is rotated by angle −  and scaled by factor  −1  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  * −1 .Now we can apply phase correlation to find the shift ( 0 ;  0 ) and image   is shifted by vector (− 0 ; − 0 ) to compensate the shift, creating image  *  which differs from  * −1 in focussed and blurred parts only.

Figure 5 : 1 S
Figure 5: The first and the fifteenth images from the series of fifteen images of blue marble (see Figure 1) displayed in supplementary pseudocolors.Real surface size 10 × 7.5 cm (b).

Figure 6 :
Figure 6: The arithmetic mean of the images from Figure 5: before registration (a) and after registration (b).Real surface size 10 × 7.5 cm (b).

Figure 7 :
Figure 7: The stairs-approximation of the data from Figure 1 constructed by means of the sharpness detector  ; see expression (28) and Figure 4(c).

Table 4 :Figure 8 :
Figure8: 3D reconstruction of the data from Figure1: preprocessing according to Section 2.5.1 (new method), the first step according to Section 2.5.2 (new method), and the second and the third steps according to Section 2.2 (present method).Expression(28) with  = 16 was used in calculation of the matrix of sharpness detectors, its columns {    };  = 1; 2; . . .;  was interpolated by the parabolic interpolation (present method).

Figure 10 :
Figure 10: Gaussian low-pass filter with dispersion  2 = 5 applied to noisily output data from Figure 9.The noise is indeed reduced (but still visible); however useful high-frequency information is already lost.

Figure 11 :
Figure11: Completely new 3D reconstruction of the data from Figure1: preprocessing according to Section 2.5.1, the first step according to Section 2.5.2-expression(28) with  = 16 was usedand the second step according to Section 2.5.3, the expected values of its columns {    };  = 1; 2; . . .;  was calculated according to expression(29) where  = 5 was used.

Figure
Figure11: Completely new 3D reconstruction of the data from Figure1: preprocessing according to Section 2.5.1, the first step according to Section 2.5.2-expression(28) with  = 16 was usedand the second step according to Section 2.5.3, the expected values of its columns {    };  = 1; 2; . . .;  was calculated according to expression(29) where  = 5 was used.

)
Example 1.One of the well-known examples (which we demonstrate in 1D for simplicity) is the series of expanding rectangular signals  *  () with constant unitary intensity on (−; );  ∈ N and zeroed elsewhere.Inverse Fourier transform gives )

Table 1 :
Parameters of the transforms indicated and applicated to separate images in the series (stated relative to the first image).

Table 2 :
Root Mean Square Error (RMSE) and Average Deviation (AD) in millimeters, Pearson Correlation Coefficient (PCC), and Difference of surface Information Entropy (DIE) (dimensionless quantities) for separate pairs of stairs-approximations  ,  , and  .

Table 3 :
Root Mean Square Error (RMSE), Average Deviation (AD), Pearson Correlation Coefficient (PCC), and Difference of surface Information Entropy (DIE) for separate pairs of parabolic interpolation approximations  ,  , and  .