1 /****************************************************************************
2 * polygon.cpp
3 *
4 * This module implements functions that manipulate polygons.
5 *
6 * This module was written by Dieter Bayer [DB].
7 *
8 * from Persistence of Vision(tm) Ray Tracer version 3.6.
9 * Copyright 1991-2003 Persistence of Vision Team
10 * Copyright 2003-2004 Persistence of Vision Raytracer Pty. Ltd.
11 *---------------------------------------------------------------------------
12 * NOTICE: This source code file is provided so that users may experiment
13 * with enhancements to POV-Ray and to port the software to platforms other
14 * than those supported by the POV-Ray developers. There are strict rules
15 * regarding how you are permitted to use this file. These rules are contained
16 * in the distribution and derivative versions licenses which should have been
17 * provided with this file.
18 *
19 * These licences may be found online, linked from the end-user license
20 * agreement that is located at http://www.povray.org/povlegal.html
21 *---------------------------------------------------------------------------
22 * This program is based on the popular DKB raytracer version 2.12.
23 * DKBTrace was originally written by David K. Buck.
24 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
25 *---------------------------------------------------------------------------
26 * $File: //depot/povray/3.6-release/source/polygon.cpp $
27 * $Revision: #3 $
28 * $Change: 3032 $
29 * $DateTime: 2004/08/02 18:43:41 $
30 * $Author: chrisc $
31 * $Log$
32 *****************************************************************************/
33
34 /****************************************************************************
35 *
36 * Explanation:
37 *
38 * Syntax:
39 *
40 * ---
41 *
42 * The "inside polygon"-test was taken from:
43 *
44 * E. Haines, "An Introduction to Ray-Tracing", ..., pp. 53-59
45 *
46 * ---
47 *
48 * May 1994 : Creation. [DB]
49 *
50 * Oct 1994 : Changed polygon structure. Polygon points are now stored
51 * in a seperate structure. This - together with the use of
52 * a transformation - allows to keep just one point definition
53 * for multiple copies of one polygon. [DB]
54 *
55 *****************************************************************************/
56
57 #include "frame.h"
58 #include "vector.h"
59 #include "matrices.h"
60 #include "objects.h"
61 #include "polygon.h"
62 #include "povray.h"
63
64 #include <algorithm>
65
66 BEGIN_POV_NAMESPACE
67
68 /*****************************************************************************
69 * Local preprocessor defines
70 ******************************************************************************/
71
72 /* Minimal intersection depth for a valid intersection. */
73 const DBL DEPTH_TOLERANCE = 1.0e-8;
74
75 /* If |x| < ZERO_TOLERANCE x is assumed to be 0. */
76 const DBL ZERO_TOLERANCE = 1.0e-10;
77
78
79
80
81 /*****************************************************************************
82 * Static functions
83 ******************************************************************************/
84
85 static int intersect_poylgon (RAY *Ray, POLYGON *Polyg, DBL *Depth);
86 static int in_polygon (int Number, UV_VECT *Points, DBL u, DBL v);
87 static int All_Polygon_Intersections (OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack);
88 static int Inside_Polygon (VECTOR point, OBJECT *Object);
89 static void Polygon_Normal (VECTOR Result, OBJECT *Object, INTERSECTION *Inter);
90 static POLYGON *Copy_Polygon (OBJECT *Object);
91 static void Translate_Polygon (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
92 static void Rotate_Polygon (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
93 static void Scale_Polygon (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
94 static void Transform_Polygon (OBJECT *Object, TRANSFORM *Trans);
95 static void Invert_Polygon (OBJECT *Object);
96 static void Destroy_Polygon (OBJECT *Object);
97 static void Compute_Polygon_BBox (POLYGON *Polyg);
98
99
100
101 /*****************************************************************************
102 * Local variables
103 ******************************************************************************/
104
105 static METHODS Polygon_Methods =
106 {
107 All_Polygon_Intersections,
108 Inside_Polygon, Polygon_Normal, Default_UVCoord,
109 (COPY_METHOD)Copy_Polygon,
110 Translate_Polygon, Rotate_Polygon,
111 Scale_Polygon, Transform_Polygon, Invert_Polygon, Destroy_Polygon
112 };
113
114
115
116 /*****************************************************************************
117 *
118 * FUNCTION
119 *
120 * All_Polygon_Intersections
121 *
122 * INPUT
123 *
124 * Object - Object
125 * Ray - Ray
126 * Depth_Stack - Intersection stack
127 *
128 * OUTPUT
129 *
130 * Depth_Stack
131 *
132 * RETURNS
133 *
134 * int - true, if a intersection was found
135 *
136 * AUTHOR
137 *
138 * Dieter Bayer
139 *
140 * DESCRIPTION
141 *
142 * Determine ray/polygon intersection and clip intersection found.
143 *
144 * CHANGES
145 *
146 * May 1994 : Creation.
147 *
148 ******************************************************************************/
149
All_Polygon_Intersections(OBJECT * Object,RAY * Ray,ISTACK * Depth_Stack)150 static int All_Polygon_Intersections(OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack)
151 {
152 DBL Depth;
153 VECTOR IPoint;
154
155 if (intersect_poylgon(Ray, (POLYGON *)Object, &Depth))
156 {
157 VEvaluateRay(IPoint, Ray->Initial, Depth, Ray->Direction);
158
159 if (Point_In_Clip(IPoint, Object->Clip))
160 {
161 push_entry(Depth, IPoint, Object, Depth_Stack);
162
163 return(true);
164 }
165 }
166
167 return(false);
168 }
169
170
171
172 /*****************************************************************************
173 *
174 * FUNCTION
175 *
176 * intersect_poylgon
177 *
178 * INPUT
179 *
180 * Ray - Ray
181 * Polyg - Polygon
182 * Depth - Depth of intersection found
183 *
184 * OUTPUT
185 *
186 * Depth
187 *
188 * RETURNS
189 *
190 * int - true if intersection found
191 *
192 * AUTHOR
193 *
194 * Dieter Bayer
195 *
196 * DESCRIPTION
197 *
198 * Determine ray/polygon intersection.
199 *
200 * CHANGES
201 *
202 * May 1994 : Creation.
203 *
204 ******************************************************************************/
205
intersect_poylgon(RAY * Ray,POLYGON * Polyg,DBL * Depth)206 static int intersect_poylgon(RAY *Ray, POLYGON *Polyg, DBL *Depth)
207 {
208 DBL x, y, len;
209 VECTOR p, d;
210
211 /* Don't test degenerate polygons. */
212
213 if (Test_Flag(Polyg, DEGENERATE_FLAG))
214 {
215 return(false);
216 }
217
218 Increase_Counter(stats[Ray_Polygon_Tests]);
219
220 /* Transform the ray into the polygon space. */
221
222 MInvTransPoint(p, Ray->Initial, Polyg->Trans);
223
224 MInvTransDirection(d, Ray->Direction, Polyg->Trans);
225
226 VLength(len, d);
227
228 VInverseScaleEq(d, len);
229
230 /* Intersect ray with the plane in which the polygon lies. */
231
232 if (fabs(d[Z]) < ZERO_TOLERANCE)
233 {
234 return(false);
235 }
236
237 *Depth = -p[Z] / d[Z];
238
239 if ((*Depth < DEPTH_TOLERANCE) || (*Depth > Max_Distance))
240 {
241 return(false);
242 }
243
244 /* Does the intersection point lie inside the polygon? */
245
246 x = p[X] + *Depth * d[X];
247 y = p[Y] + *Depth * d[Y];
248
249 if (in_polygon(Polyg->Data->Number, Polyg->Data->Points, x, y))
250 {
251 Increase_Counter(stats[Ray_Polygon_Tests_Succeeded]);
252
253 *Depth /= len;
254
255 return (true);
256 }
257 else
258 {
259 return (false);
260 }
261 }
262
263
264
265 /*****************************************************************************
266 *
267 * FUNCTION
268 *
269 * Inside_Polygon
270 *
271 * INPUT
272 *
273 * IPoint - Intersection point
274 * Object - Object
275 *
276 * OUTPUT
277 *
278 * RETURNS
279 *
280 * int - always false
281 *
282 * AUTHOR
283 *
284 * Dieter Bayer
285 *
286 * DESCRIPTION
287 *
288 * Polygons can't be used in CSG.
289 *
290 * CHANGES
291 *
292 * May 1994 : Creation.
293 *
294 ******************************************************************************/
295
Inside_Polygon(VECTOR,OBJECT *)296 static int Inside_Polygon(VECTOR, OBJECT *)
297 {
298 return(false);
299 }
300
301
302
303 /*****************************************************************************
304 *
305 * FUNCTION
306 *
307 * Polygon_Normal
308 *
309 * INPUT
310 *
311 * Result - Normal vector
312 * Object - Object
313 * Inter - Intersection found
314 *
315 * OUTPUT
316 *
317 * Result
318 *
319 * RETURNS
320 *
321 * AUTHOR
322 *
323 * Dieter Bayer
324 *
325 * DESCRIPTION
326 *
327 * Calculate the normal of the polygon in a given point.
328 *
329 * CHANGES
330 *
331 * May 1994 : Creation.
332 *
333 ******************************************************************************/
334
Polygon_Normal(VECTOR Result,OBJECT * Object,INTERSECTION *)335 static void Polygon_Normal(VECTOR Result, OBJECT *Object, INTERSECTION *)
336 {
337 Assign_Vector(Result, ((POLYGON *)Object)->S_Normal);
338 }
339
340
341
342 /*****************************************************************************
343 *
344 * FUNCTION
345 *
346 * Translate_Polygon
347 *
348 * INPUT
349 *
350 * Object - Object
351 * Vector - Translation vector
352 *
353 * OUTPUT
354 *
355 * Object
356 *
357 * RETURNS
358 *
359 * AUTHOR
360 *
361 * Dieter Bayer
362 *
363 * DESCRIPTION
364 *
365 * Translate a polygon.
366 *
367 * CHANGES
368 *
369 * May 1994 : Creation.
370 *
371 ******************************************************************************/
372
Translate_Polygon(OBJECT * Object,VECTOR,TRANSFORM * Trans)373 static void Translate_Polygon(OBJECT *Object, VECTOR, TRANSFORM *Trans)
374 {
375 Transform_Polygon(Object, Trans);
376 }
377
378
379
380 /*****************************************************************************
381 *
382 * FUNCTION
383 *
384 * Rotate_Polygon
385 *
386 * INPUT
387 *
388 * Object - Object
389 * Vector - Rotation vector
390 *
391 * OUTPUT
392 *
393 * Object
394 *
395 * RETURNS
396 *
397 * AUTHOR
398 *
399 * Dieter Bayer
400 *
401 * DESCRIPTION
402 *
403 * Rotate a polygon.
404 *
405 * CHANGES
406 *
407 * May 1994 : Creation.
408 *
409 ******************************************************************************/
410
Rotate_Polygon(OBJECT * Object,VECTOR,TRANSFORM * Trans)411 static void Rotate_Polygon(OBJECT *Object, VECTOR, TRANSFORM *Trans)
412 {
413 Transform_Polygon(Object, Trans);
414 }
415
416
417
418 /*****************************************************************************
419 *
420 * FUNCTION
421 *
422 * Scale_Polygon
423 *
424 * INPUT
425 *
426 * Object - Object
427 * Vector - Scaling vector
428 *
429 * OUTPUT
430 *
431 * Object
432 *
433 * RETURNS
434 *
435 * AUTHOR
436 *
437 * Dieter Bayer
438 *
439 * DESCRIPTION
440 *
441 * Scale a polygon.
442 *
443 * CHANGES
444 *
445 * May 1994 : Creation.
446 *
447 ******************************************************************************/
448
Scale_Polygon(OBJECT * Object,VECTOR,TRANSFORM * Trans)449 static void Scale_Polygon(OBJECT *Object, VECTOR, TRANSFORM *Trans)
450 {
451 Transform_Polygon(Object, Trans);
452 }
453
454
455
456 /*****************************************************************************
457 *
458 * FUNCTION
459 *
460 * Transform_Polygon
461 *
462 * INPUT
463 *
464 * Object - Object
465 * Trans - Transformation to apply
466 *
467 * OUTPUT
468 *
469 * Object
470 *
471 * RETURNS
472 *
473 * AUTHOR
474 *
475 * Dieter Bayer
476 *
477 * DESCRIPTION
478 *
479 * Transform a polygon by transforming all points
480 * and recalculating the polygon.
481 *
482 * CHANGES
483 *
484 * May 1994 : Creation.
485 *
486 ******************************************************************************/
487
Transform_Polygon(OBJECT * Object,TRANSFORM * Trans)488 static void Transform_Polygon(OBJECT *Object, TRANSFORM *Trans)
489 {
490 VECTOR N;
491 POLYGON *Polyg = (POLYGON *)Object;
492
493 Compose_Transforms(Polyg->Trans, Trans);
494
495 Make_Vector(N, 0.0, 0.0, 1.0);
496 MTransNormal(Polyg->S_Normal, N, Polyg->Trans);
497
498 VNormalizeEq(Polyg->S_Normal);
499
500 Compute_Polygon_BBox(Polyg);
501 }
502
503
504
505 /*****************************************************************************
506 *
507 * FUNCTION
508 *
509 * Invert_Polygon
510 *
511 * INPUT
512 *
513 * Object - Object
514 *
515 * OUTPUT
516 *
517 * Object
518 *
519 * RETURNS
520 *
521 * AUTHOR
522 *
523 * Dieter Bayer
524 *
525 * DESCRIPTION
526 *
527 * Invert a polygon.
528 *
529 * CHANGES
530 *
531 * May 1994 : Creation.
532 *
533 ******************************************************************************/
534
Invert_Polygon(OBJECT *)535 static void Invert_Polygon(OBJECT *)
536 {
537 }
538
539
540
541 /*****************************************************************************
542 *
543 * FUNCTION
544 *
545 * Create_Polygon
546 *
547 * INPUT
548 *
549 * OUTPUT
550 *
551 * RETURNS
552 *
553 * POLYGON * - new polygon
554 *
555 * AUTHOR
556 *
557 * Dieter Bayer
558 *
559 * DESCRIPTION
560 *
561 * Create a new polygon.
562 *
563 * CHANGES
564 *
565 * May 1994 : Creation.
566 *
567 ******************************************************************************/
568
Create_Polygon()569 POLYGON *Create_Polygon()
570 {
571 POLYGON *New;
572
573 New = (POLYGON *)POV_MALLOC(sizeof(POLYGON), "polygon");
574
575 INIT_OBJECT_FIELDS(New,POLYGON_OBJECT,&Polygon_Methods)
576
577 New->Trans = Create_Transform();
578
579 Make_Vector(New->S_Normal, 0.0, 0.0, 1.0);
580
581 New->Data = NULL;
582
583 return (New);
584 }
585
586
587
588 /*****************************************************************************
589 *
590 * FUNCTION
591 *
592 * Copy_Polygon
593 *
594 * INPUT
595 *
596 * Object - Object
597 *
598 * OUTPUT
599 *
600 * RETURNS
601 *
602 * void * - New polygon
603 *
604 * AUTHOR
605 *
606 * Dieter Bayer
607 *
608 * DESCRIPTION
609 *
610 * Copy a polygon structure.
611 *
612 * CHANGES
613 *
614 * May 1994 : Creation.
615 *
616 ******************************************************************************/
617
Copy_Polygon(OBJECT * Object)618 static POLYGON *Copy_Polygon(OBJECT *Object)
619 {
620 POLYGON *New, *Polyg = (POLYGON *)Object;
621
622 New = Create_Polygon ();
623
624 /* Get rid of the transformation created in Create_Polygon(). */
625
626 Destroy_Transform(New->Trans);
627
628 /* Copy polygon. */
629
630 *New = *Polyg;
631
632 New->Trans = Copy_Transform(Polyg->Trans);
633
634 New->Data->References++;
635
636 return (New);
637 }
638
639
640
641 /*****************************************************************************
642 *
643 * FUNCTION
644 *
645 * Destroy_Polygon
646 *
647 * INPUT
648 *
649 * Object - Object
650 *
651 * OUTPUT
652 *
653 * Object
654 *
655 * RETURNS
656 *
657 * AUTHOR
658 *
659 * Dieter Bayer
660 *
661 * DESCRIPTION
662 *
663 * Destroy a polygon.
664 *
665 * CHANGES
666 *
667 * May 1994 : Creation.
668 *
669 * Dec 1994 : Fixed memory leakage. [DB]
670 *
671 ******************************************************************************/
672
Destroy_Polygon(OBJECT * Object)673 static void Destroy_Polygon(OBJECT *Object)
674 {
675 POLYGON *Polyg = (POLYGON *)Object;
676
677 if (--(Polyg->Data->References) == 0)
678 {
679 POV_FREE (Polyg->Data->Points);
680
681 POV_FREE (Polyg->Data);
682 }
683
684 Destroy_Transform(Polyg->Trans);
685
686 POV_FREE (Object);
687 }
688
689
690
691 /*****************************************************************************
692 *
693 * FUNCTION
694 *
695 * Compute_Polygon
696 *
697 * INPUT
698 *
699 * Polyg - Polygon
700 * Points - 3D points describing the polygon
701 *
702 * OUTPUT
703 *
704 * Polyg
705 *
706 * RETURNS
707 *
708 * AUTHOR
709 *
710 * Dieter Bayer
711 *
712 * DESCRIPTION
713 *
714 * Compute the following things for a given polygon:
715 *
716 * - Polygon transformation
717 *
718 * - Array of 2d points describing the shape of the polygon
719 *
720 * CHANGES
721 *
722 * May 1994 : Creation.
723 *
724 ******************************************************************************/
725
Compute_Polygon(POLYGON * Polyg,int Number,VECTOR * Points)726 void Compute_Polygon(POLYGON *Polyg, int Number, VECTOR *Points)
727 {
728 int i;
729 DBL x, y, z, d;
730 VECTOR o, u, v, w, N;
731 MATRIX a, b;
732
733 /* Create polygon data. */
734
735 if (Polyg->Data == NULL)
736 {
737 Polyg->Data = (POLYGON_DATA *)POV_MALLOC(sizeof(POLYGON_DATA), "polygon points");
738
739 Polyg->Data->References = 1;
740
741 Polyg->Data->Number = Number;
742
743 Polyg->Data->Points = (UV_VECT *)POV_MALLOC(Number*sizeof(UV_VECT), "polygon points");
744 }
745 else
746 {
747 Error("Polygon data already computed.");
748 }
749
750 /* Get polygon's coordinate system (one of the many possible) */
751
752 Assign_Vector(o, Points[0]);
753
754 /* Find valid, i.e. non-zero u vector. */
755
756 for (i = 1; i < Number; i++)
757 {
758 VSub(u, Points[i], o);
759
760 if (VSumSqr(u) > EPSILON)
761 {
762 break;
763 }
764 }
765
766 if (i == Number)
767 {
768 Set_Flag(Polyg, DEGENERATE_FLAG);
769
770 Warning(0, "Points in polygon are co-linear. Ignoring polygon.");
771 }
772
773 /* Find valid, i.e. non-zero v and w vectors. */
774
775 for (i++; i < Number; i++)
776 {
777 VSub(v, Points[i], o);
778
779 VCross(w, u, v);
780
781 if ((VSumSqr(v) > EPSILON) && (VSumSqr(w) > EPSILON))
782 {
783 break;
784 }
785 }
786
787 if (i == Number)
788 {
789 Set_Flag(Polyg, DEGENERATE_FLAG);
790
791 Warning(0, "Points in polygon are co-linear. Ignoring polygon.");
792 }
793
794 VCross(u, v, w);
795 VCross(v, w, u);
796
797 VNormalize(u, u);
798 VNormalize(v, v);
799 VNormalize(w, w);
800
801 MIdentity(a);
802 MIdentity(b);
803
804 a[3][0] = -o[X];
805 a[3][1] = -o[Y];
806 a[3][2] = -o[Z];
807
808 b[0][0] = u[X];
809 b[1][0] = u[Y];
810 b[2][0] = u[Z];
811
812 b[0][1] = v[X];
813 b[1][1] = v[Y];
814 b[2][1] = v[Z];
815
816 b[0][2] = w[X];
817 b[1][2] = w[Y];
818 b[2][2] = w[Z];
819
820 MTimesC(Polyg->Trans->inverse, a, b);
821
822 MInvers(Polyg->Trans->matrix, Polyg->Trans->inverse);
823
824 /* Project points onto the u,v-plane (3D --> 2D) */
825
826 for (i = 0; i < Number; i++)
827 {
828 x = Points[i][X] - o[X];
829 y = Points[i][Y] - o[Y];
830 z = Points[i][Z] - o[Z];
831
832 d = x * w[X] + y * w[Y] + z * w[Z];
833
834 if (fabs(d) > ZERO_TOLERANCE)
835 {
836 Set_Flag(Polyg, DEGENERATE_FLAG);
837
838 Warning(0, "Points in polygon are not co-planar. Ignoring polygons.");
839 }
840
841 Polyg->Data->Points[i][X] = x * u[X] + y * u[Y] + z * u[Z];
842 Polyg->Data->Points[i][Y] = x * v[X] + y * v[Y] + z * v[Z];
843 }
844
845 Make_Vector(N, 0.0, 0.0, 1.0);
846 MTransNormal(Polyg->S_Normal, N, Polyg->Trans);
847
848 VNormalizeEq(Polyg->S_Normal);
849
850 Compute_Polygon_BBox(Polyg);
851 }
852
853
854
855 /*****************************************************************************
856 *
857 * FUNCTION
858 *
859 * Compute_Polygon_BBox
860 *
861 * INPUT
862 *
863 * Polyg - Polygon
864 *
865 * OUTPUT
866 *
867 * Polyg
868 *
869 * RETURNS
870 *
871 * AUTHOR
872 *
873 * Dieter Bayer
874 *
875 * DESCRIPTION
876 *
877 * Calculate the bounding box of a polygon.
878 *
879 * CHANGES
880 *
881 * May 1994 : Creation.
882 *
883 ******************************************************************************/
884
Compute_Polygon_BBox(POLYGON * Polyg)885 static void Compute_Polygon_BBox(POLYGON *Polyg)
886 {
887 int i;
888 VECTOR p, Puv, Min, Max;
889
890 Min[X] = Min[Y] = Min[Z] = BOUND_HUGE;
891 Max[X] = Max[Y] = Max[Z] = -BOUND_HUGE;
892
893 for (i = 0; i < Polyg->Data->Number; i++)
894 {
895 Puv[X] = Polyg->Data->Points[i][X];
896 Puv[Y] = Polyg->Data->Points[i][Y];
897 Puv[Z] = 0.0;
898
899 MTransPoint(p, Puv, Polyg->Trans);
900
901 Min[X] = min(Min[X], p[X]);
902 Min[Y] = min(Min[Y], p[Y]);
903 Min[Z] = min(Min[Z], p[Z]);
904 Max[X] = max(Max[X], p[X]);
905 Max[Y] = max(Max[Y], p[Y]);
906 Max[Z] = max(Max[Z], p[Z]);
907 }
908
909 Make_BBox_from_min_max(Polyg->BBox, Min, Max);
910
911 if (fabs(Polyg->BBox.Lengths[X]) < Small_Tolerance)
912 {
913 Polyg->BBox.Lower_Left[X] -= Small_Tolerance;
914 Polyg->BBox.Lengths[X] += 2.0 * Small_Tolerance;
915 }
916
917 if (fabs(Polyg->BBox.Lengths[Y]) < Small_Tolerance)
918 {
919 Polyg->BBox.Lower_Left[Y] -= Small_Tolerance;
920 Polyg->BBox.Lengths[Y] += 2.0 * Small_Tolerance;
921 }
922
923 if (fabs(Polyg->BBox.Lengths[Z]) < Small_Tolerance)
924 {
925 Polyg->BBox.Lower_Left[Z] -= Small_Tolerance;
926 Polyg->BBox.Lengths[Z] += 2.0 * Small_Tolerance;
927 }
928 }
929
930
931
932 /*****************************************************************************
933 *
934 * FUNCTION
935 *
936 * in_polygon
937 *
938 * INPUT
939 *
940 * Number - Number of points
941 * Points - Points describing polygon's shape
942 * u, v - 2D-coordinates of the point to test
943 *
944 * OUTPUT
945 *
946 * RETURNS
947 *
948 * int - true, if inside
949 *
950 * AUTHOR
951 *
952 * Eric Haines, 3D/Eye Inc, erich@eye.com
953 *
954 * DESCRIPTION
955 *
956 * ======= Crossings Multiply algorithm ===================================
957 *
958 * This version is usually somewhat faster than the original published in
959 * Graphics Gems IV; by turning the division for testing the X axis crossing
960 * into a tricky multiplication test this part of the test became faster,
961 * which had the additional effect of making the test for "both to left or
962 * both to right" a bit slower for triangles than simply computing the
963 * intersection each time. The main increase is in triangle testing speed,
964 * which was about 15% faster; all other polygon complexities were pretty much
965 * the same as before. On machines where division is very expensive (not the
966 * case on the HP 9000 series on which I tested) this test should be much
967 * faster overall than the old code. Your mileage may (in fact, will) vary,
968 * depending on the machine and the test data, but in general I believe this
969 * code is both shorter and faster. This test was inspired by unpublished
970 * Graphics Gems submitted by Joseph Samosky and Mark Haigh-Hutchinson.
971 * Related work by Samosky is in:
972 *
973 * Samosky, Joseph, "SectionView: A system for interactively specifying and
974 * visualizing sections through three-dimensional medical image data",
975 * M.S. Thesis, Department of Electrical Engineering and Computer Science,
976 * Massachusetts Institute of Technology, 1993.
977 *
978 *
979 * Shoot a test ray along +X axis. The strategy is to compare vertex Y values
980 * to the testing point's Y and quickly discard edges which are entirely to one
981 * side of the test ray.
982 *
983 * CHANGES
984 *
985 * -
986 *
987 ******************************************************************************/
988
in_polygon(int Number,UV_VECT * Points,DBL u,DBL v)989 static int in_polygon(int Number, UV_VECT *Points, DBL u, DBL v)
990 {
991 register int i, yflag0, yflag1, inside_flag;
992 register DBL ty, tx;
993 register DBL *vtx0, *vtx1, *first;
994
995 tx = u;
996 ty = v;
997
998 vtx0 = &Points[0][X];
999 vtx1 = &Points[1][X];
1000
1001 first = vtx0;
1002
1003 /* get test bit for above/below X axis */
1004
1005 yflag0 = (vtx0[Y] >= ty);
1006
1007 inside_flag = false;
1008
1009 for (i = 1; i < Number; )
1010 {
1011 yflag1 = (vtx1[Y] >= ty);
1012
1013 /*
1014 * Check if endpoints straddle (are on opposite sides) of X axis
1015 * (i.e. the Y's differ); if so, +X ray could intersect this edge.
1016 * The old test also checked whether the endpoints are both to the
1017 * right or to the left of the test point. However, given the faster
1018 * intersection point computation used below, this test was found to
1019 * be a break-even proposition for most polygons and a loser for
1020 * triangles (where 50% or more of the edges which survive this test
1021 * will cross quadrants and so have to have the X intersection computed
1022 * anyway). I credit Joseph Samosky with inspiring me to try dropping
1023 * the "both left or both right" part of my code.
1024 */
1025
1026 if (yflag0 != yflag1)
1027 {
1028 /*
1029 * Check intersection of pgon segment with +X ray.
1030 * Note if >= point's X; if so, the ray hits it.
1031 * The division operation is avoided for the ">=" test by checking
1032 * the sign of the first vertex wrto the test point; idea inspired
1033 * by Joseph Samosky's and Mark Haigh-Hutchinson's different
1034 * polygon inclusion tests.
1035 */
1036
1037 if (((vtx1[Y]-ty) * (vtx0[X]-vtx1[X]) >= (vtx1[X]-tx) * (vtx0[Y]-vtx1[Y])) == yflag1)
1038 {
1039 inside_flag = !inside_flag;
1040 }
1041 }
1042
1043 /* Move to the next pair of vertices, retaining info as possible. */
1044
1045 if ((i < Number-2) && (vtx1[X] == first[X]) && (vtx1[Y] == first[Y]))
1046 {
1047 vtx0 = &Points[++i][X];
1048 vtx1 = &Points[++i][X];
1049
1050 yflag0 = (vtx0[Y] >= ty);
1051
1052 first = vtx0;
1053 }
1054 else
1055 {
1056 vtx0 = vtx1;
1057 vtx1 = &Points[++i][X];
1058
1059 yflag0 = yflag1;
1060 }
1061 }
1062
1063 return(inside_flag);
1064 }
1065
1066 END_POV_NAMESPACE
1067