Creating a smooth surface from an array of points

jmurphy

Well-known member
There are two JavaScript problems I would like to solve:
  1. Create a smooth curve from a 1D array of points, with the ability to control how much smoothing is applied
  2. Create a smooth surface from a 2D array of points, with the ability to control how much smoothing is applied
Looking at the Object Explorer, it looks like there are some powerful tools available, such as:
  • SurfaceSpline
  • SurfaceSimplifier
  • ShepardInterpolation
  • GridInterpolation
Are any of these suitable, and if so, how are they used?

My current solution for creating a smooth curve for a 1D array of points works, but it is a bit clunky. This is output form the PhotometricMosaic version 2.0 which I have almost finished (currently updating the help file and tidying up the code).
GradientGraph.PNG


Thanks
John Murphy
 
Hi John,

For the first problem you can use the following JavaScript implementation of the Akima subspline interpolation algorithm. This is a JS version of the equivalent C++ class that forms part of PCL:

JavaScript:
/*
* Akima subspline interpolation algorithm.
*
* References
* ----------
*
* Hiroshi Akima, A new method of interpolation and smooth curve fitting based
* on local procedures, Journal of the ACM, Vol. 17, No. 4, October 1970,
* pages 589-602.
*
* Implementation
* --------------
*
* Our implementation is based on the book:
*
* Numerical Algorithms with C, by G. Engeln-Mullges and F. Uhlig (Springer,
* 1996), section 13.1.
*
* We properly represent corners when a data point lies between two adjacent
* straight lines with different slopes. This means that our implementation
* does not impose continuous differentiability, which deviates from the
* original work by Akima. Supporting the accurate representation of corners
* has several practical advantages in our opinion; one of them is the enhanced
* flexibility for the application of Akima interpolation to graphical
* representations of curves given by a set of prescribed x,y data points.
*
* Parameters
* ----------
*
* x  Vector of x-values:
*
*    - If x is not empty: Must be a vector of monotonically increasing,
*      distinct values: x[0] < x[1] < ... < x[n-1].
*
*    - If x is empty: The interpolation will use implicit x[i] = i for
*      i = {0,1,...,n-1}.
*
* y  Vector of function values for i = {0,1,...,n-1}.
*
* The length of the y vector (and also the length of a nonempty x vector) must
* be n >= 5. This is because the Akima algorithm requires at least 4
* interpolation subintervals.
*/
function AkimaInterpolation( x, y )
{
   if ( y.length < 5 )
      throw new Error( "AkimaInterpolation(): Less than five data points specified." );
   if ( x.length != 0 && x.length < y.length )
      throw new Error( "AkimaInterpolation(): Invalid x vector length." );

   this.x = x;
   this.y = y;

   let n = this.y.length;
   let N = n-1; // number of subintervals

   this.b = new Float32Array( N );
   this.c = new Float32Array( N );
   this.d = new Float32Array( N );

   // Chordal slopes. We need room for 4 additional prescribed slopes.
   let m = new Float32Array( N+4 ); //

   // Akima left-hand slopes to support corners.
   let tL = new Float32Array( n );

   // Calculate chordal slopes for each subinterval.
   if ( this.x.length > 0 )
      for ( let i = 0; i < N; ++i )
      {
         let h = this.x[i+1] - this.x[i];
         if ( 1 + h*h == 1 )
            throw new Error( "AkimaInterpolation(): Empty interpolation subinterval(s)." );
         m[i+2] = (this.y[i+1] - this.y[i])/h;
      }
   else
      for ( let i = 0; i < N; ++i )
         m[i+2] = this.y[i+1] - this.y[i];

   // Prescribed slopes at ending locations.
   m[  0] = 3*m[  2] - 2*m[3];
   m[  1] = 2*m[  2] -   m[3];
   m[N+2] = 2*m[N+1] -   m[N];
   m[N+3] = 3*m[N+1] - 2*m[N];

   /*
    * Akima left-hand and right-hand slopes.
    * Right-hand slopes are just the interpolation coefficients b[i].
    */
   for ( let i = 0; i < n; ++i )
   {
      let f = Math.abs( m[i+1] - m[i  ] );
      let e = Math.abs( m[i+3] - m[i+2] ) + f;
      if ( 1 + e != 1 )
      {
         tL[i] = m[i+1] + f*(m[i+2] - m[i+1])/e;
         if ( i < N )
            this.b[i] = tL[i];
      }
      else
      {
         tL[i] = m[i+1];
         if ( i < N )
            this.b[i] = m[i+2];
      }
   }

   /*
    * Interpolation coefficients b[i], c[i], d[i].
    * a[i] coefficients are the y[i] ordinate values.
    */
   for ( let i = 0; i < N; ++i )
   {
      this.c[i] = 3*m[i+2] - 2*this.b[i] - tL[i+1];
      this.d[i] = this.b[i] + tL[i+1] - 2*m[i+2];
      if ( this.x.length > 0 )
      {
         let h = this.x[i+1] - this.x[i];
         this.c[i] /= h;
         this.d[i] /= h*h;
      }
   }

   this.evaluate = function( x )
   {
      /*
       * Find the subinterval i0 such that this.x[i0] <= x < this.x[i0+1].
       * Find the distance dx = x - this.x[i], or dx = x - i for implicit
       * x = {0,1,...n-1}.
       */
      let i0;
      let dx;
      if ( this.x.length > 0 )
      {
         i0 = 0;
         let i1 = this.x.length - 1;
         while ( i1-i0 > 1 )
         {
            let im = (i0 + i1) >> 1;
            if ( x < this.x[im] )
               i1 = im;
            else
               i0 = im;
         }
         dx = x - this.x[i0];
      }
      else
      {
         if ( x <= 0 )
            return this.y[0];
         if ( x >= this.y.length-1 )
            return this.y[this.y.length-1];
         i0 = Math.trunc( x );
         dx = x - i0;
      }

      /*
       * Use a Horner scheme to calculate b[i]*dx + c[i]*dx^2 + d[i]*dx^3.
       */
      return this.y[i0] + dx*(this.b[i0] + dx*(this.c[i0] + dx*this.d[i0]));
   };
}

Akima is the default interpolation algorithm used by the CurvesTransformation tool. I have used the above JS implementation in different occasions with good results. Its main advantages, for example with respect to cubic interpolation, is its much higher resilience to oscillations and its ability to represent corners accurately, i.e. points where the function is not differentiable. This is very nice for interactive curve manipulation, or for representations of functions with strong variations.

For the smooth surface representation problem, if you need accuracy and control on the degree of interpolation smoothness, then you need thin plate splines, also known as two-dimensional surface splines. The core PJSR object SurfaceSpline is functionally equivalent to the C++ template class SurfaceSpline<>, which is fully documented here:


I'll post an example with the SurfaceSpline JavaScript object later, but the above documentation provides you with everything you need to use it. Let me know if you need more information.
 
  • Like
Reactions: dld
Thanks Juan
I have updated my code to use Akima Interpolation, and it works really well. A big improvement!
I will have a look at SurfaceSpline over the next few days.
Thanks
John
 
Back
Top