1 /*****************************************************************************
4 * This module implements the fractal sets primitive.
6 * This file was written by Pascal Massimino.
8 * from Persistence of Vision(tm) Ray Tracer
9 * Copyright 1996,1999 Persistence of Vision Team
10 *---------------------------------------------------------------------------
11 * NOTICE: This source code file is provided so that users may experiment
12 * with enhancements to POV-Ray and to port the software to platforms other
13 * than those supported by the POV-Ray Team. There are strict rules under
14 * which you are permitted to use this file. The rules are in the file
15 * named POVLEGAL.DOC which should be distributed with this file.
16 * If POVLEGAL.DOC is not available or for more info please contact the POV-Ray
17 * Team Coordinator by email to team-coord@povray.org or visit us on the web at
18 * http://www.povray.org. The latest version of POV-Ray may be found at this site.
20 * This program is based on the popular DKB raytracer version 2.12.
21 * DKBTrace was originally written by David K. Buck.
22 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
24 *****************************************************************************/
40 /*****************************************************************************
41 * Local preprocessor defines
42 ******************************************************************************/
44 #ifndef Fractal_Tolerance
45 #define Fractal_Tolerance 1e-7
50 /*****************************************************************************
52 ******************************************************************************/
56 /*****************************************************************************
58 ******************************************************************************/
60 static int All_Fractal_Intersections (OBJECT
* Object
, RAY
* Ray
, ISTACK
* Depth_Stack
);
61 static int Inside_Fractal (VECTOR IPoint
, OBJECT
* Object
);
62 static void Fractal_Normal (VECTOR Result
, OBJECT
* Object
, INTERSECTION
* Intersect
);
63 static FRACTAL
*Copy_Fractal (OBJECT
* Object
);
64 static void Translate_Fractal (OBJECT
* Object
, VECTOR Vector
, TRANSFORM
*Trans
);
65 static void Rotate_Fractal (OBJECT
* Object
, VECTOR Vector
, TRANSFORM
*Trans
);
66 static void Scale_Fractal (OBJECT
* Object
, VECTOR Vector
, TRANSFORM
*Trans
);
67 static void Transform_Fractal (OBJECT
* Object
, TRANSFORM
* Trans
);
68 static void Invert_Fractal (OBJECT
* Object
);
69 static void Destroy_Fractal (OBJECT
* Object
);
70 static void Compute_Fractal_BBox (FRACTAL
* Fractal
);
72 /*****************************************************************************
74 ******************************************************************************/
76 static METHODS Fractal_Methods
=
78 All_Fractal_Intersections
,
79 Inside_Fractal
, Fractal_Normal
,
80 (COPY_METHOD
)Copy_Fractal
,
81 Translate_Fractal
, Rotate_Fractal
,
82 Scale_Fractal
, Transform_Fractal
, Invert_Fractal
,
86 static int Allocated_Iteration_Stack_Length
= 0;
88 DBL
*Sx
= NULL
, *Sy
= NULL
, *Sz
= NULL
, *Sw
= NULL
;
92 static COMPLEX_FUNCTION_METHOD Complex_Function_List
[] =
94 /* must match STYPE list in fractal.h */
112 /*****************************************************************************
132 * Dec 1994 : Creation.
134 ******************************************************************************/
136 static int All_Fractal_Intersections(OBJECT
*Object
, RAY
*Ray
, ISTACK
*Depth_Stack
)
138 int Intersection_Found
;
141 DBL Depth
, Depth_Max
;
142 DBL Dist
, Dist_Next
, Len
;
144 VECTOR IPoint
, Mid_Point
, Next_Point
, Real_Pt
;
145 VECTOR Real_Normal
, F_Normal
;
147 FRACTAL
*Fractal
= (FRACTAL
*) Object
;
149 Increase_Counter(stats
[Ray_Fractal_Tests
]);
151 Intersection_Found
= FALSE
;
152 Precision
= Fractal
->Precision
;
154 /* Get into Fractal's world. */
156 if (Fractal
->Trans
!= NULL
)
158 MInvTransDirection(Direction
, Ray
->Direction
, Fractal
->Trans
);
159 VDot(Len
, Direction
, Direction
);
168 Len
= 1.0 / sqrt(Len
);
169 VScaleEq(Direction
, Len
);
172 Assign_Vector(New_Ray
.Direction
, Direction
);
173 MInvTransPoint(New_Ray
.Initial
, Ray
->Initial
, Fractal
->Trans
);
177 Assign_Vector(Direction
, Ray
->Direction
);
184 if (!F_Bound(&New_Ray
, Fractal
, &Depth
, &Depth_Max
))
189 if (Depth_Max
< Fractal_Tolerance
)
194 if (Depth
< Fractal_Tolerance
)
196 Depth
= Fractal_Tolerance
;
199 /* Jump to starting point */
201 VScale(Next_Point
, Direction
, Depth
);
202 VAddEq(Next_Point
, New_Ray
.Initial
);
204 CURRENT
= D_Iteration(Next_Point
, Fractal
, &Dist
);
206 /* Light ray starting inside ? */
210 VAddScaledEq(Next_Point
, 2.0 * Fractal_Tolerance
, Direction
);
212 Depth
+= 2.0 * Fractal_Tolerance
;
214 if (Depth
> Depth_Max
)
219 CURRENT
= D_Iteration(Next_Point
, Fractal
, &Dist
);
224 while (Depth
< Depth_Max
)
227 * Get close to the root: Advance with Next_Point, keeping track of last
228 * position in IPoint...
233 if (Dist
< Precision
)
240 if (Depth
> Depth_Max
)
242 if (Intersection_Found
)
244 Increase_Counter(stats
[Ray_Fractal_Tests_Succeeded
]);
247 return (Intersection_Found
);
250 Assign_Vector(IPoint
, Next_Point
);
251 VAddScaledEq(Next_Point
, Dist
, Direction
);
253 NEXT
= D_Iteration(Next_Point
, Fractal
, &Dist_Next
);
257 /* Set surface was crossed... */
264 Dist
= Dist_Next
; /* not reached */
268 /* then, polish the root via bisection method... */
270 while (Dist
> Fractal_Tolerance
)
273 VAddScaled(Mid_Point
, IPoint
, Dist
, Direction
);
275 Last
= Iteration(Mid_Point
, Fractal
);
279 Assign_Vector(IPoint
, Mid_Point
);
283 if (Depth
> Depth_Max
)
285 if (Intersection_Found
)
287 Increase_Counter(stats
[Ray_Fractal_Tests_Succeeded
]);
290 return (Intersection_Found
);
295 if (CURRENT
== FALSE
) /* Mid_Point isn't inside the set */
297 VAddScaledEq(IPoint
, Dist
, Direction
);
301 Iteration(IPoint
, Fractal
);
307 Iteration(IPoint
, Fractal
);
311 if (Fractal
->Trans
!= NULL
)
313 MTransPoint(Real_Pt
, IPoint
, Fractal
->Trans
);
314 Normal_Calc(Fractal
, F_Normal
);
315 MTransNormal(Real_Normal
, F_Normal
, Fractal
->Trans
);
319 Assign_Vector(Real_Pt
, IPoint
);
320 Normal_Calc(Fractal
, Real_Normal
);
323 if (Point_In_Clip(Real_Pt
, Object
->Clip
))
325 VNormalize(Real_Normal
, Real_Normal
);
326 push_normal_entry(Depth
* Len
, Real_Pt
, Real_Normal
, Object
, Depth_Stack
);
327 Intersection_Found
= TRUE
;
329 /* If fractal isn't used with CSG we can exit now. */
331 if (!(Fractal
->Type
& IS_CHILD_OBJECT
))
337 /* Start over where work was left */
339 Assign_Vector(IPoint
, Next_Point
);
345 if (Intersection_Found
)
347 Increase_Counter(stats
[Ray_Fractal_Tests_Succeeded
]);
349 return (Intersection_Found
);
352 /*****************************************************************************
372 * Dec 1994 : Creation.
374 ******************************************************************************/
376 static int Inside_Fractal(VECTOR IPoint
, OBJECT
*Object
)
378 FRACTAL
*Fractal
= (FRACTAL
*) Object
;
380 static VECTOR New_Point
;
382 if (Fractal
->Trans
!= NULL
)
384 MInvTransPoint(New_Point
, IPoint
, Fractal
->Trans
);
386 Result
= Iteration(New_Point
, Fractal
);
390 Result
= Iteration(IPoint
, Fractal
);
393 if (Fractal
->Inverted
)
403 /*****************************************************************************
423 * Dec 1994 : Creation.
425 ******************************************************************************/
427 static void Fractal_Normal(VECTOR Result
, OBJECT
*Object
, INTERSECTION
*Intersect
)
429 Assign_Vector(Result
, Intersect
->INormal
);
432 /*****************************************************************************
452 * Dec 1994 : Creation.
454 ******************************************************************************/
456 static void Translate_Fractal(OBJECT
*Object
, VECTOR Vector
, TRANSFORM
*Trans
)
458 Transform_Fractal(Object
, Trans
);
461 /*****************************************************************************
481 * Dec 1994 : Creation.
483 ******************************************************************************/
485 static void Rotate_Fractal(OBJECT
*Object
, VECTOR Vector
, TRANSFORM
*Trans
)
487 Transform_Fractal(Object
, Trans
);
490 /*****************************************************************************
510 * Dec 1994 : Creation.
512 ******************************************************************************/
514 static void Scale_Fractal(OBJECT
*Object
, VECTOR Vector
, TRANSFORM
*Trans
)
516 Transform_Fractal(Object
, Trans
);
519 /*****************************************************************************
539 * Dec 1994 : Creation.
540 * Mar 1996 : Moved call to Recompute_BBox to Compute_Fractal_BBox() (TW)
542 ******************************************************************************/
544 static void Transform_Fractal(OBJECT
*Object
, TRANSFORM
*Trans
)
546 FRACTAL
*Fractal
= (FRACTAL
*) Object
;
548 if (((FRACTAL
*) Object
)->Trans
== NULL
)
550 ((FRACTAL
*) Object
)->Trans
= Create_Transform();
553 Compose_Transforms(Fractal
->Trans
, Trans
);
555 Compute_Fractal_BBox((FRACTAL
*) Object
);
558 /*****************************************************************************
578 * Dec 1994 : Creation.
580 ******************************************************************************/
582 static void Invert_Fractal(OBJECT
*Object
)
584 ((FRACTAL
*) Object
)->Inverted
^= TRUE
;
587 /*****************************************************************************
607 * Dec 1994 : Creation.
608 * Mar 1996 : Added call to recompute_BBox() to bottom (TW)
610 ******************************************************************************/
612 static void Compute_Fractal_BBox(FRACTAL
*Fractal
)
616 switch (Fractal
->Algebra
)
618 case QUATERNION_TYPE
:
620 R
= 1.0 + sqrt(Sqr(Fractal
->Julia_Parm
[X
]) + Sqr(Fractal
->Julia_Parm
[Y
]) + Sqr(Fractal
->Julia_Parm
[Z
]) + Sqr(Fractal
->Julia_Parm
[T
]));
621 R
+= Fractal_Tolerance
; /* fix bug when Julia_Parameter exactly 0 */
628 Fractal
->Exit_Value
= Sqr(R
);
632 case HYPERCOMPLEX_TYPE
:
637 Fractal
->Exit_Value
= 16.0;
642 Fractal
->Radius_Squared
= Sqr(R
);
644 Fractal
->Inverted
= FALSE
;
646 Make_BBox(Fractal
->BBox
, -R
, -R
, -R
, 2.0 * R
, 2.0 * R
, 2.0 * R
);
648 Recompute_BBox(&Fractal
->BBox
, Fractal
->Trans
);
651 /*****************************************************************************
671 * Dec 1994 : Creation.
673 ******************************************************************************/
675 FRACTAL
*Create_Fractal()
679 New
= (FRACTAL
*) POV_MALLOC(sizeof(FRACTAL
), "Fractal Set");
681 INIT_OBJECT_FIELDS(New
, BASIC_OBJECT
, &Fractal_Methods
)
685 Make_Vector(New
->Center
, 0.0, 0.0, 0.0);
687 New
->Julia_Parm
[X
] = 1.0;
688 New
->Julia_Parm
[Y
] = 0.0;
689 New
->Julia_Parm
[Z
] = 0.0;
690 New
->Julia_Parm
[T
] = 0.0;
696 New
->SliceDist
= 0.0;
698 New
->Exit_Value
= 4.0;
702 New
->Precision
= 1.0 / 20.0;
704 New
->Inverted
= FALSE
;
706 New
->Algebra
= QUATERNION_TYPE
;
708 New
->Sub_Type
= SQR_STYPE
;
714 New
->Normal_Calc_Method
= NULL
;
715 New
->Iteration_Method
= NULL
;
716 New
->D_Iteration_Method
= NULL
;
717 New
->F_Bound_Method
= NULL
;
718 New
->Complex_Function_Method
= NULL
;
720 New
->Radius_Squared
= 0.0;
721 New
->exponent
.x
= 0.0;
722 New
->exponent
.y
= 0.0;
727 /*****************************************************************************
747 * Dec 1994 : Creation.
749 ******************************************************************************/
751 static FRACTAL
*Copy_Fractal(OBJECT
*Object
)
755 New
= Create_Fractal();
757 *New
= *((FRACTAL
*) Object
);
759 New
->Trans
= Copy_Transform(((FRACTAL
*) Object
)->Trans
);
764 /*****************************************************************************
784 * Dec 1994 : Creation.
786 ******************************************************************************/
788 static void Destroy_Fractal(OBJECT
*Object
)
790 Destroy_Transform(((FRACTAL
*) Object
)->Trans
);
794 /*****************************************************************************
814 * Dec 1994 : Creation.
816 ******************************************************************************/
818 void SetUp_Fractal(FRACTAL
*Fractal
)
820 switch (Fractal
->Algebra
)
822 case QUATERNION_TYPE
:
824 switch(Fractal
->Sub_Type
)
827 Fractal
->Iteration_Method
= Iteration_z3
;
828 Fractal
->Normal_Calc_Method
= Normal_Calc_z3
;
829 Fractal
->D_Iteration_Method
= D_Iteration_z3
;
832 Fractal
->Iteration_Method
= Iteration_Julia
;
833 Fractal
->D_Iteration_Method
= D_Iteration_Julia
;
834 Fractal
->Normal_Calc_Method
= Normal_Calc_Julia
;
837 Error("illegal function: quaternion only supports sqr and cube");
839 Fractal
->F_Bound_Method
= F_Bound_Julia
;
843 case HYPERCOMPLEX_TYPE
:
845 switch (Fractal
->Sub_Type
)
847 case RECIPROCAL_STYPE
:
849 Fractal
->Iteration_Method
= Iteration_HCompl_Reciprocal
;
850 Fractal
->Normal_Calc_Method
= Normal_Calc_HCompl_Reciprocal
;
851 Fractal
->D_Iteration_Method
= D_Iteration_HCompl_Reciprocal
;
852 Fractal
->F_Bound_Method
= F_Bound_HCompl_Reciprocal
;
872 Fractal
->Iteration_Method
= Iteration_HCompl_Func
;
873 Fractal
->Normal_Calc_Method
= Normal_Calc_HCompl_Func
;
874 Fractal
->D_Iteration_Method
= D_Iteration_HCompl_Func
;
875 Fractal
->F_Bound_Method
= F_Bound_HCompl_Func
;
876 Fractal
->Complex_Function_Method
= Complex_Function_List
[Fractal
->Sub_Type
];
882 Fractal
->Iteration_Method
= Iteration_HCompl_z3
;
883 Fractal
->Normal_Calc_Method
= Normal_Calc_HCompl_z3
;
884 Fractal
->D_Iteration_Method
= D_Iteration_HCompl_z3
;
885 Fractal
->F_Bound_Method
= F_Bound_HCompl_z3
;
889 default: /* SQR_STYPE or else... */
891 Fractal
->Iteration_Method
= Iteration_HCompl
;
892 Fractal
->Normal_Calc_Method
= Normal_Calc_HCompl
;
893 Fractal
->D_Iteration_Method
= D_Iteration_HCompl
;
894 Fractal
->F_Bound_Method
= F_Bound_HCompl
;
903 Error("Algebra unknown in fractal.");
906 Allocate_Iteration_Stack(Fractal
->n
);
908 Compute_Fractal_BBox(Fractal
);
911 /*****************************************************************************
931 * Dec 1994 : Creation.
933 ******************************************************************************/
935 void Allocate_Iteration_Stack(int n
)
937 if (n
> Allocated_Iteration_Stack_Length
)
939 Sx
= (DBL
*) POV_REALLOC(Sx
, (n
+ 1) * sizeof(DBL
), "x iteration stack");
940 Sy
= (DBL
*) POV_REALLOC(Sy
, (n
+ 1) * sizeof(DBL
), "y iteration stack");
941 Sz
= (DBL
*) POV_REALLOC(Sz
, (n
+ 1) * sizeof(DBL
), "w iteration stack");
942 Sw
= (DBL
*) POV_REALLOC(Sw
, (n
+ 1) * sizeof(DBL
), "z iteration stack");
944 Allocated_Iteration_Stack_Length
= n
;
948 /*****************************************************************************
968 * Dec 1994 : Creation.
970 ******************************************************************************/
972 void Free_Iteration_Stack()
987 Allocated_Iteration_Stack_Length
= 0;