more eol-style
[xonotic/netradiant.git] / radiant / winding.cpp
1 /*
2 Copyright (C) 1999-2007 id Software, Inc. and contributors.
3 For a list of contributors, see the accompanying CONTRIBUTORS file.
4
5 This file is part of GtkRadiant.
6
7 GtkRadiant is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 GtkRadiant is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with GtkRadiant; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
20 */
21
22
23
24 #include "stdafx.h"
25 #include <assert.h>
26 #include "winding.h"
27
28 #define BOGUS_RANGE (g_MaxWorldCoord+1)
29
30 /*
31 =============
32 Plane_Equal
33 =============
34 */
35 #define NORMAL_EPSILON  0.0001
36 #define DIST_EPSILON    0.02
37
38 int Plane_Equal(plane_t *a, plane_t *b, int flip)
39 {
40         vec3_t normal;
41         float dist;
42
43         if (flip) {
44                 normal[0] = - b->normal[0];
45                 normal[1] = - b->normal[1];
46                 normal[2] = - b->normal[2];
47                 dist = - b->dist;
48         }
49         else {
50                 normal[0] = b->normal[0];
51                 normal[1] = b->normal[1];
52                 normal[2] = b->normal[2];
53                 dist = b->dist;
54         }
55         if (
56            fabs(a->normal[0] - normal[0]) < NORMAL_EPSILON
57         && fabs(a->normal[1] - normal[1]) < NORMAL_EPSILON
58         && fabs(a->normal[2] - normal[2]) < NORMAL_EPSILON
59         && fabs(a->dist - dist) < DIST_EPSILON )
60                 return true;
61         return false;
62 }
63
64 /*
65 ============
66 Plane_FromPoints
67 ============
68 */
69 int Plane_FromPoints(vec3_t p1, vec3_t p2, vec3_t p3, plane_t *plane)
70 {
71         vec3_t v1, v2;
72
73         VectorSubtract(p2, p1, v1);
74         VectorSubtract(p3, p1, v2);
75         //CrossProduct(v2, v1, plane->normal);
76         CrossProduct(v1, v2, plane->normal);
77         if (VectorNormalize(plane->normal, plane->normal) < 0.1) return false;
78         plane->dist = DotProduct(p1, plane->normal);
79         return true;
80 }
81
82 /*
83 =================
84 Point_Equal
85 =================
86 */
87 int Point_Equal(vec3_t p1, vec3_t p2, float epsilon)
88 {
89         int i;
90
91         for (i = 0; i < 3; i++)
92         {
93                 if (fabs(p1[i] - p2[i]) > epsilon) return false;
94         }
95         return true;
96 }
97
98
99 /*
100 =================
101 Winding_BaseForPlane
102 =================
103 */
104 //#define DBG_WNDG
105 winding_t *Winding_BaseForPlane (plane_t *p)
106 {
107         int             i, x;
108         vec_t   max, v;
109         vec3_t  org, vright, vup;
110         winding_t       *w;
111         
112         // find the major axis
113 #ifdef DBG_WNDG
114   Sys_Printf("Winding_BaseForPlane %p\n",p);
115 #endif
116
117   max = -BOGUS_RANGE;
118         x = -1;
119         for (i=0 ; i<3; i++)
120         {
121                 v = fabs(p->normal[i]);
122                 if (v > max)
123                 {
124                         x = i;
125                         max = v;
126                 }
127         }
128         if (x==-1)
129                 Error ("Winding_BaseForPlane: no axis found");
130                 
131         VectorCopy (vec3_origin, vup);  
132         switch (x)
133         {
134         case 0:
135         case 1:
136                 vup[2] = 1;
137                 break;          
138         case 2:
139                 vup[0] = 1;
140                 break;          
141         }
142
143
144         v = DotProduct (vup, p->normal);
145         VectorMA (vup, -v, p->normal, vup);
146         VectorNormalize (vup, vup);
147                 
148         VectorScale (p->normal, p->dist, org);
149         
150         CrossProduct (vup, p->normal, vright);
151         
152         VectorScale (vup, BOGUS_RANGE, vup);
153         VectorScale (vright, BOGUS_RANGE, vright);
154
155   // project a really big       axis aligned box onto the plane
156         w = Winding_Alloc (4);
157         
158         VectorSubtract (org, vright, w->points[0]);
159         VectorAdd (w->points[0], vup, w->points[0]);
160         
161         VectorAdd (org, vright, w->points[1]);
162         VectorAdd (w->points[1], vup, w->points[1]);
163         
164         VectorAdd (org, vright, w->points[2]);
165         VectorSubtract (w->points[2], vup, w->points[2]);
166         
167         VectorSubtract (org, vright, w->points[3]);
168         VectorSubtract (w->points[3], vup, w->points[3]);
169         
170         w->numpoints = 4;
171
172         return w;       
173 }
174
175 // macro to compute winding size
176 #define WINDING_SIZE(pt) (sizeof(int)*2+sizeof(float)*5*(pt))
177
178 /*
179 ==================
180 Winding_Alloc
181 ==================
182 */
183 winding_t *Winding_Alloc (int points)
184 {
185         winding_t       *w;
186         int                     size;
187         
188         if (points > MAX_POINTS_ON_WINDING)
189                 Error ("Winding_Alloc: %i points", points);
190         
191 //      size = (int)((winding_t *)0)->points[points];
192   size = WINDING_SIZE(points);
193         w = (winding_t*) malloc (size);
194         memset (w, 0, size);
195         w->maxpoints = points;
196         
197         return w;
198 }
199
200 void Winding_Free (winding_t *w)
201 {
202         free(w);
203 }
204
205 /*
206 ==================
207 Winding_Clone
208 ==================
209 */
210 winding_t *Winding_Clone(winding_t *w)
211 {
212         int                     size;
213         winding_t       *c;
214         
215 //      size = (int)((winding_t *)0)->points[w->numpoints];
216   size = WINDING_SIZE(w->numpoints);
217         c = (winding_t*)qmalloc (size);
218         memcpy (c, w, size);
219         return c;
220 }
221
222 /*
223 ==================
224 ReverseWinding
225 ==================
226 */
227 winding_t *Winding_Reverse(winding_t *w)
228 {
229         int                     i;
230         winding_t       *c;
231
232         c = Winding_Alloc(w->numpoints);
233         for (i = 0; i < w->numpoints; i++)
234         {
235                 VectorCopy (w->points[w->numpoints-1-i], c->points[i]);
236         }
237         c->numpoints = w->numpoints;
238         return c;
239 }
240
241 /*
242 ==============
243 Winding_RemovePoint
244 ==============
245 */
246 void Winding_RemovePoint(winding_t *w, int point)
247 {
248         if (point < 0 || point >= w->numpoints)
249                 Error("Winding_RemovePoint: point out of range");
250
251         if (point < w->numpoints-1)
252         {
253                 memmove(&w->points[point], &w->points[point+1], (int)((winding_t *)0)->points[w->numpoints - point - 1]);
254         }
255         w->numpoints--;
256 }
257
258 /*
259 =============
260 Winding_InsertPoint
261 =============
262 */
263 winding_t *Winding_InsertPoint(winding_t *w, vec3_t point, int spot)
264 {
265         int i, j;
266         winding_t *neww;
267
268         if (spot > w->numpoints)
269         {
270                 Error("Winding_InsertPoint: spot > w->numpoints");
271         } //end if
272         if (spot < 0)
273         {
274                 Error("Winding_InsertPoint: spot < 0");
275         } //end if
276         neww = Winding_Alloc(w->numpoints + 1);
277         neww->numpoints = w->numpoints + 1;
278         for (i = 0, j = 0; i < neww->numpoints; i++)
279         {
280                 if (i == spot)
281                 {
282                         VectorCopy(point, neww->points[i]);
283                 }
284                 else
285                 {
286                         VectorCopy(w->points[j], neww->points[i]);
287                         j++;
288                 }
289         }
290         return neww;
291 }
292
293 /*
294 ==============
295 Winding_IsTiny
296 ==============
297 */
298 #define EDGE_LENGTH     0.2
299
300 int Winding_IsTiny (winding_t *w)
301 {
302         int             i, j;
303         vec_t   len;
304         vec3_t  delta;
305         int             edges;
306
307         edges = 0;
308         for (i=0 ; i<w->numpoints ; i++)
309         {
310                 j = i == w->numpoints - 1 ? 0 : i+1;
311                 VectorSubtract (w->points[j], w->points[i], delta);
312                 len = VectorLength (delta);
313                 if (len > EDGE_LENGTH)
314                 {
315                         if (++edges == 3)
316                                 return false;
317                 }
318         }
319         return true;
320 }
321
322 /*
323 ==============
324 Winding_IsHuge
325 ==============
326 */
327 int Winding_IsHuge(winding_t *w)
328 {
329         int             i, j;
330
331         for (i=0 ; i<w->numpoints ; i++)
332         {
333                 for (j=0 ; j<3 ; j++)
334                         if (w->points[i][j] < -BOGUS_RANGE+1 || w->points[i][j] > BOGUS_RANGE-1)
335                                 return true;
336         }
337         return false;
338 }
339
340 /*
341 =============
342 Winding_PlanesConcave
343 =============
344 */
345 #define WCONVEX_EPSILON         0.2
346
347 int Winding_PlanesConcave(winding_t *w1, winding_t *w2,
348                                                          vec3_t normal1, vec3_t normal2,
349                                                          float dist1, float dist2)
350 {
351         int i;
352
353         if (!w1 || !w2) return false;
354
355         // check if one of the points of winding 1 is at the back of the plane of winding 2
356         for (i = 0; i < w1->numpoints; i++)
357         {
358                 if (DotProduct(normal2, w1->points[i]) - dist2 > WCONVEX_EPSILON) return true;
359         }
360         // check if one of the points of winding 2 is at the back of the plane of winding 1
361         for (i = 0; i < w2->numpoints; i++)
362         {
363                 if (DotProduct(normal1, w2->points[i]) - dist1 > WCONVEX_EPSILON) return true;
364         }
365
366         return false;
367 }
368
369 /*
370 ==================
371 Winding_Clip
372
373 Clips the winding to the plane, returning the new winding on the positive side
374 Frees the input winding.
375 If keepon is true, an exactly on-plane winding will be saved, otherwise
376 it will be clipped away.
377 ==================
378 */
379 winding_t *Winding_Clip (winding_t *in, plane_t *split, qboolean keepon)
380 {
381         vec_t   dists[MAX_POINTS_ON_WINDING];
382         int             sides[MAX_POINTS_ON_WINDING];
383         int             counts[3];
384         vec_t   dot;
385         int             i, j;
386         vec_t   *p1, *p2;
387         vec3_t  mid;
388         winding_t       *neww;
389         int             maxpts;
390         
391         counts[0] = counts[1] = counts[2] = 0;
392
393         // determine sides for each point
394         for (i=0 ; i<in->numpoints ; i++)
395         {
396                 dot = DotProduct (in->points[i], split->normal);
397                 dot -= split->dist;
398                 dists[i] = dot;
399                 if (dot > ON_EPSILON)
400                         sides[i] = SIDE_FRONT;
401                 else if (dot < -ON_EPSILON)
402                         sides[i] = SIDE_BACK;
403                 else
404                 {
405                         sides[i] = SIDE_ON;
406                 }
407                 counts[sides[i]]++;
408         }
409         sides[i] = sides[0];
410         dists[i] = dists[0];
411         
412         if (keepon && !counts[0] && !counts[1])
413                 return in;
414                 
415         if (!counts[0])
416         {
417                 Winding_Free (in);
418                 return NULL;
419         }
420         if (!counts[1])
421                 return in;
422         
423         maxpts = in->numpoints+4;       // can't use counts[0]+2 because
424                                                                 // of fp grouping errors
425         neww = Winding_Alloc (maxpts);
426                 
427         for (i=0 ; i<in->numpoints ; i++)
428         {
429                 p1 = in->points[i];
430                 
431                 if (sides[i] == SIDE_ON)
432                 {
433                         VectorCopy (p1, neww->points[neww->numpoints]);
434                         neww->numpoints++;
435                         continue;
436                 }
437         
438                 if (sides[i] == SIDE_FRONT)
439                 {
440                         VectorCopy (p1, neww->points[neww->numpoints]);
441                         neww->numpoints++;
442                 }
443                 
444                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
445                         continue;
446                         
447                 // generate a split point
448                 p2 = in->points[(i+1)%in->numpoints];
449                 
450                 dot = dists[i] / (dists[i]-dists[i+1]);
451                 for (j=0 ; j<3 ; j++)
452                 {       // avoid round off error when possible
453                         if (split->normal[j] == 1)
454                                 mid[j] = split->dist;
455                         else if (split->normal[j] == -1)
456                                 mid[j] = -split->dist;
457                         else
458                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
459                 }
460                         
461                 VectorCopy (mid, neww->points[neww->numpoints]);
462                 neww->numpoints++;
463         }
464         
465         if (neww->numpoints > maxpts)
466                 Error ("Winding_Clip: points exceeded estimate");
467                 
468         // free the original winding
469         Winding_Free (in);
470         
471         return neww;
472 }
473
474 /*
475 =============
476 Winding_SplitEpsilon
477
478   split the input winding with the plane
479   the input winding stays untouched
480 =============
481 */
482 void Winding_SplitEpsilon (winding_t *in, vec3_t normal, double dist, 
483                                 vec_t epsilon, winding_t **front, winding_t **back)
484 {
485         vec_t   dists[MAX_POINTS_ON_WINDING+4];
486         int             sides[MAX_POINTS_ON_WINDING+4];
487         int             counts[3];
488         vec_t   dot;
489         int             i, j;
490         vec_t   *p1, *p2;
491         vec3_t  mid;
492         winding_t       *f, *b;
493         int             maxpts;
494         
495         counts[0] = counts[1] = counts[2] = 0;
496
497         // determine sides for each point
498         for (i = 0; i < in->numpoints; i++)
499         {
500                 dot = DotProduct (in->points[i], normal);
501                 dot -= dist;
502                 dists[i] = dot;
503                 if (dot > epsilon)
504                         sides[i] = SIDE_FRONT;
505                 else if (dot < -epsilon)
506                         sides[i] = SIDE_BACK;
507                 else
508                 {
509                         sides[i] = SIDE_ON;
510                 }
511                 counts[sides[i]]++;
512         }
513         sides[i] = sides[0];
514         dists[i] = dists[0];
515         
516         *front = *back = NULL;
517
518         if (!counts[0])
519         {
520                 *back = Winding_Clone(in);
521                 return;
522         }
523         if (!counts[1])
524         {
525                 *front = Winding_Clone(in);
526                 return;
527         }
528
529         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
530                                                                 // of fp grouping errors
531
532         *front = f = Winding_Alloc (maxpts);
533         *back = b = Winding_Alloc (maxpts);
534                 
535         for (i = 0; i < in->numpoints; i++)
536         {
537                 p1 = in->points[i];
538                 
539                 if (sides[i] == SIDE_ON)
540                 {
541                         VectorCopy (p1, f->points[f->numpoints]);
542                         f->numpoints++;
543                         VectorCopy (p1, b->points[b->numpoints]);
544                         b->numpoints++;
545                         continue;
546                 }
547         
548                 if (sides[i] == SIDE_FRONT)
549                 {
550                         VectorCopy (p1, f->points[f->numpoints]);
551                         f->numpoints++;
552                 }
553                 if (sides[i] == SIDE_BACK)
554                 {
555                         VectorCopy (p1, b->points[b->numpoints]);
556                         b->numpoints++;
557                 }
558
559                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
560                         continue;
561                         
562                 // generate a split point
563                 p2 = in->points[(i+1)%in->numpoints];
564                 
565                 dot = dists[i] / (dists[i]-dists[i+1]);
566                 for (j = 0; j < 3; j++)
567                 {
568                         // avoid round off error when possible
569                         if (normal[j] == 1)
570                                 mid[j] = dist;
571                         else if (normal[j] == -1)
572                                 mid[j] = -dist;
573                         else
574                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
575                 }
576                         
577                 VectorCopy (mid, f->points[f->numpoints]);
578                 f->numpoints++;
579                 VectorCopy (mid, b->points[b->numpoints]);
580                 b->numpoints++;
581         }
582         
583         if (f->numpoints > maxpts || b->numpoints > maxpts)
584                 Error ("Winding_Clip: points exceeded estimate");
585         if (f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING)
586                 Error ("Winding_Clip: MAX_POINTS_ON_WINDING");
587 }
588
589 /*
590 =============
591 Winding_TryMerge
592
593 If two windings share a common edge and the edges that meet at the
594 common points are both inside the other polygons, merge them
595
596 Returns NULL if the windings couldn't be merged, or the new winding.
597 The originals will NOT be freed.
598
599 if keep is true no points are ever removed
600 =============
601 */
602 #define CONTINUOUS_EPSILON      0.005
603
604 winding_t *Winding_TryMerge(winding_t *f1, winding_t *f2, vec3_t planenormal, int keep)
605 {
606         vec_t           *p1, *p2, *p3, *p4, *back;
607         winding_t       *newf;
608         int                     i, j, k, l;
609         vec3_t          normal, delta;
610         vec_t           dot;
611         qboolean        keep1, keep2;
612         
613
614         //
615         // find a common edge
616         //      
617         p1 = p2 = NULL; // stop compiler warning
618         j = 0;                  // 
619         
620         for (i = 0; i < f1->numpoints; i++)
621         {
622                 p1 = f1->points[i];
623                 p2 = f1->points[(i+1) % f1->numpoints];
624                 for (j = 0; j < f2->numpoints; j++)
625                 {
626                         p3 = f2->points[j];
627                         p4 = f2->points[(j+1) % f2->numpoints];
628                         for (k = 0; k < 3; k++)
629                         {
630                                 if (fabs(p1[k] - p4[k]) > 0.1)//EQUAL_EPSILON) //ME
631                                         break;
632                                 if (fabs(p2[k] - p3[k]) > 0.1)//EQUAL_EPSILON) //ME
633                                         break;
634                         } //end for
635                         if (k==3)
636                                 break;
637                 } //end for
638                 if (j < f2->numpoints)
639                         break;
640         } //end for
641         
642         if (i == f1->numpoints)
643                 return NULL;                    // no matching edges
644
645         //
646         // check slope of connected lines
647         // if the slopes are colinear, the point can be removed
648         //
649         back = f1->points[(i+f1->numpoints-1)%f1->numpoints];
650         VectorSubtract (p1, back, delta);
651         CrossProduct (planenormal, delta, normal);
652         VectorNormalize (normal, normal);
653         
654         back = f2->points[(j+2)%f2->numpoints];
655         VectorSubtract (back, p1, delta);
656         dot = DotProduct (delta, normal);
657         if (dot > CONTINUOUS_EPSILON)
658                 return NULL;                    // not a convex polygon
659         keep1 = (qboolean)(dot < -CONTINUOUS_EPSILON);
660         
661         back = f1->points[(i+2)%f1->numpoints];
662         VectorSubtract (back, p2, delta);
663         CrossProduct (planenormal, delta, normal);
664         VectorNormalize (normal, normal);
665
666         back = f2->points[(j+f2->numpoints-1)%f2->numpoints];
667         VectorSubtract (back, p2, delta);
668         dot = DotProduct (delta, normal);
669         if (dot > CONTINUOUS_EPSILON)
670                 return NULL;                    // not a convex polygon
671         keep2 = (qboolean)(dot < -CONTINUOUS_EPSILON);
672
673         //
674         // build the new polygon
675         //
676         newf = Winding_Alloc (f1->numpoints + f2->numpoints);
677         
678         // copy first polygon
679         for (k=(i+1)%f1->numpoints ; k != i ; k=(k+1)%f1->numpoints)
680         {
681                 if (!keep && k==(i+1)%f1->numpoints && !keep2)
682                         continue;
683                 
684                 VectorCopy (f1->points[k], newf->points[newf->numpoints]);
685                 newf->numpoints++;
686         }
687         
688         // copy second polygon
689         for (l= (j+1)%f2->numpoints ; l != j ; l=(l+1)%f2->numpoints)
690         {
691                 if (!keep && l==(j+1)%f2->numpoints && !keep1)
692                         continue;
693                 VectorCopy (f2->points[l], newf->points[newf->numpoints]);
694                 newf->numpoints++;
695         }
696
697         return newf;
698 }
699
700 /*
701 ============
702 Winding_Plane
703 ============
704 */
705 void Winding_Plane (winding_t *w, vec3_t normal, double *dist)
706 {
707         vec3_t v1, v2;
708         int i;
709
710         //find two vectors each longer than 0.5 units
711         for (i = 0; i < w->numpoints; i++)
712         {
713                 VectorSubtract(w->points[(i+1) % w->numpoints], w->points[i], v1);
714                 VectorSubtract(w->points[(i+2) % w->numpoints], w->points[i], v2);
715                 if (VectorLength(v1) > 0.5 && VectorLength(v2) > 0.5) break;
716         }
717         CrossProduct(v2, v1, normal);
718         VectorNormalize(normal, normal);
719         *dist = DotProduct(w->points[0], normal);
720 }
721
722 /*
723 =============
724 Winding_Area
725 =============
726 */
727 float Winding_Area (winding_t *w)
728 {
729         int             i;
730         vec3_t  d1, d2, cross;
731         float   total;
732
733         total = 0;
734         for (i=2 ; i<w->numpoints ; i++)
735         {
736                 VectorSubtract (w->points[i-1], w->points[0], d1);
737                 VectorSubtract (w->points[i], w->points[0], d2);
738                 CrossProduct (d1, d2, cross);
739                 total += 0.5 * VectorLength ( cross );
740         }
741         return total;
742 }
743
744 /*
745 =============
746 Winding_Bounds
747 =============
748 */
749 void Winding_Bounds (winding_t *w, vec3_t mins, vec3_t maxs)
750 {
751         vec_t   v;
752         int             i,j;
753
754         mins[0] = mins[1] = mins[2] = 99999;
755         maxs[0] = maxs[1] = maxs[2] = -99999;
756
757         for (i=0 ; i<w->numpoints ; i++)
758         {
759                 for (j=0 ; j<3 ; j++)
760                 {
761                         v = w->points[i][j];
762                         if (v < mins[j])
763                                 mins[j] = v;
764                         if (v > maxs[j])
765                                 maxs[j] = v;
766                 }
767         }
768 }
769
770
771 /*
772 =================
773 Winding_PointInside
774 =================
775 */
776 int Winding_PointInside(winding_t *w, plane_t *plane, vec3_t point, float epsilon)
777 {
778         int i;
779         vec3_t dir, normal, pointvec;
780
781         for (i = 0; i < w->numpoints; i++)
782         {
783                 VectorSubtract(w->points[(i+1) % w->numpoints], w->points[i], dir);
784                 VectorSubtract(point, w->points[i], pointvec);
785                 //
786                 CrossProduct(dir, plane->normal, normal);
787                 //
788                 if (DotProduct(pointvec, normal) < -epsilon) return false;
789         }
790         return true;
791 }
792
793 /*
794 =================
795 Winding_VectorIntersect
796 =================
797 */
798 int Winding_VectorIntersect(winding_t *w, plane_t *plane, vec3_t p1, vec3_t p2, float epsilon)
799 {
800         float front, back, frac;
801         vec3_t mid;
802
803         front = DotProduct(p1, plane->normal) - plane->dist;
804         back = DotProduct(p2, plane->normal) - plane->dist;
805         //if both points at the same side of the plane
806         if (front < -epsilon && back < -epsilon) return false;
807         if (front > epsilon && back > epsilon) return false;
808         //get point of intersection with winding plane
809         if (fabs(front-back) < 0.001)
810         {
811                 VectorCopy(p2, mid);
812         }
813         else
814         {
815                 frac = front/(front-back);
816                 mid[0] = p1[0] + (p2[0] - p1[0]) * frac;
817                 mid[1] = p1[1] + (p2[1] - p1[1]) * frac;
818                 mid[2] = p1[2] + (p2[2] - p1[2]) * frac;
819         }
820         return Winding_PointInside(w, plane, mid, epsilon);
821 }
822