Sign in to follow this  

Simpson's Rule in 3Space

This topic is 4056 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

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
   {
      floatF=ArcLength(t)-s;
      if(Abs(F)<epsilon)  //‘epsilon’isapplication-specified
         return t;

      floatDF=Speed(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

This topic is 4056 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

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