// Persistence of Vision Ray Tracer version 3.5 Include File // File: Cubic_Splines // By: Christopher D. Johnson // Last updated: 06.11.2002 #ifndef(Cubic_Splines_Temp) // if include file has already been included, do nothing #declare Cubic_Splines_Temp = version; #ifdef(View_POV_Include_Stack) #debug "including cubic_splines.inc\n" #end // ****** INCLUDES ****** #include "rand.inc" #include "math.inc" // ****** DECLARES ****** #declare Data = array [3] {0,0,0}; // [0] = time, [1] = segement length, [2] = k_error #declare CGOLD = ((3-sqrt(5))/2); #declare GOLD = 2 - CGOLD; #declare FPP = 1 / (10000); // ****** UTILITY MACROS ****** // *** Set : just makes an array out of the three arguments #macro Set(X,Y,Z) #local temp = array [3] {X,Y,Z}; temp #end // *** SHIFT : a simple housekeeping macro #macro SHIFT (A,B,C,D) #declare A = B; #declare B = C; #declare C = D; #end // *** SIGN : macro used by MIN/MAX routines #macro SIGN (A,B) #if (B >= 0) #local F = abs(A); #else #local F = -abs(A); #end F #end // *** subdivide : used by seg_length to find distances #macro subdivide (SPLINE,BEGIN,END,BEGIN_POINT,END_POINT,SEED,ERROR) #local MID_TIME = RRand (BEGIN,END,SEED); // m #local MID_POINT = ( SPLINE (MID_TIME)); // vm // determine if more subdivision is required #local ORIG_LEN = (vlength ((END_POINT - BEGIN_POINT))); #local MID_LEN = (vlength ((MID_POINT - BEGIN_POINT)) + vlength ((END_POINT - MID_POINT))); #if ((max(ORIG_LEN,MID_LEN) - min(ORIG_LEN,MID_LEN)) <= ERROR) // no more subdivision (ORIG_LEN) #else #local LENGTH = (subdivide (SPLINE,BEGIN,MID_TIME,BEGIN_POINT,MID_POINT,SEED,ERROR) + subdivide (SPLINE,MID_TIME,END,MID_POINT,END_POINT,SEED,ERROR)); (LENGTH) #end #end // *** segment_length : finds the length along spline between // begin and end to a precision of error. // uses randomized subdivision for midpoints #macro segment_length (SPLINE,BEGIN,END,ERROR) #local len_seed = seed(1); #local BEGIN_POINT = (SPLINE (BEGIN)); // va #local END_POINT = (SPLINE (END)); // vb #local LENGTH = (subdivide (SPLINE,BEGIN,END,BEGIN_POINT,END_POINT,len_seed,ERROR)); (LENGTH) #end // *** Inter_Length : Uses a reference spline to find the length to time. // This uses the precomputed refernce spline instead // of actually calculating th length. It's only as accurate // as the reference spline. #macro Inter_Length (SPLINE,A,ERROR) #local Low = 0; #local High = 1; #local Done = false; #while (!Done) #local Mid = Low + .5 * (High - Low); #local fMid = (SPLINE(Mid)).x - A; #if ((fMid != 0) | ((High - Low) > ERROR)) #if (fMid < 0) #local Low = Mid; #else #local High = Mid; #end #else #local Done = true; #end #debug concat ("! :",str(Mid,0,10),"\n") #end Mid #end // *** Create_Spline : Given an array of points and an array of times, // this creates a cubic spline. #macro Create_Spline (POINTS,TIMES) #local number_of_points = dimension_size (POINTS,1); #local this_point = 0; // Generate spline spline { cubic_spline #while (this_point < number_of_points) TIMES [ this_point ] , POINTS [ this_point ] #local this_point = this_point + 1; #end } #end // *** Estimate_Cubic_Spline : Given an array of points, it creates a cubic spline // by createing times proportional to linear distances. // It also adds the first and last points. I'll fix this later // since it's not really used by anything anymore. #macro Estimate_Cubic_Spline (POINTS) #debug "estim in\n" #local number_of_points = dimension_size (POINTS,1); #local DISTANCES = array[ (number_of_points - 1) ]; #local TIMES = array[ number_of_points ]; #local linear_length = 0; #local this_point = 1; #while (this_point < number_of_points) // Calculate linear distances #local this_length = vlength ((POINTS[this_point] - POINTS[ (this_point - 1) ])); #local linear_length = linear_length + this_length; #local DISTANCES[ (this_point - 1) ] = this_length; #local this_point = this_point + 1; #end // Calculate estimated time values #local TIMES [ 0 ] = 0; #local TIMES [ (number_of_points - 1) ] = 1; #local time_sum = 0; #if (number_of_points > 4) #local this_point = 1; #while (this_point < (number_of_points - 1)) #local time_sum = time_sum + ( DISTANCES [ (this_point - 1) ] / linear_length); #local TIMES [ this_point ] = time_sum; #local this_point = this_point + 1; #end #end // Add extra pre & post control points and copy into WORKING arrays #local WORKING_POINTS = array[ number_of_points + 2 ]; #local WORKING_TIMES = array[ number_of_points + 2 ]; #local WORKING_POINTS [0] = POINTS[0] - ( POINTS[1] - POINTS[0] ); #local WORKING_TIMES [0] = 0 - ( TIMES[1] - TIMES[0] ); #local WORKING_POINTS [ (number_of_points + 1) ] = POINTS[ (number_of_points - 1) ] + ( POINTS [ (number_of_points - 1) ] - POINTS [ (number_of_points - 2) ]); #local WORKING_TIMES [ (number_of_points + 1) ] = 1 + ( TIMES [ (number_of_points - 1) ] - TIMES [ (number_of_points - 2) ]); #local this_point = 1; #while (this_point <= number_of_points) #local WORKING_POINTS [ this_point ] = POINTS [ (this_point - 1) ]; #local WORKING_TIMES [ this_point ] = TIMES [ (this_point - 1) ]; #local this_point = this_point + 1; #end // Create working spline Create_Spline (WORKING_POINTS,WORKING_TIMES) #declare WORK = WORKING_TIMES; #debug "estim out\n" #end // *** FUNC_A : used to find the error between the reference spline and actual distances // A = point to be tested (instance of "Data") // SPLINE = actual spline , K_SPLINE = reference spline // LENGTH = length of SPLINE #macro FUNC_A (A,SPLINE,K_SPLINE,LENGTH) #local F = A; #if (F[0] < 1) #local F[1] = segment_length (SPLINE,0,F[0],FPP); #local F[2] = ((F[0]) - (K_SPLINE((F[1]/LENGTH))).x); #else #local F[1] = LENGTH; #local F[2] = 0; #end array [3] {F[0],F[1],F[2]} #end // *** FUNC_AB : used to find the error between the reference spline and actual distances // A = last point tested , B = point to be tested (both are instances of "Data") // SPLINE = actual spline , K_SPLINE = reference spline // LENGTH = length of SPLINE // this is just a little faster than FUNC_A but does the same thing #macro FUNC_AB (A,B,SPLINE,K_SPLINE,LENGTH) #local F = B; #if (F[0] < 1) #local E = A; #if (E[0] < 0) #local E = Set(0,0,0); #end #if (E[0] > 1) #local E = Set(1,LENGTH,0); #end #local temp = segment_length (SPLINE,E[0],F[0],FPP); //#local temp = Fast_Segment(F[0]) - Fast_Segment(E[0]); #local F[1] = (E[1]) + SIGN(temp,(F[0] - E[0])); #local F[2] = ((F[0])-(K_SPLINE((F[1]/LENGTH))).x); #else #local F[1] = LENGTH; #local F[2] = 0; #end array [3] {F[0],F[1],F[2]} #end // *** Lagrange2 : Lagrange interpolation with two points // A & B are existing points (vectors x,y) // X is value to be solved for // returned L = f(X) #macro Lagrange2 (A,B,X) #local Ax = (((X-B.x)/(A.x-B.x))*A.y); #local Bx = (((X-A.x)/(B.x-A.x))*B.y); #local L = Ax+Bx; L #end // *** Lagrange3 : Lagrange interpolation with three points // A & B & C are existing points (vectors x,y) // X is value to be solved for // returned L = f(X) #macro Lagrange3 (A,B,C,X) #local Ax = ((X-B.x) * (X-C.x)) / ((A.x - B.x) * (A.x - C.x)) * A.y; #local Bx = ((X-A.x) * (X-C.x)) / ((B.x - A.x) * (B.x - C.x)) * B.y; #local Cx = ((X-A.x) * (X-B.x)) / ((C.x - A.x) * (C.x - B.x)) * C.y; #local L = Ax+Bx+Cx; L #end // *** BRACKET : part of the MIN/MAX routines (roughly finds the min/max) // I'll re-write these later on as more generic funtions // right now they are tailored for this function (sort of anyway) // I'm still not happy with something here, just don't know what yet #macro BRACKET (AX,BX,SPLINE,LENGTH,K_SPLINE) // original from numerical methods in C // altered CDJ <04.28.2002> // x = time , y = segment length , z = k_error #local GLIMIT = 100; #local TINY = (1/(pow(10,7))); #local DONE = false; #declare CX = Data; #if (AX[2] > BX[2]) #local MIN = true; // finding minimum #else #local MIN = false; // finding maximum #declare AX[2] = -AX[2]; #declare BX[2] = -BX[2]; #end #declare CX[0] = (BX[0] + (GOLD * (BX[0] - AX[0]))); #declare CX = FUNC_AB(BX,CX,SPLINE,K_SPLINE,LENGTH); #declare CX[2] = (MIN ? CX[2] : -CX[2]); #local U = Data; #while ((BX[2] > CX[2]) & (!DONE)) #local mr = (BX[0] - AX[0]) * (BX[2] - CX[2]); #local mq = (BX[0] - CX[0]) * (BX[2] - AX[2]); #local Num = ((BX[0] - AX[0]) * mr - (BX[0] - CX[0]) * mq); #local Den = 2 * (SIGN ((max (( abs ((mr-mq))),TINY)),(mr-mq))); #local U[0] = BX[0] - (Num / Den); #local ulim = CX[0] + (GOLD * (CX[0]-BX[0])); #if (((BX[0]-U[0])*(U[0]-CX[0])) > 0) #local U = FUNC_AB(BX,U,SPLINE,K_SPLINE,LENGTH); #local U[2] = (MIN ? U[2] : -U[2]); #if (U[2] < CX[2]) #declare AX = BX; #declare BX = U; #local DONE = true; #else #if (U[2] > BX[2]) #declare CX = U; #local DONE = true; #else #local U[0] = CX[0] + GOLD * (CX[0] - BX[0]); #local U = FUNC_AB(BX,U,SPLINE,K_SPLINE,LENGTH); #local U[2] = (MIN ? U[2] : -U[2]); #end #end #else #if (((CX[0] - U[0]) * (U[0] - ulim)) > 0) #local U = FUNC_AB(BX,U,SPLINE,K_SPLINE,LENGTH); #local U[2] = (MIN ? U[2] : -U[2]); #if (U[2] < CX[2]) #declare AX = BX #declare BX = CX; #local U[0] = CX[0] + GOLD * (CX[0] - BX[0]); #local U = FUNC_AB(BX,U,SPLINE,K_SPLINE,LENGTH); #local U[2] = (MIN ? U[2] : -U[2]); #end #else #if (((CX[0] - ulim) * (ulim - U[0])) >= 0) #local U[0] = ulim; #local U = FUNC_AB(BX,U,SPLINE,K_SPLINE,LENGTH); #local U[2] = (MIN ? U[2] : -U[2]); #else #local U[0] = CX[0] + GOLD * (CX[0] - AX[0]); #local U = FUNC_AB(BX,U,SPLINE,K_SPLINE,LENGTH); #local U[2] = (MIN ? U[2] : -U[2]); #end #end #end #if (!DONE) #if ((U[2] = 0) & (CX[2] = 0)) #local MIN = true; #local DONE = true; #declare BX[0] = -1; #else #declare AX = BX; #declare BX = CX; #declare CX = U; #end #end #end //while #declare AX[2] = (MIN ? AX[2] : -AX[2]); #declare BX[2] = (MIN ? BX[2] : -BX[2]); #declare CX[2] = (MIN ? CX[2] : -CX[2]); #end // *** BRENT : the other MIN/MAX routine (pinpoints the MIN/MAX) // I'll also re-write a generic version of this later as these // two functions are fantastic at finding MIN's and MAX's // for a function. // I just wish I had a faster way to compute accurate distances since // that's what slows all this down #macro BRENT (AX,BX,CX,SPLINE,K_SPLINE,LENGTH,TOL) // original from numerical methods in C // altered CDJ <04.28.2002> // x = time , y = segment length , z = k_error #local ITMAX = 100; #local ZEPS = (1/pow(10,10)); #local E = 0; #local D = 0; #local U = Data; #declare A = AX; #declare B = CX; #if (A[2] > BX[2]) #local MIN = true; // finding minimum #else #local MIN = false; // finding maximum #declare A[2] = -A[2]; #declare B[2] = -B[2]; #declare BX[2] = -BX[2]; #end #local V = BX; #local W = BX; #local X = BX; #local ITER = 0; #local DONE = false; #while ((ITER <= ITMAX) & (!DONE)) #local XM = 0.5*(A[0] + B[0]); #local TOL1 = TOL * abs(X[0]) + ZEPS; #local TOL2 = 2.0 * TOL1; #if ((abs((X[0] - XM))) <= (TOL2 - 0.5 * (B[0] - A[0]))) #local DONE = true; #declare BRNT = X; #declare BRNT[2] = (MIN ? BRNT[2] : -BRNT[2]); // would rather just have it return the array X #declare AX = A; #declare CX = B; #else #if ((abs(E)) > TOL1) #local R = (X[0] - W[0]) * (X[2] - V[2]); #local Q = (X[0] - V[0]) * (X[2] - W[2]); #local P = (X[0] - V[0]) * Q - (X[0] - W[0]) * R; #local Q = 2.0 * (Q - R); #if (Q > 0.0) #local P = -P; #end #local Q = abs(Q); #local Etemp = E; #local E = D; #if ((abs(P) >= abs(0.5 * Q * Etemp)) | (P <= (Q * (A[0] - X[0]))) | (P >= (Q * (B[0] - X[0])))) #if (X[0] >= XM) #local E = A[0] - X[0]; #else #local E = B[0] - X[0]; #end #local D = CGOLD * E; #else #local D = P / Q; #local U[0] = X[0] + D; #if (((U[0] - A[0]) < TOL2) | ((B[0] - U[0]) < TOL2)) #local D = SIGN(TOL1,(XM - X[0])); #end #end #else #if (X[0] >= XM) #local E = A[0] - X[0]; #else #local E = B[0] - X[0]; #end #local D = CGOLD * E; #end #if ((abs(D)) >= TOL1) #local U[0] = X[0] + D; #else #local U[0] = X[0] + SIGN(TOL1,D); #end #local U = FUNC_AB (AX,U,SPLINE,K_SPLINE,LENGTH); #local U[2] = (MIN ? U[2] : -U[2]); #if (U[2] <= X[2]) #if (U[0] >= X[0]) #declare A[0] = X[0]; #else #declare B[0] = X[0]; #end SHIFT(V,W,X,U) #else #if (U[0] < X[0]) #declare A[0] = U[0]; #else #declare B[0] = U[0]; #end #if ((U[2] <= W[2]) | (W[0] = X[0])) #local V = W; #local W = U; #else #if ((U[2] <= V[2]) | (V[0] = X[0]) | (V[0] = W[0])) #local V = U; #end #end #end #end #local ITER = ITER + 1; #end // while #if (!DONE) #debug concat ("Too many iterations [",str(ITER,0,0),"], in brent/n") #declare BRNT = X; #declare BRNT[2] = (MIN ? BRNT[2] : -BRNT[2]); #declare AX = A; #declare CX = B; #end #declare AX[2] = (MIN ? AX[2] : -AX[2]); #declare BX = BRNT; #declare CX[2] = (MIN ? CX[2] : -CX[2]); #end // BRENT // *** KINETIC_SPLINE : this is where it all comes together // SPLINE is the original spline // TOL is how accurate you want the reference to be // BAILOUT is the maximum allowable passes to be made // this gets slow with higher TOL & BAILOUT #macro KINETIC_SPLINE (SPLINE,TOL,BAILOUT) #local AX = Data; #local BX = Data; #local CX = Data; #local HIGH = 0; #local BAILOUT = ((BAILOUT >= 1)?BAILOUT:1); #declare FPP = 1 / (10000); // working variables initialization #local _TIME = 0; #local _DISTANCE = 1; #local _DEVIATION = 2; #local MAX_ERROR = TOL; #local LENGTH = segment_length (SPLINE,0,1,FPP); #debug concat (str(LENGTH,0,8),"\n") #local DATA = array [4] {<-0.5,-0.5,0>,<0,0,0>,<1,1,0>,<1.5,1.5,0>} #local TIME = array [4] {-0.5,0,1,1.5} // process spline #local PASSES = 1; #while (PASSES <= BAILOUT) #debug concat ("pass ",str(PASSES,0,0),"\n") // allocate space for time and data arrays #local OLD_SIZE = dimension_size(TIME,1); #if (OLD_SIZE = 4) #local NEW_SIZE = 500; #else #local NEW_SIZE = (((OLD_SIZE - 3) * 3) + OLD_SIZE); #end #local NEW_DATA = array [NEW_SIZE] #local NEW_TIME = array [NEW_SIZE] // create working K spline #local WORKING = Create_Spline (DATA,TIME) // find statistics #local STEPS = 50; #local tstep = 1 / 50; #local iter = 0; #local pos = Data; #local lpos = Set (0,0,0); #local STATS = array [(STEPS+1)] #while (iter <= STEPS) #local pos = Set ((iter * tstep),0,0); #local F_eval = FUNC_AB (lpos,pos,this_path,WORKING,LENGTH); #local STATS[iter] = abs (F_eval[2]); #local lpos = F_eval; #local iter = iter + 1; #end GetStats (STATS) #declare EXP = 0.3 + (abs (log (StatisticsArray[1]))); #declare FPP = 1 / pow(10,EXP); //#declare FPP = 1 / (10 ^ 2); #debug concat ("FPP : ",str(FPP,2,10),"\n") //#local LIMIT = ( (PASSES = 1) ? -1 : StatisticsArray[1]); #local LIMIT = StatisticsArray[1]; #declare AX[0] = 0; #declare AX = FUNC_A (AX,SPLINE,WORKING,LENGTH); #declare BX[0] = AX[0] + 0.005; #declare BX = FUNC_A (BX,SPLINE,WORKING,LENGTH); // setup variables for this pass #local NEW_TIME [0] = TIME [0]; #local NEW_DATA [0] = DATA [0]; #local NEW_TIME [1] = TIME [1]; #local NEW_DATA [1] = DATA [1]; #local DONE = false; #local NEW_INDEX = 2; #local OLD_INDEX = 2; #while (!DONE) BRACKET(AX,BX,SPLINE,LENGTH,WORKING) #if (BX[0] > 0) BRENT(AX,BX,CX,SPLINE,WORKING,LENGTH,FPP) #if ((abs(BRNT[2])) > (LIMIT*0)) #local This_T = BRNT[1] / LENGTH; #local This_D = BRNT[0]; #while (This_T > TIME[OLD_INDEX]) #local NEW_TIME[NEW_INDEX] = TIME[OLD_INDEX]; #local NEW_DATA[NEW_INDEX] = DATA[OLD_INDEX]; #debug concat("[",str(OLD_INDEX,0,0),"]",str(NEW_INDEX,0,0)," \n") #local NEW_INDEX = NEW_INDEX + 1; #local OLD_INDEX = OLD_INDEX + 1; #end #local NEW_TIME[NEW_INDEX] = This_T; #local NEW_DATA[NEW_INDEX] = ; #debug concat(str(NEW_INDEX,0,0)," \n ") #local NEW_INDEX = NEW_INDEX + 1; #if ((abs(BRNT[2])) > HIGH) #local HIGH = abs(BRNT[2]); #end #end #declare AX = BRNT; #declare BX[0] = AX[0] + 0.005; #declare BX = FUNC_A (BX,SPLINE,WORKING,LENGTH); #else #local DONE = true; #end #end // add last points #while (OLD_INDEX < OLD_SIZE) #local NEW_TIME[ NEW_INDEX ] = TIME [ OLD_INDEX ]; #local NEW_DATA[ NEW_INDEX ] = DATA [ OLD_INDEX ]; #local NEW_INDEX = NEW_INDEX + 1; #local OLD_INDEX = OLD_INDEX + 1; #end // compact and transfer arrays for further refinement #undef DATA #undef TIME #local DATA = array [NEW_INDEX] #local TIME = array [NEW_INDEX] #local this = 0; #while (this < NEW_INDEX) #local DATA[this] = NEW_DATA[this]; #local TIME[this] = NEW_TIME[this]; #local this = this + 1; #end #undef NEW_DATA #undef NEW_TIME #if (HIGH > TOL) #local PASSES = PASSES + 1; #else #local PASSES = BAILOUT + 1; #end #end // OUTPUT OPTIONS // Create Mapping Spline #local this_point = 0; spline { cubic_spline #while (this_point < NEW_INDEX) TIME [ this_point ] , DATA [ this_point ] #debug concat ("TIME [",str(TIME[this_point],3,6),"] , DATA [",str(DATA[this_point].x,3,6),"]\n") #local this_point = this_point + 1; #end } #end #end // CUBIC_SPLINES.INC