Jump to content
  • Advertisement
Sign in to follow this  
Sevans

Simpson's Rule in 3Space

This topic is 4275 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
Advertisement
Sign in to follow this  

  • Advertisement
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

We are the game development community.

Whether you are an indie, hobbyist, AAA developer, or just trying to learn, GameDev.net is the place for you to learn, share, and connect with the games industry. Learn more About Us or sign up!

Sign me up!