Mathematical Principles of Object 3D Reconstruction by Shape-from-Focus Methods

mathematics Review Mathematical Principles of Object 3D Reconstruction by Shape-from-Focus Methods Dalibor Martišek * and Karel Mikulášek   Citation: Martišek, D.; Mikulášek, K. Mathematical Principles of Object 3D Reconstruction by Shape-from-Focus Methods. Mathematics 2021,9, 2253. https://doi.org/10.3390/ math9182253 Academic Editors: Akemi Galvez Tomida, Lihua You, Hassan Ugail, Andres Iglesias Prieto and Alexander Malyshev Received: 26 July 2021 Accepted: 7 September 2021 Published: 14 September 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Institute of Mathematics, Faculty of Mechanical Engineering, Brno University of Technology, 61669 Brno, Czech Republic; [email protected].cz *Correspondence: [email protected] Abstract: Shape-from-Focus (SFF) methods have been developed for about twenty years. They able to obtain the shape of 3D objects from a series of partially focused images. The plane to which the microscope or camera is focused intersects the 3D object in a contour line. Due to wave properties of light and due to finite resolution of the output device, the image can be considered as sharp not only on this contour line, but also in a certain interval of height—the zone of sharpness. SSFs are able to identify these focused parts to compose a fully focused 2D image and to reconstruct a 3D profile of the surface to be observed. Keywords: 3D reconstruction; shape-from-focus; optical cut; multifocal image; Fourier transform; phase correlation; focusing criteria 1. Introduction Three-dimensional reconstruction of general surfaces has an important role in a number of fields: the morphological analysis of fracture surfaces, for example, reveals information on the mechanical properties of natural or construction materials. There are more techniques capable of producing digital three-dimensional (3D) replicas of solid surfaces. In mechanical engineering, contacting electronic profilometers can be used to determine digital two-dimensional (2D) profiles to be combined into 3D surface profiles—see [ 1 , 2 ], for example. The contacting mode of atomic force microscopes is actually in this mechanical category [ 3 ]. In addition to the mechanical tools, optical devices exist in diverse modifications [ 2 ], light section microscopy [ 4 , 5 ], coherence scanning interferometry [ 6 ], speckle metrology [ 7 ], stereo projection [ 8 ], photogrammetry [ 9 ], and various types of light measurement of profiles [10], to mention some of them. 3D laser scanning techniques are among other ways of obtaining 3D data. They have also been tested in some rock engineering projects, such as 3D digital fracture mapping [11–13]. These devices are, however, not of universal use, with each of them having its own technical limitations [ 14 , 15 ]. Very rough surfaces, for instance, can hardly be measured by atomic force microscopes working in the nano-regions. On the other hand, it is possible to measure plane surfaces with microscopically small irregularities using a microscopic sectional technique with confocal microscopes [ 16 – 19 ]. Confocal microscopes, however, are not always suitable for technical purposes due to the small size of their visual fields (maximal visual field is about 2 cm [5,20–22]). The present paper summarized the existing methods of 3D reconstruction of objects by the Shape-From-Focus (SFF) method. This is a method for recovering depth from an image series of the same object taken with different focus settings, referred to as a multifocal image. It consists of these steps: (a) Data acquisition (confocal microscope in the standard mode, CCD camera, or standard camera)—Section 2. Mathematics 2021,9, 2253. https://doi.org/10.3390/math9182253 https://www.mdpi.com/journal/mathematics Mathematics 2021,9, 2253 2 of 20 (b) Image registration (if necessary)—Section 3. (c) Choice of the focusing criteria—Section 4. (d) The 2D and 3D reconstructions—Section 5. See [23] for one of the first papers on this subject. 2. Data Acquisition 2.1. Parallel Projection In technical practice, the confocal microscope serves as a standard instrument imaging microscopic three-dimensional surfaces. It has a very small depth of the optical field with its advanced hardware being capable of removing non-sharp points from the images. The points of the object situated close to the focal plane can be seen as sharp points. The parts lying further above or beneath the focal plane (out of the sharpness zone) are invisible, being represented as black regions if the confocal mode is on. In this way, a so-called optical cut is obtained. In the case of the non-confocal (standard) mode, the areas lying outside the sharpness zone are displayed as blurred, as they would be with a standard camera. With a confocal microscope or CCD camera, one can assume that the field of view is small with the projection used being parallel. In this case, all images are provided in the field of view with identical sizes and the corresponding pixels having the same coordinates in separate partial focused images (see Figure 1). Mathematics 2021, 9, x FOR PEER REVIEW 2 of 20 (a) Data acquisition (confocal microscope in the standard mode, CCD camera, or standard camera)—Section 2. (b) Image registration (if necessary)—Section 3. (c) Choice of the focusing criteria—Section 4. (d) The 2D and 3D reconstructions—Section 5. See [23] for one of the first papers on this subject. 2. Data Acquisition 2.1. Parallel Projection In technical practice, the confocal microscope serves as a standard instrument imaging microscopic three-dimensional surfaces. It has a very small depth of the optical field with its advanced hardware being capable of removing non-sharp points from the images. The points of the object situated close to the focal plane can be seen as sharp points. The parts lying further above or beneath the focal plane (out of the sharpness zone) are invisible, being represented as black regions if the confocal mode is on. In this way, a so-called optical cut is obtained. In the case of the non-confocal (standard) mode, the areas lying outside the sharpness zone are displayed as blurred, as they would be with a standard camera. With a confocal microscope or CCD camera, one can assume that the field of view is small with the projection used being parallel. In this case, all images are provided in the field of view with identical sizes and the corresponding pixels having the same coordinates in separate partial focused images (see Figure 1). (a) (b) Figure 1. Optical cut of fracture surface of hydrated cement paste acquired by confocal microscope Olympus LEXT 1000. Confocal mode (a), standard mode (b). Taken from [24]). However, the confocal microscope with its small visual angle is hardly suitable for technical purposes due to the small size of its visual field (maximal visual field is about 2 cm [5,20,21,24]). 2.2. Central Projection The same output (Figure 1b) can be obtained by the classical microscope or (in a wider field) common camera. The difference between the microscope or CCD camera and standard camera is given by the central projection that varies the scaling of partial images in a series and, further, by the non-sharp regions, displayed by the classic camera, while missing when taken by a confocal microscope in the confocal mode (see [25] for more information). The different image scalings, however, require subsequent corrections (including shifts and rotations). Figure 1. Optical cut of fracture surface of hydrated cement paste acquired by confocal microscope Olympus LEXT 1000. Confocal mode (a), standard mode (b). Taken from [24]. However, the confocal microscope with its small visual angle is hardly suitable for technical purposes due to the small size of its visual field (maximal visual field is about 2 cm [5,20,21,24]). 2.2. Central Projection The same output (Figure 1b) can be obtained by the classical microscope or (in a wider field) common camera. The difference between the microscope or CCD camera and standard camera is given by the central projection that varies the scaling of partial images in a series and, further, by the non-sharp regions, displayed by the classic camera, while missing when taken by a confocal microscope in the confocal mode (see [ 25 ] for more information). The different image scalings, however, require subsequent corrections (including shifts and rotations). 2.3. Multifocal Image To create a 2D or 3D reconstruction, it is necessary to obtain a series of images of an identical object, each of them with different focusing with each object point being focused Mathematics 2021,9, 2253 3 of 20 in one of the images (in the ideal case, this is referred to as a multifocal image). For the acquisition of a large multifocal image, the camera must be mounted on a stand so that it can be moved in the direction approximately orthogonal to the surface with a controlled step. Different transformations and different sharp parts must be identified and composed in a 2D or 3D model. 3. Image Registration With a confocal microscope or CCD camera, we can assume that the field of view is small and the projection used is parallel. In this case, all images are provided in the field of view with identical sizes with the corresponding pixels having the identical coordinates in separate partial focused images. However, this assumption does not hold for larger samples; then, the angle of the projection lines is not negligible with the view fields (and the coordinates of the corresponding pixels) being clearly different for each image (see Figures 2and 3). Mathematics 2021, 9, x FOR PEER REVIEW 3 of 20 2.3. Multifocal Image To create a 2D or 3D reconstruction, it is necessary to obtain a series of images of an identical object, each of them with different focusing with each object point being focused in one of the images (in the ideal case, this is referred to as a multifocal image). For the acquisition of a large multifocal image, the camera must be mounted on a stand so that it can be moved in the direction approximately orthogonal to the surface with a controlled step. Different transformations and different sharp parts must be identified and composed in a 2D or 3D model. 3. Image Registration With a confocal microscope or CCD camera, we can assume that the field of view is small and the projection used is parallel. In this case, all images are provided in the field of view with identical sizes with the corresponding pixels having the identical coordinates in separate partial focused images. However, this assumption does not hold for larger samples; then, the angle of the projection lines is not negligible with the view fields (and the coordinates of the corresponding pixels) being clearly different for each image (see Figures 2 and 3). (a) (b) Figure 2. Different scaling and different sharp and non-sharp regions in images acquired by the classic camera positioned at different distances from the 3D relief of the first (a) and the sixteenth (b) image of a series of sixteen images, blue marble, locality (Nedvedice, Czech Republic, photo Pavel Starha). Taken from [26]. (a) (b) Figure 3. Different scaling and different sharp and non-sharp regions in images acquired by the classic camera positioned at different distances from the 3D relief—the first (a) and the fourth (b) image of a series of eight images, limestone, locality Brno (Hady, Czech Republic, photo Tomas Ficker). Taken from [26]. Before reconstruction, all geometric transformations must be identified in the image series to be eliminated. The images are being analyzed as geometrically similar. Generally, Figure 2. Different scaling and different sharp and non-sharp regions in images acquired by the classic camera positioned at different distances from the 3D relief of the first ( a ) and the sixteenth ( b ) image of a series of sixteen images, blue marble, locality (Nedvedice, Czech Republic, photo Pavel Starha). Taken from [26]. Mathematics 2021, 9, x FOR PEER REVIEW 3 of 20 2.3. Multifocal Image To create a 2D or 3D reconstruction, it is necessary to obtain a series of images of an identical object, each of them with different focusing with each object point being focused in one of the images (in the ideal case, this is referred to as a multifocal image). For the acquisition of a large multifocal image, the camera must be mounted on a stand so that it can be moved in the direction approximately orthogonal to the surface with a controlled step. Different transformations and different sharp parts must be identified and composed in a 2D or 3D model. 3. Image Registration With a confocal microscope or CCD camera, we can assume that the field of view is small and the projection used is parallel. In this case, all images are provided in the field of view with identical sizes with the corresponding pixels having the identical coordinates in separate partial focused images. However, this assumption does not hold for larger samples; then, the angle of the projection lines is not negligible with the view fields (and the coordinates of the corresponding pixels) being clearly different for each image (see Figures 2 and 3). (a) (b) Figure 2. Different scaling and different sharp and non-sharp regions in images acquired by the classic camera positioned at different distances from the 3D relief of the first (a) and the sixteenth (b) image of a series of sixteen images, blue marble, locality (Nedvedice, Czech Republic, photo Pavel Starha). Taken from [26]. (a) (b) Figure 3. Different scaling and different sharp and non-sharp regions in images acquired by the classic camera positioned at different distances from the 3D relief—the first (a) and the fourth (b) image of a series of eight images, limestone, locality Brno (Hady, Czech Republic, photo Tomas Ficker). Taken from [26]. Before reconstruction, all geometric transformations must be identified in the image series to be eliminated. The images are being analyzed as geometrically similar. Generally, Figure 3. Different scaling and different sharp and non-sharp regions in images acquired by the classic camera positioned at different distances from the 3D relief—the first ( a ) and the fourth ( b ) image of a series of eight images, limestone, locality Brno (Hady, Czech Republic, photo Tomas Ficker). Taken from [26]. Before reconstruction, all geometric transformations must be identified in the image series to be eliminated. The images are being analyzed as geometrically similar. Generally, similarity is achieved by composing rotation, scale-change, shift, and axial symmetry. Axial symmetry is not possible in this particular case. Mathematics 2021,9, 2253 4 of 20 If we consider scale-change only assuming that the image size is proportional to the camera shifting (see Figure 4a), different image scaling can be obtained using elementary mathematics. This approach was used in [ 25 ]. After the initial elementary reconstruction, a subsequent 3D reconstruction is shown in Figure 5. Huge artefacts caused by inaccurate registration were blurred by brutal low-pass filters with much useful high-frequency information lost. Mathematics 2021, 9, x FOR PEER REVIEW 4 of 20 similarity is achieved by composing rotation, scale-change, shift, and axial symmetry. Axial symmetry is not possible in this particular case. If we consider scale-change only assuming that the image size is proportional to the camera shifting (see Figure 4a), different image scaling can be obtained using elementary mathematics. This approach was used in [25]. After the initial elementary reconstruction, a subsequent 3D reconstruction is shown in Figure 5. Huge artefacts caused by inaccurate registration were blurred by brutal low-pass filters with much useful high-frequency information lost. (a) (b) Figure 4. The central projection of a large sample—ideal case (a) can be solved by elementary mathematics, real case (b) necessitates sophisticated mathematical tools. Taken from [26]. Figure 5. 3D reconstruction of the data from Figure 3 after elementary registration by Figure 4a. Reconstruction taken from [27]. Software developed by the author. Figure 4. The central projection of a large sample—ideal case ( a ) can be solved by elementary mathematics, real case (b) necessitates sophisticated mathematical tools. Taken from [26]. Mathematics 2021, 9, x FOR PEER REVIEW 4 of 20 similarity is achieved by composing rotation, scale-change, shift, and axial symmetry. Axial symmetry is not possible in this particular case. If we consider scale-change only assuming that the image size is proportional to the camera shifting (see Figure 4a), different image scaling can be obtained using elementary mathematics. This approach was used in [25]. After the initial elementary reconstruction, a subsequent 3D reconstruction is shown in Figure 5. Huge artefacts caused by inaccurate registration were blurred by brutal low-pass filters with much useful high-frequency information lost. (a) (b) Figure 4. The central projection of a large sample—ideal case (a) can be solved by elementary mathematics, real case (b) necessitates sophisticated mathematical tools. Taken from [26]. Figure 5. 3D reconstruction of the data from Figure 3 after elementary registration by Figure 4a. Reconstruction taken from [27]. Software developed by the author. Figure 5. 3D reconstruction of the data from Figure 3after elementary registration by Figure 4a. Reconstruction taken from [27]. Software developed by the author. Mathematics 2021,9, 2253 5 of 20 In practice, the situation may be more sophisticated. The images may differ not just in the scale used but in the content displayed as well (different parts being focused in different images). Due to mechanical inaccuracies, the step in the z axis may be not fully constant, and the images can also be mutually shifted along the xor y-axis or rotated. Image registration is also complicated by the non-planarity of samples (see Figure 4on the right). Therefore, sophisticated pre-processing of the image series may be necessary. A method suitable tool for this is the Fourier transform and phase correlation. 3.1. Continuous Two-Dimensional Fourier Transform and Inverse Transform A continuous standard Fourier transform of a function f(x):R→Cis the function F[f(x)] =F(ξ)= ∞ Z −∞ f(x)e−ixξdx(1) (provided that this integral exists and is finite). A continuous standard Fourier transform of a function f(x;y):R2→C is the function F[f(x;y)] =F(ξ;η)=x R2 f(x;y)e−i(xξ+yη)dxdy(2) (provided that this integral exists and is finite). Function F is also referred to as the Fourier spectrum of function f . It is possible to obtain the function ffrom its Fourier spectrum Fby the inverse Fourier transform. A continuous standard inverse Fourier transform of a function G(ξ):C→C is the function F−1[F(ξ)] =f(x)=1 2π ∞ Z −∞ F(ξ)eixξdξ(3) (provided that the integral exists and is finite). A continuous standard inverse Fourier transform of a function G(ξ;η):C2→C is the function F−1[F(ξ;η)] =f(x;y)=1 4π2x R2 F(ξ;η)ei(xξ+yη)dξdη(4) (provided that this integral exists and is finite). 3.2. Discrete Two-Dimensional Fourier Transform and Inverse Transform A discrete standard Fourier transform of a function f(x;y) : {0; 1; . . . ; M−1}×{0; 1; . . . ; N−1}→Cis the function F[f(x;y)] =F(ξ;η)= M−1 ∑ x=0 N−1 ∑ y=0 f(x;y)e−2πi(xξ M+yη N)(5) A discrete standard inverse Fourier transform of a function G(ξ;η) : {0; 1; . . . ; M−1}×{0; 1; . . . ; N−1}→Cis the function F−1[G(ξ;η)] =g(x;y)=1 MN M−1 ∑ ξ=0 N−1 ∑ η=0 G(ξ;η)e2πi(xξ M+yη N)(6) If the function f(x;y)in (5) is real, then F(N−ξ;N−η)=F(ξ;η)(7) (where the bar denotes complex conjugation) and the Fast Fourier Transform (FFT) algorithm can be employed to calculate the discrete Fourier transform—see [28] for more information. Mathematics 2021,9, 2253 6 of 20 The discrete Fourier transform can be used for the image registration applied to functions that are or are assumed periodic. Generally, an image may not have the same values on the edges. Thus, by periodizing an image, the resulting function may have jumps at the edges of the original image. Such jumps are often the most contrasted structures in the function and may lead to incorrect registration. Therefore, such edges used for the shift estimation must be removed from the images. This is carried out by multiplying the image by a suitable function referred to as a window function. Its values must equal zero or almost zero at the image edges and one on a large part of the image. Primarily, the Gaussian and Hanning window functions can be used: Let σ∈R+ be a given number and sets R=h−a;ai × h−b;bi;a;b∈R+ 0 C=n(x;y)∈R2x2+y2≤ro;r∈R+ 0 Let ρ(X; S)be the distance of point X=(x;y)from set S; i.e., ρ(X; S)=inf{d∈R|d=kXYk;Y∈S}(8) Functions gGR(x;y)=e−ρ2(X;R) σ2or gGC(x;y)=e−ρ2(X;C) σ2(9) are called rectangular or circular Gaussian window functions. Functions gHR(x;y)=(1 2+1 2cos πρ(X;R) σif (X;R)≤σ 0 if (X;R)>σ or gHC(x;y)=(1 2+1 2cos πρ(X;C) σif (X;C)≤σ 0 if (X;C)>σ(10) are called rectangular or circular Hanning window functions. 3.3. δ-Distribution A one-dimensional δ -distribution δ(x) is a limit of a function sequence δn(x) ; n∈N for which (a)lim n→∞ ∞ Z −∞ δn(x)dx= (b)lim n→∞ δn(x0) lim x→0δn(x)=0; x0∈R−{0}(11) A two-dimensional δ -distribution δ(x;y) is a limit of a function sequence δn(x;y);n∈N for which lim n→∞ ∞ Z −∞ ∞ Z −∞ δn(x;y)dxdy=1; lim n→∞ δn(x0;y0) lim (x;y)→(0;0)δn(x;y)=0; (x0;y0)∈R2−{(0; 0)}(12) Example: A well-known example (here demonstrated in 1 D for simplicity) is a series of expanding rectangular signals δ∗ n(ξ) with a unitary intensity constant on (−n;n) ; n∈N and zeroed elsewhere. The inverse Fourier transform gives F−1(δ∗ n)(ξ)=δn(x)=1 2π ∞ R −∞ eixξdξ=1 2π n R −n eixξdξ=1 2πheixξ ix in ξ=−n =eixn−e−ixn 2πix =sin nx πx Mathematics 2021,9, 2253 7 of 20 As R∞ −∞sin nx πxdx= 1 for each n∈N (see [ 28 ] for proof), condition (a) is fulfilled. Since lim x→0 sin nx πx=n πlim x→0 sin nx nx =n πfor each n∈N, wehave lim n→∞ δn(x0) lim x→0δn(x)=lim n→∞ sin nx0 nx0 =0 for each x06=0. This means that condition (b) holds as well and the limit δ(x)=lim n→∞δn(x) is the (one-dimensional) δ -distribution. The expanded unitary signal δ∗ 5(ξ) with its Fourier transform is illustrated in Figure 6. Mathematics 2021, 9, x FOR PEER REVIEW 7 of 20 F −1(𝛿𝑛∗)(𝜉)=𝛿𝑛(𝑥)=1 2𝜋∫e𝑖𝑥𝜉 ∞ − ∞ d𝜉=1 2𝜋∫e𝑖𝑥𝜉 𝑛 −𝑛 d𝜉=1 2𝜋[e𝑖𝑥𝜉 𝑖𝑥]𝜉=−𝑛 𝑛 =e𝑖𝑥𝑛−e−𝑖𝑥𝑛 2𝜋𝑖𝑥 =sin𝑛𝑥 𝜋𝑥 As ∫sin𝑛𝑥 𝜋𝑥 d𝑥 ∞ −∞ =1 for each 𝑛∈ℕ (see [28] for proof), condition (a) is fulfilled. Since lim 𝑥→0sin𝑛𝑥 𝜋𝑥 =𝑛𝜋lim 𝑥→0sin𝑛𝑥 𝑛𝑥 =𝑛𝜋 for each 𝑛∈ℕ,we have lim 𝑛→∞ 𝛿𝑛(𝑥0) lim 𝑥→0𝛿𝑛(𝑥)=lim 𝑛→∞sin𝑛𝑥0 𝑛𝑥0=0 for each 𝑥0≠0. This means that condition (b) holds as well and the limit 𝛿(𝑥)=lim 𝑛→∞𝛿𝑛(𝑥) is the (one-dimensional) 𝛿−distribution. The expanded unitary signal 𝛿5∗(𝜉) with its Fourier transform is illustrated in Figure 6. Figure 6. Fourier transforms of the fifth term of the series of expanding rectangular signals. The series 𝜹𝒏∗(𝝃) converges to the 𝛅-distribution. Taken from [26]. 3.4. Phase Correlation Image processing requires the images transformed for the structures studied to be at the same position in all of them. The transformation is found by image registration. In some applications, it is possible to assume shift only while, in others, shift, rotation and scale change (i.e., similarity), general linear transformation, or even general transformations may all be present. The methods used for registration depend on the expected transformation and on the image structures. Some methods, after using the corresponding structures or points in the images, find a global transformation by measuring the positions of the structures or points [29–31]. For these methods to be applicable, the structures must be clearly visible. Other correlation-based methods work with the image as a whole. The phase correlation has proved to be a powerful tool (not only) for the registration of particular focused images. For functions 𝑓1; 𝑓2, it is defined as 𝑃𝑓1;𝑓2(𝑥;𝑦)=ℱ−1{𝐹1(𝜉;𝜂)∙𝐹2(𝜉;𝜂) |𝐹1(𝜉;𝜂)|∙|𝐹2(𝜉;𝜂)|} (13) with its modification being Figure 6. Fourier transforms of the fifth term of the series of expanding rectangular signals. The series δ∗ n(ξ)converges to the δ-distribution. Taken from [26]. 3.4. Phase Correlation Image processing requires the images transformed for the structures studied to be at the same position in all of them. The transformation is found by image registration. In some applications, it is possible to assume shift only while, in others, shift, rotation and scale change (i.e., similarity), general linear transformation, or even general transformations may all be present. The methods used for registration depend on the expected transformation and on the image structures. Some methods, after using the corresponding structures or points in the images, find a global transformation by measuring the positions of the structures or points [ 29 – 31 ]. For these methods to be applicable, the structures must be clearly visible. Other correlation-based methods work with the image as a whole. The phase correlation has proved to be a powerful tool (not only) for the registration of particular focused images. For functions f1;f2, it is defined as Pf1;f2(x;y)=F−1F1(ξ;η)·F2(ξ;η) |F1(ξ;η)|·|F2(ξ;η)|(13) with its modification being Pf1;f2;;p;q(x;y)=F−1H(ξ;η)·F1(ξ;η)·F2(ξ;η) (|F1(ξ;η)|+p)·(|F2(ξ;η)|+q)(14) where the bar denotes complex conjugation, H(ξ;η) is a bounded real function such that H(ξ;η)=H(−ξ;−η) and p ; q> 0 are arbitrary constants. It is not difficult to prove that, for real functions f1 , f2 , the phase-correlation function is real [ 32 ]. This is very useful since the extremes of the phase-correlation function can be searched for. Mathematics 2021,9, 2253 8 of 20 3.5. Identical Images Let F be the infinity periodic expansion of an image. Denote a+bi its value F(ξ;η) at (ξ;η),a+bi 6=0. Clearly, the value of the phase correlation of the Fwith itself is F(ξ;η)·F(ξ;η) |F(ξ;η)·F(ξ;η)|=(a+bi)(a−bi) |(a+bi)|·|a+bi|=a2+b2 a2+b2=1 (15) By the example in Section 3.3, we have Pf;f(x;y)=F−1F(ξ;η)·F(ξ;η) |F(ξ;η)|·|F(ξ;η)|=F−1{1}=δ(x;y)(16) which means that the inverse Fourier transform of the correlation of two identical images is the two-dimensional δ-distribution δ(x;y) 3.6. Shifted Images If two functions are shifted in arguments, that is, f2(x;y)=f1(x−x0;y−y0) , their Fourier transforms are shifted in phase: F2(ξ;η)=F1(ξ;η)·exp(−i(ξx0+ηy0)) (17) with their phase-correlation function being the δ -distribution shifted in arguments by the opposite shift vector Pf1;f2(x;y)=F−1{exp(i(ξx0+ηy0))}=δ(x+x0;y+y0)(18) This is the principal idea of phase correlation. Using phase correlation, rather than finding a shift between two images, we can just find the only non-zero point in a matrix. If the images are not identical (up to a shift), i.e., if the images are not ideal, the phasecorrelation function is more complex, but still has a global maximum at the point whose coordinates correspond to the shift vector. 3.7. Rotated Images The phase-correlation function can also be used for estimating the image rotation and rescale. Let f2be function f1rotated and shifted in arguments, i.e., f2(x;y)=f1(xcos θ−ysin θ−x0;xsin θ+ycos θ−y0)(19) Their Fourier spectra F1;F2and amplitude spectra A1;A2are related as follows: F2(ξ;η)=exp(−i(ξx0+ηy0))·F1(ξcos θ−ηsin θ;ξsin θ+ηcos θ)(20) A2(ξ;η)=A1(ξcos θ−ηsin θ;ξsin θ+ηcos θ)(21) The shift results in a phase shift and the spectra are rotated the same way as the original functions. A crucial step here is the transformation of the amplitude spectra into polar coordinates to obtain functions Ap 1;Ap 2:R+ 0×0; 2π)→R+ 0 such that Ap 1(ρ;ϕ)=Ap 2(ρ;ϕ+θ) . The rotation about an unknown centre has been transformed into a shift. This shift is estimated by the standard phase correlation (see the previous paragraph) after a reverse rotation by the angle measured, the shift (x0;y0) is then measured in another computation of the phase correlation. 3.8. Scaled Images Let f2be function f1rotated, shifted, and scaled in arguments, i.e., f2(x;y)=f1(α(xcos θ−ysin θ)−x0;α(xsin θ+ycos θ)−y0)(22) Mathematics 2021,9, 2253 9 of 20 Their Fourier spectra and amplitude spectra are related as follows: F2(ξ;η)=1 α2exp(−i(ξx0+ηy0))·F11 α(ξcos θ−ηsin θ);1 α(ξsin θ+ηcos θ)(23) A2(ξ;η)=1 α2A11 α(ξcos θ−ηsin θ);1 α(ξsin θ+ηcos θ)(24) The shift results in a phase shift, the spectra are rotated the same way as the original functions and scaled with a reciprocal factor. A crucial step here is the transformation of the amplitude spectra into the logarithmic-polar coordinates exp ρ=qx2+y2;x=exp ρcos ϕ;y=exp ρsin ϕ(25) to obtain Ap 1;Ap 2:R+ 0× h0; 2π)→R+ 0such that A1p 2(ρ;ϕ)=A1p 2(ρ−ln α;ϕ+θ). Both rotation and scale change have been transformed to a shift. The unknown angle θ and unknown factor α can be estimated by means of the phase correlation applied to the amplitude spectra in the logarithmic-polar coordinates A1p 1 ; A1p 2 . After reverse rotating function f2 by the estimated angle θ and scaling by the factor, the shift vector (x0;y0) can be estimated by means of the standard phase correlation. 3.9. Multifocal Registration Let {P1;P2; . . . ; Pn} be the image series to be registered with image P1 acquired by means of the biggest angle of view. This image will be not transformed or (formally) it will be transformed by the identity mapping into the image P∗ 1 . Now we must find the transform P2→P∗ 1 to obtain image P∗ 2 which only differs from P∗ 1 in focussed and blured parts. In the same way, transforms P3→P∗ 2 ; . . . ; Pk→P∗ k−1; . . . ; Pn→P∗ n−1 must be found. After multiplying both images Pk ; P∗ k−1 by the chosen window function, rotation angle θk and scale factor αk will be determined by the method described in Section 3.8. Then, image Pk is rotated by the angle −θk and scaled by the factor α−1 k to compensate for the rotation and scale-change found by the phase correlation, creating image Pk . Between images Pk . and P∗ k−1 , only the shifted and different focused and blurred parts remain. Now we can apply phase correlation to find the shift (x0;y0) shifting image Pk by the vector (−x0;−y0) to compensate for the shift, creating image P∗ k which only differs from P∗ k−1 in the focused and blurred parts. 4. Focusing Criteria The detectors of blurred areas (sometimes referred to as focusing criteria or sharpness detectors) can be based on different principles. Probably the first attempts to carry out non-confocal reconstructions date back to 1970’s and 1980’s [33–36]. Tenebaum [ 37 ] developed the gradient magnitude maximization method for optimizing the focus quality using the sharpness of edges. Jarvis [ 38 ] proposed a sum-modulusdifference computed by summing the first intensity differences between neighbouring pixels along a scan-line using it as a focus quality benchmark. Schlag et al. [ 39 ] implemented and tested various self-focusing algorithms. Recently, Krotkov [ 40 ] evaluated and compared the performance of different focus criterion functions. In [ 41 ], he also proposed a method for estimating the depth of an image area. Pentland [ 42 ] suggested evaluating the image blur to determine the depth of image points. Grossmann in [ 43 ] proposed estimating the depth of edge points by analyzing the blur of the edges due to defocusing. Darrell and Wohn in [ 29 ] developed a depth-from-focus method by which an image sequence can be obtained through varying the focus level using Laplacian and Gaussian pyramids to calculate the depth. Subbarao in [ 30 ] suggested changing the intrinsic camera parameters to recover the depth map of a scene. Ohta et al. [ 31 ] and Kaneda et al. [ 44 ] used images corresponding to different focus levels to obtain a single level of high focus quality. Mathematics 2021,9, 2253 16 of 20 Mathematics 2021, 9, x FOR PEER REVIEW 16 of 20 Figure 14. A 3D echelon approximation of the limestone sample by the optical cuts used in Figure 13 smoothed by 3D low-pass filters. The software was written by the first author. Taken from [45]. 6.4. 2D and 3D Reconstructions We can see optical cuts detected in the data of Figure 2 using focusing criteria (28) in Figure 11, 2D reconstruction (sharp 2D Image) of the same data using the same criterion as in Figure 12. Echelon approximation is a simple method for constructing a rough 3D model of the object, where all points belonging to the same optical cut have the same height—the height of the corresponding zone of sharpness—see Figure 13. We can also generalize the notion of low−pass filters used in image processing. They can be used for smoothing an echelon approximation. This approximation looks much better than the echelon one—see Figure 14. In Figures 15 and 16, we can see reconstructions of the limestone and blue marble sample by data registration according to Section 3, with focusing criterion (28) and profile height calculation (29) and (30). We can compare a single pore of hydrated Portland cement paste reconstructed by Olympus factory software (Figure 17) and the same pore constructed using focusing criterion (28), probability distribution function (29) and (30) for r = 2 (Figure 18). Two examples of SFF method results can be downloaded at https://dmartisek.cz/Veda/Two_Results_of_SFF_Method.avi. Figure 14. A 3D echelon approximation of the limestone sample by the optical cuts used in Figure 13 smoothed by 3D low-pass filters. The software was written by the first author. Taken from [45]. 6.4. 2D and 3D Reconstructions We can see optical cuts detected in the data of Figure 2using focusing criteria (28) in Figure 11, 2D reconstruction (sharp 2D Image) of the same data using the same criterion as in Figure 12. Echelon approximation is a simple method for constructing a rough 3D model of the object, where all points belonging to the same optical cut have the same height—the height of the corresponding zone of sharpness—see Figure 13. We can also generalize the notion of low − pass filters used in image processing. They can be used for smoothing an echelon approximation. This approximation looks much better than the echelon one—see Figure 14. In Figures 15 and 16, we can see reconstructions of the limestone and blue marble sample by data registration according to Section 3, with focusing criterion (28) and profile height calculation (29) and (30). We can compare a single pore of hydrated Portland cement paste reconstructed by Olympus factory software (Figure 17) and the same pore constructed using focusing criterion (28), probability distribution function (29) and (30) for r = 2 (Figure 18). Mathematics 2021,9, 2253 17 of 20 Mathematics 2021, 9, x FOR PEER REVIEW 17 of 20 Figure 15. A 3D reconstruction of the limestone sample by data registration according to Section 4, with focusing criterion 28 and profile height calculations 29 and 30 (compare with Figure 5). The software was written by the first author. Taken from [45]. Figure 16. A 3D reconstruction of the blue marble sample by data registration according to Section 4, with focusing criterion (28) and profile height calculations (29) and (30). The software was written by the first author. Taken from [45]. Figure 15. A 3D reconstruction of the limestone sample by data registration according to Section 4, with focusing criterion 28 and profile height calculations 29 and 30 (compare with Figure 5). The software was written by the first author. Taken from [45]. Mathematics 2021, 9, x FOR PEER REVIEW 17 of 20 Figure 15. A 3D reconstruction of the limestone sample by data registration according to Section 4, with focusing criterion 28 and profile height calculations 29 and 30 (compare with Figure 5). The software was written by the first author. Taken from [45]. Figure 16. A 3D reconstruction of the blue marble sample by data registration according to Section 4, with focusing criterion (28) and profile height calculations (29) and (30). The software was written by the first author. Taken from [45]. Figure 16. A 3D reconstruction of the blue marble sample by data registration according to Section 4, with focusing criterion (28) and profile height calculations (29) and (30). The software was written by the first author. Taken from [45]. Mathematics 2021,9, 2253 18 of 20 Mathematics 2021, 9, x FOR PEER REVIEW 18 of 20 Figure 17. A confocal 3D relief of a single pore of hydrated Portland cement paste; 47 optical cuts with a vertical stepping of 1.2 μm. Olympus LEXT 1000, confocal mode, Olympus factory software. Taken from [27]. Figure 18. The same single pore of hydrated Portland cement paste as in Figure 17; 47 optical cuts with a vertical stepping of 1.2 μm. Olympus LEXT 1000 again. Non−confocal mode. A 3D reconstruction by data registration according to Section 4, focusing criterion (28), and profile height calculations (29) and (30). The software was written by the first author. Taken from [27]. 7. Conclusions The SFF method based on the Fourier transform can provide correct 3D replicas of rough surfaces. In the case of small samples, a qualified user of this method can obtain results similar to or even better than reconstructions from a confocal microscope. For larger objects, 3D scanners and similar significantly more expensive devices can be simulated by sophisticated mathematical instruments and advanced programming techniques. Author Contributions:Conceptualization, D.M.; methodology, D.M.; investigation D.M.; software, D.M.; resources, D.M.; data curation, D.M.; writing—original draft preparation, D.M.; visualization, D.M.; project administration, D.M.; funding acquisition, D. M.; validation, K.M.; formal analysis, Figure 17. A confocal 3D relief of a single pore of hydrated Portland cement paste; 47 optical cuts with a vertical stepping of 1.2 µ m. Olympus LEXT 1000, confocal mode, Olympus factory software. Taken from [27]. Mathematics 2021, 9, x FOR PEER REVIEW 18 of 20 Figure 17. A confocal 3D relief of a single pore of hydrated Portland cement paste; 47 optical cuts with a vertical stepping of 1.2 μm. Olympus LEXT 1000, confocal mode, Olympus factory software. Taken from [27]. Figure 18. The same single pore of hydrated Portland cement paste as in Figure 17; 47 optical cuts with a vertical stepping of 1.2 μm. Olympus LEXT 1000 again. Non−confocal mode. A 3D reconstruction by data registration according to Section 4, focusing criterion (28), and profile height calculations (29) and (30). The software was written by the first author. Taken from [27]. 7. Conclusions The SFF method based on the Fourier transform can provide correct 3D replicas of rough surfaces. In the case of small samples, a qualified user of this method can obtain results similar to or even better than reconstructions from a confocal microscope. For larger objects, 3D scanners and similar significantly more expensive devices can be simulated by sophisticated mathematical instruments and advanced programming techniques. Author Contributions:Conceptualization, D.M.; methodology, D.M.; investigation D.M.; software, D.M.; resources, D.M.; data curation, D.M.; writing—original draft preparation, D.M.; visualization, D.M.; project administration, D.M.; funding acquisition, D. M.; validation, K.M.; formal analysis, Figure 18. The same single pore of hydrated Portland cement paste as in Figure 17; 47 optical cuts with a vertical stepping of 1.2 µ m. Olympus LEXT 1000 again. Non − confocal mode. A 3D reconstruction by data registration according to Section 4, focusing criterion (28), and profile height calculations (29) and (30). The software was written by the first author. Taken from [27]. Two examples of SFF method results can be downloaded in Supplementary Material. 7. Conclusions The SFF method based on the Fourier transform can provide correct 3D replicas of rough surfaces. In the case of small samples, a qualified user of this method can obtain results similar to or even better than reconstructions from a confocal microscope. For larger objects, 3D scanners and similar significantly more expensive devices can be simulated by sophisticated mathematical instruments and advanced programming techniques. Mathematics 2021,9, 2253 19 of 20 Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/math9182253/s1. Author Contributions: Conceptualization, D.M.; methodology, D.M.; investigation D.M.; software, D.M.; resources, D.M.; data curation, D.M.; writing—original draft preparation, D.M.; visualization, D.M.; project administration, D.M.; funding acquisition, D.M.; validation, K.M.; formal analysis, K.M.; writing—review and editing, K.M.; supervision, K.M. All authors have read and agreed to the published version of the manuscript. Funding: This research has been funded by Private Institute of Applied Mathematics, Slapanice, Czech Republic. Data Availability Statement: Data supporting reported results can be downloaded at https:// dmartisek.cz/Veda/Blue_Marble_Data.zip;https://dmartisek.cz/Veda/Limestone_Data.zip;https: //dmartisek.cz/Veda/Two_Results_of_SFF_Method.avi. Conflicts of Interest: The authors declare no conflict of interest. The funder has had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. References 1. Halling, J. Introduction to Tribology; John Wiley & Sons: London, UK, 1976. 2. Bennett, J.M.; Matton, L. Introduction to Surface Roughness and Scattering; Optical Society of America: Washington, DC, USA, 1999. 3. Bowen, W.R. Atomic Force Microscopy in Process Engineering; Hilal, N., Ed.; Butterworth-Heinemann: Oxford, UK, 2009. 4. Tolansky, S. A light-profile microscope for surface studies. Z. Elektrochem. 1952,56, 263–267. 5. Thiery, V.; Green, D.I. The multifocus imaging technique in petrology. CompGeosci 2012,45, 131–138. [CrossRef] 6. De Groot, P. Principles of interference microscopy for the measurement of surface topography. Adv. Opt. Photonics 2015 ,7, 1–65. [CrossRef] 7. Kaufmann, G.H. Advances in Speckle Metrology and Related Techniques; Wiley: Weinheim, Germany, 2011. 8. Mettänen, M.; Hirn, U. A comparison of five optical surface topography measurement methods. TAPPI J. 2015 ,14, 27–38. [CrossRef] 9. Bertin, S.; Friedrich, H.; Dekmas, P.; Chan, E.; Gimel’farb, G. Digital stereo photogrammetry for grain-scale monitoring offluvialsurfaces: Error evaluation and work flow optimization. ISPRS J. 2015,101, 193–208. 10. Tang, S.; Zhang, X.; Tu, D. Micro-phase measuring profilometry: Its sensitivity analysis and phase unwrapping. Opt. Lasers Eng. 2015,72, 47–57. [CrossRef] 11. Feng, Q. Novel Methods for 3-D Semi-Automatic Mapping of Fracture Geometry at Exposed Rock Surfaces. Ph.D. Thesis, KTH, Stockholm, Sweden, 2001. 12. Slob, S.; Hack, H.R.G.K.; Van Knapen, B.; Turner, K.; Kemeny, J. A method for automated discontinuity analysis of rock slopes with three-dimensional laser scanning. Transp. Res. Rec. J. Transp. Res. Board. 2005,1913, 187–194. [CrossRef] 13. Slob, S.; Hack, H.R.G.K. 3D terrestrial laser scanning as a new field measurement and monitoring technique. In Engineering Geology for Infrastructure Planning in Europe. A European Perspective; Azzam, R.H.R.a., Charlier, R., Eds.; Springer: Berlin/Heidelberg, Germany, 2004; pp. 179–190. 14. Pawlus, P.; Wieczorowski, M.; Mathia, T. The Errors of Stylus Methods in Surface Topography Measurements; ZAPOL: Szczecin, Poland, 2014. 15. Hoła, J.; Sadowski, Ł.; Reiner, J.; Stach, S. Usefulness of 3D surface roughness parameters for nondestructive evaluation of pull-off adhesion of concrete layers. Constr. Build. Mater. 2015,84, 111–120. [CrossRef] 16. Agard, D.A.; Hiraoka, Z.; Shaw, P.; Sedat, J. Fluorescence microscopy in three dimensions, in Methods in Cell Biology. In Fluorescence Microscopy of Living Cells in Culture: Part B: Quantitative Fluorescence Microcopy-Imaging and Spectroscopy; Taylor, D.L., Wang, Y., Eds.; Academic Press: San Diego, CA, USA, 1989; Volume 30, pp. 359–362. 17. Wilson, T. (Ed.) Confocal Microscopy; Academic Press Limited: London, UK, 1990. 18. Pawley, J.B. Handbook of Confocal Microscopy; Plenum Press: New York, NY, USA, 1990. 19. Logali, N. Confocal Laser Microscopy-Principles, Applications in Medicine, Biology, and the Food Sciences; InTech-open access publisher: Rijeka, Croatia, 2013. 20. Lange, D.; Jennings, H.M.; Shah, S.P. Analysis of surface roughness using confocal microscopy. J. Mater. Sci. 1993 ,28, 3879–3884. [CrossRef] 21. Ichikawa, Y.; Toriwaki, J.-I. Confocal Microscope 3d Visualizing Method for Fine Surface Characterization of Microstructures; International Society for Optics and Photonics: Denver, CO, USA, 1996. 22. Nadolny, K. Confocal laser scanning microscopy for characterization of surface micro discontinuities of vitrified bonded abrasive tools. Int. J. Mech. Eng. Robot. Res. 2012,1, 14–29. 23. Martišek, D. The 2D and 3D processing of images provided by conventional microscopes. Scanning 2002 ,24, 284–296. [CrossRef] [PubMed] Mathematics 2021,9, 2253 20 of 20 24. Ficker, T.; Martišek, D. Digital fracture surfaces and their roughness analysis: Applications to cement-based materials. Cem. Concreate Res. 2012,42, 827–833. [CrossRef] 25. Ficker, T. Sectional techniques for 3D imaging of microscopic and macroscopic objects. Optik 2017,144, 289–299. [CrossRef] 26. Martišek, D. 3D Reconstruction of the Surface Using a Standard Camera. Math. Probl. Eng. 2017,2017, 1–11. [CrossRef] 27. Martišek, D. Fast Shape-From-Focus method for 3D object reconstruction. Optik 2018,169, 16–26. [CrossRef] 28. Sorensen, H.; Jones, D.; Heideman, M.; Burrus, C. Real-valued fast Fourier transform algorithms. In Proceedings of the IEEE Transactions on Acoustics, Speech, and Signal Processing; IEEE: Grenoble, France, 1987; Volume 35, pp. 849–863. [CrossRef] 29. Darrell, T.; Wohn, K. Pyramid Based Depth from Focus. In Proceedings of the CVPR’88: The Computer Society Conference on Computer Vision and Pattern Recognition, Ann Arbor, MI, USA, 5–9 June 1988; Volume 2, pp. 504–509. 30. Subbarao, M. Direct Recovery of Depth Map 2: A New Robust Approach Technical Report 87-03; State University of New York: Stony Brook, NY, USA, 1987. 31. Ohta, T.; Sugihara, K.; Sugie, N. A Method for Image Composition Using Image Variance. Trans. IECE J66-D 2016 ,66, 1245–1246. 32. Druckmullerova, H. Phase-Correlation Based Image Registration. Master’s Thesis, Brno University of Technology, Brno-stˇred, Czech Republic, 2010. 33. Gillespie, J.; King, R. The use of self-entropy as a focus measure in digital holography. Patt. Rec. Lett. 1989,9, 19–25. [CrossRef] 34. Brenner, J.F.; Dew, B.S.; Horton, J.B.; King, T.; Neurath, P.W.; Selles, W.D. An automated microscope for cytologic research a Preliminary evaluation. J. Histochem. Cytochem. 1976,24, 100–111. [CrossRef] 35. Pieper, R.J.; Korpel, A. Image processing for extended depth of field. Appl. Opt. 1983,22, 1449–1453. [CrossRef] [PubMed] 36. Sugimoto, S.A.; Ichioka, Y. Digital composition of images with increased depth of focus considering depth information. Appl. Opt. 1985,24, 2076–2080. [CrossRef] [PubMed] 37. Tenebaum, J.M. Accomodation in Computer Vision. Ph.D. Thesis, Stanford University, Stanford, CA, USA, 1970. 38. Jarvis, R.A. Focus optimization criteria for computer image processing. Microscope 1976,24, 163–180. 39. Schlag, J.F.; Sanderson, A.C.; Neumann, C.P.; Wimberly, F.C. Implementation of Automatic Focusing Algorithms for a Computer Vision System with Camera Control; CMU-RI-TR-83-14; Carnegie Mellon University: Pittsburgh, PA, USA, 1983. 40. Krotkov, E. Focusing. Int. J. Comput. Vis. 1987,1, 223–237. [CrossRef] 41. Krotkov, E. Exploratory Visual Sensing with an Agile Camera. Ph.D. Thesis, TR-87-29. University of Pennsylvania, Philadelphia, PA, USA, 1987. 42. Pentland, A. A new sense for depth of field. IJCAI 1985,PAMI-9, 988–994. [CrossRef] [PubMed] 43. Grossmann, P. Depth from focus. Pattern Recognit. Lett. 1987,5, 63–69. [CrossRef] 44. Kaneda, K.; Wakasu, Y.; Nakamae, E.; Tazawa, E. A method of pan-focused and stereoscopic display using a series of optical microscopic images. Proc. of Fourth Sym. 1988, 189–194. 45. Martišek, D.; Procházková, J. The analysis of Rock Surface Asperities. Mendel 2018,24, 135–142. [CrossRef]