core/vnl/algo/vnl_adaptsimpson_integral.cxx
Go to the documentation of this file.
00001 #include "vnl_adaptsimpson_integral.h"
00002 #include <vcl_iostream.h>
00003 #include <vcl_cmath.h>
00004 
00005 double vnl_adaptsimpson_integral::int_fnct_(double* x)
00006 {
00007   return  pfnct_->f_(*x);
00008 }
00009 
00010 double vnl_adaptsimpson_integral::integral(vnl_integrant_fnct* f, double a,
00011                                            double b, double acury)
00012 {
00013   //set the function
00014   pfnct_ = f;
00015 
00016   return adaptivesimpson(&vnl_adaptsimpson_integral::int_fnct_, a, b, acury, 0, depth_);
00017 }
00018 
00019 double vnl_adaptsimpson_integral::adaptivesimpson(double(*f)(double*),
00020                                                   double a, double b, double eps,
00021                                                   int level, int level_max)
00022 {
00023   double c, d, e, h, result;
00024   double one_simpson, two_simpson;
00025   double left_simpson, right_simpson;
00026 
00027   h = b-a;
00028   c = 0.5*(a+b);
00029   one_simpson = h*(f(&a)+4.0*f(&c)+f(&b))/6.0;
00030   d = 0.5*(a+c);
00031   e = 0.5*(c+b);
00032   two_simpson = h*(f(&a)+4.0*f(&d)+2.0*f(&c)+4.0*f(&e)+f(&b))/12.0;
00033   /* Check for level */
00034   if (level+1 >= level_max) {
00035     result = two_simpson;
00036     vcl_cerr<< "Maximum level reached\n";
00037   }
00038   else {
00039     /* Check for desired accuracy */
00040     if (vcl_fabs(two_simpson-one_simpson) < 15.0*eps)
00041       result = two_simpson + (two_simpson-one_simpson)/15.0;
00042     /* Divide further */
00043     else {
00044       left_simpson = adaptivesimpson(f,a,c,eps/2.0,level+1,level_max);
00045       right_simpson = adaptivesimpson(f,c,b,eps/2.0,level+1,level_max);
00046       result =  left_simpson + right_simpson;
00047     }
00048   }
00049   return result;
00050 }
00051