Automatic inspection of mechanical parts using geometric models and laser range finder data

Autom.atic inspection of mechanical parts using geometric models and laker range finder data A D Marshall, R R Martin and D Hutber*

Current industrial 3D inspection machines involve tactile sensing, using co-ordinate measuring machines (CMM). These CMMs are slow in operation because a single probe has to take many readings. They also require detailed individual programming to conduct each inspection task of a different object. Programes are very time consuming to produce and are prone to errors. Thus if small quantities of an object need inspecting this approach is not very cost effective. The alternative of manual inspection is also time-consuming and inefficient. Hence the need for an automated approach capable of inspecting a wide range of objects is real and necessary. Most CMMs in current use in industry are programmed as fotlows: by knowing the object’s initial orientation in the CMM the location of certain reference points or faces can be measured. By then specifying where to find other features and measuring their positions accurately, the correctness of the object can be verified. It should be noted that all the geometric information comes from engineering drawings and this information is explicitly given to the CMM by a programmer. At no stage does any geometric reasoning occur within the CMM. Keywords: automatic inspection, 3D vision, laser/range Thus even though tactile inspection provides accufinder data rate analysis of an object, for reasons of speed, flexibility and that tactile sensing is not always desirable, the CMM approach is not suited to a totally automated method of inspection. A logical alternative There exists a great demand in manufacturing industry approach is to use an accurate 3D vision system to for a general purpose three dimensional automatic obtain views of an object and perform inspection’. inspection system. A wide variety of three dimensiona Most general purpose vision-based inspection machobjects ranging from VLSI chip packages to Iarge metal ines to-date have been built to inspect two dimensional castings need inspecting, often in small volumes. objects, such as PCB boards and sheet metal components, but none have been built for general purpose Department of Computing Mathematics, University of Wales, PO inspection of three dimensional objects. Inspection in Box 916, Cardiff CF2 4YW, UK 2D is made a lot easier by the fact that information *British Aerospace, Sowerby Research Centre, PO Box 5, Filton, provided by a single camera is itself two dimensional in Bristol BS12 7QW, UK nature, so there exists a direct one-to-one correspondPaper receivad: IO August 1990; revised paper received: 12 October ence between features in the model and the scene. 1990 Ltd 0262-88561911006385-21 @ 1991 Butte~orth-Heinemann 385 voI9 no 6 december 1991 This paper addresses the problems associated in processing arrays of depth data in order to achieve the goal of automatic inspection of mechanical parts, i.e. developing general model-based inspection strategies that can be applied to a range of objects. The main problems in processing this data are segmenting out reliable primitives from the data and matching these primitives to those in a stored geometric model of an object. The ability of a 30 vision system to provide depth data accurate enough to perform automatic inspection tasks was until recently only possible at a short range from an object, typically a few centimetres. However it is now possible to produce dense data from a vision system situated further from the object, typically half a metre to a metre. Such a system is outlined. Some current modelbased matching techniques are assessed for their suitability for employment in inspection type tasks. One approach is adopted and modifications that improve the efficiency and accuracy of the method for inspection purposes are presented. Finally, an inspection strategy is outlined and its performance assessed. Results are presented on both art~~cial and real depth data.

However, in 3D vision there exists a many-to-one correspondence between pairs of features which implies that extra information about these features is required. Two cameras mya be used to obtain a three dimensional description of a scene by combining the two cameras’ different 2D views. Difficulties arise in correlating the two views with each other which makes the matching of features more complex. Also, due to camera inaccuracies, errors in the calculation of 3D coordinates when the two views are combined, as shown below can be large enough to prevent accurate inspection from being possible. Until recently this was a problem, but with developments in camera technology it has now become possible to produce data accurate enough to enable inspection from a realistic range. This paper examines the practical ability of a vision system to perform automatic inspection tasks. First, an approach using a 3D vision system is outlined; such a system is briefly described, and errors in acquiring the depth data in the system analysed. Second, a comparison of existing matching and positioning strategies is given, and their performance for inspection is assessed. The adopted strategy is then detailed, and some improvements presented. Finally, the strategy’s ability to enable inspection to be performed is assessed. It should be noted that although the general problem studied is that of automatic inspection, the methods described may have useful results that contribute to many other related fields of study. These fields include industrial assembly, general computer integrated manufacturing (CIM), object recognition, and even computer navigation and machine perception.

APPROACH SYSTEM

- OUTLINE

OF THE VISION

Many 3D visiorrrstems have been developed over the last few years . Some produce a relatively small number of depth readings over a large region, for a given area of interest, some produce only details of edges present in a scene, while others produce a complete scan of the area of the scene producing an array of depth readings. Such an array of depth readings is called a depth map. If the array contains a large number of readings for a relatively small area of interest the depth map is said to be dense. A typical dense depth map may have individual depth readings (referred to as pixels) at every millimetre for a 300 x 300 mm2 area. Inspection tasks not only need to be able to perform accurate local measuring tasks (e.g. diameter of holes) but also infer that features of the inspected object are positioned and oriented accurately relative to each other globally (e.g. position wrt some datum face). Clearly, depth maps are capable of providing such a complete picture of a scene with both local and global knowledge of features. Information about features is more likely to be accurate given a large number of measured depth readings per unit area, because least squares fitting techniques are typically employed to derive this information from the data (see below). Thus the use of dense depth maps seems to be more desirable than the other possibilities mentioned above’. The Sowerby Research Centre (SRC) of British Aerospace (BAe) has recently developed a vision 386

View of Object L

Vision System

Errors due to camera inaccuracy L

Segmentation

Errors due to \ Y region fitting

Errors due to numerical methods A,

Matching

Inspection

Figure 1. Components of the system and data requirements.

system capable of providing accurate dense depth data for this purpose. Clearly, the inspection of every possible three dimensional object is impossible, and so a sub-task has been chosen, requiring that the part must lie within a cube of side OSm. The further constraint that the part must consist of planar, cylindrical and/or spherical faces has also been adopted. This is realistic since approximately 80-90% of manufactured parts that BAe wishes to consider for inspection are estimated to consist solely of such features’. Approach outline The outline of the elements comprising of our system for inspection and their respective data are summarized in Figure 1. The different stages are briefly discussed below. Vision acquisition

system

The vision system consists of a matched pair of cameras and a laser scanner mounted on an optical bench to reduce vibration (see Figure 2). The cameras are a high sensitivity CCD type, each mounted on turntables capable of independent rotation. The cameras are

Figure 2. SRC vision system image and vision computing

interfaced to a Sun workstation via a frame store which is able to grab frames from both cameras. The turntables are also controlled by the workstation. The scanner consists of two mirrors that can be moved independently to steer a spot of laser light. The mirrors are manipulated in such a way as to project a vertical stripe of laser light onto the object and move the stripe across it. The whole system is controlled by the workstation which also has access to a database of solid models. Geometric models of boundary representation type are used8-14. The explicit representation of surface types and connectivity is of direct usefulness when matching the model to the observed data. Some relevant aspects to the operation of the vision system are now described further.

Projective

“‘“%-‘y

Image plane

I

World cmrdinate ?.yslem

Figure 3. Camera calibration model

Camera calibration and depth extraction Initially the cameras of the system have to be calibrated. Camera calibration is probably the most important phase in the depth data acquisition stage, since without accurate calbration there is no hope of performing depth readings to a high degree of accuracy. Unfortunately, the camera calibration problem is also one of the most difficult problems in stereo vision. This is due to the fact that it is difficult to find the stereo correspondence between a point in the image of the right camera and to that of the left camera. The calibration technique employed by the described vision system is based on the techniques proposed by Tsai’57’6. Camera calibration is the process of determining the internal geometric and optical characteristics of the cameras (internal parameters) and the threedimensional position and orientation of the cameras relative to some world coordinate system (external parameters).

Tsai’s technique was designed to be accurate, efficient, versatile and to work with ‘off the shelf’ cameras and lenses. Most other existing calibration techniques do not satisfy these criteria, requiring either a special camera and/or processing equipment - prohibiting large scale automation and versatility (it may also be uneconomical), non-linear optimization - the determination of the internal and external camera parameters can involve non-linear estimations requiring iterative searches, which are not efficient and normally require good initial guesses to achieve reliable estimations. It is better if the parameters can be expressed linearly; or ignoring/assuming important camera properties - some linear estimations of the parameters require that lense distortions and uncertainties in image scale factors are ignored. Others assume certain camera information (e.g. focal length). This leads to inaccurate results. The intrinsic parameters that need calibrating are: 1 The effective focal length f or the image plane to projective centre distance 2 The lens distortion parameters kI and k2. These parameters model the radial lens distortion in the x and y directions. Only radial distortions are considered. The other kind of lens distortion, tangential distortion, is not included in the model, since any gain in accuracy is negligible compared to sensor ~019 no 6 december 1991

quantization, and also it would cause numerical instability in method’5T’6. The uncertainty scale factor s, between computed coordinates for pixels and actual pixel (CCD element) coordinates. That is to say, errors occur in the digitization. This arises due to slight hardware timing mismatches between image acquisition hardware and camera scanning hardware, or imprecisions in the hardware itself’“. The computer image coordinates of the origin in the image plane (c,, cY). The extrinsic parameters l

that need calibrating are:

The 3 X 3 rotation matrix R and translation vector t that define the transformation between world and camera coordinates.

The relationship between some of these parameters is illustrated in the camera model of Figure 3. From this model a system of equations can be set upI that describe the relationships between these parameters. The calibration of each camera is performed as follows: A calibration chart is placed in the view of the cameras. The calibration chart contains a lattice of identical black squares of known dimension, the corners of which are taken as calibration points. Since all points on the calibration chart are coplanar, the world coordinate system can be chosen so that the plane of the calibration chart is such that z, = 0. In order to maximize the efficiency of the method, best results are obtained when the world coordinate system is also chosen such that the origin is not close to the centre of view of the y axis of the camera coordinate system l6 . This ensures that not all the parameters of the translation vector t are not all zero simultaneously, and thus the computation procedure can be made more unified and simpler. Having obtained the mage of the calibration chart, the locations of the calibration points need detecting. The dark areas in the image are each iteratively grown in all directions to produce regions of connected dark pixels. The resulting dark regions correspond to the black squares on the chart. The calibration points (corners of squares) are then produced by tracing around the edges of the dark regions. 387

Knowing the position of each calibration point in the computer frame store and camera and frame specifications, then parameter values for every real image value can be approximated by from the derived system of equations (see Tsai’6). The basic method for the calibration of camera parameters can be extended to use multiple views of the calibration chart by moving the camera about on the turntables. Using multiple views over-determines the equations derived above, and least squares methods can be used to solve them giving more accurate results. The successive rotation matrices can be recalculated knowing each rotation of the turntables enabling the initial rotation matrix R to be estimated. Having calibrated the vision system, depth maps can then be extracted from the scene. The system operates by projecting and moving the stripe across the scene to obtain a series of vertical stripes of pixels which are processed to give the required dense depth map. The depth of a point is measured as being the distance from one of the cameras, adopted as the master camera. Knowing the (calibrated) relevant geometry and optical properties of the cameras the depth map is constructed as follows: 1 For each vertical stripe of laser light, form an image of the stripe in the pair of frames associated with each camera. For each row in the image from the master camera, search until the stripe is found at point P(ij), say. Form a line 1 passing through the centre C, of the muster camera and P(i,j) (see Figure 4). Construct the epipofar line which is the projection of the line 1 onto the image from the other camera (see Figure 4). This is achieved by projecting two arbitrary points Pi and Pz onto the image, and constructing a line between the two projections. Search the epipolar line for the laser stripe. If the laser stripe is found, the (x, y, Z) coordinates are calculated from the computed position P, of the point on line 1(see Figure 4), and written in the depth map. Otherwise no value is written.

First, the effects of spatial quantisation need to be evaluated. This arises due to the optical system of the cameras giving the position of a point of light correct to, typically, 0.2 pixel. This accuracy is based on observed differences in obtaining laser stripe images under identical images. It clearly depends on surface reflectance properties of objects, the lenses and camera system, and the framestore. It also depends on the working range of the vision system, as illustrated below. It should also be noted the laser spot size is finite, and this implies that the depth reading is spatially averaged over the spot. Thus an ‘area of uncertainty’ in the precise location of the sensed data point arises, as illustrated in Figure 5. Both the angular and depth resolution can be readily evaluated. Consider the schematic view of the system (see Figure 5). Here: y = field of view of cameras (radians). n = number of resolvable points across field of view. d = separation of cameras (m). B = angle between optical axis and image point (radians). z = depth (m). 6 = angular resolution of the camera. 8, = depth resolution. Then the angular resolution may be expressed by: S = 0.5 tan-’ ((tan -y)tn)

(I)

and the depth resolution by: Sz = 2d tan (6) (z*/&) (tan* l3+ (d/s) tan I3+ 1)

(2)

Thus, the error in the depth reading is proportional the square of the depth value.

to

Errors in a stereo system In order to assess the performance of subsequent matching, segmentation and inspection algorithms, the typical error in acquiring a depth value using the above arrangement needs to be evaluated:

d master

camera

Figure 4. Measuring a depth value 388

Figure 5. Uncertainty in stereo calculation of vision system image and vision computing

Inserting typical values, for this system, in equations (I) and (2):

acy of 0.01 inch or 0.254mm. Thus, to determine the depth of a point to this accuracy, allowing for depth resolution errors, the depth reading must be determined to within O-254-2 x 0.052 = 0.15 mm under the above assumptions. The system’s apparatus can achieve this, and is intended to do so under normal working conditions. However, to assess the respective performance of various algorithms their behaviour for other accuracies also needs examination (situations such as this may arise due to better/worse equipment or more/

L3= 0.0 radians, d=0.3m, y = 0.08 radians, n=512x5=2560, 2 = 0.7 m gives 6 = 4.9 X 1OW‘5 radians and 8, = 0.052mm. Typical inspection measurements require an accur-

...-.... ..._._ ._._ .....__. MI

. . . . . .

-.“..*~.___,

“...._._ _......”

e*. ;.-

b

d

. _-

less rigorous inspection requirements). Thus artificial depth maps are created for the purpose of this study. Under normal working conditions the above accuracies are taken as the system’s standard accuracy. Depth maps are simulated as an error in depth reading by adding Gaussian noise. The standard deviation of the Gaussian distribution is taken as to the size of the error. Thus, for standard accuracy 1 sd corresponds to 0.15 mm. Differing range can also be modelled by noting that the depth error is proportional to (depth)’ so increasing the range by a factor of %’ 2 corresponds to doubling the standard deviation of the added noise. Some examples of depth maps, with and without noise, are shown in Figure 6. Further details of the system and its calibration may be found elsewhere’4-‘8.

Segmentation, matching and inspection Having acquired a depth map, how do we best proceed to interrogate this huge amount of information in order to achieve an accurate match to a geometric model, stored in the computer? The term segmentation is used in the following discussion to mean processing the ‘raw’ depth data to yield some higher order description. Processing the depth data occurs in two stages: 1 Some means of segmenting the data in order to derive some higher level primitives than depth points to provide a scene description. 2 Adopting an appropriate matching strategy to match this scene description to a geometric model description, and to give the location and orientation of the scene relative to the model. Naturally the method adopted in (2) depends on the data acquired from (1). Typically a depth map is segmented to yield one or more of: 1 Oriented surface points (points in space with an associated surface normal). 2 Edges of the object. 3 Surface primitive descriptions such as planes, cylinders, spheres, quadrics, etc. Having obtained a description of the location and orientation of the viewed object relative to the model, the actual positioning of various features can then be verified to a high precision. This is the inspection part of the system. For reasons described later, inspection is best carried out on multiple views of an object. However, to justify the feasibility of this overall approach to inspection, only a single scene need be considered. Therefore, for the purposes of this paper inspection is limited to the inspection of a single view of an object. Many tests have been performed on artificially produced depth maps with added Gaussian noise, so that the behaviour of the method could be predicted independently of the actual hardware used. COWIPARISON ALGORITHMS

OF

EXISTING

Grimson Lozano-Perez method2’-25 This method was originally designed to match oriented surface points21-23, but has also been extended to match edges 6324-27 . Implementations of both methods exist. We describe the method in terms of matching edges, but the discussion is equally valid for the case of oriented surface points. At each level of the tree, one (or more) of the edges from the scene is matched with each of the N possible edges in the model. Each node has N children representing taking the match proposed so far together with all possible matches for the current scene edge. The result is a partially complete tree, where complete branches correspond to complete matchings of object and scene edges. In the Figure 7 example, the model has N edges, and three edges have been found in the scene. Thus the pairing (2,1,1) corresponds to: edge 1 in the scene matching edge 2 in the model, edge 2 in the scene matching edge 1 in the model, and edge 3 in the scene matching edge 1 in the model. (Due to errors in segmentation and other causes, e.g. partial occlusion, it is possible for more than one edge in the scene to correctly match a single edge in the model.) The interpretation tree is searched in depth first manner, and pruned by rejecting edges that fail a potential match test. The test is achieved by employing local geometric constraints on the distance, angle and direction between matching edges. If the test is unsuccessful the node and all subsequent child nodes may be removed as candidates for a match, as none of these ‘children’ can possibly correct the current mismatch. In this way the potentially large search space is controlled. The second stage involves interrogating the sets of possible interpretations which are left after all scene Faces

D

MATCHING

Over the last few years, many model-based matching strategies have been developed. These have been 390

primarily concerned with object recognition purposes, although most produce information about the position and orientation of object. Not many of them have been applied to more detailed applications, such as inspection A few methods are briefly described here, and their suitability for inspection to take place is discussed. Those chosen reflect basic trends in matching strategies employing different primitive segmentations. A more detailed discussion of current strategies may be found in Best and Jain” and Chin and Dye?‘.

PAihNG

( 241,I )

Figure 7. Example interpretation tree image and vision computing

primitives have been considered, to estimate the transformations (rotation, scale and translation) between scene and model coordinates. This is achieved by approximating the transformation by fitting the scene’s description to the matched model edges. (Details of the mathematics involved can be found in Grimson and Lozano-Perez2i.) Alternatively, an approach similar to the estimation of rotation and translation presented in later sections may be employed.

Faugeras Hebert method28-34 This method matches primitive surfaces extracted from a scene to ones in the model. In the original approach the method is described using either planar or quadric surfaces. The basic approach is to try to match a primitive surface from the model, say mi, to one segmented from the scene si to obtain a match pair (Si, mi). This can be achieved using a tree structure similar to the interpretation tree above. However, the estimation of the transformation (rotation and translation) describing the mapping from model to scene is used as a constraint to control the search as it proceeds, in contrast to the previous method where the transformation is only estimated after matching is completed. From a small set of match pairs, a match list, an estimate of the transformation can be made. Pairs of primitives from the scene and model are then added to the list if they are compatible with the match so far. Each time a successful match is found a new estimate of the transformation is calculated. The process is continued until all the scene primitives have been matched. If at any stage no possible match pair can be found, then a backtracking mechanism is used. This is easily visualized by regarding the search as a tree search. A match is successful if: 1 It is consistent with the existing transformation. 2 The error of the updated transformation is within the allowed tolerance. Thus the method exploits the global constraint of rigidity under the transformation to control the potentially combinatorially explosive search space. At any level the estimation of the transformation is found by solving a least squares minimization problem. Also, a cheap local constraint is employed at each level in the search in order to eliminate grossly mismatching pairs more quickly. This local test can be either or both of: The angle between a matched primitive in the scene and the current scene primitive is equal to the angle between a corresponding matched primitive in the model and the potentially matching model primitive, to within a given tolerance. The angle between the current scene primitive and the current model primitive rotated by the current estimation of the rotation is equal to zero, to within a given tolerance. This method will be examined in greater detail later, as it forms the basis of the inspection system proposed. ~019 no 6 december

I991

Oshima Shirai method35 This method operates in two stages: First, a learning phase takes place in which known objects are individually shown to the system in different poses. A description of the object is built up, in terms of properties of surfaces and their relationships, to create a model of the object. The description is dependent on the segmentation used, but typical properties include surface type, number of adjacent regions, area, perimeter, and many other statistical measures. Second, the matching phase operates by matching objects in a scene to a database of learned object descriptions. The best two surfaces are chosen from the scene and are used as a kernel for the match. Having matched two model surfaces to these kernel surfaces, further candidate surfaces from the model are hypothesized and tested for in the scene. This process continues until ail scene surfaces are matched.

Relaxation/labelling techniques3G38 In these approaches matching is posed as a labelling problem, thus a model primitive Mi is labelled with (i.e. matched to) a scene primitive &. Region-based measurements are normally used (corresponding to numeric and topological features similar to those employed by the Oshima Shirai method). Each primitive is given a quality measurement (normally a probability) for its possible labelling to model primitive. Starting from these initial measurements, the goal of the relaxation technique is to reduce iteratively the ambiguity and incoherence of the initial labels by employing a coherence measure of a set of matches. This may be achieved by ~onside~ng an iteratively increasing neighbourhood of each primitive, and the compatibility of the primitive with its neighbourhood.

Hough transform techniques3g42 The basic idea of these techniques is to quantize the space of the relevant transformations (Hough space) relating primitives (points, lines or surfaces) in the model to the scene. Each match between primitives from the scene and model corresponds to a family of viable transformations. An array of accumulator cells is set up, each cell corresponding to a possible transformation, and for each possible transformation the corresponding accumulator cell is incremented by one. The best matches occur at the cells containing the highest count.

PRACTICAL COMPARISON OF EXISTING MATCHING TECHNIQUES WITH RESPECT TO INSPECTION The methods described above are now to be assessed for their ability to provide sufficiently accurate results to achieve the goal of inspection.

Preliminary discussion It is apparent that the last three methods described are

391

not suited to inspection serious drawbacks:

tasks,

as they all suffer

from

The Oshima Shirai approach is data driven, and therefore requires a stable segmentation method. For example, a surface viewed from different views must not have a different description (e.g. number of adjacent regions) to provide reliable matches. At present there is no guarantee of a stable segmentation. This means that an arbitrary view of an object that may have been incorrectly segmented (or even omitted) during the learning process may lead to difficulty in matching, because some surfaces or their properties may have been missed. Also, some surface properties do not hold when occluded views of surfaces occur”. The relaxation techniques are very sensitive to initial matches, since local rather than global interaction of features is considered. Thus convergence towards a best match is not always guaranteed3*, and the method may stop at a local rather than global minimum. The Hough space is six dimensional (six degrees of freedom in the transformation - three rotational and three translational), and the accuracy of the Hough transform can only be guaranteed if the Hough space is divided up into small partitions (requiring a large 6 dimensional array). This is computationally very expensive (both in time and memory) for the level of accuracy required for inspection. Since this method is completely global, some information gets misplaced, which makes it difficult for the method to converge on the correct solution. Also, the Hough transform does not have an efficient and simple method of employing constraints to reject grossly mismatching primitives3*. For the above reasons, only the Grimson Lozano-Perez and Faugeras Hebert methods were implemented and examined further.

Segmentation Before examining these methods in practice, a brief discussion on the segmentation methods employed to produce the respective primitives is merited. As previously stated, the Grimson Lozano-Perez method can be employed to match oriented model points or edges, while the Faugeras Hebert method can use planes and We have not implemented the quadric surfaces. Faugeras Hebert method to use quadric surfaces as reliable segmentations are difficult to achieve3”. The various approaches used to segment these primitives are described below: Oriented model points - these are data points on an object in the scene together with an associated surface normal, assumed to point away from the object. This is easily estimated from the depth map by considering a small region of pixels around the point, typically 3 x 3, and approximating a plane to this set of pixels. The plane normal is then taken to be the surface normal at the oint3y74w5. Edges - an edge is taken to occur when a marked change in gradient occurs across a set of points. Thus 392

an edge detector that uses gradient ators has been imp1emented40*46.

different

oper-

Planes - a region growing method is employed which, after an initial grouping of perfectly fitting points to (small) planes, merges planes that satisfy similarity constraints. This is described further in Appendix A.

Analysis of the Faugeras Lozano-Perez methods

Hebert

and Grimson

Implementations of both these methods were examined using artificial depth maps with added Gaussian noise of different amounts, modelled on the standard accuracy given above. Tests to find which of the two methods gave the most accurate transformation from the same view (with the same noise) were carried out on relatively simple planar objects, some of whose depth maps are illustrated in Figure 6. All images are 256 x 256 depth maps and are of a cube - side 200mm (see Figure 6a), a widget - max. sides 150 x 90 x 60mm (see Figure 6b), and a widget with a cylindrical hole - radius 15mm drilled perpendicularly through a face (see Figure 6~). The Faugeras Hebert method always yielded the most accurate results - to within 0.05mm and 0.01” at standard system operation (depth accuracy of 0.152 mm), and to within 0.03” and 0.2 mm translation under moderate noise (depth accuracies of about 0.5 mm-O.6mm). The method was more stable to noise, and was also quicker in reaching the given results, typically taking a few seconds on a Sun 3/50. The Grimson Lozano-Perez method using oriented surface points was too combinatorially explosive to be satisfactory. Even with a view of three faces of a cube (Figure 6a) and using a sample of nine evenly distributed points per face, the time to achieve a match was more than 20 minutes. Using fewer points gave inaccurate estimates of the transformation. These observations can easily be explained in that more primitives are being matched in this method. However, the Grimson Lozano-Perez method using edges as primitives gave quite satisfactory results under low noise - to within 1” and 1 mm translation at standard system accuracy, but became increasingly worse with increasing noise to within 3-4”, l-2mm translation, with depth errors of 0.3mm to 0.6mm. This was due to increasingly poor edge estimates from the segmentations. The results for the matching of a more complicated test object, such as the widget (Figure 6b), using both approaches, are summarized in Figure 8. The fact that better results for the transformation are obtained from the Faugeras Hebert method is intuitively reasonable, since here all of the data for a given piece of surface are used, not just those at edges at vertices, which are furthermore highly prone to noise. These results also seem to agree with the various critiques recent algorithms in 1 From the results obtained it was concluded that planar primitives are much better than either edges or points to drive the matching, and that the Faugeras Hebert method is satisfactory - it is robust in the image and vision computing

0.3 -

/

06 0

I’

d’ 0.1

I

a w

0.2

0.3

Error in depth reading (mm)

Figure 8. Comparison of the Crimson and Faugeras Hebert methods

Lozano-Perez

presence of noise well beyond the level of accuracy required for inspection to be possible. Furthermore, it is easy to use non-planar primitives in the Faugeras Hebert method. This is not the case with the Grimson Lozano-Perez method and other methods, which were designed only for planar objects. As well as accuracy issues, another problem that occurs when using laser range-finding techniques is that the depth information obtained from edges is less accurate than face data. This is because scattering of the laser light occurs at these edges’,4,‘y. Thus it can be shown that real edge data will provided less accurate results than the simulated data presented here. For these reasons, the Faugeras Hebert approach2ss4 was adopted as the best strategy to adopt for inspection.

FAUGERAS HEBERT METHOD IMPROVEMENTS

AND

The basic approach of the method28-34 was discussed above. More details of the method are presented, and then improvements to the efficiency of the method are described. The extension of the method to non-planar surfaces that provide potentially reliable matching is also discussed. Finally, results, including the effects of the improvements and segmentation method, are presented.

More detailed description of the Faugeras Hebert method As mentioned before, the method may be used to match planar and quadric surfaces. The case for matching planar primitives is discussed below. Details for quadric surfaces may be found in Faugeras et

a/.2Y%_ 2

The approach is to obtain a list of matching pairs of planes-one from the scene, the other from the model that contains all of the planes extracted from the scene. Before each match pair is added to the list it is first ~019 no 6 december 1991

tested for correctness under the current estimated transformation of the match. If this match is successful then the transformation is updated using the new match pair. Two non-parallel match pairs are initially required to estimate the rotation, and three non-parallel match pairs are required for the translation estimation. In practice, in the original method, only the rotation component of the translation is used to control most of the search. The translation part is only employed in the lower levels of the search to find the best solution amongst several potential matchings. The reason for this is that the rotation provides the more powerful constraint, and employing only this constraint saves time. The estimation of both the translation and rotation represent the majority of the computation time in the algorithm, and in order to increase efficiency these calculations should be kept to a minimum. Grossly mismatching planes can be quickly eliminated by using local constraints as described above. Most other mismatching planes can be elimated by simply employing the rotation part of the transformation, with only certain symmetries in objects not being eliminated by this ‘rotational inconsistency’. The estimation of this transformation is now considered.

Estimation

of the transformation

At any level the estimation is treated as a least squares minimization problem. Let the transformation T be expressed as t*R, where t = translation and R = rotation. This implies that the rotation is applied first, the translation second. A plane Pi is represented by two parameters Vi and d,, where Vi is a unit vector normal to P, and d; is the signed distance of P, to the origin. If the transformation T is applied to a plane P, (vi, dJ, a plane Pli(vli,dl,) is obtained where: ~1, =

Rvi

and

dl; = vl;t + di. Consider the scene to be segmented into N planes and the model to contain N’ planes, where N’ 2 N. Let {Pi (vi, di)},=l..N be the set of scene planes and {P~(v’i,d’i)},=l..N. be the set of model planes. Consider a matching M, where M = { P,(v,,d,), P’; of N pairs. Then (v’lyd’i))i=l..N is a matching values of t and R are required that minimize:

where W is a weight factor, and W’i = Rv’;. This can be decomposed into two terms: the rotation part - a sum over IIvi- W’iJ12- and the translation part - a sum over Id,-d’;w’;tl*. Hence the rotation part is minimized separately from the translation. The weight factor is chosen so that the errors given by both minimizations are of the same order. The value was chosen based on preliminary results with simple objects. The estimation of the rotation part is not a straightforward minimization, and quarternions may be used to reduce the problem to a ‘classical’ least squares problem2”. The translation part is solved using the pseudoinverse method2’. 393

matrixes U, V and Z would introduce similar round off errors to those discussed previously. Thus A is first reduced to bi-diagonal form using Householder transformations, and secondly by the C.2: algorithm is used find the singular values of

Matrices A and Z are formed, such that: A is the NX 3 matrix : [w’, .,__,,w’,JT, where

w’; = Rv’; and Z is the N X 1 matrix : [di - d’[ ,,_.__,dN-

4vjT,

Then the error in t is minimized when t has the value : tmin = (ATA)-’ ATZ with corresponding error : emi” = ZT(Z - At,i”).

(2) Let fn;U

g+nd y, both N vectors, be defined by: Z

y=S’.t. Then, we may write equation (4) as:

In practice, instead of computing the inverse of ATA, the linear system of equations formed by (ATA).t,i” = AT.Z is solved.

min I

/[At-Z

t

Improvement

= min I/ CV’t-g

Although the Faugeras Hebert method basically provide satisfactory results, various modifications have been made to improve the method, with respect to the inspection application, etc., so that it is more accurate and efficient. These are described below.

Improvements in accuracy The improvements described here are concerned with the estimate of the translation part of the matching process. Replacing the pseudo-inverse

method

The Pseudo-inverse method used for minimising the translation part of equation (3): c/d,--d:-

11 = min 11UCVTt - Ug I]

w;tl 2

was found not to be very accurate and unstable in the presence noise. The method is also unnecessarily prone to rounding errors. This is a well known problem in numerical analysis, and arises since the conditioning of ATA is worse than that of A4’. This frequently leads to difficulties which cannot be remedied once ATA has been formed. Detailed descriptions and error analysis of the pseudo-inverse method are well documented-%-so. This problem is illustrated in Appendix B. Various methods exist for solving this type of minimization problem, and reducing these round off errors4s52. The method preferred for minimizing the expression in equation (4) is described below; in order to obtain a solution N must be at least 3, so in the subsequent calculations N is assumed to be at least 3. This problem may be solved as follows:

= min 11Cy-g Y

/I

(5)

11

(3) Let y, be the unique solution of (5) then: ZY, = g. (4) Then ail solutions of (4) are of the form: t=vyw. This method has two advantages over the alternative methods4s5’. First, the only iterative part of this method is in the calculation of the SVD; both the other methods apply iteration over the whole solution. In particular, Golub4s50 uses iteration to converge to the solution, and at each stage in the iteration calculates the SVD by iteration. The adopted method therefore finds the solution more quickly, and with no loss in accuracy. The other advantage of this approach is that it works for cases when A is rank deficient5’. This case will arise when parallel primitives occur in a scene. Using the centres of gravity of regions

If the rotation estimate inaccurately estimates the position of the plane then the estimate of the translation will be in error as well. Figure 9 illustrates this fact. In this figure the rotation estimate has resulted in the positioning of the model plane normaf vri, due to tis pairing with the scene plane normal Vi being in error. The effects of this (rotation) error on the translation, sccnc

ROtiMd

Plane

PlLUl‘Z

Factorize A using Singular Value Decomposition (SVD)” so that:

k = tiCVT where: UTU = VTV = VVT = 1, the 3 x 3 identity matrix and: C = diagonal matrix (o,, a2, g3) In fact, the N x 3 matrix U consists of the orthonormalized eigenvectors of AA*, and the 3 x 3 matrix V consists of the o~honormalized eigenvectors of ATA. The diagonal elements of Z are the nonnegative square roots of the eigenvalues of ATA (the singular values). Obviously, any straightforward evaluation of the 394

Translated

1 d -dcos13,

I I I

_c

Model

t

I I

dl

d cos Ol

I I

I

I

I

I I

I I

t

Figure 9. Errors in translation estimate image and vision computing

relative to the camera coordinates (shown in diagram), are now considered. Let 8 be the error in the matching of scene and model plane normals. Then the error in the z axis placement of the model plane is dtan 8 and the error in the x axis is d(l-cos 0). Given that d can be of the order of 500-700mm large, values of 8 are not required for the error to become significant. The least squares approximation of a plane to a set of points always passes through the centre of gravity of set of points. Thus if the model plane is transformed such that it passes through the centre of gravity of the scene plane it is being matched to, then the errors in the translation estimate are kept to a minimum. This means that minimization for the translation now minimizes the perpendicular distance between the centre of gravity of the scene plane and the model plane in the direction of the model plane, i.e.: min C 1w:Ci-d{-wit I2 (6) I where Cj = the centre of gravity of the scene plane. Equation (6) is then minimized using singular decomposition, as described above. The centre of gravity is calculated during the segmentation. Improvements

in the efficiency of the method

Some improvements on the efficiency These are discussed separately.

were made.

Order of matching pairs

The order in which feasible matching pairs are taken to estimate the transformation was found to have an effect on the speed and eventual accuracy of the result. Consider the case where three primitives from the scene sl, s2 and s3, say, are matched to three primitives from the model ml, m2, m3 to give a set of three match pairs ((~1, ml), (~2, m2h (~3, m3)>. If the first match {(s,, ml), (sz, m2)} is poor then the estimate of the rotation (the only component of the transformation evaluated at this stage) from this pair will also be poor. Since the estimate of the transformation for the next possible pair depends on the previous estimate, no matter how good a match is found the new estimate will also be poor. Thus such errors are propagated through the whole match. Again further cl;t$ of the process are described in Faugeras et . . It has been found that if the best possible matching

pairs are used first in the valuation, then the subsequent final estimate of the transformation is more accurate and the search space can be more efficiently constrained. What is needed is a quick, reliable way of finding the best possible scene primitives for matching. By taking a measure for each surface of its surface area, the quality of fit of the data points to the surface, the type of primitive, the number of adjacent surfaces and the connectivity of the surface, a value can be assigned to each primitive. Larger values indicate a better possibility of a good fit. Thus, at each stage of the search the scene primitive with the largest value is chosen as the next candidate for matching. The approach is similar to those employed by Oshima and Shirai35, Sheffield’s AIVRU”.2h,27 and Murray24.25. ~019 no 6 december 1991

Each of the measures has a weighting factor associated with it. Experimental results show that the weighting factors should be chosen to reflect the fact that: The face of the largest area should be taken first so long as it is not subject to worse segmentation errors than a face of lesser area. Planar, cylindrical and spherical faces should be taken in that order as long as the face of simpler type does not have a worse fit than a more complex face. The face with the most neighbouring faces should be chosen first between like faces, since it is more likely to provide a satisfactory match. The value of a surface is calculated as follows: V(a,e,t,c)=w,a+w21e+w3)t+w4c

(7)

where wl, w2, w3 and w4 are weighting factors, a = area of face, e = quality of fit to face, t = surface type ( = 1 if surface is a plane, = 2 if surface is cylindrical, = 5 if surface is spherical, = 10 if other surface type), c = number of connected primitives. The quality of fit of the face e is the largest absolute distance of any point fitted to the face to the approximation of that face. The values of the weighting factors were chosen based on experiments on simple objects to meet the criteria above. In practice, w3 and w4 are set to one, w1 corresponds to the reciprocal of the area of the largest segmented face, and w2 corresponds to the maximum value of e for any segmented face. The measures can be calculated during segmentation at little extra cost. It was also found that best results are obtained when the transformation is calculated completely early on in the search. If this occurs then the search space can be sufficiently reduced since mismatching planes can be more effectively eliminated once the complete transformation is defined. For this complete transformation to be defined the first two or three primitives chosen should be linearly independent. This is best when the chosen primitives are as near orthogonal as possible. Additional local constraint

It was noted that even on a relatively simple object such as the widget (see Figure 6b) that certain impossible solutions were followed for some length of the search path before being rejected. This is due to the fact that since only the estimate of the rotation is employed in the constraints for most of the search certain symmetries cannot be eliminated (see above). The main problem arises when two parallel planes which occlude each other occur in the model. This is a common occurrence, particularly when symmetries occur in an object and ambiguities about 180” arise. Each plane could be potentially matched to a scene plane, but obviously only one of the two planes is visible, and hence a candidate for the match. A match with the wrong parallel plane causes the subsequent rotation estimate to be 180” out. Thus, some time may be wasted trying to match other planes under this incorrect assumption. Similar mismatching problems also arise when parallel planes which do not or partially occlude each other occur, as illustrated in the (2D) Figure 10. 395

sight Paralle.1 Planes

53I Figure 10. Examples

of non-occluding

parallel planes

The solution adopted to help solve this problem, and hence increase the speed of the search, is to use the area of the respective faces to aid in discarding a wrong match. This is done at the same time as the other local constraint that eliminates grossly mismatching faces (see above) is tested. The areas of the faces in the scene are already known, and the areas of the model faces can easily be precomputed. If a scene face is greater in area than that of a mode1 face then there is no possibility of them matching. The exact constraint of the areas being equal cannot be used since occlusion of scene faces may occur. It should be noted that calculation of the area of scene faces is not trivial, as they are seen in projection. However, this constraint is very cheap to implement and does not greatly increase the computation required by the method. Extension to primitives other than planes

By restricting attention to just planar surfaces the potential for inspecting or even just identifying objects is severely limited. Segmenting scenes into higher order primitives that are reliable is often very difficult and time consuming; as shown by Hebert3” and Faugeras3’*32 for quadric surfaces. Since these are the ‘next simplest’ type of genera1 surfaces after planes, it is assumed that any approach to inspection that uses a general description of a surface of this type, or any other higher order surfaces, will be unreliable. The main problem in fitting general higher order algebraic surfaces to a set of points is that there are a larger number of ‘free’ parameters that need determining, i.e. a plane (z = ax+by +c) has three’ free parameters, namely u, b and c, whereas a general quadric

that an overdetermined system of equations is formed. In that case quadrics more than eight points are needed while a plane requires only three points. This leads to two problems. First, the time required to solve the system of equations is increased. Consequently, the time to segment the scene is vastly increased since very many systems of equations have to be solved. Second, larger errors in segmentation will result. Therefore, an approach is required that can reduce the number of free parameters to be determined. One method is to use the principal directions of a surface to control the match. By looking for curved surfaces of specific types the number of free parameters is reduced, and the number of natural constraints that can be employed is increased. Thus we restrict our attention to quadrics of a certain type and match each type independently. Since a vast majority of the parts that require inspecting can be modelled using planes, cylinders and spheres’ it is natural to explore the possibilities of using these as primitives in the matching process. The case for spherical and cylindrical surfaces is described: Spheres: If the centre of a sphere C can be estimated, and this centre is not at the origin of the transformation, then the vector from the centre to the origin may be used in the matching process in a similar way to the plane normal in plane matching. Spheres of different radii r can be eliminated by a simple compa~son of their respective radii. it can easily be seen that a sphere has four free parameters when expressed in the form (X- a)* + (y b)2 + (z -c)” = r2. Cylinders: a cylinder can be represented as follows (see Figure 1lb):

a unit vector n parallel to the axis of symmetry; a vector d from the origin to the closest point on the axis of symmetry; the radius of the cylinder; two scalars, A, and A2, denoting the start and points of the cylinder along the axis of symmetry, i.e. stati_vector = d + AIn. An infinite cylinder (ignoring h, and A2) has six free parameters (two for the direction of n, the vector d and. the radius). By using the vectors d and n in a similar way to the plane normal in plane matching, cylindrical surfaces can be matched.

planes we can emptoy natural interpretations of the planes parameters to derive an expression that can be minimized in the least squares sense in order that an approximation of the plane be made. For example, when a plane is expressed as ax+by+cz+d = 0 we can interpret (a, b, c) as the normal vector and d as the minimum distance from a point to a plane. However, in the case of quadric surfaces the parameters a, 6, . . . j have no simple geometric meaning, and so no natural interpretations exist to form the least squares minimization. Various artificial inte~retations need to be designed to constrain the minimization32. Using least squares techniques it is necessary that more points than the number of free parameters are present so 396

0 r

(ax’ + by2 + cz2 + 2dxy + 2exz + 2fyz + gx + hy + iz +j = 0) has nine free parameters, e.g. a/j . . . io. In the case of

C

a

d

I

40

b Figure II. (a) Representation of a sphere; (b) representation of a cylinder image and vision computing

Since the fit to a non-planar primitive is likely to be less accurate than a planar primitive, cylindrical and spherical p~mitives are given a lower priority in the hy~thesis stage. This is the purpose of the surface type parameter outlined above. This approach is easily extensible to include other primitives such as cones and toroids. The ability to match cylinders and spheres has been added to the matching algorithm and works well on theoretical data. Although reliable segmentation of these primitives is not available as of yet, work is being carried out to achieve this segmentation

2 A small modification to the segmentation procedure is required. Attempts to segment a curved surface may result in instead many small planar facets. Since any inclusion of these facets in the matching process could at worst cause failure or otherwise introduce spurious errors, it is important that they are rejected early on. We therefore insist that a plane must be larger than a certain specified area for it to be considered. If for some reason a larger planar facet avoids this detection, it is anticipated that the resulting error will be so large that it will fail the first cut-off constraint.

.

Match when possible primitives

using only reliably segmented

For reasons given below it was found that the location and orientation of the object is well defined after most of the segmented primitives of good fit have been considered for matching. Using primitives with a worse fit lower down in the search just ~ont~buties unnecessa~ errors to the estimate of the transformation. Thus a method that can optimize the fitting of only the good primitives will not only produce better transformation estimates, but will enable the method to be quicker, since poorly fitting primitives need not be considered for finding the transformation. As mentioned above, only reliable segmentations of planes are currently available. Thus we insist that the range of inspectable parts is restricted further in that at least three good planar primitives can be extracted from the scene to estimate the transformation. This appears to be no real drawback in that most of the objects that are considered for inspection at the moment fall into this class. The term ‘good’ is applied in the sense described above. However, to accommodate all possible views of objects of the classes described previously, the need for matching and locating using non-planar primitives must still be accommodated. After matching and location of the object has been achieved by using, as far as possible, only reliable segmented planes, the location of unmatched primitives is then posed as an inspection problem (see below). Only minor modifications to the existing matching method are needed to implement the strategy of matching only using reliable planes: 1 A cut-off is introduced when the order of matching pairs is evaluated, described in above, such that if there exist more than three places which provide a potential estimation of the transformation then no non-planar primitives need be considered. This cut-off can be extended further to that poorly segmented planes need not be considered to estimate the transfo~ation. It should be noted, however, that as many ‘good primitives as possible should be considered for the estimation so that the system of linear equations required to estimate the transformations is as overdetermined as possible to optimize the least squares approximations. This cut off is simple to implement in practice once more than three scene planes have been matched; subsequent planes are also matched if the error in their fit is below a certain threshold. We have found that satisfactory results are obtained when the value of the threshold is 0.1 mm for normal inspection accuracy. voE9 no 6 december 1991

Alluring for some mismatching faces to be rejected from the match

If a ‘rogue’ face is obtained from the scene and a match is sought with the model, no suitable match may be found. The source of rogue faces could be some background planes having depth values assigned to them, some other equipment in the scene (e.g. the robot or the turntables) rendering depth readings or, under less controlled conditions, the segmentation could be the source of errors. Another possibility is that the object being viewed has been incorrectly manufactured. In this case it is still useful if the matching process is continued and does not just stop at this point, so that other features can be inspected and the whole object not be immediately rejected as a bad match. The incorrect faces will clearly be rejected during the inspection phase. Accommodating such mismatching faces in the matching strategy is straightforward. Instead of having ra branches, where n is the number of model faces, coming from each node an additional null branch is added to each node. The tree search proceeds as described before, except that if no model primitive matches the scene primitive then the scene primitive is assigned to the null branch, and the search continues. Clearly, it is possible that a total mismatch of model and object could occur, resulting in a completely traversed set of null branches throughout the tree. This would be inefficient in both computational time and space considerations, and so the set of null branches in the search tree is not allowed to grow more than half the total possible depth of the search (i.e. half the number of scene primitives). This is not unreasonable, since if more than half the scene primitives have not been matched then a satisfactory match is unlikely. Attempt to recover matched planes

position

information

from

two

As mentioned previously, three linearly-independent planes are required to estimate the translation, and thus completely define the transformation from model to world coordinate frames. However, due to the relatively narrow field of view of the cameras, it is possible for the scene of the viewed object to contain a set of planes that are not linearly-independent in three dimensions. This is illustrated in Figure 12 of the widget, where each individual camera can see three linearly-independent planes, but the combined view (where only points common to both views have depth estimates) does not, However, so long as two linearly independent primitives are visible then a cruder estimate of the translation is possible. While this estimate may not be sufficiently accurate for performing inspection, it is still useful to recover some idea of the position of the object 397

3 A3 is estimated by minimizing the distance between the respective centre of gravities of between matched scene and model primitives along the direction of e3.

Results View from master camera

View from other camera

Common view of both cameras Figure 12. Example of the gathering linearly independent plane

of three non-

for many applications. In particular, for inspection applications-” the object can be repositioned to a new view likely to provide a more accurate match with some degree of certainty rather some arbitrary or random position that may not guarantee a better view. The singular value decomposition (SVD) of the matrix A, formed by rotating matched model primitives into the camera coordinate frame, gives an orthogonal basis for the set of formed vectors in A - the eigenvectors of ATA (Matrix V in the SVD - see above) and the corresponding eigenvalues (Matrix C). If there exists a zero eigenvalue in C then there is no vector present in A that is linearly-dependent to the corresponding eigenvector. Therefore, by counting the number of non-zero eigenvalues in Z the number of linearly independent vectors of A can be determined. If there are three non-zero eigenvalues then there exits a linearly independent set of vectors and the translation is determined as before. If there is only one non-zero eigenvalue then an error must have occurred since the estimation of the rotation requires two linearly independent vectors. If there are two non-zero eigenvalues then the set of vectors of A are do not have any vectors linearly-dependent in one of the eigenvector directions. In this case the translation is estimated as follows: Find the eigenvector that is linearly independent of all the vectors in A - the eigenvector corresponding to the zero eigenvalue. Work out the translation of the vectors in the two remaining eigenvector directions using the SVD of A projected into the two dimensional space defined by these two eigenvectors. This amounts to simply dropping the third eigenvector from the above calculations, since the translation vector may be expressed as: t=h,e,

+h2e2+A3e3

03)

where el, e2, e3 are the three eigenvalues and A,, AZ, A3 require estimating. Assuming e3 is the eigenvector linearly independent of all vectors of A, then A, and A2 estimated from the SVD and A3 is estimated separately. 398

The results from adding the improvements are assessed using artificial and real data. Before the matching algorithms are assessed the results from the segmentation algorithms are briefly discussed. Brief analysis of segmentation The se mentation algorithms are detailed elsewhere54-9 ‘, and are summarized in Appendix A. However, some important results from studies of these algorithsm on artificial and real data are reviewed. Results in this section and subsequent sections are presented on artificial and real data. The artificial data was produced by adding Gaussian noise to synthetic images (see Figure 6), as described above. With both types of data the segmentation algorithms proved to be robust, and provided good approximations to actual model regions. The rms fit of data to each region was low and proportional to the added noise. The average rms errors in overall (model to scene) rotational and translational position, using artificial data, on several samples proved to be proportional to the added noise. In particular, at large noise values (more than 1 mm) the error in rotational positioning was at worst 0.32”. With real data the average rms error of the planes fitted to the depth data is 0.28mm. However, the average error in the relative angles between the segmented planes on the widget is 1.2”. This error in angles is larger than the error in relative angles of about 0.32” produced in artificial depth data at the largest tolerable noise levels (added noise > 1 .Omm). There is clearly a discrepancy between these two results when the rms error in the position error corresponds to added artificial errors of about 0.3mm but the angle errors correspond to errors more than three times as large. Since the segmentation can be shown5h*57to be stable for large artificial errors and the errors in the real depth data are small, it is reasonable to assume that the errors in the relative errors of the angles of segmented planes are not due to failure of the segmentation. It has already been shown above that calibration errors in the camera set up can lead to errors in locating known points. It is difficult to analytically assess the calibration accuracy, since all subsequent measurements depend on the calibration. It is, however, possible to assess how well the calibration performs in the real world. Two measures are possible: the position of the calibration points sensed after the initial calibration (hence a radius of ambiguity in the position of these points can be calculated); how well the calibrated system can sense known points in real world. The current camera set up has an rms error of 0.5 mm in the positioning of known points in world coordinates, and a radius of ambiguity in pixel placement of 0.4mm. It is not unreasonable to assume that calibration errors account for the errors in the relative angles of planes observed above. This fact is analysed further elsewhere”. image and vision computing

Comparison of original Faugeras Hebert and modified matching algorithms on segmented artificial data The results presented in this section correspond to taking the segmentations produced above, and matching to geometric models stored in the computer. The output of the matching routines is a list of matched planes and a transformation (a rotation matrix and a translation vector) relating the object model to the real object. Since the object’s pose was known in the production of the artificial depth maps, the exact transformation is known. The errors in the estimated transformation are easily calculated, and are shown in Figures 13a and 17b. The modified matching algorithm with all the constraints was applied to the segmentations and models, as was the Faugeras Hebert algorithm. No cylind~cal primitives were considered for the match. Since both methods adopt the same strategy for estimation of the rotation no comparison can be made between them. The average errors in estimating the rotation are summarized in Figure 13. Since the estimates of the translation parts in both algorithms are dependent on good estimates of the rotation, it is important to note that errors in the matching of model surface normals to scene surface normals do not grow too large until the added depths errors have grown much larger than the expected hardware errors. The errors in the translation are now considered. The errors in the matching of model to scene planes for the widget are summarized in Figure 13. The range of the widget was 707mm.

Average error in normals 0

1.00

a

Error in depth

reading

(mm)

Error in translation

0.0

b

io.2

0.4

Error in depth

0.6

0.6

1.0

reading (mm)

Figure 13. (a) Errors in the estimation of rotation in matching; (b) errors in translation of matching the widget ~019 no 6 december 1991

As can be seen from the results, the errors in the matching algorithm remain fairly constant for a given depth reading error, even though the range of the cube is greater than that of the widget. It can easily be seen by comparing the corresponding errors in rotation and translation for a given depth reading that the errors in the translation for the Faugeras Hebert method are approximately pro~rtional to &an 8, where d is the perpendicular distance of a plane to the origin and 8 is the error in the angle. Another important result here is that the predicted error in the position of an object is within the desired inspection tolerance of 0.254mm. Even at larger than expected hardware errors, the modified matching algorithm proved capable of providing such data. At expected system errors and expected working range ~0.5-0.7m) the Faugeras Hebert method also proved satisfactory, although not as accurate as the modified matching algorithm. However, after about added depth errors of about 0.5 mm the errors in the Faugeras Hebert translation exceed this tolerance. This implies that the modified matching algorithm will be more reliable should smaller inspection tolerance limits be required. The time savings employing all the modified constraints compared to the Faugeras Hebert method were about 50%. Analysis of matching on segmented real data Results in this section are presented on the same real world views of the widget that were used for the segmentation analysis on real data. Once again, it is dif~cult to infer exactly how well the matching methods have performed on real data, since the exact positions of objects are not known. However, two analyses are possible. First, the fit of the model (under complete transformation) to the data in the camera frame can be compared. The average error in the estimation of the rotation was 0.9” and. the average error in the translation estimate was 0.2mm for the modified matching algorithm and about 4mm for the Faugeras Hebert method (once more - dtan @). It has been shown that there was a discrepancy in the relative angles between segmented faces. This accounts for the large error in the rotation and the translation as well. It is anticipated that, as illustrated by artificial data, should the rotation error be reduced the translation error will be substantilly reduced, thus enabling more accurate measurements. Errors of this type are illustrated in Figure 14, which shows the visible surfaces from the model of the widget positioned, under the matching transformation, in camera coordinates and superimposed onto a depth map of the widget. With the picture suitably enlarged slight mismatches may be observed especially at the edges. The second analysis of the matching results may be performed in world coordinates. The widget was placed on the ground plane (z = 0) and so the z position of the vertices of the model (transformed into world coordinates) can be analysed, since their exact positions may be assumed. The average error in the z positioning of the model was 0.5 mm. This is consistent with the errors calculated for the calibration analysis, where the position ground plane was found to be in error by 0.5mm. This once more shows that the most likely 399

inspection tolerances is not quite possible. However, it has been shown by considering segmentation and matching results that the errors in the calibration of the cameras greatly contribute to this problem. These errors should also be resolvable, and are currently being investigated at the SRC. From the artificial data analysis it can be assumed that if the calibration errors can be sufficiently reduced then matching to within adequate inspection tolerances should be achievable.

INSPECTION First the inspection task is discussed. The effects that errors arising from the various components in the system have on the outcome of inspection are briefly mentioned. Finally, results are presented on the adopted inspection strategy.

Inspection task Figure 14. Widget with transformed imposed onto depth map

model

super-

Given that the location and orientation of the object can be found to within a small tolerance, then known visible features extracted from that particular view of the model can be used to verify the position of features in the scene, and so some inspection tasks are possible. Moreover, once the position of the object is known, by using accurate turntables the object can be rotated to give different views of the object in which the position of the object is already known. Therefore, after the initial match, primitives that should be visible in the new view can be hypothesized from the object model. Subsequent matching of views is simplified to just testing (inspecting) the segmented scene for correct occurrences of the hypothesized primitives. The inspection tasks that we are concerned with are mainly due to tolerance errors5875y. That is to say, determining: 1 That the dimensions of the features are correct. 2 That the positioning of the features are correct. 3 That the relative positioning of features is correct, e.g. certain surfaces are parallel, perpendicular, etc.

Figure 1.5. Transformed model superimposed view of widget

on actual

explanation for this misplacement of the model is due to calibration errors and not matching (or segmentation) errors. Errors of this type are illustrated in Figure 15, which shows the model of the widget positioned, under the matching transformation, in camera coordinates, and superimposed onto the observed view of the widget from the master camera.

Conclusions By adopting and improving the strategy of planar segmentation and matching on reliable predicted errors will provide results accurate enough to enable the inspection of an object using the approach described. With the current vision system set up it has been shown that, using real data, the matching to within 400

In this paper we limit the inspection task inspection of a single scene. Thus only limited tion can take place, i.e.:

to the inspec-

1 Are the dimensions of visible primitives correct to the desired tolerance? 2 Are the locations of visible primitives correct to the desired tolerance? 3 Are the positions of other visible features relative to a visible primitive correct? The scene is inspected

as follows:

From the matching of the model with the scene a list of matched primitives and the transformation are given. This transformation is then applied to the model, and a list of visible primitives is compiled. Knowing the hypothesized position of a visible model primitive and the dimensions of the depth map, the corresponding position of this primitive in the scene can be postulated. Inspection takes place by testing the points in the depth map contained within each primitive to be within tolerance with respect to the primitives. image and vision computing

100

View of Object

90

Vision system

80 70 Database of Geometric Object Models

Lit of primitive suffaces and surface properties

% of

60

points within tolerance

so 40 30 20 10

Inspection

I

D3aiIs

Figure 16. Errors system

0

* Tramformation ( Rotation,

’

0

0.08

0.15

0.31

0.61

0.92

Errors in depth readings (mm)

Translation )

introduced

0.04

by components

of the

By posing the inspection of a single scene like this the three inspection tasks above are all clearly satisfied.

Errors in the system The errors introduced at each stage of the system are summarized in Figure 16. Each error that occurs will also cont~bute to the error produced by subsequent stages. Thus the outcome is affected by a combination of camera, segmentation and matching errors. The inspection process can therefore also be viewed as a good test for the accuracy of the whole system. As shown in previous sections, the errors due to matching can be kept to a minimum. Also, knowing the vision system camera set up allows the errors in the segmentation to be predicted.

Figure 17. Analysis

of inspection

of a widget with a hole

Results on real data Results in this section are presented on real data gathered from several unique views of the widget. In some cases a one penny coin (height 2.5 mm) was placed on a face in order to simulate a defect in the object. An example of a penny placed on the widget may be seen in Figures 14, 15 and 18. It has already been shown that errors in the position-

PttfiiiS

IN

LRRDA

Results Results on artificial data The results presented in this section correspond to artificial depth maps, with added noise being run through all the processes described to produce an estimate of location and orientation of a model. This estimate is then applied to the geometric model, and the visible primitives are inspected in the original depth map. This allows the effects of noise both inherent in the system and due to the camera range to be studied. Tests have been performed on the same data as in the previous sections. Results are presented using three different tolerances, namely 0.25 mm (0.01 inch) 0.124 mm (0.005 inch) and 0.025 mm (0.001 inch) to reflect the fact that typical inspection tolerances are 0.25 mm and detailed inspection tolerances are usually no less than 0.025 mm. Typical results are summarized in Figure 17, As can be seen from the results, inspection is possible to within the normal inspection tolerance of 0.254mm with the expected system errors (at least 99% of the object within tolerance with depth errors of 0.15mm and working at the typical range of 0.7m) for the system described. However, moving the object further away reduces the number of points within tolerance, limiting the possibility of successful inspection. ~019 no 6 december 1991

Figure 18. Inspection of objects to tolerance

401

It has also been shown that using the current camera set up inspection is possible with position tolerances of 0.8 mm. The camera calibration inaccuracies have been shown to limit the accuracy to which inspection can be performed. Should further developments in the calibration accuracy be achieved, then the artificial results show that inspection to more realistic tolerances can be expected. It is also hoped that this paper provides useful information in that the accuracy to which a vision system has to be constructed and calibrated in order to perform similar matching and inspection tasks has been demonstrated. Figure 19. Segmentation of a widget

ing of objects under the current system set up arise due to inaccurate camera calibration. Clearly, these errors have a great effect on the success of the inspection. The testing for tolerance in position of an object is clearly dependent on the matching transformation, and is affected by these errors. In all cases, all data points can be measured to within a position tolerance of between 1 mm and about 0.8mm, but as this tolerance is reduced further points begin to fall out of tolerance. At a position tolerance of 0.75 mm 2% points are out of tolerance, and at 0.5 mm 10% of points are outside tolerance. The points that fall out of tolerance are initially at the edges of segmented primitives and grow inwards towards the centre of the primitives, as illustrated in Figure 9e. This effect can be explained by the fact that the normals of the segmented primitives are in error by around l”, and consequently affect the matching results. Since the centres of gravity of the segmented planes are matched to the model planes, in order to minimize the translation error, all points close to the centre od gravity are at a minimum distance from the plane. If the plane is orientated at an incorrect angle then the further from the centre of gravity a point lies, the larger the error in the distance from the plane the point is.

FUTURE

WORK

The ability to inspect a single scene only has been demonstrated. Intended future work consists of developing strategies to fully inspect the object by taking data from several scenes so that complete inspection may be achieved. Some of this work may be found elsewhere53*57. The need to fully incorporate non-planar surfaces into this work also requires further development. ACKNOWLEDGEMENTS We would like to thank Michael Gay and Andy Page, for providing the segmentation algorithms and many helpful suggestions. Dr Marshall would also like to thank all the staff of the Advanced Information Processing Dept., SRC, British Aerospace for providing the initial and continuing support for this project. REEFERENCES Davies, P and Hutber, D A survey of possible methods for automatic inspection of three dimensional objects BAe, S.R.C. internal report no. JS

10631 (June 1986) D and Laurendeau, D ‘3D sensing for industrial computer vision’ Advances in Machine Vision, Architectures and Applications, SpringerVerlag, NY (1987) Doemens, G. et al. ‘A fast 3D sensor with high dynamic range for industrial applications’, Proc. Poussart,

CONCLUSIONS We have demonstrated segmenting planes, and a restricted range of non-planar surfaces from dense depth maps provides data upon which geometric model-based matching and inspection may be performed. The best existing matching approach that uses such data, the Faugeras Hebert method, has been justified to be suitable, after improvement, for inspection purposes. We have shown that using this matching strategy and the given segmentation methods that general inspection is possible. Inspection to tighter tolerances than 0.25 mm may be possible with the current system, but very detailed inspection, to within 0.025 mm is probably not possible, as only 25% of the depth points fall within tolerance. In addition, an overall approach to complete object inspection has also been outlined. These results imply that future research into developing such a system is worthwhile. It has been shown that under normal expected working conditions inspection to normal tolerance limits is possible, but that detailed inspection to smaller tolerances may not be possible using the outlined hardware and software. 402

ROVISEC Conf. (1986) Hutber, D ‘Automatic inspection of 3-D objects using stereo’ SPIE vol. 850: optics, illumination ilmm_ig;ge sensing for machine vision II (1987) pp Boulanger, P et al. Object input for CAD/CAM using a 3D laser scanner, National Research

Council of Canada Lab. for Intelligent NRCC Memo No. 25446 (1985)

Systems,

Porrill, J, Pollard, S B, Pridmore, T P, Mayhew, J E W and Frisby, J P TZNA: The Sheffield AZVRU vision system, Sheffield University

AIV-RU No 031, UK (1987) J Private communication. Quality Assurance Dept., British Aerospace (1988) Brown, C M ‘Some mathematical and representational aspects of solid modelling’, IEEE Trans. PAMZ, Vol 3 No 4 (July 1981) pp 444-453 Baer, A, Eastman, C and Henrion, M ‘Geometric Bremner,

image and vision computing

10

11

12

13

14

15

modelling: a survey’, Comput. Aided Design, Vol 11 No 5 (September 1979) pop 253-272 Henderson, T C ‘Efficient 3-D object representations for industrial vision systems’, IEEE Trans. PAMI, Vol5 No 6 (November 1983) pp d609-617 Requicha, A A G M~themafi~aL modes of rigid solid objects, University of Rochester, Production Automation Project, Technical memorandum 28 (November 1977) Requicha, A A G and Voelcker, H B ‘Solid modeling: a historical summary and contemporary assessment’, IEEE Comput. Graphics Appl., Vol2 No 2 (March 1982) pp 9-24 Requicha, A A G and Voelcker, H B ‘Solid modeling: current status and research directions’, IEEE Comput. Graphics Appl., Vol 3 No 7 (October 1983) pp 25-37 Silva, C E Alternative de~nitions of faces in boundary representation of solids, University of Rochester, Production Automation Project, Technical memorandum 36 (November 1981) Tsai, R Y A versatile camera calibration technique for high accuracy machine metrology using off the shelf TV cameras and Ienses, IBM Research report

RC 51342 (October 1, 1985) 16 Tsai, R Y ‘An efficient and accurate calibration technique for 3-D machine

camera vision’,

27

28

29

30

31

32

33

34

Automated inspection and high speed vision architectures (1987) pp 11-17

18 Huther, D and Wright, S M ‘A digital camera and real-time geometric correction for images’, Proc. Alvey Vision Conf, Cambridge, UK (September 1987) pp 151-157 19 Besl, P J and Jain, R K ‘Three-dimensional object Recognition’, ACM Comput. Surv., Vol 17 No 1 (March 1985) pp 75-145 20 Chin, R T and Dyer, C R ‘Model based recognition in robot vision’, ACM Comput. Surv., Vol 18 No 1 (March 1986) pp 67-108 21 Grimson, W E L and Lozano-Perez, T Model based recognition and localization from sparse range or tactile data, MIT AI Memo no. 738, MIT, MA,

22

USA (August 1983) Grimson, W E L The Combinatorics

23

MIT, MA, USA (April 1984) Lozano-Perez, T and Grimson, W E L Recognition

of local constraints in model-bused recognition and localization from sparse data, MIT, AI Memo no. 763,

and localization of overlapping parts from sparse data, MIT, AI Memo no. 842, MIT, MA, USA

24

25

(1985) Murray, D W ‘Model-based recognition structure from motion’, Image & Vision Vol 5 No 2 (May 1986) pp 85-90 Murray, D W ‘Model-based recognition shape alone’. ~omput. Vision, Graph. Process (1987)

26

using 3-D

vol9 no 6 december 1991

35

36

37 38

39

40

41 42

Comput.,

43 using 3-D & Image

Pollard, S B, Porrill, J, Mayhew, J E W and Frisby, J P ‘Matching geometric descriptions in three space’, Proc. 2nd Aivey Vision Conf., (Expanded Version Sheffield Univ. AIVRU No 022) (1986)

Faugeras, 0 D and Hebert, M ‘A 3-D recognition and positioning algorithm using geometric matching between primitive surfaces’, Proc ZJCAZ (1983) pp 996- 1002 Faugeras, 0 D and Hebert, M ‘Object representation, identification and positioning from range data’, Proc. 1st Int. Symp. Robotics Res. (1984) pp 42.5-446 Faugeras, 0 D aud Hebert, M ‘The representation, recognition and positioning of 3D shapes from range data’, in Rosenfeld, A (ed), Machine Intelligence and Pattern Recognition: Techniques for 3-D machine perception, North-Holland, Amsterdam

IEEE Int. Conf. Comput. Vision & Patt. Recogn., (1986) pp 364-374

17 Hutbcr, D ‘Improvement of CCD camera resolution using a jittering technique’, SPIE vol. 849:

Pollard, S B, Fridmore, T P, Porrill, J, Mayhew, J E W and Frisby, J P Geometric model~ing from multiple stereo views, Sheffield University AIVRU No 028, UK (1987) Faugeras, 0 D, Germain, F, Kryze, G, Boissonnat, J D, Hebert, M, Ponce, J, Pauchon, E and Ayache, N ‘Towards a flexible vision system’, in Pugh, A y;) l??tot Vision, Springer-Verlag, NY (1983) pp

44

45

(1986) pp 13-53 Faugeras, 0 D and Hebert, M ‘The representation, recognition and locating of 3-D objects’, Znt. J. Robotics Res., Vol 5 No 3 (Fall 1986) pp 27-52 Hebert, M and Kanade, T ‘The 3-D profile method for object recognition’, Proc. IEEE Comput. Vision & P. H. Recogn. (1985) pp 458-463 Ayache, N, Faugeras, 0 D, Faverjou, B and Toscani, G ~Matching depth maps obtained by passive stereo’, Proc. 3rd Workshop on Comput. Vision: Representation and Control (1985) pp 197-204 Oshima, M and Shirai, Y ‘Object recognition using 3-D information’, IEEE Trans. PAMZ, Vol5 No 4 (July 1983) pp 353-361 Rosenfield, A, Hummel, R and Zucker, S ‘Scene labelling by relaxation operations’, IEEE Truns Systems Men & Cybern., Vol6 (1979) pp 420-433 Hummel, R and Zucker, S ‘On the foundations of relaxation labelling processes’, ZEEE Trans. PAMZ, Vol 5 No 3 (May 1983) pp 267-268 Bhanu, B ‘Representation and shape matching of 3-D objects’, IEEE Trans. PAMI, Vol 6 No 3 (May 1984) pp 340-351 Silberberg, T M, Harwood, D A and Davis, L S ‘Object recognition using oriented model points’, Proc. 1st Artif. Intell. Applic. Conf., Denver, CO (December 5-7 1984) pp 645-651 Boyter, B A and Aggarwal, J K ‘Recognition of polyhedra from range data’, IEEE Expert (Spring 1986) pp 47-59 Ballard, D H and Sabbah, D ‘Viewer independent shape recognition’, IEEE Trans. PAMI, Vol 5 No 6 (November 1983) pp 653-659 Ballard, D H and Brown, C M Computer V~ion, Prentice Hall, NJ (1982) Milgram, D L and Bjorkland, C M ‘Range image processing: planar surface extraction’, Proc. 5th int. Conf. on Patt. Recogn., Miami, FL (December 1-4 1980) pp 912-929 Duda, R 0, Nitzan, D and Barret, P ‘Use of range and reflection data to find planar surface regions’, gF_;7kurzs. PA MI, Vol 1 No 3 (July 1979) pp Hebert, M and Ponce, J ‘A new method

for 403

segmenting 3-D scenes into primitives’, Proc. 6th Munich, West Germany (October 19-22 1982) pp 836-838 46 Besl, P J and Jain, R K ‘Image segmentation through variable-order surface fitting’, Proc. Vision ‘87 Conf., Detroit, MI (June 8-11 1987) 47 Osborne, E E ‘On least squares solutions of linear equations’, J. ACM, Vol 8 (1961) pp 628-636 48 Golub, G and Kahan, W ‘Calculating the singular values and pseudo-inverse a matrix’, J. SIAM Numer. Anal, Vol 2 No 2 (1965) pp 205-224 49 Golub, G ‘Numerical methods for solving linear least-squares problems’, Numerische Mathematik, Vol 7 (1965) pp 206-216 50 Golub, G and Reinsch, C ‘Singular value decomposition and least squares solutions’, Numerische ~athematik Vol 14 (1970) pp 403-420 51 Bauer, F L ‘Elimination with weighted row combinations for solving linear equations and least squares problems’, Numerische Mathematik Vol 7 (1965) pp 338-352 52 Hurley, S Solving certain least squares problems with rank deficiency, MSc dissertation, University of Reading, UK (September 1985) 53 Marshall, A D ‘Automatically inspecting gross features of machined objects using threedimensional depth data’, Proc. SPIE OEI Int. Conf. on Patt. Recogn.,

54

55

56 57

BOSTON 90: Symposium on Advances in Intelligent Systems Boston, MA (4-9 November 1990) Page, A Segmentation a~gor~th~, BAe, Sowerby

Research Centre Internal Report No. AOI/TRl BASR880201 (March 1988) Gay, M Private communications. AIP Dept, Sowerby Research Centre, British Aerospace, Bristol (1988/89). Gay, M ‘Segmentation using region merging with edges’, Proc. 5th Alvey Vision Conf., Reading, UK (September 1989) pp 115-120 Marshall, A D The Automatic Inspection of Machined Parts using Three-Dimensional Range Data and Model Based Matching Techniques, PhD

Thesis, University of Wales College of Cardiff, UK (1989) 58 Requicha, A A G ‘Toward a theory of geometric tolerancing’, Znt. .I. Robotics Res., Vol 2 No 4 (1983) pp 45-60 59 Truslove, K C E ‘The implications of tolerencing for computer-aided design’, Comput. -Aided Eng. J. (April 1988) pp 79-85

APPENDIX

A: SEGMENTATION

A brief description of the segmentation method is included for the sake of completeness. Further details of the segmentation may be obtained from Page54 and Gay 55s6 An example of the segmentation is shown in Figure 18. The segmentation process may be divided up into four distinct phases: 1 An initial stage where a mask of variable dimension is passed over the entire image and small planar approximations are made. 404

A coalescing stage: the COALESCE algorithm quickly reduces the intially large number of regions (potentially more than 65,000) to approximately 3,000 regions, which can be more rigorously tested by subsequent stages. The method is based on the physical phenomenon of condensation in which regions are gradually grown over the scene as a whole. The basic procedure of the algorithm is to iteratively group regions together under the condition that the difference between regions is less than some allowable size, which is increased on each iteration until a specified limit is reached. A forcing stage: the FORCE algorithm uses a more rigorous merging procedure. The basic idea is to examine pairs of regions and consider the ‘forces of attraction’ and ‘repulsion’ between them. The attractive forces are based on the similarity of the two regions, whilst the repulsive forces are based on the properties of the edge boundaries between the regions 54s5. Only regions with positive resultant forces of attraction are merged together. The merging process is based on simultaneous region growing. A featuring stage: various features (e.g. plane equation, region area, error in plane fit, number of connected regions, centre of gravity, etc.) are derived from the segmentation. Throughout the segmentation use is made of a region map. A region map is an array of numbers where each array element (pixel) is associated with its corresponding pixel from the depth map (the region map has the same dimensions as its corresponding depth map). The number associated with each (region map) pixel corresponds to a unique region obtained in the segmentation process.

APPENDIX B: EXAMPLE OF PSEUDO-INVERSE METHOD

ERRORS

IN

The pseudo-inverse method for least squares minimization can be prone to rounding errors and noise. This is a well known problem in numerical analysis and arises since the conditioning of ATA is worse than that of A4’. This frequently leads to difficulties which cannot be remedied once ATA has been formed. Detailed descriptions and error analysis of the pseudo-inverse method are well documented47-50. This problem is illustrated by the following example: Consider, when:

for the purpose of illustration,

1

1

1

&

0

0

0

&

0

0

E

A= [ 0

the case

1

where F is an error, then:

I image and vision computing

Now, if E is such that E’
ATA=

1

1

1

1

1

1

[ 1

1

1

1

In this case ATA is singular and, therefore,

~019 no 6 december 1991

no inverse

of ATA exists. Thus no solution to the minimization can be found. It is important to note that errors, of the form in the example, are common, since errors arising from the estimation of the rotation are introduced to A, because A = [w’], . . . . . . wINIT, where W’i= Rv’i, i= 1 .,. N and elements in R may contain errors. However, an error may be introduced to each element of w’~, and hence each element of A474”.