Compute the integral numerically using something like an adaptive version of Simpson's rule. That's probably what I would do.

Here is what I am currently using:

void ApproximateCatMullRomLength3D(double *Length,const D3DXVECTOR3* p0, const D3DXVECTOR3* p1, const D3DXVECTOR3* p2, const D3DXVECTOR3* p3){
	D3DXVECTOR3 v0,v1;
	v0=*p1;
	*Length=0;
	for (int t=1;t<1000;t++){
		float time=t/1000.0f;
		CatmullRom_Evaluate(&v1,p0,p1,p2,p3,time);//a simple CMR function
		*Length+=double(D3DXVec3Length(&(v1-v0)));
		v0=v1;
	}
}

All it does is evaluate 1000 points along the curve and sums up the lengths between each one..... I may not need that many samples, but I'm using 1000 until I get something that works right in my AutoPilot function.

Ok. Picture time:

This illustrates what I am trying to do (as best as I can draw it). I originally have all the "waypoints" calculated for avoidance of an object represented here as "start", "wp1", "wp2", and "end". I then calculate the CMR points by shooting vectors from the WP to its neighbors (or start/end). I'm doing the calc for WP1 in this illustration..... The inner control points (p1 and p2) are determined by some max half-way length and placed in-line between their respective WPs. The outer control points (p0 and p3) are simply shot out using the same vectors but the length is a multiple of the "max half-way length" I determined earlier. The multiple I use gives varying results of sometimes strange curves. I think I'm doing it wrong..... You may be thinking I should just use the waypoints as the control points, but that won't work in most cases for my application......

I think I didn't explain my second image correctly. There are 4 separate examples in it. Each one is a different example of the multiplier I use. Here is the code. I have since modified it and it works a little better with the mod:

(In the code "p#" is "CV#" representing the control points)

void AUTOPILOT::CalculateCurves(SHIP *Ship){
	if (WP.size()<2) return;//no need to do the function
	for (UINT wp=0;wp<WP.size()-1;wp++){
		//find the max distance. max is shortest WP away from WP[wp]
		D3DXVECTOR3 dir1,dir2;
		double dist;
		if (wp==0){//wp-1 is the ship's location
			HUGEVECTOR3 hvDir1=Ship->V_Location-WP[wp].vWP;
			dir1=hvDir1.Normalize();
			dist=hvDir1.Length/2.50;
		}
		else {
			HUGEVECTOR3 hvDir1=WP[wp-1].vWP-WP[wp].vWP;
			dir1=hvDir1.Normalize();
			dist=hvDir1.Length/2.50;
		}
		HUGEVECTOR3 hvDir2=WP[wp+1].vWP-WP[wp].vWP;
		dir2=hvDir2.Normalize();
		double dist2=hvDir2.Length/2.0;
		if (dist2<dist) dist=dist2;
		if (dist>WP[wp].MaxAllowedVel*10.0) dist=WP[wp].MaxAllowedVel*10.0;

		WP[wp].CV1=dir1;
		WP[wp].CV1*=float(dist);
		WP[wp].CV0=dir1*6.0f*float(dist);
		WP[wp].CV2=dir2;
		WP[wp].CV2*=float(dist);
		WP[wp].CV3=dir2*6.0f*float(dist);//this one WAS *8.0
              //***************
              //this is part of the modification
		dir1=WP[wp].CV0-WP[wp].CV3;
		dir1/=6.0f;
		WP[wp].CV0-=dir1;
		WP[wp].CV3+=dir1;
              //********************

		ApproximateCatMullRomLength3D(&WP[wp].CurveSize,&WP[wp].CV0,&WP[wp].CV1,&WP[wp].CV2,&WP[wp].CV3);
	}
}

I annotated where I changed the position of CV0 & CV3 and my "multiplier" is x6. This gives me a more gradual curve between CV1 &CV2.

.