#version 3.7;
global_settings { assumed_gamma 1.0 }

//------------------------------------------
// SDL for making a Bezier spline - given 4 data points, 
// interpolate control points p1 and p2 for 3 segments - 6 control points total
// Bill Walker - 2016
// Adapted code to POV-Ray from :
// https://www.particleincell.com/2012/bezier-splines/
// https://www.particleincell.com/wp-content/uploads/2012/06/circles.svg
// replace var with #declare or #local, function with #macro, edit #for() - #end syntax, array[] syntax
//------------------------------------------

#include "colors.inc"
#include "debug.inc"
	Set_Debug (true)
#include "math.inc"

light_source { <0, 50, -50>  color rgb <1, 1, 1>}
camera {location  <350, 150, -800> look_at   <350, 150, 0>}



//------------------------------
// begin macro code

	
//##############################################################################################################
#macro InterpolateControlPointsFrom (InputPoints)
//  macro receives an array of 4 data points (vectors - <x,y,z>) with 2 empty spaces between each
//  macro returns an array with 6 new control points filled in to make a continuous spline
                            

//#declare S = array[4];

#local K = array[4];				// working array for control point x, y, and z components
#local P = array[6];				// output array of component interpolation macro
#local CPX = array[6];				// intermediate collection array for control point x components
#local CPY = array[6];				// intermediate collection array for control point y components
#local CPX = array[6];				// intermediate collection array for control point z components
#local Computed_spline = array [10];

// Split vector data point array into separate x, y, and z component arrays
// data points are at array elements |0--3--6--9|
#local X = array[4] {InputPoints[0].x, InputPoints[3].x, InputPoints[6].x, InputPoints[9].x};
#local Y = array[4] {InputPoints[0].y, InputPoints[3].y, InputPoints[6].y, InputPoints[9].y};
#local Z = array[4] {InputPoints[0].z, InputPoints[3].z, InputPoints[6].z, InputPoints[9].z};





	#macro computeControlPoints(K)
	//computes control points given knots K, this is the brain of the operation
		#local p1 = array[4];
		#local p2 = array[4];
		#local n = dimension_size (K, 1) - 1;
			
			//rhs vector
			#local a = array[n];
			#local b = array[n];
			#local c = array[n];
			#local r = array[n];
			
			//left most segment
			#local a[0] = 0;
			#local b[0] = 2;
			#local c[0] = 1;
			#local r[0] = K[0] + 2 * K[1];
			
			//middle segment
			#for (i, 0, n-2)
				#local a[i] = 1;
				#local b[i] = 4;
				#local c[i] = 1;
				#local r[i] = 4 * K[i] + 2 * K[i+1];
			#end
					
			//right segment
			#local a[n-1] = 2;					
			#local b[n-1] = 7;
			#local c[n-1] = 0;
			#local r[n-1] = 8 * K[n-1] + K[n];
			
			//------------------------------------------------------------
			
			//solves Ax=b with the Thomas algorithm (from Wikipedia)
			#for (i, 1, n-1)
				#local m = a[i]/b[i-1];
				#local b[i] = b[i] - m * c[i - 1];
				#local r[i] = r[i] - m*r[i-1];
			#end
		 
			#local p1[n-1] = r[n-1]/b[n-1];
			#debug concat( "Defining p1.   n-1 = ", str((n-1), 3, 1), "\n")
	
	
			
			#for (i, n-2, 0, -1) 				//#for (i = n - 2; i >= 0; --i)
				#local p1[i] = (r[i] - c[i] * p1[i+1]) / b[i];
				#debug concat( "Defining p1.   i = ", str(i, 3, 1), "  Defining p1.   i+1 = ", str((i+1), 3, 1),"\n")
			#end
	
	
	
			//we have p1, now compute p2
			#for (i, 0, n-2) 				//#for (i=0;i<n-1;i++)
				#debug concat( "Defining p2   i = ", str(i, 3, 1), "  i+1 = ", str((i+1), 3, 1),"\n")
				#local p2[i] = 2 * K[i+1] - p1[i+1];
			#end
			
			#local p2[n-1] = 0.5 * (K[n] + p1[n-1]);
	
	
			// Interpolated control points for left segment | Interpolated control points for middle segment | Interpolated control points for right segment
			#debug concat( "p1[0] = ", str(p1[0], 3, 1), 	"  p1[1] = ", str(p1[1], 3, 1), 				"  p1[2] = ", str(p1[2], 3, 1), "\n")
			#debug concat( "p2[0] = ", str(p2[0], 3, 1), 	"  p2[1] = ", str(p2[1], 3, 1), 				"  p2[2] = ", str(p2[2], 3, 1), "\n")
			//  first calculated for X components, then for Y components, then for Z.
	
			//return {p1:p1, p2:p2};
			#local P[0] = p1[0]; #local P[1] = p2[0]; #local P[2] = p1[1]; #local P[3] = p2[1]; #local P[4] = p1[2]; #local P[5] = p2[2];
			
			(P)	// return components stored in the array P[] to outer macro
	#end
	
// call macro to Interpolate control point components for Bezier spline segments
// Sets of p1 and p2 for Right, Middle, and left segments - 6 components in each array
#if (Verbose = true)
	#debug "Interpolating Control Point x components \n"
#end
#local CPX = computeControlPoints(X);
#for (i, 0, 5) // for i = 0 to 5
	#debug concat( "i =", str(i, 3, 1))
	#debug concat( "CPX [", str(i, 3, 1), "] = ", str(CPX[i], 3, 1), "     ")
#end
#debug " \n"


#if (Verbose = true)
	#debug "Interpolating Control Point y components \n")
#end
#local CPY = computeControlPoints(Y);
#for (i, 0, 5) // for i = 0 to 5
	#debug concat( "CPY [", str(i, 3, 1), "] = ", str(CPY[i], 3, 1), "     ")
#end
#debug " \n"


#if (Verbose = true)
	#debug "Interpolating Control Point z components \n")
#end
#local CPZ = computeControlPoints(Z);
#for (i, 0, 5) // for i = 0 to 5
	#debug concat( "CPZ [", str(i, 3, 1), "] = ", str(CPZ[i], 3, 1), "     ")
#end
#debug " \n"


#declare Computed_spline = {	InputPoints[0], <CPX[0], CPY[0], CPZ[0]>, <CPX[1], CPY[1], CPZ[1]>, 
						InputPoints[3], <CPX[2], CPY[2], CPZ[2]>, <CPX[3], CPY[3], CPZ[3]>, 
						InputPoints[6], <CPX[4], CPY[4], CPZ[4]>, <CPX[5], CPY[5], CPZ[5]>, InputPoints[9]};
						
Computed_spline	// return spline from macro to SDL

#end // end main macro		

//#################################################################################################################################################################



// continue scene SDL

#declare Given_data = array[10] {<60, 60, 0>, <0,0,0>, <0,0,0>, <220, 300, 0>, <0,0,0>, <0,0,0>, <420, 300, 0>, <0,0,0>, <0,0,0>, <700, 240, 0>};

#declare Verbose = true;
// call interpolation macro
#declare Interpolated_data = InterpolateControlPointsFrom (Given_data)

#declare Sphere_rad = 10;

// render data points and control points
#declare Markers = union{
	sphere { Interpolated_data[0] Sphere_rad pigment {Yellow} }
	sphere { Interpolated_data[1] Sphere_rad pigment {Red} }
	sphere { Interpolated_data[2] Sphere_rad pigment {Red} }
	sphere { Interpolated_data[3] Sphere_rad pigment {Yellow} }
	sphere { Interpolated_data[4] Sphere_rad pigment {Green} }
	sphere { Interpolated_data[5] Sphere_rad pigment {Green} }
	sphere { Interpolated_data[6] Sphere_rad pigment {Yellow} }
	sphere { Interpolated_data[7] Sphere_rad pigment {Blue} }
	sphere { Interpolated_data[8] Sphere_rad pigment {Blue} }
	sphere { Interpolated_data[9] Sphere_rad pigment {Yellow} }
}

#declare Bezier_spline = union {
#for (i, 0, 6, 3)
	bicubic_patch {type 1 flatness 0 u_steps 4 v_steps 4
		Interpolated_data[i+0] Interpolated_data[i+1] Interpolated_data[i+2] Interpolated_data[1+3]
		Interpolated_data[i+0] Interpolated_data[i+1] Interpolated_data[i+2] Interpolated_data[1+3]
		Interpolated_data[i+0] Interpolated_data[i+1] Interpolated_data[i+2] Interpolated_data[1+3]
		Interpolated_data[i+0] Interpolated_data[i+1] Interpolated_data[i+2] Interpolated_data[1+3]
		pigment{Red} scale z*Sphere_rad scale y*-1
	} // end bicubic-patch
#end // end for-loop
} // end union

object {Markers}
object {Bezier_spline}

// end interpolation testing SDL