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

Abstract

Shape-from-Focus (SFF) methods have been developed for about twenty years. They able to obtain the shape of 3D objects from a series of partially focused images. The plane to which the microscope or camera is focused intersects the 3D object in a contour line. Due to wave properties of light and due to finite resolution of the output device, the image can be considered as sharp not only on this contour line, but also in a certain interval of height—the zone of sharpness. SSFs are able to identify these focused parts to compose a fully focused 2D image and to reconstruct a 3D profile of the surface to be observed.

1. Introduction

Three-dimensional reconstruction of general surfaces has an important role in a number of fields: the morphological analysis of fracture surfaces, for example, reveals information on the mechanical properties of natural or construction materials.

There are more techniques capable of producing digital three-dimensional (3D) replicas of solid surfaces. In mechanical engineering, contacting electronic profilometers can be used to determine digital two-dimensional (2D) profiles to be combined into 3D surface profiles—see [1,2], for example. The contacting mode of atomic force microscopes is actually in this mechanical category [3]. In addition to the mechanical tools, optical devices exist in diverse modifications [2], light section microscopy [4,5], coherence scanning interferometry [6], speckle metrology [7], stereo projection [8], photogrammetry [9], and various types of light measurement of profiles [10], to mention some of them.

3D laser scanning techniques are among other ways of obtaining 3D data. They have also been tested in some rock engineering projects, such as 3D digital fracture mapping [11,12,13].

These devices are, however, not of universal use, with each of them having its own technical limitations [14,15]. Very rough surfaces, for instance, can hardly be measured by atomic force microscopes working in the nano-regions. On the other hand, it is possible to measure plane surfaces with microscopically small irregularities using a microscopic sectional technique with confocal microscopes [16,17,18,19]. Confocal microscopes, however, are not always suitable for technical purposes due to the small size of their visual fields (maximal visual field is about 2 cm [5,20,21,22]).

The present paper summarized the existing methods of 3D reconstruction of objects by the Shape-From-Focus (SFF) method. This is a method for recovering depth from an image series of the same object taken with different focus settings, referred to as a multifocal image.

It consists of these steps:

(a)

Data acquisition (confocal microscope in the standard mode, CCD camera, or standard camera)—Section 2.

(b)

Image registration (if necessary)—Section 3.

(c)

Choice of the focusing criteria—Section 4.

(d)

The 2D and 3D reconstructions—Section 5.

See [23] for one of the first papers on this subject.

2. Data Acquisition

2.1. Parallel Projection

In technical practice, the confocal microscope serves as a standard instrument imaging microscopic three-dimensional surfaces. It has a very small depth of the optical field with its advanced hardware being capable of removing non-sharp points from the images. The points of the object situated close to the focal plane can be seen as sharp points. The parts lying further above or beneath the focal plane (out of the sharpness zone) are invisible, being represented as black regions if the confocal mode is on. In this way, a so-called optical cut is obtained. In the case of the non-confocal (standard) mode, the areas lying outside the sharpness zone are displayed as blurred, as they would be with a standard camera. With a confocal microscope or CCD camera, one can assume that the field of view is small with the projection used being parallel. In this case, all images are provided in the field of view with identical sizes and the corresponding pixels having the same coordinates in separate partial focused images (see Figure 1).

However, the confocal microscope with its small visual angle is hardly suitable for technical purposes due to the small size of its visual field (maximal visual field is about 2 cm [5,20,21,24]).

2.2. Central Projection

The same output (Figure 1b) can be obtained by the classical microscope or (in a wider field) common camera. The difference between the microscope or CCD camera and standard camera is given by the central projection that varies the scaling of partial images in a series and, further, by the non-sharp regions, displayed by the classic camera, while missing when taken by a confocal microscope in the confocal mode (see [25] for more information). The different image scalings, however, require subsequent corrections (including shifts and rotations).

2.3. Multifocal Image

To create a 2D or 3D reconstruction, it is necessary to obtain a series of images of an identical object, each of them with different focusing with each object point being focused in one of the images (in the ideal case, this is referred to as a multifocal image). For the acquisition of a large multifocal image, the camera must be mounted on a stand so that it can be moved in the direction approximately orthogonal to the surface with a controlled step. Different transformations and different sharp parts must be identified and composed in a 2D or 3D model.

3. Image Registration

With a confocal microscope or CCD camera, we can assume that the field of view is small and the projection used is parallel. In this case, all images are provided in the field of view with identical sizes with the corresponding pixels having the identical coordinates in separate partial focused images. However, this assumption does not hold for larger samples; then, the angle of the projection lines is not negligible with the view fields (and the coordinates of the corresponding pixels) being clearly different for each image (see Figure 2 and Figure 3).

Before reconstruction, all geometric transformations must be identified in the image series to be eliminated. The images are being analyzed as geometrically similar. Generally, similarity is achieved by composing rotation, scale-change, shift, and axial symmetry. Axial symmetry is not possible in this particular case.

If we consider scale-change only assuming that the image size is proportional to the camera shifting (see Figure 4a), different image scaling can be obtained using elementary mathematics. This approach was used in [25]. After the initial elementary reconstruction, a subsequent 3D reconstruction is shown in Figure 5. Huge artefacts caused by inaccurate registration were blurred by brutal low-pass filters with much useful high-frequency information lost.

In practice, the situation may be more sophisticated. The images may differ not just in the scale used but in the content displayed as well (different parts being focused in different images). Due to mechanical inaccuracies, the step in the z axis may be not fully constant, and the images can also be mutually shifted along the x- or y-axis or rotated. Image registration is also complicated by the non-planarity of samples (see Figure 4 on the right). Therefore, sophisticated pre-processing of the image series may be necessary. A method suitable tool for this is the Fourier transform and phase correlation.

3.1. Continuous Two-Dimensional Fourier Transform and Inverse Transform

A continuous standard Fourier transform of a function $f\left(x\right):ℝ\to ℂ$ is the function

$$\mathcal{F}\left[f\left(x\right)\right]=F\left(\xi \right)=\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}f\left(x\right)e^{-ix\xi }\mathrm{d}x$$

(1)

(provided that this integral exists and is finite).

A continuous standard Fourier transform of a function $f\left(x;y\right):{ℝ}^2\to ℂ$ is the function

$$\mathcal{F}\left[f\left(x;y\right)\right]=F\left(\xi ;\eta \right)=\underset{{ℝ}^2}{\overset{ }{{\displaystyle \iint }}}f\left(x;y\right)e^{-i\left(x\xi +y\eta \right)}\mathrm{d}x\mathrm{d}y$$

(2)

(provided that this integral exists and is finite).

Function $F$ is also referred to as the Fourier spectrum of function $f.$ It is possible to obtain the function $f$ from its Fourier spectrum $F$ by the inverse Fourier transform.

A continuous standard inverse Fourier transform of a function $G\left(\xi \right):ℂ\to ℂ$ is the function

$${\mathcal{F}}^{-1}\left[F\left(\xi \right)\right]=f\left(x\right)=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}F\left(\xi \right)e^{ix\xi }\mathrm{d}\xi$$

(3)

(provided that the integral exists and is finite).

A continuous standard inverse Fourier transform of a function $G\left(\xi ;\eta \right):{ℂ}^2\to ℂ$ is the function

$${\mathcal{F}}^{-1}\left[F\left(\xi ;\eta \right)\right]=f\left(x;y\right)=\frac{1}{4{\pi }^2}{{\displaystyle \iint }}_{{ℝ}^2}F\left(\xi ;\eta \right)e^{i\left(x\xi +y\eta \right)}\mathrm{d}\xi \mathrm{d}\eta$$

(4)

(provided that this integral exists and is finite).

3.2. Discrete Two-Dimensional Fourier Transform and Inverse Transform

A discrete standard Fourier transform of a function $f\left(x;y\right):$ $\left\{0;1;\dots ;M-1\right\}\times \left\{0;1;\dots ;N-1\right\}\to ℂ$ is the function

$$\mathcal{F}\left[f\left(x;y\right)\right]=F\left(\xi ;\eta \right)={\displaystyle \sum }_{x=0}^{M-1}{\displaystyle \sum }_{y=0}^{N-1}f\left(x;y\right)e^{-2\pi i\left(\frac{x\xi }{M}+\frac{y\eta }{N}\right)}$$

(5)

A discrete standard inverse Fourier transform of a function $G\left(\xi ;\eta \right):$$\left\{0;1;\dots ;M-1\right\}\times \left\{0;1;\dots ;N-1\right\}\to ℂ$ is the function

$${\mathcal{F}}^{-1}\left[G\left(\xi ;\eta \right)\right]=g\left(x;y\right)=\frac{1}{MN}{\displaystyle \sum }_{\xi =0}^{M-1}{\displaystyle \sum }_{\eta =0}^{N-1}G\left(\xi ;\eta \right)e^{2\pi i\left(\frac{x\xi }{M}+\frac{y\eta }{N}\right)}$$

(6)

If the function $f\left(x;y\right)$ in (5) is real, then

$$F\left(N-\xi ;N-\eta \right)=\overline{F}\left(\xi ;\eta \right)$$

(7)

(where the bar denotes complex conjugation) and the Fast Fourier Transform (FFT) algorithm can be employed to calculate the discrete Fourier transform—see [28] for more information.

The discrete Fourier transform can be used for the image registration applied to functions that are or are assumed periodic. Generally, an image may not have the same values on the edges. Thus, by periodizing an image, the resulting function may have jumps at the edges of the original image. Such jumps are often the most contrasted structures in the function and may lead to incorrect registration. Therefore, such edges used for the shift estimation must be removed from the images. This is carried out by multiplying the image by a suitable function referred to as a window function. Its values must equal zero or almost zero at the image edges and one on a large part of the image. Primarily, the Gaussian and Hanning window functions can be used: Let $\sigma \in {ℝ}^{+}$ be a given number and sets

$$\mathcal{R}=\langle -a;a\rangle \times \langle -b;b\rangle ; a;b\in {ℝ}_0^{+}$$

$$\mathcal{C}=\left\{\left(x;y\right)\in {ℝ}^2|x^2+y^2\le r\right\}; r\in {ℝ}_0^{+}$$

Let $\rho \left(X;\mathrm{S}\right)$ be the distance of point $X=\left(x;y\right)$ from set $\mathrm{S}$; i.e.,

$$\rho \left(X;\mathrm{S}\right)=\inf \left\{d\in ℝ|d=\Vert XY\Vert ;Y\in \mathrm{S}\right\}$$

(8)

Functions

$$g_{GR}\left(x;y\right)=e^{-\frac{{\rho }^2\left(X;\mathcal{R}\right)}{{\sigma }^2}} \mathrm{or} g_{GC}\left(x;y\right)=e^{-\frac{{\rho }^2\left(X;\mathcal{C}\right)}{{\sigma }^2}}$$

(9)

are called rectangular or circular Gaussian window functions. Functions

$$g_{HR}\left(x;y\right)=\{\begin{array}{c}\frac{1}{2}+\frac{1}{2}\cos \frac{\pi \rho \left(X;\mathcal{R}\right)}{\sigma } \mathrm{if} \left(X;\mathcal{R}\right)\le \sigma \\ 0 \mathrm{if} \left(X;\mathcal{R}\right)>\sigma \end{array}$$

or

$$g_{HC}\left(x;y\right)=\{\begin{array}{c}\frac{1}{2}+\frac{1}{2}\cos \frac{\pi \rho \left(X;\mathcal{C}\right)}{\sigma } \mathrm{if} \left(X;\mathcal{C}\right)\le \sigma \\ 0 \mathrm{if} \left(X;\mathcal{C}\right)>\sigma \end{array}$$

(10)

are called rectangular or circular Hanning window functions.

3.3. δ-Distribution

A one-dimensional $\delta$-distribution $\delta \left(x\right)$ is a limit of a function sequence ${\delta }_n\left(x\right);n\in \mathbb{N}$ for which

$$\left(a\right)\underset{n\to \infty }{\lim }\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}{\delta }_n\left(x\right)\mathrm{d}x=$$

$$\left(b\right)\underset{n\to \infty }{\lim }\frac{{\delta }_n\left(x_0\right)}{\underset{x\to 0}{\lim }{\delta }_n\left(x\right)}=0; x_0\in ℝ-\left\{0\right\}$$

(11)

A two-dimensional $\delta$-distribution $\delta \left(x;y\right)$ is a limit of a function sequence ${\delta }_n\left(x;y\right);n\in \mathbb{N}$ for which

$$\underset{n\to \infty }{\lim }\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}{\delta }_n\left(x;y\right)\mathrm{d}x\mathrm{d}y=1;\underset{n\to \infty }{\lim }\frac{{\delta }_n\left(x_0;y_0\right)}{\underset{\left(x;y\right)\to \left(0;0\right)}{\lim }{\delta }_n\left(x;y\right)}=0;\left(x_0;y_0\right)\in ℝ{ }^2-\left\{\left(0;0\right)\right\}$$

(12)

Example: A well-known example (here demonstrated in $1\mathrm{D}$ for simplicity) is a series of expanding rectangular signals ${\delta }_n^{*}\left(\xi \right)$ with a unitary intensity constant on $\left(-n;n\right);n\in \mathbb{N}$ and zeroed elsewhere. The inverse Fourier transform gives

$$\begin{gathered} {\mathcal{F}}^{-1}\left({\delta }_n^{*}\right)\left(\xi \right)={\delta }_n\left(x\right)=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}{\mathrm{e}}^{ix\xi }\mathrm{d}\xi =\frac{1}{2\pi }\underset{-n}{\overset{n}{{\displaystyle \int }}}{\mathrm{e}}^{ix\xi }\mathrm{d}\xi =\frac{1}{2\pi }{\left[\frac{{\mathrm{e}}^{ix\xi }}{ix}\right]}_{\xi =-n}^n \\ =\frac{{\mathrm{e}}^{ixn}-{\mathrm{e}}^{-ixn}}{2\pi ix}=\frac{\sin nx}{\pi x} \end{gathered}$$

As ${\int }_{-\infty }^{\infty }\frac{\sin nx}{\pi x}\mathrm{d}x=1$ for each $n\in \mathbb{N}$ (see [28] for proof), condition (a) is fulfilled. Since $\underset{x\to 0}{\lim }\frac{\sin nx}{\pi x}=\frac{n}{\pi }\underset{x\to 0}{\lim }\frac{\sin nx}{nx}=\frac{n}{\pi }$ for each $n\in \mathbb{N}, \mathrm{we}\mathrm{have}$

$$\underset{n\to \infty }{\lim }\frac{{\delta }_n\left(x_0\right)}{\underset{x\to 0}{\lim }{\delta }_n\left(x\right)}=\underset{n\to \infty }{\lim }\frac{\sin n{x}_0}{n{x}_0}=0 \mathrm{for} \mathrm{each} x_0\ne 0$$

.

This means that condition (b) holds as well and the limit $\delta \left(x\right)=\underset{n\to \infty }{\lim }{\delta }_n\left(x\right)$ is the (one-dimensional) $\delta$-distribution. The expanded unitary signal ${\delta }_5^{*}\left(\xi \right)$ with its Fourier transform is illustrated in Figure 6.

3.4. Phase Correlation

Image processing requires the images transformed for the structures studied to be at the same position in all of them. The transformation is found by image registration. In some applications, it is possible to assume shift only while, in others, shift, rotation and scale change (i.e., similarity), general linear transformation, or even general transformations may all be present.

The methods used for registration depend on the expected transformation and on the image structures. Some methods, after using the corresponding structures or points in the images, find a global transformation by measuring the positions of the structures or points [29,30,31]. For these methods to be applicable, the structures must be clearly visible. Other correlation-based methods work with the image as a whole. The phase correlation has proved to be a powerful tool (not only) for the registration of particular focused images. For functions $f_1; f_2$, it is defined as

$$P_{f_1;f_2}\left(x;y\right)={ℱ}^{-1}\left\{\frac{F_1\left(\xi ;\eta \right)·{\overline{F}}_2\left(\xi ;\eta \right)}{\left|F_1\left(\xi ;\eta \right)\right|·\left|F_2\left(\xi ;\eta \right)\right|}\right\}$$

(13)

with its modification being

$$P_{f_1;f_{2;};p;q}\left(x;y\right)={ℱ}^{-1}\left\{H\left(\xi ;\eta \right)·\frac{F_1\left(\xi ;\eta \right)·{\overline{F}}_2\left(\xi ;\eta \right)}{\left(\left|F_1\left(\xi ;\eta \right)\right|+p\right)·\left(\left|F_2\left(\xi ;\eta \right)\right|+q\right)}\right\}$$

(14)

where the bar denotes complex conjugation, $H\left(\xi ;\eta \right)$ is a bounded real function such that $H\left(\xi ;\eta \right)=H\left(-\xi ;-\eta \right)$ and $p;q>0$ are arbitrary constants. It is not difficult to prove that, for real functions $f_1$, $f_2$, the phase-correlation function is real [32]. This is very useful since the extremes of the phase-correlation function can be searched for.

3.5. Identical Images

Let $F$ be the infinity periodic expansion of an image. Denote $a+bi$ its value $F\left(\xi ;\eta \right)$ at $\left(\xi ;\eta \right)$, $a+bi\ne 0$. Clearly, the value of the phase correlation of the $F$ with itself is

$$\frac{F\left(\xi ;\eta \right)·\overline{F}\left(\xi ;\eta \right)}{\left|F\left(\xi ;\eta \right)·F\left(\xi ;\eta \right)\right|}=\frac{\left(a+bi\right)\left(a-bi\right)}{\left|\left(a+bi\right)\right|·\left|a+bi\right|}=\frac{a^2+b^2}{a^2+b^2}=1$$

(15)

By the example in Section 3.3, we have

$$P_{f;f}\left(x;y\right)={ℱ}^{-1}\left\{\frac{F\left(\xi ;\eta \right)·\overline{F}\left(\xi ;\eta \right)}{\left|F\left(\xi ;\eta \right)\right|·\left|F\left(\xi ;\eta \right)\right|}\right\}={ℱ}^{-1}\left\{1\right\}=\delta \left(x;y\right)$$

(16)

which means that the inverse Fourier transform of the correlation of two identical images is the two-dimensional $\delta$-distribution $\delta \left(x;y\right)$

3.6. Shifted Images

If two functions are shifted in arguments, that is, $f_2\left(x;y\right)=f_1\left(x-x_0;y-y_0\right)$, their Fourier transforms are shifted in phase:

$$F_2\left(\xi ;\eta \right)=F_1\left(\xi ;\eta \right)·\exp \left(-i\left(\xi x_0+\eta y_0\right)\right)$$

(17)

with their phase-correlation function being the δ-distribution shifted in arguments by the opposite shift vector

$$P_{f_1;f_2}\left(x;y\right)={ℱ}^{-1}\left\{\exp \left(i\left(\xi x_0+\eta y_0\right)\right)\right\}=\delta \left(x+x_0;y+y_0\right)$$

(18)

This is the principal idea of phase correlation. Using phase correlation, rather than finding a shift between two images, we can just find the only non-zero point in a matrix. If the images are not identical (up to a shift), i.e., if the images are not ideal, the phase-correlation function is more complex, but still has a global maximum at the point whose coordinates correspond to the shift vector.

3.7. Rotated Images

The phase-correlation function can also be used for estimating the image rotation and rescale. Let $f_2$ be function $f_1$ rotated and shifted in arguments, i.e.,

$$f_2\left(x;y\right)=f_1\left(x\cos \theta -y\sin \theta -x_0;x\sin \theta +y\cos \theta -y_0\right)$$

(19)

Their Fourier spectra $F_1$; $F_2$ and amplitude spectra $A_1$; $A_2$ are related as follows:

$$F_2\left(\xi ;\eta \right)=\exp \left(-i\left(\xi x_0+\eta y_0\right)\right)·F_1\left(\xi \cos \theta -\eta \sin \theta ;\xi \sin \theta +\eta \cos \theta \right)$$

(20)

$$A_2\left(\xi ;\eta \right)=A_1\left(\xi \cos \theta -\eta \sin \theta ;\xi \sin \theta +\eta \cos \theta \right)$$

(21)

The shift results in a phase shift and the spectra are rotated the same way as the original functions. A crucial step here is the transformation of the amplitude spectra into polar coordinates to obtain functions $A_1^p;A_2^p:{ℝ}_0^{+}\times 0;2\pi )\to {ℝ}_0^{+}$ such that $A_1^p\left(\rho ;\phi \right)=A_2^p\left(\rho ;\phi +\theta \right)$. The rotation about an unknown centre has been transformed into a shift. This shift is estimated by the standard phase correlation (see the previous paragraph) after a reverse rotation by the angle measured, the shift $\left(x_0;y_0\right)$ is then measured in another computation of the phase correlation.

3.8. Scaled Images

Let $f_2$ be function $f_1$ rotated, shifted, and scaled in arguments, i.e.,

$$f_2\left(x;y\right)=f_1\left(\alpha \left(x\cos \theta -y\sin \theta \right)-x_0;\alpha \left(x\sin \theta +y\cos \theta \right)-y_0\right)$$

(22)

Their Fourier spectra and amplitude spectra are related as follows:

$$F_2\left(\xi ;\eta \right)=\frac{1}{{\alpha }^2}\exp \left(-i\left(\xi x_0+\eta y_0\right)\right)·F_1\left(\frac{1}{\alpha }\left(\xi \cos \theta -\eta \sin \theta \right);\frac{1}{\alpha }\left(\xi \sin \theta +\eta \cos \theta \right)\right)$$

(23)

$$A_2\left(\xi ;\eta \right)=\frac{1}{{\alpha }^2}{A}_1\left(\frac{1}{\alpha }\left(\xi \cos \theta -\eta \sin \theta \right);\frac{1}{\alpha }\left(\xi \sin \theta +\eta \cos \theta \right)\right)$$

(24)

The shift results in a phase shift, the spectra are rotated the same way as the original functions and scaled with a reciprocal factor. A crucial step here is the transformation of the amplitude spectra into the logarithmic-polar coordinates

$$\exp \rho =\sqrt{x^2+y^2};x=\exp \rho \cos \phi ;\hspace{1em}y=\exp \rho \sin \phi$$

(25)

to obtain $A_1^p;A_2^p:{ℝ}_0^{+}\times \langle 0;2\pi )\to {ℝ}_0^{+}$ such that $A_2^{1p}\left(\rho ;\phi \right)=A_2^{1p}\left(\rho -\ln \alpha ;\phi +\theta \right)$.

Both rotation and scale change have been transformed to a shift. The unknown angle θ and unknown factor α can be estimated by means of the phase correlation applied to the amplitude spectra in the logarithmic-polar coordinates $A_1^{1p};A_2^{1p}$. After reverse rotating function $f_2$ by the estimated angle $\theta$ and scaling by the factor, the shift vector $\left(x_0;y_0\right)$ can be estimated by means of the standard phase correlation.

3.9. Multifocal Registration

Let $\left\{P_1;P_2;\dots ;P_n\right\}$ be the image series to be registered with image $P_1$ acquired by means of the biggest angle of view. This image will be not transformed or (formally) it will be transformed by the identity mapping into the image $P_1^{*}$. Now we must find the transform $P_2\to P_1^{*}$ to obtain image $P_2^{*}$ which only differs from $P_1^{*}$ in focussed and blured parts. In the same way, transforms $P_3\to P_2^{*}$;…;$P_k\to P_{k-1}^{*};$…; $P_n\to P_{n-1}^{*}$ must be found.

After multiplying both images $P_k;P_{k-1}^{*}$ by the chosen window function, rotation angle ${\theta }_k$ and scale factor ${\alpha }_k$ will be determined by the method described in Section 3.8. Then, image $P_k$ is rotated by the angle $-{\theta }_k$ and scaled by the factor ${\alpha }_k^{-1}$ to compensate for the rotation and scale-change found by the phase correlation, creating image $\overline{P_k}$. Between images $\overline{P_k}$. and $P_{k-1}^{*}$, only the shifted and different focused and blurred parts remain. Now we can apply phase correlation to find the shift $\left(x_0;y_0\right)$ shifting image $\overline{P_k}$ by the vector $\left(-x_0;-y_0\right)$ to compensate for the shift, creating image $P_k^{*}$ which only differs from $P_{k-1}^{*}$ in the focused and blurred parts.

4. Focusing Criteria

The detectors of blurred areas (sometimes referred to as focusing criteria or sharpness detectors) can be based on different principles. Probably the first attempts to carry out non-confocal reconstructions date back to 1970’s and 1980’s [33,34,35,36].

Tenebaum [37] developed the gradient magnitude maximization method for optimizing the focus quality using the sharpness of edges. Jarvis [38] proposed a sum-modulus-difference computed by summing the first intensity differences between neighbouring pixels along a scan-line using it as a focus quality benchmark. Schlag et al. [39] implemented and tested various self-focusing algorithms. Recently, Krotkov [40] evaluated and compared the performance of different focus criterion functions. In [41], he also proposed a method for estimating the depth of an image area. Pentland [42] suggested evaluating the image blur to determine the depth of image points. Grossmann in [43] proposed estimating the depth of edge points by analyzing the blur of the edges due to defocusing. Darrell and Wohn in [29] developed a depth-from-focus method by which an image sequence can be obtained through varying the focus level using Laplacian and Gaussian pyramids to calculate the depth. Subbarao in [30] suggested changing the intrinsic camera parameters to recover the depth map of a scene. Ohta et al. [31] and Kaneda et al. [44] used images corresponding to different focus levels to obtain a single level of high focus quality.

As follows from the above, the statistical range or the variance of the standard Fourier transform of a certain neighbourhood of pixel $\left[i;j\right]$ can serve as sharpness detectors (focussing criteria). A neighbourhood of pixel $\left[i;j\right]$ is a square of $s\times s$ pixels with $s$ equalling ten to twenty. In the case of the standard Fourier transform, the FFT algorithm is used. Since it requires a square with a side of $s=2^n$, $s=8$ or $s=16$ is used in this case [25,26]. However, because the standard Fourier transform is burdened with jumps along the square edges, it is necessary to use a suitable window function like in Section 3.2. For this reason, the cosine Fourier transform is preferable. It is obtained from the standard Fourier transform applied to an even extension of the neighbourhood to be processed. It is illustrated in Figure 7. This extension eliminates jumps on the edges. Since the sine frequencies in (5) are zeroed, there is no need to apply a window function.

Whereas the low frequencies in the amplitude spectrum detect the blurred parts of the image, the very high ones only mean noise. Therefore, a suitable weight must be assigned to each frequency when the sharpness detector is calculated. In our software, the following detectors may be used:

$${}_{a}P{}_{ijk}={\displaystyle \sum }_{m=-s}^s{\displaystyle \sum }_{n=-s}^s{}_{a}p{}_{i+m;j+n;k}={\displaystyle \sum }_{m=-s}^s{\displaystyle \sum }_{n=-s}^s\left(\left|m\right|+\left|n\right|\right)·\left|F_{i+m;j+n;k}\right|$$

(26)

$${}_{b}P{}_{ijk}={\displaystyle \sum }_{m=-s}^s{\displaystyle \sum }_{n=-s}^s{}_{b}p{}_{i+m;j+n;k}=\frac{1}{\mathrm{S}\left({\mathrm{A}}_{ij}\right)}{\displaystyle \sum }_{\left(m;n\right)\in {\mathrm{A}}_{ij}}\left|F_{mnk}\right|$$

(27)

$${}_{c}P{}_{ijk}={\displaystyle \sum }_{m=-s}^s{\displaystyle \sum }_{n=-s}^s{}_{c}p{}_{i+m;j+n;k}={\displaystyle \sum }_{m=-s}^s{\displaystyle \sum }_{n=-s}^s\left|F_{i+m;j+n;k}\right|·{\sin }^2\left(\frac{\pi }{s}\sqrt{m^2+n^2}\right)$$

(28)

where $F$ is the cosine spectrum of the pixel neighborhood, $\mathrm{S}\left({\mathrm{A}}_{ij}\right)$ is the volume of annulus ${\mathrm{A}}_{ij}$ with the center at $\left(i;j\right)$. The elements ${}_{a}c;{}_{b}c$; ${}_{c}c$ summed in (26)–(28) are illustrated as pixel values in Figure 8.

5. 2D and 3D Reconstructions

A 2D reconstruction consists of composing a new image that contains only the focused parts of a registered multifocal image. The registration and detection of the focused parts were described in Section 3 and Section 4.

The principal deficiency of most of the current methods of 3D reconstruction is that they assume the profile height at a point to be precisely determined by the value of the given sharpness detector. Based on such an unrealistic hypothesis, then, these values are interpolated. Parabolic interpolation [16,34] or Gaussian interpolation [35,36] is used.

This conclusion, however, is false. We can use the series $\left\{P_{ijk}\right\};k=1;2;\cdots ;n$ to assess the height of pixel $\left(i;j\right)$. Being of a random rather than deterministic nature, it cannot be interpolated, but must be processed by a statistical method. One of such methods is a regression analysis, which, however, is rather complicated. Direct calculation of the mean value is much easier.

For each pixel $\left(i;j\right),$ in k-th image in the series virtually infinitely many probability distribution functions $p_{ij}^{\left(r\right)}$ can be constructed using different exponents $r$ applied to the series terms $P_{ijk}$:

$$p_{ij}^{\left(r\right)}\left(k\right)=\frac{P_{ijk}^r}{{{\displaystyle \sum }}_{s=1}^n{P}_{ijs}^r}$$

(29)

The mean values of the random variables $P_{ij}^{\left(r\right)}$ given by these probablity distribution functions estimate the height $h_{ij}^{\left(r\right)}$ of the surface in its pixel $\left(i;j\right):$

$$h_{ij}^{\left(r\right)}=E\left(P_{ij}^{\left(r\right)}\right)={\displaystyle \sum }_{k=1}^{n}k·p_{ij}^{\left(r\right)}\left(k\right)={\displaystyle \sum }_{k=1}^{n}k·\frac{P_{ijk}^r}{{{\displaystyle \sum }}_{s=1}^{n}C{P}_{ijs}^r}$$

(30)

6. Results and Discussion

6.1. Data Acquisition

The following data were used for the purposes of this paper:

• A series of fifteen partially focused images of blue marble—Olympus camera;
• A series of eight partially focused limestone images—Canon camera;
• A series of thirty partially focused images of the surface of a hydrated cement paste —the Olympus confocal microscope in standard mode.

6.2. Image Registration

Figure 9 shows the photos of Figure 2 but this time displayed in supplementary pseudo-colours. If these images were totally identical, then the arithmetic mean of the blue-green image on the left and orange image on the right would be “perfectly gray”. The arithmetic mean of these images is shown in Figure 10a. Clearly, the components of this mean value are very different—the values of the orange image are bigger in the yellow parts of the mean value while the values of the blue-green image are bigger in the blue-violet parts of the mean value.

In Figure 10b, the same construction was used after the registration. The very low colour saturation of the arithmetic mean (of the images) testifies to a very good conformity. Of course, the arithmetic mean cannot be a ”perfect gray“ in this particular case because the original images also differ in the sharp parts.

Table 1. Parameters of the transforms indicated and applied to separate images in the series of blue marble series (stated relative to the first image). Taken from [26].

	Transforms Indicated				Transforms Applied
Image No.	Scale	Rotation	Shift Vector		Scale	Rotation	Shift Vector
		(Arcmin.)	x (Pixels)	y (Pixels)		(Arcmin.)	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

summarizes the indicated and applied transforms in separate images of blue marble (see Figure 2, Figure 9 and Figure 10). All transformations, detected with a sub-pixel precision, are listed with a precision of one thousandth of a pixel. Obviously, the scaling has the most important role; the shifts, however, cannot be neglected either. The rotation angle between the first and the last images is larger than five arcminutes, that is, about one pixel on the image periphery (the data resolution used was $1180 \mathrm{by} 885$ pixels). Such a transformation is marginal in this particular case.

6.3. Focusing Criteria

As implied by (26)–(28), and Figure 9, detector ${}_{a}C$ increased the high frequencies (right bottom corner of squares in Figure 9) too much, which means that it is very noise sensitive. Too sharp cuts removing low and high frequencies are among the disadvantages of the detector ${}_{b}C$. Although visually minuscule, the differences between detectors ${}_{a}C;$${}_{b}C$; and ${}_{c}C$ do exist.

For their quantification, Root Mean Square Error,

$$RMSE\left(h^{\left(1\right)}; h^{\left(2\right)}\right)=\sqrt{\frac{1}{WH}·{\displaystyle \sum }_{i=0}^{W-1}{\displaystyle \sum }_{j=0}^{H-1}{\left(h_{ij}^{\left(1\right)}-h_{ij}^{\left(2\right)}\right)}^2}$$

(31)

Average Deviation

$$AD\left(h^{\left(1\right)}; h^{\left(2\right)}\right)=\frac{1}{WH}·{\displaystyle \sum }_{i=0}^{W-1}{\displaystyle \sum }_{j=0}^{H-1}\left|h_{ij}^{\left(1\right)}-h_{ij}^{\left(2\right)}\right|$$

(32)

Pearson Correlation Coefficient

$$PCC\left(h^{\left(1\right)}; h^{\left(2\right)}\right)=\frac{{{\displaystyle \sum }}_{i=0}^{W-1}{{\displaystyle \sum }}_{j=1}^{H-1}\left(h_{ij}^{\left(1\right)}-\overline{h^{\left(1\right)}}\right)\left(h_{ij}^{\left(2\right)}-\overline{h^{\left(2\right)}}\right)}{\sqrt{{{\displaystyle \sum }}_{i=0}^{W-1}{{\displaystyle \sum }}_{j=0}^{H-1}{\left(h_{ij}^{\left(1\right)}-\overline{h^{\left(1\right)}}\right)}^2{{\displaystyle \sum }}_{i=0}^{W-1}{{\displaystyle \sum }}_{j=0}^{H-1}{\left(h_{ij}^{\left(2\right)}-\overline{h^{\left(2\right)}}\right)}^2}}$$

(33)

and Difference of surface Information Entropy

$$DIE\left(h^{\left(1\right)}; h^{\left(2\right)}\right)=\left|IE\left(h^{\left(1\right)}\right)-IE\left(h^{\left(2\right)}\right)\right|$$

(34)

are used. Here, $W; H$ denote the width and height of the surface domain in pixels, $h_{ij}^{\left(1\right)}$; $h_{ij}^{\left(2\right)}$ stand for the height of the first (second) surface in the pixel $\left(i;j\right)$; $\overline{h^{\left(1\right)}}$; $\overline{h^{\left(2\right)}}$ is the average height of the first (second) surface and

$$IE\left(h\right)=-{\displaystyle \sum }_{i=1}^{W-1}{\displaystyle \sum }_{j=1}^{H-1}\left(\frac{h_{ij}}{{{\displaystyle \sum }}_{m=1}^{W-1}{{\displaystyle \sum }}_{n=1}^{H-1}{h}_{mn}}{\log }_2\frac{h_{ij}}{{{\displaystyle \sum }}_{m=1}^{W-1}{{\displaystyle \sum }}_{n=1}^{H-1}{h}_{mn}}\right)$$

(35)

is the Information Entropy of the surface $h$.

Values of these characteristics are summarized in

Table 2. Root Mean Square Error ($RMSE)$, Average Deviation ($AD$), Pearson Correlation Coefficient ($PCC$) and and Difference of surface Information Entropy ($DIE$) for separate pairs of echelon−approximations ${}_{a}C$, ${}_{b}C$, and ${}_{c}C$ for the data of Figure 3. Taken from [26].

	RMSE	aC	bC	cC		AD	aC	bC	cC
PCC					DIE
aC		-	0.34152	0.46420	aC		-	0.10172	0.20114
bC		0.99320	-	0.37476	bC		0.00505	-	0.13145
cC		0.98740	0.99158	-	cC		0.00812	0.00307	-

for data from Figure 3. Note that we have $RMSE=AD=DIE=0$; $PCC=1$ for a pair of identical surfaces.

A similar summary is given by

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 ${}_{a}C$, ${}_{b}C$, and ${}_{c}C$ fo For data from Figure 2, Figure 8 and Figure 9. Taken from [26].

	RMSE	aC	bC	cC		AD	aC	bC	cC
PCC					DIE
aC		-	0.18500	0.25909	aC		-	0.10066	0.18512
bC		0.99850	-	0.16511	bC		0.00464	-	0.12144
cC		0.99627	0.99810	-	cC		0.00714	0.00250	-

for the present reconstruction methods: parabolic interpolation of the echelon−approximations ${}_{a}C$, ${}_{b}C$, and ${}_{c}C$ and in Tab 4 for parabolic, hyperbolic, and Gaussian interpolation of the same echelon−approximation ${}_{c}C$.

We can see the optical cuts detected in the data of Figure 2 using focusing criteria (29) (30) in Figure 11, 2D reconstruction (sharp 2D Image) of the same data using the same criterion in Figure 12. Echelon approximation is a simple method for constructing a rough 3D model of the object with all points belonging to the same optical cut having an identical height—the height of the corresponding zone of sharpness—see Figure 13. We can also generalize the notion of low−pass filters used in image processing. These filters are good for smoothing the echelon approximation. This approximation looks much better than the echelon one—see Figure 14.

6.4. 2D and 3D Reconstructions

We can see optical cuts detected in the data of Figure 2 using focusing criteria (28) in Figure 11, 2D reconstruction (sharp 2D Image) of the same data using the same criterion as in Figure 12. Echelon approximation is a simple method for constructing a rough 3D model of the object, where all points belonging to the same optical cut have the same height—the height of the corresponding zone of sharpness—see Figure 13.

We can also generalize the notion of low−pass filters used in image processing. They can be used for smoothing an echelon approximation. This approximation looks much better than the echelon one—see Figure 14.

In Figure 15 and Figure 16, we can see reconstructions of the limestone and blue marble sample by data registration according to Section 3, with focusing criterion (28) and profile height calculation (29) and (30). We can compare a single pore of hydrated Portland cement paste reconstructed by Olympus factory software (Figure 17) and the same pore constructed using focusing criterion (28), probability distribution function (29) and (30) for r = 2 (Figure 18).

Two examples of SFF method results can be downloaded in Supplementary Material.

7. Conclusions

The SFF method based on the Fourier transform can provide correct 3D replicas of rough surfaces. In the case of small samples, a qualified user of this method can obtain results similar to or even better than reconstructions from a confocal microscope. For larger objects, 3D scanners and similar significantly more expensive devices can be simulated by sophisticated mathematical instruments and advanced programming techniques.