DOCTORA L T H E S I S
Department of Engineering Sciences and Mathematics Division of Mathematical Sciences
Computational Methods for
Shape Verification of
Free-Form Surfaces
Per Bergström
ISSN: 1402-1544 ISBN 978-91-7439-301-9 Luleå University of Technology 2011
Per
Bergström
Computational
Methods
for
Shape
Ver
ification
of
Fr
ee-F
or
m
Surf
aces
Computational Methods for
Shape Verification of
Free-Form Surfaces
Per Bergstr¨
om
Department of Engineering Sciences and Mathematics
Division of Mathematical Sciences
Lule˚
a University of Technology
SE-971 87 Lule˚
a, Sweden
Printed by Universitetstryckeriet, Luleå 2011 ISSN: 1402-1544
ISBN 978-91-7439-301-9 Luleå 2011
Thesis for the degree of Doctor of Philosophy (sv. Teknologie doktor) Research subject: Scientific Computing
Supervisors:
Associate professor Inge S¨oderkvist Senior lecturer Ove Edlund
The supervisors are active at: Mathematical Science,
Department of Engineering Sciences and Mathematics, Lule˚a University of Technology, Sweden.
Abstract
Computational methods for shape verification of free-form surfaces are the main content in this doctoral thesis. A common property of these methods is that they enable on-line shape verification directly in the production line. Hence, the considered methods must be fast and robust.
One of the problems in shape verification of free-form surfaces is registration. That is the problem of matching data points in 3D-space, representing the measured shape, to a CAD-object by applying a rigid body transformation. A method to accomplish the registration fast and robust is developed. It is a development of the iterative closest point (ICP) method. We are pre-processing the CAD-object by creating a data structure, which enables fast closest point search. Initially, much time is spent on creating the data structure for making each registration fast. The robust registration approach is based on theories from robust statistics by using iteratively re-weighted least squares in combination with the ICP method. It gives a fast registration method that is insensitive to deviating data.
The registration method is used in an application aimed for finding deviations be-tween the shape of an object and its ideal shape. The ideal shape is known and is given by a CAD-object. An optical shape measurement method, projected fringes with a single pattern recording, is used to acquire the data points. This method is fast and insensitive to vibrations but the data points might contain errors in some regions, which are handled in the registration.
An inverse problem that arises in many optical shape measurement methods is phase unwrapping. We are introducing a novel phase unwrapping method with regularization using shape information from a CAD-object, which is describing the shape of the mea-sured object. The shape measurement method that we are using here is based on dual wavelength holography. Our phase unwrapping method is insensitive for discontinuities but the measured object can not deviate too much from the shape of the CAD-object. A procedure for obtaining the required shape information fast from the CAD-object is also developed.
In order to acquire suitable shape information from data points, a parametric curve or a parametric surface, e.g. NURBS, can be fitted to these points. A subproblem appearing when fitting NURBS using the Gauss-Newton method is studied. Computa-tional aspects for obtaining a search direction are considered. We are also considering methods for NURBS fitting that are based on techniques for separable nonlinear least squares problems. These techniques use projections in order to separate the computation of the linear parameters from the computation of the nonlinear parameters.
Acknowledgment
First of all I would like to thank my supervisors, Associate Professor Inge S¨oderkvist and Senior Lecturer Ove Edlund. They have both been a great support during my Ph.D. studies. Inge and Ove have shared their scientific proficiency which was very instructive for me. Other closely related persons I would like to thank is Ph.D. student Sara Rosendahl, Senior Lecturer Per Gren and Professor Mikael Sj¨odahl, all three active at the Division of Fluid and Experimental Mechanics, Department of Engineering Sciences and Mathematics, Lule˚a University of Technology (LTU). We have been involved in the same research projects during my time as a Ph.D. student. Sara, Per Gren and Mikael have been working with optical shape measurement methods and have done all the experimental setups. I have used their obtained measurement data in processing.
Henrik Nerg˚ard, Division of Innovation and Design, Department of Business Admin-istration, Technology and Social Sciences, LTU, is also a person I would like to thank. He has created CAD-objects which I have used in the research. In addition, he has showed me some fundamentals about creation of CAD-objects using CAD-software which was instructive for me. Gestamp HardTech has shared files of CAD-objects representing real parts of the body of a car, the B-pillars. I have learned a lot about data specification in CAD by studying these files.
Research of considerable proportions is not possible without financial support. That is why I also would like to thank VINNOVA (The Swedish Governmental Agency for Innovation Systems) for the financial support. Two research projects, OPTIPRO (sv. optisk teknik f¨or verifiering av produktdata, en. optical measurement technology for verification of product data) and HOLOPRO (sv. holografisk teknik f¨or verifiering av produktdata, en. holographic product verification), have financed my Ph.D. studies. In OPTIPRO, lots of different companies have been involved and I would also like to thank them for their contribution.
In the end of this acknowledgement I want to thank Vetenskapsr˚adet (the Swedish Research Council) and the Swedish National Infrastructure for Computing (SNIC) for giving me the opportunity to attend the Swedish National Graduate School in Scien-tific Computing (NGSSC). The courses given by NGSSC have been located at different Swedish universities so I have got an insight in other universities than my home univer-sity. These courses have given me new knowledge and much experience.
List of Papers
This doctoral thesis is a compilation thesis constituent of seven papers, Papers I–VII. The thesis starts with an introductory part giving the context of the papers followed by a summary of them and a discussion about ideas for further improvements.
Paper I
Repeated Surface Registration for On-Line Use
Per Bergstr¨om, Ove Edlund, and Inge S¨oderkvist
International Journal of Advanced Manufacturing Technology, vol. 54, no. 5,
pp. 677–689, 2011
Paper II
Robust Surface Registration for Optical
Single-Shot Measurements
Per Bergstr¨om and Ove Edlund To be submitted
Paper III
Shape Verification Aimed for
Manufacturing Process Control
Per Bergstr¨om, Sara Rosendahl, and Mikael Sj¨odahl
Optics and Lasers in Engineering, vol. 49, no. 3, pp. 403–409, 2011
Paper IV
Perspective Depth Extraction of Points on a Surface under
an Instantaneous Rigid Body Transformation
Per Bergstr¨om
Yellow series, Department of Engineering Sciences and Mathematics,
Paper V
Shape
Verification
using
Dual-Wavelength
Holographic
In-terferometry
Per Bergstr¨om, Sara Rosendahl, Per Gren, and Mikael Sj¨odahl
Optical Engineering, vol. 50, no. 10, October 2011
Paper VI
Efficient Computation of the Gauss-Newton Direction when
Fitting NURBS Using ODR
Per Bergstr¨om, Ove Edlund, and Inge S¨oderkvist Submitted to: BIT
Paper VII
Fitting NURBS Using Separable Least Squares Techniques
Per Bergstr¨om and Inge S¨oderkvist Submitted to: BIT
Additional Scientifical Contribution
Dual-wavelength holographic shape measurement with
itera-tive phase unwrapping
Sara Rosendahl, Per Bergstr¨om, Per Gren, and Mikael Sj¨odahl
Optical Measurement Systems for Industrial Inspection VII, Proc. SPIE 8082, 80820B, 2011
Table of Contents
1. Introduction . . . 7
2. Applied Shape Verification . . . . 10
2.1. Motive . . . 10
2.2. A Brief Review . . . 10
2.3. Our Approach . . . 11
3. Optical Shape Measurement Methods . . . 12
3.1. Projected Fringes . . . 12
3.2. Holography . . . 14
4. Phase Unwrapping . . . . 16
4.1. Ideal Univariate Phase Unwrapping . . . 16
4.2. Applied Bivariate Phase Unwrapping . . . 18
5. Registration. . . . 19
5.1. Problem Description . . . 19
5.2. The ICP Algorithm . . . 19
5.2.1. Rigid Body Movement Problem . . . 21
5.2.2. Point-to-Plane Distance Minimization . . . 21
6. Robust Estimation . . . 24 6.1. Ordinary Regression . . . 24 6.2. Robust Regression . . . 25 6.3. IRLS . . . 27 7. NURBS . . . 29 7.1. Curves . . . 29 7.2. Surfaces . . . 31
7.3. Trimmed NURBS Surfaces . . . 33
8. Data Specification in CAD . . . 36
8.1. The IGES Data Format . . . 36
9. Model Surface Computations . . . 39 9.1. Closest-Point Calculations . . . 39 9.2. Line Intersections . . . 40 10. Unconstrained Optimization . . . . 42 10.1. Steepest Descent . . . 42 10.2. Newton’s Method . . . 43 10.3. Gauss-Newton Method . . . 44 10.4. Levenberg-Marquardt Method . . . 45 10.5. General Remarks . . . 46
11. Applied Curve and Surface Parameter Estimation . . . 47
11.1. Unconstrained Data Fitting . . . 47
11.2. Surface Reconstruction in Reverse Engineering . . . 48
12. Data Fitting Using Separable Nonlinear Least Squares . . . . . 50
12.1. Problem Description . . . 50
12.2. Separation of Variables Using Variable Projection . . . 50
12.3. Remarks . . . 52 References . . . . 53 Summary of Papers . . . 58 Paper I . . . 58 Paper II . . . 58 Paper III . . . 59 Paper IV . . . 59 Paper V . . . 60 Paper VI . . . 60 Paper VII . . . 61
Paper I
.Repeated surface registration for on-line use
Abstract . . . 67
1. Introduction . . . 67
2. Problem Definition and the ICP algorithm . . . 68
3. Finding the Closest Point . . . 71
3.1. Related Work . . . 71
3.2. The Distance Varying Grid Tree . . . 72
4. The Model Points and Linear Approximations . . . 76
5. Numerical Experiments . . . 78
5.1. Time Comparison . . . 82
5.2. Experimental Convergence . . . 83
6. Results and Conclusion . . . 87
7. Possible Improvements . . . 88
Acknowledgment . . . 88
References . . . 89
Paper II
. Robust Surface Registration for Optical Single-Shot Measurements Abstract . . . 951. Introduction . . . 95
2. Related Work . . . 96
3. Problem Formulation . . . 97
4. Conditions for Optimality and IRLS . . . 99
5. Barrier Points . . . 101
6. The ICP Algorithm with IRLS . . . 102
7. Robust Estimators and their Parameters . . . 103
8. Numerical Experiments . . . 105
9. Conclusion . . . 114
10. Discussion . . . 115
Acknowledgment . . . 116
Paper III
.Shape Verification Aimed for Manufacturing Process Control
Abstract . . . 123
1. Introduction . . . 123
2. Shape Measurement Method . . . 125
2.1. Stereo Camera System . . . 125
2.2. Projected Fringes . . . 126
2.3. Shape Determination . . . 128
3. Pre-Processing the Design Model . . . 130
4. Surface Matching . . . 130
4.1. The ICP Algorithm . . . 131
4.2. Robust Matching . . . 132
5. Experiment . . . 132
5.1. General Description . . . 132
5.2. Experimental Setup . . . 134
5.3. Optical Shape Measurement . . . 134
5.4. Surface Matching and Results . . . 135
6. Discussions and Conclusions . . . 136
Acknowledgment . . . 137
References . . . 138
Paper IV
. Perspective Depth Extraction of Points on a Surface under an Instantaneous Rigid Body Transformation Abstract . . . 1451. Introduction . . . 145
2. Line-Surface Intersection . . . 147
3. Problem Solving Approach . . . 147
3.1. Fundamental Concepts . . . 147
3.2. Pre-Processed Shape Information . . . 149
3.3. Explicit Formulation . . . 152 4. Two Algorithms . . . 153 4.1. Algorithm IV-1 . . . 153 4.2. Algorithm IV-2 . . . 153 4.3. General Comments . . . 156 5. Results . . . 156
6. Conclusions . . . 158
7. Discussion . . . 159
References . . . 160
Paper V
. Shape Verification using Dual-Wavelength Holographic Interferometry Abstract . . . 1651. Introduction . . . 165
2. Holographic Shape Measurement . . . 167
2.1. Dual Wavelength Holographic Recording . . . 167
2.2. Experimental Setup . . . 167
3. Geometry and Pre-Processing . . . 170
3.1. Perspective Geometry . . . 170
3.2. Pre-Processing the CAD-Model . . . 171
4. Iterative Phase Unwrapping . . . 172
4.1. Theory . . . 172
4.2. Results . . . 175
5. Discussions and Conclusions . . . 177
Acknowledgment . . . 178
References . . . 179
Paper VI
. Efficient Computation of the Gauss-Newton Direction when Fitting NURBS Using ODR Abstract . . . 185 1. Introduction . . . 185 2. Problem Description . . . 187 3. Splitting of Variables . . . 189 4. Matrix Properties . . . 192 5. Computational Experiments . . . 194 5.1. General Information . . . 1945.2. Accuracy of the Output . . . 196
5.3. Time Comparison . . . 199
6. Discussion . . . 200
7. Summary . . . 201
8. Appendix . . . 202
References . . . 204
Paper VII
. Fitting NURBS Using Separble Least Squares Techniques Abstract . . . 2091. Introduction . . . 209
2. The Variable Projection Method . . . 211
3. Variable Projection for NURBS-Fitting . . . 214
3.1. The One Dimensional NURBS-Fitting Problem . . . . 215
3.2. The NURBS-Fitting Problem of Arbitrary Dimension . 217 3.3. Computational Details and Comparison of Work Load . 218 4. Numerical Tests . . . 219
4.1. Generation of Problems . . . 220
4.2. Convergence . . . 221
4.3. Elapsed Real Time . . . 226
5. Summary . . . 227
1
. Introduction
All research presented in this doctoral thesis considers methods for solving different computational problems arising in metrology. Metrology is here referring to optical methods for measuring the shape of objects. To being able to solve our considered computational problems as efficiently as possible, we have to make the best use of the capacity of computers with respect to processor speed and available memory. The limited precision in floating point arithmetic is also considered.
The presented research is within the field of scientific computing. Numerics, applied mathematics, computational geometry, computational metrology, geometric modeling, optimization and computational statistics are some disciplines in scientific computing that can be used to classify the content of this thesis.
Shape verification is done in three stages; measurement, registration and comparison. The object, whose shape that is to be verified, has first to be measured. This measure-ment results in a set of data points representing the measured shape. The measured shape is to be compared to its ideal shape, represented by a CAD-object (Computer Aided Design). With the compounded word “CAD-object”, we refer to the digital ob-ject formed by a given free-form surface specified in a CAD system. Since the measured object is not perfectly fixtured, the obtained data points must be transformed under a rigid body transformation before a comparison with the CAD-object can be done. The procedure of finding this transformation is called registration. When the registration is done, the measured shape can be compared to the CAD-object. A schematic illustration of the process of shape verification is given in Figure 1. After this introduction, in
Chap-Measurement Registration Comparison
Figure 1. Schematic illustration of the process of shape verification
ter 2, we are giving the motive and a brief review of shape verification in manufacturing processes. We are also describing our approach to the problem and the purpose of our computational methods.
In the shape verification, the CAD-object of the measured object is used as an ideal shape. The CAD-object consists of a composite surface constituent of several parametric surface patches. Those parametric surfaces are usually represented by trimmed NURBS (Non-Uniform Rational B-Splines) surfaces. Depending of the complexity of the object, the number of trimmed NURBS surfaces is usually in the order of hundreds up to several thousands. In Chapter 7 we are giving a description of NURBS and in Chapter 8 we
are mention something about data specification in CAD, that is how the geometrical entities of the CAD-object are represented and stored.
The optical measurement methods that we are using are based on detecting intensity of light. From received light intensity, adequate shape information is computed. The processing of data in the optical measurement methods can be classified in different levels of abstraction. Physical and mathematical models are used to convert detected intensity of radiation in the light detector to useful shape information. The majority of the content in this thesis considers problems of processing given data points. Hence, a representation of the measured shape is given and the task is to obtain useful information from it.
This doctoral thesis is a compilation thesis constituent of seven papers and an intro-ductory part with twelve chapters. The introintro-ductory part presents the context of the research, and the research itself is presented in the papers. Three categories of research topics are presented. These are:
• phase unwrapping; • registration; • NURBS fitting.
The problem of phase unwrapping occurs in optical shape measurements, registration is already mentioned, and NURBS fitting can be used for processing data points and obtain useful shape information. A brief introductory description of the topics and our approach to them follows.
The phase unwrapping problem is an inverse problem. A wrapped phase function is given, and the original phase function is to be recovered. In Chapter 4 a more detailed description of phase unwrapping is found. We are presenting our considered phase unwrapping problem in Paper V.
In our considered problem we have a bivariate wrapped phase function with range [−π, π]. The unknown and unwrapped original phase function corresponds to the depth of a measured surface. The surface of the measured object has discontinuities, which imply discontinuities in the unwrapped phase function that is to be recovered. Such discontinuities results in ambiguities in the phase unwrapping, and to handle that, we are using shape information from the CAD-object of the measured object. Hence, our inverse problem of phase unwrapping is solved with regularization.
The used shape measurement method is optical, based on digital holography using two different wavelengths of coherent light. See Section 3.2 for some further details about this technique. Interference and detection of light intensity gives a two-dimensional digitalized wrapped phase function, called a phase map.
For each pixel in the obtained phase map, we have to acquire a corresponding point on the object. That is done for finding the phase values on the surface of the CAD-object, which is used in the phase unwrapping. In order to obtain the shape information
fast, we are using pre-computed sampled shape information where the CAD-object is in a fixed position. A description of this method is found in Paper IV.
Another handled research topic considers registration, see Papers I–III. Registration is the problem of finding a rigid body transformation, a rotation and translation, of a set of points so their orientation fit a given shape. This problem is solved using the ICP (Iterative Closest Point) method, see Chapter 5. It is an iterative method that solves two problems in each iteration, a closest point problem and a rigid body estimation problem.
We have developed the ICP method so it makes the registration process fast and robust. It is fast because we are using pre-processed shape information, see Chapter 9 and Paper I, and the robustness is developed by applying methods from robust statistics, see Chapter 6 and Paper II. The developed method is used in an application of shape verification, where the data points are obtained from an optical measurement method using projected fringes, see Section 3.1 and Paper III.
The fast and robust registration method is developed aimed for being suitable to perform shape verification on-line in the production line. The condition of on-line per-formance in the production line implies that the method must return the output very quickly and it must be robust for strongly deviating data. If the shape verification takes too long time, it will be a bottleneck and constitute an obstacle for the production. If the shape verification is sensitive for errors, the output will not be useful in the production either.
The last presented research topic considers the problem of fitting NURBS to a set of data points. Two problems are considered. In one approach we are studying how a search direction can be efficiently computed when fitting NURBS in ODR (Orthogonal Distance Regression) sense using the Gauss-Newton algorithm, see Paper VI. The other problem, considered in Paper VII, is about separation between the linear and the non-linear parameters, which are to be estimated. In Chapter 11 we are giving some details about curve and surface fitting, and in Chapter 12 we describe the idea behind variable separation. We are also describing some basic optimization methods in Chapter 10, which is used in the NURBS fitting papers.
2
. Applied Shape Verification
The main content in thesis is about computational methods for shape verification of free-form surfaces. We will focus on the computational methods but here we will mention the fundamental motive of shape verification and give a short review of related work. We will also describe our approach to the problem.
2.1.
Motive
In the early part of the industrial revolution all the manufacturing was performed by workers who were operating the machinery. The process was controlled by humans in a great extent since the machinery devices were simple. Assembling was also done by humans. The manufactured parts were put together in a handicraft trade. The shapes of a product from mass production differed a lot, but at that time they were good enough. As time is getting on, the requirements of manufactured parts have changed. Prod-ucts must be produced at low costs and have a high quality. In the automotive industry, the requirements are very high. Nowadays the assembling of cars is performed by robots in great extent. The robots are not like humans, they are just doing exactly what they are told to do. If a part has some defect in the shape, the robot is just mount it as it was correct with the consequence that the resulting worked up product has a defect. Discarding products late in the production line is expensive and should be avoided. An-other example where shape verification is of great importance is in machine production. A typical example is production of turbine blades. Their shape can not deviate too much from an ideal shape. They should rotate at very high speed so very small deformations can affect their function significantly.
2.2.
A Brief Review
The need of shape verification is substantial in industrial manufacturing processes. Dif-ferent applications need different solutions. Researches in the area have been presented in many different types of journals considering e.g. manufacturing, computer vision, computer aided design, robotics and computer integrated engineering. Li and Gu have written a review article about free-form surface inspection, [1]. They distinguish the measurement methods between non-contact and contact. The non-contact methods use optics with electromagnetic radiation, and the contact methods use a moving mechanical tactile probe. The method that is suitable for an application is dependent of conditions and requirements.
The contact methods, performed by a coordinate measuring machine (CMM), takes longer time than the non-contact methods, but all surfaces can be measured independent of optical properties. A CMM measures 3D-coordinates with high precision using a tac-tile probe, at one point in the moment. The path of the probe must be careful prepared since it takes long time to use, [2]. The probe might also be bulky for some applications where narrowed regions are to be measured. It needs some space to remove. In next chapter, we will describe two non-contact methods. These optical shape measurement methods measure a surface in full-field, that is a whole area of the surface is measured at once.
In [3] a process for inspection process of turbine blades is described. A CMM is used to perform the measurement. The inspected part is closely attached to a fixture. After measurement, a registration to align the data points to the CAD-object must be done. Zhang et al. [4] describes a method for surface inspection of windshields. They are using an optical measurement system to measure the shape. Also they have to align the data points to the CAD-object.
In all the referenced methods, the surface is first measured and then, a workpiece localization problem of registration (see [5, 6] and Chapter 5) has to be solved. A rigorous fixturing of the object that is to be measured takes time, which is the reason for using software to localize the measured object. Since the measured object is not perfectly fixtured within the accuracy of the measurement, the data points obtained from the measurement must first be aligned to the CAD-object before a comparison can be done and the measured shape can be verified. Some workpieces must be at least partially fixtured since they are not completely rigid, making shape verification a bit more ambiguous.
2.3.
Our Approach
In contrast to many other approaches of shape verification, where the shape is measured and compared offline in industrial processes apart from the production line, we are developing methods for performing the shape verification online. That is, the methods we are describing in Papers I–V solves arisen problems under the conditions that the shape verification is being adapted for on-line functionality in the production line. Hence, the methods must be fast and robust so they can be used in fully automated processes. On-line performance makes it possible to detect defective objects without any unnecessary delays. During delays several other objects might be produced with the same defect.
The measurement methods we are using are optical, see next chapter. We are per-forming the measurement in a single full-field shot, which enables fast measurements. Hence, the measurements are less sensitive for vibrations. In industrial production en-vironments, vibrations are an apparent feature, and are usually disturbing shape mea-surements of high quality.
3
. Optical Shape Measurement
Methods
Optical shape measurement methods are techniques based on illumination and detection of light. The object to be measured is illuminated by some type of light and the intensity of the reflected light is detected. From the detected intensity, a shape representation of the measured object is acquired. Here we are giving two examples of optical shape measurement methods.
3.1.
Projected Fringes
Projected fringes [7–10] is an optical structured-light shape measurement method. We are using this method in Paper III. It is based on projection of a fringe pattern from one view onto an object whose shape that is to be measured. Commonly, this fringe pattern has a sinusoidal intensity but other patterns are also possible. The projected fringe pattern is captured from two cameras, on opposite sides of the projector, in a stereo-camera view. Triangulation and phase unwrapping are used to acquire a digital surface representation of the measured surface. This surface representation is a set of data points whose shape forms a surface.
Figure 2 illustrates a schematic set-up of the projected fringes method. A relatively large distance B, result in high accuracy in the triangulation, but low accuracy due to fringe pattern interpretation. When B is large, parts of the object might be concealed or be in the shadow of the illumination.
In Figure 3 we have used two cameras to capture fringe images of reflected light. The projected pattern of fringes has a sinusoidal intensity with fringe period of approximately 1 mm. A large fringe period results in high accuracy in the fringe pattern interpretation but low precision due to insufficient quality of detected intensity. The width of the captured fringes tells about the slope of the surface perpendicular to the direction of the fringe projection. The fringe pattern, viewed from the stereo-camera device, tells about the shape.
Regions that are in the shadow of the projections or concealed for the cameras will not be measured. That result in errors and absent data in the representation of the measured surface. When the measurements are done in an industrial production line, the industrial environment might cause errors due to disturbing reflections of light and dirt and dust on the optical system.
Lens Camera Camera Projector Lens Lens θ L Lp B Measurement volume Measurement volume
Figure 2. A schematic setup for a measurement system using projected fringes (Figure III-5)
Image creator: Sara Rosendahl
(a) (b)
Figure 3. Captured fringe images from (a) the left camera and (b) the right camera. The fringe period is approximately 1 mm in both images (Figure III-2(a),(b))
3.2.
Holography
Holography [11–14] is another optical method that can be used for shape measurements. This measurement technique is used in Paper V. Dual-wavelength holography is based on optical interferometry by illuminating an object, whose shape that is to be measured, with coherent light of two different wavelengths,λ1 andλ2. An example of a schematic
set-up for a dual-wavelength holographic shape measurement is given in Figure 4. The light beams from the laser sources are split by beam splitters to enable correlation between the object light and the reference light. The beams in the object light have different optical path lengths dependent of the shape of the object and interference gives information about the shape.
Dual-wavelength holographic shape measurement results in a phase map with a fringe pattern from detected light intensity at a camera, originating from interference of light. The phase map, with phase values on the interval [−π, π], is used for the acquirement of a resulting digital representation of the measured surface. The phase values corre-spond to a wrapped depth of the measured surface, with depth in the geometry of the measurement. The fringe pattern in the phase map appears as contour lines of the measured surface. One fringe period, with a phase difference of 2π, corresponds to a depth difference,δ, given by (V-3). The depth difference is proportional to a synthetic
wavelength, (V-4). Small values ofδ are obtained if the difference |λ1− λ2| is large and
the wavelengths are small. For a fixed depth difference of the measured surface, small values ofδ will give more dense fringes. That gives a high precision in the measurement but low accuracy due to ambiguities in fringe pattern interpretation.
In Figure 5, an example of an obtained interference fringe pattern from a holographic measurement is shown. Each fringe period correspond to a depth difference of 0.20 mm. In the measurement, a perspective projection geometry is used to describe the light beam propagation, Figure V-4, because of divergent light of illumination. If a flat surface is measured that is parallel to the detector, the fringe pattern will be circular because of differences in optical path lengths caused by the divergent light. That explains something of the pattern in Figure 5. The shown pattern corresponds to a measured surface of a milled trace on a cylinder, where the depth of the trace is 1 mm.
Unlike the method of projected fringes, we do not have the same troubles with shadowed regions, since illumination direction and view point direction is the same. However, the problem with dirt and dust on the optical system remains. This problem is common for all optical methods and cause errors. Irregularities in the surface of the measured object cause scattered reflections of incoming light. When using coherent light, these scattered reflections give rise to an interference pattern (speckles, see e.g. [10]), which has an undesired influence of the precision in the measurement.
In the both presented optical shape measurement methods, projected fringes and holography, a fringe pattern has to be interpreted in order to obtain a representation of the measured shape. A phase unwrapping problem arises in this pattern interpretation.
OD M1 M2 L Holographic setup Reference light Object light Laser diode HeNe-laser Camera Object Detector Michelson interferometer
Figure 4. Experimental setup of a dual-wavelength holographic measurement system and Michelson interferometer (Figure V-1)
Image creator: Sara Rosendahl
u
v
Figure 5. Measured phase map. Each fringe correspond to a depth difference
4
. Phase Unwrapping
In this chapter we are describing the problem of phase unwrapping. It occurs in various types of fields, e.g. signal processing. Connected to shape measurement, projected fringes analysis and optical interferometry are areas where phase unwrapping occur.
4.1.
Ideal Univariate Phase Unwrapping
Phase unwrapping is the inverse problem of recovering the original phase functionϕ(u), from a given wrapped phase functionϕw(u) with range (−π, π] or [−π, π). We are using
the notationW(ϕ) = ϕw, to denote the wrapping operator, which is
W(ϕ) = atan2(sin (ϕ), cos (ϕ)), W(ϕ) ∈ (−π, π] , (1) where atan2 is the four-quadrant arctangent operator. Another similar form of (1) is
W(ϕ) = {(ϕ + π) modulo 2π} − π, W(ϕ) ∈ [−π, π) . (2) The outputs of W in (1) and (2) are equal except where ϕ = π + n 2π, n ∈ Z. At these values of ϕ, W(ϕ) differ 2π between (1) and (2), giving an output of π and −π respectively. The task of phase unwrapping is to find a multiplen of 2π for each value ofu, so when adding n 2π to ϕw(u), the original phase function ϕ(u) is obtained. That
is, we have to find the integern = n(u), such that
ϕ(u) = n 2π + ϕw(u) . (3)
If we assume that the original phase functionϕ has jump discontinuities of magnitude less thanπ, we can recover the original phase function from ϕwusing phase unwrapping
without any ambiguities. This is performed by integrating the phase differences over
u, and when the magnitude of the phase difference is greater than π, 2π is added or
subtracted to make the magnitude of the jump less than π. An initial value of the unwrapped phase is necessary. This can be the value of the wrapped phase function at the start. However, the unwrapped phase function might still differ with a multiple of 2π from the original phase function, but this difference is independent of u. The phase values ofϕ is lost when it has jump discontinuities of magnitude greater or equal to π. In Figure 6, we have plotted three univariate phase functions; the original one, the wrapped one, and the unwrapped one. The unwrapped phase function is recovered as previously described, by starting at u = 0 and assigning the initial value of the unwrapped phase function to ϕw(0). We can see that the unwrapped phase function
u ϕ 2π 3π 4π 5π 6π π 0 -π Original phase Wrapped phase Unwrapped phase
Figure 6. Example of phase functions and successful unwrapping
u ϕ 2π 3π 4π 5π 6π π 0 -π Original phase Wrapped phase Unwrapped phase
Figure 7. Example of phase functions and unsuccessful unwrapping
and the original phase function are the same, so the phase unwrapping succeed in this example.
Ambiguities occur in the phase unwrapping whenϕ has jump discontinuities of mag-nitude greater or equal to π, which follows from (2). If we only know ϕw, we cannot
recover the original one, since important information is lost. This is illustrated in Fig-ure 7, where the same method of unwrapping was used. It is clear that the original phase function is not the same as the unwrapped one, so the phase unwrapping failed here.
4.2.
Applied Bivariate Phase Unwrapping
In the previous section we saw two examples of phase unwrapping of univariate phase functions. One of the phase functions was continuous while the other had jump discon-tinuities of magnitude greater thanπ. The jump discontinuities resulted in ambiguities, so the phase unwrapping did not succeeded. In fringe pattern analysis and optical in-terferometry, the wrapped phase function is bivariate and digitalized. It is represented by a discrete set of pixels with phase values, called phase map. Because of undesired circumstances, the phase values in the phase map have a noise. The noise makes the unwrapping procedure much more complicated than in previous section, see e.g. [15, 16]. In some regions in the phase map, the noise can be dominating so the phase values are completely blurred.
Two types of phase unwrapping techniques for bivariate phase functions are:
• Path-following methods; • Minimum-norm methods.
Both methods use the phase gradient, which is defined to be the wrapped phase difference of adjacent pixels in the phase map.
Path-following methods are methods that perform the unwrapping of the phase map by integrating the phase gradients over a path in the phase map. These methods are similar to the univariate phase unwrapping methods since the unwrapping is performed along a curve. Because of noise and discontinuities, the unwrapped values can have different values dependent of path of integration. The path-following methods are fast but not robust.
Minimum-norm methods are phase unwrapping methods that perform the unwrap-ping in terms of minimization of a sum of norms of the phase gradients. They do not rely on a path of integration. Methods using the minimum-norm technique are known to be robust but also more computational intensive. The 2-norm minimization method is the fastest one, since there exists a direct solution to the considered minimization problem.
In Paper V we are proposing a method to perform a bivariate phase unwrapping with regularization. We are assuming that the measured shape is similar to the shape of the CAD-object, and using this shape information in the phase unwrapping. It makes it possible to unwrap phases with discontinuities. Unfortunately, the measured shape can not deviate too much from the shape of the CAD-object. Our method is a type of minimum-norm method, since we are minimizing an objective function of a sum of norms.
5
. Registration
One of the main computational problems in shape verification is registration. That is the problem of finding a rigid body transformation, so the transformed data points fit the CAD-object in some error metric. In many applications, the term “best-fit” is used for the output of the registration. We are considering registration in Papers I–III and Paper V.
5.1.
Problem Description
The surface representation obtained from a shape measurement consists of a set ofNP data points, P = {pk}, and the surface of the CAD-object is represented by a set of
model points,X. We will give more details about the model points in Chapter 9. The problem of registration is to find a rigid body transformation of the data points, so that they fit to the model points inX, see Figure 8. This transformation consists of a rotation R∈ Ω, and a translation t ∈ R3, where
Ω ={R ∈ R3×3|RTR= I3, det (R) = +1}.
The transformation is found by solving min R∈Ω,t NP k=1 d(Rpk+ t, X)2, (4) where d(p, X) = min x∈Xx − p2.
That is, the operator d(p, X) is the minimum distance between an arbitrary point, p,
and the model points inX.
5.2.
The ICP Algorithm
A common way to solve problem (4) is to use the ICP algorithm, see e.g. [17–20], developed in the early nineties. A simple form of the ICP algorithm is presented in Al-gorithm 1. Here is the identity transformation used as initial transformation. In applied surface registration, there might exist better choices of initial transformation. The initial transformation should be chosen such that the expected transformation obtained from
Model Data points
Figure 8. Data points and model. The arrows indicate the rigid body transformation to be determined (Figure III-3)
the ICP algorithm is as small as possible. The ICP-algorithm is re-iteratively solving a closest point problem and a rigid body movement problem.
An operator C is used to denote the closest point search in Algorithm 1. For an arbitrary point, p, the closest point
y=C(p, X) ∈ X, (5)
satisfies
y − p2= d(p, X).
The closest-point search and the handling of the data points is the most time consuming part in the ICP algorithm. It must be done in each iteration for all used data points. In Paper I, we are considering a method to perform a closest point search, appropriate for usage in the ICP algorithm.
Algorithm 1ICP Require: X, P R= I3, t = 0 repeat yk=C(Rpk+ t, X), k = 1, . . . , NP [R, t] = arg min R∗∈Ω,t∗ NP k=1 R∗p k+ t∗− yk22 untilconvergence return R, t
5.2.1. Rigid Body Movement Problem
Once the closest model points to the data points are found in the ICP algorithm, the minimization problem min R∗∈Ω,t∗ NP k=1 R∗p k+ t∗− yk 2 2, (6)
has to be solved. Problem (6) is a rigid body movement problem, where a rotation matrix and a translation vector have to be found. For properties of this problem see e.g. [21–25]. The solution is computed from a singular value decomposition (SVD) of a cross-covariance matrix. There is also a solution based on quaternions, but we will not use that approach.
For the solution, we need the centroids of the two sets of points Y = {yk} and P , which are μY = 1 NP NP k=1 yk, μP = 1 NP NP k=1 pk, (7)
respectively. Using these centroids, the cross-covariance matrix, ΣY,P, of the two sets of points,Y and P , is ΣY,P = NP k=1 (yk− μY) (pk− μP)T= NP k=1 ykpTk − NPμYμTP. (8) Let the SVD of ΣY,P be
ΣY,P = UΣVT. (9) Given this SVD, the optimal rotation matrix R and translation vector t to the problem (6) are R= U ⎡ ⎣ 10 01 00 0 0 det(UVT) ⎤ ⎦ VT , t=μY − RμP. (10) In most common situations, the determinant satisfies det(UVT) = 1, so R = UVT.
5.2.2. Point-to-Plane Distance Minimization
An alternative to the point ICP algorithm, Algorithm 1, is to use a point-to-plane ICP algorithm, see [26, 27], stated in Algorithm 2. Instead of solving a rigid body movement problem (6), the problem
min R∗∈Ω,t∗ NP k=1 nTk (R∗p˜k+ t∗− yk) 2, (11) is solved. Here are the surface normals, nk, n2 = 1, at the model points, yk, also
Algorithm 2ICP, point-to-plane distance minimization Require: X, P R= I3, t = 0 repeat ˜ pk = Rpk+ t, [yk, nk] =C(˜pk, X), k = 1, . . . , NP [R0, t0] = arg min R∗∈Ω,t∗ NP k=1 nTk (R∗p˜k+ t∗− yk) 2 R= R0R, t = R0t+ t0 untilconvergence return R, t
distances. Unlike (6), it does not exist any simple solution to (11). Fortunately, there is an approximate solution to (11), which is simply obtained.
Consider the rotation matrix
R∗= ⎡
⎣ cos (ω1) cos (ω2)cos (ω1) sin (ω2) − cos (ω0) sin (ω2)+sin (ω0) sin (ω1) cos (ω2)cos (ω0) cos (ω2)+sin (ω0) sin (ω1) sin (ω2) − sin (ω0) cos (ω2)+cos (ω0) sin (ω1) sin (ω2)sin (ω0) sin (ω2)+cos (ω0) sin (ω1) cos (ω2)
− sin (ω1) sin (ω0) cos (ω1) cos (ω0) cos (ω1)
⎤ ⎦ ,
with anglesω0,ω1, andω2. For small angles, we can do the approximation
nT kR∗p˜k≈ nTk ⎡ ⎣ ω12 −ω12 −ωω10 −ω1 ω0 1 ⎤ ⎦ ˜pk = nT kp˜k+ (˜pk× nk) Tω, where ω = ω0 ω1 ω2 T
, and × denote the vector cross product. Using this
ap-proximation, problem (11) can be written as min ω,t∗ NP k=1 (˜pk× nk)Tω + nT kt∗− nTk(yk− ˜pk) 2 , (12)
which itself can be written as a linear least-squares problem min
θ Aθ − b
2
where θ = ω t∗ , A= ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ (˜p1× n1) T nT1 (˜p2× n2) T nT 2 .. . ... (˜pNP× nNP)T nT NP ⎤ ⎥ ⎥ ⎥ ⎥ ⎦∈ R NP×6, b= ⎡ ⎢ ⎢ ⎢ ⎣ nT 1 (y1− ˜p1) nT 2 (y2− ˜p2) .. . nT NP(yNP − ˜pNP) ⎤ ⎥ ⎥ ⎥ ⎦∈ RNP.
Problem (13) is a linear least squares problem and can be solved using e.g. the normal equations or QR-factorization.
For point sets that are quite well aligned, Algorithm 2 has usually faster convergence than Algorithm 1, see e.g. [28]. Initially, Algorithm 1 can be used but after some iterations, the registration method can be switched to Algorithm 2. An advantage with Algorithm 1 is that it is proven to be globally and monotonically convergent to a local minimum, [17].
We have not been considering the point-to-plane distance minimization in any of the papers, Papers I–VII, where other problems are analyzed. This subsection is a supplement to the main content, but point-to-plane distance minimization is important for surface registration in general. It constitutes a good alternative in real applications of registration, where the normals of the surface at the model points are known.
6
. Robust Estimation
The aim of robust estimation is to reduce the effect of deviating data values when fitting the data to a model. In Paper II, a robust ICP-method is used to reduce the influence of deviating data points in surface registration. The ideas originate from robust statistics, see e.g. [29–31].
In Sections 6.1 and 6.2, examples of linear regression are presented for better under-standing and justification of the concept “robust estimation”. In Section 6.3 we present a brief description of the IRLS method, which is used in Paper II.
6.1.
Ordinary Regression
Here follows an example of ordinary linear regression. Let the independent variable be x, and let the dependent variable be y. Further, let β0 and β1 be two parameters
for modeling the linear relation y =β0+β1x+, where the residual is a stochastic
variable.
The task in regression is to estimate values of β0 and β1 from given observations.
LetxiandyibeN observations of x and y respectively. A linear model for the observed values is
yi=b0+b1xi+ei, i = 1, . . . , N , (14)
and the problem is to find the estimates ˆb0 and ˆb1 of β0 and β1 respectively. If
originates from a normal distribution and ei =yi− (b0+b1xi) are observations from
this distribution, the values according to ˆ b0, ˆb1 = arg min b0,b1 N i=1 e2 i, (15)
will give the best estimation ofβ0 and β1 with respect to their expected values.
Ob-servations that are strongly deviating from the linear model in (14) will have a high influence of the parameter estimation. An example of linear regression, where a sum of squared residuals is minimized and the set of data points, (xi, yi), contains an outlier, is plotted in Figure 9. It is clear from the figure that the linear model does not fit the set of data points very well. The outlier has a high leverage in the estimation.
y
x
Figure 9. Example of linear regression by minimizing a sum of squared residuals
6.2.
Robust Regression
It is possible to overcome the problem of strong influence from deviating data points with the use of robust estimation methods. Instead of minimizing a sum of squared resid-uals (15), a functionρ is introduced and the estimated parameter values are obtained according to ˆb 0, ˆb1 = arg min b0,b1 N i=1 ρ(ei). (16)
Aρ-function should be continuous, symmetric, non-decreasing for positive values, and
not growing as fast as the square functions do for values of great magnitude. The last property makes strongly deviating values not so dominant in the estimation (16). If a non-convex ρ-function is used, new local minima are introduced, but the influence of strongly deviating data values are effectually reduced. No additional local minima are introduced for the convexρ-functions, which is an advantage in numerical computations. In statistic literature, some commonfunctions are used. Huber’s and Tukey’s ρ-functions are two of them. These are defined and plotted in Figure 10. The Huber
ρ-function is convex and grows quadratically for small residuals and it grows linearly for
larger residuals. Tukey’sρ-function is non-convex and grows approximately quadrati-cally for small residuals and it is constant for large residuals. An interpretation of the constant behavior of the ρ-function is that data with residuals in the constant
inter-val are rejected. The Huber’sρ-function and the Tukey’s ρ-function are both useful in different ways. Tukey’s ρ-function is less sensitive for outliers but Huber’s ρ-function has other advantages. The behavior of theseρ-functions is appropriate if the residual
ρ
r
0
Huber Tukeys Biweight
Huber (kHu> 0) Tukey’s Bi-weight (kT u> 0)
ρ(r) = r2/2, |r| ≤ kHu kHu|r| − k2Hu/2, |r| > kHu ρ(r) = (k2T u/6)[1 − {1 − (r/kT u)2}3], |r| ≤ kT u k2T u/6, |r| > kT u
Figure 10. Twoρ-functions. Huber is convex, while Biweight is not (Figure II-1)
y
x
Least−Square Huber
Figure 11. Regression examples: Linear robust regression using the Huberρ-function, solid line. Linear least squares regression, dotted line
is approximately normally distributed, but some of the residuals are strongly deviating. The choice ofρ-function depends on the application and the expected residuals.
If the estimated parameters are obtained from (16) instead of (15) using a robust
ρ-function, the outlier in Figure 9 will not have so much influence. The linear model
(14) for the same set of data points, (xi, yi), is shown in Figure 11 when the Huber
ρ-function is used and parameters are estimated according to (16). We can see that the
linear model obtained from the robust method fits the observed points, (xi, yi), except the outlier, much better than the linear model obtained from least squares regression. The robust parameter estimation gives in this case one distinct point that is strongly deviating from the others.
Robust estimation is used in data analysis to fit a model to the bulk of the data points. It uses a majority of the data points that fits the model well and leaves the deviating ones out. Robust estimation can be used to detect outliers in a set of data points. If minimization of least squares is used, the deviating points will have a significant influence of the estimation, which results in that many data points appears as outliers. Hence, the detection of outlier is more ambiguous than if a robust method were used. That is shown in Figure 9, where many data points are deviating from the estimated linear model. It is more tricky to conclude which ones of them that are outliers, than using the linear model shown in Figure 11. In multi dimensional parameter estimation, detection of outliers is even more complicated.
6.3.
IRLS
The method of iteratively re-weighted least squares (IRLS) is widely used for solving the minimization problem in (16). In each iteration, it finds the optimum of a weighted least squares problem, which is computational easier to obtain than the optimum of the original problem. Let a weight functionw(r) be defined as
w(r) =
1
rρ(r) if r = 0,
ρ(r) ifr = 0,
(17) where and denote the first and second derivative respectively. The w-function is derived from the conditions of optimality of the minimization problem in (16). The IRLS algorithm is stated in Algorithm 3 and solves a weighted least squares problem with weights from (17). The initial values of ˆb0 and ˆb1are required and these values are
obtained using for example least squares estimation, (15). In [32] it is shown that the IRLS method converges to the optimum of the minimization problem in (16) for convex
Algorithm 3IRLS
Require: xi, yi, i = 1, . . . , N and initial values of ˆb0 and ˆb1
repeat ei =yi− ˆb0− ˆb1xi, i = 1, . . . , N ˆb0, ˆb1= arg min b0,b1 N i=1 w(ei) yi− b0− b1xi 2 untilconvergence return ˆb0, ˆb1
There are also other methods than the IRLS method that can be used to solve problems like (16). The Newton’s method, see for example Chapter 10, is another method that can be used for finding the solution. Which method to chose is dependent of application and requirements. For some types of problems, the IRLS method is a good choice, for some other types of problems, another method is better.
In Paper II, the problem of registration is solved with a robust method based on IRLS. A transformation of a measured surface containing errors is estimated, such that it fits a model surface. The registration-derived weighted least squares problem (II-13), is a weighted least-squares rigid-body movement problem (see e.g. Section 5.2.1 for the non-weighted case). There exists a direct method for solving (II-13), which makes it possible to combine the IRLS method with the ICP algorithm for obtaining a robust method with (almost) no loss of computational speed. This method, the IRLS-ICP method, is presented in Algorithm II-1.
7
. NURBS
NURBS is an abbreviation for non-uniform rational B-splines. Here follows a brief summary of its concept with some properties included. A careful description of NURBS is found in e.g. [33–35]. We are giving a description of NURBS curves, NURBS surfaces and trimmed NURBS surfaces, which are all frequently occurring in CAD applications. The dimensiond in the considered space, Rd, is usually 2 or 3. In all our research we
are using NURBS.
7.1.
Curves
A NURBS curve, c = c(u), of degree φ is represented by
c(u) = n i=0Ni,φ (u)wipi n i=0Ni,φ (u)wi , 0≤ u ≤ 1, (18)
where{Ni,φ} is a set of B-spline basis functions, which are non-negative and dependent of u, and {pi ∈ Rd} is a set of control points forming a control polygon, and {wi} is
a set of weights corresponding to the control points. If all weights,wi, are equal, the denominator in (18) is constant and the curve c is “just” a B-spline curve. The basis functions,Ni,φ, with degreeφ, are defined over a knot vector
U = {0, . . . , 0, μφ+1, μφ+2, . . . , μk, μk+1, . . . , μn, 1, . . . 1}, (19)
where the knots are ordered nondecreasingly, i.e. μk ≤ μk+1. We have μk = 0, k ≤ φ andμk = 1, k > n. With n + 1 control points and a φth-degree NURBS curve, the knot vectorU has n + φ + 2 elements. The knots do not have to be equally spaced on [0, 1], wherefrom the word non-uniform originates.
In Figure 12 there is an example of a NURBS curve. The control points of the curve are also drawn connected with dashed lines forming a control polygon.
For simplicity we write the knots and u-values on the interval [0, 1]. The B-spline basis functions,Ni,φ, are invariant under an affine transformation,αμk+β, of the knots
and the u-parameter space. It is just the relative differences between the knots that affect the values of the basis functions.
The control points enter linearly to the NURBS curve representation and they keep the curve in place inRd. An affine transformation applied on all the control points will
Figure 12. Example of NURBS curve with control points
do the same transformation on the curve. Hence, such a transformation on the curve is computational easy to perform.
We can see from the definition (18) that a rescaling of the weights,αw1, α = 0, has no
effect of the curve. It is just their ratios that make sense. The weights have the property of distributing influence of the control points. A relatively large value of a weight pulls the curve towards the corresponding control point. Weights,wi, with different signs are allowed but usually it is avoided because it might give zero in the denominator. With positive weights, the curve c(u), is always in the convex hull spanned by the control points. That property is lost when using negative weights.
Some important properties of{Ni,φ(u)} (see [33], Ch. 2) are:
1. Ni,φ(u) = 0 if u /∈ [μi, μi+φ+1) (local support property);
2. in any given knot span, [μi∗, μi∗+1), at mostφ + 1 of the Ni,φ are nonzero, namely
the functionsNi∗−φ,φ, . . . , Ni∗,φ; 3. Ni,φ(u) ≥ 0 for all i, φ and u; 4. ni=0Ni,φ(u) = 1 for all φ and u.
This is the same as saying that for all feasible values of u, at most φ + 1 consecutive basis functions are non-zero and none of them are negative. We also know that the sum of the basis functions is one.
For simplified computations, a convenient way to rewrite the NURBS form (18) is to write it on B-spline form
cw(u) = n i=0 Ni,φ(u)pwi , (20) where pw i = wipi wi , c[k] = 1 cw[d]c w[k], k = 0, 1, . . . , d − 1. (21)
The notation of a vector together with a pair of square brackets with an integer inside represents an element in the given vector. Three steps to compute a point on a B-spline curve at a fixed parameter value,u∗, are required [33]. These steps are:
1. find the knot span [μk, μk+1) in whichu∗ lies; 2. evaluate the nonzero basis functions;
3. multiply the values of the nonzero basis functions with the corresponding control points.
For NURBS computations we also have to do a division according to (21).
One of the greatest advantages of NURBS is their capability of precisely representing conic sections, i.e. circles, ellipses, parabolas and hyperbolas, as well as more general free-form curves. In CAD applications, the circles and ellipses together with general free-form curves have a fundamental role. B-splines, that is NURBS with constant weights, do not have that property. Allowing different values of the weights gives more freedom and flexibility of the shape. Some more data has to be stored and some more computations have to be done in applications of NURBS compared to B-splines, but the advantages dominate.
7.2.
Surfaces
NURBS surfaces have similar properties as NURBS curves. A NURBS surface, s =
s(u, v), is represented by s(u, v) = n i=0 m j=0
Ni,φ(u)Nj,ω(v)wi,jpi,j n i=0 m j=0
Ni,φ(u)Nj,ω(v)wi,j
, 0≤ u, v ≤ 1, (22)
where {pi,j ∈ R3} is a set of control points forming a bidirectional control net, and
{wi,j} is a set of weights corresponding to the control points. If all weights are equal,
the denominator is constant and the surface is a B-spline surface. The two sets of basis functions{Ni,φ(u)} and {Nj,ω(v)}, with degree φ and ω, are the B-spline basis functions defined on the knot vectors
U = {0, . . . , 0, μφ+1, μφ+2, . . . , μk, μk+1, . . . , μn, 1, . . . , 1},
V = {0, . . . , 0, νω+1, νω+2, . . . , νk, νk+1, . . . , νn, 1, . . . , 1},
(23) respectively. The knots in the knot vectors are ordered nondecreasingly, i.e. μk ≤ μk+1
andνk ≤ νk+1. An example of a NURBS surface with control points is given in Figure 13. The control points, lying in a bidirectional control net, are also drawn connected with dashed lines.
Figure 13. Example of NURBS surface with control points in a bidirectional control net
Many of the properties for NURBS surfaces are similar to the properties of NURBS curves, which is expected because of the similarity between (22) and (18). It is just the product basis functions that differ. A rescaling of the weights has no effect of the surface, it is just their internal ratios that make sense. With positive weights, the surface is in the convex hull spanned by the bidirectional control net. An isoparametric curve on a NURBS surface is a NURBS curve.
Some important properties of the product basis functions,Ni,φ(u)Nj,ω(v), (see [33], Ch. 3) are stated below.
1. Ni,φ(u)Nj,ω(v) = 0 if (u, v) /∈ [μi, μi+φ+1)× [νj, νj+ω+1);
2. In any given knot span rectangle, [μi∗, μi∗+1)×[νj∗, νj∗+ω+1), at most (φ+1)(ω+1)
of the basis functions are nonzero, in particular the Ni,φ(u)Nj,ω(v) for i∗− φ ≤
i ≤ i∗ andj∗− ω ≤ j ≤ j∗;
3. Nonnegativity: Ni,φ(u)Nj,ω(v) ≥ 0 for all i, j, φ, ω and (u, v);
4. Partition of unity: ni=0mj=0Ni,φ(u)Nj,ω(v) = 1 for all φ, ω and (u, v).
These properties are consequences of the properties for the single basis function,Ni,φ(u),
given previous. When the NURBS degrees are low and the number of control points is large, many of the product basis functions,Ni,φ(u)Nj,ω(v), are zero at a specific (u, v). That is utilized in NURBS surface computations.
The NURBS surface (22) can be written as a B-spline, sw(u, v) = n i=0 m j=0
Ni,φ(u)Nj,ω(v)pwi,j, (24)
where pw i,j= wi,jpi,j wi,j , s[k] = 1 sw[d]s w[k], k = 0, 1, . . . , d − 1, (25)
which is convenient for computations. Five steps to compute a point on a B-spline surface at fixed parameter values, (u∗, v∗), are required [33]. These steps are:
1. find the knot span [μk, μk+1) in whichu∗ lies;
2. compute the nonzero basis functionsNk−φ,φ(u∗), . . . , Nk,φ(u∗); 3. find the knot span [νk, νk+1) in whichv∗ lies;
4. compute the nonzero basis functionsNk−ω,ω(v∗), . . . , Nk,ω(v∗);
5. multiply the values of the nonzero basis functions with the corresponding control points.
For NURBS computations we also have to do a division according to (25).
NURBS surfaces provide a unified mathematical basis for representing both analytic shapes, such as surfaces of revolution as well as more general surfaces. Many important types of surfaces of revolution are exactly represented by NURBS. These surfaces are very useful in CAD, for example to represent spheres, circular cylinders and other rotational symmetric surfaces. NURBS can also be used to represent bilinear surfaces and ruled surfaces.
7.3.
Trimmed NURBS Surfaces
In CAD-applications it is important that the surface representations are able to fit an arbitrary shape. A single NURBS surface is not able to fit an arbitrary surface with sharp edges so well. However, a composite surface of several trimmed NURBS surfaces is suitable for doing that. That is the reason why trimmed NURBS surfaces are normally used in CAD.
The NURBS surface in (22), has a rectangular domain with sides parallel to the coordinate axis. That is not the case for trimmed NURBS surfaces [35–37], where the domain is just an arbitrarily connected region with no other restrictions. It does not have to be convex and it can contain holes with inner boundaries where the surface parameter values are not defined. The domain is bounded by an outer closed curve and any possible inner closed curves that may exist. An example of a domain of a NURBS
