1 /*******************************************************************************
2 * quatern.cpp
3 *
4 * This module implements Quaternion algebra julia fractals.
5 *
6 * This file was written by Pascal Massimino.
7 * Revised and updated for POV-Ray 3.x by Tim Wegner
8 *
9 * ---------------------------------------------------------------------------
10 * Persistence of Vision Ray Tracer ('POV-Ray') version 3.7.
11 * Copyright 1991-2013 Persistence of Vision Raytracer Pty. Ltd.
12 *
13 * POV-Ray is free software: you can redistribute it and/or modify
14 * it under the terms of the GNU Affero General Public License as
15 * published by the Free Software Foundation, either version 3 of the
16 * License, or (at your option) any later version.
17 *
18 * POV-Ray is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU Affero General Public License for more details.
22 *
23 * You should have received a copy of the GNU Affero General Public License
24 * along with this program. If not, see <http://www.gnu.org/licenses/>.
25 * ---------------------------------------------------------------------------
26 * POV-Ray is based on the popular DKB raytracer version 2.12.
27 * DKBTrace was originally written by David K. Buck.
28 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
29 * ---------------------------------------------------------------------------
30 * $File: //depot/public/povray/3.x/source/backend/math/quatern.cpp $
31 * $Revision: #1 $
32 * $Change: 6069 $
33 * $DateTime: 2013/11/06 11:59:40 $
34 * $Author: chrisc $
35 *******************************************************************************/
36
37 // frame.h must always be the first POV file included (pulls in platform config)
38 #include "backend/frame.h"
39 #include "povray.h"
40 #include "backend/math/vector.h"
41 #include "backend/shape/fractal.h"
42 #include "backend/math/quatern.h"
43 #include "backend/shape/spheres.h"
44
45 // this must be the last file included
46 #include "base/povdebug.h"
47
48 namespace pov
49 {
50
51 /*****************************************************************************
52 * Local preprocessor defines
53 ******************************************************************************/
54
55 #define Deriv_z2(n1,n2,n3,n4) \
56 { \
57 tmp = (n1)*x - (n2)*y - (n3)*z - (n4)*w; \
58 (n2) = (n1)*y + x*(n2); \
59 (n3) = (n1)*z + x*(n3); \
60 (n4) = (n1)*w + x*(n4); \
61 (n1) = tmp; \
62 }
63
64 #define Deriv_z3(n1,n2,n3,n4) \
65 { \
66 dtmp = 2.0*((n2)*y + (n3)*z + (n4)*w); \
67 dtmp2 = 6.0*x*(n1) - dtmp; \
68 (n1) = ( (n1)*x3 - x*dtmp )*3.0; \
69 (n2) = (n2)*x4 + y*dtmp2; \
70 (n3) = (n3)*x4 + z*dtmp2; \
71 (n4) = (n4)*x4 + w*dtmp2; \
72 }
73
74
75 /*****************************************************************************
76 * Local typedefs
77 ******************************************************************************/
78
79
80
81 /*****************************************************************************
82 * Static functions
83 ******************************************************************************/
84
85
86
87 /*****************************************************************************
88 * Local variables
89 ******************************************************************************/
90
91
92
93 /*****************************************************************************
94 *
95 * FUNCTION
96 *
97 * INPUT
98 *
99 * OUTPUT
100 *
101 * RETURNS
102 *
103 * AUTHOR
104 *
105 * Pascal Massimino
106 *
107 * DESCRIPTION
108 *
109 * -
110 *
111 * CHANGES
112 *
113 * Dec 1994 : Creation.
114 *
115 ******************************************************************************/
116
Iteration_z3(const VECTOR point,const Fractal * Julia,DBL ** IterStack)117 int Iteration_z3(const VECTOR point, const Fractal *Julia, DBL **IterStack)
118 {
119 int i;
120 DBL x, y, z, w;
121 DBL d, x2, tmp;
122 DBL Exit_Value;
123
124 IterStack[X][0] = x = point[X];
125 IterStack[Y][0] = y = point[Y];
126 IterStack[Z][0] = z = point[Z];
127 IterStack[W][0] = w = (Julia->SliceDist
128 - Julia->Slice[X]*x
129 - Julia->Slice[Y]*y
130 - Julia->Slice[Z]*z)/Julia->Slice[T];
131
132 Exit_Value = Julia->Exit_Value;
133
134 for (i = 1; i <= Julia->Num_Iterations; ++i)
135 {
136 d = y * y + z * z + w * w;
137
138 x2 = x * x;
139
140 if ((d + x2) > Exit_Value)
141 {
142 return (false);
143 }
144
145 tmp = 3.0 * x2 - d;
146
147 IterStack[X][i] = x = x * (x2 - 3.0 * d) + Julia->Julia_Parm[X];
148 IterStack[Y][i] = y = y * tmp + Julia->Julia_Parm[Y];
149 IterStack[Z][i] = z = z * tmp + Julia->Julia_Parm[Z];
150 IterStack[W][i] = w = w * tmp + Julia->Julia_Parm[T];
151 }
152
153 return (true);
154 }
155
156
157
158 /*****************************************************************************
159 *
160 * FUNCTION
161 *
162 * INPUT
163 *
164 * OUTPUT
165 *
166 * RETURNS
167 *
168 * AUTHOR
169 *
170 * Pascal Massimino
171 *
172 * DESCRIPTION
173 *
174 * -
175 *
176 * CHANGES
177 *
178 * Dec 1994 : Creation.
179 *
180 ******************************************************************************/
181
Iteration_Julia(const VECTOR point,const Fractal * Julia,DBL ** IterStack)182 int Iteration_Julia(const VECTOR point, const Fractal *Julia, DBL **IterStack)
183 {
184 int i;
185 DBL x, y, z, w;
186 DBL d, x2;
187 DBL Exit_Value;
188
189 IterStack[X][0] = x = point[X];
190 IterStack[Y][0] = y = point[Y];
191 IterStack[Z][0] = z = point[Z];
192 IterStack[W][0] = w = (Julia->SliceDist
193 - Julia->Slice[X]*x
194 - Julia->Slice[Y]*y
195 - Julia->Slice[Z]*z)/Julia->Slice[T];
196
197 Exit_Value = Julia->Exit_Value;
198
199 for (i = 1; i <= Julia->Num_Iterations; ++i)
200 {
201 d = y * y + z * z + w * w;
202
203 x2 = x * x;
204
205 if ((d + x2) > Exit_Value)
206 {
207 return (false);
208 }
209
210 x *= 2.0;
211
212 IterStack[Y][i] = y = x * y + Julia->Julia_Parm[Y];
213 IterStack[Z][i] = z = x * z + Julia->Julia_Parm[Z];
214 IterStack[W][i] = w = x * w + Julia->Julia_Parm[T];
215 IterStack[X][i] = x = x2 - d + Julia->Julia_Parm[X];;
216
217 }
218
219 return (true);
220 }
221
222
223
224 /*****************************************************************************
225 *
226 * FUNCTION
227 *
228 * INPUT
229 *
230 * OUTPUT
231 *
232 * RETURNS
233 *
234 * AUTHOR
235 *
236 * Pascal Massimino
237 *
238 * DESCRIPTION
239 *
240 * D_Iteration puts in *Dist a lower bound for the distance from *point to the
241 * set
242 *
243 * CHANGES
244 *
245 * Dec 1994 : Creation.
246 *
247 ******************************************************************************/
248
249 /*----------- Distance estimator + iterations ------------*/
250
D_Iteration_z3(const VECTOR point,const Fractal * Julia,const VECTOR &,DBL * Dist,DBL ** IterStack)251 int D_Iteration_z3(const VECTOR point, const Fractal *Julia, const VECTOR &, DBL *Dist, DBL **IterStack)
252 {
253 int i, j;
254 DBL Norm, d;
255 DBL xx, yy, zz;
256 DBL x, y, z, w;
257 DBL tmp, x2;
258 DBL Exit_Value;
259 DBL Pow;
260
261 x = IterStack[X][0] = point[X];
262 y = IterStack[Y][0] = point[Y];
263 z = IterStack[Z][0] = point[Z];
264 w = IterStack[W][0] = (Julia->SliceDist
265 - Julia->Slice[X]*x
266 - Julia->Slice[Y]*y
267 - Julia->Slice[Z]*z)/Julia->Slice[T];
268
269 Exit_Value = Julia->Exit_Value;
270
271 for (i = 1; i <= Julia->Num_Iterations; i++)
272 {
273 d = y * y + z * z + w * w;
274
275 x2 = x * x;
276
277 if ((Norm = d + x2) > Exit_Value)
278 {
279 /* Distance estimator */
280
281 x = IterStack[X][0];
282 y = IterStack[Y][0];
283 z = IterStack[Z][0];
284 w = IterStack[W][0];
285
286 Pow = 1.0 / 3.0;
287
288 for (j = 1; j < i; ++j)
289 {
290 xx = x * IterStack[X][j] - y * IterStack[Y][j] - z * IterStack[Z][j] - w * IterStack[W][j];
291 yy = x * IterStack[Y][j] + y * IterStack[X][j] - z * IterStack[W][j] + w * IterStack[Z][j];
292 zz = x * IterStack[Z][j] + y * IterStack[W][j] + z * IterStack[X][j] - w * IterStack[Y][j];
293 w = x * IterStack[W][j] - y * IterStack[Z][j] + z * IterStack[Y][j] + w * IterStack[X][j];
294
295 x = xx;
296 y = yy;
297 z = zz;
298
299 Pow /= 3.0;
300 }
301
302 *Dist = Pow * sqrt(Norm / (x * x + y * y + z * z + w * w)) * log(Norm);
303
304 return (false);
305 }
306
307 tmp = 3.0 * x2 - d;
308
309 IterStack[X][i] = x = x * (x2 - 3.0 * d) + Julia->Julia_Parm[X];
310 IterStack[Y][i] = y = y * tmp + Julia->Julia_Parm[Y];
311 IterStack[Z][i] = z = z * tmp + Julia->Julia_Parm[Z];
312 IterStack[W][i] = w = w * tmp + Julia->Julia_Parm[T];
313 }
314
315 *Dist = Julia->Precision;
316
317 return (true);
318 }
319
320
321
322 /*****************************************************************************
323 *
324 * FUNCTION
325 *
326 * INPUT
327 *
328 * OUTPUT
329 *
330 * RETURNS
331 *
332 * AUTHOR
333 *
334 * Pascal Massimino
335 *
336 * DESCRIPTION
337 *
338 * -
339 *
340 * CHANGES
341 *
342 * Dec 1994 : Creation.
343 *
344 ******************************************************************************/
345
D_Iteration_Julia(const VECTOR point,const Fractal * Julia,const VECTOR &,DBL * Dist,DBL ** IterStack)346 int D_Iteration_Julia(const VECTOR point, const Fractal *Julia, const VECTOR &, DBL *Dist, DBL **IterStack)
347 {
348 int i, j;
349 DBL Norm, d;
350 DBL Exit_Value;
351 DBL x, y, z, w;
352 DBL xx, yy, zz, x2;
353 DBL Pow;
354
355 x = IterStack[X][0] = point[X];
356 y = IterStack[Y][0] = point[Y];
357 z = IterStack[Z][0] = point[Z];
358 w = IterStack[W][0] = (Julia->SliceDist
359 - Julia->Slice[X]*x
360 - Julia->Slice[Y]*y
361 - Julia->Slice[Z]*z)/Julia->Slice[T];
362
363 Exit_Value = Julia->Exit_Value;
364
365 for (i = 1; i <= Julia->Num_Iterations; i++)
366 {
367 d = y * y + z * z + w * w;
368
369 x2 = x * x;
370
371 if ((Norm = d + x2) > Exit_Value)
372 {
373 /* Distance estimator */
374
375 x = IterStack[X][0];
376 y = IterStack[Y][0];
377 z = IterStack[Z][0];
378 w = IterStack[W][0];
379
380 Pow = 1.0 / 2.0;
381
382 for (j = 1; j < i; ++j)
383 {
384 xx = x * IterStack[X][j] - y * IterStack[Y][j] - z * IterStack[Z][j] - w * IterStack[W][j];
385 yy = x * IterStack[Y][j] + y * IterStack[X][j] + w * IterStack[Z][j] - z * IterStack[W][j];
386 zz = x * IterStack[Z][j] + z * IterStack[X][j] + y * IterStack[W][j] - w * IterStack[Y][j];
387 w = x * IterStack[W][j] + w * IterStack[X][j] + z * IterStack[Y][j] - y * IterStack[Z][j];
388
389 x = xx;
390 y = yy;
391 z = zz;
392
393 Pow /= 2.0;
394 }
395
396 *Dist = Pow * sqrt(Norm / (x * x + y * y + z * z + w * w)) * log(Norm);
397
398 return (false);
399 }
400
401 x *= 2.0;
402
403 IterStack[Y][i] = y = x * y + Julia->Julia_Parm[Y];
404 IterStack[Z][i] = z = x * z + Julia->Julia_Parm[Z];
405 IterStack[W][i] = w = x * w + Julia->Julia_Parm[T];
406 IterStack[X][i] = x = x2 - d + Julia->Julia_Parm[X];
407
408 }
409
410 *Dist = Julia->Precision;
411
412 return (true);
413 }
414
415 /*****************************************************************************
416 *
417 * FUNCTION
418 *
419 * INPUT
420 *
421 * OUTPUT
422 *
423 * RETURNS
424 *
425 * AUTHOR
426 *
427 * Pascal Massimino
428 *
429 * DESCRIPTION
430 *
431 * Provided the iterations sequence has been built, perform the computation of
432 * the Normal
433 *
434 * CHANGES
435 *
436 * Dec 1994 : Creation.
437 *
438 ******************************************************************************/
439
Normal_Calc_z3(VECTOR Result,int N_Max,const Fractal *,DBL ** IterStack)440 void Normal_Calc_z3(VECTOR Result, int N_Max, const Fractal *, DBL **IterStack)
441 {
442 DBL
443 n11 = 1.0, n12 = 0.0, n13 = 0.0, n14 = 0.0,
444 n21 = 0.0, n22 = 1.0, n23 = 0.0, n24 = 0.0,
445 n31 = 0.0, n32 = 0.0, n33 = 1.0, n34 = 0.0;
446
447 DBL x, y, z, w;
448 int i;
449 DBL tmp, dtmp, dtmp2, x2, x3, x4;
450
451 x = IterStack[X][0];
452 y = IterStack[Y][0];
453 z = IterStack[Z][0];
454 w = IterStack[W][0];
455
456 for (i = 1; i <= N_Max; i++)
457 {
458 tmp = y * y + z * z + w * w;
459
460 x2 = x * x;
461 x3 = x2 - tmp;
462 x4 = 3.0 * x2 - tmp;
463
464 Deriv_z3(n11, n12, n13, n14);
465 Deriv_z3(n21, n22, n23, n24);
466 Deriv_z3(n31, n32, n33, n34);
467
468 x = IterStack[X][i];
469 y = IterStack[Y][i];
470 z = IterStack[Z][i];
471 w = IterStack[W][i];
472 }
473
474 Result[X] = n11 * x + n12 * y + n13 * z + n14 * w;
475 Result[Y] = n21 * x + n22 * y + n23 * z + n24 * w;
476 Result[Z] = n31 * x + n32 * y + n33 * z + n34 * w;
477 }
478
479
480
481 /*****************************************************************************
482 *
483 * FUNCTION
484 *
485 * INPUT
486 *
487 * OUTPUT
488 *
489 * RETURNS
490 *
491 * AUTHOR
492 *
493 * Pascal Massimino
494 *
495 * DESCRIPTION
496 *
497 * -
498 *
499 * CHANGES
500 *
501 * Dec 1994 : Creation.
502 *
503 ******************************************************************************/
504
Normal_Calc_Julia(VECTOR Result,int N_Max,const Fractal *,DBL ** IterStack)505 void Normal_Calc_Julia(VECTOR Result, int N_Max, const Fractal *, DBL **IterStack)
506 {
507 DBL
508 n11 = 1.0, n12 = 0.0, n13 = 0.0, n14 = 0.0,
509 n21 = 0.0, n22 = 1.0, n23 = 0.0, n24 = 0.0,
510 n31 = 0.0, n32 = 0.0, n33 = 1.0, n34 = 0.0;
511 DBL tmp;
512 DBL x, y, z, w;
513 int i;
514
515 x = IterStack[X][0];
516 y = IterStack[Y][0];
517 z = IterStack[Z][0];
518 w = IterStack[W][0];
519
520 for (i = 1; i <= N_Max; i++)
521 {
522 Deriv_z2(n11, n12, n13, n14);
523 Deriv_z2(n21, n22, n23, n24);
524 Deriv_z2(n31, n32, n33, n34);
525
526 x = IterStack[X][i];
527 y = IterStack[Y][i];
528 z = IterStack[Z][i];
529 w = IterStack[W][i];
530 }
531
532 Result[X] = n11 * x + n12 * y + n13 * z + n14 * w;
533 Result[Y] = n21 * x + n22 * y + n23 * z + n24 * w;
534 Result[Z] = n31 * x + n32 * y + n33 * z + n34 * w;
535 }
536
537
538
539 /*****************************************************************************
540 *
541 * FUNCTION
542 *
543 * INPUT
544 *
545 * OUTPUT
546 *
547 * RETURNS
548 *
549 * AUTHOR
550 *
551 * Pascal Massimino
552 *
553 * DESCRIPTION
554 *
555 * -
556 *
557 * CHANGES
558 *
559 * Dec 1994 : Creation.
560 *
561 ******************************************************************************/
562
F_Bound_Julia(const Ray & ray,const Fractal * fractal,DBL * Depth_Min,DBL * Depth_Max)563 int F_Bound_Julia(const Ray& ray, const Fractal *fractal, DBL *Depth_Min, DBL *Depth_Max)
564 {
565 return (Sphere::Intersect(ray, fractal->Center, fractal->Radius_Squared, Depth_Min, Depth_Max));
566 }
567
568 }
569