# 3D Reconstruction of the Surface Using a Standard Camera

Authors: Dalibor Martišek
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

This paper deals with3D 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 3D 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.


---

## Body

## 1. Introduction

The three-dimensional reconstruction of general surfaces plays an important role in many branches; for example, the morphological analysis of fracture surfaces reveals information on mechanical properties of natural or construction materials.There are more techniques capable of producing digital three-dimensional (3D) replicas of solid surfaces. Mechanical engineers can use contacting electronic profilometers to determine digital two-dimensional (2D) profiles that can be combined into 3D surface profiles; see [1, 2], for example. The contacting mode of atomic force microscopes actually belongs to this mechanical category [3]. Besides the mechanical tools, there exist different optical devices [2], light section microscopy [4, 5], coherence scanning interferometry [6], speckle metrology [7], stereo projection [8], photogrammetry [9], and various kinds of light profilometry [10], to mention some of them.3D laser scanning techniques are the next possibilities of how to obtain the3D 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 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–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 Figure1). Different image scaling (including possibly shift and rotation) may be also quantified and corrected. This problem has not been solved in available sources.Figure 1
Different scaling and different sharp and nonsharp regions in images acquired by classic camera placed in different distances from the3D 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.


### 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 Cijk is stated. Its element Cijk determines statistical range, variance, or standard Fourier transform of certain neighborhood of the pixel i;j in kth image in the series of n images [29].The Second Step. Maxima maxk⁡Cijk in the columns Cij1;Cij2;…;Cijn are found. Height hij is assigned to the pixel i;j if and only if the maximum is detected on the hth 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 Ci;j;h-1; Ci;j;h+1: inverse proportionality [29, 30], parabolic fits [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) 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 the3D reconstruction method which is able to process a series of noisy images acquired using central projection with variable center, that is, in various scaling (various shift and rotation eventually).


### 2.3. The Fourier Transform and Phase Correlation

#### 2.3.1. Standard Fourier Transform and Inverse Transform

Standard (continuous) Fourier transform of functionfx:R→C is function(1)Ffx=Fξ=∫-∞∞fxe-ixξdx (if this integral exists and is finite).Standard (continuous) Fourier transform of functionfx;y:R2→C is function(2)Ffx;y=Fξ;η=∬R2fx;ye-ixξ+yηdxdy (if this integral exists and is finite).Inverse (continuous) Fourier transform of functionGξ:C→C is function(3)F-1Gξ=gx=12π∫-∞∞Gξeixξdξ (if this integral exists and is finite).Inverse (continuous) Fourier transform of functionGξ;η:C2→C is function(4)F-1Gξ;η=gx;y=14π2∬R2Gξ;ηeixξ+yηdξdη (if this integral exists and is finite).Standard discrete Fourier transform of functionfx;y:0;1;…;N-12→C is function (5)Ffξ;η=Fξ;η=∑x=0N-1∑y=0N-1fx;ye-2πi/Nxξ+yη.Function F is also called the Fourier spectrum of function f. It is possible to obtain the function f from its Fourier spectrum Fξ;η using inverse discrete Fourier transform(6)F-1Fx;y=fx;y=1N2∑ξ=0N-1∑η=0N-1Fξ;ηe2πi/Nxξ+yη (see [38] for proof). Function Aξ;η=Fξ;η is called amplitude spectrum of fx;y.


#### 2.3.2.δ-Distribution

One-dimensional
δ
-distribution
δ
x is a limit of a sequence of functions δnx; n∈N such that(7)alimn→∞⁡∫-∞∞δnxdx=1;blimn→∞δnx0limx→0⁡δnx=0;x0∈R-0.Two-dimensional
δ
-distribution
δ
x
;
yis a limit of a sequence of functions δnx;y; n∈N such that(8)limn→∞⁡∫-∞∞∫-∞∞δnx;ydxdy=1;limn→∞δnx0;y0limx;y→0;0⁡δnx;y=0;x0;y0∈R2-0;0.Example 1.
One of the well-known examples (which we demonstrate in1D for simplicity) is the series of expanding rectangular signals δn∗ξ with constant unitary intensity on -n;n; n∈N and zeroed elsewhere. Inverse Fourier transform gives(9)F-1δn∗ξ=δnx=12π∫-∞∞eixξdξ=12π∫-nneixξdξ=12πeixξixξ=-nn=eixn-e-ixn2πix=sin⁡nxπx.Because ∫-∞∞((sin⁡nx)/πx)dx=1 for each n∈N (see [36] for proof), condition (a) is fulfilled. Because limx→0⁡((sin⁡nx)/πx)=n/πlimx→0⁡((sin⁡nx)/nx)=n/π for each n∈N, there is (10)limn→∞δnx0limx→0⁡δnx=limn→∞sin⁡nx0nx0=0foreachx0≠0.It means that also condition (b) holds and limit δx=limn→∞⁡δnx is the (one-dimensional) δ-distribution. Expanded unitary signal and its Fourier transform are illustrated in Figure 2.Figure 2
Fourier transforms of the fifth member of the series of expanding rectangular signals. The series converges toδ-distribution.


#### 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 f1; f2, it is defined as (11)Pf1;f2x;y=F-1F1ξ;η·F-2ξ;ηF1ξ;η·F2ξ;ηand its modification as(12)Pf1;f2;;p;qx;y=F-1Hξ;η·F1ξ;η·F-2ξ;ηF1ξ;η+p·F2ξ;η+q,where strip means complex conjugation and Hξ;η is a bounded real function such that Hξ;η=H-ξ;-η and p;q>0 are arbitrary constants. It can be proved that for real functions f1; f2 the phase-correlation function is real [42]. This is of great value, since it enables us to search for extremes of the phase-correlation function.


### 2.4. Image Transformation

Identical Images. Let F be the infinity periodic expansion of an image. Denote a+bi as its value Fξ;η in the point ξ;η, a+bi≠0. It is evident that the value of the phase correlation of the F with itself is (13)Fξ;η·F-ξ;ηFξ;η·Fξ;η=a+bia-bia+bi·a+bi=a2+b2a2+b2=1.According to Example 1 in Section 2.3.2, one has(14)Pf;fx;y=F-1Fξ;η·F-ξ;ηFξ;η·Fξ;η=F-11=δx;y.It means that the inverse Fourier transform of the convolution of two identical images is the two-dimensional δ-distribution δx;y.Shifted Images. If two functions are shifted in arguments, that is, f2x;y=f1x-x0;y-y0, their Fourier transforms are shifted in phase; that is,(15)F2ξ;η=F1ξ;η·exp⁡-iξx0+ηy0,and their phase-correlation function is the δ-distribution shifted in arguments by the opposite shift vector(16)Pf1;f2x;y=F-1expiξx0+ηy0=δx+x0;y+y0.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 f2 be function f1 rotated and shifted in arguments; that is,(17)f2x;y=f1xcos⁡θ-ysin⁡θ-x0;xsin⁡θ+ycos⁡θ-y0.Their Fourier spectra F1; F2 and amplitude spectra A1; A2 are related as follows:(18)F2ξ;η=exp-iξx0+ηy0·F1ξcos⁡θ-ηsin⁡θ;ξsin⁡θ+ηcos⁡θ,A2ξ;η=A1ξcos⁡θ-ηsin⁡θ;ξsin⁡θ+ηcos⁡θ.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 A1p;A2p:R0+×0;2π→R0+ such that A1pρ;φ=A2pρ;φ+θ. 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 x0;y0 is measured with another computation of the phase correlation.Scaled Images. Let f2 be function f1 rotated, shifted, and scaled in arguments; that is,(19)f2x;y=f1αxcos⁡θ-ysin⁡θ-x0;αxsin⁡θ+ycos⁡θ-y0.Their Fourier spectra and amplitude spectra are related as follows:(20)F2ξ;η=1α2exp-iξx0+ηy0·F11αξcos⁡θ-ηsin⁡θ;1αξsin⁡θ+ηcos⁡θ,A2ξ;η=1α2A11αξcos⁡θ-ηsin⁡θ;1αξsin⁡θ+ηcos⁡θ.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(21)exp⁡ρ=x2+y2;x=exp⁡ρcos⁡φ;y=exp⁡ρsin⁡φto obtain A1p;A2p:R0+×0;2π→R0+ such that A21pρ;φ=A21pρ-ln⁡α;φ+θ.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 A11p; A21p. After rotating function f2 back by the estimated angle θ and scaling by factor α, the shift vector x0;y0.


### 2.5. Proposed Method of Surface Approximation

#### 2.5.1. Preprocessing: Image Registration

Before3D 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(22)R=-a;a×-b;b;a;b∈R0+,C=x;y∈R2∣x2+y2≤r;r∈R0+.Let ρX;S be the distance of point X=x;y from set S; that is,(23)ρX;S=inf⁡d∈R∣d=XY;Y∈S.Functions(24)gGRx;y=e-ρ2X;R/σ2orgGCx;y=e-ρ2X;C/σ2are called rectangular or circular Gaussian window function. Functions(25)gHRx;y=12+12cos⁡πρX;RσifX;R≤σ0ifX;R>σorgHCx;y=12+12cos⁡πρX;CσifX;C≤σ0ifX;C>σare called rectangular or circular Hanning window function.LetI1;I2;…;In be the image series to be registered; image I1 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 I1∗. Now we must find transform I2→I1∗ to obtain image I2∗ which differs from I1∗ in focussed and blurred parts only. In the same way, transforms I3→I2∗;…;Ik→Ik-1∗;…;In→In-1∗ must be found.After the multiplying both imagesIk; Ik-1∗ by the chosen window function, we determine the rotation angle θk and scale factor αk by means of the method described in Section 2.4. Then, image Ik is rotated by angle -θk and scaled by factor αk-1 to compensate the rotation and scale change found by the phase correlation, creating image Ik¯. Only shift and different focussed and blurred parts remain between image Ik¯ and image Ik-1∗. Now we can apply phase correlation to find the shift x0;y0 and image Ik¯ is shifted by vector -x0;-y0 to compensate the shift, creating image Ik∗ which differs from Ik-1∗ in focussed and blurred parts only.


#### 2.5.2. The First Step: Matrix of Sharpness Detectors

As was said in Section2.2, statistical range, variance, or standard Fourier transform of certain neighborhood of the pixel i;j may be used as the sharpness detectors (focusing criteria). The neighborhood of the pixel i;j is a square of s×s pixels where s is obviously ten to twenty. In case of the standard Fourier transform, algorithm FFT is used. It requires the square with the side s=2n; that is, s=8 or s=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.Figure 3
Even extension of the32×32 neighborhood (framed) of the processing pixel (cross).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)Caijk=∑m=-ss∑n=-sscai+m;j+n;k=∑m=-ss∑n=-ssm+n·Fi+m;j+n;k,(27)Cbijk=∑m=-ss∑n=-sscbi+m;j+n;k=1SAij∑m;n∈AijFmnk,(28)Ccijk=∑m=-ss∑n=-sscci+m;j+n;k=∑m=-ss∑n=-ssFi+m;j+n;k·sin2⁡πsm2+n2,where F is the cosine spectrum of the pixel neighborhood and SAij is volume of the annulus Aij with the center i;j. Elements ca; cb; cc 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.Figure 4
Graphical representation of sharpness detectorsCa (a); Cb (b); and Cc (c).
(a)
(b)
(c)


#### 2.5.3. The Second Step: Profile Heights Calculation

The main imperfection of current methods is the hypothesis that the profile height in given pixel is exactly determined by values of chosen sharpness detector. This hypothesis implies that these values can be interpolated. However, this conclusion is quite false. We have available the seriesCijk; k=1;2;…;n for assessment of the height of pixel i;j. 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 pixeli;j, we can construct theoretically infinitely many probability distribution functions pij(r) using different exponents r applied to series members Cijk:(29)pijrk=Cijkr∑s=1nCijsr.Expected values of random variables Pij(r) given by these probability functions estimate the height hij(r) of surface in its pixel i;j:(30)hijr=EPijr=∑k=1nk·pijrk=∑k=1nk·Cijkr∑s=1nCijsr.


## 2.1. Equipment

In technical practice the confocal microscope is used as a standard tool for imaging microscopic three-dimensional surfaces. The depth of the optical field of the microscope is very small and its advanced hardware is capable of removing nonsharp points from the images. The points of the object close to the focal plane are visible as sharp points. The parts lying above or beneath the focal plane are invisible and represented as black regions. To create a 2D and 3D reconstruction, it is necessary to obtain a series of images of the same object, each of them with different focusing, and each point of the object focussed in one of the images (in the ideal case). The sharp parts are identified and composed of a 2D and 3D model.However, the confocal microscope is often not suitable for technical purposes due to the small size of the visual field (maximal visual field is approximately 2 cm [5, 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 Figure1). Different image scaling (including possibly shift and rotation) may be also quantified and corrected. This problem has not been solved in available sources.Figure 1
Different scaling and different sharp and nonsharp regions in images acquired by classic camera placed in different distances from the3D 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.


## 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 Cijk is stated. Its element Cijk determines statistical range, variance, or standard Fourier transform of certain neighborhood of the pixel i;j in kth image in the series of n images [29].The Second Step. Maxima maxk⁡Cijk in the columns Cij1;Cij2;…;Cijn are found. Height hij is assigned to the pixel i;j if and only if the maximum is detected on the hth 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 Ci;j;h-1; Ci;j;h+1: inverse proportionality [29, 30], parabolic fits [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) 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 the3D reconstruction method which is able to process a series of noisy images acquired using central projection with variable center, that is, in various scaling (various shift and rotation eventually).


## 2.3. The Fourier Transform and Phase Correlation

### 2.3.1. Standard Fourier Transform and Inverse Transform

Standard (continuous) Fourier transform of functionfx:R→C is function(1)Ffx=Fξ=∫-∞∞fxe-ixξdx (if this integral exists and is finite).Standard (continuous) Fourier transform of functionfx;y:R2→C is function(2)Ffx;y=Fξ;η=∬R2fx;ye-ixξ+yηdxdy (if this integral exists and is finite).Inverse (continuous) Fourier transform of functionGξ:C→C is function(3)F-1Gξ=gx=12π∫-∞∞Gξeixξdξ (if this integral exists and is finite).Inverse (continuous) Fourier transform of functionGξ;η:C2→C is function(4)F-1Gξ;η=gx;y=14π2∬R2Gξ;ηeixξ+yηdξdη (if this integral exists and is finite).Standard discrete Fourier transform of functionfx;y:0;1;…;N-12→C is function (5)Ffξ;η=Fξ;η=∑x=0N-1∑y=0N-1fx;ye-2πi/Nxξ+yη.Function F is also called the Fourier spectrum of function f. It is possible to obtain the function f from its Fourier spectrum Fξ;η using inverse discrete Fourier transform(6)F-1Fx;y=fx;y=1N2∑ξ=0N-1∑η=0N-1Fξ;ηe2πi/Nxξ+yη (see [38] for proof). Function Aξ;η=Fξ;η is called amplitude spectrum of fx;y.


### 2.3.2.δ-Distribution

One-dimensional
δ
-distribution
δ
x is a limit of a sequence of functions δnx; n∈N such that(7)alimn→∞⁡∫-∞∞δnxdx=1;blimn→∞δnx0limx→0⁡δnx=0;x0∈R-0.Two-dimensional
δ
-distribution
δ
x
;
yis a limit of a sequence of functions δnx;y; n∈N such that(8)limn→∞⁡∫-∞∞∫-∞∞δnx;ydxdy=1;limn→∞δnx0;y0limx;y→0;0⁡δnx;y=0;x0;y0∈R2-0;0.Example 1.
One of the well-known examples (which we demonstrate in1D for simplicity) is the series of expanding rectangular signals δn∗ξ with constant unitary intensity on -n;n; n∈N and zeroed elsewhere. Inverse Fourier transform gives(9)F-1δn∗ξ=δnx=12π∫-∞∞eixξdξ=12π∫-nneixξdξ=12πeixξixξ=-nn=eixn-e-ixn2πix=sin⁡nxπx.Because ∫-∞∞((sin⁡nx)/πx)dx=1 for each n∈N (see [36] for proof), condition (a) is fulfilled. Because limx→0⁡((sin⁡nx)/πx)=n/πlimx→0⁡((sin⁡nx)/nx)=n/π for each n∈N, there is (10)limn→∞δnx0limx→0⁡δnx=limn→∞sin⁡nx0nx0=0foreachx0≠0.It means that also condition (b) holds and limit δx=limn→∞⁡δnx is the (one-dimensional) δ-distribution. Expanded unitary signal and its Fourier transform are illustrated in Figure 2.Figure 2
Fourier transforms of the fifth member of the series of expanding rectangular signals. The series converges toδ-distribution.


### 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 f1; f2, it is defined as (11)Pf1;f2x;y=F-1F1ξ;η·F-2ξ;ηF1ξ;η·F2ξ;ηand its modification as(12)Pf1;f2;;p;qx;y=F-1Hξ;η·F1ξ;η·F-2ξ;ηF1ξ;η+p·F2ξ;η+q,where strip means complex conjugation and Hξ;η is a bounded real function such that Hξ;η=H-ξ;-η and p;q>0 are arbitrary constants. It can be proved that for real functions f1; f2 the phase-correlation function is real [42]. This is of great value, since it enables us to search for extremes of the phase-correlation function.


## 2.3.1. Standard Fourier Transform and Inverse Transform

Standard (continuous) Fourier transform of functionfx:R→C is function(1)Ffx=Fξ=∫-∞∞fxe-ixξdx (if this integral exists and is finite).Standard (continuous) Fourier transform of functionfx;y:R2→C is function(2)Ffx;y=Fξ;η=∬R2fx;ye-ixξ+yηdxdy (if this integral exists and is finite).Inverse (continuous) Fourier transform of functionGξ:C→C is function(3)F-1Gξ=gx=12π∫-∞∞Gξeixξdξ (if this integral exists and is finite).Inverse (continuous) Fourier transform of functionGξ;η:C2→C is function(4)F-1Gξ;η=gx;y=14π2∬R2Gξ;ηeixξ+yηdξdη (if this integral exists and is finite).Standard discrete Fourier transform of functionfx;y:0;1;…;N-12→C is function (5)Ffξ;η=Fξ;η=∑x=0N-1∑y=0N-1fx;ye-2πi/Nxξ+yη.Function F is also called the Fourier spectrum of function f. It is possible to obtain the function f from its Fourier spectrum Fξ;η using inverse discrete Fourier transform(6)F-1Fx;y=fx;y=1N2∑ξ=0N-1∑η=0N-1Fξ;ηe2πi/Nxξ+yη (see [38] for proof). Function Aξ;η=Fξ;η is called amplitude spectrum of fx;y.


## 2.3.2.δ-Distribution

One-dimensional
δ
-distribution
δ
x is a limit of a sequence of functions δnx; n∈N such that(7)alimn→∞⁡∫-∞∞δnxdx=1;blimn→∞δnx0limx→0⁡δnx=0;x0∈R-0.Two-dimensional
δ
-distribution
δ
x
;
yis a limit of a sequence of functions δnx;y; n∈N such that(8)limn→∞⁡∫-∞∞∫-∞∞δnx;ydxdy=1;limn→∞δnx0;y0limx;y→0;0⁡δnx;y=0;x0;y0∈R2-0;0.Example 1.
One of the well-known examples (which we demonstrate in1D for simplicity) is the series of expanding rectangular signals δn∗ξ with constant unitary intensity on -n;n; n∈N and zeroed elsewhere. Inverse Fourier transform gives(9)F-1δn∗ξ=δnx=12π∫-∞∞eixξdξ=12π∫-nneixξdξ=12πeixξixξ=-nn=eixn-e-ixn2πix=sin⁡nxπx.Because ∫-∞∞((sin⁡nx)/πx)dx=1 for each n∈N (see [36] for proof), condition (a) is fulfilled. Because limx→0⁡((sin⁡nx)/πx)=n/πlimx→0⁡((sin⁡nx)/nx)=n/π for each n∈N, there is (10)limn→∞δnx0limx→0⁡δnx=limn→∞sin⁡nx0nx0=0foreachx0≠0.It means that also condition (b) holds and limit δx=limn→∞⁡δnx is the (one-dimensional) δ-distribution. Expanded unitary signal and its Fourier transform are illustrated in Figure 2.Figure 2
Fourier transforms of the fifth member of the series of expanding rectangular signals. The series converges toδ-distribution.


## 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 f1; f2, it is defined as (11)Pf1;f2x;y=F-1F1ξ;η·F-2ξ;ηF1ξ;η·F2ξ;ηand its modification as(12)Pf1;f2;;p;qx;y=F-1Hξ;η·F1ξ;η·F-2ξ;ηF1ξ;η+p·F2ξ;η+q,where strip means complex conjugation and Hξ;η is a bounded real function such that Hξ;η=H-ξ;-η and p;q>0 are arbitrary constants. It can be proved that for real functions f1; f2 the phase-correlation function is real [42]. This is of great value, since it enables us to search for extremes of the phase-correlation function.


## 2.4. Image Transformation

Identical Images. Let F be the infinity periodic expansion of an image. Denote a+bi as its value Fξ;η in the point ξ;η, a+bi≠0. It is evident that the value of the phase correlation of the F with itself is (13)Fξ;η·F-ξ;ηFξ;η·Fξ;η=a+bia-bia+bi·a+bi=a2+b2a2+b2=1.According to Example 1 in Section 2.3.2, one has(14)Pf;fx;y=F-1Fξ;η·F-ξ;ηFξ;η·Fξ;η=F-11=δx;y.It means that the inverse Fourier transform of the convolution of two identical images is the two-dimensional δ-distribution δx;y.Shifted Images. If two functions are shifted in arguments, that is, f2x;y=f1x-x0;y-y0, their Fourier transforms are shifted in phase; that is,(15)F2ξ;η=F1ξ;η·exp⁡-iξx0+ηy0,and their phase-correlation function is the δ-distribution shifted in arguments by the opposite shift vector(16)Pf1;f2x;y=F-1expiξx0+ηy0=δx+x0;y+y0.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 f2 be function f1 rotated and shifted in arguments; that is,(17)f2x;y=f1xcos⁡θ-ysin⁡θ-x0;xsin⁡θ+ycos⁡θ-y0.Their Fourier spectra F1; F2 and amplitude spectra A1; A2 are related as follows:(18)F2ξ;η=exp-iξx0+ηy0·F1ξcos⁡θ-ηsin⁡θ;ξsin⁡θ+ηcos⁡θ,A2ξ;η=A1ξcos⁡θ-ηsin⁡θ;ξsin⁡θ+ηcos⁡θ.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 A1p;A2p:R0+×0;2π→R0+ such that A1pρ;φ=A2pρ;φ+θ. 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 x0;y0 is measured with another computation of the phase correlation.Scaled Images. Let f2 be function f1 rotated, shifted, and scaled in arguments; that is,(19)f2x;y=f1αxcos⁡θ-ysin⁡θ-x0;αxsin⁡θ+ycos⁡θ-y0.Their Fourier spectra and amplitude spectra are related as follows:(20)F2ξ;η=1α2exp-iξx0+ηy0·F11αξcos⁡θ-ηsin⁡θ;1αξsin⁡θ+ηcos⁡θ,A2ξ;η=1α2A11αξcos⁡θ-ηsin⁡θ;1αξsin⁡θ+ηcos⁡θ.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(21)exp⁡ρ=x2+y2;x=exp⁡ρcos⁡φ;y=exp⁡ρsin⁡φto obtain A1p;A2p:R0+×0;2π→R0+ such that A21pρ;φ=A21pρ-ln⁡α;φ+θ.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 A11p; A21p. After rotating function f2 back by the estimated angle θ and scaling by factor α, the shift vector x0;y0.


## 2.5. Proposed Method of Surface Approximation

### 2.5.1. Preprocessing: Image Registration

Before3D 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(22)R=-a;a×-b;b;a;b∈R0+,C=x;y∈R2∣x2+y2≤r;r∈R0+.Let ρX;S be the distance of point X=x;y from set S; that is,(23)ρX;S=inf⁡d∈R∣d=XY;Y∈S.Functions(24)gGRx;y=e-ρ2X;R/σ2orgGCx;y=e-ρ2X;C/σ2are called rectangular or circular Gaussian window function. Functions(25)gHRx;y=12+12cos⁡πρX;RσifX;R≤σ0ifX;R>σorgHCx;y=12+12cos⁡πρX;CσifX;C≤σ0ifX;C>σare called rectangular or circular Hanning window function.LetI1;I2;…;In be the image series to be registered; image I1 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 I1∗. Now we must find transform I2→I1∗ to obtain image I2∗ which differs from I1∗ in focussed and blurred parts only. In the same way, transforms I3→I2∗;…;Ik→Ik-1∗;…;In→In-1∗ must be found.After the multiplying both imagesIk; Ik-1∗ by the chosen window function, we determine the rotation angle θk and scale factor αk by means of the method described in Section 2.4. Then, image Ik is rotated by angle -θk and scaled by factor αk-1 to compensate the rotation and scale change found by the phase correlation, creating image Ik¯. Only shift and different focussed and blurred parts remain between image Ik¯ and image Ik-1∗. Now we can apply phase correlation to find the shift x0;y0 and image Ik¯ is shifted by vector -x0;-y0 to compensate the shift, creating image Ik∗ which differs from Ik-1∗ in focussed and blurred parts only.


### 2.5.2. The First Step: Matrix of Sharpness Detectors

As was said in Section2.2, statistical range, variance, or standard Fourier transform of certain neighborhood of the pixel i;j may be used as the sharpness detectors (focusing criteria). The neighborhood of the pixel i;j is a square of s×s pixels where s is obviously ten to twenty. In case of the standard Fourier transform, algorithm FFT is used. It requires the square with the side s=2n; that is, s=8 or s=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.Figure 3
Even extension of the32×32 neighborhood (framed) of the processing pixel (cross).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)Caijk=∑m=-ss∑n=-sscai+m;j+n;k=∑m=-ss∑n=-ssm+n·Fi+m;j+n;k,(27)Cbijk=∑m=-ss∑n=-sscbi+m;j+n;k=1SAij∑m;n∈AijFmnk,(28)Ccijk=∑m=-ss∑n=-sscci+m;j+n;k=∑m=-ss∑n=-ssFi+m;j+n;k·sin2⁡πsm2+n2,where F is the cosine spectrum of the pixel neighborhood and SAij is volume of the annulus Aij with the center i;j. Elements ca; cb; cc 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.Figure 4
Graphical representation of sharpness detectorsCa (a); Cb (b); and Cc (c).
(a)
(b)
(c)


### 2.5.3. The Second Step: Profile Heights Calculation

The main imperfection of current methods is the hypothesis that the profile height in given pixel is exactly determined by values of chosen sharpness detector. This hypothesis implies that these values can be interpolated. However, this conclusion is quite false. We have available the seriesCijk; k=1;2;…;n for assessment of the height of pixel i;j. 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 pixeli;j, we can construct theoretically infinitely many probability distribution functions pij(r) using different exponents r applied to series members Cijk:(29)pijrk=Cijkr∑s=1nCijsr.Expected values of random variables Pij(r) given by these probability functions estimate the height hij(r) of surface in its pixel i;j:(30)hijr=EPijr=∑k=1nk·pijrk=∑k=1nk·Cijkr∑s=1nCijsr.


## 2.5.1. Preprocessing: Image Registration

Before3D 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(22)R=-a;a×-b;b;a;b∈R0+,C=x;y∈R2∣x2+y2≤r;r∈R0+.Let ρX;S be the distance of point X=x;y from set S; that is,(23)ρX;S=inf⁡d∈R∣d=XY;Y∈S.Functions(24)gGRx;y=e-ρ2X;R/σ2orgGCx;y=e-ρ2X;C/σ2are called rectangular or circular Gaussian window function. Functions(25)gHRx;y=12+12cos⁡πρX;RσifX;R≤σ0ifX;R>σorgHCx;y=12+12cos⁡πρX;CσifX;C≤σ0ifX;C>σare called rectangular or circular Hanning window function.LetI1;I2;…;In be the image series to be registered; image I1 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 I1∗. Now we must find transform I2→I1∗ to obtain image I2∗ which differs from I1∗ in focussed and blurred parts only. In the same way, transforms I3→I2∗;…;Ik→Ik-1∗;…;In→In-1∗ must be found.After the multiplying both imagesIk; Ik-1∗ by the chosen window function, we determine the rotation angle θk and scale factor αk by means of the method described in Section 2.4. Then, image Ik is rotated by angle -θk and scaled by factor αk-1 to compensate the rotation and scale change found by the phase correlation, creating image Ik¯. Only shift and different focussed and blurred parts remain between image Ik¯ and image Ik-1∗. Now we can apply phase correlation to find the shift x0;y0 and image Ik¯ is shifted by vector -x0;-y0 to compensate the shift, creating image Ik∗ which differs from Ik-1∗ in focussed and blurred parts only.


## 2.5.2. The First Step: Matrix of Sharpness Detectors

As was said in Section2.2, statistical range, variance, or standard Fourier transform of certain neighborhood of the pixel i;j may be used as the sharpness detectors (focusing criteria). The neighborhood of the pixel i;j is a square of s×s pixels where s is obviously ten to twenty. In case of the standard Fourier transform, algorithm FFT is used. It requires the square with the side s=2n; that is, s=8 or s=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.Figure 3
Even extension of the32×32 neighborhood (framed) of the processing pixel (cross).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)Caijk=∑m=-ss∑n=-sscai+m;j+n;k=∑m=-ss∑n=-ssm+n·Fi+m;j+n;k,(27)Cbijk=∑m=-ss∑n=-sscbi+m;j+n;k=1SAij∑m;n∈AijFmnk,(28)Ccijk=∑m=-ss∑n=-sscci+m;j+n;k=∑m=-ss∑n=-ssFi+m;j+n;k·sin2⁡πsm2+n2,where F is the cosine spectrum of the pixel neighborhood and SAij is volume of the annulus Aij with the center i;j. Elements ca; cb; cc 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.Figure 4
Graphical representation of sharpness detectorsCa (a); Cb (b); and Cc (c).
(a)
(b)
(c)


## 2.5.3. The Second Step: Profile Heights Calculation

The main imperfection of current methods is the hypothesis that the profile height in given pixel is exactly determined by values of chosen sharpness detector. This hypothesis implies that these values can be interpolated. However, this conclusion is quite false. We have available the seriesCijk; k=1;2;…;n for assessment of the height of pixel i;j. 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 pixeli;j, we can construct theoretically infinitely many probability distribution functions pij(r) using different exponents r applied to series members Cijk:(29)pijrk=Cijkr∑s=1nCijsr.Expected values of random variables Pij(r) given by these probability functions estimate the height hij(r) of surface in its pixel i;j:(30)hijr=EPijr=∑k=1nk·pijrk=∑k=1nk·Cijkr∑s=1nCijsr.


## 3. Results and Discussion

### 3.1. The Preprocessing: Image Registration

In Figure5, we can see the photos from Figure 1 but they are displayed in so-called supplementary pseudocolors. If these images were completely identical, then arithmetic mean of the blue-green image on Figure 5(a) and orange image on Figure 5(b) would be “perfectly grey.” Arithmetic mean of these images is constructed in Figure 6(a). It is clear that components of this mean are very different, values of the orange image are bigger in yellow parts of the mean, and values of the blue-green image are bigger in blue-violet parts of the mean.Figure 5
The first and the fifteenth images from the series of fifteen images of blue marble (see Figure1) displayed in supplementary pseudocolors. Real surface size 10×7.5 cm (b).
(a)
(b)Figure 6
The arithmetic mean of the images from Figure5: before registration (a) and after registration (b). Real surface size 10×7.5 cm (b).
(a)
(b)In Figure6(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 Table1, 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.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
x (pixels)
y (pixels)
x (pixels)
y (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


### 3.2. The First Step: Sharpness Detectors

As is clear from (26), (27), and (28) and Figure 4, detector Ca 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 Cb are too sharp cuts which remove low and height frequencies.Differences between the sharpness detectors stated in Section2.5.2 are best demonstrated by their maxima. The graphical representation of detector Cc (i.e., corresponding stair-approximation of the surface) is shown in Figure 7.Figure 7
The stairs-approximation of the data from Figure1 constructed by means of the sharpness detector Cc; see expression (28) and Figure 4(c).Differences between representationsCa; Cb; and Cc are visually quite miniscule (therefore, Ca and Cb are not illustrated); however, differences exist. For their quantification, Root Mean Square Error is as follows:(31)RMSEh1;h2=1WH·∑i=0W-1∑j=0H-1hij1-hij22.Average Deviation is as follows:(32)ADh1;h2=1WH·∑i=0W-1∑j=0H-1hij1-hij2.Pearson Correlation Coefficient is as follows:(33)PCCh1;h2=∑i=0W-1∑j=1H-1hij1-h1¯hij2-h2¯∑i=0W-1∑j=0H-1hij1-h1¯2∑i=0W-1∑j=0H-1hij2-h2¯2and Difference of surface Information Entropy is as follows:(34)DIEh1;h2=IEh1-IEh2,where W; H are width and height of surface domain in pixels; hij(1); hij2 are height of first (second) surface in the pixel (i;j); h(1)¯; h(2)¯ are the average height of the first (second) surface; and(35)IEh=-∑i=1W-1∑j=1H-1hij∑m=1W-1∑n=1H-1hmnlog2hij∑m=1W-1∑n=1H-1hmnis the Information Entropy of the surface h.Values of these characteristics are summarized in Table2. Note that it is RMSE=AD=DIE=0; PCC=1 for a pair of identical surfaces.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 Ca, Cb, and Cc.
PCC (DL)
RMSE (mm)
DIE (DL)
AD (mm)
C
a
C
b
C
c
C
a
C
b
C
c
C
a
∖
0.51228
0.69630
C
a
∖
0.15258
0.30171
C
b
0.99320
∖
0.56214
C
b
0.00505
∖
0.13145
C
c
0.98740
0.99158
∖
C
c
0.00812
0.00307
∖Similar summary is made in Table3 for the present reconstruction methods: parabolic interpolation of the stair-approximations Ca, Cb, and Cc. In Table 4, there are summarized RMSE, AD, PCC, and DIE for parabolic and Gaussian interpolation of the same stair-approximation Cc and proposed method of the expected values of columns Ccijk; k=1;2;…;n.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 Ca, Cb, and Cc.
PCC (DL)
RMSE (mm)
DIE (DL)
AD (mm)
C
a
C
b
C
c
C
a
C
b
C
c
C
a
∖
0.27750
0.38644
C
a
∖
0.15066
0.27180
C
b
0.99850
∖
0.24776
C
b
0.00464
∖
0.18216
C
c
0.99627
0.99810
∖
C
c
0.00714
0.00250
∖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 Cc. The expected values of columns Ccijk; k=1;2;…;n was calculated according to expression (31) where r=5 was 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
∖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 Cc is shown in Figure 8; Gaussian interpolation is illustrated in Figure 9. Practically the same and large noise is evident in these reconstructions.Figure 8
3
D 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 s=16 was used in calculation of the matrix of sharpness detectors, its columns Ccijk; k=1;2;…;n was interpolated by the parabolic interpolation (present method).Figure 9
3
D 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 s=16 was used in calculation of the matrix of sharpness detectors, its columns Ccijk; k=1;2;…;n was interpolated by the Gaussian interpolation (present method).Low-pass filters are commonly used to reduce the noise, but these filters are not ableto differentiate whether high-frequency information is caused by noise or by small details of useful information. Therefore, the loss of useful information is inevitable as shown in Figure10.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.In Figure11, there is illustrated significantly better result of proposed method applied to the same matrix Cc. There is no visible noise at the output, and high-frequency useful information (small surface details) is retained.Figure 11
Completely new3D reconstruction of the data from Figure 1: preprocessing according to Section 2.5.1, the first step according to Section 2.5.2—expression (28) with s=16 was used— and the second step according to Section 2.5.3, the expected values of its columns Ccijk; k=1;2;…;n was calculated according to expression (29) where r=5 was used.


## 3.1. The Preprocessing: Image Registration

In Figure5, we can see the photos from Figure 1 but they are displayed in so-called supplementary pseudocolors. If these images were completely identical, then arithmetic mean of the blue-green image on Figure 5(a) and orange image on Figure 5(b) would be “perfectly grey.” Arithmetic mean of these images is constructed in Figure 6(a). It is clear that components of this mean are very different, values of the orange image are bigger in yellow parts of the mean, and values of the blue-green image are bigger in blue-violet parts of the mean.Figure 5
The first and the fifteenth images from the series of fifteen images of blue marble (see Figure1) displayed in supplementary pseudocolors. Real surface size 10×7.5 cm (b).
(a)
(b)Figure 6
The arithmetic mean of the images from Figure5: before registration (a) and after registration (b). Real surface size 10×7.5 cm (b).
(a)
(b)In Figure6(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 Table1, 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.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
x (pixels)
y (pixels)
x (pixels)
y (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


## 3.2. The First Step: Sharpness Detectors

As is clear from (26), (27), and (28) and Figure 4, detector Ca 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 Cb are too sharp cuts which remove low and height frequencies.Differences between the sharpness detectors stated in Section2.5.2 are best demonstrated by their maxima. The graphical representation of detector Cc (i.e., corresponding stair-approximation of the surface) is shown in Figure 7.Figure 7
The stairs-approximation of the data from Figure1 constructed by means of the sharpness detector Cc; see expression (28) and Figure 4(c).Differences between representationsCa; Cb; and Cc are visually quite miniscule (therefore, Ca and Cb are not illustrated); however, differences exist. For their quantification, Root Mean Square Error is as follows:(31)RMSEh1;h2=1WH·∑i=0W-1∑j=0H-1hij1-hij22.Average Deviation is as follows:(32)ADh1;h2=1WH·∑i=0W-1∑j=0H-1hij1-hij2.Pearson Correlation Coefficient is as follows:(33)PCCh1;h2=∑i=0W-1∑j=1H-1hij1-h1¯hij2-h2¯∑i=0W-1∑j=0H-1hij1-h1¯2∑i=0W-1∑j=0H-1hij2-h2¯2and Difference of surface Information Entropy is as follows:(34)DIEh1;h2=IEh1-IEh2,where W; H are width and height of surface domain in pixels; hij(1); hij2 are height of first (second) surface in the pixel (i;j); h(1)¯; h(2)¯ are the average height of the first (second) surface; and(35)IEh=-∑i=1W-1∑j=1H-1hij∑m=1W-1∑n=1H-1hmnlog2hij∑m=1W-1∑n=1H-1hmnis the Information Entropy of the surface h.Values of these characteristics are summarized in Table2. Note that it is RMSE=AD=DIE=0; PCC=1 for a pair of identical surfaces.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 Ca, Cb, and Cc.
PCC (DL)
RMSE (mm)
DIE (DL)
AD (mm)
C
a
C
b
C
c
C
a
C
b
C
c
C
a
∖
0.51228
0.69630
C
a
∖
0.15258
0.30171
C
b
0.99320
∖
0.56214
C
b
0.00505
∖
0.13145
C
c
0.98740
0.99158
∖
C
c
0.00812
0.00307
∖Similar summary is made in Table3 for the present reconstruction methods: parabolic interpolation of the stair-approximations Ca, Cb, and Cc. In Table 4, there are summarized RMSE, AD, PCC, and DIE for parabolic and Gaussian interpolation of the same stair-approximation Cc and proposed method of the expected values of columns Ccijk; k=1;2;…;n.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 Ca, Cb, and Cc.
PCC (DL)
RMSE (mm)
DIE (DL)
AD (mm)
C
a
C
b
C
c
C
a
C
b
C
c
C
a
∖
0.27750
0.38644
C
a
∖
0.15066
0.27180
C
b
0.99850
∖
0.24776
C
b
0.00464
∖
0.18216
C
c
0.99627
0.99810
∖
C
c
0.00714
0.00250
∖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 Cc. The expected values of columns Ccijk; k=1;2;…;n was calculated according to expression (31) where r=5 was 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
∖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 Cc is shown in Figure 8; Gaussian interpolation is illustrated in Figure 9. Practically the same and large noise is evident in these reconstructions.Figure 8
3
D 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 s=16 was used in calculation of the matrix of sharpness detectors, its columns Ccijk; k=1;2;…;n was interpolated by the parabolic interpolation (present method).Figure 9
3
D 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 s=16 was used in calculation of the matrix of sharpness detectors, its columns Ccijk; k=1;2;…;n was interpolated by the Gaussian interpolation (present method).Low-pass filters are commonly used to reduce the noise, but these filters are not ableto differentiate whether high-frequency information is caused by noise or by small details of useful information. Therefore, the loss of useful information is inevitable as shown in Figure10.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.In Figure11, there is illustrated significantly better result of proposed method applied to the same matrix Cc. There is no visible noise at the output, and high-frequency useful information (small surface details) is retained.Figure 11
Completely new3D reconstruction of the data from Figure 1: preprocessing according to Section 2.5.1, the first step according to Section 2.5.2—expression (28) with s=16 was used— and the second step according to Section 2.5.3, the expected values of its columns Ccijk; k=1;2;…;n was calculated according to expression (29) where r=5 was used.


## 4. Conclusion

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



---

*Source: 4642397-2017-11-29.xml*