1 /****************************************************************************
2 * splines.cpp
3 *
4 * Contributed by ???
5 *
6 * This module implements splines.
7 *
8 * from Persistence of Vision(tm) Ray Tracer version 3.6.
9 * Copyright 1991-2003 Persistence of Vision Team
10 * Copyright 2003-2004 Persistence of Vision Raytracer Pty. Ltd.
11 *---------------------------------------------------------------------------
12 * NOTICE: This source code file is provided so that users may experiment
13 * with enhancements to POV-Ray and to port the software to platforms other
14 * than those supported by the POV-Ray developers. There are strict rules
15 * regarding how you are permitted to use this file. These rules are contained
16 * in the distribution and derivative versions licenses which should have been
17 * provided with this file.
18 *
19 * These licences may be found online, linked from the end-user license
20 * agreement that is located at http://www.povray.org/povlegal.html
21 *---------------------------------------------------------------------------
22 * This program is based on the popular DKB raytracer version 2.12.
23 * DKBTrace was originally written by David K. Buck.
24 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
25 *---------------------------------------------------------------------------
26 * $File: //depot/povray/3.6-release/source/splines.cpp $
27 * $Revision: #3 $
28 * $Change: 3032 $
29 * $DateTime: 2004/08/02 18:43:41 $
30 * $Author: chrisc $
31 * $Log$
32 *****************************************************************************/
33
34 #include "frame.h"
35 #include "userio.h"
36 #include "splines.h"
37
38 BEGIN_POV_NAMESPACE
39
40 /*****************************************************************************
41 * Global Variables
42 ******************************************************************************/
43
44
45 /*****************************************************************************
46 * Local preprocessor defines
47 ******************************************************************************/
48
49
50 /*****************************************************************************
51 * Local typedefs
52 ******************************************************************************/
53
54
55 /*****************************************************************************
56 * Local variables
57 ******************************************************************************/
58
59
60 /*****************************************************************************
61 * Local functions
62 ******************************************************************************/
63
64 DBL linear_interpolate(SPLINE_ENTRY * se, int i, int k, DBL p);
65 DBL quadratic_interpolate(SPLINE_ENTRY * se, int i, int k, DBL p);
66 DBL natural_interpolate(SPLINE_ENTRY * se, int i, int k, DBL p);
67 DBL catmull_rom_interpolate(SPLINE_ENTRY * se, int i, int k, DBL p);
68 int findt(SPLINE * sp, DBL Time);
69 void mkfree(SPLINE * sp, int i);
70 void Precompute_Cubic_Coeffs(SPLINE *sp);
71
72 /*****************************************************************************
73 *
74 * FUNCTION
75 *
76 * Precompute_Cubic_Coeffs
77 *
78 * INPUT
79 *
80 * sp : a pointer to the spline to compute interpolation coefficients for
81 *
82 * OUTPUT
83 *
84 * RETURNS
85 *
86 * AUTHOR
87 *
88 * Mark Wagner
89 *
90 * DESCRIPTION
91 *
92 * Computes the coefficients used in cubic_interpolate.
93 *
94 * CHANGES
95 *
96 ******************************************************************************/
97
Precompute_Cubic_Coeffs(SPLINE * sp)98 void Precompute_Cubic_Coeffs(SPLINE *sp)
99 {
100 int i, k;
101 DBL *h;
102 DBL *b;
103 DBL *u;
104 DBL *v;
105
106 h = (DBL *)POV_MALLOC(sp->Number_Of_Entries*sizeof(DBL), "Spline coefficient storage");
107 b = (DBL *)POV_MALLOC(sp->Number_Of_Entries*sizeof(DBL), "Spline coefficient storage");
108 u = (DBL *)POV_MALLOC(sp->Number_Of_Entries*sizeof(DBL), "Spline coefficient storage");
109 v = (DBL *)POV_MALLOC(sp->Number_Of_Entries*sizeof(DBL), "Spline coefficient storage");
110
111 for(k = 0; k < 5; k++)
112 {
113 for(i = 0; i <= sp->Number_Of_Entries - 2; i++)
114 {
115 h[i] = sp->SplineEntries[i+1].par - sp->SplineEntries[i].par;
116 b[i] = (sp->SplineEntries[i+1].vec[k] - sp->SplineEntries[i].vec[k])/h[i];
117 }
118 u[1] = 2*(h[0]+h[1]);
119 v[1] = 6*(b[1]-b[0]);
120 for(i = 2; i <= sp->Number_Of_Entries - 2; i++)
121 {
122 u[i] = 2*(h[i]+h[i-1]) - (h[i-1]*h[i-1])/u[i-1];
123 v[i] = 6*(b[i]-b[i-1]) - (h[i-1]*v[i-1])/u[i-1];
124 }
125 sp->SplineEntries[sp->Number_Of_Entries-1].coeff[k] = 0;
126 for(i = sp->Number_Of_Entries-2; i > 0; i--)
127 {
128 sp->SplineEntries[i].coeff[k] = (v[i] - h[i]*sp->SplineEntries[i+1].coeff[k])/u[i];
129 }
130 sp->SplineEntries[0].coeff[k] = 0;
131 }
132 sp->Coeffs_Computed = true;
133
134 POV_FREE(h);
135 POV_FREE(b);
136 POV_FREE(u);
137 POV_FREE(v);
138 }
139
140
141 /*****************************************************************************
142 *
143 * FUNCTION
144 *
145 * linear_interpolate
146 *
147 * INPUT
148 *
149 * se : a pointer to the entries in the spline
150 * i : the first point to interpolate between
151 * k : which dimension of the 5D vector to interpolate in
152 * p : the parameter to interpolate the value for
153 *
154 * OUTPUT
155 *
156 * RETURNS
157 *
158 * The value of the kth dimension of the vector linearly interpolated at p
159 *
160 * AUTHOR
161 *
162 * Wolfgang Ortmann
163 *
164 * DESCRIPTION
165 *
166 * CHANGES
167 *
168 ******************************************************************************/
169
linear_interpolate(SPLINE_ENTRY * se,int i,int k,DBL p)170 DBL linear_interpolate(SPLINE_ENTRY * se, int i, int k, DBL p)
171 {
172 DBL p1, p2, v1, v2;
173 p1 = se[i].par;
174 p2 = se[i+1].par;
175 v1 = se[i].vec[k];
176 v2 = se[i+1].vec[k];
177 return (p-p1)*(v2-v1)/(p2-p1)+v1;
178 }
179
180
181 /*****************************************************************************
182 *
183 * FUNCTION
184 *
185 * quadratic_interpolate
186 *
187 * INPUT
188 *
189 * se : a pointer to the entries in the spline
190 * i : the second point of three to interpolate between
191 * k : which dimension of the 5D vector to interpolate in
192 * p : the parameter to interpolate the value for
193 *
194 * OUTPUT
195 *
196 * RETURNS
197 *
198 * The value of the kth dimension of the quadratic interpolation of the
199 * vector at p
200 *
201 * AUTHOR
202 *
203 * Wolfgang Ortmann
204 *
205 * DESCRIPTION
206 *
207 * CHANGES
208 *
209 ******************************************************************************/
210
quadratic_interpolate(SPLINE_ENTRY * se,int i,int k,DBL p)211 DBL quadratic_interpolate(SPLINE_ENTRY * se, int i, int k, DBL p)
212 {
213 /* fit quadratic function to three point*/
214 DBL n;
215 DBL p1, p2, p3,
216 v1, v2, v3,
217 a, b, c;
218 /* assignments to make life easier */
219 p1=se[i-1].par; p2=se[i].par; p3=se[i+1].par;
220 v1=se[i-1].vec[k]; v2=se[i].vec[k]; v3=se[i+1].vec[k];
221
222 n=(p2-p1)*(p3-p1)*(p3-p2);
223
224 /* MWW NOTE: I'm assuming these are correct. I didn't write them */
225 a=(-p2*v1+p3*v1
226 +p1*v2-p3*v2
227 -p1*v3+p2*v3) /n;
228 b=( p2*p2*v1 - p3*p3*v1
229 -p1*p1*v2 + p3*p3*v2
230 +p1*p1*v3 - p2*p2*v3) /n;
231 c=(-p2*p2*p3*v1+p2*p3*p3*v1
232 +p1*p1*p3*v2-p1*p3*p3*v2
233 -p1*p1*p2*v3+p1*p2*p2*v3) /n;
234 return (a*p+b)*p+c; /* Fast way of doing ap^2+bp+c */
235 }
236
237
238 /*****************************************************************************
239 *
240 * FUNCTION
241 *
242 * natural_interpolate
243 *
244 * INPUT
245 *
246 * se : a pointer to the entries in the spline
247 * i : the first point in the interpolation interval
248 * k : which dimension of the 5D vector to interpolate in
249 * p : the parameter to interpolate the value for
250 *
251 * OUTPUT
252 *
253 * RETURNS
254 *
255 * The value of the kth dimension of the natural cubic interpolation
256 * of the vector at p
257 *
258 * AUTHOR
259 *
260 * Mark Wagner
261 *
262 * DESCRIPTION
263 *
264 * CHANGES
265 *
266 ******************************************************************************/
267
natural_interpolate(SPLINE_ENTRY * se,int i,int k,DBL p)268 DBL natural_interpolate(SPLINE_ENTRY * se, int i, int k, DBL p)
269 {
270 DBL h, tmp;
271 h = se[i+1].par - se[i].par;
272 tmp = se[i].coeff[k]/2.0 + ((p - se[i].par)*(se[i+1].coeff[k] - se[i].coeff[k]))/(6.0*h);
273 tmp = -(h/6.0)*(se[i+1].coeff[k] + 2.0*se[i].coeff[k]) + (se[i+1].vec[k] - se[i].vec[k])/h + (p - se[i].par)*tmp;
274 return se[i].vec[k] + (p - se[i].par)*tmp;
275 }
276
277
278 /*****************************************************************************
279 *
280 * FUNCTION
281 *
282 * catmull_rom_interpolate
283 *
284 * INPUT
285 *
286 * se : a pointer to the entries in the spline
287 * i : the first point in the interpolation interval
288 * k : which dimension of the 5D vector to interpolate in
289 * p : the parameter to interpolate the value for
290 *
291 * OUTPUT
292 *
293 * RETURNS
294 *
295 * The value of the kth dimension of the Catmull-Rom interpolation of the
296 * vector at p
297 *
298 * AUTHOR
299 *
300 * Mark Wagner
301 *
302 * DESCRIPTION
303 *
304 * CHANGES
305 *
306 ******************************************************************************/
307
catmull_rom_interpolate(SPLINE_ENTRY * se,int i,int k,DBL p)308 DBL catmull_rom_interpolate(SPLINE_ENTRY * se, int i, int k, DBL p)
309 {
310 DBL dt = (se[i+1].par - se[i].par); /* Time between se[i] and se[i+1] */
311 DBL u = (p - se[i].par)/dt; /* Fractional time from se[i] to p */
312 DBL dp0 = ((se[i].vec[k] - se[i-1].vec[k])/(se[i].par - se[i-1].par) + (se[i+1].vec[k] - se[i].vec[k])/(se[i+1].par - se[i].par))/2.0 * dt;
313 DBL dp1 = ((se[i+2].vec[k] - se[i+1].vec[k])/(se[i+2].par - se[i+1].par) + (se[i+1].vec[k] - se[i].vec[k])/(se[i+1].par - se[i].par))/2.0 * dt;
314
315 return (se[i].vec[k] * (2*u*u*u - 3*u*u + 1) +
316 se[i+1].vec[k] * (3*u*u - 2*u*u*u) +
317 dp0 * (u*u*u - 2*u*u + u) +
318 dp1 * (u*u*u - u*u) );
319 }
320
321
322 /*****************************************************************************
323 *
324 * FUNCTION
325 *
326 * findt
327 *
328 * INPUT
329 *
330 * sp : a pointer to a spline
331 * Time : The parameter to search for
332 *
333 * OUTPUT
334 *
335 * RETURNS
336 *
337 * The first spline entry with a parameter greater than Time
338 *
339 * AUTHOR
340 *
341 * Wolfgang Ortmann
342 *
343 * DESCRIPTION
344 *
345 * CHANGES
346 *
347 * Mark Wagner 6 Nov 2000 : Changed from linear search to binary search
348 *
349 ******************************************************************************/
350
findt(SPLINE * sp,DBL Time)351 int findt(SPLINE * sp, DBL Time)
352 {
353 int i, i2;
354 SPLINE_ENTRY * se;
355 se = sp->SplineEntries;
356 if(sp->Number_Of_Entries == 0) return 0;
357
358 if(Time <= se[0].par) return 0;
359
360 if(Time >= se[sp->Number_Of_Entries-1].par) return sp->Number_Of_Entries;
361
362 i = sp->Number_Of_Entries / 2;
363 /* Bracket the proper entry */
364 if( Time > se[i].par ) /* i is lower, i2 is upper */
365 {
366 i2 = sp->Number_Of_Entries-1;
367 while(i2 - i > 1)
368 {
369 if(Time > se[i+(i2-i)/2].par)
370 i = i+(i2-i)/2;
371 else
372 i2 = i+(i2-i)/2;
373 }
374 return i2;
375 }
376 else /* i is upper, i2 is lower */
377 {
378 i2 = 0;
379 while(i - i2 > 1)
380 {
381 if(Time > se[i2+(i-i2)/2].par)
382 i2 = i2+(i-i2)/2;
383 else
384 i = i2+(i-i2)/2;
385 }
386 return i;
387 }
388 }
389
390
391 /*****************************************************************************
392 *
393 * FUNCTION
394 *
395 * mkfree
396 *
397 * INPUT
398 *
399 * sp : a pointer to the entries in the spline
400 * i : the index of the entry to make available
401 *
402 * OUTPUT
403 *
404 * RETURNS
405 *
406 * AUTHOR
407 *
408 * Wolfgang Ortmann
409 *
410 * DESCRIPTION
411 *
412 * Makes the spline entry at index i available for inserting a new entry
413 * by moving all entries starting with that index down by one slot in the
414 * array.
415 *
416 * CHANGES
417 *
418 ******************************************************************************/
419
mkfree(SPLINE * sp,int i)420 void mkfree(SPLINE * sp, int i)
421 {
422 int j;
423 SPLINE_ENTRY * se;
424 se = sp->SplineEntries;
425
426 for (j=sp->Number_Of_Entries; j>i; j--)
427 se[j] = se[j-1];
428 }
429
430
431 /*****************************************************************************
432 *
433 * FUNCTION
434 *
435 * Create_Spline
436 *
437 * INPUT
438 *
439 * Type : the type of the new spline
440 *
441 * OUTPUT
442 *
443 * RETURNS
444 *
445 * A pointer to the newly-created spline
446 *
447 * AUTHOR
448 *
449 * Wolfgang Ortmann
450 *
451 * DESCRIPTION
452 *
453 * CHANGES
454 *
455 * Mark Wagner 6 Nov 2000 : Added support for dynamic resizing of the
456 * SplineEntry array.
457 * Mark Wagner 25 Aug 2001 : Added support for pre-computing coefficients
458 * of cubic splines
459 *
460 ******************************************************************************/
461
Create_Spline(int Type)462 SPLINE * Create_Spline(int Type)
463 {
464 SPLINE * New;
465 New = (SPLINE *)POV_MALLOC(sizeof(SPLINE), "spline");
466 New->SplineEntries = (SPLINE_ENTRY *)POV_MALLOC(INIT_SPLINE_SIZE*sizeof(SPLINE_ENTRY), "spline entry");
467 New->Max_Entries = INIT_SPLINE_SIZE;
468 New->Number_Of_Entries = 0;
469 New->Type = Type;
470 New->Coeffs_Computed = false;
471 New->Terms = 2;
472 New->Cache_Valid = false;
473 int i;
474 for (i=0; i< New->Max_Entries; i++)
475 {
476 New->SplineEntries[i].par=-1e6; //this should be a negative large number
477 }
478
479 return New;
480 }
481
482
483 /*****************************************************************************
484 *
485 * FUNCTION
486 *
487 * Copy_Spline
488 *
489 * INPUT
490 *
491 * Old : A pointer to the old spline
492 *
493 * OUTPUT
494 *
495 * RETURNS
496 *
497 * A pointer to the new spline
498 *
499 * AUTHOR
500 *
501 * Wolfgang Ortmann
502 *
503 * DESCRIPTION
504 *
505 * CHANGES
506 *
507 * Mark Wagner 6 Nov 2000 : Added support for dynamic resizing of the
508 * SplineEntry array
509 * Mark Wagner 25 Aug 2001 : Added support for pre-computing coefficients
510 * of cubic splines
511 *
512 ******************************************************************************/
513
Copy_Spline(SPLINE * Old)514 SPLINE * Copy_Spline(SPLINE * Old)
515 {
516 SPLINE * New;
517 New = (SPLINE *)POV_MALLOC(sizeof(SPLINE), "spline");
518
519 New->SplineEntries = (SPLINE_ENTRY *)POV_MALLOC(Old->Number_Of_Entries*sizeof(SPLINE_ENTRY), "spline entry");
520 POV_MEMCPY(New->SplineEntries, Old->SplineEntries, Old->Number_Of_Entries*sizeof(SPLINE_ENTRY));
521
522 New->Max_Entries = Old->Number_Of_Entries;
523 New->Number_Of_Entries = Old->Number_Of_Entries;
524 New->Type = Old->Type;
525 New->Coeffs_Computed = Old->Coeffs_Computed;
526 New->Terms = Old->Terms;
527 New->Cache_Valid = false; // we don't copy the cache so mark it as invalid
528
529 return New;
530 }
531
532
533 /*****************************************************************************
534 *
535 * FUNCTION
536 *
537 * Destroy_Spline
538 *
539 * INPUT
540 *
541 * Spline : A pointer to the spline to delete
542 *
543 * OUTPUT
544 *
545 * RETURNS
546 *
547 * AUTHOR
548 *
549 * Wolfgang Ortmann
550 *
551 * DESCRIPTION
552 *
553 * Deletes the SplineEntry array
554 *
555 * CHANGES
556 *
557 ******************************************************************************/
558
Destroy_Spline(SPLINE * Spline)559 void Destroy_Spline(SPLINE * Spline)
560 {
561 POV_FREE(Spline->SplineEntries);
562 POV_FREE(Spline);
563 }
564
565
566 /*****************************************************************************
567 *
568 * FUNCTION
569 *
570 * Insert_Spline_Entry
571 *
572 * INPUT
573 *
574 * sp : a pointer to the spline to insert the new value in
575 * p : the value of the parameter at that point
576 * v : a 5D vector that is the value of the spline at that parameter
577 *
578 * OUTPUT
579 *
580 * RETURNS
581 *
582 * AUTHOR
583 *
584 * Wolfgang Ortmann
585 *
586 * DESCRIPTION
587 *
588 * Inserts a value into the given spline, sorting the entries in order by
589 * increasing p
590 * If a value of the parameter already exists, that value is overwritten
591 *
592 * CHANGES
593 *
594 * Mark Wagner 8 Nov 2000 : Added dynamic sizing of the SplineEntry array
595 * If a value of the parameter already exists, that value is overwritten
596 * Mark Wagner 25 Aug 2001 : Modified for compatibility with new method of
597 * computing cubic splines.
598 * Mark Wagner 30 Aug 2001 : Fixed a bug with over-writing of parameter
599 * values.
600 *
601 ******************************************************************************/
602
Insert_Spline_Entry(SPLINE * sp,DBL p,EXPRESS v)603 void Insert_Spline_Entry(SPLINE * sp, DBL p, EXPRESS v)
604 {
605 int i, k;
606
607 /* Reset the Coeffs_Computed flag. Inserting a new point invalidates
608 * pre-computed coefficients */
609 sp->Coeffs_Computed = false;
610 sp->Cache_Valid = false;
611 /* If all space is used, reallocate */
612 if(sp->Number_Of_Entries >= sp->Max_Entries)
613 {
614 sp->Max_Entries += INIT_SPLINE_SIZE;
615 sp->SplineEntries = (SPLINE_ENTRY *)POV_REALLOC(sp->SplineEntries, sp->Max_Entries * sizeof(SPLINE_ENTRY), "Temporary Spline Entries");
616 for (i = sp->Number_Of_Entries; i < sp->Max_Entries; i++)
617 {
618 sp->SplineEntries[i].par=-1e6;
619 }
620 }
621 i = findt(sp, p);
622 /* If p is already in spline, replace */
623 /* The clause after the || is needed because findt returns sp->Number_Of_Entries
624 * if p is greater than OR EQUAL TO the highest par in the spline */
625 if(sp->Number_Of_Entries != 0 && ((sp->SplineEntries[i].par == p) || (i == sp->Number_Of_Entries && sp->SplineEntries[i-1].par == p)))
626 {
627 for(k=0; k<5; k++)
628 sp->SplineEntries[i].vec[k] = v[k];
629 }
630 else
631 {
632 mkfree(sp, i);
633 sp->SplineEntries[i].par = p;
634
635 for(k=0; k<5; k++)
636 sp->SplineEntries[i].vec[k] = v[k];
637
638 sp->Number_Of_Entries += 1;
639 }
640 }
641
642
643 /*****************************************************************************
644 *
645 * FUNCTION
646 *
647 * Get_Spline_Value
648 *
649 * INPUT
650 *
651 * sp : a pointer to the spline to interpolate
652 * p : the parameter to interpolate the value for
653 * v : a pointer to a 5D vector to store the interpolated values in
654 *
655 * OUTPUT
656 *
657 * v : a pointer to a 5D vector to store the interpolated values in
658 *
659 * RETURNS
660 *
661 * The value of the first element of the interpolated vector
662 *
663 * AUTHOR
664 *
665 * Wolfgang Ortmann
666 *
667 * DESCRIPTION
668 *
669 * Interpolates the spline at the given point. If the point is outside the
670 * range of parameters of the spline, returns the value at the
671 * appropriate endpoint (ie. no extrapolation). If there are not
672 * enough points in the spline to do the desired type of
673 * interpolation, does the next best type (cubic->quadratic->linear).
674 *
675 * CHANGES
676 *
677 * Mark Wagner 8 Nov 2000 : Complete overhaul. I am apalled at the
678 * problems I found in the original implementation
679 * Mark Wagner 25 Aug 2001 : Changed interpolation method to pre-compute
680 * coefficients for cubic splines.
681 * Mark Wagner 24 Feb 2002 : Added support for Catmull-Rom interpolation
682 * Re-arranged the code to make future additions cleaner by moving
683 * more of the code into the "switch" statement
684 *
685 ******************************************************************************/
686
Get_Spline_Val(SPLINE * sp,DBL p,EXPRESS v,int * Terms)687 DBL Get_Spline_Val(SPLINE *sp, DBL p, EXPRESS v, int *Terms)
688 {
689 int i, k;
690 int last;
691 SPLINE_ENTRY * se;
692
693 *Terms = sp->Terms;
694
695 if(!sp->Coeffs_Computed)
696 {
697 switch(sp->Type)
698 {
699 case NATURAL_SPLINE:
700 Precompute_Cubic_Coeffs(sp);
701 break;
702 default:
703 break;
704 }
705 }
706
707 // check if the value is in the cache
708 if((sp->Cache_Point == p) && (sp->Cache_Type == sp->Type))
709 {
710 if(sp->Cache_Valid == true) // doing this here is more efficient as it is rarely false [trf]
711 {
712 Assign_Express(v, sp->Cache_Data);
713 return sp->Cache_Data[0];
714 }
715 }
716
717 // init some cache data
718 sp->Cache_Valid = false;
719 sp->Cache_Type = sp->Type;
720 sp->Cache_Point = p;
721
722 last = sp->Number_Of_Entries-1;
723 se = sp->SplineEntries;
724
725 if(last == 0)
726 {/* if only one entry then return this */
727 for(k=0; k<5; k++)
728 v[k] = se[0].vec[k];
729 return se[0].vec[0];
730 }
731
732 /* Find which spline segment we're in. i is the control point at the end of the segment */
733 i = findt(sp, p);
734
735 switch(sp->Type)
736 {
737 case LINEAR_SPLINE:
738 for(k=0; k<5; k++)
739 {
740 /* If outside spline range, return first or last point */
741 if(i == 0)
742 v[k] = se[0].vec[k];
743 else if(i > last)
744 v[k] = se[last].vec[k];
745 /* Else, normal case */
746 else
747 v[k] = linear_interpolate(se, i-1, k, p);
748 }
749 break;
750 case QUADRATIC_SPLINE:
751 for(k=0; k<5; k++)
752 {
753 /* If outside the spline range, return the first or last point */
754 if(i == 0)
755 v[k] = se[0].vec[k];
756 else if(i > last)
757 v[k] = se[last].vec[k];
758 /* If not enough points, reduce order */
759 else if(last == 1)
760 v[k] = linear_interpolate(se, i-1, k, p);
761 /* Normal case: between the second and last points */
762 else if(i > 1)
763 {
764 v[k] = quadratic_interpolate(se, i-1, k, p);
765 }
766 else /* Special case: between first and second points */
767 {
768 v[k] = quadratic_interpolate(se, i, k, p);
769 }
770 }
771 break;
772 case NATURAL_SPLINE:
773 for(k=0; k<5; k++)
774 {
775 /* If outside the spline range, return the first or last point */
776 if(i == 0)
777 v[k] = se[0].vec[k];
778 else if(i > last)
779 v[k] = se[last].vec[k];
780 /* Else, normal case. cubic_interpolate can handle the case of not enough points */
781 else
782 v[k] = natural_interpolate(se, i-1, k, p);
783 }
784 break;
785 case CATMULL_ROM_SPLINE:
786 for(k=0; k<5; k++)
787 {
788 /* If only two points, return their average */
789 if(last == 1)
790 v[k] = (se[0].vec[k] + se[1].vec[k])/2.0;
791 /* Catmull-Rom: If only three points, return the second one */
792 /* Catmull-Rom: Can't interpolate before second point or after next-to-last */
793 else if(i < 2)
794 v[k] = se[1].vec[k];
795 else if(i >= last)
796 v[k] = se[last-1].vec[k];
797 /* Else, normal case */
798 else
799 v[k] = catmull_rom_interpolate(se, i-1, k, p);
800 }
801 break;
802 default:
803 Error("Unknown spline type %d found.\n", sp->Type);
804
805 }
806
807 // put data in cache
808 Assign_Express(sp->Cache_Data, v);
809 sp->Cache_Valid = true;
810
811 return v[0];
812 }
813
814 END_POV_NAMESPACE
815