1. Introduction

Trace detection is common in criminal investigation^[1], archaeology^[2], and other fields. Over the last two decades, reflectance transformation imaging (RTI) has played a prominent role in trace detection since it was proposed^[1,3]. RTI can provide a powerful enhancement algorithm that interactively reconstructs two-dimensional (2D) images of objects in any virtual light direction; however, it cannot display three-dimensional (3D) shapes and cannot provide the trace size. Later, the combination of 3D reconstruction and 2D detection has been attempted. Artal-Isbran et al.^[4] used two surface examination methods, RTI and 3D laser scanning confocal microscopy, to evaluate the relief line and the contour line on Greek red-figure vases, which was time-consuming. Elfarargy et al.^[5] used the diffuse map and normal map to generate a height map, but the detailed texture of the surface was lost.

In past decades, 3D imaging has been the interest of researchers. Frankot and Chellappa^[6] proposed a global integral method, but with high accumulative errors. The resolutions of the structured light method^[7,8] and integral imaging^{[9–12" target="_self" style="display: inline;">–12]} were subject to wavelength. Li et al.^[13] developed a dynamic virtual camera method, where enough feature points were desired. Besides, combining meta-surfaces and lenses could realize high-performance 3D imaging, however, it was still in the experimental stage^[14,15].

However, most of these common 3D methods focus on the reconstruction of macroscopic objects. There are few reports on the 3D reconstruction of faint traces of millimeters and microns. Therefore, in this paper, we use RTI to enhance the weak traces on the front and back of metal and obtain the surface normal data at the same time. Based on these normal data, the 3D morphologies on both sides of the metal are reconstructed, and the depth information is captured. This method can obtain not only the 2D image of weak traces clearly but also the 3D morphology and 3D information, which provides a powerful tool for archaeology and criminal investigation.

2. Method

The reflection characteristics of the object directly affect its realistic rendering, so the brightness L of each pixel (u,v) of the object image is related to the direction of the incident light. RTI models this dependence by a biquadratic polynomial: $$\begin{gathered} L=a0(u,v)lu2+a1(u,v)lv2+a2(u,v)lulv \\ +a3(u,v)lu+a4(u,v)lv+a5(u,v), \end{gathered}$$(1)where lu,lv are the projections of the normalized direction of incident light on pixel coordinates, L is the illumination of each pixel, and a0(u,v)–a5(u,v) are the six coefficients of the polynomial calculated and stored for each pixel (u,v). The difference in the polynomial coefficient stored in each pixel means that the albedo of the object surface is also different. Due to the curvilinear nature of the biquadratic, it is expected that the high frequencies of the trace will be slightly smoothed, but the impact is acceptable.

Figure 1 shows the image acquisition device, which is composed of a fixed digital single lens reflex (DSLR) camera (Canon EOS 5D Mark IV) with a resolution of 6720×4480 pixels directly above the detection plane and a motor-driven RTI half-arc with a radius of 0.5 m. Nine LEDs (Cree XLamp XM-L), lit in turn, facing the right center of the detection plane, are fixed on the RTI arc at equal angle intervals (15°). The RTI arc rotates six times by 30° in turn from the angle 15° with respect to the detection plane. Simultaneously, the camera (1/25 s, F 8.0, ISO 250) captures 54 images with a small black ball to record the direction of the incident lights (red box in Fig. 1), which can then be calculated from the relative positions of the highlights. Then, we can solve the coefficients of the biquadratic polynomial for each pixel by singular value decomposition (SVD): [lu02lv02lu0lv0lu0lv01lu12lv12lu1lv1lu1lv11⋮⋮⋮⋮⋮⋮lu532lv532lu53lv53lu53lv531][a0a1⋮a5]=[L0L1⋮L53].(2)

Fig. 1. Image acquisition device and the distribution of the direction of the incident light.

下载图片 查看所有图片

With the polynomial coefficients of each pixel, we can reconstruct the image under any virtual illumination direction, providing convenience in capturing the panorama of traces. After deriving Eq. (1), we can get the normalized surface normal N→: {lu0=a2a4−2a1a34a0a1−a22lv0=a2a3−2a0a44a0a1−a22N⇀=(lu0,lv0,1−lu02−lv02).(3)

These noramls are stored in a normal map, where the RGB values of each texel represent the normal components X, Y and Z. X and Y are in the range [−1,1], and Z is in the range [0,1]. To facilitate storage, these data are mapped to [0,255], forming a normal map.

Depth recovery by gradient integration first assumes that the surface depth of the object is Z(x,y), and then the gradients p, q and the partial derivative of the gradients px, qy are defined as follows^[16]: {Zx=∂Z∂x=p,Zxx=∂Zx∂x=pxZy=∂Z∂y=q,Zyy=∂Zy∂y=qy.(4)

Then, we can obtain the error equation: $$\begin{gathered} E(Z;p,q)=(Zx−p)2+(Zy−q)2 \\ +λ[(Zxx−px)2+(Zyy−qy)2]. \end{gathered}$$(5)

Therefore, the problem of solving the optimal depth information Z(x,y) is shifted into finding the minimum solution of the continuous function: const=∬E(Z;p,q)dxdy.(6)

Conventional algorithms focus on restoring the 3D shape of a macroscopic object. However, the traces detected are only a few hundreds or dozens of microns in size, and the noise introduced in the measurement has a great influence on the recovery of 3D morphology. Therefore, we proposed a data processing method, which implements the surface normal data into the surface gradient data, magnifies the micro change of the gradient, and reduces the height of the accumulated error.

Firstly, we convert the normal n→=(nx,ny,nz) of each pixel as follows: {p=−nxnz+w×(−nxnz+nx)q=−nynz+w×(−nynz+ny).(7)

In this step, the normals are converted into gradient data of the surface height to facilitate depth recovery; moreover, slight surface changes can be amplified by the weight coefficient w.

Secondly, to avoid some gross errors, we select threshold T to screen p and q. If |p|>T or |q|>T, we set them to zero.

Then, we choose a local window of size d with each pixel of the gradient map as the center to calculate the mean value of the window as a reference value and amplify the gap between the traces and the background to enhance the visual effect of 3D reconstruction: {p′=p+k[p−1d2(∑dpi)]q′=q+k[q−1d2(∑dqi)].(8)

After that, we obtained the new gradient maps P and Q of the surface with the size of N×M, which we extended to avoid the periodicity problem of the Fourier transform: {P˜=[P(x,y)−P(N−x,y)P(x,M−y)−P(N−x,M−y)]Q˜=[Q(x,y)Q(N−x,y)−Q(x,M−y)−Q(N−x,M−y)].(9)

Finally, the depth information of trace is recovered by a global integral^[6]. We used a high-precision DetakXT profilometer (Bruker, Billerica, MA, USA) (128 µm of scanning range, 1% accuracy) to scan the dimensions of a series of traces, and approximated these data to the truth value of these dimensions. According to the ratio of the data from the profilometer to the data from RTI, we determined the fixed relationship between physical size and pixel size. The algorithm is implemented under MATLAB R2019b on a personal computer (PC) with 2.6 GHz CPU (Intel Core i5-4210M), 8G RAM and 2G Graphics card (NIVIDIA GTX960M).

3. Results and Discussions

A metal sheet with Chinese characters printed on its surface was selected as the experimental sample. Because of oxidation and corrosion, it is difficult to get the information of these characters directly and clearly; there is almost no trace on its back [Figs. 2(a) and 2(c)]. The normal maps of these two surfaces are shown in Fig. 3. Since the metal sheet is a plane, and its normal almost points to the direction of the Z axis, the normal maps are nearly all blue, and part of the color change represents the change of the normal, which indicates the trace information. These maps also show that the changes in normal of traces are weak. RTI with an enhanced algorithm helps capture the trace information easily and accurately [Figs. 2(b) and 2(d)]. The selected red box in Fig. 2(b) shows the most complicated structure of the Chinese character clearly. Then, after this part of the normal data is extracted, the depth recovery algorithm is tested.

Fig. 2. Sample and RTI reconstructions. (a) Original photo of the front. (b) RTI reconstruction of the front with specular enhancement. (c) Original photo of the back. (d) RTI reconstruction of the back with specular enhancement.

下载图片 查看所有图片

Fig. 3. Normal maps of the two sides of the metal in Fig. 2. (a) The normal map of the front. (b) The normal map of the back.

下载图片 查看所有图片

Figure 4 is the result of the 3D reconstruction of the trace in the red box in Fig. 2. Figure 4(a) is the result from the common photometric stereo (PS) algorithm without filtering normal data, where there are some isolated peaks, indicating that there are some gross errors in the gradient data that need to be eliminated. Figure 4(b) shows the recovery using the PS algorithm (λ=0) after setting the threshold T=1 to screen the normal data. The burr noise is eliminated, but the whole surface is not flat enough. Then, as shown in Fig. 4(c), we set w=0.5 to convert the normal data and set λ=0 again for 3D recovery, which improves the overall face shape and edges. Finally, as shown in Fig. 4(d), we get a better shape, clearer edge, and lower height error with our method.

Fig. 4. 3D recovery results of the trace in the red box in Fig. 2(b). (a) Results without filtering normal data. (b) Recovery by PS algorithm after data filtering (T = 1, w = 0, k = 0, λ = 0). (c) Recovery using PS algorithm after data filtering and data conversion (T = 1, w = 0.5, k = 0, λ = 0). (d) The result of our algorithm (T = 1, w = 0.5, k = 0, λ = 1). (e) The planes taken from the same parts marked by different colors of the above reconstructions and MSEs of these parts. (f) The contours of the above reconstructions along the direction indicated by the colored arrows.

下载图片 查看所有图片

Meanwhile, as shown in Fig. 4(e), we intercepted the same flat parts of the recovery results, marked them with different colors, and calculated their mean square error (MSE) to indicate flatness. It is not difficult to find that the flatness of the surface recovered by our method is better, and the accumulated height error is also lower. Moreover, we also captured the relative depth of the same edges of the traces recovered by the above methods along the directions indicated by the colored arrows, as shown in Fig. 4(f). The result confirms the noise in Fig. 4(a) and illustrates that our method can effectively reconstruct the contour of the trace and reduce the accumulated height error. In short, proper w and λ can improve the surface and edge of reconstruction.

In experiments, we found that w=0.5 and λ=1 were good enough to rebuild the 3D shape, so we kept both settings for the rest of the experiment. Meanwhile, by adjusting the magnification factor k and window size d, we can get different 3D effects and contrast. For Fig. 2(d), as different k and d were applied, the effects are shown in Fig. 5. It is not difficult to find that proper enlargement of k and d can effectively help us observe the 3D structure of faint traces. Figure 5(a) shows the 3D shape without enhancement, where the trace and background are partly blended. As we increased k and d, the trace and background were efficiently separated, and the edge of the trace was significantly magnified [Figs. 5(b)–5(d)]. The 3D reconstruction of the metal back is shown in Fig. 6(a); in contrast with Fig. 2(c), some traces on the surface are faintly revealed, whereas, since the metal is much larger than traces, these traces are somewhat submerged. Instead, we only focused on the traces on the surface and enhanced the 3D reconstruction after capturing the regions of the traces from the 2D enhancement results. As shown in Fig. 6(b), the 3D information of these traces is effectively extracted.

Fig. 5. Different enhanced 3D topographies of the trace marked by red arrows in Fig. 2(d). (a) The 3D reconstruction without enhancement (k = 0). (b) The 3D topography enhanced by k = 15 and d = 15. (c) The 3D topography enhanced by k = 25 and d = 15. (d) The 3D topography enhanced by k = 25 and d = 25.

下载图片 查看所有图片

Fig. 6. (a) Metal back 3D reconstruction (w = 0.5, k = 25, d = 25, λ = 1). (b) Trace 3D reconstruction on the back (w = 0.5, k = 25, d = 25, λ = 1).

下载图片 查看所有图片

Moreover, along the direction of the red arrows in Figs. 2(b) and 2(d), we scanned both parts by a DetakXT profilometer, and the results are shown as black curves in Fig. 7. Given the fixed relationship between pixel size and physical size, we set k=0 for 3D reconstruction of these parts, obtaining the altitude lines of their longitudinal sections. Due to the high precision of the DetakXT profilometer, we can evaluate the depth recovery by calculating the degree of coincidence between the depth recovery curve and the scanned curve. Therefore, we refer to the R2 value: R2=1−∑i(zi−z^i)2∑i(zi−zI¯)2,(10)where zi is the point on the red curve of the depth recovery, z^i is the corresponding point on the black curve, and zi¯ is the mean of the black curve. It is generally believed that the closer R2 value is to one, the more similar the two curves are. We calculated that the R2 value of Fig. 7(a) is equal to 0.8772, and the R2 value in Fig. 7(b) is equal to 0.9387.

Fig. 7. Scanning results with the DetakXT profilometer (128 µm scanning range, 1% accuracy) and depth recovery results of the same part. The scanning position and direction are shown as red arrows in Fig. 2. (a) Scanning results and depth recovery of Fig. 2(b), R² = 0.8772. (b) Scanning results and depth recovery of Fig. 2(b), R² = 0.9387.

下载图片 查看所有图片

As shown in Fig. 7, the recovery results of the parts with gentle changes were basically identical with the scanning results, while there were differences in the parts with abrupt changes, which is exactly what we predicted. Because of the properties of biquadratic polynomials, the high-frequency parts will be partially smoothed. Some other reasons, such as the approximate parallelism of beams and the limited number of sampling points, may lead to errors. Fortunately, only in the sharp parts is the error slightly larger, but also within 10 µm, and the depth information of interest is still available.

4. Conclusion

First of all, we describe the basic principle of RTI, and, then, for the faint traces below the millimeter level, we propose an optimized 3D reconstruction algorithm, which can effectively suppress the noise, improve the edge and surface shape, and enhance the 3D morphology. Next, the oxidized metal sheet was tested using the RTI device, and the traces on the metal surface were clearly visible in the 2D image with the specular enhancement. Finally, 3D reconstruction and depth information calculation are carried out on the weak traces. The results show that our method can realize the 3D reconstruction of tens of micron depth traces, and the error of depth recovery is less than 10 µm. Promisingly, our method has an important application prospect in the fields of cultural relic protection and criminal investigation.