Sign in to follow this  

Simpson's Rule in 3Space

Recommended Posts

Hi, I am implementing the following algorithm:

Point Y(float t);   //Theparameterizationofthecurve,a<=t<=b.
Point DY(float t);  //ThederivativeofY(t),a<=t<=b.
float Speed(float t) { return Length(DY(t)); }
float ArcLength(float t) {return Integral(a, t, Length(DY()); }

float GetT(floats)   //0<=s<=total_arc_length
   floatL=ArcLength(b);   //precomputeLifyoufrequentlycallGetT
   floatt=a+s*(b-a)/L;    //precompute(b-a)/LifyoufrequentlycallGetT
   for(inti=0;i<max;i++)  //‘max’isapplication-specified
      if(Abs(F)<epsilon)  //‘epsilon’isapplication-specified
         return t;

      t -= F/DF; // warning: DF = 0 might be a problem

   //maximum iterations exceeded, you might want to trap this
   return t;

/* The function Integral(a,t,f()) is any numerical integrator that computes the integral of f(t) over the interval[a,t].

I have tried this algorithm using simpson's rule as the numerical Integrator and have had no such luck. I believe this is because I do not have a good understanding of simpson's rule when doing vector math :). I have been using the following simpson's rule algorithm (coded in python I believe):
def simpson_rule(f, a, b):
  "Approximate the definite integral of f from a to b by Simpson's rule."
  c = (a + b) / 2.0
  h3 = abs(b - a) / 6.0
  return h3 * (f(a) + 4.0*f(c) + f(b))

# Calculates integral of sin(x) from 0 to 1
from math import sin
print simpson_rule(sin, 0, 1)

It appears very straight forward, here is my implementation of it in c++ using my function getWorldCoordsAt() instead of sin. getWorldCoordsAt() finds the world coordinates of the position on the line that corresponds to the normalized point passed in:
 * Finds the arc length to a given normalized point on the curve 
 * who's segment begins at the passed in control point. Thus we will 
 * measure u as being passed the location of the previous curve control
 * point, measuring as if the last control point is u. Thus getWorldCoordsAt()
 * is much like my f(x) for the curve.
 * @param prevCtrlPt The beginning of this segment of our overall line
 *                   it is an index into a list of points that I use as 
 *                   control points when calculating and drawing my curve.
 * @param u The normalized point on our curve that we are trying to find 
 *          the arc length to.
GLdouble Curve::arcLengthTo( int prevCtrlPt, GLdouble u ){
        // for me a (in the python algo above) is 0
        // and b is u, a value from 0 to 1
	GLdouble c = (u+0) / 2.0;
	GLdouble h = ABS(u-0) / 2.0;

        // gets the world coordinates at the 
        // distance of the normalized coordinate 0 to 1, 
        // on the curve segment that is anchored at the prevCtrlPt
        // curve control point. This is basically just running our 
        // point u through the "function" (matricies) that compute
        // our line points. 
        // so run the different values through our line function
	Pointd fa = getWorldCoordsAt( prevCtrlPt, 0 );
	Pointd fu = getWorldCoordsAt( prevCtrlPt, u );
	Pointd fc = getWorldCoordsAt( prevCtrlPt, c );

        // total up each of the function call values, f(?)
        // this is like adding up the inside set of (   ) 
	Pointd fT;
	fT.x = fa.x + 4*fc.x + fu.x;
	fT.y = fa.y + 4*fc.y + fu.y;
	fT.z = fa.z + 4*fc.z + fu.z;

        // multiply h through
	fT.x *= h;
	fT.y *= h;
	fT.z *= h;

        // divide by 3 through
	fT.x /= 3.0;
	fT.y /= 3.0;
	fT.z /= 3.0;

        // Since we are in 3Space i just took a wild guess
        // and found the length of the vector to see if that
        // is the correct value. It does not appear to be!
	GLdouble length = MatrixOps::vectorLength( fT );
        // The python algo:
	// return (h * ( f(0.0) + 4.0 * f(c) + f(u) ) / 3.0);
	return length;

So I really do not know whats wrong but im guessing it has something to do with using the simpsons algorithm in 3 space that i do not understand! ANY help would be appreciated! -sevans

Share this post

Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this