centripetal parametrization for curve fitting or interpolation

very useful to compute parameter for curve fitting or interpolation
the detail explanation is written here


and here is the code for three dimensional parametric curve

 void GetCentripetal_Param(  
       int ncp,                                                  //number of control point  
       double point[][3],                                        //control points of curve  
       double t[]                                                //array of t parameter (output) 
       )  
 {  
      double *dcp = (double*)malloc(sizeof(double) * (ncp-1));   //CP0-CP1, CP2-CP1, ..., CP[ncp-1][ncp]  
      double *sum = (double*)malloc(sizeof(double) * ncp);  
             
      int i, j;  
      double dcptotal = 0.0;  
      for(i=0; i<ncp-1; i++)  
           dcp[i] = 0.0;  
        
      sum[0] = 0.0;  
      for(i=0; i<ncp-1; i++)  
      {  
        for(j=0; j<3; j++)  
          dcp[i] += (point[i+1][j]-point[i][j]) * (point[i+1][j]-point[i][j]);        
       
        dcp[i] = sqrt(sqrt(dcp[i]));  
        dcptotal += dcp[i];  
        sum[i+1] = dcptotal;  
      }  
      sum[0] = 0.0;  
      double inv_dcptotal = 1 / dcptotal;  
      for(i=0; i<ncp; i++)  
         t[i] = (sum[i] * inv_dcptotal);  
        
      if(dcp) free(dcp);  
      if(sum) free(sum);  
 }     

comparison of secant and bisection method

   
 #include<stdio.h>  
 #include<math.h>  
 #include<time.h>  
   
 //example of function f(t), you can make any function you like  
 float f(float t)  
 {  
   return(((48.2291-30.6208*t-43) * (48.2291-30.6208*t-43)) + 35*t - 25.0);   
 }  
   
 void main()  
 {  
     
   float a = 0.1, b = 0.3, c = 0, e = 7.5e-3;  
   //a, b interval  
   //e error  
   
   int count = 1, n = 50; //n=maximum iteration  
     
   //  
   //secant  
   //  
   clock_t begin = clock();  
   do  
   {      
     c = (a*f(b) - b*f(a)) / (f(b) - f(a));  
     //c = b - (f(b)*(b-a) / (f(b) - f(a)));  
     a=b;  
     b=c;  
     printf("i No-%d  t=%f  d=%f\n",count,c, f(c));  
     count++;  
     if(count==n) break;  
   } while(fabs(f(c))&gt;e);  
   
   clock_t end = clock();  
   double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;  
   printf("time:%f", time_spent);  
   
   printf("\nThe required solution is %f\n",c);  
   printf("d = %f\n", f(c) + 25);  
   
   //  
   //bisection  
   //  
   double fa, fb, fc;  
   a = 0.1; b = 0.3; c = 0; e = 7.5e-3;  
   count = 1;  
   clock_t begin2 = clock();  
   do{  
     c = (a+b)/2;  
     fa = f(a);  
     fb = f(b);  
     fc = f(c);  
   
     if(fa*fc &lt; 0){  
       b = c;  
     }  
     else{  
       a = c;  
     }  
     printf("i No-%d  t=%f  d=%f\n", count, c, fc);  
     count++;  
     if(count==n) break;  
   } while(fabs(fc)&gt;e);  
   
   clock_t end2 = clock();  
   double time_spent2 = (double)(end2 - begin2) / CLOCKS_PER_SEC;  
   printf("time:%f", time_spent2);  
   
   printf("\nThe required solution is %f\n",c);  
   printf("d = %f\n", fc + 25);  
   
    
 }