Axial stereo-photogrammetry for 360° measurement on tubular samples

ARTICLE IN PRESS

Optics and Lasers in Engineering 45 (2007) 637–650

Axial stereo-photogrammetry for 3601 measurement on tubular samples Katia Genovesea,, Carmine Pappalettereb a

Universita` degli Studi della Basilicata, Dipartimento di Ingegneria e Fisica dell’Ambiente, Viale dell’Ateneo Lucano 10, 85100, Potenza, Italy b Politecnico di Bari, Dipartimento di Ingegneria Meccanica e Gestionale, Viale Japigia 182, 70126, Bari, Italy Available online 10 October 2006

Abstract This work presents a stereo-photogrammetry (SP) based procedure to perform whole-body measurements on tubular samples. Such a system is designed for future applications to the study of vascular wall mechanics. The use of a concave conical mirror surrounding the specimen makes it possible to capture the reflected 3601 surface with a single snapshot moving neither cameras nor object. Then, according to 3-D computer vision principles, a stereo camera system retrieves control points depth information from image-pairs of the investigated surface. An axial-SP arrangement is selected since is more suitable for this specific application than the more popular lateralstereo model. In this paper, particular emphasis is given to a formulation taking into account even small camera misalignments. A calibration process based on optimization concepts is used together with a feature-based matching algorithm to efficiently find correspondence between highly distorted images reflected by the conical mirror. Both theoretical and experimental analyses on calibrated samples demonstrated feasibility and accuracy of the proposed procedure. r 2006 Elsevier Ltd. All rights reserved. Keywords: 3601 measurements; Stereo-photogrammetry; Tubular samples

1. Introduction Understanding and modeling vascular wall mechanics is a primary issue in the study of circulatory diseases. Although theoretical and numerical studies on arterial compliance are increasing in numbers, relatively little work has been documented on the use of non-invasive full-field imaging techniques for monitoring 3-D vascular wall deformations. Usually, 2-D video dimension analyzer (VDA) systems recover length and diameter variations during inflation/extension tests by tracking position changes of two couples of markers applied longitudinally and transversely onto blood vessel surface [1]. Strain determination then relies on the assumption of axisymmetric deformations. However, because of the highly anisotropic and inhomogeneous structure of vascular tissue, full-field evaluation of blood vessel deformation map should be achieved to properly model its mechanical response. Bi-plane video systems for 3-D measurements on Corresponding author. Tel.: +39 0971 205146; fax: +39 0971 205160.

E-mail address: [email protected] (K. Genovese). 0143-8166/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.optlaseng.2006.08.009

vascular segments appeared just very recently [2]. However, still no whole-body reconstruction has yet been reported. To address issues mentioned above, a stereo-photogrammetry (SP) based procedure for reconstructing whole 3-D shape of tubular samples is proposed in this paper. Besides other interesting applications, in fact, the present system can indeed be used for recovering whole deformation map of arterial wall under inflation-extension tests. SP falls in the variety of full-field non-contact optical methods developed for 3-D topographical data measurements [3]. Most commonly, when 3601 profiling of an object is required, measurements are carried out in two steps: first, object images are grabbed from different views (using a multi-camera system or a camera-turntable arrangement); then, all views are merged together into a main point cloud [4]. In this research, principles of SP and radial metrology are combined together into a system capable to perform 3-D whole-body measurements with a single snapshot. A concave conical mirror coaxial to the tubular sample captures the whole reflected surface of the specimen moving neither acquisition system nor object. Images are

ARTICLE IN PRESS 638

K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

recorded with a stereo-camera system and then processed, following SP principles, to locate the virtual image of a pattern applied onto the sample surface. Finally, once cone mirror geometry and position have been assessed through calibration procedure, a straightforward specular transformation recovers object shape from virtual point location map. Similar arrangements based on radial metrology concepts have already been used for obtaining panoramic views over a wide range of resolutions. Panoramic annular lens [5] and concave mirrors [6], for example, can be used for highly sensitive measurements on internal and external axi-symmetric surfaces with digital correlation and speckle interferometry, respectively. At lower spatial resolution, curved mirrors are extensively used for building SP catadioptric sensors in robot-vision applications [7]. However, on authors’ knowledge, a radial-SP system for measuring large deformations on tubular sample has not yet been reported in literature. To perform SP measurements [8] it is necessary to properly process image-pairs of the area of interest as captured by two cameras. The procedure for estimating depth from homologous points coordinates is defined by the specific camera model. In this study, we choose an axial-stereo (AS) arrangement (where the two cameras share a common optical axis). In fact, this set-up has, especially for this specific application, some significant advantages over the most widely used lateral-stereo (LS) model [9]. The paper describes as camera model and image processing algorithms are selected for optimally dealing with the highly distorted and in-homogeneous images reflected by the conical mirror. A formulation of AS model accounting cameras unknown misplacements is derived and implemented within a calibration procedure relying on optimization concepts. Both simulations and experimental tests served to evaluate validity and accuracy of the proposed method. System proved itself able to provide very accurate results in spite of inherent drawbacks of axial-stereo configuration. Whole 3601 shape of calibrated samples have been reconstructed from single snapshot acquisitions thus avoiding time-consuming multiple acquisitions and merging procedures. 2. Theoretical This section briefly outlines the AS model and the issues connected to its use in omnidirectional view systems. 2.1. The 3601 stereo-photogrammetric system 2.1.1. The axial- stereo model Stereo-vision [8] includes the following main phases:





image acquisition; features extraction;





image matching (solving the correspondence problem); depth determination from image-pairs according to the camera model equations.

In SP, the object is imaged by two cameras [10]. Many different camera models have been used for specific applications: converging or parallel camera optical axes; camera translation along optical axis; and rotating camera optical axis. The most widely used stereo camera arrangement found in literature is the lateral-model which simulates the human visual system. There is a variety of techniques based on this camera model, see [8,11] for a survey. Once the camera model has been defined, and two images of the same scene obtained, it is necessary to locate the projections of the same physical point onto camera sensors (matching problem). In fact, according to the formulation of the specific camera model, depth information is derived from homologous points disparities (i.e. differences in the xy coordinates of corresponding points in the images planes). Image matching techniques can be divided into three major categories: area-based, featurebased and hybrid. Area-based algorithms use the local gray level distribution information to find the best match between a pair of images [12,13]. Feature-based algorithms compare specific characteristics extracted from imagepairs, such as centroids, edges, lines, corners and other regular shapes [14]. Hybrid approaches properly combine both area-based and feature-based algorithms. Fig. 1a shows the optical set-up including a concave cone mirror, a beam splitter and two CCD cameras proposed in this research to perform 3601 measurements. To image the sample mounted coaxially to the conical mirror, an SP rig with two cameras sharing a common optical axis is adopted. Fig. 1b illustrates the rationale behind AS model. The concave mirror surrounding the specimen forms a perfect virtual image point P behind mirror for each world point P0 on the sample surface within the field of view of the acquisition system. According to the scheme of Fig. 1b, the two cameras should be aligned with the axis of the cone mirror. However, to avoid view occlusion, one camera is placed away by using a beam splitter yet practically maintaining the alignment with optical axis. It can clearly be understood from Fig. 1b, how image point position on each camera sensor changes according to different viewpoints and focal lengths and is directly related to the world point location. Such disparity permits to calculate world coordinates of point P from local coordinates information on point-pairs PN and PF in the two camera views. Once virtual location map is reconstructed based on SP equations, a reflection allows to recover position of the object surface point P0 . The AS arrangement is chosen for dealing with the highly distorted images reflected by conical mirror. Briefly, for an AS model, the two views are easier to be compared than those obtained from the more common LS

ARTICLE IN PRESS K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

639

Fig. 1. Schematic of the optical arrangement for 3601 measurements on tubular samples (a); axial-stereo model (b).

arrangement. This is because equivalent points fall in a smaller ROI and undergo similar distortion. Conversely, the adoption of the AS model requires special cares to be taken in selecting a robust matching algorithm for treating differently sized images. As it will be illustrated for this special application, feature-based approaches appear more appropriate than area-based algorithms. Retrieval of depth information from 2-D images captured by near and far cameras (whose relative parameters will be indicated in the rest of the paper with the N and F subscripts, respectively) is usually done by referring to the simple scheme sketched in Fig. 1b (see, for example, [7,9]). A projective geometry (PG) model is adopted to describe image formation process: cameras are considered as ideal pin-hole cameras. Image planes related to projection centers C N and C F are indicated as  iN  and  iF  respectively; K N and K F are the centers of the two sensors. The global coordinate system center coincides with the near projection centre C N . The two cameras are separated by the system baseline D in the z-direction. Distances C N K N and C F K F will be indicated as f N and f F , respectively. It should be noted that, when object distance from projection centre is large, f N and f F can be considered equal to nominal focal lengths. If this condition is not satisfied, more rigorous equations should be used to evaluate distance between lens system center and image plane. PN ðxN ; yN Þ and PF ðxF ; yF Þ are the projections of a point PðxP ; yP ; zP Þ in the local systems  iN and  iF ; PPN ðxPN ; yPN ; 0Þ and PPF ðxPF ; yPF ; DÞ, respectively, are their counterparts in the global system. If perfect alignment of cameras is achieved, it will hold xk ¼ xPk and yk ¼ yPk ðk ¼ N; F Þ. From similarity between triangles in Fig. 1b: 8 x y < xPNP ¼ yPNP ¼ fzNP ; (1) : xxPFP ¼ yyP ¼ ðzPfþDÞ : PF

F

This leads to obtain: 8 > > > <

xP ¼ 

D fN fF xPF xPN

 ; yP ¼ 

D fN fF yPF yPN

;

; zP ¼ zP_X ;Y ;d ¼    D   > > xPN yPN d PN fF > : xPF ; yPF ; d PF f N 1

(2)

where d Pk ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2Pk þ y2Pk with k ¼ N; F .

By combining relationships (1) we obtain yPF y ¼ PN . xPF xPN

(3)

This equation states that lines l N and l F , respectively, connecting K N ! PN and K F ! PF are parallel. This is a particular feature of the AS model: each point lying on the plane passing through P, CN and CF (called epipolar plane) will fall on two parallel lines l N and l F called epipolar lines, respectively, passing by points J N K N and J F K F (called epipoles and obtained as the intersection of C N C F with image planes). We will return later on AS epipolar geometry showing as it can be used to simplify and enhance matching algorithm. It should be stressed that conditions (2) occur only when a perfect alignment of cameras is ensured. In fact, next paragraph will demonstrate how even small misplacements of camera sensors may strongly affect measurement accuracy. For this reason, some generalized AS model such that proposed in this study should always be adopted. 2.1.2. The generalized axial-stereo model A generalized AS model is represented in Fig. 2. Each sensor  ik  whose coordinate system is originally aligned with the global system, is rotated first by the tilt angle ak about the x-axis, then by the pan angle bk about the

ARTICLE IN PRESS K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

640

Fig. 2. Schematic of the generalized axial-stereo model (special case withaN ¼ aF , bN ¼ bF andgN ¼ gF ).

rotated y-axis and, finally, by the swing angle gk about the twice rotated z-axis. In the most general case, no information is given on cameras position. Since the proposed system is specifically designed for measuring shape and deformation without requiring any data on absolute position, we can arbitrarily locate global coordinate system without losing in generality. Here, we choose the 3-D world coordinate system where z-axis coincides with the C N C F and coordinate origin is located in the projection centre C N of the near camera. Therefore, any translation of cameras with respect to the global system are null (except for O0 z-coordinate). Hence, relationships between world and local coordinate systems become much simpler. Likewise the ideal case, PG relationships between global coordinates ðxP ; yP ; zP Þ of point P and its projections PPN and PPF into camera systems  iN and  iF can be expressed as (see Fig. 2): 8 <

xP xPN

P ¼ yyP ¼ jzzPN j; PN

: yyP ¼ xxPFP ¼ ðjzðzPFP þjjDjDjÞjÞ :

(4)

PF

In Eq. (4) the terns ðxPN ; yPN ; zPN Þ and ðxPF ; yPF ; zPF Þ, respectively, represent the global coordinates of points PN ðxN ; yN Þ and PF ðxF ; yF Þ. In this work a scale factor of ffi0.006 mm/pixels has been adopted for a 7.7 mm  6.2 mm camera sensor with 1280  1024 pixels. PN and PF coordinates can be calculated as: xPk ¼ xk cos bk cos gk þ yk cos bk sin gk  f k sin bk , yPk ¼  xk ðsin ak sin bk cos gk þ sin gk cos ak Þ

þ yk ðcos ak cos gk  sin ak sin bk sin gk Þ  f k sin ak cos bk , zPk ¼ xk ðsin ak sin gk  cos ak sin bk cos gk Þ  yk ðsin ak cos gk þ cos ak sin bk sin gk Þ  f k cos ak cos bk  D,

ð5Þ

where D ¼ 0 if k ¼ N. Hence, it follows: 8 < xP ¼ : zP ¼

D

ðððjzPF jjDjÞ=xPF ÞðjzPN j=xPN ÞÞ zP_X ;Y ;d ¼

;

yP ¼

D

ðððjzPF jjDjÞ=yPF ÞjzPN j=yPN Þ

jDj

ðððjzPF jjDjÞ=jzPN jÞðxPN =xPF Þ;ðyPN =yPF Þ;ðd PN =d PF Þ1Þ

:

(6) It should be noted that, whatever reasonable accuracy in system alignment is reached, it is practically impossible to avoid even very small misplacements with respect to the ideal scheme of Fig. 1b. Therefore, when Eq. (2) is used in place of Eq. (6) errors in topography reconstruction will occur. To evaluate the degree of inaccuracy entailed by simplified formulation, virtual image pairs PN ðxN ; yN Þ and PF ðxF ; yF Þ of a tilted plane have been reconstructed for an AS system with ak ¼ bk ¼ gk ¼ 0:01 , f k ¼ 50 mm and D ¼ 300 mm following a ray-tracing procedure. Simulations served to evaluate influence of system geometry on topography reconstruction errors without including the effect of optical aberrations and matching procedure. Computed control points of synthetic plane have been graphed (from left to right) in Figs. 3–4:






green plot: ðxP ; yP ; zP_X Þ obtained from Eq. (2); blue plot: ðxP ; yP ; zP_X Þ or equivalently ðxP ; yP ; zP_Y Þ and ðxP ; yP ; zP_d Þ obtained from Eq. (6); red plot: ðxP ; yP ; zP_Y Þ obtained from Eq. (2).

ARTICLE IN PRESS K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

Fig. 3. 3-D representation of virtually reconstructed planes for an axialstereo system with ak ¼ bk ¼ gk ¼ 0:01 , f k ¼ 50 mm, D ¼ 300 mm by using Eq. (2) (green and red plots) or Eq. (6) (blue plot) (tails have been cut and surfaces have been relatively shifted along z-axis for the sake of clarity).

1025

z (mm)

1015 1005 995 985 975 -40

-30

-20

-10

0

10

20

30

40

x (mm) ZP vs. X

ZP_X vs. X ZP_Y vs. X

Fig. 4. x–z view of control points for surfaces in Fig. 3.

Figs. 3–4 provide an overall representation of reconstruction error trend depending from all parameters included in (6). It is to be noted that using Eq. (2), when ak ; bk ; gk a0, erroneously leads to have zP_X azP_Y azP_d . In particular, Fig. 4 shows as the errors in evaluating zP_X and zP_Y increase when x and y tend to zero, respectively. It will be shown later in the paper that the condition zP_X =zP_Y ¼ 1 can effectively be used to build an objective function for evaluating ak ,bk and gk within the optimization-based procedure for system calibration. 2.1.3. The epipolar geometry for the generalized axial stereo model Computational procedures for stereo matching often take advantage of a geometric constraint stating that corresponding points of two image pairs will always fall on the above-defined epipolar lines. The use of the epipolar constraint reduces the search space from two dimensions to one, hence producing a tremendous saving in the computation time required to find the matching solution. Usually,

641

to avoid image rectification (the operation that transforms the original image pair into new digital images where epipolar lines become the rows) and hence to simplify matching procedure, stereo cameras are arranged in a lateral-parallel configuration. This scheme generates horizontal epipolar lines thus reducing stereo search space to one dimension along the horizontal axis. However, in general, use of epipolar constraint will require the precise knowledge of the relative positions of the two sensors. Such information can only be obtained by performing system calibration. The matching algorithm developed in this work takes advantage of the information of epipolar lines. One of the major advantages of the present setup, in fact, is the simple epipolar geometry. In the ideal AS configuration, when each component of the system is perfectly aligned (Fig. 5), projections of points belonging to the pencil of planes through the z-axis will fall on radial lines passing through epipoles J N K N and J F K F . Since epipolar lines are the traces of the plane through P, C N and C F in the image planes, in the ideal AS configuration angles between epipolar lines on each camera sensor are preserved. Using this constraint in the general situation portrayed in Fig. 2 is less straightforward then in the ideal case since, because of rotations ak ,bk andgk , corresponding families of radial lines result differentially shifted (i.e., J N aK N and J F aK F ) and rotated in the sensor plane. To retrieve epipolar geometry of the generalized AS model we can start by finding the equation of a generic epipolar plane pP passing through the world points PðxP ; yP ; zP Þ,C N and C F . In the local system  ik , equation of plane pP : yP ¼ xP tan w can be written as pkP : xk ðsin ak sin bk cos gk þ sin gk cos ak þ cos bk cos gk tan wÞ þ yk ð cos ak cos gk þ sin ak sin bk sin gk þ cos bk sin gk tan wÞ þ zk ð sin a cos b þ sin b tan wÞ þ f ðsin ak cos bk  sin bk tan wÞ ¼ 0.

ð7Þ

Fig. 5. Schematic of an axial-stereo system showing epipolar geometry (ideal case with ak ¼ bk ¼ gk ¼ 0).

ARTICLE IN PRESS K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

642

The line of intersection between pkP and  ik : zk ¼ 0 represents the epipolar line rkP where projection of point P will fall: rkP : xk ðsin ak sin bk cos gk þ sin gk cos ak þ cos bk cos gk tan wÞ þ yk ð cos ak cos gk þ sin ak sin bk sin gk þ cos bk sin gk tan wÞ þ f ðsin ak cos bk  sin bk tan wÞ ¼ 0.

ð8Þ

From Eq. (8) it is possible to see that slopes of corresponding epipolar lines on the two camera sensors are different in the general case where aN aaF , bN abF andgN agF . In fact it holds tan dk ¼

sin ak sin bk cos gk þ cos ak sin gk þ cos bk cos gk tan w , cos ak cos gk  sin ak sin bk sin gk  cos bk sin gk tan w

(9) where tan dk ffi tan g if ak ; bk ; gk ! 0. Furthermore, epipolesJ N e J F (i.e. points with xP ¼ yP ¼ zk ¼ 0) are not centered in the camera planes as it occurs for the ideal case but have local coordinates: f tan ak f tan ak k sin gk ; k cos gk J k ¼ f k tan bk cos gk  k cos bk cos bk

þf k tan bk sin gk ; 0 . ð10Þ In conclusion, once angles ak , bk and gk are known from calibration, Eqs. (9) and (10) allow us to recover epipolar geometry of the AS system to be used for simplifying and speeding up the image matching procedure. 2.2. Implication of the use of an axial stereo rig into an omnidirectional system 2.2.1. Strength points and limitations From the technical point of view, the AS model possesses some significant advantages over the traditional LS model. A comprehensive description of differences between stereo camera models can be found in Ref. [9]. Briefly, it has been proven that, for the same range of depth values, correspondence ROI in the AS model is always less than half of its LS model counterpart. Moreover, in the AS system, the common field of view of the two cameras coincides with the field of view of the camera near. Thus, when image features are selected from the front image, projections of the same scene features are guaranteed to exist in the common view space. Moreover, for the AS model, as camera baseline D increases, depth measurement error decreases at a faster rate than for the LS model. Therefore, while in the lateral arrangement better depth accuracy is achieved at the cost of a larger number of ambiguous matches, the camera displacement in the AS model may be increased for smaller depth measurement errors, without affecting the number of correct matches. However, the main reason considered in this study behind the choice of an AS model instead of the more common LS model comes from the use of a conical mirror

Fig. 6. Right image grabbed in a lateral-stereo configuration. Left image is symmetric with respect to the dashed axis. The picture shows also illumination non-uniformity for the lateral model.

for obtaining the 3601 vision. Fig. 6 shows the virtual image of a cylindrical surface coaxial to the mirror with a regular pattern of dots applied onto it. The figure shows how differently distorted images are observed if LS arrangement is adopted. In fact, left and right images exhibit an higher pattern frequency in the circumferential direction inside the region closest to the camera and vice versa. This results in very differently sized search areas to be compared. Therefore, possible failure of matching procedure, especially in case of area-based algorithms, can occur. Moreover, since for improving accuracy baseline should be increased the most as possible (yet trying to keep low the number of unambiguous matches), it can clearly be understood that problem of occlusion arises for LS arrangement with the sample in the present position. Similar limitations, although with less significant occlusions and with the advantage to have epipolar lines parallel to the camera sensors, would be encountered with the parallel-lateral camera configuration. Finally, since images form beyond a very reflective surface, illumination source position is more critical for a LS model (see Fig. 6) with respect to the AS configuration where an annular light source coaxial to the system axis is enough to ensure uniform illumination and contrast (see Fig. 12a). However, some limitations also exist for the AS model: 2.2.1.1. Multiresolution. An entity is not represented by the same number of pixels in both images. In fact, although the two views are easier to compare that for more common lateral vision modes since searching area is smaller, special care has to be taken in selecting a robust matching algorithm for treating differently sized images without resampling. Unfortunately, accuracy in depth measurement increases when relative displacement between

ARTICLE IN PRESS K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

643

Fig. 7. Resampling of images for solving correspondence problem with area-based correlation algorithms. A speckle pattern is applied on a cylindrical surface coaxial to the conical mirror (not visible in the cropped images) (a); distortion along the circumferential direction of a regular pattern made up of square markers (b).

cameras also increases. This fact implies a greater difference in size for the two image pairs. Hence, a tradeoff must be found between overall accuracy and computational efficiency of the matching algorithm in order to obtain reliable results. As it can be easily understood, for AS photogrammetry, and especially for this particular application where control patterns on the sample surface are reflected by a conical mirror, feature-based algorithms usually perform better than area-based methods. The reasons for this statement can be explained as follows. When measuring topography, the primary task to accomplish is to compare two corresponding areas that have different sizes in the two image pairs. If the largest image (the front image captured from the near camera) is resampled to be comparable with the smallest one (the far image), preservation of gray-level distribution (Fig. 7a) is not always ensured. This could lead to failure of areabased matching approach. Should even this problem be brilliantly overcome by using highly efficient correlation algorithms together with additional geometric and scene constraints, additional problems may arise when deformations instead of topography must be measured with the omnidirectional system proposed in this work. As it can be seen from Fig. 7b, in fact, the degree of circumferential distortion of the image of a marker applied on the sample depends on the axial position of the marker itself (radial direction on the image). Thus, when deformation has to be measured, the control area in the undeformed image not only covers a different number of pixels than in the deformed image but it could also be distorted differently. It has been proven [15] that accuracy in measurements using area-based digital correlation algorithms drops rapidly when correlation is lost. This means that, to get the best measurements as possible, it is important to compare the most similar image fields as possible. In this case, due to

presence of large distortions, finding the correspondence could be very cumbersome. 2.2.1.2. Central invalid area. There always will be a set of corresponding points within a certain radius from the center of the front image that leads to an erroneous depth value evaluation [9]. Referring to the scheme of Fig. 1, in fact, for any given error e in pixels, it is possible to find the radius rN of the circular area centered in the image  iN  whose corresponding points possess a disparity less than  and must thus be considered invalid for depth estimation. Radius rN can usefully be computed in the experiment planning phase to evaluate the maximum coordinate z¯ P of points belonging to the object in order to obtain an invalid area within the smaller base of conical mirror (area A in Fig. 14a). From Eq. (1), assuming f N ¼ f F for simplicity, it is possible to calculate disparity as: dN  dF ¼

Dd 2N . fd P þ Dd N

(11)

Then, considering only points at z ¼ z¯ P , it holds: dP ¼

z¯ P d N . f

(12)

By combining Eqs. (11) and (12), it is possible to obtain the minimum radiusrN of the invalid area: rN ¼ 

ð¯zP þ DÞ . D

(13)

Hence, in the model, image points in the image  iN  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi for which 0p x2N þ y2N prN must not be considered in the image matching process. As example, considering a relative displacement between cameras D ¼ 254 mm, an object distance from C N equal to z¯ P ¼ 1000 mm and a camera

ARTICLE IN PRESS 644

K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

sensor with 0:006 mm=pixel, for  ¼ 1 pixel, the invalid area radius is rN ¼ 0:0296 mm ffi 5 pixels. It should be noted that this can be seen as the theoretical limit for invalid area radius. In fact, it will be shown in Section 4 that, even outside of the invalid area, the measurement efficiency drops if x; y ! 0 (i.e. towards the center of the image) since, for a given error in pixels in establishing point location, the smaller the x or y coordinate, the greater the percentage of relative error in disparity calculation. 3. Experimental

reflecting surface in the global coordinate system. In fact, knowledge of position and orientation of mirror surface is necessary to perform the specular transformation when 3601 measurements on tubular samples are done. Since positions of system components are no further moved after this test, this calibration step allowed also to extract the system epipolar geometry to be used in the last set of experiments. 3601 measurement of a calibrated cylindrical sample: Whole-body topography measurements on a cylinder with a nominal diameter f ¼ 12:7 mm are performed to finally assess validity of the proposed procedure.

3.1. Experimental apparatus

3.3. Matching algorithm

The experimental apparatus of the proposed omnidirectional system is comprised of two CCD cameras (B/W, 1280  1024 pixels), a beam splitter, multi-axial stages for fine positioning, a PC-based data acquisition system equipped by an 8-bit frame grabber, an annular fluorescent lamp and a conical mirror. In this work, since real-time measurement is not entailed, only one camera is used and object has been moved back and forth along a rail to capture the two views. The conical mirror has been manufactured in steel by an ultra-precision turning machine. It is designed to possess an outer diameter of 110 mm and an inner diameter of 22 mm in order to fit the set of diameters and heights of samples to be measured. A 451 cone aperture is chosen for its peculiar optical properties. In fact, the image of a cylinder coaxial to this mirror is a flat disc. Hence, even with small depth of focus it is easy to obtain well-focused virtual images of smoothly deformed shapes of the sample. Furthermore, fairly accurate measurements can be made without need to compensate for variation in magnification at points with different depth across the image. Both camera and conical mirror are positioned on multiaxial stages for fine position adjustment. The adoption of an annular lamp allowed to obtain uniform and similar illumination conditions for both grabbed images.

In photogrammetry and remote sensing, matching can be defined as setting a correspondence between different data sets. According to the specific formulation behind the adopted stereo system, image matching serves to reconstruct the 3-D object surface from projections onto camera sensors. A great deal of effort has been spent by the computer vision community on this problem and several approaches have been reported in literature in the last decades [16]. Performing rigorous comparisons between efficiency of various algorithms for this special application is beyond the scope of this work. However, as it is explained in Section 2, we have chosen a feature-based approach for solving the correspondence problem. A feature-based algorithm compares specific characteristics extracted from image pairs, such as edges, lines, vertices and other regular shapes. In this work, a straightforward point-to-point matching algorithm is adopted. For each point of the front image, its corresponding point in the back image can be found using a centroid search criterion implemented by a Matlabs routine. This approach is particularly efficient for images containing polyhedral objects such as in the case of the adopted square dot pattern applied onto the plane and conical surface (see Figs. 8 and 12a, respectively). When 3601 measurements on tubular sample are performed, since conical mirror is used, it can clearly be understood that further problems arise when correspondence task is afforded. First, illumination conditions certainly affect efficiency of matching procedure even if feature-based approaches are rather insensitive to photometric variations. Fortunately, using an annular light source coaxial to the conical mirror (Fig. 12a) was enough for obtaining well uniformly contrasted images for both cameras. However, the main problem is still the circumferential distortion of the marker geometry generated by the specular surface (see Fig. 7b). Moving towards the most external portion of the mirror, regular geometry of reflected pattern starts to be distorted and contrast is lost. Since preliminary tests carried out by implementing simple point-to-point matching algorithms led to poor results, an ad hoc developed procedure has been implemented in order

3.2. Samples To evaluate the validity of the above-proposed system, the following tests have been performed: Surface topography measurement of a tilted plane: To calculate system parameters a perfectly flat plate is reconstructed. This preliminary test serves to verify applicability of the generalized AS formulation in the context of an optimization based procedure for system calibration. Topography measurement of conical mirror surface: After setting up the optical system the conical mirror is positioned on a multi-axial stage and a regular dot pattern printed on a thin sheet of paper is adjusted to adhere to qthe cone inner surface. Measurements of conical mirror geometry serves to retrieve the exact position of the

ARTICLE IN PRESS K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

Fig. 8. Image grabbed by near camera showing the path followed for sorting point locations onto calibration target.

to deal properly with the highly deformed virtual images of tubular sample surface. In particular, the matching algorithm developed for the last set of experiments exploits epipolar geometry. Since using of the epipolar constraint reduces search space from two dimensions to one, a tremendous saving in computation time is achieved. Moreover, in this special case, reduction of the search area to a line has the further great advantage to increase reliability of the matching algorithm in finding circumferential position of homologous point pairs in the two highly distorted virtual images. Remarkably, by adopting this approach, it is sufficient to mark the tubular sample with a pattern made up of circumferential lines instead of dots (Fig. 14). Then, by processing distribution of pixel gray levels along epipolar lines, it is easy to locate the middle points of each circumferential line section (see detail in Fig. 14) and using these coordinates together with their counterpart pairs for SP calculations. In this way, the matching algorithm is greatly enhanced by reducing of one dimension ambiguity in finding point location. It should be noted that, even if control points are sparse, measurement spatial resolution can be refined by simply increasing the number of epipolar lines investigated. However, we showed in Section 2 that the use of epipolar constraint requires a precise knowledge of the relative position of cameras. This information can be obtained only by performing a system calibration: that is evaluating parameters contained in Eqs. (9)–(10). Before running 3601 measurements it is hence necessary to calibrate the system. 3.4. System calibration The problem of estimating the 3-D position of a camera from a number of certain matches gathered considerable attention from scientific community working in the computer-vision field. A relevant example of the state-of-

645

art of camera calibration procedures can be found in Ref. [17]. In this paper we propose a straightforward approach for calibrating AS system based on optimization procedures that showed to be adequate to fit the scope of this work. Calibration has to find the unknown camera misalignments to correctly recover object topography using Eq. (6). To test calibration procedure developed here, a first set of tests are run on a tilted plane. Optically measured data are used within an optimization procedure where angles ak , bk and gk are the design variables. The goal is to find angle values that ensure the best fit between theoretical and measured geometry of the calibration target. The optimization problem has been formulated and solved with the Sequential Quadratic Programming (SQP) routine available in Matlabs. SQP [18] is a powerful optimization technique which employs quadratic approximations of cost function and linear approximations of constraints. The higher computational cost associated to each approximate subproblem is however compensated by the smaller number of design iterations required in the optimization process. In the first set of tests, optimization problem is to minimize the difference between optically measured points and the ideal plane. The problem has hence been formulated as: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi3 2 u N u1 X 4 Min Cplane ðak ; bk ; gk Þ ¼ t ðaxj þ byj þ czj þ dÞ2 5; N j¼1 k ¼ N; F ,

ð14Þ

where N is the number of control points, a, b, c, d are the coefficients of the ideal plane equation and xj ; yj ; zj are the measured coordinates of the jth control point. Very large bounds should be imposed on the optimization variables because problem (14) is actually unconstrained. To avoid numerical instabilities that may arise during the computations it is useful to introduce additional information on the problem and constraints to easily estimate initial values of the unknown parameters. Then, since the exact position of the object plane is actually unknown in the global coordinate system (i.e. a, b, c, d are trial values), optimization has been restarted from design variable values previously found. The new objective is now expressed by the condition (see Section 2.1.2): zP_X ¼1 zP_Y

(15)

which holds true only when right values of camera tilts are introduced in Eq. (6). Obviously, as it happens for any selfcalibration procedure, the method proposed here relies on the assumption that the matching algorithm finds the exact correspondence between analogous points. Finally, it should be remarked that problem (14) can effectively be built for any shape whose theoretical geometry is known. This fact allowed us to utilize the perfectly conical internal mirror surface as calibration

ARTICLE IN PRESS 646

K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

target for the 3601 measurements, following the procedure described in the next section. 3.5. Measurement procedure The whole experimental procedure proposed to perform 3601 measurements on tubular samples can be summarized as follows: (i) Camera alignment: To this aim, it is important to create a reference point by putting a marker on the object, better if in correspondence of the cone axis (see the dot on the tip of the cylindrical sample in Fig. 14a). Object positions are adjusted until the image of the marker falls approximately in the center of both views. Fine adjustment is then done by moving camera on a micrometric stage. (ii) System calibration: After the system component has been positioned, a regular dot pattern printed on a thin sheet of paper is adjusted in order to adhere to the inner surface of conical mirror. Topography assessment of the cone geometry serves to retrieve the exact position of the reflecting surface into the global coordinate system. Moreover, since positions of system components are no further moved, the calibration procedure performed on the cone allows to extract parameters ak , bk and gk . Optimization problem is built around the Ccone function in fashion of:

(v) 3601 measurements: A pattern made up of circumferential lines is applied on tubular sample surface (the sample investigated in this research is a calibrated cylinder with f ¼ 12:7 mm). The virtual image of this pattern is a series of concentric circles (Fig. 14a). Using epipolar constraint allows to simplify matching of the two highly distorted images by searching homologous points on epipolar lines (in this case, control points are individuated as the middle point of dark regions; see detail in Fig. 14a). (vi) Recovering of the 3-D shape of tubular sample: Measurement performed in step (ii) allows to retrieve both position and orientation of conical mirror surface. This permits to perform specular transformation for recovering 3-D whole external surface of the tubular sample from virtual map obtained in step (v). 4. Results and discussion 4.1. Surface topography estimation of a tilted plane The first set of experiments consisted in measuring topography of a perfectly flat inclined plane. This preliminary measurement served to test the optimizationbased procedure for calibrating the system with the AS generalized formulation. In Fig. 8 the image of the square dot pattern printed on the calibration target as it is captured by near camera is

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx cos Z  z sen Z cos x  y sen Z sen x  x0 Þ2 þ ðy cos x  z sen x  y0 Þ2 z cos x cos Z þ y sen x cos Z þ x sen x  z0

(16)

shown with the path followed to sort control point coordinates superimposed on it. As expected, experimentally measured data show the trend predicted in Section 2.1.2 (see the graph in the top of Fig. 9). The calibration 10000

Z (mm)

where x0 ; y0 ; z0 are the vertex coordinates, j is the cone aperture, and Z and x, respectively, are the tilt and the pan angles of the cone axis. (iii) Recovering of epipolar geometry: Once system parameters are determined it is possible to recover epipolar geometry of the system by using Eqs. (9)–(10). That is, it is possible to know the exact position of epipoles J k in the two images along with the value dk of epipolar lines slopes. (iv) Refinement of calibration: Once epipolar geometry is determined, it could be possible to use it for establishing a new set of correspondences by imposing epipolar constraint. More precisely, the algorithm consists of three steps: (a) establish initial correspondence using a classical feature-based matching algorithm (e.g., centroidbased criterion); (b) estimate epipolar geometry; (c) establish a more reliable new set of control point pairs imposing epipolar constraint. In this way, for each iteration, it is possible to reduce the error fraction caused by the matching process out of the global error made in calibrating the system.

 tan j ¼ 0,

1000 20

30

40

50

60

70

80

90

100

X (mm) camera tilt correction

no correction

Fig. 9. Plots of zP_X vs. xP calculated from Eq. (2) (blue) or Eq. (6)(purple) (data series have been shifted into x–z plane to be plotted into a semi-logarithmic diagram and between them to avoid cluttering).

ARTICLE IN PRESS K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

orientation have been found using the optimization procedure described in Section 3.5. This test served also to assess properly the system parameters and epipolar geometry of the optical arrangement to be used in the next set of tests. 4.3. 3601 measurement on a calibrated cylindrical sample Whole-body shape measurements on a calibrated cylinder (f ¼ 12:7 mm) mounted in order to be coaxial to the conical mirror have performed to finally verify feasibility and validity of the proposed procedure. Once AS system parameters are determined, paper sheet can be removed from specular surface and cameras are focused on the virtual image of the patterned adhesive label applied onto tubular sample. A set of circumferential thick lines is chosen as control pattern to implement a simple featurebased algorithm utilizing epipolar constraint. The matching algorithm, coded in Matlabs environment, takes as the control points the middle points of each dark region along 1150 1050 950 850

Z (mm)

procedure described in Section 3.4 is able to correctly reconstruct plane geometry (see bottom of Fig. 11) even if residual errors still occurred. It appears, by observing the blue curve in Fig. 10, that before calibration, each time the control point possesses at least one coordinate close to zero (i.e., each time the control path crosses the y- or x- axes), the zP_X =zP_Y ratio strongly deviates from 1. After calibration, (i.e., when correct values of camera position parameters are introduced in Eq. (6)), plane topography results is reconstructed more correctly. However, residual error are yet higher for points closer to the axes (green curve in Fig. 10). This fact is evident from Fig. 11 where trends of zP_X and zP_Y vs. xP are plotted. It can be easily understood that, if the optimization process found the best solution available, this residual error indicates lower efficiency of the matching algorithm when x; y ! 0 (i.e., towards the axes of the image). In fact, for a given error in pixels in establishing point location, as x or y coordinates get smaller, the percentage of relative error in disparity calculation increases. Disparity field thus obtained, in fact, is an integer-pixel dimension deriving from the discrete nature of digital images. Efficiency of matching algorithm can be increased only entailing sub-pixel interpolation schemes [19,20]. Randomly distributed error for points close to axes can be reduced by calculating zP as zP_d . Obviously, when both x and y are very small, zP is unavoidable affected by error. However, this inaccuracy is related to a small central area of the image (in this particular application useless for the 3601 measurements) and it can eventually be reduced by adopting high resolution cameras and more accurate matching algorithms.

750 650 550

4.2. Topography measurement of conical mirror surface

ZP_X/ZP_Y

Figs. 12b and 13 show results of measurement carried out to retrieve topography of the inner surface of conical mirror (Fig. 12a). A very straightforward feature-based algorithm resulted enough to achieve a fairly accurate shape reconstruction even if features (in this case square dots) are sparse and fitting procedures are needed to obtain the whole surface. Cone mirror vertex position and axis 7 6 5 4 3 2 1 0 -1 0 -2 -3

647

-35

-25

450 -5

-15

5

15

25

35

X (mm) ZP_D vs.X

ZP_Y vs. X

ZP_X vs.X

Fig. 11. Plots of zP_X , zP_Y and zP_d vs. xP for the tilted plane. Note that random errors due to a loss of efficiency of the matching algorithm for x ! 0and y ! 0can be reduced by considering zP_d . CALIBRATION PARAMETERS

100

200

300

400

500

600

700

800

αN

1.494’

βN

0.66’

γN

0.24’

αF

0.18’’

βF

-0.36’’

γF

-0.36’

point number Fig. 10. Ratio zP_X =zP_Y computed with Eq. (2) (blue) or Eq. (6)(green) (green data have been shifted by +3 for the sake of clarity). System parameter values calculated through the optimization-based calibration are listed in the table.

ARTICLE IN PRESS K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

648

Fig. 12. Cone mirror calibration: a thin sheet of paper with a regular dots pattern is applied onto the inner cone surface (a); results of the optimization procedure: points from SP measurements interpolated by the best-fitting cone (b).

1525 1520

z (mm)

1515 1510 1505 1500 1495 1490 -50

-40

-30

-20

-10

0

10

20

30

40

50

x (mm) Fig. 13. x–z profiles of the 3-D plot in Fig. 12b (measured points are in green, calculated points are in red).

epipolar lines. However, it should be remarked that, in this case, since virtual points lie on a plane z ¼ cost, near and far images are homothetic and magnification does not depend on position of each point across both image fields. Thus, it is expected that, when the sample is fairly a cylinder, also area-based matching algorithm can perform efficiently since near and far images are similarly distorted. Fig. 14 shows the near image (an in-focus image of the cylindrical sample has been pasted onto area A), along with distribution of gray level for pixels located along the yellow epipolar line before (blue plot) and after (green plot) binarization. Dots also show position of middle points (control points) for each dark region. Matching process can achieve a sub-pixel accuracy if images are preliminarily filtered (e.g., with a Gaussian filter—radius 1.5) and then up-sampled by a factor of four with a bi-cubic interpolation function. Fig. 15 shows the error map for the virtual image of the cylindrical surface computed as the difference between measured z-coordinate and its theoretical counterpart calculated through a best fitting procedure. The virtual

Fig. 14. Near camera image with superimposed epipolar lines.

plane reconstructed experimentally shows some inclination as it should certainly be expected in view of the assessed misalignment of the system. Finally, Fig. 16 shows the error made on cylinder radius with respect to the nominal value. Data are reported for each control point belonging to the 3601 surface which has been recovered from virtual plane of Fig. 15 through a specular transformation. It appears that measurement procedure left a residual error periodically distributed. 5. Concluding remarks This work presented preliminary results for a SP-based procedure to achieve 3601 topography reconstruction of

ARTICLE IN PRESS K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

649

Fig. 15. x–z view of the z-error map for the virtual image of the cylinder geometry.

Fig. 16. Map of the radius reconstruction error for the investigated cylinder surface.

tubular samples. Applications of this system can be manifold and of great interest such as full-field deformation measurement on vascular segment under inflation/ extension tests. The whole-body measurement ability of the proposed system is ensured by the combination of radial metrology and SP. Accuracy of the measurement is ensured by correcting errors caused by optics misalignments through a an optimization-based calibration procedure. The procedure was successfully verified with numerical simulations and tested on calibrated samples. Although at the present stage the system cannot be considered optimized, results are very encouraging and confirm the feasibility of the approach. In view of future developments of the system special care has to be taken for the following issues:

(i) Accuracy of the proposed method depends on many factors including quality of specular surface, signal to noise ratio in the CCD, optical system aberration and, above all, on efficiency of the matching algorithm. Sensitivity to factors as illumination, control pattern aspect and image resolution on accuracy of the matching algorithm should also be investigated. (ii) The disparity field obtained with this simple matching procedure presented here is quantized because of the discrete aspect of digital images. However, the true disparity to be estimated is, for most points, small and not integer. In this paper sub-pixel accuracy was achieved by interpolating signal in order to have an image value also for non-integer pixel positions. It would be interesting, however, to compare relative

ARTICLE IN PRESS 650

K. Genovese, C. Pappalettere / Optics and Lasers in Engineering 45 (2007) 637–650

merits of alternative approaches proposed in literature to achieve sub-pixel accuracy. (iii) Area-based algorithm works fast and well if patches to be matched does not contain too high frequencies and if geometrical and radiometric distortions are minimized. Both conditions are usually not satisfied in this case since using a cone mirror generates rather smoothly textured areas and image pairs with significant geometric differences. However, a performance comparison between conventional area-based correlation algorithms and the proposed feature-based approach could be useful to validate the qualitative statements made in this work. Future research will focus on the application of the proposed method to measuring deformation on hyperelastic tubular samples and to developing more reliable computer algorithms to improve matching efficiency. References

[7]

[8] [9]

[10]

[11] [12]

[13]

[14] [15]

[1] Hayashi K. Experimental approaches on measuring the mechanical properties and constitutive laws of arterial walls. ASME Journal of Biomechanical Engineering 1993;115:481–8. [2] Everett WN, Shih P, Humphrey JD. A bi-plane video-based system for studying the mechanics of arterial bifurcations. Experimental Mechanics 2005;45(4):377–82. [3] Optical Engineering. Special issue on optical methods for shape measurement. SPIE 2000;39(1). [4] Genovese K, Pappalettere C. Feasibility of 3D reconstruction of vascular segments via Projection Moire´. SPIE Proceedings, vol. 5776, 8th international symposium on laser metrology, 2005, Merida Yucatan, Mexico. [5] Matthys DR, Gilbert JA, Greguss P. Endoscopic measurement using radial metrology with digital correlation. Optical Engineering 1991;30(19):1140–455. [6] Albertazzi GA, Melao I, Devece E. Measurement of thermal deformation of an engine piston using a conical mirror and ESPI.

[16]

[17] [18] [19]

[20]

In: Pryputniewicz RJ, Brown GM, Jueptner WP, editors. Proceedings SPIE, Laser interferometry IX: applications. Vol. 3479, 1998, p. 274–283. Lin SS, Bajcsy R. High resolution catadioptric omni-directional stereo sensor for robot-vision, Proceedings of the 2003 IEEE international conference on robotics & automation, 2003. Barnard ST, Fishler MA. Computational stereo. ACM Computing Surveys 1982;14:553–72. Alvertos N, Brzakovic D, Gonzales RC. Camera geometries for image matching in 3D machine vision. IEEE Transactions on Pattern Analysis and Machine Intelligence 1989;11(9):897–915. Cardenas-Garcia JF, Yao HG, Zheng S. 3D Reconstruction of objects using stereo imaging. Optics and Lasers in Engineering 1995;22:193–213. Dhond UR, Aggarwal JK. Structure from stereo-a review. IEEE Transactions on Systems, Man and Cybernetics 1989;19:1489–510. Barnea DI, Silverman HF. A class of algorithms for fast digital image registration. IEEE Transactions on Computers 1972;C-21(82): 179–86. Helm JD, McNeill SR, Sutton MA. Improved 3-D Image Correlation for Surface Displacement Measurement. Optical Engineering 1996;35(7):1911–20. Grimson WEL. Computational experiments with a feature based stereo algorithm. IEEE Trans. PAMI 1985;PAMI-7(1):17–34. Synnergren P, Sjodahl M. A stereoscopic digital speckle photography system for 3-D displacement field measurements. Optics and Lasers in Engineering 1999;31:425–43. Brown MZ, Burschka D, Hager GD. Advances in computational stereo. IEEE Transactions on Pattern Analysis and Machine Intelligence 2003;25(8):993–1008. Faugeras OD. Three-dimensional computer vision: a geometric viewpoint. Cambridge, MA: MIT Press; 1993. Rao SS. Engineering optimization. New York: Wiley Interscience; 1997. Sutton MA, Cheng MQ, Peters WH, Chao YJ, McNeill SR. Application of an optimized digital correlation method to planar deformation analysis. Image and Vision Computing 1986;4(3): 143–51. Sutton MA, McNeill SR, Jang J, Babai MK. Effects of subpixel image restoration on digital correlation error estimates. Optical Engineering 1988;27(10):870–7.