No Pity for the Masses

Free form deformation(FFD) is a fun little technique in animation allowing for region specific deformations.  You specify a series of control points around the surface of your object.  Each control point asserts a certain level of regional influence along the object.  That influence is directly related to the control point’s position in the grid and the number of contending points nearby.  (If that sounds convoluted just look at the link I have above.)

In this image we see the brain model is being pulled by a single control point. Notice how only a small part of the brain is modified. This is due to the other control points remaining in their original positions.

Computing FFD’s on the fly can be computationally intensive.   The most common FFD approach was developed by Sederberg and Parry in their paper: “Free-Form Deformation of Solid Geometric Model.” Basically you follow a pretty simple series of steps:

1. Define a series of axes: S,T,U.  These don’t need to be orthogonal but in my example they are.  I also go one step further and align these axes with X,Y, and Z but scaled to the size of my object.
2. Redefine your point P in terms of your S,T, and U axes and an origin.  The origin in my system is the minimum point on my mesh.
3. Create a grid of control points.  The easiest way is to create a uniform series of points along S,T, and U.
4. Solve the vector valued  trivariate bernstein polynomial shown in the Sederberg/Parry paper
5. Draw new point

Code for my project can found here. SVN login instructions can found to your left in my page rollout. (Code will not be posted until 2/30/10 as this project has not been turned in yet.)

But let me point a few things out:

void tyrFreeFormDef::reParamVertices(vector<RectCoord> &param, vector3 &S,vector3 &T, vector3 &U){
	printf("  - Reparameterizing vertices\n");
 
	int nVtx = _ply.GetNumVertices();
	point3D_t *vtx = _ply.GetVertices();
 
	vector3 min = _ply.GetMin();
	vector3 p0 = min;
	vector3 max = _ply.GetMax();
 
	S = vector3(max.x - min.x,0.0,0.0);
	T = vector3(0.0,max.y - min.y,0.0);
	U = vector3(0.0,0.0,max.z - min.z);
 
 
	vector3 TcU = T.Cross(U);
	vector3 ScU = S.Cross(U);
	vector3 ScT = S.Cross(T);
 
	double TcUdS = TcU.Dot(S);
	double ScUdT = ScU.Dot(T);
	double ScTdU = ScT.Dot(U);
 
	for(int v = 0; v < nVtx; v++){
		vector3 diff = vtx[v].ToVector3() - p0;
 
		RectCoord tmp;
		tmp.s = TcU.Dot(diff/TcUdS);
		tmp.t = ScU.Dot(diff/ScUdT);
		tmp.u = ScT.Dot(diff/ScTdU);
		tmp.p = p0 + (tmp.s * S) + (tmp.t * T) + (tmp.u * U);
		tmp.p0 = p0;
 
		{	///Pre-calculate bernstein polynomial expansion.  It only needs to be done once per parameterization
			for(int i = 0; i <= 5; i++){
				for(int j = 0; j <= 5; j++){
 
					for(int k = 0; k <= 5; k++){
						tmp.bernPolyPack[n][k][2] = bern_poly(n,k,tmp.u);
					}
					tmp.bernPolyPack[m][j][1] = bern_poly(m,j,tmp.t);
				}
				tmp.bernPolyPack[l][i][0] = bern_poly(l,i,tmp.s);
			}
		}
		param.push_back(tmp);
		if(tmp.p.Dist(vtx[v].ToVector3()) > Math::EPSILON){
			printf("     - Warning, vtx[%d] does not match it's parameterization.\n",v);
		}
	}
	for(int i = 0; i <= 5; i++){
		for(int j = 0; j <= 5; j++){
			vector3 tK(0.0f,0.0f,0.0f);
			for(int k = 0; k <= 5; k++){
				controlPoints[i][j][k] = createControlPoint(min,i,j,k,S,T,U);
 
			}
		}
	}
}

This is my reparameterization function relating point P to an origin and distance along axes S/T/U. Notice the “pre-calculate bernstein polynomial…” section. Here each new parameterized point is also given a table of lookup bernstein polynomial values. At a cost in memory we save on recalculation later on. This becomes important if we’re trying to render our object in real-time.

double tyrFreeFormDef::binom_coeff4(int k){
	///Manually calculated binomial coefficients
	if(k==0||k==4)return 1.0f;
	else if(k==1||k==3)return 4.0f;
	else if(k==2)return 6.0f;
	else return 0.0f;
}	
double tyrFreeFormDef::bern_poly(int n, int v, double x){
	///bernstein polynomial expansion (currently only supports marker size of 4)
	return binom_coeff4(v) * pow(x,(double)v) *  pow((1.0f-x),(double)(n-v));
}

Bernstein polynomial expansion and binomial coefficient functions. Make note of the binomial coef. function. Before I implemented the bernstein poly. packing I was calculating binomial coef’s every frame, to speed up the process I wrote a simply function that assumes the total # of control points on each axis is 4. This function could easily be replaced with one that actually does a factorial expansion.

vector3 tyrFreeFormDef::getGlobalVertice(RectCoord &r,vector3 &S,vector3 &T,vector3 &U){
	double s = r.s;
	double t = r.t;
	double u = r.u;
 
	vector3 tS(0.0f,0.0f,0.0f);
	for(int i = 0; i <= l; i++){
 
		vector3 tM(0.0f,0.0f,0.0f);
		for(int j = 0; j <= m; j++){
 
			vector3 tK(0.0f,0.0f,0.0f);
			for(int k = 0; k <= n; k++){
				tK += r.bernPolyPack[n][k][2] * controlPoints[i][j][k];
			}
			tM += r.bernPolyPack[m][j][1] * tK;
		}
		tS += r.bernPolyPack[l][i][0] * tM;
	}
	return tS;
}

Last is the conversion function. This is where we apply our interpolations to the actual mesh. The return value is in global coordinate space. In other words the point coming out of this function is what you would display.