// .----------------------------------------------------------------------.
// |  Nonuniform rational B-spline (NURBS) macro for povray       v 0.01  |
// |                                                                      |
// |  calculates a point on a nurbs curve of any order                    |
// |                                                                      |
// |  2006 J.W. Stolk                       http://www.jwstolk.xs4all.nl  |
// |                                                                      |
// |  based on:                                                           |
// |  "An Introduction to NURBS" by David F. Rogers. (2000)               |
// |  http://www.nar-associates.com/nurbs/c_code.html                     |
// '----------------------------------------------------------------------'

// run:
// megapov nurbs.pov +w10 +h10 -a -j -F

// .----------------------------------------------------------------------.
// |  Why? / Usage:                                                       |
// '----------------------------------------------------------------------'
// 
// I use this to for animations. The object, camera and camera_look_at points
// points can be drawn and edited in Rhino 3D, and used in pov-ray.
//
// these 3 macros need to be called:
//
//    macro_nurbs_initialize()
//  
//    #declare nurbs_max_t = macro_nurbs_get_max_t();
// 
//    macro_nurbs_get_point(0.0 .. nurbs_max_t)
//    (result is <nurbs_p[0],nurbs_p[1],nurbs_p[2]> )
// 
// notes:
// - I used megapov 1.2.1, but it should work with the official povray
// - this macro only works for one nurbs curve at a time. (it uses global variables)
// - the examples at http://www.nar-associates.com/nurbs/c_code.html do not use
//   the first element of each array ( [0] ). but this macro does!
// - this macro can only sample nurbs, the nurbs can't be used like regular
//   splines in pov-ray. But you sample the nurbs curve at many points and
//   use those points to define a bicubic spline, you will get close :-) 
// - Rhino 3D uses "degree" which is one less than the nurbs_order used here!
//   "degree=2 in Rhino" == "nurbs_order=3 in this macro"
//   "degree=3 in Rhino" == "nurbs_order=4 in this macro", etc.
// - this macro should be able to do high orders (but it will become slow)
// - "export selected, .pov" does not work for curves in Rhino. Use 
//   "export selected, .igs" and get the points using a text editor.

#declare nurbs_vertexes = 5;
#declare nurbs_order    = 4;

// array containing the defining polygon vertices
#declare nurbs_vert_array = array[nurbs_vertexes*3] {
//x:   y:   z:  
  0.0, 0.0, 1.0,
  1.0, 2.0, 1.0,
  2.5, 0.0, 1.0,
  4.0, 2.0, 1.0,
  5.0, 0.0, 1.0
}

// array containing the homogeneous weights:
#declare nurbs_weight_array = array[nurbs_vertexes] {
  1.0,
  1.0,
  1.0,
  1.0,
  1.0,
}

// .----------------------------------------------------------------------.
// |  internal NURBS data:                                                |
// '----------------------------------------------------------------------'

// knot vectors:  
#declare nurbs_knot_nr_array = array[nurbs_vertexes+nurbs_order*2-1]
// array containing the rationalbasis functions:
#declare nurbs_nbasis_array = array[20]
#declare nurbs_temp_array = array[35]
// result:
#declare nurbs_p = array[3]

// .----------------------------------------------------------------------.
// |  Macro:                                                              |
// '----------------------------------------------------------------------'

#macro macro_nurbs_initialize()
  #local npluso = nurbs_vertexes + nurbs_order; // maximum value of the knot vector

  #local nurbs_i = 0;
  #while( nurbs_i <= npluso )
    #declare nurbs_knot_nr_array[nurbs_i] = 0.0;    // fixme: not needed ?
    #local nurbs_i = nurbs_i + 1;
  #end
  
  // generate the uniform open knot vector:
  //  generate a B-spline open knot vector with multiplicity equal to the order at the ends.
  #local nplus2 = nurbs_vertexes + 1; // index of x() for the first occurence of the maximum knot vector value
  #declare nurbs_knot_nr_array[0] = 0;
  #local nurbs_i = 1; // (next)
  #while( nurbs_i < npluso )
    #declare nurbs_knot_nr_array[nurbs_i] = nurbs_knot_nr_array[nurbs_i-1];
    #if ( nurbs_i >= nurbs_order )
      #if ( nurbs_i < nplus2 )
        #declare nurbs_knot_nr_array[nurbs_i] = nurbs_knot_nr_array[nurbs_i] + 1;
      #end      
    #end
    #local nurbs_i = nurbs_i + 1;
  #end
#end

// .----------------------------------------------------------------------.
// '----------------------------------------------------------------------'

#macro macro_nurbs_get_max_t()
  #local npluso = nurbs_vertexes + nurbs_order;
  nurbs_knot_nr_array[npluso-1]
#end

// .----------------------------------------------------------------------.
// '----------------------------------------------------------------------'

#macro macro_nurbs_get_point( nurbs_t_requested )

  #local npluso = nurbs_vertexes + nurbs_order; // maximum value of the knot vector
  
  // calculate a point on the rational B-spline curve:
  
  #if (nurbs_knot_nr_array[npluso-1] - nurbs_t_requested <  0.00005) 
    #local nurbs_t_requested = nurbs_knot_nr_array[npluso-1];
  #end
  
  // generate rational B-spline basis functions--open knot vector:
  
  // calculate the first order nonrational basis functions n[i]:
  #local nurbs_i = 0;
  #while( nurbs_i < (npluso-1) )
    #declare nurbs_temp_array[nurbs_i] = 0.0;
    #if ( nurbs_t_requested >= nurbs_knot_nr_array[nurbs_i] )
      #if ( nurbs_t_requested < nurbs_knot_nr_array[nurbs_i+1] )
        #declare nurbs_temp_array[nurbs_i] = 1.0;
      #end
    #end
    #local nurbs_i = nurbs_i + 1;
  #end

  // calculate the higher order nonrational basis functions:
  #local nurbs_k = 2;
  #while( nurbs_k <= nurbs_order )
    #local nurbs_i = 0;
    #while( nurbs_i < (npluso-nurbs_k) )
      // first term of the basis function recursion relation:
      #if ( nurbs_temp_array[nurbs_i] = 0 )
        #local nurbs_d = 0.0;
      #else
        #local nurbs_d = ((nurbs_t_requested-nurbs_knot_nr_array[nurbs_i])*nurbs_temp_array[nurbs_i])/(nurbs_knot_nr_array[nurbs_i+nurbs_k-1]-nurbs_knot_nr_array[nurbs_i]);
      #end
      // second term of the basis function recursion relation:
      #if ( nurbs_temp_array[nurbs_i+1] = 0 )
        #local nurbs_e = 0.0;
      #else
        #local nurbs_e = ((nurbs_knot_nr_array[nurbs_i+nurbs_k]-nurbs_t_requested)*nurbs_temp_array[nurbs_i+1])/(nurbs_knot_nr_array[nurbs_i+nurbs_k]-nurbs_knot_nr_array[nurbs_i+1]);
      #end
      #declare nurbs_temp_array[nurbs_i] = nurbs_d + nurbs_e;
      #local nurbs_i = nurbs_i + 1;
    #end
    #local nurbs_k = nurbs_k + 1;
  #end  
  
  #if (nurbs_t_requested = nurbs_knot_nr_array[npluso-1] ) 
    #declare nurbs_temp_array[nurbs_vertexes-1] = 1.0;    //    pick up last point
  #end

  // calculate sum for denominator of rational basis functions:
  #local nurbs_sum = 0.0;
  #local nurbs_i = 0;
  #while( nurbs_i < nurbs_vertexes )
    #local nurbs_sum = nurbs_sum + nurbs_temp_array[nurbs_i] * nurbs_weight_array[nurbs_i];
    #local nurbs_i = nurbs_i + 1;
  #end
  
  // form rational basis functions and put in r vector:
  #local nurbs_i = 0;
  #while( nurbs_i < nurbs_vertexes )
    #if( nurbs_sum = 0 )
      #declare nurbs_nbasis_array[nurbs_i] = 0.0;
    #else
      #declare nurbs_nbasis_array[nurbs_i] = (nurbs_temp_array[nurbs_i]*nurbs_weight_array[nurbs_i]) / nurbs_sum;
    #end
    #local nurbs_i = nurbs_i + 1;
  #end
   
  #local nurbs_j = 0;
  #while( nurbs_j <= 2 ) //(do for x, y and z)
    #local jcount = nurbs_j;
    #declare nurbs_p[nurbs_j] = 0.0;
    #local nurbs_i = 0;
    #while( nurbs_i < nurbs_vertexes )
      #declare nurbs_p[nurbs_j] = nurbs_p[nurbs_j] + nurbs_nbasis_array[nurbs_i]*nurbs_vert_array[jcount];
      #local jcount = jcount + 3;
      #local nurbs_i = nurbs_i + 1;
    #end
    #local nurbs_j = nurbs_j + 1;
  #end
  // (result is <nurbs_p[0],nurbs_p[1],nurbs_p[2]> )
#end

// .----------------------------------------------------------------------.
// |  Test:                                                               |
// '----------------------------------------------------------------------'

#debug " -=- -=- -=- -=- NURBS curve macro test: -=- -=- -=- -=-\n"
  
macro_nurbs_initialize()
#declare nurbs_max_t = macro_nurbs_get_max_t();

#declare calculate_points = 21;

#local nurbs_t = 0.0;
#while( nurbs_t <= nurbs_max_t )
  macro_nurbs_get_point(nurbs_t)
  #debug concat("result for nurbs_t = ",str(nurbs_t,0,5)," : ")
  #debug concat("<",str(nurbs_p[0],0,5))
  #debug concat(",",str(nurbs_p[1],0,5))
  #debug concat(",",str(nurbs_p[2],0,5),">\n")
  #local nurbs_t = nurbs_t + (nurbs_max_t/(calculate_points-1));
#end

#debug " -=- -=- -=- -=- -=- -=- -=- -=- -=- -=- -=- -=- -=- -=-\n"
  
/*
results of C code: (result of this macro should be the same with same input)

Polygon points:
 0.000000 0.000000 1.000000
 1.000000 2.000000 1.000000
 2.500000 0.000000 1.000000
 4.000000 2.000000 1.000000
 5.000000 0.000000 1.000000

Homogeneous weighting vector:
 1.000000  1.000000  1.000000  1.000000  1.000000

Curve points:
 0.000000 0.000000 1.000000
 0.292750 0.514000 1.000000
 0.572000 0.872000 1.000000
 0.839250 1.098000 1.000000
 1.096000 1.216000 1.000000
 1.343750 1.250000 1.000000
 1.584000 1.224000 1.000000
 1.818250 1.162000 1.000000
 2.048000 1.088000 1.000000
 2.274750 1.026000 1.000000
 2.500000 1.000000 1.000000
 2.725251 1.026000 1.000000
 2.952000 1.088000 1.000000
 3.181751 1.162000 1.000000
 3.416000 1.224000 1.000000
 3.656251 1.250000 1.000000
 3.904001 1.216000 1.000000
 4.160751 1.098000 1.000000
 4.428001 0.871999 1.000000
 4.707251 0.513999 1.000000
 5.000000 0.000000 1.000000
*/