updated with latest N-D natural B-Spline fit model and example code

Author: superchromix

Gpufit/constants.h

@@ -20,7 +20,12 @@ enum ModelID {
     SPLINE_4D = 13,
     SPLINE_4D_MULTICHANNEL = 14,
     SPLINE_5D = 15,
-    NATURAL_BSPLINE_1D = 16
+    NATURAL_BSPLINE_1D = 16,
+    NATURAL_BSPLINE_2D = 17,
+    NATURAL_BSPLINE_3D = 18,
+    NATURAL_BSPLINE_4D = 19,
+    NATURAL_BSPLINE_5D = 20,
+    NATURAL_BSPLINE_6D = 21
 };
 
 // estimator ID

Gpufit/models/models.cuh

@@ -19,6 +19,7 @@
 #include "spline_4d_multichannel.cuh"
 #include "spline_5d.cuh"
 #include "natural_bspline_1d.cuh"
+#include "natural_bspline_nd.cuh"
 
 __device__ void calculate_model(
     ModelID const model_id,
@@ -86,6 +87,21 @@ __device__ void calculate_model(
     case NATURAL_BSPLINE_1D:
         calculate_natural_bspline1d(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size);
         break;
+    case NATURAL_BSPLINE_2D:
+        calculate_natural_bspline_nd(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size);
+        break;
+    case NATURAL_BSPLINE_3D:
+        calculate_natural_bspline_nd(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size);
+        break;
+    case NATURAL_BSPLINE_4D:
+        calculate_natural_bspline_nd(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size);
+        break;
+    case NATURAL_BSPLINE_5D:
+        calculate_natural_bspline_nd(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size);
+        break;
+    case NATURAL_BSPLINE_6D:
+        calculate_natural_bspline_nd(parameters, n_fits, n_points, value, derivative, point_index, fit_index, chunk_index, user_info, user_info_size);
+        break;
     default:
         assert(0); // unknown model ID
     }
@@ -112,6 +128,11 @@ void configure_model(ModelID const model_id, int & n_parameters, int & n_dimensi
     case SPLINE_4D_MULTICHANNEL:         n_parameters = 6; n_dimensions = 5; break;
     case SPLINE_5D:             n_parameters = 7; n_dimensions = 5; break;
     case NATURAL_BSPLINE_1D:    n_parameters = 3; n_dimensions = 1; break;
+    case NATURAL_BSPLINE_2D:    n_parameters = 4; n_dimensions = 2; break;
+    case NATURAL_BSPLINE_3D:    n_parameters = 5; n_dimensions = 3; break;
+    case NATURAL_BSPLINE_4D:    n_parameters = 6; n_dimensions = 4; break;
+    case NATURAL_BSPLINE_5D:    n_parameters = 7; n_dimensions = 5; break;
+    case NATURAL_BSPLINE_6D:    n_parameters = 8; n_dimensions = 6; break;
     default: throw std::runtime_error("unknown model ID");
     }
 }

Gpufit/models/natural_bspline_1d.cuh

@@ -18,7 +18,7 @@
 *
 * parameters: An input vector of concatenated sets of model parameters.
 *             p[0]: amplitude
-*             p[1]: center coordinate
+*             p[1]: x-shift
 *             p[2]: offset
 *
 * n_fits: The number of fits.
@@ -60,30 +60,40 @@ __device__ void calculate_natural_bspline1d(
     REAL const * ui = (REAL*)user_info;
 
     int const num_coeff = static_cast<int>(ui[0]);
+    int const N_tpl = num_coeff - 2;   // template length
     int const num_knots = num_coeff + 4; // cubic
     REAL const * knots = ui + 1;
     REAL const * coeff = ui + 1 + num_knots;
 
-    // Model parameters: [amp, center, offset]
+    // Model parameters: [amp, shift, offset]
     REAL const * p = parameters;
     REAL amp    = p[0];
-    REAL center = p[1];
+    REAL shift  = p[1];
     REAL offset = p[2];
 
     // Data point coordinate
     REAL x = static_cast<REAL>(point_index);
-    REAL xq = x - center;
+    REAL xq = x - p[1]; 
+
+    bool clamped = false;
+    if (xq < (REAL)0)           { xq = (REAL)0;           clamped = true; }
+    if (xq > (REAL)(N_tpl - 1)) { xq = (REAL)(N_tpl - 1); clamped = true; }
 
     // Find knot span as in host code
     const int k = 3;
-    int N = num_coeff - 2;
     int span;
-    if (xq <= 0.0)
+    if (xq <= (REAL)0) 
+    {               
         span = k;
-    else if (xq >= REAL(N - 1))
+    }
+    else if (xq >= (REAL)(N_tpl - 1))
+    {
         span = num_coeff - 1;
+    }
     else
-        span = int(xq) + k;
+    {
+        span = (int)xq + k;
+    }
 
     // Evaluate basis and derivative (device function!)
     REAL basis[4], dbasis[4];
@@ -102,13 +112,15 @@ __device__ void calculate_natural_bspline1d(
         }
     }
 
+    if (clamped) spline_dx = (REAL)0;
+
     // Write value
     value[point_index] = amp * spline_val + offset;
 
-    // Write derivatives [amp, center, offset]
+    // Write derivatives [amp, shift, offset]
     REAL * der = derivative + point_index;
     der[0 * n_points] = spline_val;         // d/d(amp)
-    der[1 * n_points] = -amp * spline_dx;   // d/d(center)
+    der[1 * n_points] = -amp * spline_dx;   // d/d(shift)
     der[2 * n_points] = 1;                  // d/d(offset)
 }