/*
Include file: curve.inc
Spline and curve interpolation for many variants.

Based upon the works of Steve H. Noskowicz !

http://web.archive.org/web/20151002232205/http://home.comcast.net/~k9dci/site/?/page/Piecewise_Polynomial_Interpolation/
Get the book PDF and the appendix.

Author : Ingo Janssen
 
Version: 0.1 Alpha, 2019-06-16
         initial
Version: 0.2 Alpha, 2021-05-24
         complete re-write from earlier attempt.
         cubic Lagrange still not right.
Version: 0.3 Alpha, 2021-06-06
         yet again a complete re-write. Simplification.
         cubic Lagrange still not right, left out.

Run the file without including to render the tests at the end of the file.
These also show how to use the macro's. 

Curve Interpolation macros:
  CIlinear(VectorArray)
  CIsimple3(VectorArray)
  CIsimple3ESC(VectorArray, EndSlope)
  CIsimple3MSC(VectorArray, MidSlope)
  CIsimple3IESC(VectorArray, EndSlope, EndSlope)
  
  CIbezier2(VectorArray)
  CIbezier2EVC(VectorArray, Velocity)
  CIbezier3(VectorArray)
  CIbezier3EVC(VectorArray, Velocity)
  
  CIbspline2(VectorArray)
  CIbspline3(VectorArray)
  
  CIgolden(VectorArray)
  CI4p3(VectorArray)
  CIhermite(VectorArray)
  CIcatmul(VectorArray)
  CIlagrange2(VectorArray, T2)

  CIbeta(VectorArray, Tension, Bias)
  CItcb(VectorArray, Tension, Continuity, Bias)
  CIbezier3EVC(VectorArray, Velocity)
  CIpalmer(VectorArray, Tension, Bias)

#declare BzCurve = CIbezier3(array[n]{....});
#declare PointOnBz = CIget(BzCurve, 0.5);

Curve data structure:

Curve = array[4]{
  _cu_ctrl: array[6]{_cu_dimseg, _cu_step, _cu_end, _cu_segments, _cu_curseg},
  _cu_point: array[n]{<>,<>, ...},
  _cu_curmult: array[len matrix]{<>, <>, ...},
  _cu_def: array[]{
      _cu_var:array[1]{_cu_div},
      _cu_matrix:array[n]{
        array[n]{martrix parameters},
        array[n]{..},
        array[n]{..},
        array ...
      }
  },
 }
*/

#version 3.7;
#ifndef(CURVE_INC_TEMP)
  #declare CURVE_INC_TEMP = version;
  #ifdef(View_POV_Include_Stack)
    #debug "including curve.inc\n"
  #end

//enum CICurve items
#declare _cu_ctrl = 0;
#declare _cu_point = 1;
#declare _cu_curmult = 2; //store for per section precalcultated vector values
#declare _cu_def = 3;
//enum curve ctrl items
#declare _cu_dimseg = 0;
#declare _cu_step = 1;
#declare _cu_end = 2;    //not used in calculations
#declare _cu_len = 3;
#declare _cu_segments = 4;
#declare _cu_curseg = 5;
//enumerate secondary curve type definitions
#declare _cu_var = 0;
#declare _cu_matrix = 1; 
//enumerate variables in _cu_var
#declare _cu_div = 0;


#macro _arrFill(Len, Val)
  /*:Create and fill an array of length Len with value Val*/
  #local Arr = array[Len];
  #for(i,0,Len-1)
    #local Arr[i] = Val;
  #end
  Arr
#end 

#macro _CIinit(CIarr, DimSeg, Step, End, Div, Def)
  #local Len = dimension_size(CIarr, 1);
  #debug concat("len : ",str(Len,0,-1),"\n")
  #local Segments = (Len-1) / Step;
  #debug concat("seg : ",str(Segments,0,-1),"\n")
  #if(Len < DimSeg)
    #error "Too few points in curve array"
  #end
  #if (Segments - int(Segments) != 0)
    #error "Not the right amount of points in curve array"
  #end
  #local CurrSeg = -1;
  #local CICurve = array[4];
  #local CICurve[_cu_ctrl] = array[6]{DimSeg, Step, End, Len, Segments, CurrSeg};
  #local CICurve[_cu_point] = CIarr; 
  #local CICurve[_cu_curmult] = _arrFill(dimension_size(Def[_cu_matrix],1), <0,0,0>);
  #local CICurve[_cu_def] = Def;
  CICurve
#end 

#macro CIlinear(CIarr)
  /*:Linear interpolation*/
  #local DimSeg = 2;
  #local Step = 1;
  #local End = 1;
  #local Div = 1;
  #local Def = array[2]{
    array[1]{Div},              
    array[2]{
      array[2]{-1, 1},
      array[1]{ 1}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIsimple3(CIarr)
  /*:Simple Cubic, linear sections*/
  #local DimSeg = 2;
  #local Step = 1;
  #local End = 1;
  #local Div = 1;
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[2]{ 2,-2},
      array[2]{-3, 3},
      array[2]{ 0, 0},
      array[1]{ 1}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIsimple3ESC(CIarr, ESC)
  /*:Simple Cubic, Endpoint Slope Control*/
  #local DimSeg = 2;
  #local Step = 1;
  #local End = 1;
  #local Div = 1;
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[2]{ 2 - (2 * ESC), (2 * ESC) - 2},
      array[2]{ (3 * ESC) - 3, 3 - (3 * ESC)},
      array[2]{-ESC, ESC},
      array[1]{ 1}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIsimple3MSC(CIarr, MSC)
  /*:Simple Cubic, Midpoint Slope Control*/
  #local DimSeg = 2;
  #local Step = 1;
  #local End = 1;
  #local Div = 1;
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[2]{ (4 * MSC) - 4, 4 - (4 * MSC)},
      array[2]{ 6 - (6 * MSC), (6 * MSC) - 6},
      array[2]{ (2 * MSC) - 3, 3 - (2 * MSC)},
      array[1]{ 1}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIsimple3IESC(CIarr, ESC0, ESC1)
  /*:Simple Cubic, Independent Endpoint Slope Control*/
  #local DimSeg = 2;
  #local Step = 1;
  #local End = 1;
  #local Div = 1;
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[2]{ 2 - ESC0 - ESC1, ESC0 + ESC1 - 2},
      array[2]{ (2 * ESC0) + ESC1 - 3, 3 - (2 * ESC0) - ESC1},
      array[2]{ -ESC0, ESC0},
      array[1]{ 1}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIbezier2(CIarr)
  //:Quadratic Bezier curve
  #local DimSeg = 3;
  #local Step = 2;
  #local End = 2;
  #local Div = 1;
  #local Def = array[2]{
    array[1]{Div},              
    array[3]{
      array[3]{ 1,-2, 1},
      array[2]{-2, 2},
      array[1]{ 1},
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIbezier2EVC(CIarr, EVC)
  //:Quadratic Bezier curve with End Velocity Control
  #local DimSeg = 3;
  #local Step = 2;
  #local End = 2;
  #local Div = 1;
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[3]{-EVC, 0.0, EVC},
      array[3]{ 1 + (2 * EVC), -2 - EVC, 1 - EVC},
      array[2]{-2 - EVC, 2 + EVC},
      array[1]{ 1},
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIbezier3(CIarr)
  //:Cubic Bezier curve
  #local DimSeg = 4;
  #local Step = 3;
  #local End = 3;
  #local Div = 1;
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[4]{-1, 3,-3, 1},
      array[3]{ 3,-6, 3},
      array[2]{-3, 3},
      array[1]{ 1}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIbezier3EVC(CIarr, EVC)
  //:Cubic Bezier curve with End Velocity Control
  #local DimSeg = 4;
  #local Step = 3;
  #local End = 3;
  #local Div = 1;
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[4]{-1.0 - EVC, 3.0 + EVC, -3.0 - EVC, 1.0 + EVC},
      array[3]{ 3.0 + (2 * EVC), -6.0 - (2 * EVC), 3.0 + EVC},
      array[2]{-3.0 - EVC, 3.0 + EVC},
      array[1]{ 1}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIbspline2(CIarr)
  //:Cubic b-spline
  #local DimSeg = 3;
  #local Step = 1;
  #local End = 1;
  #local Div = 2;
  #local Def = array[2]{
    array[1]{Div},              
    array[3]{
      array[3]{ 1,-2, 1},
      array[2]{-2, 2},
      array[2]{ 1, 1}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIbspline3(CIarr)
  //:Cubic b-spline
  #local DimSeg = 4;
  #local Step = 1;
  #local End = 2;
  #local Div = 6;
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[4]{-1, 3,-3, 1},
      array[3]{ 3,-6, 3},
      array[3]{-3, 0, 3},
      array[3]{ 1, 4, 1}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIgolden(CIarr)
  //:Golden curve
  #local DimSeg = 3;
  #local Step = 2;
  #local End = 2;
  #local Div = 1;
  #local Def = array[2]{
    array[1]{Div},              
    array[3]{
      array[3]{ 2,-4, 2},
      array[3]{-3, 4,-1},
      array[1]{1},
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CI4p3(CIarr)
  //:Four Point Cubic
  #local DimSeg = 4;
  #local Step = 4;
  #local End = 4;
  #local Div = 2;
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[4]{ -9, 27,-27, 9},
      array[4]{ 18,-45, 36,-9},
      array[4]{-11, 18, -9, 2},
      array[1]{  2}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIhermite(CIarr)
  //:Cubic Hermite spline
  #local DimSeg = 4;
  #local Step = 3;
  #local End = 3;
  #local Div = 1;
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[4]{ 2, 1, 1,-2},
      array[4]{-3,-2,-1, 3},
      array[2]{ 0, 1},
      array[1]{ 1}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIcatmull(CIarr)
  //:Catmull Rom spline
  #local DimSeg = 4;
  #local Step = 1;
  #local End = 2;
  #local Div = 2;
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[4]{-1, 3,-3, 1},
      array[4]{ 2,-5, 4,-1},
      array[3]{-1, 0, 1},
      array[2]{0,2}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIlagrange2(CIarr, T2)
  //:Quadratic Lagrange curve
  #local DimSeg = 3;
  #local Step = 2;
  #local End = 2;
  #local Div = 1;
  #local N1 = 1/T2;
  #local N2 = 1/(T2 * (T2 - 1));
  #local N3 = 1/(1 - T2);     
  #local Def = array[2]{
    array[1]{Div},              
    array[3]{
      array[3]{ N1, N2, N3},
      array[3]{-N1 - 1, -N2, -N3 * T2},
      array[1]{ 1},
     }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIbeta(CIarr, T, B)
  //:Beta-spline
  #local DimSeg = 4;
  #local Step = 1;
  #local End = 2;
  #local Div = (T + (2 * pow(B, 3)) + (4 * pow(B,2)) + (4 * B) + 2);
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[4]{
        -2 * pow(B, 3),
         2 * (T + pow(B, 3) + pow(B, 2) + B), 
        -2 * (T + pow(B, 2) + B + 1),         
         2        
      },
      array[3]{
         6 * pow(B, 3),                             
        -3 * (T + (2 * pow(B, 3) + (2 * pow(B, 2)))),
         3 * (T + (2 * pow(B, 2)))
      },
      array[3]{
        -6 * pow(B, 3),     
         6 * (pow(B, 3) - B),
         6 * B  
      },
      array[3]{
         2 * pow(B, 3),            
         T + (4 * (pow(B, 2) + B)),
         2
      }
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CItcb(CIarr, T, C, B)
  //:Kochanek Bartels aka TCB spline, Tension, Continuity, Bias
  #local DimSeg = 4;
  #local Step = 1;
  #local End = 2;
  #local Div = 2;
  #local A = (1 - T) * (1 + C) * (1 + B);
  #local B = (1 - T) * (1 - C) * (1 - B);
  #local C = (1 - T) * (1 - C) * (1 + B);
  #local D = (1 - T) * (1 + C) * (1 - B);
  #local Def = array[2]{
    array[1]{Div},              
    array[4]{
      array[4]{
        -A,          
         4 + A - B - C,
        -4 + A + C - D,
         D       
      },
      array[4]{
         2 * A,
        -6 - (2 * A) + (2 * B) + C,
         6 - (2 * B) - C + D,
        -D
      },
      array[3]{
        -A,
         A - B,
         B      
      },
      array[2]{0,2}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end

#macro CIpalmer(CIarr, T, B)
  //:Palmer
  #local DimSeg = 3;
  #local Step = 2;
  #local End = 2;
  #local Div = 1;
  #local N0 = 1 / B;
  #local N1 = 1 / (B * (1 - B));
  #local N2 = 1 / (1 - B);
  #local Def = array[2]{
    array[1]{Div},              
      array[4]{
        array[3]{
          -2 * T, 
          0, 
          2 * T,
        },
        array[3]{
           N0 + (3 * T),
          -N1,
           N2 - (3 * T),
        },
        array[3]{
          -1 - N0 - T,
           N1,
           T - (N2 / N0),
        },
        array[1]{
          1
        }
      }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end
/*template
#macro CI*******(CIarr)
  //:****************************
  #local DimSeg = *;
  #local Step = *;
  #local End = *;
  #local Div = *;
  #local Def = array[2]{
    array[1]{Div},              
    array[*]{
      array[*]{*,*,*,*,},
      array[*]{*,*,*,*},
      array[*]{*,*,*,*},
      array[*]{*,*,*,*}
    }
  }
  _CIinit(CIarr, DimSeg, Step, End, Div, Def)
#end
*/


#macro CIget(CICurve, T)
  /*: Gets the point on the curve at T
  
  parameter:
    CICurve (array, vector): an array of vectors and other parameters as created via one 
                              of the spline initialisation macros. 
    T (float, [0,1]): position at which to interpolate the curve.
  
  return: 
    (vector) returns a point on the curve at T.
  
  use example:
    #declare TheCurve = CIlinear(PointArray);
    #declare P = CIget(TheCurve, 0.1);    
  */
  
  /* T goes from 0 to 1 for the whole curve. From current T find the proper curve segment.
     Tseg also goes from 0 to 1. Once the proper segment is found, get the proper Tseg.
  */
  #local Len = CICurve[_cu_ctrl][_cu_len];
  #local Segments = CICurve[_cu_ctrl][_cu_segments];
  #local CImult = CICurve[_cu_curmult];
  #if(T = 1)
    #local Segment = Segments-1;
    #local Tseg = 1;
  #else
    #local Segment = floor(T * Segments);
    #local Tseg =  (T * Segments) - Segment;    
  #end
  #if(Segment != CICurve[_cu_ctrl][_cu_curseg])
    #local CICurve[_cu_ctrl][_cu_curseg] = Segment;  
    #local CurrPos = Segment * CICurve[_cu_ctrl][_cu_step];
    #local CImatrix = CICurve[_cu_def][_cu_matrix];
    #local Div = CICurve[_cu_def][_cu_var][_cu_div];  
    #local MaxPos = Len - CICurve[_cu_ctrl][_cu_dimseg] + CICurve[_cu_ctrl][_cu_step];
    #if(CurrPos < MaxPos)
      #for(i, 0, dimension_size(CImatrix, 1) - 1)
        #local CImult[i] = <0, 0, 0>;
        #for(j, 0, dimension_size(CImatrix[i], 1) - 1)
          #local jSegment = j + CurrPos;
          #local CImult[i] = CImult[i] + (CImatrix[i][j] / Div) * CICurve[_cu_point][jSegment];    
        #end
      #end
    #end
    #local CICurve[_cu_curmult] = CImult; 
  #end
  #local F = CImult[0];
  #for(i, 1, dimension_size(CImult, 1) - 1)
    #local F = (F * Tseg) + CImult[i];
  #end
  F
#end



#if(input_file_name="curve.inc")
  //tests
  global_settings{ assumed_gamma 1.0 }
  #default{ finish{ ambient 0.1 diffuse 0.9 emission 0}} 

  #declare Lin = CIlinear(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Simp3 = CIsimple3(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Bspl2 = CIbspline2(
    array[4]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Bez2 = CIbezier2(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Gold = CIgolden(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Bez3 = CIbezier3(
    array[4]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Four3 = CI4p3(
    array[5]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>,<2,0,0>}
  );
  #declare Herm = CIhermite(
    array[4]{<0, 0, 0>, <0, 2, 0>, <1, 2, 0>, <1, 0, 0>}
  );
  #declare Bspl3 = CIbspline3(
    array[4]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>}
  );
  #declare Catmull = CIcatmull(
    array[6]{<-1, 0, 0>,<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>,<2,0,0>}
  );
  #declare ESC = CIsimple3ESC(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>},
    3 //End point Slope control  parameter
  );
  #declare MSC = CIsimple3MSC(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>},
    0.5 //Mid point Slope control  parameter
  );          
  #declare IESC = CIsimple3IESC(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>},
    -.5, 2  //Independent End point Slope control parameters
  );
  #declare Bez2EVC = CIbezier2EVC(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>},
    -3 //End point Velocity Control parameter
  );
  #declare Lag2 = CIlagrange2(
    array[5]{<0, 0, 0>, <0, 1, 0>, <0.5, 0.5, 0>, <1, 1, 0>, <1, 0, 0>},
    0.3 
  ); 
  #declare Beta = CIbeta(
    array[6]{<0,0,0>, <0,1,0>, <1,1,0>, <1,0,0>, <2,0,0>, <2,1,0>}
    0, 1  //Tension, Bias
  )
  #declare TCB = CItcb(
    array[6]{<0,0,0>, <0,1,0>, <1,1,0>, <1,0,0>, <2,0,0>, <2,1,0>}
    0, 0, 0 //Tension, Continuity, Bias
  )
  #declare Bez3EVC = CIbezier3EVC(
    array[7]{<0, 0, 0>, <0, 1, 0>, <1, 1, 0>, <1, 0, 0>, <1,-1,0>, <2,-1,0>, <2,0,0>},
    2 //End point Velocity Control parameter
  );
   #declare Palmer = CIpalmer(
    array[5]{<0, 0, 0>, <0, 1, 0>, <1, 0.5, 0>, <2,1,0>,<2,0,0>},
    0, 0.5 //Tension, Bias 
  );


  
  // Set curve for test render //
  #declare Test =  Palmer;
  
  //get the points on the curve
  #for(T, 0, 1, 0.01)
    sphere{CIget(Test, T), 0.02 pigment{rgb 1}}
  #end
  
  //control points
  sphere{Test[1][0] 0.03 pigment{rgb <1,0,0>}}
  #for(i, 1, dimension_size(Test[1],1) - 1)
    sphere{Test[1][i] 0.03 pigment{rgb <1,0,0>}}
    cylinder{Test[1][i-1], Test[1][i], 0.01 pigment{rgb <1,0,0>}}
  #end
  
  camera {
    perspective angle 95 
    location  <0.5, 0.5,-2>
    right     x*image_width/image_height
    look_at   <0.5, 0.5, 0.0>
  }
  
  light_source{<0, 0,-3000> color rgb 1}  

#end

#version CURVE_INC_TEMP;
#end