Among the main advantages of interferometric techniques, the possibility to carry out full-field measurements is certainly one of the most appealing. When interferometric techniques, like photoelasticity [1
1. M. M. Frocht, Photoelasticity, Vol. I (John Wiley and Sons, 1941).
], moiré interferometry [2
2. D. Post, B. Han, and P. Ifju, High sensitivity moiré (Springer Verlag1997).
], speckle interferometry [3
3. R. K. Erf, Speckle metrology (Academic Press1978).
] or digital holography [4
4. U. Schnars and W. Jueptner, Digital holography (Springer2005).
] are applied, the experimental data typically consist of fringe patterns which represent the contour map of a particular quantity. If a real full-field approach is desired it is necessary the application of a phase-shifting procedure [5
5. K. Creath, “Temporal phase measurement methods,” in Interferogram analysis, D.W. Robinson and G.T. Reid, ed. (Institute of Physics Publishing, 1993).
6. M. Kujawinska, “Spatial phase measurement methods,” in Interferogram analysis, D.W. Robinson and G.T. Reid, ed. (Institute of Physics Publishing, 1993).
], in order to obtain the phase maps by means of a number of fringe patterns. Actually the phase maps provide point-wise measurements, although an unwrapping procedure becomes necessary in order to remove the phase jumps always present when the optical path difference exceeds one wavelength.
Aside from the drawbacks (many of which are nowadays overtaken) that arise by unwrapping a phase map particularly noisy [7
7. D. C. Ghiglia and M. D. Pritt, Two-dimensional phase unwrapping (John Wiley & Sons, 1998).
], a residual noise on the measurements is still present. The entity and the origin of this noise is strictly dependent on the type of technique and on the operative experimental conditions.
Many are the attempts available in the scientific literature to smooth the noise either in the space domain or in the spatial frequency domain [8–12
8. H. A. Aebischer and S. Waldner, “A simple and effective method for filtering speckle-interferometric phase fringe patterns,” Opt. Commun. 162, 205–210, (1999). [CrossRef]
], but, as it is shown in the papers dealing with these approaches, it is not possible to obtain a smooth distribution of the phase without fluctuations.
The only way to obtain a smooth phase distribution consists in using an approach based on a fitting operation [13
13. J. Novak and A. Miks, “Least-squares fitting of wavefront using rational function,” Opt. Lasers Eng. 44, 40–51, (2005).
]. By this approach, very popular in the field of reverse engineering [14–15
14. C. L. Bajaj and G. Xu, “Spline approximations of real algebraic surfaces,” J. Symbolic Comput. 23, 315–333, (1997). [CrossRef]
], it is necessary to define an analytical expression based on a number of parameters which must be evaluated by a proper optimization procedure. If a reliable analytical model is available for the problem under investigation the fitting can be carried out with high probability of success, otherwise the choice of the analytical model becomes a very crucial step of the whole fitting procedure.
A typical analytical model is based on the use of polynomial functions which are particularly easy to manage and allow to obtain always a system of linear equations. On the other hand the need of higher modeling capability brings to increase the order of the polynomial and it implies unwanted fluctuations of the model, which does not follow the physical nature of the phenomena represented by the phase maps. The use of particular polynomial functions, like orthogonal Legendre or Zernike [16
16. F. Furgiuele, M. Muzzupappa, and L. Pagnotta, “A full-field procedure for evaluating the elastic properties of advanced ceramics,” Exp. Mech. 37, 285–291, (1997). [CrossRef]
] polynomials, can improve the performance of the polynomial fitting, but the instabilities can be only lowered or shifted outside the domain of interest, not removed. Further improvements can be obtained by adopting different functions or nonlinear models [17
17. D. M. Bates and D. G. Watts, Nonlinear regression analysis and its applications (John Wiley & Sons, 1988). [CrossRef]
] that lead to use particular solving algorithms. Finally it is important to say that, although most of mono-dimensional fitting is straightforward, when we handle interferometric data the problem becomes bi-dimensional with consequent increase of the calculation complexity.
The present paper reports a bi-dimensional fitting procedure based on the use of B-spline functions [18
18. V. B. Anand, Computer graphics and geometric modeling for engineers (John Wiley & Sons, 1993).
] applicable to differently shaped domains. According to this model a number of piecewise polynomial functions, usually called blending functions, are defined and their existence is limited to a portion of the total domain. The number of the functions and the order of the polynomials can be chosen independently from each other; in this way a complex phase map can be fitted simply by increasing the number of functions without increasing the order of polynomials. The advantages of this approach are the linearity of the model to be optimized, the local control and the possibility to keep low the order of polynomial functions involved in the definition of the model. On the other hand the piecewise formulation implies a more complex mathematical model to be implemented, although it is easily achievable by the actual software environments like MatLab
19. A. Knight, Basics of Matlab and beyond (Chapman & Hall/CRC, 2000).
] or Mathematica
20. S. Wolfram, The Mathematica book, 5th edition (Wolfram Media inc, 2003).
]. Moreover the application of the procedure to a non rectangular domain requires a space transformation consisting of further mathematical manipulations.
The fitting procedure proposed in the present paper was initially calibrated on analytical data by considering as case study the deformation field of a thin square plate supported along the edges and loaded by a uniform pressure. The out-of-plane displacements, calculated by an analytical solution of the problem, were fitted by varying the number and the order of the blending functions. The choice of these operative parameters was performed by the standard deviation evaluated on the difference between the analytical and the fitted data. Subsequently, the capability of the procedure to be noise tolerant was evaluated by introducing an additive and multiplicative noise (with uniform and gaussian distribution) of the same order of magnitude of the typical interferometric techniques.
The results obtained for the analytical case study were applied to experimental data obtained by phase shifting speckle interferometry. The out-of-plane displacements, measured by a Michelson interferometer with a speckle reference beam, were fitted by 49 (7×7) bicubic blending functions. Again the performance of the fitting procedure was tested by evaluating the standard deviation on the difference between the best fit surface and the experimental data. The reference values for standard deviation were calculated by imposing a uniform tilting of the reference beam, that can be considered as known out-of-plane distributions.
2. Mathematical formulation
The mathematical formulation of the fitting procedure is based on the well-assessed B-spline formulation [18
18. V. B. Anand, Computer graphics and geometric modeling for engineers (John Wiley & Sons, 1993).
] according to which the generic point P
) of a surface can be represented by the following equation:
where Ni,h(x) and Nj,k(y) are the mono-dimensional B-spline blending functions, Wi,j are the weights to be evaluated by the optimization procedure, (n+1) and (m+1) are the number of the control points in x- and y- direction, h and k are the orders of the polynomial functions in x-and y- direction. The blending functions can be defined by the following recursive equation:
where ξi is the i-th element of the knot vector whose number of elements is equal to (n+k). By this approach a piecewise polynomial function of (k-1) degree is obtained, that is C
k-2 continuous, and it spreads throughout (k+1) control points.
By operating on the knot vector it is possible to obtain a non-periodic and non-uniform B-spline. The non-periodic formulation is based on the repeated knots at the ends of the knot vector; by repeating these knots k times we obtain a curve starting and finishing from and to the ends knots. In the non-uniform formulation the intervals of the knot vector are not equispaced, and the internal knots can be even repeated.
reports an example of the B-spline blending functions assuming a non-periodic uniform formulation on a rectangular domain; in this example parabolic functions (h
=3) and 4x4 control points (n
=4 and m
=4) were used. It can be seen that each function affects only a limited portion of the whole domain.
2.1 The fitting procedure
The optimization of the model is achieved by evaluating the weights Wi,j. The advantage of this approach consists in reducing the number of parameters to be optimized, and the low degree of the polynomial functions, that avoids unwanted fluctuations.
The data obtained by an experimental technique are arranged on a bi-dimensional array, the dimensions of this array are the resolution of the sensors, typical values are: 640×480, 768×576, 1024×768, 1280×1024 and also other formats are available. It is not important the exact resolution, what matters is that we have to deal with hundreds of thousands of experimental information, out of which any tens of parameters must be evaluated.
The optimization procedure was carried out according to the following steps:
Calibration of the spatial coordinates of the input data;
Calculation of the blending functions;
Construction of the model;
Evaluation of the best fit parameters.
In the first step of the fitting procedure the indices of the input matrix are transformed into the spatial coordinates of the object under investigation, or vice versa. Maybe it is better to work in the object reference system, but sometimes it could be convenient to use the indices of the input matrix. However, if there are not distortion effects, this step is very easy to carry out because it consists of an origin translation and a stretching (or a shrinkage) of the axes.
In the second step the blending functions, like those reported in Fig. 1
, are analytically evaluated in the specified domain. At the moment a rectangular domain is assumed, in the next sections this constrain will be removed.
Fig. 1. Blending functions of a B-spline non-periodic uniform surface with h=k=3 (parabolic approximation), n+1=4 and m+1=4 (16 control points in all).The z-values of the functions are represented as phase maps just for a better visualization.
In the third step the model to be optimized is built as an over-determined system of linear equations by means of eq. (1
). In this way a rectangular matrix M
is obtained with a very high number of rows (as many as the experimental data) and a reduced number of columns (as many as the parameters to be evaluated, the weights Wi,j
). The elements of this matrix are the values assumed by the blending functions at each experimental datum location. The known vector V
is simply the experimental data flattened in a column vector. So assuming the following form for the problem:
where each equation of the system is relative to the l-th datum, and p is the total number of the experimental data. According to this approach the unknown vector W is calculated as the pseudoinverse of M, that can be denoted as M*, multiplied by the known vector V. Another advantage of this formulation is that if more than one set of experimental data is carried out on the same object, the matrix M must be calculated and inverted only the first time, whereas the unknown vectors for each set of data will be evaluated by a simple multiplication between M* and V.
2.2 Data distributed on a circular domain
The blending functions can be easily modified in a cylindrical coordinate system simply by acting on the knot vectors. If for a rectangular domain is rather natural a non-periodic formulation in order to dispose of functions starting and finishing from and to fixed points, for a circular domain a periodic formulation along the circumferential coordinate is necessarily required, while along the radial coordinate a non-periodic approach remains more convenient.
reports the blending functions defined on a circular domain with 6 control points along the radial coordinate and 16 control points along the circumferential coordinate. In the example reported in Fig. 2
cubic functions were used in both directions (h
=4), hence the knot vectors are of 10 and 20 elements in radial and circumferential directions, respectively. A non-periodic formulation was used for the radial coordinate while a periodic formulation is necessary for the circumferential coordinate. In Fig. 2
only the 6 blending functions for a fixed circumferential control point were reported, because in the periodic formulation these functions are simply translated (and eventually rotated in cylindrical coordinates) when moving from a control point to another. In the example reported in Fig. 2
the blending functions relative to the next or to the previous control points are obtained by a rotation of π/8 (=2π/16).
It must be pointed out that the case of the circular domain can be applied also in presence of a hole in the center of the surface, in this circumstance the radial knot vector does not start from zero but from a finite value.
2.3 Data distributed on a domain of any shape
It is not unusual to deal with experimental data distributed on non regular domains, like a rectangle or a circle, in this case the application of the procedure is no more straightforward, but still possible.
In order to achieve the purpose it is necessary to add one step more between steps 1 and 2 of the procedure defined in section 2.1. Another coordinate transformation must be defined with the aim of transforming the initial any shaped real domain into a regular auxiliary domain. This transformation must be invertible, in order to bring back to the real domain the data after the application of the fitting procedure. This operation is not always possible, but in most practical situations it is. In the following the way by which two cases could be faced is presented.
Fig. 2. Cubic blending functions (h=k=4) of a B-spline surface defined on a circular domain. Non-periodic formulation and 6 control points were used in the radial direction, while periodic formulation and 16 control points were used in the circumferential direction, n+1=6 and m+1=16 (96 control points in all). Again the z-values of the functions are represented as phase maps just for a better visualization.
In the first example a generic quadrangular domain is considered, in this case it is easy to define a bilinear transformation function to pass from a domain to the other. Thus, as reported in Fig. 3(a)
, it is possible to calculate the functions f
that allow the transformation between the two coordinate systems by simple geometric considerations as those reported in Fig. 3(b)
As a numeric example let us consider the real domain reported in Fig. 3(a)
and we suppose that ξ, and η vary between 0 and 1. By assuming the following values for the end point coordinates of the real domain:
the transformation functions assume the following expressions:
It is possible to recognize the coordinates of the real domain in the second set of equations (for transforming the auxiliary coordinates (ξ,η) into the real coordinates (x,y
)), while the first set of equations (for transforming the real coordinates (x,y
) into the auxiliary coordinates (ξ,η)) hides the geometric data used to generate them. This fact is due to the complexity introduced by the geometry of the problem, which however do not implies a complex algebraic formula, as it can be seen from eq. (5
Fig. 3. A generic quadrangular domain: a) the geometric transformation of the coordinate system; b) the geometric constructions for the calculation of the transformation equations.
In the second numeric example let us consider a convex domain whose boundary is a polyline, as that reported in Fig. 4(a)
. In this conditions it is not possible to define a single function able to transform all points inside the domain, but a single transformation can be defined on a number of sub-regions. This domain can be divided into triangle simply by fixing a point inside the domain and joining this point with each vertex, as shown in Fig. 4(a)
; in this way N
triangles are obtained if N
is the number of the edges of the domain. Finally, according to the Fig. 4(b)
each triangle can be transformed into a circular sector by the following general formula:
i) and (x,y) are the coordinates of the triangle and the generic coordinates in the real domain respectively, (r,θ) are the generic auxiliary coordinates that vary in the ranges [0,1] and [θ
Fig. 4. A generic convex domain whose boundary is a polyline: a) the geometric transformation of the coordinate system; b) the geometric constructions for the calculation of the transformation equations.
It is worthwhile to mention the possibility to apply the same approach in the case of the same kind of domain with a hole. In this condition, by approximating the hole with a polyline with a number of edges equal to the number of the boundary edges, it is possible to obtain again a number of sub-regions each of which can be treated either as a rectangular regular region or as a circular sector starting form a non zero radius. This case is shown in Fig. 5
, the formula are not reported for the sake of brevity.
Fig. 5. A generic convex domain whose boundary is a polyline with a hole: a) partition into sub-regions of the domain in the real coordinate system; b) transformation of a single sub-region into a regular rectangular region; c) transformation of a single sub-region into a regular circular sector.
3. Application of the procedure
The present section reports two examples of application of the fitting procedure described in the paper. In the first example the procedure was applied on analytical data in order to identify some guidelines for the choice of the order and of the number of control points of the B-spline surface. Subsequently the same case study was used with the aim of evaluating the capability of the procedure to work in presence of noise. In the second example the procedure was applied on experimental data obtained by phase shifting speckle interferometry.
Both the examples consider the out-of-plane displacements of a plate subjected to flexural loading condition, because often for these mechanical configurations these displacement components are enough to give an exhaustive description of the stress and strain fields; it is obvious that the fitting procedure can be applied to any displacement components or even to any bi-dimensional field. For the sake of brevity rectangular domains were considered in both examples.
3.1 Application to analytical data
The first case study is a square thin plate supported along the edges and loaded by a uniform pressure. For this problem the out-of-plane displacements can be evaluated analytically by the general following expression [21
21. R. Szilard, Theory and analysis of plates: classical and numerical methods (Prentice-Hall, 1974).
) is the out-of-plane displacement at the generic coordinate, p
is the pressure applied on the plate, E
are the Young’s modulus and Poisson’s ratio of the material, a
are the length of the edges, t
is the thickness, n
are two indices used to sweep the odd integer numbers between 0 and infinity. The loading configuration and a typical phase map of the out-of-plane displacements are reported in Fig. 6
The first task of this section consists in defining a procedure for choosing the order of polynomials (h,k
) and the number of control points (n
+1) necessary to fit accurately the simulated experimental data. By eq. (7
) the displacements were evaluated on a grid of 300x300 elements, that is a proper number of information gathered by an experimental test. Then these data were fitted by varying the parameters (h,k
) and (n
+1) and for each values assumed by the fitting parameters the maximum error (Emax
) and the standard deviation (σ) were evaluated. Because of the symmetry of the problem the same parameters were assumed along the x
- and y
- directions (h
). Furthermore, as said in section 2, because of each blending function spreads throughout (k
+1) control points the order of polynomials fixes the minimum number of control points.
reports the values obtained for Emax
and σ and provides quantitative information to choose the fitting parameters. For example if a maximum error less than 0.1% is desired it is necessary to use 10 control points for h
=3 or 7 control point for h
=4, while more than 12 control points are necessary if a linear approximation (h
=2) is adopted.
Geometry of the loading configuration: a) square thin plate (L/t=100, ν
=0.3) simply supported along the edges and loaded by a uniform pressure; b) simulated wrapped phase map of the out-of-plane displacements obtained by eq. (7
It is worthwhile making some further remarks about the choice of the degree of the polynomials. The degree to be adopted depends on the shape of the 2D field to approximate: if the field does not present change of concavity a parabolic approximation should be enough, while the presence of flexes requires the use of at least a cubic approximation. If more than one flex is present and they are close it is necessary to increase the degree of polynomials or to insert at least one control point between a pair of flexes. Furthermore a linear approximation (k=2) is not in general a good choice, especially when a displacement component must be fitted, because it provides very poor strain fields which must be evaluated by derivative operations.
The second step consists in evaluating the capability of the procedure to be noise tolerant. To perform this task an additive and a multiplicative noise with flat and gaussian distribution were superimposed to the data generated analytically. Then the fitting procedure was applied to the noisy data and the performance of the fitting was evaluated by calculating the standard deviations reported in Tab. 2
In particular the standard deviations are evaluated on: a) the difference between the fitted surface and the data without noise, σ1; b) the difference between the fitted surface and the data with noise, σ2; c) the noise as is, σ3. In the simulations the level of noise was increased by increasing the standard deviation of the random data, which was varied between 0.02/30.5 and 0.10/30.5 by 5 steps. These values, calculated as r/30.5, are adopted in order to have a flat distribution in the range [-r, r] which results to have the same standard deviation of gaussian distribution. Obviously for the additive noise the mean value of both distributions is 0, while for the multiplicative noise the mean value is 1. The parameter r is expressed as a fraction of the maximum displacement observed, in the simulations of the noise its value is varied until 10%, which represents a particularly severe noise.
By observing the standard deviation reported in Tab. 2
it is possible to state that the fitting procedure is able to filter the noise very efficiently. The differences between the fitted surface and the data without noise are very small, in fact the value σ1
is more than one order of magnitude less than σ2
. For the additive noise the value σ2
is practically equal to σ3
, and it means that the noise is not able to bring to fail the procedure. For the multiplicative noise this fact does not happen because in this case it is necessary to take into account the distribution of the displacement; nevertheless the standard deviation σ1
is still very small.
Table 1. Maximum error and standard deviation (expressed as percentage of the maximum displacements) obtained by varying the order of the polynomial (h,k) and the number of control points (n+1,m+1).
| | |
Table 2. Standard deviation evaluated on the data with various type and levels of noise generated numerically.
| | |
3.2 Application to experimental data
The experimental case study is an isotropic plate subjected to the flexural loading configuration reported in Fig. 7(a)
: the plate is supported on the three points S1
and loaded by a punctual force F. The experimental data are acquired only on a reduced area of the specimen because the loading fixture creates patches of shade. The specimen is a square with an edge of 50 mm, while the observed area is a square chosen in the middle of the specimen with an edge of about 35 mm, the experimental data consists of 400x400 pixels. The specimen is inserted along one arm of a Michelson interferometer capable of measuring out-of-plane displacements, according to the layout reported in Fig. 7(b)
. The reference beam is obtained by a rough surface fixed to a three degrees of freedom PZT actuator able to perform a uniform translation along its axis (necessary for phase shifting) and two tilting around an axis belonging to a plane parallel to the reference surface.
Fig. 7. a) Geometry of the loading configuration: S1, S2 and S3 are the support points, F is a punctual load applied in the middle of the specimen. The area observed by the interferometer is emphasized by projecting on it an experimental fringe pattern. b) Layout of Michelson interferometer for measuring out-of-plane displacements.
As for the analytic case study it is necessary to define an index of performance of the fitting algorithm when it is applied to experimental data. Unfortunately a known deformation field it is difficult to generate straightforwardly, apart from elementary loading configurations for which a complex fitting operation would be unnecessary.
Another kind of reference solution consists of a numerical simulation of a linear elastic body subjected to mechanical loads, in fact the solutions of these problems are nowadays well assessed and perfectly reliable. Nevertheless a very accurate simulation requires the exact knowledge of the mechanical configuration, i.e. location of constraints and loads, geometry, properties of material. For example the assumption of a punctual load or a punctual constrain is an approximation which can lead to substantial differences in the proximity of these zones between numeric and experimental data. Even if all these information are available with precision other sources of uncertainty are present that can produce systematic errors. Except for the electronic noise, which introduce a random error easily filtered by the fitting procedure, a rigid body motion is always present because of the compliance of the loading device. In fact as far as the deformation field is known it is difficult fix exactly the deformation, because also the constraints, which typically are zero displacement points, are compliant and spatially extended.
For all these reasons, excluding the numerical simulation, the only way to dispose of a known deformation field as accurate as possible consists of a uniform tilting of the specimen obtained by imposing a tilting to the reference surface fixed on the PZT actuator. According to this approach for an imposed tilting given to the reference surface we obtain a phase map whose exact approximation is a plane. By fitting this phase map with a plane of any orientation three parameters are obtained out of the number of information gathered experimentally (usually hundreds of thousands, 400x400 in the example reported in the present section). If the error superimposed to the measurements has a zero mean distribution, as the errors due to the electronic noise and to the decorrelation of the speckle patterns can be assumed, the best fit plane can be considered the exact deformation field. By this plane and by the experimental data the standard deviation can be evaluated and considered as a synthetic index of performance, that can be compared with the standard deviation evaluated by a best fit surface and the corresponding experimental data. The standard deviation obtained on the deformation field of known shape can be seen as a synthetic parameter for identifying the best fitting attainable, if a fitting performed on a deformation field of unknown shape is characterized by the same standard deviation it means that the two experimental conditions (known and unknown deformation field) are the same in term of noise (essentially due, for a phase shifting speckle interferometry technique, to the electronic noise and to the decorrelation).
By this approach two different levels of load were applied to the specimen and a bi-cubic (h
=4) B-spline surface with 7x7 control points were used, results are reported in Fig. 8
and Fig. 9
. Figures 8(a)
report the experimental data without any manipulation or fitting operation, while Fig. 8(b)
and Fig. 9(b)
are the fitted data represented as phase maps in order to be compared more easily with the experimental data. Figures 8(c)
are the histograms of the error evaluated as the difference between the experimental and the fitted data for the known (continuous line) and unknown (dots) deformation field. In both cases the two error distributions, which are symmetric and have a mean equal to zero, are slightly different but they have the same standard deviation. In the upper-right corner of Fig. 8(c)
and Fig. 9(c)
there are the experimental phase maps obtained by tilting the reference surface which determine the same standard deviation of the fitted surface, while in the upper-left corner are reported the standard deviations of known (σk
) and unknown (σu
) deformation field. These quantities are expressed in radians because all the quantities are evaluated as phase variation, as also the errors used for plotting the histograms of Fig. 8(c)
and Fig. 9(c)
; for the scope of the present paper it is not necessary the exact calibration in terms of the physical units of the measured quantity.
Fig. 8. Low level of load: a) experimental data; b) fitted data; c) histogram of the error for known (continuous line) and unknown (dots) deformation field, experimental phase maps of the known deformation field in the upper-right corner, standard deviation of the known (σk) and unknown (σu) deformation field in the upper-left corner.
Fig. 9. High level of load: a) experimental data; b) fitted data; c) histogram of the error for known (continuous line) and unknown (dots) deformation field, experimental phase maps of the known deformation field in the upper-right corner, standard deviation of the known (σk) and unknown (σu) deformation field in the upper-left corner.
The paper presents a fitting procedure able to provide a smooth approximation of a bi-dimensional distribution of a general quantity. The procedure is based on a global approach according to which a linear model was defined by the geometry of the surface and by the approximation functions called blending functions. Then this model is inverted and applied on the experimental data to evaluate the best fit parameters, which represent the weights of the blending functions.
The proposed procedure, that can be applied on domains of several shapes, is based on a B-spline formulation by which the blending functions are built by piecewise polynomials defined in sub-regions of whole domain; the order of polynomials and the number of blending functions can be chosen independently from each other. In the paper a procedure to evaluate these parameters and the equations for studying some type of domains are proposed.
Finally the fitting procedure was tested on two sets of experimental data obtained by an out-of-plane speckle interferometer. The performance of the procedure was evaluated by calculating the standard deviation of the difference between the experimental and the fitted data. This standard deviation was compared with that calculated by fitting a known deformation field, obtained by imposing a uniform tilting to the reference beam of the interferometer.