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
FromExtrusionOf(SBezier * sb,Vector t0,Vector t1)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
IsExtrusion(SBezier * of,Vector * alongp)27 bool SSurface::IsExtrusion(SBezier *of, Vector *alongp) {
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
IsCylinder(Vector * axis,Vector * center,double * r,Vector * start,Vector * finish)54 bool SSurface::IsCylinder(Vector *axis, Vector *center, double *r,
55 Vector *start, Vector *finish)
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
FromRevolutionOf(SBezier * sb,Vector pt,Vector axis,double thetas,double thetaf)66 SSurface SSurface::FromRevolutionOf(SBezier *sb, Vector pt, Vector axis,
67 double thetas, double thetaf)
68 {
69 SSurface ret = {};
70
71
72 ret.degm = sb->deg;
73 ret.degn = 2;
74
75 double dtheta = fabs(WRAP_SYMMETRIC(thetaf - thetas, 2*PI));
76
77 // We now wish to revolve the curve about the z axis
78 int i;
79 for(i = 0; i <= ret.degm; i++) {
80 Vector p = sb->ctrl[i];
81
82 Vector ps = p.RotatedAbout(pt, axis, thetas),
83 pf = p.RotatedAbout(pt, axis, thetaf);
84
85 Vector ct;
86 if(ps.Equals(pf)) {
87 // Degenerate case: a control point lies on the axis of revolution,
88 // so we get three coincident control points.
89 ct = ps;
90 } else {
91 // Normal case, the control point sweeps out a circle.
92 Vector c = ps.ClosestPointOnLine(pt, axis);
93
94 Vector rs = ps.Minus(c),
95 rf = pf.Minus(c);
96
97 Vector ts = axis.Cross(rs),
98 tf = axis.Cross(rf);
99
100 ct = Vector::AtIntersectionOfLines(ps, ps.Plus(ts),
101 pf, pf.Plus(tf),
102 NULL, NULL, NULL);
103 }
104
105 ret.ctrl[i][0] = ps;
106 ret.ctrl[i][1] = ct;
107 ret.ctrl[i][2] = pf;
108
109 ret.weight[i][0] = sb->weight[i];
110 ret.weight[i][1] = sb->weight[i]*cos(dtheta/2);
111 ret.weight[i][2] = sb->weight[i];
112 }
113
114 return ret;
115 }
116
FromPlane(Vector pt,Vector u,Vector v)117 SSurface SSurface::FromPlane(Vector pt, Vector u, Vector v) {
118 SSurface ret = {};
119
120 ret.degm = 1;
121 ret.degn = 1;
122
123 ret.weight[0][0] = ret.weight[0][1] = 1;
124 ret.weight[1][0] = ret.weight[1][1] = 1;
125
126 ret.ctrl[0][0] = pt;
127 ret.ctrl[0][1] = pt.Plus(u);
128 ret.ctrl[1][0] = pt.Plus(v);
129 ret.ctrl[1][1] = pt.Plus(v).Plus(u);
130
131 return ret;
132 }
133
FromTransformationOf(SSurface * a,Vector t,Quaternion q,double scale,bool includingTrims)134 SSurface SSurface::FromTransformationOf(SSurface *a,
135 Vector t, Quaternion q, double scale,
136 bool includingTrims)
137 {
138 SSurface ret = {};
139
140 ret.h = a->h;
141 ret.color = a->color;
142 ret.face = a->face;
143
144 ret.degm = a->degm;
145 ret.degn = a->degn;
146 int i, j;
147 for(i = 0; i <= 3; i++) {
148 for(j = 0; j <= 3; j++) {
149 ret.ctrl[i][j] = a->ctrl[i][j];
150 ret.ctrl[i][j] = (ret.ctrl[i][j]).ScaledBy(scale);
151 ret.ctrl[i][j] = (q.Rotate(ret.ctrl[i][j])).Plus(t);
152
153 ret.weight[i][j] = a->weight[i][j];
154 }
155 }
156
157 if(includingTrims) {
158 STrimBy *stb;
159 for(stb = a->trim.First(); stb; stb = a->trim.NextAfter(stb)) {
160 STrimBy n = *stb;
161 n.start = n.start.ScaledBy(scale);
162 n.finish = n.finish.ScaledBy(scale);
163 n.start = (q.Rotate(n.start)) .Plus(t);
164 n.finish = (q.Rotate(n.finish)).Plus(t);
165 ret.trim.Add(&n);
166 }
167 }
168
169 if(scale < 0) {
170 // If we mirror every surface of a shell, then it will end up inside
171 // out. So fix that here.
172 ret.Reverse();
173 }
174
175 return ret;
176 }
177
GetAxisAlignedBounding(Vector * ptMax,Vector * ptMin)178 void SSurface::GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin) {
179 *ptMax = Vector::From(VERY_NEGATIVE, VERY_NEGATIVE, VERY_NEGATIVE);
180 *ptMin = Vector::From(VERY_POSITIVE, VERY_POSITIVE, VERY_POSITIVE);
181
182 int i, j;
183 for(i = 0; i <= degm; i++) {
184 for(j = 0; j <= degn; j++) {
185 (ctrl[i][j]).MakeMaxMin(ptMax, ptMin);
186 }
187 }
188 }
189
LineEntirelyOutsideBbox(Vector a,Vector b,bool segment)190 bool SSurface::LineEntirelyOutsideBbox(Vector a, Vector b, bool segment) {
191 Vector amax, amin;
192 GetAxisAlignedBounding(&amax, &amin);
193 if(!Vector::BoundingBoxIntersectsLine(amax, amin, a, b, segment)) {
194 // The line segment could fail to intersect the bbox, but lie entirely
195 // within it and intersect the surface.
196 if(a.OutsideAndNotOn(amax, amin) && b.OutsideAndNotOn(amax, amin)) {
197 return true;
198 }
199 }
200 return false;
201 }
202
203 //-----------------------------------------------------------------------------
204 // Generate the piecewise linear approximation of the trim stb, which applies
205 // to the curve sc.
206 //-----------------------------------------------------------------------------
MakeTrimEdgesInto(SEdgeList * sel,int flags,SCurve * sc,STrimBy * stb)207 void SSurface::MakeTrimEdgesInto(SEdgeList *sel, int flags,
208 SCurve *sc, STrimBy *stb)
209 {
210 Vector prev = Vector::From(0, 0, 0);
211 bool inCurve = false, empty = true;
212 double u = 0, v = 0;
213
214 int i, first, last, increment;
215 if(stb->backwards) {
216 first = sc->pts.n - 1;
217 last = 0;
218 increment = -1;
219 } else {
220 first = 0;
221 last = sc->pts.n - 1;
222 increment = 1;
223 }
224 for(i = first; i != (last + increment); i += increment) {
225 Vector tpt, *pt = &(sc->pts.elem[i].p);
226
227 if(flags & AS_UV) {
228 ClosestPointTo(*pt, &u, &v);
229 tpt = Vector::From(u, v, 0);
230 } else {
231 tpt = *pt;
232 }
233
234 if(inCurve) {
235 sel->AddEdge(prev, tpt, sc->h.v, stb->backwards);
236 empty = false;
237 }
238
239 prev = tpt; // either uv or xyz, depending on flags
240
241 if(pt->Equals(stb->start)) inCurve = true;
242 if(pt->Equals(stb->finish)) inCurve = false;
243 }
244 if(inCurve) dbp("trim was unterminated");
245 if(empty) dbp("trim was empty");
246 }
247
248 //-----------------------------------------------------------------------------
249 // Generate all of our trim curves, in piecewise linear form. We can do
250 // so in either uv or xyz coordinates. And if requested, then we can use
251 // the split curves from useCurvesFrom instead of the curves in our own
252 // shell.
253 //-----------------------------------------------------------------------------
MakeEdgesInto(SShell * shell,SEdgeList * sel,int flags,SShell * useCurvesFrom)254 void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, int flags,
255 SShell *useCurvesFrom)
256 {
257 STrimBy *stb;
258 for(stb = trim.First(); stb; stb = trim.NextAfter(stb)) {
259 SCurve *sc = shell->curve.FindById(stb->curve);
260
261 // We have the option to use the curves from another shell; this
262 // is relevant when generating the coincident edges while doing the
263 // Booleans, since the curves from the output shell will be split
264 // against any intersecting surfaces (and the originals aren't).
265 if(useCurvesFrom) {
266 sc = useCurvesFrom->curve.FindById(sc->newH);
267 }
268
269 MakeTrimEdgesInto(sel, flags, sc, stb);
270 }
271 }
272
273 //-----------------------------------------------------------------------------
274 // Compute the exact tangent to the intersection curve between two surfaces,
275 // by taking the cross product of the surface normals. We choose the direction
276 // of this tangent so that its dot product with dir is positive.
277 //-----------------------------------------------------------------------------
ExactSurfaceTangentAt(Vector p,SSurface * srfA,SSurface * srfB,Vector dir)278 Vector SSurface::ExactSurfaceTangentAt(Vector p, SSurface *srfA, SSurface *srfB,
279 Vector dir)
280 {
281 Point2d puva, puvb;
282 srfA->ClosestPointTo(p, &puva);
283 srfB->ClosestPointTo(p, &puvb);
284 Vector ts = (srfA->NormalAt(puva)).Cross(
285 (srfB->NormalAt(puvb)));
286 ts = ts.WithMagnitude(1);
287 if(ts.Dot(dir) < 0) {
288 ts = ts.ScaledBy(-1);
289 }
290 return ts;
291 }
292
293 //-----------------------------------------------------------------------------
294 // Report our trim curves. If a trim curve is exact and sbl is not null, then
295 // add its exact form to sbl. Otherwise, add its piecewise linearization to
296 // sel.
297 //-----------------------------------------------------------------------------
MakeSectionEdgesInto(SShell * shell,SEdgeList * sel,SBezierList * sbl)298 void SSurface::MakeSectionEdgesInto(SShell *shell,
299 SEdgeList *sel, SBezierList *sbl)
300 {
301 STrimBy *stb;
302 for(stb = trim.First(); stb; stb = trim.NextAfter(stb)) {
303 SCurve *sc = shell->curve.FindById(stb->curve);
304 SBezier *sb = &(sc->exact);
305
306 if(sbl && sc->isExact && (sb->deg != 1 || !sel)) {
307 double ts, tf;
308 if(stb->backwards) {
309 sb->ClosestPointTo(stb->start, &tf);
310 sb->ClosestPointTo(stb->finish, &ts);
311 } else {
312 sb->ClosestPointTo(stb->start, &ts);
313 sb->ClosestPointTo(stb->finish, &tf);
314 }
315 SBezier junk_bef, keep_aft;
316 sb->SplitAt(ts, &junk_bef, &keep_aft);
317 // In the kept piece, the range that used to go from ts to 1
318 // now goes from 0 to 1; so rescale tf appropriately.
319 tf = (tf - ts)/(1 - ts);
320
321 SBezier keep_bef, junk_aft;
322 keep_aft.SplitAt(tf, &keep_bef, &junk_aft);
323
324 sbl->l.Add(&keep_bef);
325 } else if(sbl && !sel && !sc->isExact) {
326 // We must approximate this trim curve, as piecewise cubic sections.
327 SSurface *srfA = shell->surface.FindById(sc->surfA),
328 *srfB = shell->surface.FindById(sc->surfB);
329
330 Vector s = stb->backwards ? stb->finish : stb->start,
331 f = stb->backwards ? stb->start : stb->finish;
332
333 int sp, fp;
334 for(sp = 0; sp < sc->pts.n; sp++) {
335 if(s.Equals(sc->pts.elem[sp].p)) break;
336 }
337 if(sp >= sc->pts.n) return;
338 for(fp = sp; fp < sc->pts.n; fp++) {
339 if(f.Equals(sc->pts.elem[fp].p)) break;
340 }
341 if(fp >= sc->pts.n) return;
342 // So now the curve we want goes from elem[sp] to elem[fp]
343
344 while(sp < fp) {
345 // Initially, we'll try approximating the entire trim curve
346 // as a single Bezier segment
347 int fpt = fp;
348
349 for(;;) {
350 // So construct a cubic Bezier with the correct endpoints
351 // and tangents for the current span.
352 Vector st = sc->pts.elem[sp].p,
353 ft = sc->pts.elem[fpt].p,
354 sf = ft.Minus(st);
355 double m = sf.Magnitude() / 3;
356
357 Vector stan = ExactSurfaceTangentAt(st, srfA, srfB, sf),
358 ftan = ExactSurfaceTangentAt(ft, srfA, srfB, sf);
359
360 SBezier sb = SBezier::From(st,
361 st.Plus (stan.WithMagnitude(m)),
362 ft.Minus(ftan.WithMagnitude(m)),
363 ft);
364
365 // And test how much this curve deviates from the
366 // intermediate points (if any).
367 int i;
368 bool tooFar = false;
369 for(i = sp + 1; i <= (fpt - 1); i++) {
370 Vector p = sc->pts.elem[i].p;
371 double t;
372 sb.ClosestPointTo(p, &t, false);
373 Vector pp = sb.PointAt(t);
374 if((pp.Minus(p)).Magnitude() > SS.ChordTolMm()/2) {
375 tooFar = true;
376 break;
377 }
378 }
379
380 if(tooFar) {
381 // Deviates by too much, so try a shorter span
382 fpt--;
383 continue;
384 } else {
385 // Okay, so use this piece and break.
386 sbl->l.Add(&sb);
387 break;
388 }
389 }
390
391 // And continue interpolating, starting wherever the curve
392 // we just generated finishes.
393 sp = fpt;
394 }
395 } else {
396 if(sel) MakeTrimEdgesInto(sel, AS_XYZ, sc, stb);
397 }
398 }
399 }
400
TriangulateInto(SShell * shell,SMesh * sm)401 void SSurface::TriangulateInto(SShell *shell, SMesh *sm) {
402 SEdgeList el = {};
403
404 MakeEdgesInto(shell, &el, AS_UV);
405
406 SPolygon poly = {};
407 if(el.AssemblePolygon(&poly, NULL, true)) {
408 int i, start = sm->l.n;
409 if(degm == 1 && degn == 1) {
410 // A surface with curvature along one direction only; so
411 // choose the triangulation with chords that lie as much
412 // as possible within the surface. And since the trim curves
413 // have been pwl'd to within the desired chord tol, that will
414 // produce a surface good to within roughly that tol.
415 //
416 // If this is just a plane (degree (1, 1)) then the triangulation
417 // code will notice that, and not bother checking chord tols.
418 poly.UvTriangulateInto(sm, this);
419 } else {
420 // A surface with compound curvature. So we must overlay a
421 // two-dimensional grid, and triangulate around that.
422 poly.UvGridTriangulateInto(sm, this);
423 }
424
425 STriMeta meta = { face, color };
426 for(i = start; i < sm->l.n; i++) {
427 STriangle *st = &(sm->l.elem[i]);
428 st->meta = meta;
429 if(st->meta.color.alpha != 255) sm->isTransparent = true;
430 st->an = NormalAt(st->a.x, st->a.y);
431 st->bn = NormalAt(st->b.x, st->b.y);
432 st->cn = NormalAt(st->c.x, st->c.y);
433 st->a = PointAt(st->a.x, st->a.y);
434 st->b = PointAt(st->b.x, st->b.y);
435 st->c = PointAt(st->c.x, st->c.y);
436 // Works out that my chosen contour direction is inconsistent with
437 // the triangle direction, sigh.
438 st->FlipNormal();
439 }
440 } else {
441 dbp("failed to assemble polygon to trim nurbs surface in uv space");
442 }
443
444 el.Clear();
445 poly.Clear();
446 }
447
448 //-----------------------------------------------------------------------------
449 // Reverse the parametrisation of one of our dimensions, which flips the
450 // normal. We therefore must reverse all our trim curves too. The uv
451 // coordinates change, but trim curves are stored as xyz so nothing happens
452 //-----------------------------------------------------------------------------
Reverse(void)453 void SSurface::Reverse(void) {
454 int i, j;
455 for(i = 0; i < (degm+1)/2; i++) {
456 for(j = 0; j <= degn; j++) {
457 swap(ctrl[i][j], ctrl[degm-i][j]);
458 swap(weight[i][j], weight[degm-i][j]);
459 }
460 }
461
462 STrimBy *stb;
463 for(stb = trim.First(); stb; stb = trim.NextAfter(stb)) {
464 stb->backwards = !stb->backwards;
465 swap(stb->start, stb->finish);
466 }
467 }
468
ScaleSelfBy(double s)469 void SSurface::ScaleSelfBy(double s) {
470 int i, j;
471 for(i = 0; i <= degm; i++) {
472 for(j = 0; j <= degn; j++) {
473 ctrl[i][j] = ctrl[i][j].ScaledBy(s);
474 }
475 }
476 }
477
Clear(void)478 void SSurface::Clear(void) {
479 trim.Clear();
480 }
481
482 typedef struct {
483 hSCurve hc;
484 hSSurface hs;
485 } TrimLine;
486
MakeFromExtrusionOf(SBezierLoopSet * sbls,Vector t0,Vector t1,RgbaColor color)487 void SShell::MakeFromExtrusionOf(SBezierLoopSet *sbls, Vector t0, Vector t1, RgbaColor color)
488 {
489 // Make the extrusion direction consistent with respect to the normal
490 // of the sketch we're extruding.
491 if((t0.Minus(t1)).Dot(sbls->normal) < 0) {
492 swap(t0, t1);
493 }
494
495 // Define a coordinate system to contain the original sketch, and get
496 // a bounding box in that csys
497 Vector n = sbls->normal.ScaledBy(-1);
498 Vector u = n.Normal(0), v = n.Normal(1);
499 Vector orig = sbls->point;
500 double umax = 1e-10, umin = 1e10;
501 sbls->GetBoundingProjd(u, orig, &umin, &umax);
502 double vmax = 1e-10, vmin = 1e10;
503 sbls->GetBoundingProjd(v, orig, &vmin, &vmax);
504 // and now fix things up so that all u and v lie between 0 and 1
505 orig = orig.Plus(u.ScaledBy(umin));
506 orig = orig.Plus(v.ScaledBy(vmin));
507 u = u.ScaledBy(umax - umin);
508 v = v.ScaledBy(vmax - vmin);
509
510 // So we can now generate the top and bottom surfaces of the extrusion,
511 // planes within a translated (and maybe mirrored) version of that csys.
512 SSurface s0, s1;
513 s0 = SSurface::FromPlane(orig.Plus(t0), u, v);
514 s0.color = color;
515 s1 = SSurface::FromPlane(orig.Plus(t1).Plus(u), u.ScaledBy(-1), v);
516 s1.color = color;
517 hSSurface hs0 = surface.AddAndAssignId(&s0),
518 hs1 = surface.AddAndAssignId(&s1);
519
520 // Now go through the input curves. For each one, generate its surface
521 // of extrusion, its two translated trim curves, and one trim line. We
522 // go through by loops so that we can assign the lines correctly.
523 SBezierLoop *sbl;
524 for(sbl = sbls->l.First(); sbl; sbl = sbls->l.NextAfter(sbl)) {
525 SBezier *sb;
526 List<TrimLine> trimLines = {};
527
528 for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {
529 // Generate the surface of extrusion of this curve, and add
530 // it to the list
531 SSurface ss = SSurface::FromExtrusionOf(sb, t0, t1);
532 ss.color = color;
533 hSSurface hsext = surface.AddAndAssignId(&ss);
534
535 // Translate the curve by t0 and t1 to produce two trim curves
536 SCurve sc = {};
537 sc.isExact = true;
538 sc.exact = sb->TransformedBy(t0, Quaternion::IDENTITY, 1.0);
539 (sc.exact).MakePwlInto(&(sc.pts));
540 sc.surfA = hs0;
541 sc.surfB = hsext;
542 hSCurve hc0 = curve.AddAndAssignId(&sc);
543
544 sc = {};
545 sc.isExact = true;
546 sc.exact = sb->TransformedBy(t1, Quaternion::IDENTITY, 1.0);
547 (sc.exact).MakePwlInto(&(sc.pts));
548 sc.surfA = hs1;
549 sc.surfB = hsext;
550 hSCurve hc1 = curve.AddAndAssignId(&sc);
551
552 STrimBy stb0, stb1;
553 // The translated curves trim the flat top and bottom surfaces.
554 stb0 = STrimBy::EntireCurve(this, hc0, false);
555 stb1 = STrimBy::EntireCurve(this, hc1, true);
556 (surface.FindById(hs0))->trim.Add(&stb0);
557 (surface.FindById(hs1))->trim.Add(&stb1);
558
559 // The translated curves also trim the surface of extrusion.
560 stb0 = STrimBy::EntireCurve(this, hc0, true);
561 stb1 = STrimBy::EntireCurve(this, hc1, false);
562 (surface.FindById(hsext))->trim.Add(&stb0);
563 (surface.FindById(hsext))->trim.Add(&stb1);
564
565 // And form the trim line
566 Vector pt = sb->Finish();
567 sc = {};
568 sc.isExact = true;
569 sc.exact = SBezier::From(pt.Plus(t0), pt.Plus(t1));
570 (sc.exact).MakePwlInto(&(sc.pts));
571 hSCurve hl = curve.AddAndAssignId(&sc);
572 // save this for later
573 TrimLine tl;
574 tl.hc = hl;
575 tl.hs = hsext;
576 trimLines.Add(&tl);
577 }
578
579 int i;
580 for(i = 0; i < trimLines.n; i++) {
581 TrimLine *tl = &(trimLines.elem[i]);
582 SSurface *ss = surface.FindById(tl->hs);
583
584 TrimLine *tlp = &(trimLines.elem[WRAP(i-1, trimLines.n)]);
585
586 STrimBy stb;
587 stb = STrimBy::EntireCurve(this, tl->hc, true);
588 ss->trim.Add(&stb);
589 stb = STrimBy::EntireCurve(this, tlp->hc, false);
590 ss->trim.Add(&stb);
591
592 (curve.FindById(tl->hc))->surfA = ss->h;
593 (curve.FindById(tlp->hc))->surfB = ss->h;
594 }
595 trimLines.Clear();
596 }
597 }
598
599
600 typedef struct {
601 hSSurface d[4];
602 } Revolved;
603
MakeFromRevolutionOf(SBezierLoopSet * sbls,Vector pt,Vector axis,RgbaColor color,Group * group)604 void SShell::MakeFromRevolutionOf(SBezierLoopSet *sbls, Vector pt, Vector axis, RgbaColor color, Group *group)
605 {
606 SBezierLoop *sbl;
607
608 int i0 = surface.n, i;
609
610 // Normalize the axis direction so that the direction of revolution
611 // ends up parallel to the normal of the sketch, on the side of the
612 // axis where the sketch is.
613 Vector pto;
614 double md = VERY_NEGATIVE;
615 for(sbl = sbls->l.First(); sbl; sbl = sbls->l.NextAfter(sbl)) {
616 SBezier *sb;
617 for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {
618 // Choose the point farthest from the axis; we'll get garbage
619 // if we choose a point that lies on the axis, for example.
620 // (And our surface will be self-intersecting if the sketch
621 // spans the axis, so don't worry about that.)
622 Vector p = sb->Start();
623 double d = p.DistanceToLine(pt, axis);
624 if(d > md) {
625 md = d;
626 pto = p;
627 }
628 }
629 }
630 Vector ptc = pto.ClosestPointOnLine(pt, axis),
631 up = (pto.Minus(ptc)).WithMagnitude(1),
632 vp = (sbls->normal).Cross(up);
633 if(vp.Dot(axis) < 0) {
634 axis = axis.ScaledBy(-1);
635 }
636
637 // Now we actually build and trim the surfaces.
638 for(sbl = sbls->l.First(); sbl; sbl = sbls->l.NextAfter(sbl)) {
639 int i, j;
640 SBezier *sb;
641 List<Revolved> hsl = {};
642
643 for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {
644 Revolved revs;
645 for(j = 0; j < 4; j++) {
646 if(sb->deg == 1 &&
647 (sb->ctrl[0]).DistanceToLine(pt, axis) < LENGTH_EPS &&
648 (sb->ctrl[1]).DistanceToLine(pt, axis) < LENGTH_EPS)
649 {
650 // This is a line on the axis of revolution; it does
651 // not contribute a surface.
652 revs.d[j].v = 0;
653 } else {
654 SSurface ss = SSurface::FromRevolutionOf(sb, pt, axis,
655 (PI/2)*j,
656 (PI/2)*(j+1));
657 ss.color = color;
658 if(sb->entity != 0) {
659 hEntity he;
660 he.v = sb->entity;
661 hEntity hface = group->Remap(he, Group::REMAP_LINE_TO_FACE);
662 if(SK.entity.FindByIdNoOops(hface) != NULL) {
663 ss.face = hface.v;
664 }
665 }
666 revs.d[j] = surface.AddAndAssignId(&ss);
667 }
668 }
669 hsl.Add(&revs);
670 }
671
672 for(i = 0; i < sbl->l.n; i++) {
673 Revolved revs = hsl.elem[i],
674 revsp = hsl.elem[WRAP(i-1, sbl->l.n)];
675
676 sb = &(sbl->l.elem[i]);
677
678 for(j = 0; j < 4; j++) {
679 SCurve sc;
680 Quaternion qs = Quaternion::From(axis, (PI/2)*j);
681 // we want Q*(x - p) + p = Q*x + (p - Q*p)
682 Vector ts = pt.Minus(qs.Rotate(pt));
683
684 // If this input curve generate a surface, then trim that
685 // surface with the rotated version of the input curve.
686 if(revs.d[j].v) {
687 sc = {};
688 sc.isExact = true;
689 sc.exact = sb->TransformedBy(ts, qs, 1.0);
690 (sc.exact).MakePwlInto(&(sc.pts));
691 sc.surfA = revs.d[j];
692 sc.surfB = revs.d[WRAP(j-1, 4)];
693
694 hSCurve hcb = curve.AddAndAssignId(&sc);
695
696 STrimBy stb;
697 stb = STrimBy::EntireCurve(this, hcb, true);
698 (surface.FindById(sc.surfA))->trim.Add(&stb);
699 stb = STrimBy::EntireCurve(this, hcb, false);
700 (surface.FindById(sc.surfB))->trim.Add(&stb);
701 }
702
703 // And if this input curve and the one after it both generated
704 // surfaces, then trim both of those by the appropriate
705 // circle.
706 if(revs.d[j].v && revsp.d[j].v) {
707 SSurface *ss = surface.FindById(revs.d[j]);
708
709 sc = {};
710 sc.isExact = true;
711 sc.exact = SBezier::From(ss->ctrl[0][0],
712 ss->ctrl[0][1],
713 ss->ctrl[0][2]);
714 sc.exact.weight[1] = ss->weight[0][1];
715 (sc.exact).MakePwlInto(&(sc.pts));
716 sc.surfA = revs.d[j];
717 sc.surfB = revsp.d[j];
718
719 hSCurve hcc = curve.AddAndAssignId(&sc);
720
721 STrimBy stb;
722 stb = STrimBy::EntireCurve(this, hcc, false);
723 (surface.FindById(sc.surfA))->trim.Add(&stb);
724 stb = STrimBy::EntireCurve(this, hcc, true);
725 (surface.FindById(sc.surfB))->trim.Add(&stb);
726 }
727 }
728 }
729
730 hsl.Clear();
731 }
732
733 for(i = i0; i < surface.n; i++) {
734 SSurface *srf = &(surface.elem[i]);
735
736 // Revolution of a line; this is potentially a plane, which we can
737 // rewrite to have degree (1, 1).
738 if(srf->degm == 1 && srf->degn == 2) {
739 // close start, far start, far finish
740 Vector cs, fs, ff;
741 double d0, d1;
742 d0 = (srf->ctrl[0][0]).DistanceToLine(pt, axis);
743 d1 = (srf->ctrl[1][0]).DistanceToLine(pt, axis);
744
745 if(d0 > d1) {
746 cs = srf->ctrl[1][0];
747 fs = srf->ctrl[0][0];
748 ff = srf->ctrl[0][2];
749 } else {
750 cs = srf->ctrl[0][0];
751 fs = srf->ctrl[1][0];
752 ff = srf->ctrl[1][2];
753 }
754
755 // origin close, origin far
756 Vector oc = cs.ClosestPointOnLine(pt, axis),
757 of = fs.ClosestPointOnLine(pt, axis);
758
759 if(oc.Equals(of)) {
760 // This is a plane, not a (non-degenerate) cone.
761 Vector oldn = srf->NormalAt(0.5, 0.5);
762
763 Vector u = fs.Minus(of), v;
764
765 v = (axis.Cross(u)).WithMagnitude(1);
766
767 double vm = (ff.Minus(of)).Dot(v);
768 v = v.ScaledBy(vm);
769
770 srf->degm = 1;
771 srf->degn = 1;
772 srf->ctrl[0][0] = of;
773 srf->ctrl[0][1] = of.Plus(u);
774 srf->ctrl[1][0] = of.Plus(v);
775 srf->ctrl[1][1] = of.Plus(u).Plus(v);
776 srf->weight[0][0] = 1;
777 srf->weight[0][1] = 1;
778 srf->weight[1][0] = 1;
779 srf->weight[1][1] = 1;
780
781 if(oldn.Dot(srf->NormalAt(0.5, 0.5)) < 0) {
782 swap(srf->ctrl[0][0], srf->ctrl[1][0]);
783 swap(srf->ctrl[0][1], srf->ctrl[1][1]);
784 }
785 continue;
786 }
787
788 if(fabs(d0 - d1) < LENGTH_EPS) {
789 // This is a cylinder; so transpose it so that we'll recognize
790 // it as a surface of extrusion.
791 SSurface sn = *srf;
792
793 // Transposing u and v flips the normal, so reverse u to
794 // flip it again and put it back where we started.
795 sn.degm = 2;
796 sn.degn = 1;
797 int dm, dn;
798 for(dm = 0; dm <= 1; dm++) {
799 for(dn = 0; dn <= 2; dn++) {
800 sn.ctrl [dn][dm] = srf->ctrl [1-dm][dn];
801 sn.weight[dn][dm] = srf->weight[1-dm][dn];
802 }
803 }
804
805 *srf = sn;
806 continue;
807 }
808 }
809
810 }
811
812 }
813
MakeFromCopyOf(SShell * a)814 void SShell::MakeFromCopyOf(SShell *a) {
815 MakeFromTransformationOf(a,
816 Vector::From(0, 0, 0), Quaternion::IDENTITY, 1.0);
817 }
818
MakeFromTransformationOf(SShell * a,Vector t,Quaternion q,double scale)819 void SShell::MakeFromTransformationOf(SShell *a,
820 Vector t, Quaternion q, double scale)
821 {
822 booleanFailed = false;
823
824 SSurface *s;
825 for(s = a->surface.First(); s; s = a->surface.NextAfter(s)) {
826 SSurface n;
827 n = SSurface::FromTransformationOf(s, t, q, scale, true);
828 surface.Add(&n); // keeping the old ID
829 }
830
831 SCurve *c;
832 for(c = a->curve.First(); c; c = a->curve.NextAfter(c)) {
833 SCurve n;
834 n = SCurve::FromTransformationOf(c, t, q, scale);
835 curve.Add(&n); // keeping the old ID
836 }
837 }
838
MakeEdgesInto(SEdgeList * sel)839 void SShell::MakeEdgesInto(SEdgeList *sel) {
840 SSurface *s;
841 for(s = surface.First(); s; s = surface.NextAfter(s)) {
842 s->MakeEdgesInto(this, sel, SSurface::AS_XYZ);
843 }
844 }
845
MakeSectionEdgesInto(Vector n,double d,SEdgeList * sel,SBezierList * sbl)846 void SShell::MakeSectionEdgesInto(Vector n, double d,
847 SEdgeList *sel, SBezierList *sbl)
848 {
849 SSurface *s;
850 for(s = surface.First(); s; s = surface.NextAfter(s)) {
851 if(s->CoincidentWithPlane(n, d)) {
852 s->MakeSectionEdgesInto(this, sel, sbl);
853 }
854 }
855 }
856
TriangulateInto(SMesh * sm)857 void SShell::TriangulateInto(SMesh *sm) {
858 SSurface *s;
859 for(s = surface.First(); s; s = surface.NextAfter(s)) {
860 s->TriangulateInto(this, sm);
861 }
862 }
863
IsEmpty(void)864 bool SShell::IsEmpty(void) {
865 return (surface.n == 0);
866 }
867
Clear(void)868 void SShell::Clear(void) {
869 SSurface *s;
870 for(s = surface.First(); s; s = surface.NextAfter(s)) {
871 s->Clear();
872 }
873 surface.Clear();
874
875 SCurve *c;
876 for(c = curve.First(); c; c = curve.NextAfter(c)) {
877 c->Clear();
878 }
879 curve.Clear();
880 }
881
882