Quick Overview of the Algorithm

In order to perform a free-form deformation on models composed of point sets (such as polygonal or curve surface models composed of vertices or control points respectively), we can surround a model (or just a portion of a model) with an oriented bounding box (OBB, as opposed to axis-aligned bounding box - AABB) that is defined by a center point, three axes (X,Y,Z), and a half-width or radius along each axis. The bounding box is then subdivided along each of its axes to form a lattice structure. Each point in the model point set (either vertices or control points) is given 3D normalized coordinates (s,t,u) in the bounding box's space. The control points of the lattice structure can then be repositioned and each point's new position can be calculated using a trivariate tensor product Bernstein polynomial function.

Preparing the Surrounding Lattice Structure

The given bounding box (OBB) becomes the all-positive octant in the 3D cartesian coordinate system where each contained point can be defined by normalized coordinates (s,t,u) where 0 <= s,t,u <= 1.

The OBB is given to us as follows:

• 3D center point (Cntr)
• 3 3D normalized vectors corresponding the axes X, Y, and Z
• 3 values corresponding to the half-dimensions or radius along each axis (HalfDimX, HalfDimY, HalfDimZ)

Throughout this paper we will be assuming a right-handed coordinate system with the positive x-axis to the right, y-axis up, and the z-axis towards the viewer (we are looking down the z-axis in the negative direction). In order to make the OBB correspond to the all-positive octant, we need to make the left-bottom-back (-x, -y, -z) vertex of the OBB the origin of the lattice system. The eight corners of the OBB will define the extents within this octant.

Given the surrounding OBB, the origin of the lattice (the left-bottom-back vertex) is defined as follows:

LatticeOrigin = Cntr - X*HalfDimX - Y*HalfDimY - Z*HalfDimZ

In order to find the points of the lattice structure, we must be given the number of intervals to subdivide along each axis: NumIntX, NumIntY, NumIntZ. These numbers also correspond to the number of times we must step along each axis starting at the origin in order to find the lattice points. We now need to find the step vector along each axis to each subsequent lattice control points.

Given number of intervals along each axis and the lattice origin point, we can find the remaining lattice points:

Allocate a 3D array of lattice control points: 

LatticeCtrlPts[NumIntX+1][NumIntY+1][NumIntZ+1]

Calculate actual dimension of the OBB (will be used throughout):

DimX = HalfDimX*2.0
DimY = HalfDimY*2.0
DimZ = HalfDimZ*2.0

Calculate step vectors along each axis:

StepInX = (X*DimX) / NumIntX
StepInY = (Y*DimY) / NumIntY
StepInZ = (Z*DimZ) / NumIntZ

Find lattice control points by looping through all three dimensions and stepping, starting at the origin:

TempX = LatticeOrigin;
for (int i=0; i<=NumIntX; i++)
{
  TempY = TempX;
  for (int j=0; j<=NumIntY; j++)
  {
    TempZ = TempY;
    for (int k=0; k<=NumIntZ; k++)
    {
      LatticeCtrlPts[i][j][k] = TempZ;
      TempZ += StepInZ;
    }
    TempY += StepInY;
  }
  TempX += StepInX;
}

The code to find the lattice points may at first be a bit confusing since we have to use so many temp variables. These lattice points are scanned out much like filling in a 2D array; the temp vectors are used to ensure that the next row or complete level starts back at the left or bottom-right respectively.

Finding Point (s,t,u) Coordinates in the Surrounding Lattice Structure

Now that we have the origin of the lattice (LatticeOrigin) and the dimensions and axes of the surrounding OBB, we can easily calculate the (s,t,u) lattice coordinates of any (x,y,z) vertex or control point as follows:

xyz = (x,y,z) - LatticeOrigin
(s,t,u) = { (xyz*X)/DimX, (xyz*Y)/DimY, (xyz*Z)/DimZ }

xyz is the (x,y,z) point is the world vector within the OBB. This vector is projected onto the OBB axes (dot product) and then normalized to obtain lattice space (s,t,u) coordinates.

Moving the Lattice Control Points

Move the lattice control points to deform the contained object. The lattice forms a 3D control mesh analagous to the control polygon for a 2D bezier curve. Deformations of the control mesh will result in a corresponding deformation in the object when we remap the (s,t,u) point coordinates back into world (x,y,z) coordinates.

Mapping (s,t,u) Coordinates back into World (x,y,z) Coordinates

1D FFD:

Since we are allowing for an arbitrary number of subdivisions along each axis, we must be able to handle arbitrary degree trivariate Bezier hyperpatches in order to map the space of (s,t,u) coordinates back into (x,y,z) coordinates as a function of the lattice control points. These hyperpatches are based on the Berstein polynomial blending functions:

B_i,n(t) = n!/(i!*(n-i)!) * tⁱ * (1-t)^n-i

Where n is the degree of the polynomial and n+1 is the order. The n+1 nth degree Berstein blending functions correspond to each of the n+1 control points required for the degree n Bezier space curve; the ith Berstein polynomial, for all i : 0 <= i <= n, corresponds to the ith control point. For example, for a 2D cubic Bezier curve defined with four control points (P0,P1,P2,P3), we require the set of third degree polynomials (n=3) which blend the control points to find a particular point along the curve at time t : 0<=t<=1 (parametric "time") as follows:

PtOnCurve(t) = P0*B_0,3(t) + P1*B_1,3(t) + P2*B_2,3(t) + P3*B_3,3(t)

Each Bernstein polynomial weights each control point at a particular time t. The set of Bernstein polynomials for any degree n always sum to unity (one) for any value of t in [0,1], Also, B_0,n(0)=1 and B_n,n(1)=1 so the endpoints are interpolated at the start and end of the time interval since the other blending functions are zero at t=0 and t=1 respectively. In general the 2D space curve is defined as:

PtOnCurve(t) = SUM_i=0,n { P_i * B_i,n(t) }

In relation to FFD's, the 2D Bezier space curve is a 1D FFD. We can define 1D t-points along a line represented by a collinear set of control points; the number of control points n+1 refers to n interval subdivisions along the single 1D axis. We can then move the control points and find the t-point's new position. This will be a little clearer when we define a 3D Bezier surface which can be though of as a 2D FFD.

2D FFD:

Given an n+1 x m+1 rectangular array of control points P_i,j where 0<=i<=n, 0<=j<=m, we can define a 3D Bezier surface as follows:

PtOnSurface(u,v) = SUM_i=0,n { B_i,n(u) * SUM_j=0,m { B_j,m(v) * P_i,j } }

Even though the control points can be of arbitrary degree (in this case, mostly 3D for a 3D surface), we will restrict the 3D Bezier surface to be composed of 2D control points. Points on the surface are just 2D points in a mesh of control points. This is analagous to our previous FFD discussion, where this mesh of control points forms a lattice that subdivides the bounded 2D rectangular region. This 2D FFD can be used to warp or "morph" 2D points that are mapped into the (u,v) space of the lattice. If we give a set of points normalized (u,v) coordinates in this lattice, we can move the control points, apply the PtOnSurface function, and obtain new "warped" (x,y) points. This strategy extends to the 3D FFD quite easily.

3D FFD:

Now finally, given an m+1 x n+1 x o+1 (NumIntX+1 by NumIntY+1 by NumIntZ+1) cubical array of 3D (x,y,z) control points, P_i,j,k (LatticeCtrlPoints[i][j][k]), we can map point (s,t,u) coordinates back into world (x,y,z) coordinates as follows:

PtInLattice(s,t,u) = SUM_i=0,m { B_i,m(s) * SUM_j=0,n { B_j,n(t) * SUM_k=0,o { B_k,o(u) * P_i,j,k } } }

where m=NumIntX, n=NumIntY, o=NumIntZ, P=LatticeCtrlPoints.

Algorithm for the (s,t,u) to (x,y,z) Mapping

As seen above, the PtInLattice function is a nice, clean, and elegant formulation that proceeded directly from the simpler 1D and 2D cases, but how do we implement such a formulation efficiently? Normally, if we were to simply find points on the curve or surface, we would utilize many efficient techniques that take advantage of discrete, uniformly sampled points in the domain (t, (u,v), and (s,t,u) points in the 1D, 2D and 3D cases respectively) such as the use of lookup tables for the sample values of the Bernstein polynomials, recursive subdivision to quickly find a point, and forward-differencing that takes advantage of a uniform stepsize. All of these techniques are much more efficient than the direct evaluation of the function; however, all but recursive subdivision, are not well-suited for finding a single point given an arbitrary point in the domain. In order to utilize recursive subdivision, we will have to gain some intuition on how these multivariate Bezier functions operate. This may be a bit tricky, but the effort will be worthwhile since the resulting algorithm will be very efficient and will not require the direct evaluation of any of the Bernstein polynomials.

Using the previously defined space curve function, PtOnCurve, but for 3D control points P_i, we can think of the PtInLattice function as a series of space curve evaluations. So, looking at the PtInLattice function, the innermost k-summation is evaluating the 3D space curve along the i,j row of control points for the parametric value u. The second j-summation is evaluating the ith 3D surface for the ith cross-section for the parametric domain coordinates (t,u). The outermost summation evaluates the blending of all of the (t,u) points on the cross-sections to give the single point (x,y,z) for the single parametric domain (s,t,u) coordinates. What this basically means is that we can think of the evaluation of the lattice point as being a blending of (t,u) points on the 3D surface cross-sections perpendicular to the X-axis. Each of these surface points come from an evaluation of the space curve point blending the t-points on the space curves in the j-direction perpendicular to the Y-axis. The points on these j-space curves result from evaluating the space curve defined by lattice points along the k direction for the given i,j. As we can clearly see here, it is a blending of space curve evaluations. So this means that if we write a single space curve recursive subdivision evaluation routine for a given set of 3D lattice control points, we can find a point in the 3D lattice.

If we let Nx, Ny, and Nz be the number of lattice control points in each direction (NumInX+1, NumInY+1, NumInZ+1 respectively), we can write the overall lattice evaluation routine as follows:

1. Evaluate the (t,u) points tuPt_i for each of the i cross-sections perpendicular to the X-axis, where 0<=i<Nx (there are Nx cross-sections).
2. Evaluate the final xyzPt as the resulting point from PtOnCurve(s, tuPt_i, Nx) where s is the parametric "position" along the curve, tuPt_i is the array of (t,u) points evaluated on each of the i cross-sections and Nx is the number of cross-sections.

We must now create a routine to evaluate the (t,u) point for the ith cross-section. Remember, the i cross-sections are perpendicular to the X-axis and are composed of an Ny x Nz lattice control point mesh. This is exactly like a normal 3D Bezier surface control mesh. We need only evaluate (t,u) on this surface as follows:

1. Evaluate the (u) points uPt_i for each of the j space curves perpendicular to X and Y-axes, where 0<=j<Ny (there are Ny space curves).
2. Evaluate the final tuPt as the resulting point from PtOnCurve(t, uPt_i, Ny) where t is parametric position along the curve, uPti is the array of u-points evaluated on each of the j space curves and Ny is the number of those space curves.

The two only remaining parts is to evaluate the u-points and to implement PtOnCurve; these two are actually one in the same. We can evaluate the u-points as follows:

1. For a particular ith cross-section, and a jth space curve, we have an an array of k-points along a particular i,j row. These k-points can be copied into the kPts_i array from LatticePoints[i][j][k] where k : 0<=k<Nz. The reasons for copying them will become apparent in our implementation of PtOnCurve.
2. Now to evaluate the u-point along this space curve, we need only call PtOnCurve(u, kPts_i, Nz).

Evaluation of a Point on the Space Curve using Recursive Subdivision: PtOnCurve implementation

We have developed a hierarchy of evaluations that are dependent only on one 3D Bezier space curve evaluation routine - PtOnCurve(). We will implement this routine using the DeCasteljau algorithm (refer to Foley & Van Dam for a proper treatment of this algorithm).

Given:
       t: parametric "position" along the curve on which to evaluate
     Pts: an array of 3D control points
  NumPts: the number of control points

Returns:
  Point on the curve evaluated at parametric position t.
  NOTE: Pts array is overwritten and Pts[0] is the resulting evaluated point.

Vec3f PtOnCurve(float t, Vec3f *Pts, int NumPts)
{
  for (int i=NumPts-1; i>0; i--)         // FOR EACH OF THE LEVELS OF RECURSION
    for (int j=0; j<i; j++)              // FOR EACH OF THE CNTRL PT PAIRS ON THE LEVEL
      Pts[j]=Pts[j]+(Pts[j+1]-Pts[j])*t; // EVALUATE IN-BETWEEN POINTS FOR EACH PAIR
  return(Pts[0]);                        // FIRST PT IN ARRAY BECOME EVALUATED PT
}

This results in an efficient implementation of PtOnCurve as opposed to the previous direct evaluation:

PtOnCurve(t,Pts,NumPts) = SUM_i=0,NumPts { Pts[i] * n!/(i!*(n-i)!) * tⁱ * (1-t)^n-i }

Evaluation of a Point on a 3D Bezier Surface: implementation of PtOnSurface for the ith cross section

The first step in our lattice point evaluation routine involves finding the (t,u) point on the ith cross-section. This amounts to nothing more than evaluating a point on a 3D Bezier surface defined by an N x O array of control points; the N direction refers to the t parametric direction, and O is the u direction. We will evaluate this (t,u) point using our PtOnCurve routine. In this case, the control points will NOT be overwritten:

Given:
  t,u: parametric 2D "position" in the domain of the surface on which to evaluate
  Pts: an NxO array of control points
  N,O: number of control points in the t and u directions respectively

Returns:
  Point on the surface evaluated at (t,u)

Vec3f PtOnSurface(float t, float u, Vec3f **Pts, int N, int O)
{
  Vec3f *jPts = new Vec3f[O];       // J-TH SPACE CURVE CTRL PTS TO EVAL IN u-DIR
  Vec3f *kPts = new Vec3f[N];       // RESULTING CTRL PTS TO EVAL IN t-DIRECTION
  for (int j=0; j<N; j++)           // FOR EACH J SPACE CURVE
  {
    for (int k=0; k<O; k++)         // COPY THE K-PTS
      jPts[k] = Pts[j][k];
    kPts[j] = PtOnCurve(u,jPts,O);  // EVAL THE K-PT ALONG THE JTH SPACE CURVE
  }
  Vec3f tuPt = PtOnCurve(t,kPts,N); // EVAL THE Pt on the kPts SPACE CURVE
  delete[] jPts;
  delete[] kPts;
  return (tuPt);
}

Evaluation of the Point in the Lattice: implementation of PtInLattice for an (s,t,u) point

Given:
           s,t,u: parametric 3D "position" or normalized coordinates inside the lattice structure
  LatticeCtrlPts: NxMxO 3D array of lattice control points
           N,M,O: number of control pts in the lattice in each axis direction

Returns:
  Point in the lattice at the position (s,t,u)

Vec3f PtInLattice(float s, float t, float u, Vec3f ***LatticeCtrlPts, int N, int M, int O)
{
  Vec3f *tuPts = new Vec3f[N];        // ith CROSS-SECTION (t,u) PTS
  for (int i=0; i<N; i++)             // FOR EACH i CROSS-SECTION
    tuPts[i] = PtOnSurface(t,u,LatticeCtrlPts[i],M,O); // EVAL (t,u) PT
  Vec3f xyzPt = PtOnCurve(s,tuPts,N); // EVAL FINAL LATTICE POINT
  delete[] tuPts;
  return (xyzPt);
}

Complete Final source Code

Source