1 //******************************************************************************
2 ///
3 /// @file core/shape/polynomial.cpp
4 ///
5 /// Implementation of the 3-variable polynomial geometric primitive.
6 ///
7 /// @author Alexander Enzmann
8 ///
9 /// @copyright
10 /// @parblock
11 ///
12 /// Persistence of Vision Ray Tracer ('POV-Ray') version 3.8.
13 /// Copyright 1991-2018 Persistence of Vision Raytracer Pty. Ltd.
14 ///
15 /// POV-Ray is free software: you can redistribute it and/or modify
16 /// it under the terms of the GNU Affero General Public License as
17 /// published by the Free Software Foundation, either version 3 of the
18 /// License, or (at your option) any later version.
19 ///
20 /// POV-Ray is distributed in the hope that it will be useful,
21 /// but WITHOUT ANY WARRANTY; without even the implied warranty of
22 /// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 /// GNU Affero General Public License for more details.
24 ///
25 /// You should have received a copy of the GNU Affero General Public License
26 /// along with this program. If not, see <http://www.gnu.org/licenses/>.
27 ///
28 /// ----------------------------------------------------------------------------
29 ///
30 /// POV-Ray is based on the popular DKB raytracer version 2.12.
31 /// DKBTrace was originally written by David K. Buck.
32 /// DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
33 ///
34 /// @endparblock
35 ///
36 //******************************************************************************
37
38 // Unit header file must be the first file included within POV-Ray *.cpp files (pulls in config)
39 #include "core/shape/polynomial.h"
40
41 #include "core/bounding/boundingbox.h"
42 #include "core/math/matrix.h"
43 #include "core/math/polynomialsolver.h"
44 #include "core/render/ray.h"
45 #include "core/scene/tracethreaddata.h"
46
47 // this must be the last file included
48 #include "base/povdebug.h"
49
50 /*
51 * Basic form of a quartic equation:
52 *
53 * a00*x^4 + a01*x^3*y + a02*x^3*z + a03*x^3 + a04*x^2*y^2 +
54 * a05*x^2*y*z+ a06*x^2*y + a07*x^2*z^2 + a08*x^2*z + a09*x^2 +
55 * a10*x*y^3 + a11*x*y^2*z + a12*x*y^2 + a13*x*y*z^2 + a14*x*y*z +
56 * a15*x*y + a16*x*z^3 + a17*x*z^2 + a18*x*z + a19*x + a20*y^4 +
57 * a21*y^3*z + a22*y^3 + a23*y^2*z^2 + a24*y^2*z + a25*y^2 + a26*y*z^3 +
58 * a27*y*z^2 + a28*y*z + a29*y + a30*z^4 + a31*z^3 + a32*z^2 + a33*z + a34
59 *
60 */
61
62 namespace pov
63 {
64
65 /*****************************************************************************
66 * Local preprocessor defines
67 ******************************************************************************/
68
69 const DBL DEPTH_TOLERANCE = 1.0e-4;
70 const DBL INSIDE_TOLERANCE = 1.0e-4;
71 const DBL ROOT_TOLERANCE = 1.0e-4;
72 const DBL COEFF_LIMIT = 1.0e-20;
73 // const int BINOMSIZE = 40;
74
75
76
77 /*****************************************************************************
78 * Local variables
79 ******************************************************************************/
80
81 /* The following table contains the binomial coefficients up to 35 */
82 #if (MAX_ORDER > 35)
83 #error "Poly.cpp code would break loose due to a too short pascal triangle table."
84 #endif
85 /* this is a pascal's triangle : [k][j]=[k-1][j-1]+[k-1][j] ;
86 [k][0]=1, [k][k]=1
87 */
88 /* Max value in [35] is 0x8B18014C, so unsigned int is enough
89 * If you want to go down to from [36] to [69],
90 * a 64 bits unsigned long must be used
91 */
92 const unsigned int binomials[35][35] =
93 {
94
95 { 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
96 { 1u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
97 { 1u, 2u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
98 { 1u, 3u, 3u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
99 { 1u, 4u, 6u, 4u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
100 { 1u, 5u, 10u, 10u, 5u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
101 { 1u, 6u, 15u, 20u, 15u, 6u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
102 { 1u, 7u, 21u, 35u, 35u, 21u, 7u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
103 { 1u, 8u, 28u, 56u, 70u, 56u, 28u, 8u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
104 { 1u, 9u, 36u, 84u, 126u, 126u, 84u, 36u, 9u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
105 { 1u, 10u, 45u, 120u, 210u, 252u, 210u, 120u, 45u, 10u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
106 { 1u, 11u, 55u, 165u, 330u, 462u, 462u, 330u, 165u, 55u, 11u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
107 { 1u, 12u, 66u, 220u, 495u, 792u, 924u, 792u, 495u, 220u, 66u, 12u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
108 { 1u, 13u, 78u, 286u, 715u, 1287u, 1716u, 1716u, 1287u, 715u, 286u, 78u, 13u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
109 { 1u, 14u, 91u, 364u, 1001u, 2002u, 3003u, 3432u, 3003u, 2002u, 1001u, 364u, 91u, 14u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
110 { 1u, 15u, 105u, 455u, 1365u, 3003u, 5005u, 6435u, 6435u, 5005u, 3003u, 1365u, 455u, 105u, 15u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
111 { 1u, 16u, 120u, 560u, 1820u, 4368u, 8008u, 11440u, 12870u, 11440u, 8008u, 4368u, 1820u, 560u, 120u, 16u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
112 { 1u, 17u, 136u, 680u, 2380u, 6188u, 12376u, 19448u, 24310u, 24310u, 19448u, 12376u, 6188u, 2380u, 680u, 136u, 17u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
113 { 1u, 18u, 153u, 816u, 3060u, 8568u, 18564u, 31824u, 43758u, 48620u, 43758u, 31824u, 18564u, 8568u, 3060u, 816u, 153u, 18u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
114 { 1u, 19u, 171u, 969u, 3876u, 11628u, 27132u, 50388u, 75582u, 92378u, 92378u, 75582u, 50388u, 27132u, 11628u, 3876u, 969u, 171u, 19u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
115 { 1u, 20u, 190u, 1140u, 4845u, 15504u, 38760u, 77520u, 125970u, 167960u, 184756u, 167960u, 125970u, 77520u, 38760u, 15504u, 4845u, 1140u, 190u, 20u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
116 { 1u, 21u, 210u, 1330u, 5985u, 20349u, 54264u, 116280u, 203490u, 293930u, 352716u, 352716u, 293930u, 203490u, 116280u, 54264u, 20349u, 5985u, 1330u, 210u, 21u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
117 { 1u, 22u, 231u, 1540u, 7315u, 26334u, 74613u, 170544u, 319770u, 497420u, 646646u, 705432u, 646646u, 497420u, 319770u, 170544u, 74613u, 26334u, 7315u, 1540u, 231u, 22u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
118 { 1u, 23u, 253u, 1771u, 8855u, 33649u, 100947u, 245157u, 490314u, 817190u, 1144066u, 1352078u, 1352078u, 1144066u, 817190u, 490314u, 245157u, 100947u, 33649u, 8855u, 1771u, 253u, 23u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
119 { 1u, 24u, 276u, 2024u, 10626u, 42504u, 134596u, 346104u, 735471u, 1307504u, 1961256u, 2496144u, 2704156u, 2496144u, 1961256u, 1307504u, 735471u, 346104u, 134596u, 42504u, 10626u, 2024u, 276u, 24u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
120 { 1u, 25u, 300u, 2300u, 12650u, 53130u, 177100u, 480700u, 1081575u, 2042975u, 3268760u, 4457400u, 5200300u, 5200300u, 4457400u, 3268760u, 2042975u, 1081575u, 480700u, 177100u, 53130u, 12650u, 2300u, 300u, 25u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
121 { 1u, 26u, 325u, 2600u, 14950u, 65780u, 230230u, 657800u, 1562275u, 3124550u, 5311735u, 7726160u, 9657700u, 10400600u, 9657700u, 7726160u, 5311735u, 3124550u, 1562275u, 657800u, 230230u, 65780u, 14950u, 2600u, 325u, 26u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
122 { 1u, 27u, 351u, 2925u, 17550u, 80730u, 296010u, 888030u, 2220075u, 4686825u, 8436285u, 13037895u, 17383860u, 20058300u, 20058300u, 17383860u, 13037895u, 8436285u, 4686825u, 2220075u, 888030u, 296010u, 80730u, 17550u, 2925u, 351u, 27u, 1u, 0u, 0u, 0u, 0u, 0u, 0u, 0u},
123 { 1u, 28u, 378u, 3276u, 20475u, 98280u, 376740u, 1184040u, 3108105u, 6906900u, 13123110u, 21474180u, 30421755u, 37442160u, 40116600u, 37442160u, 30421755u, 21474180u, 13123110u, 6906900u, 3108105u, 1184040u, 376740u, 98280u, 20475u, 3276u, 378u, 28u, 1u, 0u, 0u, 0u, 0u, 0u, 0u},
124 { 1u, 29u, 406u, 3654u, 23751u, 118755u, 475020u, 1560780u, 4292145u, 10015005u, 20030010u, 34597290u, 51895935u, 67863915u, 77558760u, 77558760u, 67863915u, 51895935u, 34597290u, 20030010u, 10015005u, 4292145u, 1560780u, 475020u, 118755u, 23751u, 3654u, 406u, 29u, 1u, 0u, 0u, 0u, 0u, 0u},
125 { 1u, 30u, 435u, 4060u, 27405u, 142506u, 593775u, 2035800u, 5852925u, 14307150u, 30045015u, 54627300u, 86493225u, 119759850u, 145422675u, 155117520u, 145422675u, 119759850u, 86493225u, 54627300u, 30045015u, 14307150u, 5852925u, 2035800u, 593775u, 142506u, 27405u, 4060u, 435u, 30u, 1u, 0u, 0u, 0u, 0u},
126 { 1u, 31u, 465u, 4495u, 31465u, 169911u, 736281u, 2629575u, 7888725u, 20160075u, 44352165u, 84672315u, 141120525u, 206253075u, 265182525u, 300540195u, 300540195u, 265182525u, 206253075u, 141120525u, 84672315u, 44352165u, 20160075u, 7888725u, 2629575u, 736281u, 169911u, 31465u, 4495u, 465u, 31u, 1u, 0u, 0u, 0u},
127 { 1u, 32u, 496u, 4960u, 35960u, 201376u, 906192u, 3365856u, 10518300u, 28048800u, 64512240u, 129024480u, 225792840u, 347373600u, 471435600u, 565722720u, 601080390u, 565722720u, 471435600u, 347373600u, 225792840u, 129024480u, 64512240u, 28048800u, 10518300u, 3365856u, 906192u, 201376u, 35960u, 4960u, 496u, 32u, 1u, 0u, 0u},
128 { 1u, 33u, 528u, 5456u, 40920u, 237336u, 1107568u, 4272048u, 13884156u, 38567100u, 92561040u, 193536720u, 354817320u, 573166440u, 818809200u, 1037158320u, 1166803110u, 1166803110u, 1037158320u, 818809200u, 573166440u, 354817320u, 193536720u, 92561040u, 38567100u, 13884156u, 4272048u, 1107568u, 237336u, 40920u, 5456u, 528u, 33u, 1u, 0u},
129 { 1u, 34u, 561u, 5984u, 46376u, 278256u, 1344904u, 5379616u, 18156204u, 52451256u, 131128140u, 286097760u, 548354040u, 927983760u, 1391975640u, 1855967520u, 2203961430u, 2333606220u, 2203961430u, 1855967520u, 1391975640u, 927983760u, 548354040u, 286097760u, 131128140u, 52451256u, 18156204u, 5379616u, 1344904u, 278256u, 46376u, 5984u, 561u, 34u, 1u}
130
131 };
132
133 /*
134 */
Set_Coeff(const unsigned int x,const unsigned int y,const unsigned int z,const DBL value)135 bool Poly::Set_Coeff(const unsigned int x, const unsigned int y, const unsigned int z, const DBL value)
136 {
137 int a,b,c;
138 unsigned int pos;
139 a=Order-x;
140 b=Order-x-y;
141 c=Order-x-y-z;
142 /* a bit overprotective */
143 if ((x+y+z>Order)||(a<0)||(b<0)||(c<0))
144 {
145 return false;
146 }
147 /* pos = binomials[a+2][3]+binomials[b+1][2]+binomials[c][1];
148 * rewritten to stay in bound (max is "Order", not Order+2)
149 */
150 pos =
151 // binomials[a+2][3]
152 // binomials[a+1][3]
153 binomials[a][3]+binomials[a][2]
154 // +binomials[a+1][2]
155 +binomials[a][2]+binomials[a][1]
156 // +binomials[b+1][2]
157 +binomials[b][1]+binomials[b][2]
158 +binomials[c][1];
159 /* It's magic
160 * Nah... a is the tetrahedral sum to jump to get to the power of x index (first entry)
161 * b is then the triangular sum to add to get to the power of y index (also first entry)
162 * and c is the linear sum to add to get to the power of z index (that the one we want)
163 *
164 * Notice that binomials[c][1] == c, but the formula would loose its magic use of
165 * pascal triangle everywhere.
166 * triangular sums are in the third ([2] column)
167 * tetrahedral sums are in the fourth ([3] column)
168 *
169 * (and yes, the 0 at the start of each column is useful)
170 */
171 Coeffs[pos] = value;
172 return true;
173 }
174
175
176 /*****************************************************************************
177 *
178 * FUNCTION
179 *
180 * All_Poly_Intersections
181 *
182 * INPUT
183 *
184 * OUTPUT
185 *
186 * RETURNS
187 *
188 * AUTHOR
189 *
190 * Alexander Enzmann
191 *
192 * DESCRIPTION
193 *
194 * -
195 *
196 * CHANGES
197 *
198 * -
199 *
200 ******************************************************************************/
201
All_Intersections(const Ray & ray,IStack & Depth_Stack,TraceThreadData * Thread)202 bool Poly::All_Intersections(const Ray& ray, IStack& Depth_Stack, TraceThreadData *Thread)
203 {
204 DBL Depths[MAX_ORDER];
205 DBL len;
206 Vector3d IPoint;
207 int cnt, i, j, Intersection_Found, same_root;
208 BasicRay New_Ray;
209
210 /* Transform the ray into the polynomial's space */
211
212 MInvTransRay(New_Ray, ray, Trans);
213
214 len = New_Ray.Direction.length();
215 New_Ray.Direction /= len;
216
217 Intersection_Found = false;
218
219 Thread->Stats()[Ray_Poly_Tests]++;
220
221 switch (Order)
222 {
223 case 1:
224
225 cnt = intersect_linear(New_Ray, Coeffs, Depths);
226
227 break;
228
229 case 2:
230
231 cnt = intersect_quadratic(New_Ray, Coeffs, Depths);
232
233 break;
234
235 default:
236
237 cnt = intersect(New_Ray, Order, Coeffs, Test_Flag(this, STURM_FLAG), Depths, Thread);
238 }
239
240 if (cnt > 0)
241 {
242 Thread->Stats()[Ray_Poly_Tests_Succeeded]++;
243 }
244
245 for (i = 0; i < cnt; i++)
246 {
247 if (Depths[i] > DEPTH_TOLERANCE)
248 {
249 same_root = false;
250
251 for (j = 0; j < i; j++)
252 {
253 if (Depths[i] == Depths[j])
254 {
255 same_root = true;
256
257 break;
258 }
259 }
260
261 if (!same_root)
262 {
263 IPoint = New_Ray.Evaluate(Depths[i]);
264
265 /* Transform the point into world space */
266
267 MTransPoint(IPoint, IPoint, Trans);
268
269 if (Clip.empty() || Point_In_Clip(IPoint, Clip, Thread))
270 {
271 Depth_Stack->push(Intersection(Depths[i] / len,IPoint,this));
272
273 Intersection_Found = true;
274 }
275 }
276 }
277 }
278
279 return (Intersection_Found);
280 }
281
282
283
284 /*****************************************************************************
285 *
286 * FUNCTION
287 *
288 * evaluate_linear
289 *
290 * INPUT
291 *
292 * OUTPUT
293 *
294 * RETURNS
295 *
296 * AUTHOR
297 *
298 * Alexander Enzmann
299 *
300 * DESCRIPTION
301 *
302 * -
303 *
304 * CHANGES
305 *
306 * -
307 *
308 ******************************************************************************/
309
310 /* For speedup of low order polynomials, expand out the terms
311 involved in evaluating the poly. */
312 /* unused
313 DBL evaluate_linear(const Vector3d& P, DBL *a)
314 {
315 return(a[0] * P[X]) + (a[1] * P[Y]) + (a[2] * P[Z]) + a[3];
316 }
317 */
318
319
320
321 /*****************************************************************************
322 *
323 * FUNCTION
324 *
325 * evaluate_quadratic
326 *
327 * INPUT
328 *
329 * OUTPUT
330 *
331 * RETURNS
332 *
333 * AUTHOR
334 *
335 * Alexander Enzmann
336 *
337 * DESCRIPTION
338 *
339 * -
340 *
341 * CHANGES
342 *
343 * -
344 *
345 ******************************************************************************/
346
347 /*
348 DBL evaluate_quadratic(const Vector3d& P, DBL *a)
349 {
350 DBL x, y, z;
351
352 x = P[X];
353 y = P[Y];
354 z = P[Z];
355
356 return(a[0] * x * x + a[1] * x * y + a[2] * x * z +
357 a[3] * x + a[4] * y * y + a[5] * y * z +
358 a[6] * y + a[7] * z * z + a[8] * z + a[9]);
359 }
360 */
361
362
363
364 /*****************************************************************************
365 *
366 * FUNCTION
367 *
368 * factor_out
369 *
370 * INPUT
371 *
372 * OUTPUT
373 *
374 * RETURNS
375 *
376 * AUTHOR
377 *
378 * Alexander Enzmann
379 *
380 * DESCRIPTION
381 *
382 * Remove all factors of i from n.
383 *
384 * CHANGES
385 *
386 * -
387 *
388 ******************************************************************************/
389 /*
390 int Poly::factor_out(int n, int i, int *c, int *s)
391 {
392 while (!(n % i))
393 {
394 n /= i;
395
396 s[(*c)++] = i;
397 }
398
399 return(n);
400 }
401 */
402
403
404 /*****************************************************************************
405 *
406 * FUNCTION
407 *
408 * factor1
409 *
410 * INPUT
411 *
412 * OUTPUT
413 *
414 * RETURNS
415 *
416 * AUTHOR
417 *
418 * Alexander Enzmann
419 *
420 * DESCRIPTION
421 *
422 * Find all prime factors of n. (Note that n must be less than 2^15.
423 *
424 * CHANGES
425 *
426 * -
427 *
428 ******************************************************************************/
429 #if 0
430 void Poly::factor1(int n, int *c, int *s)
431 {
432 int i,k;
433
434 /* First factor out any 2s. */
435
436 n = factor_out(n, 2, c, s);
437
438 /* Now any odd factors. */
439
440 k = (int)sqrt((DBL)n) + 1;
441
442 for (i = 3; (n > 1) && (i <= k); i += 2)
443 {
444 if (!(n % i))
445 {
446 n = factor_out(n, i, c, s);
447 k = (int)sqrt((DBL)n)+1;
448 }
449 }
450
451 if (n > 1)
452 {
453 s[(*c)++] = n;
454 }
455 }
456 #endif
457
458
459 /*****************************************************************************
460 *
461 * FUNCTION
462 *
463 * binomial
464 *
465 * INPUT
466 *
467 * OUTPUT
468 *
469 * RETURNS
470 *
471 * AUTHOR
472 *
473 * Alexander Enzmann
474 *
475 * DESCRIPTION
476 *
477 * Calculate the binomial coefficent of n, r.
478 *
479 * CHANGES
480 *
481 * -
482 *
483 ******************************************************************************/
484 #if 0
485 int Poly::binomial(int n, int r)
486 {
487 int h,i,j,k,l;
488 unsigned int result;
489 int stack1[BINOMSIZE], stack2[BINOMSIZE];
490
491 if ((n < 0) || (r < 0) || (r > n))
492 {
493 result = 0L;
494 }
495 else
496 {
497 if (r == n)
498 {
499 result = 1L;
500 }
501 else
502 {
503 if ((r < 16) && (n < 16))
504 {
505 result = (int)binomials[n][r];
506 }
507 else
508 {
509 j = 0;
510
511 for (i = r + 1; i <= n; i++)
512 {
513 stack1[j++] = i;
514 }
515
516 for (i = 2; i <= (n-r); i++)
517 {
518 h = 0;
519
520 factor1(i, &h, stack2);
521
522 for (k = 0; k < h; k++)
523 {
524 for (l = 0; l < j; l++)
525 {
526 if (!(stack1[l] % stack2[k]))
527 {
528 stack1[l] /= stack2[k];
529
530 goto l1;
531 }
532 }
533 }
534
535 // Error if we get here
536 // Debug_Info("Failed to factor\n");
537 l1:;
538 }
539
540 result = 1;
541
542 for (i = 0; i < j; i++)
543 {
544 result *= stack1[i];
545 }
546 }
547 }
548 }
549
550 return(result);
551 }
552 #endif
553
554
555 /*****************************************************************************
556 *
557 * FUNCTION
558 *
559 * inside
560 *
561 * INPUT
562 *
563 * OUTPUT
564 *
565 * RETURNS
566 *
567 * AUTHOR
568 *
569 * Alexander Enzmann
570 *
571 * DESCRIPTION
572 *
573 * -
574 *
575 * CHANGES
576 *
577 * -
578 *
579 ******************************************************************************/
580
inside(const Vector3d & IPoint,int Order,const DBL * Coeffs)581 DBL Poly::inside(const Vector3d& IPoint, int Order, const DBL *Coeffs)
582 {
583 DBL x[MAX_ORDER+1], y[MAX_ORDER+1], z[MAX_ORDER+1];
584 DBL c, Result;
585 int i, j, k, term;
586
587 x[0] = 1.0; y[0] = 1.0; z[0] = 1.0;
588 x[1] = IPoint[X]; y[1] = IPoint[Y]; z[1] = IPoint[Z];
589
590 for (i = 2; i <= Order; i++)
591 {
592 x[i] = x[1] * x[i-1];
593 y[i] = y[1] * y[i-1];
594 z[i] = z[1] * z[i-1];
595 }
596
597 Result = 0.0;
598
599 term = 0;
600
601 for (i = Order; i >= 0; i--)
602 {
603 for (j=Order-i;j>=0;j--)
604 {
605 for (k=Order-(i+j);k>=0;k--)
606 {
607 if ((c = Coeffs[term]) != 0.0)
608 {
609 Result += c * x[i] * y[j] * z[k];
610 }
611 term++;
612 }
613 }
614 }
615
616 return(Result);
617 }
618
619
620
621 /*****************************************************************************
622 *
623 * FUNCTION
624 *
625 * intersect
626 *
627 * INPUT
628 *
629 * OUTPUT
630 *
631 * RETURNS
632 *
633 * AUTHOR
634 *
635 * Alexander Enzmann
636 *
637 * DESCRIPTION
638 *
639 * Intersection of a ray and an arbitrary polynomial function.
640 *
641 * CHANGES
642 *
643 * -
644 *
645 ******************************************************************************/
646
intersect(const BasicRay & ray,int Order,const DBL * Coeffs,int Sturm_Flag,DBL * Depths,TraceThreadData * Thread)647 int Poly::intersect(const BasicRay &ray, int Order, const DBL *Coeffs, int Sturm_Flag, DBL *Depths, TraceThreadData *Thread)
648 {
649 DBL eqn_v[3][MAX_ORDER+1], eqn_vt[3][MAX_ORDER+1];
650 DBL eqn[MAX_ORDER+1];
651 DBL t[3][MAX_ORDER+1];
652 Vector3d P, D;
653 DBL val;
654 int h, i, j, k, i1, j1, k1, term;
655 int offset;
656
657 /* First we calculate the values of the individual powers
658 of x, y, and z as they are represented by the ray */
659
660 P = ray.Origin;
661 D = ray.Direction;
662
663 for (i = 0; i < 3; i++)
664 {
665 eqn_v[i][0] = 1.0;
666 eqn_vt[i][0] = 1.0;
667 }
668
669 eqn_v[0][1] = P[X];
670 eqn_v[1][1] = P[Y];
671 eqn_v[2][1] = P[Z];
672
673 eqn_vt[0][1] = D[X];
674 eqn_vt[1][1] = D[Y];
675 eqn_vt[2][1] = D[Z];
676
677 for (i = 2; i <= Order; i++)
678 {
679 for (j=0;j<3;j++)
680 {
681 eqn_v[j][i] = eqn_v[j][1] * eqn_v[j][i-1];
682 eqn_vt[j][i] = eqn_vt[j][1] * eqn_vt[j][i-1];
683 }
684 }
685
686 for (i = 0; i <= Order; i++)
687 {
688 eqn[i] = 0.0;
689 }
690
691 /* Now walk through the terms of the polynomial. As we go
692 we substitute the ray equation for each of the variables. */
693
694 term = 0;
695
696 for (i = Order; i >= 0; i--)
697 {
698 for (h = 0; h <= i; h++)
699 {
700 t[0][h] = binomials[i][h] * eqn_vt[0][i-h] * eqn_v[0][h];
701 }
702
703 for (j = Order-i; j >= 0; j--)
704 {
705 for (h = 0; h <= j; h++)
706 {
707 t[1][h] = binomials[j][h] * eqn_vt[1][j-h] * eqn_v[1][h];
708 }
709
710 for (k = Order-(i+j); k >= 0; k--)
711 {
712 if (Coeffs[term] != 0)
713 {
714 for (h = 0; h <= k; h++)
715 {
716 t[2][h] = binomials[k][h] * eqn_vt[2][k-h] * eqn_v[2][h];
717 }
718
719 /* Multiply together the three polynomials. */
720
721 offset = Order - (i + j + k);
722
723 for (i1 = 0; i1 <= i; i1++)
724 {
725 for (j1=0;j1<=j;j1++)
726 {
727 for (k1=0;k1<=k;k1++)
728 {
729 h = offset + i1 + j1 + k1;
730 val = Coeffs[term];
731 val *= t[0][i1];
732 val *= t[1][j1];
733 val *= t[2][k1];
734 eqn[h] += val;
735 }
736 }
737 }
738 }
739
740 term++;
741 }
742 }
743 }
744
745 for (i = 0, j = Order; i <= Order; i++)
746 {
747 if (eqn[i] != 0.0)
748 {
749 break;
750 }
751 else
752 {
753 j--;
754 }
755 }
756
757 if (j > 1)
758 {
759 return(Solve_Polynomial(j, &eqn[i], Depths, Sturm_Flag, ROOT_TOLERANCE, Thread->Stats()));
760 }
761 else
762 {
763 return(0);
764 }
765 }
766
767
768
769 /*****************************************************************************
770 *
771 * FUNCTION
772 *
773 * intersect_linear
774 *
775 * INPUT
776 *
777 * OUTPUT
778 *
779 * RETURNS
780 *
781 * AUTHOR
782 *
783 * Alexander Enzmann
784 *
785 * DESCRIPTION
786 *
787 * Intersection of a ray and a quadratic.
788 *
789 * CHANGES
790 *
791 * -
792 *
793 ******************************************************************************/
794
intersect_linear(const BasicRay & ray,const DBL * Coeffs,DBL * Depths)795 int Poly::intersect_linear(const BasicRay &ray, const DBL *Coeffs, DBL *Depths)
796 {
797 DBL t0, t1;
798 const DBL *a = Coeffs;
799
800 t0 = a[0] * ray.Origin[X] + a[1] * ray.Origin[Y] + a[2] * ray.Origin[Z];
801 t1 = a[0] * ray.Direction[X] + a[1] * ray.Direction[Y] +
802
803 a[2] * ray.Direction[Z];
804
805 if (fabs(t1) < EPSILON)
806 {
807 return(0);
808 }
809
810 Depths[0] = -(a[3] + t0) / t1;
811
812 return(1);
813 }
814
815
816
817 /*****************************************************************************
818 *
819 * FUNCTION
820 *
821 * intersect_quadratic
822 *
823 * INPUT
824 *
825 * OUTPUT
826 *
827 * RETURNS
828 *
829 * AUTHOR
830 *
831 * Alexander Enzmann
832 *
833 * DESCRIPTION
834 *
835 * Intersection of a ray and a quadratic.
836 *
837 * CHANGES
838 *
839 * -
840 *
841 ******************************************************************************/
842
intersect_quadratic(const BasicRay & ray,const DBL * Coeffs,DBL * Depths)843 int Poly::intersect_quadratic(const BasicRay &ray, const DBL *Coeffs, DBL *Depths)
844 {
845 DBL x, y, z, x2, y2, z2;
846 DBL xx, yy, zz, xx2, yy2, zz2, ac, bc, cc, d, t;
847 const DBL *a;
848
849 x = ray.Origin[X];
850 y = ray.Origin[Y];
851 z = ray.Origin[Z];
852
853 xx = ray.Direction[X];
854 yy = ray.Direction[Y];
855 zz = ray.Direction[Z];
856
857 x2 = x * x; y2 = y * y; z2 = z * z;
858 xx2 = xx * xx; yy2 = yy * yy; zz2 = zz * zz;
859
860 a = Coeffs;
861
862 /*
863 * Determine the coefficients of t^n, where the line is represented
864 * as (x,y,z) + (xx,yy,zz)*t.
865 */
866
867 ac = (a[0]*xx2 + a[1]*xx*yy + a[2]*xx*zz + a[4]*yy2 + a[5]*yy*zz + a[7]*zz2);
868
869 bc = (2*a[0]*x*xx + a[1]*(x*yy + xx*y) + a[2]*(x*zz + xx*z) +
870 a[3]*xx + 2*a[4]*y*yy + a[5]*(y*zz + yy*z) + a[6]*yy +
871 2*a[7]*z*zz + a[8]*zz);
872
873 cc = a[0]*x2 + a[1]*x*y + a[2]*x*z + a[3]*x + a[4]*y2 +
874 a[5]*y*z + a[6]*y + a[7]*z2 + a[8]*z + a[9];
875
876 if (fabs(ac) < COEFF_LIMIT)
877 {
878 if (fabs(bc) < COEFF_LIMIT)
879 {
880 return(0);
881 }
882
883 Depths[0] = -cc / bc;
884
885 return(1);
886 }
887
888 /*
889 * Solve the quadratic formula & return results that are
890 * within the correct interval.
891 */
892
893 d = bc * bc - 4.0 * ac * cc;
894
895 if (d < 0.0)
896 {
897 return(0);
898 }
899
900 d = sqrt(d);
901
902 bc = -bc;
903
904 t = 2.0 * ac;
905
906 Depths[0] = (bc + d) / t;
907 Depths[1] = (bc - d) / t;
908
909 return(2);
910 }
911
912
913
914 /*****************************************************************************
915 *
916 * FUNCTION
917 *
918 * normal0
919 *
920 * INPUT
921 *
922 * OUTPUT
923 *
924 * RETURNS
925 *
926 * AUTHOR
927 *
928 * Alexander Enzmann
929 *
930 * DESCRIPTION
931 *
932 * Normal to a polynomial (used for polynomials with order > 4).
933 *
934 * CHANGES
935 *
936 * -
937 *
938 ******************************************************************************/
939
normal0(Vector3d & Result,int Order,const DBL * Coeffs,const Vector3d & IPoint)940 void Poly::normal0(Vector3d& Result, int Order, const DBL *Coeffs, const Vector3d& IPoint)
941 {
942 int i, j, k, term;
943 DBL x[MAX_ORDER+1], y[MAX_ORDER+1], z[MAX_ORDER+1];
944 DBL val;
945 const DBL *a;
946
947 x[0] = 1.0;
948 y[0] = 1.0;
949 z[0] = 1.0;
950
951 x[1] = IPoint[X];
952 y[1] = IPoint[Y];
953 z[1] = IPoint[Z];
954
955 for (i = 2; i <= Order; i++)
956 {
957 x[i] = IPoint[X] * x[i-1];
958 y[i] = IPoint[Y] * y[i-1];
959 z[i] = IPoint[Z] * z[i-1];
960 }
961
962 a = Coeffs;
963
964 term = 0;
965
966 Result = Vector3d(0.0, 0.0, 0.0);
967
968 for (i = Order; i >= 0; i--)
969 {
970 for (j = Order-i; j >= 0; j--)
971 {
972 for (k = Order-(i+j); k >= 0; k--)
973 {
974 if (i >= 1)
975 {
976 val = x[i-1] * y[j] * z[k];
977 Result[X] += i * a[term] * val;
978 }
979
980 if (j >= 1)
981 {
982 val = x[i] * y[j-1] * z[k];
983 Result[Y] += j * a[term] * val;
984 }
985
986 if (k >= 1)
987 {
988 val = x[i] * y[j] * z[k-1];
989 Result[Z] += k * a[term] * val;
990 }
991
992 term++;
993 }
994 }
995 }
996 }
997
998
999
1000 /*****************************************************************************
1001 *
1002 * FUNCTION
1003 *
1004 * nromal1
1005 *
1006 * INPUT
1007 *
1008 * OUTPUT
1009 *
1010 * RETURNS
1011 *
1012 * AUTHOR
1013 *
1014 * Alexander Enzmann
1015 *
1016 * DESCRIPTION
1017 *
1018 * Normal to a polynomial (for polynomials of order <= 4).
1019 *
1020 * CHANGES
1021 *
1022 * -
1023 *
1024 ******************************************************************************/
1025
normal1(Vector3d & Result,int Order,const DBL * Coeffs,const Vector3d & IPoint)1026 void Poly::normal1(Vector3d& Result, int Order, const DBL *Coeffs, const Vector3d& IPoint)
1027 {
1028 DBL x, y, z, x2, y2, z2, x3, y3, z3;
1029 const DBL *a;
1030
1031 a = Coeffs;
1032
1033 x = IPoint[X];
1034 y = IPoint[Y];
1035 z = IPoint[Z];
1036
1037 switch (Order)
1038 {
1039 case 1:
1040
1041 /* Linear partial derivatives */
1042
1043 Result = Vector3d(a[0], a[1], a[2]);
1044
1045 break;
1046
1047 case 2:
1048
1049 /* Quadratic partial derivatives */
1050
1051 Result[X] = 2*a[0]*x+a[1]*y+a[2]*z+a[3];
1052 Result[Y] = a[1]*x+2*a[4]*y+a[5]*z+a[6];
1053 Result[Z] = a[2]*x+a[5]*y+2*a[7]*z+a[8];
1054
1055 break;
1056
1057 case 3:
1058
1059 x2 = x * x; y2 = y * y; z2 = z * z;
1060
1061 /* Cubic partial derivatives */
1062
1063 Result[X] = 3*a[0]*x2 + 2*x*(a[1]*y + a[2]*z + a[3]) + a[4]*y2 +
1064 y*(a[5]*z + a[6]) + a[7]*z2 + a[8]*z + a[9];
1065 Result[Y] = a[1]*x2 + x*(2*a[4]*y + a[5]*z + a[6]) + 3*a[10]*y2 +
1066 2*y*(a[11]*z + a[12]) + a[13]*z2 + a[14]*z + a[15];
1067 Result[Z] = a[2]*x2 + x*(a[5]*y + 2*a[7]*z + a[8]) + a[11]*y2 +
1068 y*(2*a[13]*z + a[14]) + 3*a[16]*z2 + 2*a[17]*z + a[18];
1069
1070 break;
1071
1072 case 4:
1073
1074 /* Quartic partial derivatives */
1075
1076 x2 = x * x; y2 = y * y; z2 = z * z;
1077 x3 = x * x2; y3 = y * y2; z3 = z * z2;
1078
1079 Result[X] = 4*a[ 0]*x3+3*x2*(a[ 1]*y+a[ 2]*z+a[ 3])+
1080 2*x*(a[ 4]*y2+y*(a[ 5]*z+a[ 6])+a[ 7]*z2+a[ 8]*z+a[ 9])+
1081 a[10]*y3+y2*(a[11]*z+a[12])+y*(a[13]*z2+a[14]*z+a[15])+
1082 a[16]*z3+a[17]*z2+a[18]*z+a[19];
1083 Result[Y] = a[ 1]*x3+x2*(2*a[ 4]*y+a[ 5]*z+a[ 6])+
1084 x*(3*a[10]*y2+2*y*(a[11]*z+a[12])+a[13]*z2+a[14]*z+a[15])+
1085 4*a[20]*y3+3*y2*(a[21]*z+a[22])+2*y*(a[23]*z2+a[24]*z+a[25])+
1086 a[26]*z3+a[27]*z2+a[28]*z+a[29];
1087 Result[Z] = a[ 2]*x3+x2*(a[ 5]*y+2*a[ 7]*z+a[ 8])+
1088 x*(a[11]*y2+y*(2*a[13]*z+a[14])+3*a[16]*z2+2*a[17]*z+a[18])+
1089 a[21]*y3+y2*(2*a[23]*z+a[24])+y*(3*a[26]*z2+2*a[27]*z+a[28])+
1090 4*a[30]*z3+3*a[31]*z2+2*a[32]*z+a[33];
1091 }
1092 }
1093
1094
1095
1096 /*****************************************************************************
1097 *
1098 * FUNCTION
1099 *
1100 * Inside_Poly
1101 *
1102 * INPUT
1103 *
1104 * OUTPUT
1105 *
1106 * RETURNS
1107 *
1108 * AUTHOR
1109 *
1110 * Alexander Enzmann
1111 *
1112 * DESCRIPTION
1113 *
1114 * -
1115 *
1116 * CHANGES
1117 *
1118 * -
1119 *
1120 ******************************************************************************/
1121
Inside(const Vector3d & IPoint,TraceThreadData * Thread) const1122 bool Poly::Inside(const Vector3d& IPoint, TraceThreadData *Thread) const
1123 {
1124 Vector3d New_Point;
1125 DBL Result;
1126
1127 /* Transform the point into polynomial's space */
1128
1129 MInvTransPoint(New_Point, IPoint, Trans);
1130
1131 Result = inside(New_Point, Order, Coeffs);
1132
1133 if (Result < INSIDE_TOLERANCE)
1134 {
1135 return ((int)(!Test_Flag(this, INVERTED_FLAG)));
1136 }
1137 else
1138 {
1139 return ((int)Test_Flag(this, INVERTED_FLAG));
1140 }
1141 }
1142
1143
1144
1145 /*****************************************************************************
1146 *
1147 * FUNCTION
1148 *
1149 * Poly_Normal
1150 *
1151 * INPUT
1152 *
1153 * OUTPUT
1154 *
1155 * RETURNS
1156 *
1157 * AUTHOR
1158 *
1159 * Alexander Enzmann
1160 *
1161 * DESCRIPTION
1162 *
1163 * Normal to a polynomial.
1164 *
1165 * CHANGES
1166 *
1167 * -
1168 *
1169 ******************************************************************************/
1170
Normal(Vector3d & Result,Intersection * Inter,TraceThreadData * Thread) const1171 void Poly::Normal(Vector3d& Result, Intersection *Inter, TraceThreadData *Thread) const
1172 {
1173 DBL val;
1174 Vector3d New_Point;
1175
1176 /* Transform the point into the polynomials space. */
1177
1178 MInvTransPoint(New_Point, Inter->IPoint, Trans);
1179
1180 if (Order > 4)
1181 {
1182 normal0(Result, Order, Coeffs, New_Point);
1183 }
1184 else
1185 {
1186 normal1(Result, Order, Coeffs, New_Point);
1187 }
1188
1189 /* Transform back to world space. */
1190
1191 MTransNormal(Result, Result, Trans);
1192
1193 /* Normalize (accounting for the possibility of a 0 length normal). */
1194
1195 val = Result.lengthSqr();
1196
1197 if (val > 0.0)
1198 {
1199 val = 1.0 / sqrt(val);
1200
1201 Result *= val;
1202 }
1203 else
1204 {
1205 Result = Vector3d(1.0, 0.0, 0.0);
1206 }
1207 }
1208
1209
1210
1211 /*****************************************************************************
1212 *
1213 * FUNCTION
1214 *
1215 * Translate_Poly
1216 *
1217 * INPUT
1218 *
1219 * OUTPUT
1220 *
1221 * RETURNS
1222 *
1223 * AUTHOR
1224 *
1225 * Alexander Enzmann
1226 *
1227 * DESCRIPTION
1228 *
1229 * -
1230 *
1231 * CHANGES
1232 *
1233 * -
1234 *
1235 ******************************************************************************/
1236
Translate(const Vector3d &,const TRANSFORM * tr)1237 void Poly::Translate(const Vector3d&, const TRANSFORM *tr)
1238 {
1239 Transform(tr);
1240 }
1241
1242
1243
1244 /*****************************************************************************
1245 *
1246 * FUNCTION
1247 *
1248 * Rotate_Poly
1249 *
1250 * INPUT
1251 *
1252 * OUTPUT
1253 *
1254 * RETURNS
1255 *
1256 * AUTHOR
1257 *
1258 * Alexander Enzmann
1259 *
1260 * DESCRIPTION
1261 *
1262 * -
1263 *
1264 * CHANGES
1265 *
1266 * -
1267 *
1268 ******************************************************************************/
1269
Rotate(const Vector3d &,const TRANSFORM * tr)1270 void Poly::Rotate(const Vector3d&, const TRANSFORM *tr)
1271 {
1272 Transform(tr);
1273 }
1274
1275
1276
1277 /*****************************************************************************
1278 *
1279 * FUNCTION
1280 *
1281 * Scale_Poly
1282 *
1283 * INPUT
1284 *
1285 * OUTPUT
1286 *
1287 * RETURNS
1288 *
1289 * AUTHOR
1290 *
1291 * Alexander Enzmann
1292 *
1293 * DESCRIPTION
1294 *
1295 * -
1296 *
1297 * CHANGES
1298 *
1299 * -
1300 *
1301 ******************************************************************************/
1302
Scale(const Vector3d &,const TRANSFORM * tr)1303 void Poly::Scale(const Vector3d&, const TRANSFORM *tr)
1304 {
1305 Transform(tr);
1306 }
1307
1308
1309
1310 /*****************************************************************************
1311 *
1312 * FUNCTION
1313 *
1314 * Transform_Poly
1315 *
1316 * INPUT
1317 *
1318 * OUTPUT
1319 *
1320 * RETURNS
1321 *
1322 * AUTHOR
1323 *
1324 * Alexander Enzmann
1325 *
1326 * DESCRIPTION
1327 *
1328 * -
1329 *
1330 * CHANGES
1331 *
1332 * -
1333 *
1334 ******************************************************************************/
1335
Transform(const TRANSFORM * tr)1336 void Poly::Transform(const TRANSFORM *tr)
1337 {
1338 Compose_Transforms(Trans, tr);
1339
1340 Compute_BBox();
1341 }
1342
1343
1344
1345 /*****************************************************************************
1346 *
1347 * FUNCTION
1348 *
1349 * Create_Poly
1350 *
1351 * INPUT
1352 *
1353 * OUTPUT
1354 *
1355 * RETURNS
1356 *
1357 * AUTHOR
1358 *
1359 * Alexander Enzmann
1360 *
1361 * DESCRIPTION
1362 *
1363 * -
1364 *
1365 * CHANGES
1366 *
1367 * -
1368 *
1369 ******************************************************************************/
1370
Poly(int order)1371 Poly::Poly(int order) : ObjectBase(POLY_OBJECT)
1372 {
1373 Order = order;
1374
1375 Trans = Create_Transform();
1376
1377 Coeffs = reinterpret_cast<DBL *>(POV_MALLOC(term_counts(Order) * sizeof(DBL), "coefficients for POLY"));
1378
1379 for (int i = 0; i < term_counts(Order); i++)
1380 Coeffs[i] = 0.0;
1381 }
1382
1383
1384
1385 /*****************************************************************************
1386 *
1387 * FUNCTION
1388 *
1389 * Copy_Poly
1390 *
1391 * INPUT
1392 *
1393 * OUTPUT
1394 *
1395 * RETURNS
1396 *
1397 * AUTHOR
1398 *
1399 * Alexander Enzmann
1400 *
1401 * DESCRIPTION
1402 *
1403 * Make a copy of a polynomial object.
1404 *
1405 * CHANGES
1406 *
1407 * -
1408 *
1409 ******************************************************************************/
1410
Copy()1411 ObjectPtr Poly::Copy()
1412 {
1413 Poly *New = new Poly(Order);
1414 DBL *tmpCoeffs = New->Coeffs;
1415 int i;
1416
1417 Destroy_Transform(New->Trans);
1418 *New = *this;
1419 New->Coeffs = tmpCoeffs;
1420 New->Trans = Copy_Transform(Trans);
1421
1422 for(i = 0; i < term_counts(New->Order); i++)
1423 New->Coeffs[i] = Coeffs[i];
1424
1425 return New;
1426 }
1427
1428
1429
1430 /*****************************************************************************
1431 *
1432 * FUNCTION
1433 *
1434 * Destroy_Poly
1435 *
1436 * INPUT
1437 *
1438 * OUTPUT
1439 *
1440 * RETURNS
1441 *
1442 * AUTHOR
1443 *
1444 * Alexander Enzmann
1445 *
1446 * DESCRIPTION
1447 *
1448 * -
1449 *
1450 * CHANGES
1451 *
1452 * -
1453 *
1454 ******************************************************************************/
1455
~Poly()1456 Poly::~Poly()
1457 {
1458 POV_FREE(Coeffs);
1459 }
1460
1461
1462
1463 /*****************************************************************************
1464 *
1465 * FUNCTION
1466 *
1467 * Compute_Poly_BBox
1468 *
1469 * INPUT
1470 *
1471 * Poly - Poly
1472 *
1473 * OUTPUT
1474 *
1475 * Poly
1476 *
1477 * RETURNS
1478 *
1479 * AUTHOR
1480 *
1481 * Dieter Bayer
1482 *
1483 * DESCRIPTION
1484 *
1485 * Calculate the bounding box of a poly.
1486 *
1487 * CHANGES
1488 *
1489 * Aug 1994 : Creation.
1490 *
1491 ******************************************************************************/
1492
Compute_BBox()1493 void Poly::Compute_BBox()
1494 {
1495 Make_BBox(BBox, -BOUND_HUGE/2, -BOUND_HUGE/2, -BOUND_HUGE/2, BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
1496
1497 if(!Clip.empty())
1498 BBox = Clip[0]->BBox; // FIXME - does not seem to support more than one bounding object? [trf]
1499 }
1500
Intersect_BBox(BBoxDirection,const BBoxVector3d &,const BBoxVector3d &,BBoxScalar) const1501 bool Poly::Intersect_BBox(BBoxDirection, const BBoxVector3d&, const BBoxVector3d&, BBoxScalar) const
1502 {
1503 return true;
1504 }
1505
1506 }
1507