Solvespace

1
//-----------------------------------------------------------------------------
2
// Routines for ray-casting: intersecting a line segment or an infinite line
3
// with a surface or shell. Ray-casting against a shell is used for point-in-
4
// shell testing, and the intersection of edge line segments against surfaces
5
// is used to get rough surface-curve intersections, which are later refined
6
// numerically.
7
//
8
// Copyright 2008-2013 Jonathan Westhues.
9
//-----------------------------------------------------------------------------
10
#include "solvespace.h"
11

12
// Dot product tolerance for perpendicular; this is on the direction cosine,
13
// so it's about 0.001 degrees.
14
const double SShell::DOTP_TOL = 1e-5;
15

16
extern int FLAG;
17

18

19
double SSurface::DepartureFromCoplanar() const {
20
    int i, j;
21
    int ia, ja, ib = 0, jb = 0, ic = 0, jc = 0;
22
    double best;
23

24
    // Grab three points to define a plane; first choose (0, 0) arbitrarily.
25
    ia = ja = 0;
26
    // Then the point farthest from pt a.
27
    best = VERY_NEGATIVE;
28
    for(i = 0; i <= degm; i++) {
29
        for(j = 0; j <= degn; j++) {
30
            if(i == ia && j == ja) continue;
31

32
            double dist = (ctrl[i][j]).Minus(ctrl[ia][ja]).Magnitude();
33
            if(dist > best) {
34
                best = dist;
35
                ib = i;
36
                jb = j;
37
            }
38
        }
39
    }
40
    // Then biggest magnitude of ab cross ac.
41
    best = VERY_NEGATIVE;
42
    for(i = 0; i <= degm; i++) {
43
        for(j = 0; j <= degn; j++) {
44
            if(i == ia && j == ja) continue;
45
            if(i == ib && j == jb) continue;
46

47
            double mag =
48
                ((ctrl[ia][ja].Minus(ctrl[ib][jb]))).Cross(
49
                 (ctrl[ia][ja].Minus(ctrl[i ][j ]))).Magnitude();
50
            if(mag > best) {
51
                best = mag;
52
                ic = i;
53
                jc = j;
54
            }
55
        }
56
    }
57

58
    Vector n = ((ctrl[ia][ja].Minus(ctrl[ib][jb]))).Cross(
59
                (ctrl[ia][ja].Minus(ctrl[ic][jc])));
60
    n = n.WithMagnitude(1);
61
    double d = (ctrl[ia][ja]).Dot(n);
62

63
    // Finally, calculate the deviation from each point to the plane.
64
    double farthest = VERY_NEGATIVE;
65
    for(i = 0; i <= degm; i++) {
66
        for(j = 0; j <= degn; j++) {
67
            double dist = fabs(n.Dot(ctrl[i][j]) - d);
68
            if(dist > farthest) {
69
                farthest = dist;
70
            }
71
        }
72
    }
73
    return farthest;
74
}
75

76
void SSurface::WeightControlPoints() {
77
    int i, j;
78
    for(i = 0; i <= degm; i++) {
79
        for(j = 0; j <= degn; j++) {
80
            ctrl[i][j] = (ctrl[i][j]).ScaledBy(weight[i][j]);
81
        }
82
    }
83
}
84
void SSurface::UnWeightControlPoints() {
85
    int i, j;
86
    for(i = 0; i <= degm; i++) {
87
        for(j = 0; j <= degn; j++) {
88
            ctrl[i][j] = (ctrl[i][j]).ScaledBy(1.0/weight[i][j]);
89
        }
90
    }
91
}
92
void SSurface::CopyRowOrCol(bool row, int this_ij, SSurface *src, int src_ij) {
93
    if(row) {
94
        int j;
95
        for(j = 0; j <= degn; j++) {
96
            ctrl  [this_ij][j] = src->ctrl  [src_ij][j];
97
            weight[this_ij][j] = src->weight[src_ij][j];
98
        }
99
    } else {
100
        int i;
101
        for(i = 0; i <= degm; i++) {
102
            ctrl  [i][this_ij] = src->ctrl  [i][src_ij];
103
            weight[i][this_ij] = src->weight[i][src_ij];
104
        }
105
    }
106
}
107
void SSurface::BlendRowOrCol(bool row, int this_ij, SSurface *a, int a_ij,
108
                                                    SSurface *b, int b_ij)
109
{
110
    if(row) {
111
        int j;
112
        for(j = 0; j <= degn; j++) {
113
            Vector c = (a->ctrl  [a_ij][j]).Plus(b->ctrl  [b_ij][j]);
114
            double w = (a->weight[a_ij][j]   +   b->weight[b_ij][j]);
115
            ctrl  [this_ij][j] = c.ScaledBy(0.5);
116
            weight[this_ij][j] = w / 2;
117
        }
118
    } else {
119
        int i;
120
        for(i = 0; i <= degm; i++) {
121
            Vector c = (a->ctrl  [i][a_ij]).Plus(b->ctrl  [i][b_ij]);
122
            double w = (a->weight[i][a_ij]   +   b->weight[i][b_ij]);
123
            ctrl  [i][this_ij] = c.ScaledBy(0.5);
124
            weight[i][this_ij] = w / 2;
125
        }
126
    }
127
}
128
void SSurface::SplitInHalf(bool byU, SSurface *sa, SSurface *sb) {
129
    sa->degm = sb->degm = degm;
130
    sa->degn = sb->degn = degn;
131

132
    // by de Casteljau's algorithm in a projective space; so we must work
133
    // on points (w*x, w*y, w*z, w) so create a temporary copy
134
    SSurface st;
135
    st = *this;
136
    st.WeightControlPoints();
137

138
    switch(byU ? degm : degn) {
139
        case 1:
140
            sa->CopyRowOrCol (byU, 0, &st, 0);
141
            sb->CopyRowOrCol (byU, 1, &st, 1);
142

143
            sa->BlendRowOrCol(byU, 1, &st, 0, &st, 1);
144
            sb->BlendRowOrCol(byU, 0, &st, 0, &st, 1);
145
            break;
146

147
        case 2:
148
            sa->CopyRowOrCol (byU, 0, &st, 0);
149
            sb->CopyRowOrCol (byU, 2, &st, 2);
150

151
            sa->BlendRowOrCol(byU, 1, &st, 0, &st, 1);
152
            sb->BlendRowOrCol(byU, 1, &st, 1, &st, 2);
153

154
            sa->BlendRowOrCol(byU, 2, sa,   1, sb,   1);
155
            sb->BlendRowOrCol(byU, 0, sa,   1, sb,   1);
156
            break;
157

158
        case 3: {
159
            sa->CopyRowOrCol (byU, 0, &st, 0);
160
            sb->CopyRowOrCol (byU, 3, &st, 3);
161

162
            sa->BlendRowOrCol(byU, 1, &st, 0, &st, 1);
163
            sb->BlendRowOrCol(byU, 2, &st, 2, &st, 3);
164
            st. BlendRowOrCol(byU, 0, &st, 1, &st, 2); // use row/col 0 as scratch
165

166
            sa->BlendRowOrCol(byU, 2, sa,   1, &st,  0);
167
            sb->BlendRowOrCol(byU, 1, sb,   2, &st,  0);
168

169
            sa->BlendRowOrCol(byU, 3, sa,   2, sb,   1);
170
            sb->BlendRowOrCol(byU, 0, sa,   2, sb,   1);
171
            break;
172
        }
173

174
        default: ssassert(false, "Unexpected degree of spline");
175
    }
176

177
    sa->UnWeightControlPoints();
178
    sb->UnWeightControlPoints();
179
}
180

181
//-----------------------------------------------------------------------------
182
// Find all points where the indicated finite (if segment) or infinite (if not
183
// segment) line intersects our surface. Report them in uv space in the list.
184
// We first do a bounding box check; if the line doesn't intersect, then we're
185
// done. If it does, then we check how small our surface is. If it's big,
186
// then we subdivide into quarters and recurse. If it's small, then we refine
187
// by Newton's method and record the point.
188
//-----------------------------------------------------------------------------
189
void SSurface::AllPointsIntersectingUntrimmed(Vector a, Vector b,
190
                                              int *cnt, int *level,
191
                                              List<Inter> *l, bool asSegment,
192
                                              SSurface *sorig)
193
{
194
    // Test if the line intersects our axis-aligned bounding box; if no, then
195
    // no possibility of an intersection
196
    if(LineEntirelyOutsideBbox(a, b, asSegment)) return;
197

198
    if(*cnt > 2000) {
199
        dbp("!!! too many subdivisions (level=%d)!", *level);
200
        dbp("degm = %d degn = %d", degm, degn);
201
        return;
202
    }
203
    (*cnt)++;
204

205
    // If we might intersect, and the surface is small, then switch to Newton
206
    // iterations.
207
    if(DepartureFromCoplanar() < 0.2*SS.ChordTolMm()) {
208
        Vector p = (ctrl[0   ][0   ]).Plus(
209
                    ctrl[0   ][degn]).Plus(
210
                    ctrl[degm][0   ]).Plus(
211
                    ctrl[degm][degn]).ScaledBy(0.25);
212
        Inter inter;
213
        sorig->ClosestPointTo(p, &(inter.p.x), &(inter.p.y), /*mustConverge=*/false);
214
        if(sorig->PointIntersectingLine(a, b, &(inter.p.x), &(inter.p.y))) {
215
            Vector p = sorig->PointAt(inter.p.x, inter.p.y);
216
            // Debug check, verify that the point lies in both surfaces
217
            // (which it ought to, since the surfaces should be coincident)
218
            double u, v;
219
            ClosestPointTo(p, &u, &v);
220
            l->Add(&inter);
221
        } else {
222
            // Might not converge if line is almost tangent to surface...
223
        }
224
        return;
225
    }
226

227
    // But the surface is big, so split it, alternating by u and v
228
    SSurface surf0, surf1;
229
    SplitInHalf((*level & 1) == 0, &surf0, &surf1);
230

231
    int nextLevel = (*level) + 1;
232
    (*level) = nextLevel;
233
    surf0.AllPointsIntersectingUntrimmed(a, b, cnt, level, l, asSegment, sorig);
234
    (*level) = nextLevel;
235
    surf1.AllPointsIntersectingUntrimmed(a, b, cnt, level, l, asSegment, sorig);
236
}
237

238
//-----------------------------------------------------------------------------
239
// Find all points where a line through a and b intersects our surface, and
240
// add them to the list. If seg is true then report only intersections that
241
// lie within the finite line segment (not including the endpoints); otherwise
242
// we work along the infinite line. And we report either just intersections
243
// inside the trim curve, or any intersection with u, v in [0, 1]. And we
244
// either disregard or report tangent points.
245
//-----------------------------------------------------------------------------
246
void SSurface::AllPointsIntersecting(Vector a, Vector b,
247
                                     List<SInter> *l,
248
                                     bool asSegment, bool trimmed, bool inclTangent)
249
{
250
    if(LineEntirelyOutsideBbox(a, b, asSegment)) return;
251

252
    Vector ba = b.Minus(a);
253
    double bam = ba.Magnitude();
254

255
    List<Inter> inters = {};
256

257
    // All the intersections between the line and the surface; either special
258
    // cases that we can quickly solve in closed form, or general numerical.
259
    Vector center, axis, start, finish;
260
    double radius;
261
    if(degm == 1 && degn == 1) {
262
        // Against a plane, easy.
263
        Vector n = NormalAt(0, 0).WithMagnitude(1);
264
        double d = n.Dot(PointAt(0, 0));
265
        // Trim to line segment now if requested, don't generate points that
266
        // would just get discarded later.
267
        if(!asSegment ||
268
           (n.Dot(a) > d + LENGTH_EPS && n.Dot(b) < d - LENGTH_EPS) ||
269
           (n.Dot(b) > d + LENGTH_EPS && n.Dot(a) < d - LENGTH_EPS))
270
        {
271
            Vector p = Vector::AtIntersectionOfPlaneAndLine(n, d, a, b, NULL);
272
            Inter inter;
273
            ClosestPointTo(p, &(inter.p.x), &(inter.p.y));
274
            inters.Add(&inter);
275
        }
276
    } else if(IsCylinder(&axis, &center, &radius, &start, &finish)) {
277
        // This one can be solved in closed form too.
278
        Vector ab = b.Minus(a);
279
        if(axis.Cross(ab).Magnitude() < LENGTH_EPS) {
280
            // edge is parallel to axis of cylinder, no intersection points
281
            return;
282
        }
283
        // A coordinate system centered at the center of the circle, with
284
        // the edge under test horizontal
285
        Vector u, v, n = axis.WithMagnitude(1);
286
        u = (ab.Minus(n.ScaledBy(ab.Dot(n)))).WithMagnitude(1);
287
        v = n.Cross(u);
288
        Point2d ap = (a.Minus(center)).DotInToCsys(u, v, n).ProjectXy(),
289
                bp = (b.Minus(center)).DotInToCsys(u, v, n).ProjectXy(),
290
                sp = (start. Minus(center)).DotInToCsys(u, v, n).ProjectXy(),
291
                fp = (finish.Minus(center)).DotInToCsys(u, v, n).ProjectXy();
292

293
        double thetas = atan2(sp.y, sp.x), thetaf = atan2(fp.y, fp.x);
294

295
        Point2d ip[2];
296
        int ip_n = 0;
297
        if(fabs(fabs(ap.y) - radius) < LENGTH_EPS) {
298
            // tangent
299
            if(inclTangent) {
300
                ip[0] = Point2d::From(0, ap.y);
301
                ip_n = 1;
302
            }
303
        } else if(fabs(ap.y) < radius) {
304
            // two intersections
305
            double xint = sqrt(radius*radius - ap.y*ap.y);
306
            ip[0] = Point2d::From(-xint, ap.y);
307
            ip[1] = Point2d::From( xint, ap.y);
308
            ip_n = 2;
309
        }
310
        int i;
311
        for(i = 0; i < ip_n; i++) {
312
            double t = (ip[i].Minus(ap)).DivProjected(bp.Minus(ap));
313
            // This is a point on the circle; but is it on the arc?
314
            Point2d pp = ap.Plus((bp.Minus(ap)).ScaledBy(t));
315
            double theta = atan2(pp.y, pp.x);
316
            double dp = WRAP_SYMMETRIC(theta  - thetas, 2*PI),
317
                   df = WRAP_SYMMETRIC(thetaf - thetas, 2*PI);
318
            double tol = LENGTH_EPS/radius;
319

320
            if((df > 0 && ((dp < -tol) || (dp > df + tol))) ||
321
               (df < 0 && ((dp >  tol) || (dp < df - tol))))
322
            {
323
                continue;
324
            }
325

326
            Vector p = a.Plus((b.Minus(a)).ScaledBy(t));
327

328
            Inter inter;
329
            ClosestPointTo(p, &(inter.p.x), &(inter.p.y));
330
            inters.Add(&inter);
331
        }
332
    } else {
333
        // General numerical solution by subdivision, fallback
334
        int cnt = 0, level = 0;
335
        AllPointsIntersectingUntrimmed(a, b, &cnt, &level, &inters, asSegment, this);
336
    }
337

338
    // Remove duplicate intersection points
339
    inters.ClearTags();
340
    int i, j;
341
    for(i = 0; i < inters.n; i++) {
342
        for(j = i + 1; j < inters.n; j++) {
343
            if(inters[i].p.Equals(inters[j].p)) {
344
                inters[j].tag = 1;
345
            }
346
        }
347
    }
348
    inters.RemoveTagged();
349

350
    for(i = 0; i < inters.n; i++) {
351
        Point2d puv = inters[i].p;
352

353
        // Make sure the point lies within the finite line segment
354
        Vector pxyz = PointAt(puv.x, puv.y);
355
        double t = (pxyz.Minus(a)).DivProjected(ba);
356
        if(asSegment && (t > 1 - LENGTH_EPS/bam || t < LENGTH_EPS/bam)) {
357
            continue;
358
        }
359

360
        // And that it lies inside our trim region
361
        Point2d dummy = { 0, 0 };
362
        SBspUv::Class c = (bsp) ? bsp->ClassifyPoint(puv, dummy, this) : SBspUv::Class::OUTSIDE;
363
        if(trimmed && c == SBspUv::Class::OUTSIDE) {
364
            continue;
365
        }
366

367
        // It does, so generate the intersection
368
        SInter si;
369
        si.p = pxyz;
370
        si.surfNormal = NormalAt(puv.x, puv.y);
371
        si.pinter = puv;
372
        si.srf = this;
373
        si.onEdge = (c != SBspUv::Class::INSIDE);
374
        l->Add(&si);
375
    }
376

377
    inters.Clear();
378
}
379

380
void SShell::AllPointsIntersecting(Vector a, Vector b,
381
                                   List<SInter> *il,
382
                                   bool asSegment, bool trimmed, bool inclTangent)
383
{
384
    for(SSurface &ss : surface) {
385
        ss.AllPointsIntersecting(a, b, il,
386
            asSegment, trimmed, inclTangent);
387
    }
388
}
389

390

391

392
SShell::Class SShell::ClassifyRegion(Vector edge_n, Vector inter_surf_n,
393
                           Vector edge_surf_n) const
394
{
395
    double dot = inter_surf_n.DirectionCosineWith(edge_n);
396
    if(fabs(dot) < DOTP_TOL) {
397
        // The edge's surface and the edge-on-face surface
398
        // are coincident. Test the edge's surface normal
399
        // to see if it's with same or opposite normals.
400
        if(inter_surf_n.Dot(edge_surf_n) > 0) {
401
            return Class::SURF_COINC_SAME;
402
        } else {
403
            return Class::SURF_COINC_OPP;
404
        }
405
    } else if(dot > 0) {
406
        return Class::SURF_OUTSIDE;
407
    } else {
408
        return Class::SURF_INSIDE;
409
    }
410
}
411

412
//-----------------------------------------------------------------------------
413
// Does the given point lie on our shell? There are many cases; inside and
414
// outside are obvious, but then there's all the edge-on-edge and edge-on-face
415
// possibilities.
416
//
417
// To calculate, we intersect a ray through p with our shell, and classify
418
// using the closest intersection point. If the ray hits a surface on edge,
419
// then just reattempt in a different random direction.
420
//-----------------------------------------------------------------------------
421

422
// table of vectors in 6 arbitrary directions covering 4 of the 8 octants.
423
// use overlapping sets of 3 to reduce memory usage.
424
static const double Random[8] = {1.278, 5.0103, 9.427, -2.331, 7.13, 2.954, 5.034, -4.777};
425
 
426
bool SShell::ClassifyEdge(Class *indir, Class *outdir,
427
                          Vector ea, Vector eb,
428
                          Vector p,
429
                          Vector edge_n_in, Vector edge_n_out, Vector surf_n)
430
{
431
    List<SInter> l = {};
432

433
    // First, check for edge-on-edge
434
    int edge_inters = 0;
435
    Vector inter_surf_n[2], inter_edge_n[2];
436
    for(SSurface &srf : surface) {
437
        if(srf.LineEntirelyOutsideBbox(ea, eb, /*asSegment=*/true)) continue;
438

439
        SEdgeList *sel = &(srf.edges);
440
        SEdge *se;
441
        for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {
442
            if((ea.Equals(se->a) && eb.Equals(se->b)) ||
443
               (eb.Equals(se->a) && ea.Equals(se->b)) ||
444
                p.OnLineSegment(se->a, se->b))
445
            {
446
                if(edge_inters < 2) {
447
                    // Edge-on-edge case
448
                    Point2d pm;
449
                    srf.ClosestPointTo(p,  &pm, /*mustConverge=*/false);
450
                    // A vector normal to the surface, at the intersection point
451
                    inter_surf_n[edge_inters] = srf.NormalAt(pm);
452
                    // A vector normal to the intersecting edge (but within the
453
                    // intersecting surface) at the intersection point, pointing
454
                    // out.
455
                    inter_edge_n[edge_inters] =
456
                      (inter_surf_n[edge_inters]).Cross((se->b).Minus((se->a)));
457
                }
458

459
                edge_inters++;
460
            }
461
        }
462
    }
463

464
    if(edge_inters == 2) {
465
        //! @todo make this use the appropriate curved normals
466
        double dotp[2];
467
        for(int i = 0; i < 2; i++) {
468
            dotp[i] = edge_n_out.DirectionCosineWith(inter_surf_n[i]);
469
        }
470

471
        if(fabs(dotp[1]) < DOTP_TOL) {
472
            swap(dotp[0],         dotp[1]);
473
            swap(inter_surf_n[0], inter_surf_n[1]);
474
            swap(inter_edge_n[0], inter_edge_n[1]);
475
        }
476

477
        Class coinc = (surf_n.Dot(inter_surf_n[0])) > 0 ? Class::SURF_COINC_SAME : Class::SURF_COINC_OPP;
478

479
        if(fabs(dotp[0]) < DOTP_TOL && fabs(dotp[1]) < DOTP_TOL) {
480
            // This is actually an edge on face case, just that the face
481
            // is split into two pieces joining at our edge.
482
            *indir  = coinc;
483
            *outdir = coinc;
484
        } else if(fabs(dotp[0]) < DOTP_TOL && dotp[1] > DOTP_TOL) {
485
            if(edge_n_out.Dot(inter_edge_n[0]) > 0) {
486
                *indir  = coinc;
487
                *outdir = Class::SURF_OUTSIDE;
488
            } else {
489
                *indir  = Class::SURF_INSIDE;
490
                *outdir = coinc;
491
            }
492
        } else if(fabs(dotp[0]) < DOTP_TOL && dotp[1] < -DOTP_TOL) {
493
            if(edge_n_out.Dot(inter_edge_n[0]) > 0) {
494
                *indir  = coinc;
495
                *outdir = Class::SURF_INSIDE;
496
            } else {
497
                *indir  = Class::SURF_OUTSIDE;
498
                *outdir = coinc;
499
            }
500
        } else if(dotp[0] > DOTP_TOL && dotp[1] > DOTP_TOL) {
501
            *indir  = Class::SURF_INSIDE;
502
            *outdir = Class::SURF_OUTSIDE;
503
        } else if(dotp[0] < -DOTP_TOL && dotp[1] < -DOTP_TOL) {
504
            *indir  = Class::SURF_OUTSIDE;
505
            *outdir = Class::SURF_INSIDE;
506
        } else {
507
            // Edge is tangent to the shell at shell's edge, so can't be
508
            // a boundary of the surface.
509
            return false;
510
        }
511
        return true;
512
    }
513

514
    if(edge_inters != 0) dbp("bad, edge_inters=%d", edge_inters);
515

516
    // Next, check for edge-on-surface. The ray-casting for edge-inside-shell
517
    // would catch this too, but test separately, for speed (since many edges
518
    // are on surface) and for numerical stability, so we don't pick up
519
    // the additional error from the line intersection.
520

521
    for(SSurface &srf : surface) {
522
        if(srf.LineEntirelyOutsideBbox(ea, eb, /*asSegment=*/true)) continue;
523

524
        Point2d puv;
525
        srf.ClosestPointTo(p, &(puv.x), &(puv.y), /*mustConverge=*/false);
526
        Vector pp = srf.PointAt(puv);
527

528
        if((pp.Minus(p)).Magnitude() > LENGTH_EPS) continue;
529
        Point2d dummy = { 0, 0 };
530
        SBspUv::Class c = (srf.bsp) ? srf.bsp->ClassifyPoint(puv, dummy, &srf) : SBspUv::Class::OUTSIDE;
531
        if(c == SBspUv::Class::OUTSIDE) continue;
532

533
        // Edge-on-face (unless edge-on-edge above superseded)
534
        Point2d pin, pout;
535
        srf.ClosestPointTo(p.Plus(edge_n_in),  &pin,  /*mustConverge=*/false);
536
        srf.ClosestPointTo(p.Plus(edge_n_out), &pout, /*mustConverge=*/false);
537

538
        Vector surf_n_in  = srf.NormalAt(pin),
539
               surf_n_out = srf.NormalAt(pout);
540

541
        *indir  = ClassifyRegion(edge_n_in,  surf_n_in,  surf_n);
542
        *outdir = ClassifyRegion(edge_n_out, surf_n_out, surf_n);
543
        return true;
544
    }
545

546
    // Edge is not on face or on edge; so it's either inside or outside
547
    // the shell, and we'll determine which by raycasting.
548
    int cnt = 0;
549
    for(;;) {
550
        // Cast a ray in a random direction (two-sided so that we test if
551
        // the point lies on a surface, but use only one side for in/out
552
        // testing)
553
        Vector ray = Vector::From(Random[cnt], Random[cnt+1], Random[cnt+2]);
554

555
        AllPointsIntersecting(
556
            p.Minus(ray), p.Plus(ray), &l,
557
                /*asSegment=*/false, /*trimmed=*/true, /*inclTangent=*/false);
558

559
        // no intersections means it's outside
560
        *indir  = Class::SURF_OUTSIDE;
561
        *outdir = Class::SURF_OUTSIDE;
562
        double dmin = VERY_POSITIVE;
563
        bool onEdge = false;
564
        edge_inters = 0;
565

566
        SInter *si;
567
        for(si = l.First(); si; si = l.NextAfter(si)) {
568
            double t = ((si->p).Minus(p)).DivProjected(ray);
569
            if(t*ray.Magnitude() < -LENGTH_EPS) {
570
                // wrong side, doesn't count
571
                continue;
572
            }
573

574
            double d = ((si->p).Minus(p)).Magnitude();
575

576
            // We actually should never hit this case; it should have been
577
            // handled above.
578
            if(d < LENGTH_EPS && si->onEdge) {
579
                edge_inters++;
580
            }
581

582
            if(d < dmin) {
583
                dmin = d;
584
                // Edge does not lie on surface; either strictly inside
585
                // or strictly outside
586
                if((si->surfNormal).Dot(ray) > 0) {
587
                    *indir  = Class::SURF_INSIDE;
588
                    *outdir = Class::SURF_INSIDE;
589
                } else {
590
                    *indir  = Class::SURF_OUTSIDE;
591
                    *outdir = Class::SURF_OUTSIDE;
592
                }
593
                onEdge = si->onEdge;
594
            }
595
        }
596
        l.Clear();
597

598
        // If the point being tested lies exactly on an edge of the shell,
599
        // then our ray always lies on edge, and that's okay. Otherwise
600
        // try again in a different random direction.
601
        if(!onEdge) break;
602
        cnt++;
603
        if(cnt > 5) {
604
            dbp("can't find a ray that doesn't hit on edge!");
605
            dbp("on edge = %d, edge_inters = %d", onEdge, edge_inters);
606
            SS.nakedEdges.AddEdge(ea, eb);
607
            break;
608
        }
609
    }
610

611
    return true;
612
}
613

614