3D Reconstruction of the Surface using a Standard Camera

Research Article 3D Reconstruction of the Surface Using a Standard Camera Dalibor Martišek Faculty of Mechanical Engineering, Institute of Materials Science and Engineering, NETME Center, Brno University of Technology, Technick´ a 2, 616 69 Brno, Czech Republic Correspondence should be addressed to Dalibor Martiˇ sek; martis[email protected].cz Received 5 August 2017; Accepted 10 October 2017; Published 29 November 2017 Academic Editor: Yakov Strelniker Copyright © 2017 Dalibor Martiˇ sek. 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 2Dand3D model. 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. 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 howtoobtainthe3D data. These techniques have also been tested in some rock engineering projects, such as 3D digital fracture mapping [11–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–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– 22]). 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 themicroscopeisverysmallanditsadvancedhardwareis capableofremovingnonsharppointsfromtheimages.The pointsoftheobjectclosetothefocalplanearevisibleassharp points. The parts lying above or beneath the focal plane are Hindawi Mathematical Problems in Engineering Volume 2017, Article ID 4642397, 11 pages https://doi.org/10.1155/2017/4642397 2Mathematical Problems in Engineering 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.5cm. Locality Nedvˇ edice, Czech Republic, photo Pavel ˇ Starha. 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 fortechnicalpurposesduetothesmallsizeofthevisual field (maximal visual field is approximately 2 cm [5, 20–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. 2.2. Present Methods. Probably the first attempts to carry out nonconfocal reconstructions come from the seventies and eighties of the last century [23–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. 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 𝐶𝑖;𝑗;ℎ−1; 𝐶𝑖;𝑗;ℎ+1:inverseproportionality[29,30],parabolicfits[31–33], 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) isnotuseablefornoisedatafrompointofviewofnumerical 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 3Dreconstructionmethod 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 𝑓(𝑥): R→Cis function F(𝑓)(𝑥)=𝐹(𝜉)=∫∞ −∞ 𝑓(𝑥)𝑒−𝑖𝑥𝜉d𝑥(1) (if this integral exists and is finite). Standard (continuous) Fourier transform of function 𝑓(𝑥;𝑦):R2→Cis function F(𝑓)(𝑥;𝑦)=𝐹(𝜉;𝜂)=∬R2𝑓(𝑥;𝑦)𝑒−𝑖(𝑥𝜉+𝑦𝜂)d𝑥d𝑦(2) (if this integral exists and is finite). Inverse (continuous) Fourier transform of function 𝐺(𝜉): C→Cis function F−1 (𝐺)(𝜉)=𝑔(𝑥)=1 2𝜋∫∞ −∞ 𝐺(𝜉)𝑒𝑖𝑥𝜉d𝜉(3) (if this integral exists and is finite). Inverse (continuous) Fourier transform of function 𝐺(𝜉;𝜂):C2→Cis function F−1 (𝐺)(𝜉;𝜂)=𝑔(𝑥;𝑦) =1 4𝜋2∬R2𝐺(𝜉;𝜂)𝑒𝑖(𝑥𝜉+𝑦𝜂)d𝜉d𝜂(4) (if this integral exists and is finite). Mathematical Problems in Engineering 3 Standard discrete Fourier transform of function 𝑓(𝑥;𝑦): {0;1;...;𝑁−1}2→Cis function F(𝑓)(𝜉;𝜂)=𝐹(𝜉;𝜂) =𝑁−1 ∑ 𝑥=0 𝑁−1 ∑ 𝑦=0𝑓(𝑥;𝑦)𝑒−(2𝜋𝑖/𝑁)(𝑥𝜉+𝑦𝜂).(5) Function 𝐹is also called the Fourier spectrum of function 𝑓.Itispossibletoobtainthefunction𝑓from its Fourier spectrum 𝐹(𝜉;𝜂)using inverse discrete Fourier transform F−1 (𝐹)(𝑥;𝑦)=𝑓(𝑥;𝑦) =1 𝑁2 𝑁−1 ∑ 𝜉=0 𝑁−1 ∑ 𝜂=0𝐹(𝜉;𝜂)𝑒(2𝜋𝑖/𝑁)(𝑥𝜉+𝑦𝜂) (6) (see [38] for proof). Function 𝐴(𝜉;𝜂) = |𝐹(𝜉;𝜂)|is called amplitude spectrum of 𝑓(𝑥;𝑦). 2.3.2. 𝛿-Distribution. One-dimensional 𝛿-distribution 𝛿(𝑥)is alimitofasequenceoffunctions𝛿𝑛(𝑥);𝑛∈Nsuch that (a)lim 𝑛→∞ ∫∞ −∞ 𝛿𝑛(𝑥)d𝑥=1; (b)lim 𝑛→∞ 𝛿𝑛(𝑥0) lim𝑥→0𝛿𝑛(𝑥)=0; 𝑥0∈R−{0}.(7) Two-dimensional 𝛿-distribution 𝛿(𝑥;𝑦)is a limit of a sequence of functions 𝛿𝑛(𝑥;𝑦);𝑛∈Nsuch that lim 𝑛→∞ ∫∞ −∞ ∫∞ −∞ 𝛿𝑛(𝑥;𝑦)d𝑥d𝑦=1; lim 𝑛→∞ 𝛿𝑛(𝑥0;𝑦0) lim(𝑥;𝑦)→(0;0)𝛿𝑛(𝑥;𝑦)=0; (𝑥0;𝑦0)∈R2−{(0;0)}. (8) 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 (−𝑛;𝑛);𝑛∈Nandzeroedelsewhere.InverseFourier transform gives F−1 (𝛿∗ 𝑛)(𝜉)=𝛿𝑛(𝑥)=1 2𝜋∫∞ −∞ 𝑒𝑖𝑥𝜉d𝜉 =1 2𝜋∫𝑛 −𝑛 𝑒𝑖𝑥𝜉d𝜉=1 2𝜋[𝑒𝑖𝑥𝜉 𝑖𝑥]𝑛 𝜉=−𝑛 =𝑒𝑖𝑥𝑛 −𝑒−𝑖𝑥𝑛 2𝜋𝑖𝑥 =sin 𝑛𝑥 𝜋𝑥 . (9)  1 −5 5  ∗ 5() y ℱ−1 (∗ 5)()= 5(x)=ＭＣＨ 5x x x Figure 2: Fourier transforms of the fifth member of the series of expanding rectangular signals. The series converges to 𝛿distribution. Because ∫∞ −∞((sin 𝑛𝑥)/𝜋𝑥)d𝑥=1for each 𝑛∈N(see [36] for proof), condition (a) is fulfilled. Because lim𝑥→0((sin 𝑛𝑥)/ 𝜋𝑥)=(𝑛/𝜋)lim𝑥→0((sin 𝑛𝑥)/𝑛𝑥)=𝑛/𝜋for each 𝑛∈N,there is lim 𝑛→∞ 𝛿𝑛(𝑥0) lim𝑥→0𝛿𝑛(𝑥)=lim 𝑛→∞ sin 𝑛𝑥0 𝑛𝑥0=0 for each 𝑥0=0.(10) 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. 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 [39–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 𝑃𝑓1;𝑓2(𝑥;𝑦)=F−1 {𝐹1(𝜉;𝜂)⋅𝐹2(𝜉;𝜂) 󵄨󵄨󵄨󵄨𝐹1(𝜉;𝜂)󵄨󵄨󵄨󵄨⋅󵄨󵄨󵄨󵄨𝐹2(𝜉;𝜂)󵄨󵄨󵄨󵄨}(11) and its modification as 𝑃𝑓1;𝑓2;;𝑝;𝑞 (𝑥;𝑦)=F−1 {𝐻(𝜉;𝜂) ⋅𝐹1(𝜉;𝜂)⋅𝐹2(𝜉;𝜂) (󵄨󵄨󵄨󵄨𝐹1(𝜉;𝜂)󵄨󵄨󵄨󵄨+𝑝)⋅(󵄨󵄨󵄨󵄨𝐹2(𝜉;𝜂)󵄨󵄨󵄨󵄨+𝑞)}, (12) 4Mathematical Problems in Engineering where strip means complex conjugation and 𝐻(𝜉;𝜂)is a bounded real function such that 𝐻(𝜉;𝜂) = 𝐻(−𝜉;−𝜂)and 𝑝;𝑞>0are arbitrary constants. It can be proved that for real functions 𝑓1;𝑓2the 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 (𝜉;𝜂), 𝑎+𝑏𝑖 =0.Itisevidentthatthevalueofthephasecorrelation of the 𝐹with itself is 𝐹(𝜉;𝜂)⋅𝐹(𝜉;𝜂) 󵄨󵄨󵄨󵄨𝐹(𝜉;𝜂)⋅𝐹(𝜉;𝜂)󵄨󵄨󵄨󵄨=(𝑎+𝑏𝑖)(𝑎−𝑏𝑖) |(𝑎+𝑏𝑖)|⋅|𝑎+𝑏𝑖|=𝑎2+𝑏2 𝑎2+𝑏2 =1. (13) According to Example 1 in Section 2.3.2, one has 𝑃𝑓;𝑓 (𝑥;𝑦)=F−1 {𝐹(𝜉;𝜂)⋅𝐹(𝜉;𝜂) 󵄨󵄨󵄨󵄨𝐹(𝜉;𝜂)󵄨󵄨󵄨󵄨⋅󵄨󵄨󵄨󵄨𝐹(𝜉;𝜂)󵄨󵄨󵄨󵄨}=F−1 {1} =𝛿(𝑥;𝑦). (14) 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, 𝐹2(𝜉;𝜂)=𝐹1(𝜉;𝜂)⋅exp (−𝑖(𝜉𝑥0+𝜂𝑦0)),(15) and their phase-correlation function is the 𝛿-distribution shifted in arguments by the opposite shift vector 𝑃𝑓1;𝑓2(𝑥;𝑦)=F−1 {exp (𝑖(𝜉𝑥0+𝜂𝑦0))} =𝛿(𝑥+𝑥0;𝑦+𝑦0). (16) 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 𝑓2be function 𝑓1rotated and shifted in arguments; that is, 𝑓2(𝑥;𝑦) =𝑓1(𝑥cos 𝜃−𝑦sin 𝜃−𝑥0;𝑥sin 𝜃+𝑦cos 𝜃−𝑦0). (17) Their Fourier spectra 𝐹1;𝐹2and amplitude spectra 𝐴1;𝐴2are related as follows: 𝐹2(𝜉;𝜂)=exp (−𝑖(𝜉𝑥0+𝜂𝑦0)) ⋅𝐹1(𝜉cos 𝜃−𝜂sin 𝜃;𝜉sin 𝜃+𝜂cos 𝜃), 𝐴2(𝜉;𝜂)=𝐴1(𝜉cos 𝜃−𝜂sin 𝜃;𝜉sin 𝜃+𝜂cos 𝜃). (18) Theshiftresultsinaphaseshiftandthespectraarerotated 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 𝐴𝑝 1;𝐴𝑝 2:R+ 0× ⟨0;2𝜋)→R+ 0such that 𝐴𝑝 1(𝜌;𝜑)=𝐴𝑝 2(𝜌;𝜑+𝜃).Therotation around an unknown center of rotation was transformed toashift.Thisshiftisestimatedwiththestandardphase 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𝑓2be function 𝑓1rotated, shifted, and scaled in arguments; that is, 𝑓2(𝑥;𝑦)=𝑓1(𝛼(𝑥cos 𝜃−𝑦sin 𝜃) −𝑥0;𝛼(𝑥sin 𝜃+𝑦cos 𝜃)−𝑦0). (19) Their Fourier spectra and amplitude spectra are related as follows: 𝐹2(𝜉;𝜂)= 1 𝛼2exp (−𝑖(𝜉𝑥0+𝜂𝑦0)) ⋅𝐹1(1 𝛼(𝜉cos 𝜃−𝜂sin 𝜃);1 𝛼(𝜉sin 𝜃+𝜂cos 𝜃)), 𝐴2(𝜉;𝜂)= 1 𝛼2 ⋅𝐴1(1 𝛼(𝜉cos 𝜃−𝜂sin 𝜃);1 𝛼(𝜉sin 𝜃+𝜂cos 𝜃)). (20) 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 𝐴𝑝 1;𝐴𝑝 2:R+ 0×⟨0;2𝜋)→R+ 0such that 𝐴1𝑝 2(𝜌;𝜑)= 𝐴1𝑝 2(𝜌−ln 𝛼;𝜑+𝜃). Both rotation and scale change were transformed to a shift. The unknown angle 𝜃and unknown factor 𝛼can be estimatedbymeansofthephasecorrelationappliedonthe amplitude spectra in the logarithmic-polar coordinate system 𝐴1𝑝 1;𝐴1𝑝 2. After rotating function 𝑓2back by the estimated angle 𝜃and scaling by factor 𝛼,theshiftvector(𝑥0;𝑦0). 2.5. Proposed Method of Surface Approximation 2.5.1. Preprocessing: Image Registration. Before 3Dprofile calculation, it is necessary to identify all geometric Mathematical Problems in Engineering 5 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 imageedgesandoneonalargepartoftheimage.Wecanprimarily use the Gaussian window function and the Hanning window function. Let 𝜎∈R+be a given number and set R=⟨−𝑎;𝑎⟩×⟨−𝑏;𝑏⟩;𝑎;𝑏∈R+ 0, C={(𝑥;𝑦)∈R2|𝑥2+𝑦2≤𝑟}; 𝑟∈R+ 0.(22) Let 𝜌(𝑋;S)be the distance of point 𝑋=(𝑥;𝑦)from set S; that is, 𝜌(𝑋;S)=inf {𝑑∈R|𝑑=‖𝑋𝑌‖;𝑌∈S}.(23) Functions 𝑔GR (𝑥;𝑦)=𝑒−𝜌2(𝑋;R)/𝜎2 or 𝑔GC (𝑥;𝑦)=𝑒−𝜌2(𝑋;C)/𝜎2(24) are called rectangular or circular Gaussian window function. Functions 𝑔HR (𝑥;𝑦) ={ { {12+12cos 𝜋𝜌(𝑋;R) 𝜎if (𝑋;R)≤𝜎 0if (𝑋;R)>𝜎 or 𝑔HC (𝑥;𝑦) ={ { {12+12cos 𝜋𝜌(𝑋;C) 𝜎if (𝑋;C)≤𝜎 0if (𝑋;C)>𝜎 (25) are called rectangular or circular Hanning window function. Let {𝐼1;𝐼2;...;𝐼𝑛}be the image series to be registered; image 𝐼1wasacquiredbymeansofthebiggestangleofview. This image will be not transformed or (formally) it will be transformed by identity mapping to the image 𝐼∗ 1.Nowwe must find transform 𝐼2→𝐼∗ 1to obtain image 𝐼∗ 2which differs from 𝐼∗ 1infocussedandblurredpartsonly.Inthesameway, 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 Figure 3: Even extension of the 32×32neighborhood (framed) of the processing pixel (cross). 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 𝐼𝑘.Onlyshiftand 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. 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 𝑠=2 𝑛;thatis,𝑠=8or 𝑠=16is usedinthiscase.However,standardFouriertransformsuffers 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: 𝑎𝐶𝑖𝑗𝑘 =𝑠 ∑ 𝑚=−𝑠 𝑠 ∑ 𝑛=−𝑠 𝑎𝑐𝑖+𝑚;𝑗+𝑛;𝑘 =𝑠 ∑ 𝑚=−𝑠 𝑠 ∑ 𝑛=−𝑠 (|𝑚|+|𝑛|)⋅󵄨󵄨󵄨󵄨󵄨𝐹𝑖+𝑚;𝑗+𝑛;𝑘󵄨󵄨󵄨󵄨󵄨,(26) 6Mathematical Problems in Engineering (a) (b) (c) Figure 4: Graphical representation of sharpness detectors 𝑎𝐶(a); 𝑏𝐶(b); and 𝑐𝐶(c). (a) (b) 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.5cm (b). 𝑏𝐶𝑖𝑗𝑘 =𝑠 ∑ 𝑚=−𝑠 𝑠 ∑ 𝑛=−𝑠 𝑏𝑐𝑖+𝑚;𝑗+𝑛;𝑘 =1 S(A𝑖𝑗)∑ (𝑚;𝑛)∈A𝑖𝑗 󵄨󵄨󵄨󵄨𝐹𝑚𝑛𝑘󵄨󵄨󵄨󵄨,(27) 𝑐𝐶𝑖𝑗𝑘 =𝑠 ∑ 𝑚=−𝑠 𝑠 ∑ 𝑛=−𝑠 𝑐𝑐𝑖+𝑚;𝑗+𝑛;𝑘 =𝑠 ∑ 𝑚=−𝑠 𝑠 ∑ 𝑛=−𝑠 󵄨󵄨󵄨󵄨󵄨𝐹𝑖+𝑚;𝑗+𝑛;𝑘󵄨󵄨󵄨󵄨󵄨⋅sin2(𝜋𝑠√𝑚2+𝑛2), (28) 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×16or 32×32pixels. 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 quitefalse.Wehaveavailabletheseries{𝐶𝑖𝑗𝑘};𝑘=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. Oneofthepossibilitiesisaregressionanalysisbutitwouldbe very complicated. Direct calculation of the expected value is much easier. For each pixel (𝑖;𝑗),wecanconstructtheoreticallyinfinitely many probability distribution functions 𝑝(𝑟) 𝑖𝑗 using different exponents 𝑟appliedtoseriesmembers𝐶𝑖𝑗𝑘: 𝑝(𝑟) 𝑖𝑗 (𝑘)=𝐶𝑟 𝑖𝑗𝑘 ∑𝑛 𝑠=1 𝐶𝑟 𝑖𝑗𝑠 .(29) Expected values of random variables 𝑃(𝑟) 𝑖𝑗 given by these probability functions estimate the height ℎ(𝑟) 𝑖𝑗 of surface in its pixel (𝑖;𝑗): ℎ(𝑟) 𝑖𝑗 =𝐸(𝑃(𝑟) 𝑖𝑗 )=𝑛 ∑ 𝑘=1𝑘⋅𝑝(𝑟) 𝑖𝑗 (𝑘)=𝑛 ∑ 𝑘=1𝑘⋅ 𝐶𝑟 𝑖𝑗𝑘 ∑𝑛 𝑠=1 𝐶𝑟 𝑖𝑗𝑠 .(30) 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 bluegreen image on Figure 5(a) and orange image on Figure 5(b) would be “perfectly grey.” Arithmetic mean of these images Mathematical Problems in Engineering 7 (a) (b) Figure 6: The arithmetic mean of the images from Figure 5: before registration (a) and after registration (b). Real surface size 10×7.5cm (b). 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 bluegreenimagearebiggerinblue-violetpartsofthemean. 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 themeandiffersinblurredandsharpparts. 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 pixelontheimageperiphery(useddataresolutionis1600× 1200pixels). 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(rightbottomcornerofsquaresinFigure4) toomuch;itmeansthatitisverynoisesensitive.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: RMSE (ℎ(1);ℎ(2))=√1 𝑊𝐻⋅𝑊−1 ∑ 𝑖=0 𝐻−1 ∑ 𝑗=0 (ℎ(1) 𝑖𝑗 −ℎ(2) 𝑖𝑗 )2.(31) Average Deviation is as follows: AD (ℎ(1);ℎ(2))= 1 𝑊𝐻⋅𝑊−1 ∑ 𝑖=0 𝐻−1 ∑ 𝑗=0 󵄨󵄨󵄨󵄨󵄨ℎ(1) 𝑖𝑗 −ℎ(2) 𝑖𝑗 󵄨󵄨󵄨󵄨󵄨.(32) 0.00 0.00 10.00 20.00 30.00 40.0050.00 60.00 70.00 80.00 90.00 10.0020.00 30.00 40.0050.00 60.0070.0080.0090.00 100.00 (mm) (mm) x z 7.50 15.00 22.50 30.00 37.50 45.00 52.50 60.00 67.50 (mm) y 7.50 15.00 22.50 30.00 37.50 45.00 52.50 60.00 67.50 (mm) 15.00 7.50 0.00 (mm) 15.00 7.50 0.00 (mm) 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). Pearson Correlation Coefficient is as follows: PCC (ℎ(1);ℎ(2)) =∑𝑊−1 𝑖=0 ∑𝐻−1 𝑗=1 (ℎ(1) 𝑖𝑗 −ℎ(1))(ℎ(2) 𝑖𝑗 −ℎ(2)) √∑𝑊−1 𝑖=0 ∑𝐻−1 𝑗=0 (ℎ(1) 𝑖𝑗 −ℎ(1))2∑𝑊−1 𝑖=0 ∑𝐻−1 𝑗=0 (ℎ(2) 𝑖𝑗 −ℎ(2))2 (33) and Difference of surface Information Entropy is as follows: DIE (ℎ(1);ℎ(2))=󵄨󵄨󵄨󵄨󵄨IE (ℎ(1))−IE (ℎ(2))󵄨󵄨󵄨󵄨󵄨,(34) where 𝑊;𝐻are width and height of surface domain in pixels; ℎ(1) 𝑖𝑗 ;ℎ(2) 𝑖𝑗 are height of first (second) surface in the pixel (𝑖;𝑗); ℎ(1);ℎ(2) are the average height of the first (second) surface; and IE (ℎ) =−𝑊−1 ∑ 𝑖=1 𝐻−1 ∑ 𝑗=1 (ℎ𝑖𝑗 ∑𝑊−1 𝑚=1 ∑𝐻−1 𝑛=1 ℎ𝑚𝑛 log2ℎ𝑖𝑗 ∑𝑊−1 𝑚=1 ∑𝐻−1 𝑛=1 ℎ𝑚𝑛 )(35) is the Information Entropy of the surface ℎ. 8Mathematical Problems in Engineering Table 1: Parameters of the transforms indicated and applicated to separate images in the series (stated relative to the first image). Image number Indicated transforms Applicated transforms Scale Rotation (arcmin.) Shift vector Scale Rotation (arcmin.) Shift vector 𝑥(pixels) 𝑦(pixels) 𝑥(pixels) 𝑦(pixels) (2) 1.01343 −0.554 −1.011 1.064 0.98675 0.554 1.011 −1.064 (3) 1.02351 −0.221 −0.990 1.861 0.97703 0.221 0.990 −1.861 (4) 1.03796 −0.061 −1.869 3.027 0.96343 0.061 1.869 −3.027 (5) 1.04888 0.053 −2.903 3.032 0.95340 −0.053 2.903 −3.032 (6) 1.06228 −0.409 −4.085 5.142 0.94137 0.409 4.085 −5.142 (7) 1.07055 −0.140 −4.947 4.987 0.93410 0.140 4.947 −4.987 (8) 1.08105 −0.027 −4.964 5.856 0.92503 0.027 4.964 −5.856 (9) 1.09173 0.134 −4.847 6.141 0.91598 −0.134 4.847 −6.141 (10) 1.10340 4.652 −4.988 7.003 0.90629 −4.652 4.988 −7.003 (11) 1.11426 5.215 −5.003 7.934 0.89746 −5.215 5.003 −7.934 (12) 1.12590 5.728 −5.895 8.972 0.88818 −5.728 5.895 −8.972 (13) 1.13784 5.278 −5.907 8.872 0.87886 −5.278 5.907 −8.872 (14) 1.14940 5.378 −5.928 8.980 0.87002 −5.378 5.928 −8.980 (15) 1.16066 5.275 −6.935 9.137 0.86158 −5.275 6.935 −9.137 (16) 1.17297 5.324 −8.036 9.095 0.85254 −5.324 8.036 −9.095 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 𝑐𝐶. PCC (DL) RMSE (mm) DIE (DL) AD (mm) 𝑎𝐶𝑏𝐶𝑐𝐶𝑎𝐶𝑏𝐶𝑐𝐶 𝑎𝐶\0.51228 0.69630 𝑎𝐶\0.15258 0.30171 𝑏𝐶0.99320 \0.56214 𝑏𝐶0.00505 \0.13145 𝑐𝐶0.98740 0.99158 \𝑐𝐶0.00812 0.00307 \ 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 𝑐𝐶. PCC (DL) RMSE (mm) DIE (DL) AD (mm) 𝑎𝐶𝑏𝐶𝑐𝐶𝑎𝐶𝑏𝐶𝑐𝐶 𝑎𝐶\0.27750 0.38644 𝑎𝐶\0.15066 0.27180 𝑏𝐶0.99850 \0.24776 𝑏𝐶0.00464 \0.18216 𝑐𝐶0.99627 0.99810 \𝑐𝐶0.00714 0.00250 \ Values of these characteristics are summarized in Table 2. Note that it is RMSE =AD =DIE =0;PCC=1for a pair of identical surfaces. 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. Mathematical Problems in Engineering 9 Table 4: 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, Gaussian interpolation, and the proposed method applied to the same matrix 𝑐𝐶. The expected values of columns {𝑐𝐶𝑖𝑗𝑘};𝑘=1;2;...;𝑛was calculated according to expression (31) where 𝑟=5was used. PCC (DL) RMSE (mm) DIE (DL) AD (mm) Parabolic Gaussian Mean Parabolic Gaussian Mean Parabolic \0.00592 0.49098 Parabolic \0.00477 0.34745 Gaussian 1.00000 \0.48993 Gaussian 0.00002 \0.34679 Mean 0.99377 0.99379 \Mean 0.00309 0.00468 \ Figure 8: 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 𝑠=16wasusedincalculationofthematrixofsharpnessdetectors, its columns {𝑐𝐶𝑖𝑗𝑘};𝑘=1;2;...;𝑛was interpolated by the parabolic interpolation (present method). Results illustrated in Figures 8 and 9 have been obtained by combinationofpresentandproposedmethod(theyarenot realized by present methods only due to different scaling of the input data). Visualization of the parabolic interpolation of 𝑐𝐶isshowninFigure8;Gaussianinterpolationisillustrated 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. 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 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 𝑠=16was used in calculation of the matrix of sharpness detectors, its columns {𝑐𝐶𝑖𝑗𝑘};𝑘=1;2;...;𝑛was interpolated by the Gaussian interpolation (present method). Figure 10: Gaussian low-pass filter with dispersion 𝜎2=5applied to noisily output data from Figure 9. The noise is indeed reduced (but still visible); however useful high-frequency information is already lost. 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.