Data fusion of surface normals and point coordinates for deflectometric measurements

Measuring specular surfaces can be realized by means of deflectometric measurement systems with at least two reference planes as proposed by Petz and Tutsch (2004). The results are the point coordinates and the normal direction of each valid measurement point. The typical evaluation strategy for continuous surfaces involves an integration or regularization of the measured normals. This method yields smooth results of the surface with deviations in the nanometer range but it is sensitive to systematic deviations. The measured point coordinates are robust against systematic deviations but the noise level is in the order of micrometers. As an alternative evaluation strategy a data-fusion process that combines both the normal direction and the point coordinates has been developed. A linear fitting technique is proposed to increase the accuracy of the point coordinate measurements by forming an objective functional as the mean squared misfit of the gradients with respect to the point coordinates on the one hand and to the normals on the other hand. Moreover, a constraint on the maximal change of the coordinate measurements is added to the optimization problem. To minimize to objective under the constraint a projected gradient method is used. The results show that the proposed method is able to adjust the point coordinate measurement to the measured normals and hence decrease the spatial noise level by more than an order of magnitude.


Introduction and measurement principle
The combination of close range photogrammetry and structured illumination is well established in the field of threedimensional measurements of diffusely reflecting surfaces.The so-called fringe projection technique is based on the triangulation principle and makes use of at least one electronic camera and one projector.By means of structured illumination, which especially involves phase shifting techniques, the surface under test is optically coded.Based on this spatial coding, the camera images can be evaluated in a way that delivers a three-dimensional object point for each image pixel.
This fringe projection technique depends on the optical imaging of the surface under test onto the image sensor of the camera.This requires a minimum amount of diffusely scattered light from the surface.For reflecting surfaces the amount of diffusely scattered light is very low -in the ideal case there should be no scattered light at all.Reflecting surfaces therefore are not directly visible, but their presence can affect the propagation of light in a characteristic way (Tutsch et al., 2011).
It is a common experience that an image generated by reflection at a curved mirror shows characteristic deformations.These optical deformations of an object's mirror image carry information about the geometry of the reflective surface.If the imaged object is a known regular pattern, observed deformations of the pattern can be taken as an indication of form deviations or defects of the surface under test.This technique, known as deflectometry, is well suited for the qualitative inspection of reflective surfaces.An absolute shape measurement, however, is not possible directly.As can be seen in Fig. 1, in general there are various combinations of position and orientation of a mirror element that lead to the same observation.
Published by Copernicus Publications on behalf of the AMA Association for Sensor Technology.(Tutsch et al., 2011).
Figure 2. Principle of an enhanced deflectometric approach using two reference pattern positions (Tutsch et al., 2011).
This ambiguity can only be resolved by additional information.For instance assumptions can be made concerning the properties of the mirror surface.Then the surface can be reconstructed by means of regularization (Werling et al., 2009).However, the parameters of this reconstruction process have to be determined by other means, such as approximation or additional measurements for at least one point on the surface (Li et al., 2012).Petz and Ritter (2011), Petz andTutsch (2004), andPetz (2006) therefore proposed an enhanced deflectometric approach, shown in Fig. 2, that solves the ambiguity discussed for the basic approach.The enhanced approach uses at least two different positions of the reference pattern, which can easily be realized by moving the reference structure with a linear stage.With this approach not only the camera defines a unique ray for each image point, but also the information  (Petz, 2006).
from the spatially coded reference patterns delivers an unambiguous ray.The intersection of these two rays defines the location of the object point.In addition to the threedimensional position of the object point also the surface normal in that point can be determined, as for the reflective surface the law of reflection must be valid.
A technical realization of the approach illustrated in Fig. 2 can be seen in Fig. 3 (Petz, 2006).As reference structure a common LCD-TFT display is used.The absolute twodimensional spatial coding of the display area is accomplished by a combination of phase shifting techniques and a heterodyne approach for deconvolution.The display is mounted onto a motorized linear stage and can thus be moved to various positions during the measurement process.The object under test will be mounted in the position where the figure shows a calibration mirror, so that the camera observes a mirror image of the LCD display.
As mentioned above, the enhanced deflectometric approach directly delivers three-dimensional coordinates of the observed surface points.Thus there is no need for surface reconstruction by means of regularization or integration of the surface slope.However, a closer investigation of the measurement data reveals that the direct measurement of coordinates is more sensitive to noise than the measurement of the surface normals.In earlier evaluations Petz (2006) has shown that the noise of the measured surface is approximately 3 orders of magnitude higher for the direct coordinate measurement than it is for the reconstruction based on the integration of the surface slopes.
However, the surface reconstruction based on integration is susceptible to local disturbances and systematic errors.Therefore an advantageous approach for the surface measurement in deflectometry should consider both, the absolute coordinates determined by triangulation and the observed surface normals.A correspondent approach is the subject of the present paper.

Notation and data type
In this paper the Euclidean norm • on R 3 will be used.Vectors v will be written in italic, bold letters and matrices M in bold, capital letters.Scalars c will be written in italic letters.A matrix A, when of dimension m × n × 3, is of form A = (A x , A y , A z ) with A x , A y , A z matrices of dimension m × n.The vector in the ith row and the j th column of A is denoted with the lowercase, italic, bold letter v i,j = (v x i,j , v y i,j , v z i,j ) T ∈ R 3 with its x, y and z coordinates.If any vector w or matrix A is corresponding to m i,j it will be written as w i,j or A i,j .
The type of data set being worked with is the result of a deflectometric measurement process as described earlier in Sect. 1.After postprocessing the data set contains a matrix P = (P x , P y , P z ) of measured point coordinates together with its information of neighborhood and a matrix N = (N x , N y , N z ) of measured normal vectors.The normal vectors n i,j are corresponding to the point coordinates p i,j .
If some measurement point did not yield a valid measurement, the corresponding entries in P and N will be undefined.Hence, each point coordinate has at most eight neighbors, i.e., the eight neighborhood.Point coordinates can have fewer neighbors because of invalid measurements and also at the boundary of the measured object.

Objective functional for data fusion
As mentioned in Sect. 1 the accuracy of the normal vectors is approximately 3 orders of magnitude higher than the accuracy of the measured point coordinates.The purpose of this section is to combine both measured data sets with the objective to get consistent data.
To that end an objective functional will be defined that will fuse the data in a process of optimization.According to the measurement process, it is known that the x and the y coordinates of the point coordinates are more trustworthy compared to the z coordinates.Hence, the objective is to adjust the z coordinate with respect to the orientation of the measured surface.
The first idea is to take a look at the normal vectors.For every measured point coordinate a calculated normal can be defined by the cross product with two neighboring point coordinates.Then the normals can be compared by minimizing the mean squared misfit of the measured and the calculated normals.However, this objective functional leads to an optimization problem which is nonlinear with respect to the z coordinate.This is because the mapping from the point coordinates to the normals is highly nonlinear.This kind of optimization problem is surely solvable, but costly in realization.
By rethinking about information of surface orientation beside normal vectors, gradient vectors come to mind.Gradients carry the same information about the surface structure as normals.Hence, gradients can be calculated separately with respect to each measured normal and to the z coordinates of each point coordinate.
Thus, an objective functional can be formed as the mean squared misfit of two kinds of calculated gradients.On the one hand, gradients calculated with respect to the normal vectors and, on the other hand, to the point coordinates.With this approach a quadratic objective functional will be obtained.
First the gradient vector g corresponding to the normal vectors is formed by rearranging the vectors.The first mn rows of g are the partial derivatives of the z coordinate with respect to the x coordinate and the second mn rows are the partial derivatives with respect to the y coordinate.Thus, g is given by and is of dimension 2mn.
Next the gradient with respect to the point coordinates will be calculated.In order to do that, all neighbors in the eight neighborhood of one point coordinate p i,j , for all i = 1, 2, . .., m and j = 1, 2, . .., n, will be located, and build the vectors v 1 , . .., v 8 from p i,j to its neighbors.The vectors v 1 , . .., v 8 have form The point coordinate p i,j with its possible neighbors and the vectors v k , k = 1, . .., 8, are shown in Fig. 4. Note that the grid defined by the point coordinates is usually not regular.
On the one hand the directional derivatives of p i,j in direction to its neighbors can be calculated with v 1 , . .., v 8 as follows: All directional vectors v T k can be written line by line in a matrix V i,j , which is then at most of dimension 8 × 2. All neighbors of a point coordinate are considered into the calculation to obtain a more robust solution.In order to be able to calculate the matrix V i,j at all, at least two neighbors are www.j-sens-sens-syst.net/3/281/2014/J. Sens. Sens. Syst., 3, 281-290, 2014 Figure 4. Neighborhood of p i,j with vectors v 1 , . .., v 8 pointing to all direct neighbors of p i,j .
needed.It can happen that a point coordinate has fewer than eight neighbors, e.g., at boundary points, so that V i,j has a smaller number of rows.
On the other hand, the directional derivatives in direction of the neighboring points can be calculated by taking the differences of the z coordinates of the neighbors of p i,j and p i,j itself: Putting both Eqs. ( 1) and (2) for the directional derivatives together, the result is a linear system for ∇p z i,j .
The least-squares solution of this overdetermined system for the partial derivatives ∇p z i,j is given by with (V i,j ) + the pseudo inverse of V i,j .The matrix D i,j = (V i,j ) + • C is the calculated difference matrix of the point coordinate p i,j .
Since the objective is not to calculate the partial derivatives separately but to calculate all in one, the entries of each difference matrix D i,j can be written in the proper row and column of a matrix D, which contains all differences for the partial derivatives with respect to x in its first mn rows, and all differences for the partial derivatives with respect to y in the second mn rows.So D is of dimension 2mn × mn.Note that since one point coordinate has at most eight direct neighbors, D will be a sparse matrix with at most nine nonzero entries per row.
In that way, an objective functional is given as the misfit of the gradients corresponding to the point coordinates and to the normals, which will be minimized with respect to P z .The objective functional F is quadratic in the vectorized z coordinates of the points.Note that the gradient vector g and the difference matrix D have only to be calculated once at the beginning of the process.

Linear fitting with steepest descent method
The objective functional was formed in Sect. 3 and it is quadratic in P z .In this section the focus goes to the measured data and the optimization process itself.
An optimal P z corresponding to the x and y coordinates can easily be found by shifting the z coordinates.Since the point coordinates are given from a deflectometric measurement process, the z coordinates cannot be shifted without considering the original position of the measured point coordinates and the occurring measurement uncertainties.Hence, the points can only be shifted in a given maximal range.
So the estimated accuracy of the coordinate measurement is used to add a constraint to the optimization problem on the maximal change of the measured points.We assume that the point coordinates are disturbed by noise that follows a Gaussian distribution in the z direction with known variance σ .Hence, the distribution of the sum of all squared distances in the z direction follows from a chi-squared distribution with the number of degrees of freedom equal to the number N of valid measurement points.We deduce that the expected value of this sum is equal to N σ 2 .We use this value as tolerance δ = N σ 2 and define the set in which the point coordinates are allowed to be shifted as B = { P z − P z 2 ≤ δ}, the sphere with radius δ around the measured z coordinate P z .The shifted z coordinate is given by P z .Therefore, with objective functional F (P z ) = 1 2 g − D • vec(P z ) 2 , the optimization problem becomes min P z F (P z ) s.t In order to minimize the objective functional under the constraint, a projected gradient method will be used.Therefore, the search direction for the optimum is given by the negative gradient of F , which is and, hence, the update for each optimization step k with respect to the optimization constraint is with τ k the step size for the iteration k.
The projection of the point coordinates P z k onto the set B is given by Now, a sensible choice of the step size τ k is to be made.By observing the objective functional F , it can be seen that a convex quadratic problem is to be solved.
The valid step size is The step size with respect to the gradient of F at P z k has to be calculated for each step k.
Putting all results together, the data-fusion algorithm is as follows.Now, a sensible choice of the step size τ k is to be made.By observing the objective functional F, it can be seen that a convex quadratic problem is to be solved.
The valid step size is The step size with respect to the gradient of F at P z k has to be calculated for each step k.
Putting all results together, the data-fusion algorithm is as follows: Algorithm Data-Fusion

Experimental Results and Conclusions
To verify the feasibility of the algorithm, a MATLAB implementation of the algorithm has been created.Two sets of measurement results were used.On the one hand a plane mirror and on the other hand a spherical mirror.In the first case all normal directions of the point coordinates are almost the same, since the measured object is plane.But, a calculation of normals with respect to the measured point coordinates leads to deviation up to 20 • at some points.With the sec-

Experimental results and conclusions
To verify the feasibility of the algorithm, a MATLAB implementation of the algorithm has been created.Two sets of measurement results were used.On the one hand a plane mirror and on the other hand a spherical mirror.In first case all normal directions of the point coordinates are almost the same, since the measured object is plane.But, a calculation of normals with respect to the measured point coordinates leads to deviations of up to 20 • at some points.With the second data set similar results in deviations are observed.Since problems can occur during data fusion, further examinations of the data type itself will be done.

Implementation
In practice, a point coordinate at the boundary of the measured object will have fewer neighbors than eight.Also, point coordinates within the object can have fewer neighbors, e.g., if the measured object has a hole within, has scratches or if some coordinate points could not be measured.
Thus, when the matrix V i,j is formed, there can be entries that are not defined and so V i,j (cf.Sect.3) will have a smaller number of rows.Thus, the matrices C i,j and D i,j and the vector z i,j will also be smaller during the calculation of the gradient of p z i,j .In this case, all rows and columns corresponding to the undefined point coordinates have to be deleted.
After calculating D i,j , the deleted rows and columns of the matrix will be reinserted by writing zeros in the proper row and column, thereby acquiring its original size of 2 × 8.For insertion, zeros were chosen, because these entries should not and will not be changed during the calculation process.After that the entries of D i,j can be written within the matrix D.
The matrix D is a sparse matrix, since one point coordinate has at most eight direct neighbors.This matrix is of dimension 2m × n, the first m rows are corresponding to the partial derivatives with respect to x, the second to the partial derivatives with respect to y.The principal diagonal of the first m rows and of the second m rows of D has entries corresponding to the point coordinate p z i,j and so on.So the entries of each matrix D i,j can be saved in vectors, which can be written on the diagonals of D.

Results
Different data sets were used for the tests.Here, results of the data-fusion process will be shown in direct comparison of the data before and after the calculation.The comparison will be shown by two different measured objects.These objects are a plane reference mirror and a spherical reference mirror.
In Fig. 5 the logarithm of the value of the objective functional for both measured objects is shown.In both semilogarithmic plots it can be seen that the value is falling even by 1000 iteration steps.
As an a posteriori validation of the assumption that the noise in the z coordinate of the coordinate measurement was indeed Gaussian distributed, we show histograms of the residuum of the originally measured z coordinate with the new z coordinate after data fusion in Fig. 6.It turns out that the histograms are approximately Gaussian and, since our data-fusion model in fact produces consistent data of point coordinates and normals, this verifies a posteriori that the assumption of a Gaussian distribution was valid.
In Fig. 7 the gradients of the measured plane and measured spherical reference mirror were plotted.The two plots on the left-hand side show the gradients of the plane mirror before and after the data-fusion process.The two plots on the right-hand side show the same for the spherical mirror.In every plot of this figure the gradients with respect to the measured normals are given in red and the calculated gradients with respect to the point coordinates are given in green.As one can see, the measured gradients all are pointing almost in the same direction, which is reasonable because the two measured objects did not have any visible jumps on the surface.The calculated gradients before data fusion otherwise differ considerably from the calculated ones.After data fusion the gradients with respect to the fitted point coordinates were calculated.On each of the right-hand side plots it can be seen that the new calculated gradients lie on the measured gradients.Hence, the data-fusion process leads to same information of the orientation of the measured surface.
In Fig. 8, the point coordinates before and after data fusion are shown.The surface is colored according to the deviation of the angle between the measured normal and the normal calculated with respect to the point coordinates before and after the optimization process.The color bar is given in degrees.
Before data fusion, the error between the normal vectors on the surface defined by the point coordinates and the measured normal vectors was up to 20 • at some coordinates.After optimization the maximal error was below 1 • in all coordinates.
In Fig. 9 the flatness deviation of a cross section through the plane mirror is shown.The original data obtained from triangulation shows a noise amplitude of about 20 µm.The data-fusion process reduces the noise level by approximately an order of magnitude.For the spherical mirror very similar results are obtained.Gradients after data fusion Gradients of the plane mirror . Left: measured point coordinates.Right: point coordinates after optimization.The red vector is a measured normal, the green a calculated normal.The color encodes the error between calculated normals and measured normals in degrees.

Conclusions
In this paper a novel evaluation approach for a deflectometric measurement has been developed.The results of such a measurement, when carried out with multiple reference planes and active triangulation of the incident and reflected light rays, are the three-dimensional point coordinates and normal directions for each valid pixel of the camera.Due to the local and absolute measurement, the point coordinates are very robust against systematic deviations but with stochastic deviations in the order of 10 µm.However, the commonly applied integration or regularization of the normal directions yields a smooth surface with stochastic deviations in the nanometer range but is sensitive to global form deviations due to systematic influences.The data-fusion method proposed in this work combines both data sets by adjusting the point coordinates to the measured normal directions such that the data is consistent up to the measured tolerances.This adjustment is done in a uniform way such that no systematic deviations are expected to occur.The method has been tested with real data sets for a planar and a spherical mirror.It has been shown that the spatial noise level of the resulting point cloud could be reduced by more than an order of magnitude without sacrificing any of the major advantages of a local and absolute measurement, such as the ability to measure nonsmooth surfaces.

Figure 1 .
Figure1.Ambiguity of determination of position and orientation of a mirror element in deflectometry(Tutsch et al., 2011).

Figure 5 .
Figure 5. Logarithm of functional value during iteration for a measured plane reference mirror (lower curve) and for the spherical reference mirror (upper curve).

Figure 6 .
Figure 6.Top panel: residuum of the measured z coordinate and the z coordinate after data fusion in the plane reference mirror case.Bottom panel: residuum of the measured z coordinate and the z coordinate after data fusion in the spherical reference mirror case.

Figure 7 .
Figure 7. Clipped area of plane and spherical mirror gradients before and after the data-fusion process with 1000 iterations.Measured gradients in red, calculated gradients from the point coordinates in green.

Figure 9 .
Figure 9. Flatness deviation of one row of the plane mirror before and after data fusion.