Solvespace

1
//-----------------------------------------------------------------------------
2
// Anything involving surfaces and sets of surfaces (i.e., shells); except
3
// for the real math, which is in ratpoly.cpp.
4
//
5
// Copyright 2008-2013 Jonathan Westhues.
6
//-----------------------------------------------------------------------------
7
#include "../solvespace.h"
8

9
SSurface SSurface::FromExtrusionOf(SBezier *sb, Vector t0, Vector t1) {
10
    SSurface ret = {};
11

12
    ret.degm = sb->deg;
13
    ret.degn = 1;
14

15
    int i;
16
    for(i = 0; i <= ret.degm; i++) {
17
        ret.ctrl[i][0] = (sb->ctrl[i]).Plus(t0);
18
        ret.weight[i][0] = sb->weight[i];
19

20
        ret.ctrl[i][1] = (sb->ctrl[i]).Plus(t1);
21
        ret.weight[i][1] = sb->weight[i];
22
    }
23

24
    return ret;
25
}
26

27
bool SSurface::IsExtrusion(SBezier *of, Vector *alongp) const {
28
    int i;
29

30
    if(degn != 1) return false;
31

32
    Vector along = (ctrl[0][1]).Minus(ctrl[0][0]);
33
    for(i = 0; i <= degm; i++) {
34
        if((fabs(weight[i][1] - weight[i][0]) < LENGTH_EPS) &&
35
           ((ctrl[i][1]).Minus(ctrl[i][0])).Equals(along))
36
        {
37
            continue;
38
        }
39
        return false;
40
    }
41

42
    // yes, we are a surface of extrusion; copy the original curve and return
43
    if(of) {
44
        for(i = 0; i <= degm; i++) {
45
            of->weight[i] = weight[i][0];
46
            of->ctrl[i] = ctrl[i][0];
47
        }
48
        of->deg = degm;
49
        *alongp = along;
50
    }
51
    return true;
52
}
53

54
bool SSurface::IsCylinder(Vector *axis, Vector *center, double *r,
55
                            Vector *start, Vector *finish) const
56
{
57
    SBezier sb;
58
    if(!IsExtrusion(&sb, axis)) return false;
59
    if(!sb.IsCircle(*axis, center, r)) return false;
60

61
    *start = sb.ctrl[0];
62
    *finish = sb.ctrl[2];
63
    return true;
64
}
65

66
// Create a surface patch by revolving and possibly translating a curve.
67
// Works for sections up to but not including 180 degrees.
68
SSurface SSurface::FromRevolutionOf(SBezier *sb, Vector pt, Vector axis, double thetas,
69
                                    double thetaf, double dists,
70
                                    double distf) { // s is start, f is finish
71
    SSurface ret = {};
72
    ret.degm = sb->deg;
73
    ret.degn = 2;
74

75
    double dtheta = fabs(WRAP_SYMMETRIC(thetaf - thetas, 2*PI));
76
    double w      = cos(dtheta / 2);
77

78
    // Revolve the curve about the z axis
79
    int i;
80
    for(i = 0; i <= ret.degm; i++) {
81
        Vector p = sb->ctrl[i];
82

83
        Vector ps = p.RotatedAbout(pt, axis, thetas),
84
               pf = p.RotatedAbout(pt, axis, thetaf);
85

86
        // The middle control point should be at the intersection of the tangents at ps and pf.
87
        // This is equivalent but works for 0 <= angle < 180 degrees.
88
        Vector mid = ps.Plus(pf).ScaledBy(0.5);
89
        Vector c   = ps.ClosestPointOnLine(pt, axis);
90
        Vector ct  = mid.Minus(c).ScaledBy(1 / (w * w)).Plus(c);
91

92
        // not sure this is needed
93
        if(ps.Equals(pf)) {
94
            ps = c;
95
            ct = c;
96
            pf = c;
97
        }
98
        // moving along the axis can create hilical surfaces (or straight extrusion if
99
        // thetas==thetaf)
100
        ret.ctrl[i][0] = ps.Plus(axis.ScaledBy(dists));
101
        ret.ctrl[i][1] = ct.Plus(axis.ScaledBy((dists + distf) / 2));
102
        ret.ctrl[i][2] = pf.Plus(axis.ScaledBy(distf));
103

104
        ret.weight[i][0] = sb->weight[i];
105
        ret.weight[i][1] = sb->weight[i] * w;
106
        ret.weight[i][2] = sb->weight[i];
107
    }
108

109
    return ret;
110
}
111

112
SSurface SSurface::FromPlane(Vector pt, Vector u, Vector v) {
113
    SSurface ret = {};
114

115
    ret.degm = 1;
116
    ret.degn = 1;
117

118
    ret.weight[0][0] = ret.weight[0][1] = 1;
119
    ret.weight[1][0] = ret.weight[1][1] = 1;
120

121
    ret.ctrl[0][0] = pt;
122
    ret.ctrl[0][1] = pt.Plus(u);
123
    ret.ctrl[1][0] = pt.Plus(v);
124
    ret.ctrl[1][1] = pt.Plus(v).Plus(u);
125

126
    return ret;
127
}
128

129
SSurface SSurface::FromTransformationOf(SSurface *a, Vector t, Quaternion q, double scale,
130
                                        bool includingTrims)
131
{
132
    bool needRotate    = !EXACT(q.vx == 0.0 && q.vy == 0.0 && q.vz == 0.0 && q.w == 1.0);
133
    bool needTranslate = !EXACT(t.x  == 0.0 && t.y  == 0.0 && t.z  == 0.0);
134
    bool needScale     = !EXACT(scale == 1.0);
135

136
    SSurface ret = {};
137
    ret.h = a->h;
138
    ret.color = a->color;
139
    ret.face = a->face;
140

141
    ret.degm = a->degm;
142
    ret.degn = a->degn;
143
    int i, j;
144
    for(i = 0; i <= 3; i++) {
145
        for(j = 0; j <= 3; j++) {
146
            Vector ctrl = a->ctrl[i][j];
147
            if(needScale) {
148
                ctrl = ctrl.ScaledBy(scale);
149
            }
150
            if(needRotate) {
151
                ctrl = q.Rotate(ctrl);
152
            }
153
            if(needTranslate) {
154
                ctrl = ctrl.Plus(t);
155
            }
156
            ret.ctrl[i][j] = ctrl;
157
            ret.weight[i][j] = a->weight[i][j];
158
        }
159
    }
160

161
    if(includingTrims) {
162
        STrimBy *stb;
163
        ret.trim.ReserveMore(a->trim.n);
164
        for(stb = a->trim.First(); stb; stb = a->trim.NextAfter(stb)) {
165
            STrimBy n = *stb;
166
            if(needScale) {
167
                n.start  = n.start.ScaledBy(scale);
168
                n.finish = n.finish.ScaledBy(scale);
169
            }
170
            if(needRotate) {
171
                n.start  = q.Rotate(n.start);
172
                n.finish = q.Rotate(n.finish);
173
            }
174
            if(needTranslate) {
175
                n.start  = n.start.Plus(t);
176
                n.finish = n.finish.Plus(t);
177
            }
178
            ret.trim.Add(&n);
179
        }
180
    }
181

182
    if(scale < 0) {
183
        // If we mirror every surface of a shell, then it will end up inside
184
        // out. So fix that here.
185
        ret.Reverse();
186
    }
187

188
    return ret;
189
}
190

191
void SSurface::GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin) const {
192
    *ptMax = Vector::From(VERY_NEGATIVE, VERY_NEGATIVE, VERY_NEGATIVE);
193
    *ptMin = Vector::From(VERY_POSITIVE, VERY_POSITIVE, VERY_POSITIVE);
194

195
    int i, j;
196
    for(i = 0; i <= degm; i++) {
197
        for(j = 0; j <= degn; j++) {
198
            (ctrl[i][j]).MakeMaxMin(ptMax, ptMin);
199
        }
200
    }
201
}
202

203
bool SSurface::LineEntirelyOutsideBbox(Vector a, Vector b, bool asSegment) const {
204
    Vector amax, amin;
205
    GetAxisAlignedBounding(&amax, &amin);
206
    if(!Vector::BoundingBoxIntersectsLine(amax, amin, a, b, asSegment)) {
207
        // The line segment could fail to intersect the bbox, but lie entirely
208
        // within it and intersect the surface.
209
        if(a.OutsideAndNotOn(amax, amin) && b.OutsideAndNotOn(amax, amin)) {
210
            return true;
211
        }
212
    }
213
    return false;
214
}
215

216
//-----------------------------------------------------------------------------
217
// Generate the piecewise linear approximation of the trim stb, which applies
218
// to the curve sc.
219
//-----------------------------------------------------------------------------
220
void SSurface::MakeTrimEdgesInto(SEdgeList *sel, MakeAs flags,
221
                                 SCurve *sc, STrimBy *stb)
222
{
223
    Vector prev = Vector::From(0, 0, 0);
224
    bool inCurve = false, empty = true;
225
    double u = 0, v = 0;
226

227
    int i, first, last, increment;
228
    if(stb->backwards) {
229
        first = sc->pts.n - 1;
230
        last = 0;
231
        increment = -1;
232
    } else {
233
        first = 0;
234
        last = sc->pts.n - 1;
235
        increment = 1;
236
    }
237
    for(i = first; i != (last + increment); i += increment) {
238
        Vector tpt, *pt = &(sc->pts[i].p);
239

240
        if(flags == MakeAs::UV) {
241
            ClosestPointTo(*pt, &u, &v);
242
            tpt = Vector::From(u, v, 0);
243
        } else {
244
            tpt = *pt;
245
        }
246

247
        if(inCurve) {
248
            sel->AddEdge(prev, tpt, sc->h.v, stb->backwards);
249
            empty = false;
250
        }
251

252
        prev = tpt;     // either uv or xyz, depending on flags
253

254
        if(pt->Equals(stb->start)) inCurve = true;
255
        if(pt->Equals(stb->finish)) inCurve = false;
256
    }
257
    if(inCurve) dbp("trim was unterminated");
258
    if(empty)   dbp("trim was empty");
259
}
260

261
//-----------------------------------------------------------------------------
262
// Generate all of our trim curves, in piecewise linear form. We can do
263
// so in either uv or xyz coordinates. And if requested, then we can use
264
// the split curves from useCurvesFrom instead of the curves in our own
265
// shell.
266
//-----------------------------------------------------------------------------
267
void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, MakeAs flags,
268
                             SShell *useCurvesFrom)
269
{
270
    STrimBy *stb;
271
    for(stb = trim.First(); stb; stb = trim.NextAfter(stb)) {
272
        SCurve *sc = shell->curve.FindById(stb->curve);
273

274
        // We have the option to use the curves from another shell; this
275
        // is relevant when generating the coincident edges while doing the
276
        // Booleans, since the curves from the output shell will be split
277
        // against any intersecting surfaces (and the originals aren't).
278
        if(useCurvesFrom) {
279
            sc = useCurvesFrom->curve.FindById(sc->newH);
280
        }
281

282
        MakeTrimEdgesInto(sel, flags, sc, stb);
283
    }
284
}
285

286
//-----------------------------------------------------------------------------
287
// Compute the exact tangent to the intersection curve between two surfaces,
288
// by taking the cross product of the surface normals. We choose the direction
289
// of this tangent so that its dot product with dir is positive.
290
//-----------------------------------------------------------------------------
291
Vector SSurface::ExactSurfaceTangentAt(Vector p, SSurface *srfA, SSurface *srfB, Vector dir)
292
{
293
    Point2d puva, puvb;
294
    srfA->ClosestPointTo(p, &puva);
295
    srfB->ClosestPointTo(p, &puvb);
296
    Vector ts = (srfA->NormalAt(puva)).Cross(
297
                (srfB->NormalAt(puvb)));
298
    ts = ts.WithMagnitude(1);
299
    if(ts.Dot(dir) < 0) {
300
        ts = ts.ScaledBy(-1);
301
    }
302
    return ts;
303
}
304

305
//-----------------------------------------------------------------------------
306
// Report our trim curves. If a trim curve is exact and sbl is not null, then
307
// add its exact form to sbl. Otherwise, add its piecewise linearization to
308
// sel.
309
//-----------------------------------------------------------------------------
310
void SSurface::MakeSectionEdgesInto(SShell *shell, SEdgeList *sel, SBezierList *sbl)
311
{
312
    STrimBy *stb;
313
    for(stb = trim.First(); stb; stb = trim.NextAfter(stb)) {
314
        SCurve *sc = shell->curve.FindById(stb->curve);
315
        SBezier *sb = &(sc->exact);
316

317
        if(sbl && sc->isExact && (sb->deg != 1 || !sel)) {
318
            double ts, tf;
319
            if(stb->backwards) {
320
                sb->ClosestPointTo(stb->start,  &tf);
321
                sb->ClosestPointTo(stb->finish, &ts);
322
            } else {
323
                sb->ClosestPointTo(stb->start,  &ts);
324
                sb->ClosestPointTo(stb->finish, &tf);
325
            }
326
            SBezier junk_bef, keep_aft;
327
            sb->SplitAt(ts, &junk_bef, &keep_aft);
328
            // In the kept piece, the range that used to go from ts to 1
329
            // now goes from 0 to 1; so rescale tf appropriately.
330
            tf = (tf - ts)/(1 - ts);
331

332
            SBezier keep_bef, junk_aft;
333
            keep_aft.SplitAt(tf, &keep_bef, &junk_aft);
334

335
            sbl->l.Add(&keep_bef);
336
        } else if(sbl && !sel && !sc->isExact) {
337
            // We must approximate this trim curve, as piecewise cubic sections.
338
            SSurface *srfA = shell->surface.FindById(sc->surfA);
339
            SSurface *srfB = shell->surface.FindById(sc->surfB);
340

341
            Vector s = stb->backwards ? stb->finish : stb->start,
342
                   f = stb->backwards ? stb->start : stb->finish;
343

344
            int sp, fp;
345
            for(sp = 0; sp < sc->pts.n; sp++) {
346
                if(s.Equals(sc->pts[sp].p)) break;
347
            }
348
            if(sp >= sc->pts.n) return;
349
            for(fp = sp; fp < sc->pts.n; fp++) {
350
                if(f.Equals(sc->pts[fp].p)) break;
351
            }
352
            if(fp >= sc->pts.n) return;
353
            // So now the curve we want goes from elem[sp] to elem[fp]
354

355
            while(sp < fp) {
356
                // Initially, we'll try approximating the entire trim curve
357
                // as a single Bezier segment
358
                int fpt = fp;
359

360
                for(;;) {
361
                    // So construct a cubic Bezier with the correct endpoints
362
                    // and tangents for the current span.
363
                    Vector st = sc->pts[sp].p,
364
                           ft = sc->pts[fpt].p,
365
                           sf = ft.Minus(st);
366
                    double m = sf.Magnitude() / 3;
367

368
                    Vector stan = ExactSurfaceTangentAt(st, srfA, srfB, sf),
369
                           ftan = ExactSurfaceTangentAt(ft, srfA, srfB, sf);
370

371
                    SBezier sb = SBezier::From(st,
372
                                               st.Plus (stan.WithMagnitude(m)),
373
                                               ft.Minus(ftan.WithMagnitude(m)),
374
                                               ft);
375

376
                    // And test how much this curve deviates from the
377
                    // intermediate points (if any).
378
                    int i;
379
                    bool tooFar = false;
380
                    for(i = sp + 1; i <= (fpt - 1); i++) {
381
                        Vector p = sc->pts[i].p;
382
                        double t;
383
                        sb.ClosestPointTo(p, &t, /*mustConverge=*/false);
384
                        Vector pp = sb.PointAt(t);
385
                        if((pp.Minus(p)).Magnitude() > SS.ChordTolMm()/2) {
386
                            tooFar = true;
387
                            break;
388
                        }
389
                    }
390

391
                    if(tooFar) {
392
                        // Deviates by too much, so try a shorter span
393
                        fpt--;
394
                        continue;
395
                    } else {
396
                        // Okay, so use this piece and break.
397
                        sbl->l.Add(&sb);
398
                        break;
399
                    }
400
                }
401

402
                // And continue interpolating, starting wherever the curve
403
                // we just generated finishes.
404
                sp = fpt;
405
            }
406
        } else {
407
            if(sel) MakeTrimEdgesInto(sel, MakeAs::XYZ, sc, stb);
408
        }
409
    }
410
}
411

412
void SSurface::TriangulateInto(SShell *shell, SMesh *sm) {
413
    SEdgeList el = {};
414

415
    MakeEdgesInto(shell, &el, MakeAs::UV);
416

417
    SPolygon poly = {};
418
    if(el.AssemblePolygon(&poly, NULL, /*keepDir=*/true)) {
419
        int i, start = sm->l.n;
420
        if(degm == 1 && degn == 1) {
421
            // A surface with curvature along one direction only; so
422
            // choose the triangulation with chords that lie as much
423
            // as possible within the surface. And since the trim curves
424
            // have been pwl'd to within the desired chord tol, that will
425
            // produce a surface good to within roughly that tol.
426
            //
427
            // If this is just a plane (degree (1, 1)) then the triangulation
428
            // code will notice that, and not bother checking chord tols.
429
            poly.UvTriangulateInto(sm, this);
430
        } else {
431
            // A surface with compound curvature. So we must overlay a
432
            // two-dimensional grid, and triangulate around that.
433
            poly.UvGridTriangulateInto(sm, this);
434
        }
435

436
        STriMeta meta = { face, color };
437
        for(i = start; i < sm->l.n; i++) {
438
            STriangle *st = &(sm->l[i]);
439
            st->meta = meta;
440
            st->an = NormalAt(st->a.x, st->a.y);
441
            st->bn = NormalAt(st->b.x, st->b.y);
442
            st->cn = NormalAt(st->c.x, st->c.y);
443
            st->a = PointAt(st->a.x, st->a.y);
444
            st->b = PointAt(st->b.x, st->b.y);
445
            st->c = PointAt(st->c.x, st->c.y);
446
            // Works out that my chosen contour direction is inconsistent with
447
            // the triangle direction, sigh.
448
            st->FlipNormal();
449
        }
450
    } else {
451
        dbp("failed to assemble polygon to trim nurbs surface in uv space");
452
    }
453

454
    el.Clear();
455
    poly.Clear();
456
}
457

458
//-----------------------------------------------------------------------------
459
// Reverse the parametrisation of one of our dimensions, which flips the
460
// normal. We therefore must reverse all our trim curves too. The uv
461
// coordinates change, but trim curves are stored as xyz so nothing happens
462
//-----------------------------------------------------------------------------
463
void SSurface::Reverse() {
464
    int i, j;
465
    for(i = 0; i < (degm+1)/2; i++) {
466
        for(j = 0; j <= degn; j++) {
467
            swap(ctrl[i][j], ctrl[degm-i][j]);
468
            swap(weight[i][j], weight[degm-i][j]);
469
        }
470
    }
471

472
    STrimBy *stb;
473
    for(stb = trim.First(); stb; stb = trim.NextAfter(stb)) {
474
        stb->backwards = !stb->backwards;
475
        swap(stb->start, stb->finish);
476
    }
477
}
478

479
void SSurface::ScaleSelfBy(double s) {
480
    int i, j;
481
    for(i = 0; i <= degm; i++) {
482
        for(j = 0; j <= degn; j++) {
483
            ctrl[i][j] = ctrl[i][j].ScaledBy(s);
484
        }
485
    }
486
}
487

488
void SSurface::Clear() {
489
    trim.Clear();
490
}
491

492

493