Computer-Aided Design 33 (2001) 593-603 Parameterization for surface fitting in reverse engineering L.A. Piegla,*, W. Tillerb aDepartment of Computer Science and Engineering, University of South Florida, 4202 Fowler Avenue, ENG 118, Tampa, FL 33620, USA bGeomWare, Inc., 3035 Great Oak Circle, Tyler, TX 75703, USA
Abstract Given four boundary curves and a set of random points lying on a surface patch, a method for assigning parameters to these points is presented. The algorithm uses various base surfaces to project the points onto these surfaces to recover the parameters based on the surfaces’ underlying parameterization. Several techniques for speeding up the time consuming projection process are also presented. A thinning method is introduced as well to select a subset of the points that may have been sampled at a much higher rate than necessary. The thinning is based on the geometry of the base surface and relies on a meaningful geometric tolerance.
Keywords: Surface fitting; Reverse engineering; NURBS
1 Introduction Object reconstruction from sampled data forms an important part of obtaining the design parameters of an existing object. Graphical models [1—3,7—14] as well as engineering parts[15,27,28,30—33] must be faithfully reconstructed, that is, some form of representation is sought that defines an object that does not deviate from the points more than a given error measure. The representation can be discrete, e.g. a set of triangles, or continuous, e.g. a B-rep model. The focus of this paper is to obtain a B-rep model through surface fitting to subsets of points. Fitting surfaces to randomly distributed data has been the subject of significant interest [4—6, 10, 16—25, 29, 33]. Using B-splines, the critical issue of such fitting process is the computation of parameters and the knot vectors. This paper deals with the first issue. More precisely, the following problem is dealt with. Assume that the data has been split up into rectangular regions bounded by curves. That is, given four boundary NURBS curves (Fig. 1)
a set of randomly distributed points (it is assumed that the points on the boundaries are stripped off)
are sought that best represent the points’ position within the rectangular domain, i.e. the data points are assumed in the NURBS representation at these parameters. A viable solution to this problem is to start with a base surface, fitted somehow to the boundary curves and/or (some of) the points and then projecting the data points to this surface to obtain the parameters of the projected points. Theoretically this all sounds easy to do, however, there are many technical hurdles that can render this base surface based approach completely unusable.
The issue of parameterization for curve and surface fitting has been a difficult one since the advent of computers in design engineering. There is no known best parameterization, and even after many decades of research some data sets just cannot be dealt with adequately. This paper does not claim to solve the parameterization problem either. Especially not for random surface points. The aim is to give the designer a set of tools by which he can obtain good parameters at a reasonable cost. It is our opinion that the parameterization process cannot be completely automatic. The user may have to choose different options for different data sets, and he may even have to iterate a few times to obtain better results. The contribution of this paper is twofold: (1) the introduction of various base surfaces to capture more detail inherent in the data; and (2) a methodology to gain significant speed-ups during point projection. The paper is organized as follows. Section 2 deals with base surfaces, and Section 3 is devoted to obtaining stability and speed-ups during point projection. A Section 4 closes the paper.
2 Base surfaces The general idea of obtaining parameters is to use a surface attached to the input parameters [21]. In the sections that follow, we investigate three kinds of surfaces: (1) the bilinearly blended Coons patch; (2) the bicubic Coons patch; and (3) the tensor-product surface.
2.1 Bilinear Coons The simplest base surface is the so called bilinearly blended Coons surface. Using NURBS, the surface is expressed as follows [26]:
where is a linearly lofted surface (also known as a ruled surface) in the u-direction; is a linearly lofted surface in the v-direction; and is a tensor product surface, i.e. it is a bilinear surface defined by the four corner points. Once the two ruled surfaces and the bilinear surface are constructed, they are made compatible, and the control points of are obtained by adding the control points of the ruled surfaces and subtracting those of the bilinear patch. The advantages of using a bilinear Coons patch are: it is very easy to compute; the computation is stable and fast; and the surface is well parameterized. If the data set is simple, and in many cases it is (a good reverse engineering package includes a segmentation module that subdivides the data set into very simple patches), then the bilinear Coons patch is the best base surface we can think of. However, one disadvantage with this patch is that it tends to miss internal detail, as illustrated in Fig. 2. If the segmentation process failed to do a good job, then the bilinear Coons may provide inadequate parameters.
2.2 Bicubic Coons The idea behind using a bicubic patch is to include some internal detail into the base surface. The bicubic Coons patch has the same form as the bilinear except that the lofted surfaces are not linearly but cubically blended. To do so, we need cross derivatives. The computation of such a surface is quite involved, therefore only the main ideas are outlined below. Assume now that in addition to the boundary curves some number of internal points is also selected. Fixing the number of control points and the degrees in both directions, the unknown entities are the cross derivatives
where, for example, is a cross derivative field across the boundary , i.e. it is a function of u and its vectors, evaluated at various parameters, are the derivatives pointing across the boundary curve. To compute these cross derivatives, the unknown control points
are sought. For simplicity, the surface ; at fixed parameters, are expressed using cubic Hermite polynomials
defined over the domain [0, 1]. If the surface is defined over a general domain of , then the parameter transformations
bring the domains to the unit span. Note that
Now, for a fixed parameter pair ; the bicubic surface is written as
where and are the cubic Bernstein polynomials for the bicubic tensor product patch . Since the boundary curves are known, the control points of can be expressed in terms of the surface derivatives at the corners. For example
where and denote partial derivatives with respect to u, and u and v, respectively. Using the cross derivatives, the mixed partials are also expressed as
where dus and dvs denote the first nonzero knot spans of the cross derivatives in u- and v-directions, respectively. Taking to be the average of the two expressions, and substituting the equations of , and into the equation of one gets
Similar results hold for the remaining inner control points , and . Now, substituting into the Coons formulation above, the only unknown elements are the control points of the cross derivatives. Assuming that the subset of data points
had parameters assigned, ( denotes a random selection of indexes of the set {0,…,N} ) ; an over determined system of equations
can be set up and solved for the unknown control points. The system represents an unconstrained set of equations. In order to compute the Coons patch properly, twist compatibility must be ensured at the corners. This gives rise to the following constraint equations:
where due and dve denote the last nonzero knot spans of the cross derivatives in u- and v-directions, respectively. Putting the constrained and the above unconstrained systems together, the system can be solved using Lagrange multipliers as detailed in [26]. In order to make the above scheme useful for base surface computation, the following issues must be addressed: How many points to select? Practical experience shows that the number of points needed for a base surface that follows the data reasonably faithfully, heavily depends on the data. Simple data sets need no more than 10% of the original points, whereas more complicated ones may require up to 70%. How to select the points? Since the method computes cross-derivatives, it may appear that points close to the boundaries are better than those in the interior. Again, practical experience shows that this is not always the case. In our implementation we chose randomly distributed points. How many control points are necessary? Again, this is data dependent. At first, it may be a good idea to select the numbers of control points from the boundary curves. If the surface is not adequate, this number can be increased (sometimes decreased). Where do the parameters come from? The easiest way to deal with this is to put a bilinear Coons patch to the boundary curves, select a subset of points, and compute the parameters by projecting them onto the bilinear Coons patch. How to choose the knots? The best way to select knots is to use the ones of the boundary curves. This will ensure that no additional data is generated when computing cross derivatives. The above scheme works very well on most of the data sets exhibiting the kind of simplicity one expects from good segmentation algorithms. Unfortunately, some surface patches are more complicated, in which case the bicubic Coons tends to overshoot the data, i.e. it includes more wiggles than implied by the points. The reason for the overshoot is that the magnitudes of the cross derivatives are much larger than needed. Although the surface is a perfect Coons patch, it may not be a good base surface. A simple solution to this problem is to scale the cross derivatives, if needed. Let Lu and Lv be the average arc lengths of the opposite boundary curves in u- and v-directions, respectively. To scale, for example, , let mag be the average magnitude of the cross derivatives, then multiply the curve with a constant factor as follows:
This formula first scales the magnitude to become the average arc length, and then it scales again to account for the other directional derivative. After the scaling is done, the bicubic Coons surface will violate the derivative and twist conditions, e.g. the surface’s partial derivatives along the boundaries will be different from those of the boundary curves, however, this is irrelevant so long as the base surface is well parameterized and incorporates local details. Fig. 3(a) and (b) shows bicubic base surfaces with derivative scaling. In Fig. 3(a) 30% of the original points were used, and the number of control points were inherited from the boundary curves, i.e. (n, m) = (20, 8). Fig. 3(b) shows the case with 70% of the data points used and (n, m) = (50, 25). Notice that small improvement is gained at a significant cost, i.e. in general it is not worth increasing the number of data points over 30-50%.
2.3 Tensor product In some cases the parameterization process may require a tensor product surface, i.e. a surface fit is performed to a subset of points given their parameters. The same issues as in the bicubic Coons case apply. That is, we choose a percentage of the original points based on the complexity of the data; the points are randomly picked; the number of control points can be changed; and the parameters are obtained from the bilinear Coons patch. Fig. 4(a) and (b) shows two examples. In the first 30% of the points were used and (n, m) = (20, 8), (p, q) = (2, 2). The second example had the same percentage, however, (n, m) = (30, 15), (p, q) = (2, 2). The large number of data points (3455) and control points (496) required several minutes of compute time (including parameterization and fitting) to obtain the surface shown in Fig. 4(b). Nevertheless, it shows the power behind the tensor product surface fitting. Compare the result of this figure with that of Fig. 3(b). Even though that the Coons patch had many more data points as well as control points, the surface is not as good as the tensor product surface. However, computing the bicubic Coons is much faster because a lot smaller system of equations must be solved.
To conclude the base surface construction, the following is suggested to the reader: If the data is simple, use the bilinear Coons patch. It is the absolute fastest and the most reliable method. If internal details are to be reproduced, select no more than 30% of the data points to get a bicubic Coons surface. If all else fails and intricate details must be reproduced for a good parameterization, use the tensor product surface fitter.
3 Point projection Once the base surface is computed, the points are projected onto the surface to obtain the parameters. A standard way of doing that is a two step process 1. for each point and a good initial guess; and 2. perform Newton iteration until convergence is reached. Unfortunately, this is a very expensive and error prone process which can fail quite often especially for points near the boundaries. We propose a different kind of point projection. The idea is that precise Newton iteration is not needed considering the fact that the base surface is quite arbitrary. The method works as follows: 1. decompose the base surface into quadrilaterals using a relative geometric tolerance (we used 0.5% in all test examples); 2. project the points onto the closest quadrilateral; and 3. recover the parameter from the ones of the corners of the quadrilateral. Details of this method and several speed-up possibilities are explained in the subsequent sections.
3.1 Surface decomposition Given a (base) surface and a tolerance, a quadrilateral decomposition of the surface is sought so that each quadrilateral does not deviate from the surface more than the tolerance. Such decomposition is normally done adaptively, i.e. check if the surface is at; if not , subdivide and recurse. Although this method works well, it gives rise to cracks in the subdivision, i.e. corners of some quadrilaterals may lie on sides of others. Cracks, in turn, give rise to problems in point projection in that some points may actually project into the crack. For numerically stable point projection we needed a surface decomposition method that gives a crack free blanket covering the entire surface. The following method was implemented: Initialize a parameter stack with the parameters of the four corners. Add these parameters to the u- and v-parameter lists, which must be kept in sort order. While the stack is not empty do: Pop the stack, i.e. get the parameters of a quadrilateral. Extract a surface patch at these parameters. Check if the patch is flat. If not, then: - Compute the midpoint of the parameter span either in u- or in v-direction, depending on which direction to subdivide. - If this parameter is not yet in the parameter list, add it to the list keeping the list in sort order. - Subdivide the current parameter rectangle into two, and add them to the stack. Generate the output blanket by evaluating points at the computed parameters. A few words about flatness check. Since the parameterization process is quite an expensive task, and since most surfaces in reverse engineering fitting should not be too complex (after segmentation), an inexpensive flatness test was implemented. It works as follows: Given a NURBS surface patch. Compute an approximating plane to the four corners. This can be a least-squares plane or a plane passing through the center of gravity of the corners and is parallel to the mid-rulings of the hyperbolic paraboloid determined by the corners. Call this plane the mid-plane. Get the normalized implicit equation of the plane. Compute the distances of the control points from this plane by plugging their x; y; z coordinates into the plane’s equation (this is a cheap way of computing the distance of a point from a plane requiring only three adds and three multiplications). If for some point the distance is larger then the tolerance, stop. The surface is not planar and the direction of subdivision is the direction of larger extent. Otherwise, the surface is planar. The direction of subdivision can be obtained more accurately by computing discrete curvatures, however, that would be significantly more expensive. What we needed was (1) geometry based subdivision, (2) reasonable number of quadrilaterals, and (3) very low computational cost. Fig. 5(a) shows a subdivision of the base surface (shown in Fig. 4(b)) using 0.5% relative tolerance. In Fig. 5(b) the same surface is subdivided with 0.1% relative tolerance. Notice that (1) the subdivision is geometry based, i.e. high curvature areas receive more quadrilaterals, and (2) low curvature areas contain more quadrilaterals than necessary due to the crack-preventing nature of the process. This kind of subdivision, although not economical in terms of the number of quadrilaterals, is precisely what we want for obtaining the parameters, i.e. a uniform, geometry/curvature based distribution of quadrilaterals is more important to get good parameters than the minimum number. In fact, the quality of the parameterization can be improved by decreasing the sizes of the quadrilaterals.
3.2 Fast range searching The next step in the projection process is the projection of points onto the quadrilaterals. A brute force method, i.e. each point is projected onto each quadrilateral to find the closest one, is fairly expensive. For example, the subdivision in Fig. 5(a) has 1080 quadrilaterals, and the data set contains 11,517 points, giving rise to 12,496,680 point projections. One can do significantly better than that. The idea is to put the points and the quadrilaterals in a data structure for fast range searching. Since range searching in three-dimensional (3-D) is expensive, and since the surfaces are in fact topologically two-dimensional (2-D) entities, point projection to a plane is used. It is assumed that the surface does not curve more than 1808 (the point cloud is projectable to a plane), e.g. the method does not work for surfaces such as a complete cylinder or a sphere. The segmentation module must subdivide such surfaces into simpler patches. The first idea that pops up is to use the well-known leasts squares plane. Although it is easy to compute, it tends to behave quite badly especially for thin and long data sets. A better plane is the so-called maximum area plane, i.e. the plane on which the convex hull of the projection forms a maximum area. To obtain such a plane, we use the base surface’s decomposition. The derivation below is based upon the ideas and notes of Joe Klahs, which we respectfully acknowledge. Assume that the base surface is decomposed into quadrilaterals in the form
the normal N of the maximum area plane is sought. Using the diagonals of each quadrilateral in a consistent manner, the dot product
is the area (vector) projected onto the normal (vectors are written as column vectors, i.e. their transpose (T) is needed for product computation). Collecting all triangles, yields the system
To get the maximum area, one needs the minimum projected area, that is, the minimization problem
is solved. This leads to the following eigenproblem
where λi are the eigenvalues associated with the eigenvectors Ni; i= 1, 2, 3. A quick way of solving this eigenproblem is by iteration
If one normalizes Nk+1 after each iteration, the iteration above quickly converges to the eigenvector associated with the largest eigenvalue. In all of our practical examples, the iteration converged within five steps using an angular tolerance of 1。 (which is more than adequate for practical applications). Fig. 6(a) shows the maximum area plane of the base surface and the projection of the quadrilaterals. The data points and their projections are depicted in Fig. 6(b).
Once the maximum area plane is computed and the points are projected, the plane is put in a position parallel to the [x,y] plane, and the problem is transformed into a 2-D search problem. Fig. 7 shows a 2-D view of the points and the quadrilaterals (the base surface was subdivided using a relative tolerance of 2%). The search proceeds as follows: Compute the minmax box [xl,xr]×[yb,yt] of the projected quadrilaterals and of the points. Put a grid over the data set as follows. Let nov be the number of vertices in the surface decomposition, and denote the lengths of the sides of bounding box by xL and yL: Then the size of the grid is
which gives the grid resolutions
Bin the projected points into the cells, i.e. for each projected data point [x; y], find the cell indexes
and register the point with that cell. Note that quadrilateral vertices need not be binned. Initialize a distance array for each point with a large number. For each quadrilateral in the base surface decomposition do: Find the cells covering the quadrilateral. This can easily be done by computing the indexes of the cells the vertices fall in. For each cell do: - Retrieve points registered with the cell. - Project the point onto the quadrilateral. - If the projection falls within the boundaries, and if the distance of the point from the quadrilateral is less than the current distance, then - Compute the parameters from those of the vertices of the quadrilateral. - Update the current distance.
In general each quadrilateral is covered by a few cells, e.g. 4-6, which in turn contain a small portion of the data set. Using the example of Fig. 7, the points are binned into about 160 bins, which in turn contain about 300-350 points per quadrilateral, i.e. only 2-3% of the data set is needed to process a quadrilateral.
3.3 Thinning In many practical applications data acquisition devices produce much more data than necessary for accurate surface reconstruction, i.e. a subset of the data is quite satisfactory. Such a subset is obtained by the process called thinning. Our method is based on the decomposition of the base surface. The idea is to find the nearest point for each vertex in the quadrilateral decomposition. That is, the decomposition, mentioned in Section 3.1, uses two kinds of tolerances: (1) for point projection; and (2) for thinning. The point projection tolerance (set in the system at 0.5% relative) is hidden from the user; however, the thinning tolerance is an input allowing the user to thin out as many points as necessary. Finding nearest points for each quadrilateral vertex proceeds as follows: Use the data structure set up in the previous section. For each vertex in the quadrilateral decomposition do: Find the cell the point is in. Find the distance between the point and the closest cell wall. Call it the wall distance. Set the minimum distance to a large number. While the minimum distance is greater than the wall distance do: - Find the cell ring around the cell of the vertex. - Find the closest point in all the cells in the ring. - Reset the minimum distance to be the distance between the vertex and the closest point. - Increase the wall distance by the size of the grid. This method first searches the cell the vertex is in. Then it searches rectangular rings of cells around the first cell. The search stops when the distance between the closest points is smaller than the distance between the vertex and the closest wall of the last ring. For most practical applications this algorithm hardly ever goes out of the first ring, i.e. the maximum number of cells it searches is 9. If lots of points are to be thinned out, then considering that the point distribution is much denser than the vertex distribution, only one cell, i.e. the cell the vertex is in, must be searched. Fig.8 (a) shows thinning of the input data using 0.5% relative tolerance. The original 11,517 points have been reduced to a set of only 1052. Fig.8 (b) illustrates thinning with 0.1% tolerance yielding a subset of 4647 points. Notice that the thinning is geometry based, i.e. in ¯at areas more points are thinned out than in high curvature areas.
3.4 Computing the parameters The final step in the parameterization process is to compute the parameters from the projections onto the quadrilaterals. If the points have been thinned, then only points in the subset are considered for projection, i.e. thinning is a preprocessing step where selected points are marked. In the projection phase the unmarked points are skipped whereas the marked ones are projected. A few words about projecting onto quadrilaterals. Four points in space determine a hyperbolic paraboloid, i.e. a bilinear surface. Projecting onto such surface can be done either numerically, e.g. using Newton iterations, or geometrically. Since we introduced surface decomposition to do away with Newton iteration, numerical projection is out of the question. Geometric projection is viable in the following sense. The opposite sides of the quadrilateral (when no planar) determine two plane directions. These planes are called the generating planes which can be used to slice the quadrilateral in a parallel manner to obtain its rulings. Putting the planes through the point to be projected, gives two rulings determining a point on the paraboloid which may be considered as the projection. This works quite well for points on or near the paraboloid, and the parameters can be recovered very quickly. However, if the paraboloid is quite ¯at and the point is far from it, this kind of appreciative projection gives fairly bad parameters. Therefore an appreciative method was applied, which is justifiable because the decomposition is done fine enough (in our experience 0.5% proved to be quite good); we are looking for reasonable parameters, not parameters associated with specific points with respect to a given surface, i.e. point inversion. The method is simple and works as follows: Get the mid-plane of the quadrilateral (as above). Project the vertices onto this plane. This gives a planar quadrilateral. Project the data point onto this plane. If the projection falls within the planar quadrilateral, retrieve the parameter based on those of the vertices (this can easily be done based on the two sets of planar rulings emanating from the intersections of the opposite sides). This process is very fast, reliable, i.e. no numerical computation is needed, and produces good parameters. An added feature is that it can be used for point inversion as well in essentially two ways. The one is to compute good start parameters for fast Newton iteration, and the other is to use a fine surface decomposition for approximate point inversion. Fig. 9(a) and (b) shows fitting examples using the parameters obtained by the process outlined above. Fig. 9(a) shows a surface ®t to 1052 points obtained by thinning out the original data with a tolerance of 0.5%. In Fig. 9(b) the surface was fit to 8061 points, obtained by thinning to 0.05% tolerance; the computational time was quite significant.
4 Conclusions Several ideas and methodologies for obtaining parameters for surface fitting to random data were presented in this paper. The bottom line of parameterization is to have simple surfaces with as little internal detail as possible. That is, the most critical part of the reverse engineering process is segmentation. Well segmented data sets are easy to thin out, and the thinned set along with the simple boundaries produce high quality surfaces at a reasonable cost. However, if segmentation fails to provide simple patches, and if intricate details are to be reproduced, as in the figures of this paper, then the only alternative is to use complicated base surfaces, e.g. bicubic Coons or tensor product surfaces, to interactively adjust their parameters, e.g. number of data points and control points, and to possibly iterate from one base surface to another. That is, the first set of parameters may be used to construct another base surface, which in turn is used to construct a second set of parameters, and so on. There are two major problems with this process: (1) it is extremely time consuming; and (2) if the initial base surface is poor, due to the complexity of the data, then the second one can be even worse, and the entire iterative process produces nonsensical surfaces.
Acknowledgements The authors thank Prof. C.-H. Menq, Prof. W. Ma, and Dr G. Lukacs for their many valuable comments and criticisms.
逆向工程中曲面拟合的参数化
摘要: 给定四条边界曲线和表面贴片上的一系列随机的点,那么就可以提出一种给这些点赋参数的方法。这种算法是用许多基本曲面并把这些点投影到这些面上来恢复基于这些面的基本参数的一些参数。同时也出现了一些加快这个投影过程的技术。一种稀释方法也同样被引进,这种方法选择原有点的一个子集作为取样点,速率比实际需要更快,这种稀释方法基于基本曲面的几何特征并且依赖于一个有实际使用价值的形位公差。
关键字:曲面拟合;逆向工程;非均匀有理B样条曲面
1. 引言 从取样的数据中进行物体重建是获得一个已有物体的设计参数的一种重要组成部分。图形模型[1-3,7-14]以及工程零件[15,27,28,30-33]必须被如实的重建出来,也就是说,可以寻求到一种表现形式来定义一个物体使其在给定的误差范围内不偏离这些点。这种表示法可以是离散的,例如一组三角形,或者是连续的,例如边界表示法模型。本文的中心是通过对点的子集进行曲面拟合来获得边界表示法模型。对随机分布的点进行曲面拟合就成了有重大利害关系的一个主题[4-6,10,16-25,29,33]。使用B样条曲线的话,那么这个拟合过程的关键就是参数和结点向量的计算了。本文主要是解决第一个问题。更精确地说是要解决下面的问题。假设数据被分割成由曲线界限的四边形区域。更确切地说,给定四条非均匀有理B样条曲线作为边界(Fig.1) ; 寻求一系列随机分布的点(假定边界上点是被剥离的) (在Fig.1中 );一组参数( ) 来最好地表示矩形区域中点的位置。也就是说,在这些参数方面这些数据点是假定在非均匀有理B样条曲线的表达式中的。一种解决这个问题的可行方法是开始于一个基本面,以某种方法拟合到边界曲线和一些点并且将这些数据点投影到这个面上来获得这些投影点的参数。理论上来说这些操作似乎都是做起来很简单的,然而这其中有很多技术的障碍可能使这种基于基本面的方法完全不可用。
这篇文章的贡献分两部分:(1)介绍各种基本曲面来获得更详细的内在数据;(2)一种可以在点投影过程中获得显著加速的方法。本文组织如下。第二部分主要涉及一些基本曲面,第三部分主要是在点投影过程中获得稳定性和加速,第四部分则是对本文做一总结。
2. 基本曲面 获得参数的一般方法是用曲面依附到输入的参数上[21]。在下面的这部分中,我们主要研究这三种曲面:(1)双线性混成昆氏曲面; (2) 双三次昆氏曲面;(3)张量乘积曲面。 2.1 双线性昆氏曲面 最简单的基本曲面被称为双线性混成昆氏曲面。这种曲面用非均匀有理B样条曲线可以表示为如下所示[26]:
其中: 表示一个u方向的线性举升曲面(也被称为直纹面); 表示一个v方向的线性举升曲面; 表示一个张量乘积曲面,也就是说,它是一个由四个角点定义的双线性曲面。 一旦这两个直纹面和那个双线性曲面被构建出来,它们就会被谐调。通过加上直纹面的控制点及减去双线性曲面控制点的方法来获得 的控制点。 使用双线性昆氏曲面的优点为: 计算简单; 计算过程稳定、快速; 曲面更易用参数表示。 如果数据设置非常简单(一个好的逆向工程软件要包括一个分割模块用来将数据细分成非常简单的片),那么双线性昆氏曲面就是我们所能想到的最好的基本曲面。然而,这种曲面有一个缺点就是它会丢失一些内部详细信息,如Fig.2中所示。如果分割过程不能很好地进行,那么双线性昆氏曲面就不能提供详尽的参数。
2.2 双三次昆氏曲面 使用双三次昆氏曲面的意图就是让基本曲面包含一些内部详细信息。双三次昆氏曲面除了举升曲面不是线性而是三次混成的外,其他方面都和双线性昆氏曲面相同。这样做,我们就需要截面的导函数。这样一个曲面的计算是非常复杂的,因此下面只略述一下其大意。 假定除了边界曲线之外,相当数量的内部点也被选中,固定了控制点的数量和两个方向的角度后,未知量就是截面导函数了。
其中, 是穿过边界 截面导函域,也就是说,它是u的一个函数并且它的通过多个参数评价的向量就是导函数,其指向是沿着边界曲线的。为了计算这些截面导函数,就要寻求一些未知的控制点。
简单来说,固定了参数后曲面 就可用定义在值域为[0,1]的厄密多项式表示:
如果在一般域 上来定义曲面,那么以下这些参数变换便可以使上面那些域统一跨度。
注意: 现在,对于一个固定的参数组 ,双三次昆氏曲面就可以写作:
这里, 和 是双三次张量乘积曲面 的伯恩斯坦立方体多项式。既然边界曲线已经知道了,那么 的控制点就可以根据曲面在拐点处的导函数来表示了。例如:
这里, 和 分别表示关于 的偏函数。使用截面导函数,那么混合的偏函数同样可表示为:
这里, 和 分别表示截面导函数在u和v方向的第一非零结点跨度。将 看作二者的平均并且将 , 和 的等式用 的等式代替就可以得到:
剩下的内部控制点 , 和 也可以得到相似的结果。现在把 用上面的昆氏公式代替,就仅有截面导函数的控制点是未知的了。 假定数据的子集 ( )表示在数集{0,…,N}中随机选择的数字)被指定了参数 ,那么那些未知控制点就可用一个超定系统的等式 来提出并解决。这个系统表示一组不受约束的等式。为了正确的计算昆氏曲面,必须要保证拐点处的扭曲兼容性。这就导致了下面约束等式的出现:
这里, 和 分别表示截面导函数在u和v方向上的最后一个非零结点的跨度。将上述约束等式和不受约束的系统整合,那么系统就可以用[26]中介绍的拉格朗日乘子来解决。 为了使上述方案对基本曲面的计算有所帮助,那么就必须解决下面这些问题: 选择多少点?实际经验告诉我们一个能够适度如实地遵循数据的基本曲面需要的点数在很大程度上是和数据有关的。简单的数据需要的点不超过原始点数的10%,然而复杂的可能会高达70%。 怎么选择点?既然这种方法要计算截面导函数,看起来好像靠近边界的点要比那些内部的点更好。但实际经验告诉我们事实并非总是如此。在执行的过程中,我们选择随机分布的点。 多少控制点是必须的?这也是和数据相关的。起初,从边界曲线上选择控制点可能是个比较好的主意,但如果曲面不够的话,这个数量可能会增加(有时也会减少)。 参数从哪里来?解决这个问题最简单的方法就是把一个双线性昆氏曲面放到边界上,选择一个点的子集并且通过将它们投影到双线性昆氏曲面来计算那些参数。 怎样选择结点?最好的方法就使用边界曲线。这样就会保证在计算截面导函数的时候不会产生冗余的数据。 上面这种方案在展示简单曲面的时候效果比较好。但有一些曲面片是非常复杂的,在这种情况下,双三次昆氏曲面就会超出原有的数据,也就是说,它包含比原有点中更多的摆动。超出的原因是导函数的量比实际需要的大得多。尽管这个曲面是个比较完善的昆氏曲面,但它却不是一个好的基本曲面。如果需要的话,一种比较简单的解决这个问题的方法是按比例缩放截面导函数。 和 分别为u-和v-两个反方向边界曲线的平均长度。为了缩放 ,设mag为截面导函数的平均量,那么将原有曲线乘上一个常量如下所示:
这个公式首先将量级缩入成平均弧长,然后又在另一方向的导函数上进行缩放。缩放后双三次昆氏曲面就会偏离导函数和扭曲条件,例如,沿着曲面边界的偏导数就会和边界曲线有所不同,但只要基本曲面能够很好地参数化并且表现出局部的详细信息,那么这些就是无关紧要的。Fig.3(a)和(b)显示的是导函数缩放后的双三次基本曲面。在Fig.3(a)中用到了原始数据的30%并且控制点来自于边界曲线,也就是说,(n, m)=(20,8)。Fig.3(b)显示的是用到原始数据70%的情况并且(n, m)=(50,25)。需要注意的是,用较大的成本仅获得了很小的改善,也就是一般情况下并不用增加数据点到30%—50%以上。
2.3 张量乘积曲面 在一些情况下参数化的过程可能需要一个张量乘积曲面,也就是一个表示成可以给出它们参数的点的曲面拟合。这和在双三次昆氏曲面中涉及的问题相同。更确切地说,我们根据原始数据的复杂情况选择一定比例的点;点的选择是随机的;控制点的数量是可变的;从双线性昆氏曲面上获得参数。Fig.4(a)和(b)给我们显示了两个例子。在第一例子中使用了原有数据的30%并且(n, m)=(20,8),(p, q)=(2,2)。第二个例子的取点比例和第一个相同,但(n, m)=(30,15),(p, q)=(2,2)。大量的数据点(3455)和控制点(496)需要几分钟的时间(包括参数化和拟合)来获得如Fig.4(b)中所示的曲面。然而,它体现的是张量乘积曲面拟合背后的力量。拿这个图形的结果和Fig.3(b)对比一下,尽管昆氏曲面包含更多的数据点和控制点,但它还是比不上张量乘积曲面。可是,双三次昆氏曲面的计算要快得多,因为它只需计算更少的系统方程。
为了总结一下基本曲面的构建,建议阅读一下以下这些内容: 如果数据比较简单,提倡使用双线性昆氏曲面。这绝对是最快和最可靠的方法。 如果要再生内部的详细信息,建议选择不超过30%的数据点来生成一个双三次昆氏曲面。 对于一个好的参数化过程来说,如果所有其它的方法都不行并且要再生一些复杂的详细信息,就要使用张量乘积曲面了。 3. 点的投影 一旦基本曲面计算出来了,就可以把点投影到上面来获得参数。做这个的标准过程一般要分两步。 给每个点找到一个好的初始猜想。 使用牛顿迭代法直到达到收敛状态。 但是,这是一个非常耗时和经常出错的过程,特别对于靠近边界上的点它会经常失效。在这里我们建议另外一种点投影的方法。考虑到基本曲面是比较任意的,所以并不需要精确的牛顿迭代法。这个方法的工作过程如下: 1) 使用一个相对的几何公差将基本曲面分解成一些四边形(在下面的实例中我们使用的是0.5%); 2) 将点投影到距离它最近的四边形上; 3) 从四边形的拐角处重新获得参数。 在后面的部分中我们会详细介绍一下这种方法的详细内容和一些加速的可能性。
3.1 曲面分解 给定了一个基本曲面和公差后,将基本曲面进行四边形分解,所以在这个公差的范围内每个四边形都不会偏离基本曲面。这个分解的过程通常都自适应完成的,也就是说: 检查面是否是平的; 如果不平就再次细分并且重复上一个过程。 尽管这种方法工作良好,但它在细分这步也有一些缺陷,就是,一些四边形的拐角可能位于另外一些四边形的边上。反过来这些缺陷在点投影的过程中又会产生一些问题,一些点可能会投影到这些缺陷上来。对于一个数字上稳定的点投影过程我们需要一个能够提供无缺陷的曲面分解方法来覆盖整个曲面。这种方法的执行如下所示: 用四个拐角处的参数初始化一个参数栈。将这些参数添加到u-和v-方向的参数列表中,这个列表是用类别排序的。 当栈不为空时就做以下操作: 取中栈中的数据,就是得到四边形的参数。 在这些参数处解析曲面片。 检查这些曲面片是否平坦,如果不平,就做以下操作: - 根据要细分哪一个方向,计算u-或v-任一方向的参数跨度的中点。 - 如果这个参数还没有在参数列表中,就把它添加到列表中,并还按类别排序。 - 将当前的参数长方形一分为二,并把它们添加到栈中。 通过分析在被计算处参数的点来生成总括输出。 还有一些文字关于平整度的检查。既然参数化的过程是个比较耗时的工作并且逆向工程中大多数曲面的拟合是不需要太复杂的(分割后),所以就要执行一个不太耗时的平整度检查。它工作如下所示: 给定一个非均匀有理B样条曲面片。 计算一个近似平面的四个拐角。这可以是一个最小平方平面或者是一个通过拐角的重心并且平行于由拐点决定的双曲线抛物面的半轴的平面。称这个平面为中间平面。 得到平面的规范化隐式方程。 通过把控制点的x, y, z的坐标代入平面的方程中来计算这些控制点与平面的距离(这是耗时最小的计算点与平面距离的方法,只需要三个加法和三个乘法运算)。 如果发现一些点的距离超过了误差,就停下。这个面不是平面并且再分的方向是更大跨度的方向。 否则,这个面就是平面。 可以通过计算离散的曲率来获得更精确的细分方向,但是这个过程也更加耗时。我们需要的是:(1)基与几何的再分,(2)合理的四边形数量,(3)非常低的计算成本。Fig.5 (a)中展示的就是使用0.5%的相对误差来对基本曲面(Fig.4(b)中的曲面)的细分。Fig.5(b)中的是相同曲面用0.1%的相对误差进行细分。注意(1)这个细分是基与几何的,就是,大曲率的区域有更多的四边形,因为这个过程有防止出错的特性。这种再分方法尽管在四边形的数量方面不是很经济,但它可以使我们得到精确的参数,就是说对于要得到一个好的参数来说,一个统一、基与几何或曲率的四边形分配比最小数量更重要。实际小,可以通过缩小四边形的大小来提高参数化的质量。
3.2 快速范围搜索 投影过程的第二步就是把点投影到四边形上去。一种严格的方法是把每个点投影到每个小四边形上去,从而找到最近的一个,这种方法是非常耗时的。例如,Fig.5(a)中的细分有1080个四边形,并且数据中包含11517个点,也就平生了12496680个点投影。有另外一种方法可以比上面做得更好。 这种方法就是把那些点和四边形放在一个数据结构中用来快速范围搜索。尽管在三维空间里范围搜索是非常耗时的并且实际上那些面在拓扑上是二维的,但还是使用了点投影到曲面。假设这些曲面弯曲都不超过180度(点云可以投影到曲面上),例如,这种方法不适合完全是一个圆柱或是球体的曲面。分割模块必须把这些曲面分成更为简单的面片。 首先想到的就是使用众所周知的最小二乘曲面。尽管它更容易计算,但它却表现得不尽如人意,特别是稀薄狭长的数据集。一种更好的曲面被称为最大面积曲面,就是曲面所在的投影凸起外壳构成最大面积。为了获得这样的曲面,我们要用到基本曲面的再分。下面的推导是基于Joe Klahs的思想和笔记,在这里我们给予尊敬地承认。 假设基本曲面用这种形式再分成四边形,那么就可以寻求到最大面积曲面的法向N。用一致的方式使用每个四边形的对角线,点积 就可以表示面积投影到法向(向量被写成了列向量的形式,就是点积的计算需要它们的转置)。收集所有的三角形,产生了下面这个规律:
为了得到最大面积,就需要最小投影面积,就是说,最小限度问题被解决了。这导致了下面这个特征值问题的产生, ,这里 表示特征向量 的特征值 i=1,2,3 。一种快速解决特征值问题的方法就是利用迭代:
如果每次迭代后都规格化 ,上面的迭代就可以快速收敛到最大特征值相关的特征向量。在我们所有的实例中,迭代都可以以不大于1度(这对于实际应用已经足够了)的误差在5步之内收敛。Fig.6(a)显示的是基本曲面的最大面积曲面和四边形的投影。数据点及它们的投影在Fig.6(b)中被描述。
一旦最大面积平面被计算出来并且点被投影上去后,曲面就被置于平行于[x,y]的平面上了,问题就变成了二维搜索问题。Fig.7显示的就是二维视图下的点和四边形(基本曲面被用2%的相对误差细分)。搜索过程进行如下: 计算投影四边形和投影点的最小最大盒[xl, xr]×[yb, yt]。 在数据上放一个如下的栅格。设nov为曲面再分的顶点数,并且用xL和yL表示包容盒的边长。那么这个栅格的大小为:
它给出了栅格的解析率: 把投影点放到每个单元中,就是给每个投影点(x, y)找到一个单元索引
并且用该单元注册该点。注意四边形的顶点无须放到单元中。 用一个大数来初始化每个点的距离数组。 对于基本曲面上的每个四边形做如下操作: 寻找覆盖四边形的单元。这步可以通过计算顶点落入的单元索引来很容易完成。 对于每个单元做如下操作: - 取出单元中的点。 - 把点投影到四边形上。 - 如果点投影在边界内并且点与四边形的距离小于当前的距离,那么 - 计算这个四边形顶点的参数。 - 更新当前距离。
一般说来,每个四边形都被几个单元所覆盖,如4—6,依次包含了数据集中的一小部分。拿Fig.7这个例子来说,点被分成了160份,每个四边形包含300—350个点,就是说,处理一个四边形只须用到数据集的2—3%。 3.3 稀释 在许多实际应用中数据获取设备产生了比所需的能够重建准确曲面更多的数据,就是说,它的一个子集是比较满意的。这样一个子集的获取就被称作稀释。我们的方法是基于基本曲面的再分的。这个方法就是在四边形分解时给每个顶点找到一个最近的点。也就是,在3.1中提到的分解用到两种误差:(1)用来点投影;(2)用来稀释。点投影误差(系统中设置为相对0.5%)对于用户来说是隐藏的,然而,稀释误差是一个可以输入的量从而允许用户根据需要来稀释。为四边形每个顶点找最近的点进行如下: 使用前面讲到的数据结构。 对于四边形分解中的每个项点做如下操作: 找到点所在的单元。 找出点与最近的单元墙之间的距离,称之为墙厚。 用一个大数来设置最小距离。 当最小距离大于墙厚时,做如下操作: - 找到围绕在顶点单元周围的单元环。 - 找到所有单元中最近的点。 - 将最小距离设成顶点与最近点之间的距离。 - 将墙厚增加一个栅格的大小。 这个方法首先搜索顶点所在的单元。然后它会搜索第一个单元周围的矩形环。当最近点的距离小于顶点与最近点的距离时,搜索就会停止。对于大多数实际应用来说,这个算法很少能够跳出第一环,就是说它搜索的最大单元数为9。如果有太多的点要被稀释,就要考虑到点的分布可能比顶点分布密集,只有一个单元,就是顶点所在单元必须被搜索到。Fig.8(a)显示的相对误差为0.5%时的稀释。原始11517个点减少到只剩下1052个。Fig.8(b)显示的是用0.1%的误差产生4647个点。注意这里的稀释是基于几何的,就是说在平面区域比在大曲率区域更多的点被稀释。
3.4 计算参数 参数化过程的最后一步就是从投影到四边形中计算参数。如果点被稀释过,那么只有子集中的点被考虑用来投影,就是说稀释是标记所选点的前处理步骤。在投影过程中未被标记的点就被略去而被标记的点就要被投影。 还有一些关于投影到四边形的文字。空间中的四个点决定一个双曲线抛物面,就是双线性曲面。投影到这样一下曲线既可以用数字方法来做,例如用牛顿迭代法,或者用几何方法。既然我们在介绍曲面分解时没有用牛顿迭代法,那么数字投影就不可能了。下面几何投影是可行的。四边形的反面(当不是平面时)决定两个面的方向。这些面被称作发生面,它们将四边形用平行的方式划分从而得到引导线。将面穿过被投影的点并且给出两条引导线来决定哪些双曲线上的点将被投影。这种方法对于在或者临近双曲线的点动作良好,并且参数可以被迅速恢复。然而,如果双曲线抛物面比较平并且点离它比较远,这种近似的投影就可能会产生比较不好的参数。因此,一种近似的方法被采用,它之所以正确是因为: 分解过程做得足够好(在我们的经验中证实0.5%就比较好); 我们寻求的是合理的参数,并不是和一个给定曲面上特定点相联系的参数,就是点倒置。 这种方法比较简单,流程如下: 得到四边形的中间面(如上)。 投影顶点到面上。这给出了一下平面四边形。 将点投影到这个平面上。 如果投影落到了平面四边形内,就取出基于这些顶点的参数(这很易由两条从相对边的交点散发出来的引导线来完成)。 这个过程是非常快速、可靠的,就是说不需要数字计算,并且还可以产生好的参数。还有一个附加的特性就是它也可以以两种方式被点倒置所用。一种是为牛顿迭代法计算起始参数,另外一种就是用一个比较好的曲面分解来进行近似点倒置。 Fig.9(a)和(b)显示的是用上面叙述的过程获得的参数拟合曲面的例子。Fig.9(a)显示的是用0.5%的公差和稀释原始数据后获得的1052个点拟合成的曲面。Fig.9(b)中显示是用0.05%的公差和稀释原始数据后获得的8061个点拟合成的曲面,计算时间是非常显著的。
4. 结论 在这篇文章中介绍了一些关于曲面拟合到随机数据从而获得参数的思想和方法。参数化的底线就是得到含有尽可能少的内部信息的简单曲面。也就是,逆向工程最重要的部分就是分割。好的分割数据很容易被稀释,而且稀释后的数据沿着简单的边界可以在一个合理的代价内产生高质量的曲面。然后,如果分割过程不能产生简单面片并且再生了一些复杂详细信息,就像这篇文章里的数据一样,那唯一的替代方法就是用复杂的基本曲面,例如双三次昆氏曲面或张量乘积曲面,来交互的调节参数,例如点和控制点的数量,并且尽可能从一种基本曲面迭代到另外一种。也就是说,第一组参数可能被用来构建另外一个基本曲面,它反过来又被用来产生第二组参数,如些下去。这个过程有两个主要问题:(1)它是极其耗时的;(2)如果初始基本曲面不好,由于数据的复杂性,第二个可以更不好,而且整个迭代过程会产生没有实际意义的曲面。
致谢 非常感谢Prof. C.-H. Menq, Prof. W. Ma, t和Dr G. Lukacs 的很多宝贵意见和评价。
|