1 //******************************************************************************
2 ///
3 /// @file core/shape/quadric.cpp
4 ///
5 /// Implementation of the quadric geometric primitive.
6 ///
7 /// @copyright
8 /// @parblock
9 ///
10 /// Persistence of Vision Ray Tracer ('POV-Ray') version 3.8.
11 /// Copyright 1991-2018 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 /// ----------------------------------------------------------------------------
27 ///
28 /// POV-Ray is based on the popular DKB raytracer version 2.12.
29 /// DKBTrace was originally written by David K. Buck.
30 /// DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
31 ///
32 /// @endparblock
33 ///
34 //******************************************************************************
35
36 // Unit header file must be the first file included within POV-Ray *.cpp files (pulls in config)
37 #include "core/shape/quadric.h"
38
39 #include <algorithm>
40
41 #include "core/bounding/boundingbox.h"
42 #include "core/math/matrix.h"
43 #include "core/render/ray.h"
44 #include "core/shape/plane.h"
45 #include "core/scene/tracethreaddata.h"
46
47 // this must be the last file included
48 #include "base/povdebug.h"
49
50 namespace pov
51 {
52
53 /*****************************************************************************
54 * Local preprocessor defines
55 ******************************************************************************/
56
57 /* The following defines make typing much easier! [DB 7/94] */
58
59 #define Xd ray.Direction[X]
60 #define Yd ray.Direction[Y]
61 #define Zd ray.Direction[Z]
62
63 #define Xo ray.Origin[X]
64 #define Yo ray.Origin[Y]
65 #define Zo ray.Origin[Z]
66
67 #define QA Square_Terms[X]
68 #define QE Square_Terms[Y]
69 #define QH Square_Terms[Z]
70
71 #define QB Mixed_Terms[X]
72 #define QC Mixed_Terms[Y]
73 #define QF Mixed_Terms[Z]
74
75 #define QD Terms[X]
76 #define QG Terms[Y]
77 #define QI Terms[Z]
78
79 #define QJ Constant
80
81 const DBL DEPTH_TOLERANCE = 1.0e-6;
82
83
84 /*****************************************************************************
85 *
86 * FUNCTION
87 *
88 * All_Quadric_Intersections
89 *
90 * INPUT
91 *
92 * OUTPUT
93 *
94 * RETURNS
95 *
96 * AUTHOR
97 *
98 * POV-Ray Team
99 *
100 * DESCRIPTION
101 *
102 * -
103 *
104 * CHANGES
105 *
106 * -
107 *
108 ******************************************************************************/
109
All_Intersections(const Ray & ray,IStack & Depth_Stack,TraceThreadData * Thread)110 bool Quadric::All_Intersections(const Ray& ray, IStack& Depth_Stack, TraceThreadData *Thread)
111 {
112 DBL Depth1, Depth2;
113 Vector3d IPoint;
114 bool Intersection_Found;
115
116 Intersection_Found = false;
117
118 Thread->Stats()[Ray_Quadric_Tests]++;
119 if (Intersect(ray, &Depth1, &Depth2))
120 {
121 Thread->Stats()[Ray_Quadric_Tests_Succeeded]++;
122 if ((Depth1 > DEPTH_TOLERANCE) && (Depth1 < MAX_DISTANCE))
123 {
124 IPoint = ray.Evaluate(Depth1);
125 if (Clip.empty() || Point_In_Clip(IPoint, Clip, Thread))
126 {
127 Depth_Stack->push(Intersection(Depth1, IPoint, this));
128
129 Intersection_Found = true;
130 }
131 }
132
133 if ((Depth2 > DEPTH_TOLERANCE) && (Depth2 < MAX_DISTANCE))
134 {
135 IPoint = ray.Evaluate(Depth2);
136
137 if (Clip.empty() || Point_In_Clip(IPoint, Clip, Thread))
138 {
139 Depth_Stack->push(Intersection(Depth2, IPoint, this));
140
141 Intersection_Found = true;
142 }
143 }
144 }
145
146 return(Intersection_Found);
147 }
148
149
150
151 /*****************************************************************************
152 *
153 * FUNCTION
154 *
155 * Intersect_Quadric
156 *
157 * INPUT
158 *
159 * OUTPUT
160 *
161 * RETURNS
162 *
163 * AUTHOR
164 *
165 * POV-Ray Team
166 *
167 * DESCRIPTION
168 *
169 * -
170 *
171 * CHANGES
172 *
173 * -
174 *
175 ******************************************************************************/
176
Intersect(const BasicRay & ray,DBL * Depth1,DBL * Depth2) const177 bool Quadric::Intersect(const BasicRay& ray, DBL *Depth1, DBL *Depth2) const
178 {
179 DBL a, b, c, d;
180
181 a = Xd * (QA * Xd + QB * Yd + QC * Zd) +
182 Yd * (QE * Yd + QF * Zd) +
183 Zd * QH * Zd;
184
185 b = Xd * (QA * Xo + 0.5 * (QB * Yo + QC * Zo + QD)) +
186 Yd * (QE * Yo + 0.5 * (QB * Xo + QF * Zo + QG)) +
187 Zd * (QH * Zo + 0.5 * (QC * Xo + QF * Yo + QI));
188
189 c = Xo * (QA * Xo + QB * Yo + QC * Zo + QD) +
190 Yo * (QE * Yo + QF * Zo + QG) +
191 Zo * (QH * Zo + QI) +
192 QJ;
193
194 if (a != 0.0)
195 {
196 /* The equation is quadratic - find its roots */
197
198 d = Sqr(b) - a * c;
199
200 if (d <= 0.0)
201 {
202 return(false);
203 }
204
205 d = sqrt (d);
206
207 *Depth1 = (-b + d) / (a);
208 *Depth2 = (-b - d) / (a);
209 }
210 else
211 {
212 /* There are no quadratic terms. Solve the linear equation instead. */
213
214 if (b == 0.0)
215 {
216 return(false);
217 }
218
219 *Depth1 = - 0.5 * c / b;
220 *Depth2 = MAX_DISTANCE;
221 }
222
223 return(true);
224 }
225
226
227
228 /*****************************************************************************
229 *
230 * FUNCTION
231 *
232 * Inside_Quadric
233 *
234 * INPUT
235 *
236 * OUTPUT
237 *
238 * RETURNS
239 *
240 * AUTHOR
241 *
242 * POV-Ray Team
243 *
244 * DESCRIPTION
245 *
246 * -
247 *
248 * CHANGES
249 *
250 * -
251 *
252 ******************************************************************************/
253
Inside(const Vector3d & IPoint,TraceThreadData * Thread) const254 bool Quadric::Inside(const Vector3d& IPoint, TraceThreadData *Thread) const
255 {
256 /* This is faster and shorter. [DB 7/94] */
257
258 return((IPoint[X] * (QA * IPoint[X] + QB * IPoint[Y] + QD) +
259 IPoint[Y] * (QE * IPoint[Y] + QF * IPoint[Z] + QG) +
260 IPoint[Z] * (QH * IPoint[Z] + QC * IPoint[X] + QI) + QJ) <= 0.0);
261 }
262
263
264
265 /*****************************************************************************
266 *
267 * FUNCTION
268 *
269 * Quadric_Normal
270 *
271 * INPUT
272 *
273 * OUTPUT
274 *
275 * RETURNS
276 *
277 * AUTHOR
278 *
279 * POV-Ray Team
280 *
281 * DESCRIPTION
282 *
283 * -
284 *
285 * CHANGES
286 *
287 * -
288 *
289 ******************************************************************************/
290
Normal(Vector3d & Result,Intersection * Inter,TraceThreadData * Thread) const291 void Quadric::Normal(Vector3d& Result, Intersection *Inter, TraceThreadData *Thread) const
292 {
293 DBL Len;
294
295 /* This is faster and shorter. [DB 7/94] */
296
297 Result[X] = 2.0 * QA * Inter->IPoint[X] +
298 QB * Inter->IPoint[Y] +
299 QC * Inter->IPoint[Z] +
300 QD;
301
302 Result[Y] = QB * Inter->IPoint[X] +
303 2.0 * QE * Inter->IPoint[Y] +
304 QF * Inter->IPoint[Z] +
305 QG;
306
307 Result[Z] = QC * Inter->IPoint[X] +
308 QF * Inter->IPoint[Y] +
309 2.0 * QH * Inter->IPoint[Z] +
310 QI;
311
312 Len = Result.length();
313
314 if (Len == 0.0)
315 {
316 /* The normal is not defined at this point of the surface. */
317 /* Set it to any arbitrary direction. */
318
319 Result = Vector3d(1.0, 0.0, 0.0);
320 }
321 else
322 {
323 Result /= Len;
324 }
325 }
326
327
328
329 /*****************************************************************************
330 *
331 * FUNCTION
332 *
333 * Transform_Quadric
334 *
335 * INPUT
336 *
337 * OUTPUT
338 *
339 * RETURNS
340 *
341 * AUTHOR
342 *
343 * POV-Ray Team
344 *
345 * DESCRIPTION
346 *
347 * -
348 *
349 * CHANGES
350 *
351 * -
352 *
353 ******************************************************************************/
354
Transform(const TRANSFORM * tr)355 void Quadric::Transform(const TRANSFORM *tr)
356 {
357 MATRIX Quadric_Matrix, Transform_Transposed;
358
359 Quadric_To_Matrix (Quadric_Matrix);
360
361 MTimesB (tr->inverse, Quadric_Matrix);
362 MTranspose (Transform_Transposed, tr->inverse);
363 MTimesA (Quadric_Matrix, Transform_Transposed);
364
365 Matrix_To_Quadric (Quadric_Matrix);
366
367 Recompute_BBox(&BBox, tr);
368 }
369
370
371
372 /*****************************************************************************
373 *
374 * FUNCTION
375 *
376 * Quadric_To_Matrix
377 *
378 * INPUT
379 *
380 * OUTPUT
381 *
382 * RETURNS
383 *
384 * AUTHOR
385 *
386 * POV-Ray Team
387 *
388 * DESCRIPTION
389 *
390 * -
391 *
392 * CHANGES
393 *
394 * -
395 *
396 ******************************************************************************/
397
Quadric_To_Matrix(MATRIX Matrix) const398 void Quadric::Quadric_To_Matrix(MATRIX Matrix) const
399 {
400 MZero (Matrix);
401
402 Matrix[0][0] = Square_Terms[X];
403 Matrix[1][1] = Square_Terms[Y];
404 Matrix[2][2] = Square_Terms[Z];
405 Matrix[0][1] = Mixed_Terms[X];
406 Matrix[0][2] = Mixed_Terms[Y];
407 Matrix[0][3] = Terms[X];
408 Matrix[1][2] = Mixed_Terms[Z];
409 Matrix[1][3] = Terms[Y];
410 Matrix[2][3] = Terms[Z];
411 Matrix[3][3] = Constant;
412 }
413
414
415
416 /*****************************************************************************
417 *
418 * FUNCTION
419 *
420 * Matrix_To_Quadric
421 *
422 * INPUT
423 *
424 * OUTPUT
425 *
426 * RETURNS
427 *
428 * AUTHOR
429 *
430 * POV-Ray Team
431 *
432 * DESCRIPTION
433 *
434 * -
435 *
436 * CHANGES
437 *
438 * -
439 *
440 ******************************************************************************/
441
Matrix_To_Quadric(const MATRIX Matrix)442 void Quadric::Matrix_To_Quadric(const MATRIX Matrix)
443 {
444 Square_Terms[X] = Matrix[0][0];
445 Square_Terms[Y] = Matrix[1][1];
446 Square_Terms[Z] = Matrix[2][2];
447
448 Mixed_Terms[X] = Matrix[0][1] + Matrix[1][0];
449 Mixed_Terms[Y] = Matrix[0][2] + Matrix[2][0];
450 Mixed_Terms[Z] = Matrix[1][2] + Matrix[2][1];
451
452 Terms[X] = Matrix[0][3] + Matrix[3][0];
453 Terms[Y] = Matrix[1][3] + Matrix[3][1];
454 Terms[Z] = Matrix[2][3] + Matrix[3][2];
455
456 Constant = Matrix[3][3];
457 }
458
459
460
461 /*****************************************************************************
462 *
463 * FUNCTION
464 *
465 * Translate_Quadric
466 *
467 * INPUT
468 *
469 * OUTPUT
470 *
471 * RETURNS
472 *
473 * AUTHOR
474 *
475 * POV-Ray Team
476 *
477 * DESCRIPTION
478 *
479 * -
480 *
481 * CHANGES
482 *
483 * -
484 *
485 ******************************************************************************/
486
Translate(const Vector3d &,const TRANSFORM * tr)487 void Quadric::Translate(const Vector3d&, const TRANSFORM *tr)
488 {
489 Transform(tr);
490 }
491
492
493
494 /*****************************************************************************
495 *
496 * FUNCTION
497 *
498 * Rotate_Quadric
499 *
500 * INPUT
501 *
502 * OUTPUT
503 *
504 * RETURNS
505 *
506 * AUTHOR
507 *
508 * POV-Ray Team
509 *
510 * DESCRIPTION
511 *
512 * -
513 *
514 * CHANGES
515 *
516 * -
517 *
518 ******************************************************************************/
519
Rotate(const Vector3d &,const TRANSFORM * tr)520 void Quadric::Rotate(const Vector3d&, const TRANSFORM *tr)
521 {
522 Transform(tr);
523 }
524
525
526
527 /*****************************************************************************
528 *
529 * FUNCTION
530 *
531 * Scale_Quadric
532 *
533 * INPUT
534 *
535 * OUTPUT
536 *
537 * RETURNS
538 *
539 * AUTHOR
540 *
541 * POV-Ray Team
542 *
543 * DESCRIPTION
544 *
545 * -
546 *
547 * CHANGES
548 *
549 * -
550 *
551 ******************************************************************************/
552
Scale(const Vector3d &,const TRANSFORM * tr)553 void Quadric::Scale(const Vector3d&, const TRANSFORM *tr)
554 {
555 Transform(tr);
556 }
557
558
559
560 /*****************************************************************************
561 *
562 * FUNCTION
563 *
564 * Invert_Quadric
565 *
566 * INPUT
567 *
568 * OUTPUT
569 *
570 * RETURNS
571 *
572 * AUTHOR
573 *
574 * POV-Ray Team
575 *
576 * DESCRIPTION
577 *
578 * -
579 *
580 * CHANGES
581 *
582 * -
583 *
584 ******************************************************************************/
585
Invert()586 ObjectPtr Quadric::Invert()
587 {
588 Square_Terms.invert();
589 Mixed_Terms.invert();
590 Terms.invert();
591
592 Constant *= -1.0;
593
594 return this;
595 }
596
597
598
599 /*****************************************************************************
600 *
601 * FUNCTION
602 *
603 * Create_Quadric
604 *
605 * INPUT
606 *
607 * OUTPUT
608 *
609 * RETURNS
610 *
611 * AUTHOR
612 *
613 * POV-Ray Team
614 *
615 * DESCRIPTION
616 *
617 * -
618 *
619 * CHANGES
620 *
621 * -
622 *
623 ******************************************************************************/
624
Quadric()625 Quadric::Quadric() : ObjectBase(QUADRIC_OBJECT)
626 {
627 Square_Terms = Vector3d(1.0, 1.0, 1.0);
628 Mixed_Terms = Vector3d(0.0, 0.0, 0.0);
629 Terms = Vector3d(0.0, 0.0, 0.0);
630 Constant = 1.0;
631 Automatic_Bounds = false;
632 }
633
634
635
636 /*****************************************************************************
637 *
638 * FUNCTION
639 *
640 * Copy_Quadric
641 *
642 * INPUT
643 *
644 * OUTPUT
645 *
646 * RETURNS
647 *
648 * AUTHOR
649 *
650 * POV-Ray Team
651 *
652 * DESCRIPTION
653 *
654 * -
655 *
656 * CHANGES
657 *
658 * -
659 *
660 ******************************************************************************/
661
Copy()662 ObjectPtr Quadric::Copy()
663 {
664 Quadric *New = new Quadric();
665 Destroy_Transform(New->Trans);
666 *New = *this;
667
668 return(New);
669 }
670
671
672
673 /*****************************************************************************
674 *
675 * FUNCTION
676 *
677 * Destroy_Quadric
678 *
679 * INPUT
680 *
681 * OUTPUT
682 *
683 * RETURNS
684 *
685 * AUTHOR
686 *
687 * POV-Ray Team
688 *
689 * DESCRIPTION
690 *
691 * -
692 *
693 * CHANGES
694 *
695 * -
696 *
697 ******************************************************************************/
698
~Quadric()699 Quadric::~Quadric()
700 {
701 }
702
703
704
705 /*****************************************************************************
706 *
707 * FUNCTION
708 *
709 * Compute_Quadric_BBox
710 *
711 * INPUT
712 *
713 * Quadric - Quadric object
714 *
715 * OUTPUT
716 *
717 * Quadric
718 *
719 * RETURNS
720 *
721 * AUTHOR
722 *
723 * Dieter Bayer
724 *
725 * DESCRIPTION
726 *
727 * Compute a bounding box around a quadric.
728 *
729 * This function calculates the bounding box of a quadric given in
730 * its normal form, i.e. f(x,y,z) = A*x^2 + E*y^2 + H*z^2 + J = 0.
731 *
732 * NOTE: Translated quadrics can also be bounded by this function.
733 *
734 * Supported: cones, cylinders, ellipsoids, hyperboloids.
735 *
736 * CHANGES
737 *
738 * May 1994 : Creation.
739 *
740 * Sep 1994 : Added support of hyperboloids. Improved bounding of
741 * quadrics used in CSG intersections. [DB]
742 *
743 ******************************************************************************/
744
Compute_BBox()745 void Quadric::Compute_BBox()
746 {
747 Vector3d clipMin(-1.0), clipMax(1.0);
748 Compute_BBox(clipMin, clipMax);
749 }
750
Compute_BBox(Vector3d & ClipMin,Vector3d & ClipMax)751 void Quadric::Compute_BBox(Vector3d& ClipMin, Vector3d& ClipMax)
752 {
753 DBL A, B, C, D, E, F, G, H, I, J;
754 DBL a, b, c, d;
755 DBL rx, ry, rz, rx1, rx2, ry1, ry2, rz1, rz2, x, y, z;
756 DBL New_Volume, Old_Volume;
757 Vector3d TmpMin, TmpMax, NewMin, NewMax, T1;
758 BoundingBox Old = BBox;
759
760 if(!Clip.empty())
761 {
762 /* Intersect the members bounding boxes. */
763 for(vector<ObjectPtr>::iterator it = Clip.begin(); it != Clip.end(); it++)
764 {
765 ObjectPtr p = *it;
766 if (Test_Flag(p, INVERTED_FLAG) == false)
767 {
768 if (dynamic_cast<Plane *> (p) != nullptr)
769 Compute_Plane_Min_Max(dynamic_cast<Plane *> (p), TmpMin, TmpMax);
770 else
771 Make_min_max_from_BBox(TmpMin, TmpMax, p->BBox);
772
773 ClipMin = max(ClipMin, TmpMin);
774 ClipMax = min(ClipMax, TmpMax);
775 }
776 }
777 }
778
779 /*
780 * Check for 'normal' form. If the quadric isn't in it's normal form
781 * we can't do anything (we could, but that would be to tedious!
782 * Diagonalising the quadric's 4x4 matrix, i.e. finding its eigenvalues
783 * and eigenvectors -> solving a 4th order polynom).
784 */
785
786 /* Get quadrics coefficients. */
787
788 A = Square_Terms[X];
789 E = Square_Terms[Y];
790 H = Square_Terms[Z];
791 B = Mixed_Terms[X] / 2.0;
792 C = Mixed_Terms[Y] / 2.0;
793 F = Mixed_Terms[Z] / 2.0;
794 D = Terms[X] / 2.0;
795 G = Terms[Y] / 2.0;
796 I = Terms[Z] / 2.0;
797 J = Constant;
798
799 /* Set small values to 0. */
800
801 if (fabs(A) < EPSILON) A = 0.0;
802 if (fabs(B) < EPSILON) B = 0.0;
803 if (fabs(C) < EPSILON) C = 0.0;
804 if (fabs(D) < EPSILON) D = 0.0;
805 if (fabs(E) < EPSILON) E = 0.0;
806 if (fabs(F) < EPSILON) F = 0.0;
807 if (fabs(G) < EPSILON) G = 0.0;
808 if (fabs(H) < EPSILON) H = 0.0;
809 if (fabs(I) < EPSILON) I = 0.0;
810 if (fabs(J) < EPSILON) J = 0.0;
811
812 /* Non-zero mixed terms --> return */
813
814 if ((B != 0.0) || (C != 0.0) || (F != 0.0))
815 {
816 New_Volume = (ClipMax[X] - ClipMin[X]) * (ClipMax[Y] - ClipMin[Y]) * (ClipMax[Z] - ClipMin[Z]);
817 BOUNDS_VOLUME(Old_Volume, Old);
818 if (New_Volume < Old_Volume)
819 Make_BBox_from_min_max(BBox, ClipMin, ClipMax);
820 return;
821 }
822
823 /* Non-zero linear terms --> get translation vector */
824
825 if ((D != 0.0) || (G != 0.0) || (I != 0.0))
826 {
827 if (A != 0.0)
828 {
829 T1[X] = -D / A;
830 }
831 else
832 {
833 if (D != 0.0)
834 {
835 T1[X] = -J / (2.0 * D);
836 }
837 else
838 {
839 T1[X] = 0.0;
840 }
841 }
842
843 if (E != 0.0)
844 {
845 T1[Y] = -G / E;
846 }
847 else
848 {
849 if (G != 0.0)
850 {
851 T1[Y] = -J / (2.0 * G);
852 }
853 else
854 {
855 T1[Y] = 0.0;
856 }
857 }
858
859 if (H != 0.0)
860 {
861 T1[Z] = -I / H;
862 }
863 else
864 {
865 if (I != 0.0)
866 {
867 T1[Z] = -J / (2.0 * I);
868 }
869 else
870 {
871 T1[Z] = 0.0;
872 }
873 }
874
875 /* Recalculate coefficients. */
876
877 D += A * T1[X];
878 G += E * T1[Y];
879 I += H * T1[Z];
880 J -= T1[X]*(A*T1[X] - 2.0*D) + T1[Y]*(E*T1[Y] - 2.0*G) + T1[Z]*(H*T1[Z] - 2.0*I);
881 }
882 else
883 {
884 T1 = Vector3d(0.0, 0.0, 0.0);
885 }
886
887 /* Init new bounding box. */
888
889 NewMin[X] = NewMin[Y] = NewMin[Z] = -BOUND_HUGE/2;
890 NewMax[X] = NewMax[Y] = NewMax[Z] = BOUND_HUGE/2;
891
892 /* Translate clipping box. */
893
894 ClipMin -= T1;
895 ClipMax -= T1;
896
897 // TODO FIXME - The following code disregards inside/outside information.
898 // This is a huge problem in CSG intersection bounding box computations.
899
900 /* We want A to be non-negative. */
901
902 if (A < 0.0)
903 {
904 A = -A;
905 D = -D;
906 E = -E;
907 G = -G;
908 H = -H;
909 I = -I;
910 J = -J;
911 }
912
913 /*
914 *
915 * Check for ellipsoid.
916 *
917 * x*x y*y z*z
918 * ----- + ----- + ----- - 1 = 0
919 * a*a b*b c*c
920 *
921 */
922
923 if ((A > 0.0) && (E > 0.0) && (H > 0.0) && (J < 0.0))
924 {
925 a = sqrt(-J/A);
926 b = sqrt(-J/E);
927 c = sqrt(-J/H);
928
929 NewMin[X] = -a;
930 NewMin[Y] = -b;
931 NewMin[Z] = -c;
932 NewMax[X] = a;
933 NewMax[Y] = b;
934 NewMax[Z] = c;
935 }
936
937 /*
938 *
939 * Check for cylinder (x-axis).
940 *
941 * y*y z*z
942 * ----- + ----- - 1 = 0
943 * b*b c*c
944 *
945 */
946
947 if ((A == 0.0) && (E > 0.0) && (H > 0.0) && (J < 0.0))
948 {
949 b = sqrt(-J/E);
950 c = sqrt(-J/H);
951
952 NewMin[Y] = -b;
953 NewMin[Z] = -c;
954 NewMax[Y] = b;
955 NewMax[Z] = c;
956 }
957
958 /*
959 *
960 * Check for cylinder (y-axis).
961 *
962 * x*x z*z
963 * ----- + ----- - 1 = 0
964 * a*a c*c
965 *
966 */
967
968 if ((A > 0.0) && (E == 0.0) && (H > 0.0) && (J < 0.0))
969 {
970 a = sqrt(-J/A);
971 c = sqrt(-J/H);
972
973 NewMin[X] = -a;
974 NewMin[Z] = -c;
975 NewMax[X] = a;
976 NewMax[Z] = c;
977 }
978
979 /*
980 *
981 * Check for cylinder (z-axis).
982 *
983 * x*x y*y
984 * ----- + ----- - 1 = 0
985 * a*a b*b
986 *
987 */
988
989 if ((A > 0.0) && (E > 0.0) && (H == 0.0) && (J < 0.0))
990 {
991 a = sqrt(-J/A);
992 b = sqrt(-J/E);
993
994 NewMin[X] = -a;
995 NewMin[Y] = -b;
996 NewMax[X] = a;
997 NewMax[Y] = b;
998 }
999
1000 /*
1001 *
1002 * Check for cone (x-axis).
1003 *
1004 * x*x y*y z*z
1005 * ----- - ----- - ----- = 0
1006 * a*a b*b c*c
1007 *
1008 */
1009
1010 if ((A > 0.0) && (E < 0.0) && (H < 0.0) && (J == 0.0))
1011 {
1012 a = sqrt(1.0/A);
1013 b = sqrt(-1.0/E);
1014 c = sqrt(-1.0/H);
1015
1016 /* Get radii for lower x value. */
1017
1018 x = ClipMin[X];
1019
1020 ry1 = fabs(x * b / a);
1021 rz1 = fabs(x * c / a);
1022
1023 /* Get radii for upper x value. */
1024
1025 x = ClipMax[X];
1026
1027 ry2 = fabs(x * b / a);
1028 rz2 = fabs(x * c / a);
1029
1030 ry = max(ry1, ry2);
1031 rz = max(rz1, rz2);
1032
1033 NewMin[Y] = -ry;
1034 NewMin[Z] = -rz;
1035 NewMax[Y] = ry;
1036 NewMax[Z] = rz;
1037 }
1038
1039 /*
1040 *
1041 * Check for cone (y-axis).
1042 *
1043 * x*x y*y z*z
1044 * ----- - ----- + ----- = 0
1045 * a*a b*b c*c
1046 *
1047 */
1048
1049 if ((A > 0.0) && (E < 0.0) && (H > 0.0) && (J == 0.0))
1050 {
1051 a = sqrt(1.0/A);
1052 b = sqrt(-1.0/E);
1053 c = sqrt(1.0/H);
1054
1055 /* Get radii for lower y value. */
1056
1057 y = ClipMin[Y];
1058
1059 rx1 = fabs(y * a / b);
1060 rz1 = fabs(y * c / b);
1061
1062 /* Get radii for upper y value. */
1063
1064 y = ClipMax[Y];
1065
1066 rx2 = fabs(y * a / b);
1067 rz2 = fabs(y * c / b);
1068
1069 rx = max(rx1, rx2);
1070 rz = max(rz1, rz2);
1071
1072 NewMin[X] = -rx;
1073 NewMin[Z] = -rz;
1074 NewMax[X] = rx;
1075 NewMax[Z] = rz;
1076 }
1077
1078 /*
1079 *
1080 * Check for cone (z-axis).
1081 *
1082 * x*x y*y z*z
1083 * ----- + ----- - ----- = 0
1084 * a*a b*b c*c
1085 *
1086 */
1087
1088 if ((A > 0.0) && (E > 0.0) && (H < 0.0) && (J == 0.0))
1089 {
1090 a = sqrt(1.0/A);
1091 b = sqrt(1.0/E);
1092 c = sqrt(-1.0/H);
1093
1094 /* Get radii for lower z value. */
1095
1096 z = ClipMin[Z];
1097
1098 rx1 = fabs(z * a / c);
1099 ry1 = fabs(z * b / c);
1100
1101 /* Get radii for upper z value. */
1102
1103 z = ClipMax[Z];
1104
1105 rx2 = fabs(z * a / c);
1106 ry2 = fabs(z * b / c);
1107
1108 rx = max(rx1, rx2);
1109 ry = max(ry1, ry2);
1110
1111 NewMin[X] = -rx;
1112 NewMin[Y] = -ry;
1113 NewMax[X] = rx;
1114 NewMax[Y] = ry;
1115 }
1116
1117 /*
1118 *
1119 * Check for hyperboloid (x-axis).
1120 *
1121 * x*x y*y z*z
1122 * ----- - ----- - ----- + 1 = 0
1123 * a*a b*b c*c
1124 *
1125 */
1126
1127 if ((A > 0.0) && (E < 0.0) && (H < 0.0) && (J > 0.0))
1128 {
1129 /* Get radii for lower x value. */
1130
1131 x = ClipMin[X];
1132
1133 d = 1.0 + A * Sqr(x);
1134
1135 ry1 = sqrt(-d / E);
1136 rz1 = sqrt(-d / H);
1137
1138 /* Get radii for upper x value. */
1139
1140 x = ClipMax[X];
1141
1142 d = 1.0 + A * Sqr(x);
1143
1144 ry2 = sqrt(-d / E);
1145 rz2 = sqrt(-d / H);
1146
1147 ry = max(ry1, ry2);
1148 rz = max(rz1, rz2);
1149
1150 NewMin[Y] = -ry;
1151 NewMin[Z] = -rz;
1152 NewMax[Y] = ry;
1153 NewMax[Z] = rz;
1154 }
1155
1156 /*
1157 *
1158 * Check for hyperboloid (y-axis).
1159 *
1160 * x*x y*y z*z
1161 * ----- - ----- + ----- - 1 = 0
1162 * a*a b*b c*c
1163 *
1164 */
1165
1166 if ((A > 0.0) && (E < 0.0) && (H > 0.0) && (J < 0.0))
1167 {
1168 /* Get radii for lower y value. */
1169
1170 y = ClipMin[Y];
1171
1172 d = 1.0 - E * Sqr(y);
1173
1174 rx1 = sqrt(d / A);
1175 rz1 = sqrt(d / H);
1176
1177 /* Get radii for upper y value. */
1178
1179 y = ClipMax[Y];
1180
1181 d = 1.0 - E * Sqr(y);
1182
1183 rx2 = sqrt(d / A);
1184 rz2 = sqrt(d / H);
1185
1186 rx = max(rx1, rx2);
1187 rz = max(rz1, rz2);
1188
1189 NewMin[X] = -rx;
1190 NewMin[Z] = -rz;
1191 NewMax[X] = rx;
1192 NewMax[Z] = rz;
1193 }
1194
1195 /*
1196 *
1197 * Check for hyperboloid (z-axis).
1198 *
1199 * x*x y*y z*z
1200 * ----- + ----- - ----- - 1 = 0
1201 * a*a b*b c*c
1202 *
1203 */
1204
1205 if ((A > 0.0) && (E > 0.0) && (H < 0.0) && (J < 0.0))
1206 {
1207 /* Get radii for lower z value. */
1208
1209 z = ClipMin[Z];
1210
1211 d = 1.0 - H * Sqr(z);
1212
1213 rx1 = sqrt(d / A);
1214 ry1 = sqrt(d / E);
1215
1216 /* Get radii for upper z value. */
1217
1218 z = ClipMax[Z];
1219
1220 d = 1.0 - H * Sqr(z);
1221
1222 rx2 = sqrt(d / A);
1223 ry2 = sqrt(d / E);
1224
1225 rx = max(rx1, rx2);
1226 ry = max(ry1, ry2);
1227
1228 NewMin[X] = -rx;
1229 NewMin[Y] = -ry;
1230 NewMax[X] = rx;
1231 NewMax[Y] = ry;
1232 }
1233
1234 /*
1235 *
1236 * Check for paraboloid (x-axis).
1237 *
1238 * y*y z*z
1239 * x - ----- - ----- = 0
1240 * b*b c*c
1241 *
1242 */
1243
1244 if ((A == 0.0) && (D != 0.0) && (E * H > 0.0) && (J == 0.0))
1245 {
1246 /* Get radii for lower x value. */
1247
1248 x = D * E < 0 ? max(0., ClipMin[X]) : ClipMin[X];
1249
1250 ry1 = sqrt(fabs(2.0 * D * x / E));
1251 rz1 = sqrt(fabs(2.0 * D * x / H));
1252
1253 /* Get radii for upper x value. */
1254
1255 x = D * E > 0 ? min(0., ClipMax[X]) : ClipMax[X];
1256
1257 ry2 = sqrt(fabs(2.0 * D * x / E));
1258 rz2 = sqrt(fabs(2.0 * D * x / H));
1259
1260 ry = max(ry1, ry2);
1261 rz = max(rz1, rz2);
1262
1263 if (D*E < 0) NewMin[X] = max(0., ClipMin[X]);
1264 NewMin[Y] = -ry;
1265 NewMin[Z] = -rz;
1266 if (D*E > 0) NewMax[X] = min(0., ClipMax[X]);
1267 NewMax[Y] = ry;
1268 NewMax[Z] = rz;
1269 }
1270
1271 /*
1272 *
1273 * Check for paraboloid (y-axis).
1274 *
1275 * x*x z*z
1276 * y - ----- - ----- = 0
1277 * a*a c*c
1278 *
1279 */
1280
1281 if ((E == 0.0) && (G != 0.0) && (A * H > 0.0) && (J == 0.0))
1282 {
1283 /* Get radii for lower y-value. */
1284 y = G > 0 ? ClipMin[Y] : max(0., ClipMin[Y]);
1285
1286 rx1 = sqrt(fabs(2.0 * G * y / A));
1287 rz1 = sqrt(fabs(2.0 * G * y / H));
1288
1289 /* Get radii for upper y value. */
1290
1291 y = G < 0 ? ClipMax[Y] : min(0., ClipMax[Y]);
1292
1293 rx2 = sqrt(fabs(2.0 * G * y / A));
1294 rz2 = sqrt(fabs(2.0 * G * y / H));
1295
1296 rx = max(rx1, rx2);
1297 rz = max(rz1, rz2);
1298
1299 NewMin[X] = -rx;
1300 if (G < 0) NewMin[Y] = max(0., ClipMin[Y]);
1301 NewMin[Z] = -rz;
1302 NewMax[X] = rx;
1303 if (G > 0) NewMax[Y] = min(0., ClipMax[Y]);
1304 NewMax[Z] = rz;
1305 }
1306
1307 /*
1308 *
1309 * Check for paraboloid (z-axis).
1310 *
1311 * x*x y*y
1312 * z - ----- - ----- = 0
1313 * a*a b*b
1314 *
1315 */
1316
1317 if ((H == 0.0) && (I != 0.0) && (A * E > 0.0) && (J == 0.0))
1318 {
1319 /* Get radii for lower z-value. */
1320
1321 z = I < 0 ? max(ClipMin[Z], 0.) : ClipMin[Z];
1322
1323 rx1 = sqrt(fabs(2.0 * I * z / A));
1324 ry1 = sqrt(fabs(2.0 * I * z / E));
1325
1326 /* Get radii for upper z value. */
1327
1328 z = I > 0 ? min(0., ClipMax[Z]) : ClipMax[Z];
1329
1330 rx2 = sqrt(fabs(2.0 * I * z / A));
1331 ry2 = sqrt(fabs(2.0 * I * z / E));
1332
1333 rx = max(rx1, rx2);
1334 ry = max(ry1, ry2);
1335
1336 NewMin[X] = -rx;
1337 NewMin[Y] = -ry;
1338 if (I < 0) NewMin[Z] = max(ClipMin[Z], 0.);
1339 NewMax[X] = rx;
1340 NewMax[Y] = ry;
1341 if (I > 0) NewMax[Z] = min(ClipMax[Z], 0.);
1342 }
1343
1344 /* Intersect clipping object's and quadric's bounding boxes */
1345
1346 NewMin = max(NewMin, ClipMin);
1347 NewMax = min(NewMax, ClipMax);
1348
1349 /* Use old or new bounding box? */
1350
1351 New_Volume = (NewMax[X] - NewMin[X]) * (NewMax[Y] - NewMin[Y]) * (NewMax[Z] - NewMin[Z]);
1352
1353 BOUNDS_VOLUME(Old_Volume, Old);
1354
1355 if (New_Volume < Old_Volume)
1356 {
1357 /* Add translation. */
1358 Automatic_Bounds = true;
1359
1360 NewMin += T1;
1361 NewMax += T1;
1362
1363 Make_BBox_from_min_max(BBox, NewMin, NewMax);
1364
1365 /* Beware of bounding boxes to large. */
1366
1367 if ((BBox.size[X] > CRITICAL_LENGTH) ||
1368 (BBox.size[Y] > CRITICAL_LENGTH) ||
1369 (BBox.size[Z] > CRITICAL_LENGTH))
1370 {
1371 Make_BBox(BBox, -BOUND_HUGE/2, -BOUND_HUGE/2, -BOUND_HUGE/2,
1372 BOUND_HUGE, BOUND_HUGE, BOUND_HUGE);
1373 }
1374 }
1375 }
1376
1377
1378
1379 /*****************************************************************************
1380 *
1381 * FUNCTION
1382 *
1383 * Compute_Plane_Min_Max
1384 *
1385 * INPUT
1386 *
1387 * Plane - Plane
1388 * Min, Max - Vectors containing plane's dimensions
1389 *
1390 * OUTPUT
1391 *
1392 * Min, Max
1393 *
1394 * RETURNS
1395 *
1396 * AUTHOR
1397 *
1398 * Dieter Bayer
1399 *
1400 * DESCRIPTION
1401 *
1402 * Compute min/max vectors for planes perpendicular to an axis.
1403 *
1404 * CHANGES
1405 *
1406 * May 1994 : Creation.
1407 *
1408 ******************************************************************************/
1409
Compute_Plane_Min_Max(const Plane * plane,Vector3d & Min,Vector3d & Max)1410 void Quadric::Compute_Plane_Min_Max(const Plane *plane, Vector3d& Min, Vector3d& Max)
1411 {
1412 DBL d;
1413 Vector3d P, N;
1414
1415 if (plane->Trans == nullptr)
1416 {
1417 N = plane->Normal_Vector;
1418
1419 d = -plane->Distance;
1420 }
1421 else
1422 {
1423 MInvTransNormal(N, plane->Normal_Vector, plane->Trans);
1424
1425 MInvTransPoint(P, N, plane->Trans);
1426
1427 d = -plane->Distance - P[X] * N[X] - P[Y] * N[Y] - P[Z] * N[Z] + 1.0;
1428 }
1429
1430 Min[X] = Min[Y] = Min[Z] = -BOUND_HUGE/2;
1431 Max[X] = Max[Y] = Max[Z] = BOUND_HUGE/2;
1432
1433 /* y-z-plane */
1434
1435 if ((fabs(N[Y]) < EPSILON) && (fabs(N[Z]) < EPSILON)) // [CLi] can't test for fabs(1-N[X])<EPSILON because N isn't normalized
1436 {
1437 if (N[X] > 0.0)
1438 {
1439 Max[X] = d;
1440 }
1441 else
1442 {
1443 Min[X] = -d;
1444 }
1445 }
1446
1447 /* x-z-plane */
1448
1449 if ((fabs(N[X]) < EPSILON) && (fabs(N[Z]) < EPSILON)) // [CLi] can't test for fabs(1-N[Y])<EPSILON because N isn't normalized
1450 {
1451 if (N[Y] > 0.0)
1452 {
1453 Max[Y] = d;
1454 }
1455 else
1456 {
1457 Min[Y] = -d;
1458 }
1459 }
1460
1461 /* x-y-plane */
1462
1463 if ((fabs(N[X]) < EPSILON) && (fabs(N[Y]) < EPSILON)) // [CLi] can't test for fabs(1-N[Z])<EPSILON because N isn't normalized
1464 {
1465 if (N[Z] > 0.0)
1466 {
1467 Max[Z] = d;
1468 }
1469 else
1470 {
1471 Min[Z] = -d;
1472 }
1473 }
1474 }
1475
1476 }
1477