fpParametricCurve: improvements (#24)

Author: perazz

include/fitpack_core_c.h

@@ -245,8 +245,24 @@ FP_FLAG parder_c(const FP_REAL* tx, FP_SIZE nx, const FP_REAL* ty, FP_SIZE ny, c
 
 #endif
 
+// Core C++ functions
 #ifdef __cplusplus
-    typedef vector<FP_REAL> fpPoint;
+
+typedef vector<FP_REAL> fpPoint;
+    
+// Flatten a 2D vector
+inline vector<FP_REAL> flatten_2d_vector(const vector<fpPoint>& x2d) {
+    FP_SIZE ndim = 0;
+    for (const auto& p : x2d) ndim = std::max(ndim, static_cast<FP_SIZE>(p.size()));
+
+    vector<FP_REAL> x(x2d.size() * ndim, 0.0);
+    for (FP_SIZE i = 0; i < x2d.size(); i++)
+        copy(x2d[i].begin(), x2d[i].end(), x.begin() + i * ndim);
+
+    return x;
+}
+
+    
 #endif
 
 #endif // FITPACK_CORE_C_H_INCLUDED

include/fitpack_parametric_curves_c.h

@@ -50,7 +50,7 @@ void    fitpack_parametric_curve_c_new_points            (fitpack_parametric_cur
 void    fitpack_parametric_curve_c_set_default_parameters(fitpack_parametric_curve_c *self);
 FP_FLAG fitpack_parametric_curve_c_new_fit               (fitpack_parametric_curve_c *self, FP_SIZE ndim, FP_SIZE npts, FP_REAL* x, FP_REAL* u, FP_REAL* w, FP_REAL* smoothing, FP_SIZE* order);
 FP_FLAG fitpack_parametric_curve_c_fit                   (fitpack_parametric_curve_c *self, FP_REAL* smoothing, FP_SIZE* order);
-FP_FLAG fitpack_parametric_curve_c_interpolating         (fitpack_parametric_curve_c *self);
+FP_FLAG fitpack_parametric_curve_c_interpolating         (fitpack_parametric_curve_c *self, FP_SIZE* order);
 FP_FLAG fitpack_parametric_curve_c_eval_one              (fitpack_parametric_curve_c *self, FP_REAL u, FP_REAL* y);
 //FP_REAL fitpack_parametric_curve_c_integral       (fitpack_parametric_curve_c *self, FP_REAL from, FP_REAL to);
 //void    fitpack_parametric_curve_c_fourier        (fitpack_parametric_curve_c *self, FP_SIZE nparm, const FP_REAL* alpha, FP_REAL* A, FP_REAL* B, FP_FLAG* ierr);

include/fpParametricCurve.hpp

@@ -29,37 +29,6 @@
 #include <vector>
 using std::vector;
 
-// Flatten a 2D vector
-vector<FP_REAL> flatten_2d_vector(vector<fpPoint> x2d)
-{
-
-    // Get number of points
-    FP_SIZE npts = x2d.size();
-
-    // Get max vector size
-    FP_SIZE ndim = 0;
-    for (FP_SIZE i=0; i<x2d.size(); i++)
-    {
-        ndim = x2d[i].size()>ndim? x2d[i].size() : ndim;
-    }
-
-    // Allocate flattened vector
-    FP_SIZE n = npts*ndim;
-
-    vector<FP_REAL> x(n,0.0);
-
-    for (FP_SIZE pt=0; pt<x2d.size(); pt++)
-    {
-        for (FP_SIZE idim=0; idim<x2d[pt].size(); idim++)
-        {
-            x[pt*ndim+idim] = x2d[pt][idim];
-        }
-    }
-
-    return x;
-
-}
-
 class fpParametricCurve
 {
     public:
@@ -110,7 +79,8 @@ class fpParametricCurve
         FP_FLAG fit(FP_REAL smoothing, FP_SIZE order) { return fitpack_parametric_curve_c_fit(&cptr,&smoothing,&order); }
 
         // Get the interpolating fit
-        FP_FLAG interpolate()                         { return fitpack_parametric_curve_c_interpolating(&cptr); }
+        FP_FLAG interpolate()                         { return fitpack_parametric_curve_c_interpolating(&cptr,nullptr); }
+        FP_FLAG interpolate(FP_SIZE order)            { return fitpack_parametric_curve_c_interpolating(&cptr,nullptr); }
 
         // Fit properties
         const FP_SIZE degree   () { return fitpack_parametric_curve_c_degree(&cptr); };
@@ -148,6 +118,21 @@ class fpParametricCurve
 
         }
 
+        // Get values at a range of parameter values
+        vector<fpPoint> eval(FP_REAL umin, FP_REAL umax, FP_SIZE npts = 100, FP_SIZE* ierr=nullptr)
+        {
+           vector<FP_REAL> uu(npts);
+           
+           // Generate linearly spaced u values
+           FP_REAL du = (umax - umin) / (npts - 1);
+           for (FP_SIZE i = 0; i < npts; ++i) {
+               uu[i] = umin + i * du;
+           }
+           
+           return eval(uu,ierr);
+           
+        }             
+        
         // Get single derivative at u
         fpPoint ddu(FP_REAL u, FP_SIZE order, FP_FLAG* ierr=nullptr)
         {

src/fitpack_parametric.f90

@@ -437,11 +437,15 @@ function curve_eval_many(this,u,ierr) result(x)
     end function curve_eval_many
 
     ! Interpolating curve
-    integer function interpolating_curve(this) result(ierr)
+    integer function interpolating_curve(this,order) result(ierr)
         class(fitpack_parametric_curve), intent(inout) :: this
+        integer(FP_SIZE), optional, intent(in) :: order
+
+        ! Set zero smoothing
+        this%iopt = IOPT_NEW_SMOOTHING
 
         ! Set zero smoothing
-        ierr = curve_fit_automatic_knots(this,zero)
+        ierr = curve_fit_automatic_knots(this,smoothing=zero,order=order)
 
     end function interpolating_curve
 

src/fitpack_parametric_c.f90

@@ -222,16 +222,17 @@ integer(FP_FLAG) function fitpack_parametric_curve_c_fit(this,smoothing,order) r
      end function fitpack_parametric_curve_c_fit
 
      !> Wrapper to interpolating_curve
-     integer(FP_FLAG) function fitpack_parametric_curve_c_interpolating(this) result(ierr) &
+     integer(FP_FLAG) function fitpack_parametric_curve_c_interpolating(this,order) result(ierr) &
                                 bind(c,name='fitpack_parametric_curve_c_interpolating')
         type(fitpack_parametric_curve_c), intent(inout) :: this
+        integer(FP_SIZE), optional, intent(in) :: order
 
         type(fitpack_parametric_curve), pointer :: fcurve
 
         !> Get object; allocate it in case
         call fitpack_parametric_curve_c_pointer(this, fcurve)
 
-        ierr = fcurve%interpolate()
+        ierr = fcurve%interpolate(order)
 
      end function fitpack_parametric_curve_c_interpolating
 

test/fitpack_tests.f90

@@ -2189,16 +2189,16 @@ logical function mnparc(iunit) result(success)
           end do approximations
 
          ! Format statements
-         910  format(31h0least-squares curve of degree ,i1,7h  ipar=,i1)
-         915  format(27h0smoothing curve of degree ,i1,7h  ipar=,i1)
-         920  format(20h smoothing factor s=,f7.2)
-         925  format(1x,23hsum squared residuals =,e15.6,5x,11herror flag=,i2)
-         930  format(1x,24htotal number of knots n=,i3)
-         935  format(1x,22hposition of the knots )
+         910  format('0least-squares curve of degree ',i1,'  ipar=',i1)
+         915  format('0smoothing     curve of degree ',i1,'  ipar=',i1)
+         920  format(' smoothing factor s=',f7.2)
+         925  format(1x,'sum squared residuals =',e15.6,5x,'error flag=',i2)
+         930  format(1x,'htotal number of knots n=',i3)
+         935  format(1x,'position of the knots ')
          940  format(5x,10f6.0)
-         945  format(1x,30hb-spline coefficients of sx(u))
+         945  format(1x,'b-spline coefficients of sx(u)')
          950  format(5x,8f8.4)
-         955  format(1x,30hb-spline coefficients of sy(u))
+         955  format(1x,'b-spline coefficients of sy(u)')         
          960  format(1h0,2(4x,'xi',7x,'yi',6x,'sx(ui)',3x,'sy(ui)'))
          965  format(1x,8f9.4)
         1000  format('[mnparc] parametric curve test ',i0,' failed: ',a)