From 246e984e6401d0777efccb109548fe0fbfc2165b Mon Sep 17 00:00:00 2001 From: Val Erastov Date: Thu, 19 Oct 2017 20:07:43 -0700 Subject: [PATCH] boolean / nurbs curve --- web/app/3d/craft/brep/cut-extrude.js | 1 - web/app/3d/craft/sketch/sketch-model.js | 22 +- web/app/3d/debug.js | 5 +- web/app/3d/modeler-app.js | 56 ++-- web/app/3d/scene/brep-scene-object.js | 2 +- web/app/3d/tess/triangulation.js | 2 +- web/app/brep/brep-builder.js | 10 +- web/app/brep/brep-enclose.js | 3 +- web/app/brep/geom/curve.js | 26 -- web/app/brep/geom/impl/approx.js | 126 ------- web/app/brep/geom/impl/curve/curve-tess.js | 58 ++++ web/app/brep/geom/impl/curve/curves-isec.js | 100 ++++++ web/app/brep/geom/impl/line.js | 8 +- .../geom/impl/{nurbs-impl.js => nurbs-ext.js} | 110 +++--- web/app/brep/geom/impl/nurbs.js | 312 +++++++++++++----- web/app/brep/geom/tolerance.js | 21 +- web/app/brep/operations/boolean.js | 3 +- web/app/brep/operations/evolve-face.js | 2 +- web/app/brep/topo/loop.js | 5 +- web/app/math/math.js | 2 + web/app/math/vec.js | 46 +-- web/lib/verb.js | 6 + 22 files changed, 570 insertions(+), 356 deletions(-) delete mode 100644 web/app/brep/geom/impl/approx.js create mode 100644 web/app/brep/geom/impl/curve/curve-tess.js create mode 100644 web/app/brep/geom/impl/curve/curves-isec.js rename web/app/brep/geom/impl/{nurbs-impl.js => nurbs-ext.js} (73%) diff --git a/web/app/3d/craft/brep/cut-extrude.js b/web/app/3d/craft/brep/cut-extrude.js index 0cd633a0..790bd266 100644 --- a/web/app/3d/craft/brep/cut-extrude.js +++ b/web/app/3d/craft/brep/cut-extrude.js @@ -6,7 +6,6 @@ import * as stitching from '../../../brep/stitching' import {Loop} from '../../../brep/topo/loop' import {incRefCounter} from '../../../brep/topo/topo-object' import {Line} from '../../../brep/geom/impl/line' -import {CompositeCurve} from '../../../brep/geom/curve' import {ReadSketchFromFace} from '../sketch/sketch-reader' import {Segment} from '../sketch/sketch-model' import {isCurveClass} from '../../cad-utils' diff --git a/web/app/3d/craft/sketch/sketch-model.js b/web/app/3d/craft/sketch/sketch-model.js index 33178b0f..2bc66701 100644 --- a/web/app/3d/craft/sketch/sketch-model.js +++ b/web/app/3d/craft/sketch/sketch-model.js @@ -1,6 +1,4 @@ -import {CompositeCurve} from '../../../brep/geom/curve' -import {ApproxCurve} from '../../../brep/geom/impl/approx' -import {NurbsCurve} from '../../../brep/geom/impl/nurbs' +import {NurbsCurve, NurbsCurveImpl} from '../../../brep/geom/impl/nurbs' import {Point} from '../../../brep/geom/point' import {Line} from '../../../brep/geom/impl/line' import {LUT} from '../../../math/bezier-cubic' @@ -38,7 +36,7 @@ class SketchPrimitive { if (this.inverted) { verbNurbs = verbNurbs.reverse(); } - return new NurbsCurve(verbNurbs); + return new NurbsCurve(new NurbsCurveImpl(verbNurbs)); } toVerbNurbs(plane, _3dtr) { @@ -308,3 +306,19 @@ function to3DTrFunc(surface) { return _3dTransformation.apply(v); } } + +class CompositeCurve { + + constructor() { + this.curves = []; + this.points = []; + this.groups = []; + } + + add(curve, point, group) { + this.curves.push(curve); + this.points.push(point); + this.groups.push(group); + } +} + diff --git a/web/app/3d/debug.js b/web/app/3d/debug.js index 4ae214d2..a8c4f09f 100644 --- a/web/app/3d/debug.js +++ b/web/app/3d/debug.js @@ -79,7 +79,7 @@ function addGlobalDebugActions(app) { app.viewer.render(); }, AddHalfEdge: (he, color) => { - const points = he.edge.curve.verb.tessellate().map(p => new Vector().set3(p)); + const points = he.edge.curve.tessellate(); if (he.inverted) { points.reverse(); } @@ -128,8 +128,7 @@ function addGlobalDebugActions(app) { app.viewer.render(); }, AddCurve: (curve, color) => { - __DEBUG__.AddPolyLine( curve.verb.tessellate().map(v => new Vector().set3(v)), 0xffffff); - __DEBUG__.AddPolyLine( curve.tessellate(0.5, 1), color); + __DEBUG__.AddPolyLine( curve.tessellate(), color); }, AddNurbsCorners: (nurbs) => { __DEBUG__.AddPoint(nurbs.point(0, 0), 0xff0000); diff --git a/web/app/3d/modeler-app.js b/web/app/3d/modeler-app.js index f3ed1d37..201dd452 100644 --- a/web/app/3d/modeler-app.js +++ b/web/app/3d/modeler-app.js @@ -24,7 +24,10 @@ import * as BREPBool from '../brep/operations/boolean' import {BREPValidator} from '../brep/brep-validator' import {BREPSceneSolid} from './scene/brep-scene-object' import TPI from './tpi' -import {NurbsCurve, NurbsSurface} from "../brep/geom/impl/nurbs"; +import {NurbsCurve, NurbsCurveImpl, NurbsSurface} from "../brep/geom/impl/nurbs"; +import {Circle} from "./craft/sketch/sketch-model"; +import {Plane} from "../brep/geom/impl/plane"; +import {enclose} from "../brep/brep-enclose"; // import {createSphere, rayMarchOntoCanvas, sdfIntersection, sdfSolid, sdfSubtract, sdfTransform, sdfUnion} from "../hds/sdf"; function App() { @@ -113,7 +116,16 @@ App.prototype.test1 = function() { App.prototype.cylTest = function() { const cylinder1 = BREPPrimitives.cylinder(200, 500); - const cylinder2 = BREPPrimitives.cylinder(200, 500, Matrix3.rotateMatrix(90, AXIS.Y, ORIGIN)); + + + // const cylinder2 = (function () { + // let circle1 = new Circle(-1, new Vector(0,0,0), 200).toNurbs( new Plane(AXIS.X, 500)); + // let circle2 = circle1.translate(new Vector(-1000,0,0)); + // return enclose([circle1], [circle2]) + // })(); + + + const cylinder2 = BREPPrimitives.cylinder(200, 500, Matrix3.rotateMatrix(90, AXIS.Y, ORIGIN)); let result = this.TPI.brep.bool.subtract(cylinder1, cylinder2); @@ -201,18 +213,21 @@ App.prototype.test5 = function() { const g = vx(0.73, 0.73); const h = vx(0.73, 0.33); + function fromVerb(verb) { + return new NurbsCurve(new NurbsCurveImpl(verb)); + } let shell = bb.face(srf) .loop() - .edgeTrim(a, b, srf.verb.isocurve(0.13, true)) - .edgeTrim(b, c, srf.verb.isocurve(0.9, false)) - .edgeTrim(c, d, srf.verb.isocurve(0.9, true).reverse()) - .edgeTrim(d, a, srf.verb.isocurve(0.13, false).reverse()) + .edgeTrim(a, b, fromVerb(srf.verb.isocurve(0.13, true))) + .edgeTrim(b, c, fromVerb(srf.verb.isocurve(0.9, false))) + .edgeTrim(c, d, fromVerb(srf.verb.isocurve(0.9, true).reverse())) + .edgeTrim(d, a, fromVerb(srf.verb.isocurve(0.13, false).reverse())) .loop() - .edgeTrim(e, f, srf.verb.isocurve(0.33, false)) - .edgeTrim(f, g, srf.verb.isocurve(0.73, true)) - .edgeTrim(g, h, srf.verb.isocurve(0.73, false).reverse()) - .edgeTrim(h, e, srf.verb.isocurve(0.33, true).reverse()) + .edgeTrim(e, f, fromVerb(srf.verb.isocurve(0.33, false))) + .edgeTrim(f, g, fromVerb(srf.verb.isocurve(0.73, true))) + .edgeTrim(g, h, fromVerb(srf.verb.isocurve(0.73, false).reverse())) + .edgeTrim(h, e, fromVerb(srf.verb.isocurve(0.33, true).reverse())) .build(); this.addShellOnScene(shell); @@ -220,8 +235,8 @@ App.prototype.test5 = function() { App.prototype.scratchCode = function() { // const app = this; - this.test3(); - // this.cylTest(); + // this.test5(); + this.cylTest(); return @@ -230,21 +245,22 @@ return var p1 = [-50,0,0], p2 = [100,0,0], p3 = [100,100,0], p4 = [0,100,0], p5 = [50, 50, 0]; var pts = [p1, p2, p3, p4, p5]; - let curve1 = new NurbsCurve(verb.geom.NurbsCurve.byPoints( pts, 3 )); + let curve1 = new NurbsCurve(new NurbsCurveImpl(verb.geom.NurbsCurve.byPoints( pts, 3 ))); var p1a = [-50,0,0], p2a = [50,-10,0], p3a = [150,50,0], p4a = [30,100,0], p5a = [50, 120, 0]; var ptsa = [p1a, p2a, p3a, p4a, p5a]; - let curve2 = new NurbsCurve(verb.geom.NurbsCurve.byPoints( ptsa, 3 )); + let curve2 = new NurbsCurve(new NurbsCurveImpl(verb.geom.NurbsCurve.byPoints( ptsa, 3 ))); - // __DEBUG__.AddCurve(curve1); + curve1 = curve1.splitByParam(0.6)[0]; + __DEBUG__.AddCurve(curve1); __DEBUG__.AddCurve(curve2); - // let points = curve1.intersectCurve(curve2); - // for (let p of points) { - // __DEBUG__.AddPoint(p.p0); - // } + let points = curve1.intersectCurve(curve2); + for (let p of points) { + __DEBUG__.AddPoint(p.p0); + } - app.viewer.render(); + // app.viewer.render(); }; diff --git a/web/app/3d/scene/brep-scene-object.js b/web/app/3d/scene/brep-scene-object.js index 41aa74ff..387dc555 100644 --- a/web/app/3d/scene/brep-scene-object.js +++ b/web/app/3d/scene/brep-scene-object.js @@ -44,7 +44,7 @@ export class BREPSceneSolid extends SceneSolid { for (let edge of this.shell.edges) { if (edge.data[EDGE_AUX] === undefined) { const line = new THREE.Line(undefined, WIREFRAME_MATERIAL); - const contour = edge.curve.verb.tessellate(); + const contour = edge.curve.tessellate(); for (let p of contour) { line.geometry.vertices.push(new THREE.Vector3().fromArray(p)); } diff --git a/web/app/3d/tess/triangulation.js b/web/app/3d/tess/triangulation.js index 0c87703f..6ab70f9e 100644 --- a/web/app/3d/tess/triangulation.js +++ b/web/app/3d/tess/triangulation.js @@ -118,7 +118,7 @@ export function TriangulateFace(face) { tessy.gluTessVertex(arr(e.vertexA.point), e.vertexA); } else if (e.edge.curve.verb) { tessy.gluTessVertex(arr(e.vertexA.point), e.vertexA); - let points = e.edge.curve.verb.tessellate(); + let points = e.edge.curve.tessellate(); if (e.inverted) { points.reverse(); } diff --git a/web/app/brep/brep-builder.js b/web/app/brep/brep-builder.js index 34d87740..864e9a95 100644 --- a/web/app/brep/brep-builder.js +++ b/web/app/brep/brep-builder.js @@ -40,11 +40,11 @@ export default class BrepBuilder { } edgeTrim(a, b, curve) { - let u1 = curve.closestParam(a.point.data()); - curve = curve.split(u1)[1]; - let u2 = curve.closestParam(b.point.data()); - curve = curve.split(u2)[0]; - this.edge(a, b, new NurbsCurve(curve)); + let u1 = curve.param(a.point); + curve = curve.splitByParam(u1)[1]; + let u2 = curve.param(b.point); + curve = curve.splitByParam(u2)[0]; + this.edge(a, b, curve); return this; } diff --git a/web/app/brep/brep-enclose.js b/web/app/brep/brep-enclose.js index 9895b27f..eb69b74c 100644 --- a/web/app/brep/brep-enclose.js +++ b/web/app/brep/brep-enclose.js @@ -8,7 +8,6 @@ import {NurbsSurface, NurbsCurve} from './geom/impl/nurbs' import {Plane} from './geom/impl/plane' import {Point} from './geom/point' import {BasisForPlane, Matrix3} from '../math/l3space' -import {CompositeCurve} from './geom/curve' import * as cad_utils from '../3d/cad-utils' import * as math from '../math/math' import mergeNullFace from './null-face-merge' @@ -124,7 +123,7 @@ export function createWall(curve1, curve2) { if (bothClassOf(curve1, curve2, 'Line')) { throw 'unsupported' } else if (bothClassOf(curve1, curve2, 'NurbsCurve')) { - return new NurbsSurface(verb.geom.NurbsSurface.byLoftingCurves([curve2.verb, curve1.verb], 1)); + return NurbsSurface.loft(curve2, curve1, 1); } else { throw 'unsupported'; } diff --git a/web/app/brep/geom/curve.js b/web/app/brep/geom/curve.js index a498f0a2..35acbdc2 100644 --- a/web/app/brep/geom/curve.js +++ b/web/app/brep/geom/curve.js @@ -1,30 +1,4 @@ -export class Curve { - - constructor() { - } - - isSameClass(other) { - return this.constructor.name == other.constructor.name; - } - - parametricEquation(t) { - throw 'not implemented'; - } - - translate(vector) { - throw 'not implemented'; - } - - toNurbs(domainA, domainB) { - throw 'not implemented'; - } - - approximate(resolution, from, to, path) { - } -} -Curve.prototype.isLine = false; - export class CompositeCurve { constructor() { diff --git a/web/app/brep/geom/impl/approx.js b/web/app/brep/geom/impl/approx.js deleted file mode 100644 index 06b45e6d..00000000 --- a/web/app/brep/geom/impl/approx.js +++ /dev/null @@ -1,126 +0,0 @@ -import {Surface} from '../surface' -import {Line} from '../impl/line' -import {Curve} from '../curve' -import {Matrix3, AXIS, BasisForPlane} from '../../../math/l3space' -import * as math from '../../../math/math' - -export class ApproxSurface extends Surface { - constructor(mesh) { - super(); - this.mesh = mesh; - } -} - -export class ApproxCurve extends Curve { - constructor(points, proto) { - super(); - this.points = points; - this.segments = []; - this.tShifts = [0]; - for (let i = 1; i < points.length; ++i) { - const a = points[i - 1]; - const b = points[i]; - const line = Line.fromSegment(a, b); - this.segments.push(line); - this.tShifts.push(this.tShifts[this.tShifts.length - 1] + line.t(b)); - } - this.proto = proto; - } - - - t(point) { - for (let i = 0; i < this.points.length; ++i) { - if (math.vectorsEqual(this.points[i], point)) { - return this.tShifts[i]; - } - } - - for (let i = 0; i < this.segments.length; ++i) { - const line = this.segments[i]; - const subT = line.t(point); - if (subT > 0 && subT < this.tShifts[i + 1]) { - return this.tShifts[i] + subT; - } - } - return NaN; - } - - parametricEquation(t) { - for (let i = 0; i < this.points.length; ++i) { - if (math.equal(t, this.tShifts[i])) { - return this.points[i]; - } - } - for (let i = 1; i < this.points.length; ++i) { - if (t > this.tShifts[i - 1] && t < this.tShifts[i]) { - return this.segments[i - 1].parametricEquation(t - this.tShifts[i - 1]); - } - } - return null; - } - - getChunk(p1, p2) { - const result = []; - this.getChunkImpl(p1, p2, result, true); - return result; - } - - getChunkImpl(p1, p2, result, includeBounds) { - const t1 = this.t(p1); - const t2 = this.t(p2); - if (t1 > t2) { - let tmp = p1; - p1 = p2; - p2 = tmp; - } - return this.getChunkDirectional(p1, p2, result, includeBounds); - } - - getChunkDirectional(p1, p2, result, includeBounds) { - let inState = false; - for (let i = 1; i < this.points.length; ++i) { - const a = this.points[i - 1]; - const b = this.points[i]; - const line = this.segments[i - 1]; - if (inState) { - result.push(a); - } - if (!inState) { - if (math.vectorsEqual(a, p1)) { - if (includeBounds) result.push(p1); - inState = true; - } else if (math.equal(b, p1)) { - //nothing - } else { - const t = line.t(p1); - if (t > 0 && t < this.tShifts[i] - this.tShifts[i - 1]) { - if (includeBounds) result.push(p1); - inState = true; - } - } - } - if (inState) { - if (math.vectorsEqual(b, p2)) { - if (includeBounds) result.push(p2); - break; - } else if (math.equal(a, p2)) { - //nothing, can't be here - } else { - const t = line.t(p2); - if (t > 0 && t < this.tShifts[i] - this.tShifts[i - 1]) { - if (includeBounds) result.push(p2); - break; - } - } - } - } - } - - translate(vector) { - return new ApproxCurve(this.points.map(p => p.plus(vector)), this.proto); - } - - approximate(resolution, from, to, path) { - this.getChunkImpl(from, to, path, false); - } -} diff --git a/web/app/brep/geom/impl/curve/curve-tess.js b/web/app/brep/geom/impl/curve/curve-tess.js new file mode 100644 index 00000000..be6a19cf --- /dev/null +++ b/web/app/brep/geom/impl/curve/curve-tess.js @@ -0,0 +1,58 @@ +import * as vec from "../../../../math/vec"; + +export default function curveTess(curve, min, max, tessTol, scale) { + return curveTessParams(curve, min, max, tessTol, scale).map(u => curve.point(u)); +} + +export function curveTessParams(curve, min, max, tessTol, scale) { + + let [dmin, dmax] = curve.domain(); + + let out = []; + let nSplits = curve.optimalSplits(); + + let splitStep = (dmax - dmin) / nSplits; + nSplits = Math.round((max - min) / splitStep); + splitStep = (max - min) / nSplits; + + let splits = [min]; + + for (let i = 1; i < nSplits; ++i) { + splits.push(min + i * splitStep); + } + splits.push(max); + + function refine(u1, u2, step) { + if (step < u2 - u1) { + let mid = u1 + (u2 - u1) * 0.5; + refine(u1, mid, step); + out.push(mid); + refine(mid, u2, curveStep(curve, mid, tessTol, scale)); + } + } + for (let i = 1; i < splits.length; ++i) { + let u1 = splits[i - 1]; + out.push(u1); + refine(u1, splits[i], curveStep(curve, u1, tessTol, scale)); + } + + out.push(max); + return out; +} + +export function curveStep(curve, u, tessTol, scale) { + tessTol = tessTol || 1; + scale = scale || 1; + let ders = curve.eval(u, 2); + let r1 = ders[1]; + let r2 = ders[2]; + + let r1lsq = vec.lengthSq(r1); + let r1l = Math.sqrt(r1lsq); + + let r = r1lsq * r1l / vec.length(vec.cross(r1, r2)); + let tol = tessTol / scale; + + let step = 2 * Math.sqrt(tol*(2*r - tol)) / r1l; + return step; +} \ No newline at end of file diff --git a/web/app/brep/geom/impl/curve/curves-isec.js b/web/app/brep/geom/impl/curve/curves-isec.js new file mode 100644 index 00000000..8af2b1cf --- /dev/null +++ b/web/app/brep/geom/impl/curve/curves-isec.js @@ -0,0 +1,100 @@ +import * as vec from "../../../../math/vec"; +import {TOLERANCE, TOLERANCE_SQ} from '../../tolerance'; +import * as math from '../../../../math/math' +import {fmin_bfgs} from "../../../../math/optim"; + +export default function curveIntersect(curve1, curve2, isecRange1, isecRange2, tesselator) { + + let result = []; + let segs1 = tesselator(curve1, isecRange1[0], isecRange1[1]); + let segs2 = tesselator(curve2, isecRange2[0], isecRange2[1]); + + for (let i = 0; i < segs1.length - 1; i++) { + let a1 = segs1[i]; + let b1 = segs1[i + 1]; + for (let j = 0; j < segs2.length - 1; j++) { + let a2 = segs2[j]; + let b2 = segs2[j + 1]; + + let isec = intersectSegs(a1, b1, a2, b2, TOLERANCE); + if (isec !== null) { + let {point1, point2, l1, l2} = isec; + + let u1 = curve1.param(point1); + let u2 = curve2.param(point2); + [u1, u2] = curveExactIntersection(curve1, curve2, u1, u2); + + if (!isInRange(u1, isecRange1) || !isInRange(u2, isecRange2)) { + continue; + } + + result.push({ + u0: u1, + u1: u2, + p0: point1, + p1: point2 + }); + if (math.areEqual(u1, l1, TOLERANCE )) { + i ++; + } + if (math.areEqual(u2, l2, TOLERANCE )) { + j ++; + } + } + } + } + return result; +} + +function curveExactIntersection(curve1, curve2, u1, u2) { + + function f([u1, u2]) { + return vec.lengthSq( vec.sub(curve1.point(u1), curve2.point(u2))); + } + function grad([u1, u2]) { + let d1 = curve1.eval( u1, 1); + let d2 = curve2.eval( u2, 1); + let r = vec.sub(d1[0], d2[0]); + let drdu = d1[1]; + let drdt = vec.mul(-1, d2[1]); + return [2 * vec.dot(drdu, r), 2 * vec.dot(drdt,r)]; + } + let params = [u1, u2]; + return fmin_bfgs(f, params, TOLERANCE_SQ, grad).solution; +} + +function lineLineIntersection(p1, p2, v1, v2) { + let zAx = vec.cross(v1, v2); + const n1 = vec._normalize(vec.cross(zAx, v1)); + const n2 = vec._normalize(vec.cross(zAx, v2)); + return { + u1: vec.dot(n2, vec.sub(p2, p1)) / vec.dot(n2, v1), + u2: vec.dot(n1, vec.sub(p1, p2)) / vec.dot(n1, v2), + } +} + +function intersectSegs(a1, b1, a2, b2) { + let v1 = vec.sub(b1, a1); + let v2 = vec.sub(b2, a2); + let l1 = vec.length(v1); + let l2 = vec.length(v2); + vec._div(v1, l1); + vec._div(v2, l2); + + let {u1, u2} = lineLineIntersection(a1, a2, v1, v2); + let point1 = vec.add(a1, vec.mul(v1, u1)); + let point2 = vec.add(a2, vec.mul(v2, u2)); + let p2p = vec.lengthSq(vec.sub(point1, point2)); + let eq = (a, b) => math.areEqual(a, b, TOLERANCE); + if (u1 !== Infinity && u2 !== Infinity && math.areEqual(p2p, 0, TOLERANCE_SQ) && + ((u1 >0 && u1 < l1) || eq(u1, 0) || eq(u1, l1)) && + ((u2 >0 && u2 < l2) || eq(u2, 0) || eq(u2, l2)) + ) { + return {point1, point2, u1, u2, l1, l2} + } + return null; +} + +function isInRange(p, range) { + return p >= range[0] && p <= range[1]; +} diff --git a/web/app/brep/geom/impl/line.js b/web/app/brep/geom/impl/line.js index aed2ca8d..d1e59740 100644 --- a/web/app/brep/geom/impl/line.js +++ b/web/app/brep/geom/impl/line.js @@ -1,9 +1,7 @@ -import {Curve} from '../curve' -export class Line extends Curve { +export class Line { constructor(p0, v) { - super(); throw 'only nurbs for now' this.p0 = p0; this.v = v; @@ -27,7 +25,7 @@ export class Line extends Curve { return super.intersectCurve(curve, surface); } - parametricEquation(t) { + point(t) { return this.p0.plus(this.v.multiply(t)); } @@ -39,7 +37,7 @@ export class Line extends Curve { let point = this._pointsCache.get(surface); if (!point) { const t = this.intersectSurface(surface); - point = this.parametricEquation(t); + point = this.point(t); this._pointsCache.set(surface, point); } return point; diff --git a/web/app/brep/geom/impl/nurbs-impl.js b/web/app/brep/geom/impl/nurbs-ext.js similarity index 73% rename from web/app/brep/geom/impl/nurbs-impl.js rename to web/app/brep/geom/impl/nurbs-ext.js index 8cf2fac7..f855fc2a 100644 --- a/web/app/brep/geom/impl/nurbs-impl.js +++ b/web/app/brep/geom/impl/nurbs-ext.js @@ -1,6 +1,6 @@ import * as vec from "../../../math/vec"; import * as math from '../../../math/math' -import {TOLERANCE, TOLERANCE_SQ} from '../tolerance'; +import {eqEps, TOLERANCE, TOLERANCE_SQ} from '../tolerance'; import {fmin_bfgs} from "../../../math/optim"; export function curveStep(curve, u, tessTol, scale) { @@ -24,7 +24,7 @@ export function curveDomain(curve) { return [curve.knots[0], curve.knots[curve.knots.length - 1]]; } -export function curveParts(curve) { +export function distinctKnots(curve) { let out = [curve.knots[0]]; for (let i = 1; i < curve.knots.length; ++i) { if (out[out.length - 1] !== curve.knots[i]) { @@ -34,19 +34,22 @@ export function curveParts(curve) { return out; } -export function curveTessellateToParams(curve, tessTol, scale) { - let domain = curveDomain(curve); +export function curveTessellate(curve, min, max, tessTol, scale) { if (curve.degree === 1) { - return domain; + return distinctKnots(curve); } + let domain = curveDomain(curve); - let [min, max] = domain; + let [dmin, dmax] = domain; let out = []; let nSplits = curve.knots.length - 1; - let splitStep = (max - min) / nSplits + let splitStep = (dmax - dmin) / nSplits; + nSplits = Math.round((max - min) / splitStep); + splitStep = (max - min) / nSplits; + let splits = [min]; for (let i = 1; i < nSplits; ++i) { @@ -70,40 +73,6 @@ export function curveTessellateToParams(curve, tessTol, scale) { out.push(max); return out; - // let out = []; - // function tessRange(begin, end) { - // let u = begin; - // while (u < end) { - // out.push(u); - // u += curveStep(curve, u, tessTol, scale ); - // } - // } - - // let parts = curveParts(curve); - // for (let i = 1; i < parts.length; ++i) { - // let begin = parts[i - 1]; - // let end = parts[i]; - // tessRange(begin, end); - // } - // out.push(parts[parts.length - 1]); - // return out; -} - -export function curveTessellate(curve, tessTol, scale) { - let params = curveTessellateToParams(curve, tessTol, scale); - let out = []; - if (params.length === 0) { - return out; - } - out.push(curvePoint(curve, params[0])); - - for (let i = 1; i < params.length; ++i) { - let p = curvePoint(curve, params[i]); - if (!math.areVectorsEqual3(out[out.length - 1], p, TOLERANCE)) { - out.push(p); - } - } - return out; } export function curvePoint(curve, u) { @@ -115,19 +84,38 @@ export function curveClosestParam(curve, point) { } export function surfaceIntersect(surface0, surface1) { - const tess1 = verb.eval.Tess.rationalSurfaceAdaptive(surface0); - const tess2 = verb.eval.Tess.rationalSurfaceAdaptive(surface1); - const resApprox = verb.eval.Intersect.meshes(tess1,tess2); + const tess0 = verb.eval.Tess.rationalSurfaceAdaptive(surface0); + const tess1 = verb.eval.Tess.rationalSurfaceAdaptive(surface1); + + function fixTessNaNPoitns(s, tess) { + for (let i = 0; i < tess.points.length; i++) { + let pt = tess.points[i]; + if (Number.isNaN(pt[0]) || Number.isNaN(pt[1]) || Number.isNaN(pt[2])) { + let [u, v] = tess.uvs[i]; + tess.points[i] = verb.eval.Eval.rationalSurfacePoint(s, u, v); + } + } + } + + fixTessNaNPoitns(surface0, tess0); + fixTessNaNPoitns(surface1, tess1); + + const resApprox = verb.eval.Intersect.meshes(tess0,tess1); const exactPls = resApprox.map(function(pl) { return pl.map(function(inter) { return verb.eval.Intersect.surfacesAtPointWithEstimate(surface0,surface1,inter.uv0,inter.uv1,TOLERANCE); }); }); + + + //temporary workaround + return exactPls.map(pl => verb.eval.Make.polyline(pl.map(ip => ip.point))); + return exactPls.map(function(x) { return verb.eval.Make.rationalInterpCurve(x.map(function(y) { return y.point; }), surfaceMaxDegree(surface0) === 1 && surfaceMaxDegree(surface1) === 1 ? 1 : x.length - 1); - }).map(cd => new verb.geom.NurbsCurve(cd)); + }); } export function surfaceMaxDegree(surface) { @@ -159,8 +147,8 @@ export function curveIntersect(curve1, curve2) { result.push({ u0: u1, u1: u2, - point0: point1, - point1: point2 + p0: point1, + p1: point2 }); if (math.areEqual(u1, l1, TOLERANCE )) { i ++; @@ -223,6 +211,30 @@ function intersectSegs(a1, b1, a2, b2) { return null; } -function dist(p1, p2) { - return math.distance3(p1[0], p1[1], p1[2], p2[0], p2[1], p2[2]); +export function normalizeCurveParametrization(curve) { + let [min, max] = curveDomain(curve); + let d = max - min; + for (let i = 0; i < curve.knots.length; i++) { + let val = curve.knots[i]; + if (eqEps(val, min)) { + curve.knots[i] = 0; + } else if (eqEps(val, max)) { + curve.knots[i] = 1; + } else { + curve.knots[i] = (val - min) / d; + } + } + return curve; +} + +export function normalizeCurveParametrizationIfNeeded(curve) { + let [min, max] = curveDomain(curve); + if (min !== 0 || max !== 1) { + normalizeCurveParametrization(curve) + } +} + +export function curveInvert(curve) { + let reversed = verb.eval.Modify.curveReverse(curve); + return reversed; } \ No newline at end of file diff --git a/web/app/brep/geom/impl/nurbs.js b/web/app/brep/geom/impl/nurbs.js index 17c6ba4e..08efb23f 100644 --- a/web/app/brep/geom/impl/nurbs.js +++ b/web/app/brep/geom/impl/nurbs.js @@ -3,146 +3,263 @@ import * as math from '../../../math/math' import {Point} from '../point' import {Surface} from "../surface"; import Vector from "../../../math/vector"; -import {Curve} from "../curve"; -import * as impl from "./nurbs-impl"; +import * as ext from "./nurbs-ext"; +import {EPSILON, eqEps, TOLERANCE} from "../tolerance"; +import curveIntersect from "./curve/curves-isec"; +import curveTess from "./curve/curve-tess"; +import {areEqual} from "../../../math/math"; -export class NurbsCurve extends Curve { + +class ParametricCurve { + + domain() { } + + degree() { } + + degree1Tess() {} + + eval(u, num) { } + + point(param) { } + + param(point) { } + + transform(tr) { } + + optimalSplits() { } + + normalizeParametrization() { } + + invert() { } +} + +export class NurbsCurveImpl { //TODO: rename to NurbsCurve implements ParametricCurve constructor(verbCurve) { - super(); this.verb = verbCurve; this.data = verbCurve.asNurbs(); } - translate(vector) { - const tr = new Matrix3().translate(vector.x, vector.y, vector.z).toArray(); - return new NurbsCurve(this.verb.transform(tr)); - } - - tangentAtPoint(point) { - return pt(this.verb.tangent(this.verb.closestParam(point.data())))._normalize(); + domain() { + return ext.curveDomain(this.data); } - tangentAtParam(param) { - return pt(this.verb.tangent(param ))._normalize(); + degree1Tess() { + return ext.distinctKnots(this.data); } - - closestDistanceToPoint(point) { - const closest = this.verb.closestPoint(point.data()); - return math.distance3(point.x, point.y, point.z, closest[0], closest[1], closest[2]); + + degree() { + return this.data.degree; + } + + transform(tr) { + return new NurbsCurveImpl(this.verb.transform(tr)); + } + + point(u) { + return this.verb.point(u); + } + + param(point) { + return this.verb.closestParam(point); + } + + eval(u, num) { + return verb.eval.Eval.rationalCurveDerivatives( this.data, u, num ); + } + + optimalSplits() { + return this.data.knots.length - 1; + } + + invert() { + + let inverted = ext.curveInvert(this.data); + ext.normalizeCurveParametrizationIfNeeded(inverted); + // let [min, max] = curveDomain(curve); + // for (let i = 0; i < reversed.knots.length; i++) { + // if (eqEps(reversed.knots[i], max)) { + // reversed.knots[i] = max; + // } else { + // break; + // } + // } + // for (let i = reversed.knots.length - 1; i >= 0 ; i--) { + // if (eqEps(reversed.knots[i], min)) { + // reversed.knots[i] = min; + // } else { + // break; + // } + // } + + return new NurbsCurveImpl(newVerbCurve(inverted)); + } + + split(u) { + let split = verb.eval.Divide.curveSplit(this.data, u); + split.forEach(n => ext.normalizeCurveParametrization(n)); + return split.map(c => new NurbsCurveImpl(newVerbCurve(c))); + } +} + +export class NurbsCurve { //TODO: rename to BrepCurve + + constructor(_impl, uMin, uMax) { + let [iMin, iMax] = _impl.domain(); + if (iMin !== 0 || iMax !== 1) { + throw 'only normalized(0..1) parametrization is supported'; + } + this.impl = _impl; + // if (uMin === undefined || uMax === undefined) { + // [uMin, uMax] = this.impl.domain(); + // } + // this.uMin = uMin; + // this.uMax = uMax; + this.uMin = 0; + this.uMax = 1; + } + + translate(vector) { + const tr = new Matrix3().translate(vector.x, vector.y, vector.z); + return new NurbsCurve(this.impl.transform(tr.toArray()), this.uMin, this.uMax); + } + + tangentAtPoint(point) { + let u = this.impl.param(point.data()); + if (areEqual(u, this.uMax, 1e-3)) { // we don't need much tolerance here + //TODO: + // let cps = this.impl.data.controlPoints; + // return pt(cps[cps.length - 1])._minus(pt(cps[cps.length - 2]))._normalize(); + u -= 1e-3; + } + return this.tangentAtParam(u); + } + + tangentAtParam(u) { + const dr = this.impl.eval(u, 1); + return pt(dr[1])._normalize(); + } + + param(point) { + return this.impl.param(point.data()); } split(point) { - return this.splitByParam(this.verb.closestParam(point.data())); + return this.splitByParam(this.param(point)); } splitByParam(u) { - let split = verb.eval.Divide.curveSplit(this.data, u); - split.forEach(n => { - let min = n.knots[0]; - let max = n.knots[n.knots.length - 1]; - let d = max - min; - for (let i = 0; i < n.knots.length; i++) { - let val = n.knots[i]; - if (val === min) { - n.knots[i] = 0; - } else if (val === max) { - n.knots[i] = 1; - } else { - n.knots[i] = (val - min) / d; - } - } - }); - split = split.map(c => new verb.geom.NurbsCurve(c)); + if (u < this.uMin || u > this.uMax) { + throw 'illegal splitting parameter ' + u; + } + let split = this.impl.split(u); + const splitCheck = (split) => { return ( - math.equal(this.verb.closestParam(split[0].point(1)), this.verb.closestParam(split[1].point(0))) && - math.equal(this.verb.closestParam(split[0].point(0)), 0) && - math.equal(this.verb.closestParam(split[0].point(1)), u) && - math.equal(this.verb.closestParam(split[1].point(0)), u) && - math.equal(this.verb.closestParam(split[1].point(1)), 1) + math.equal(this.impl.param(split[0].point(1)), this.impl.param(split[1].point(0))) && + math.equal(this.impl.param(split[0].point(0)), 0) && + math.equal(this.impl.param(split[0].point(1)), u) && + math.equal(this.impl.param(split[1].point(0)), u) && + math.equal(this.impl.param(split[1].point(1)), 1) ) }; if (!splitCheck(split)) { throw 'wrong split'; } - // if (!splitCheck(split)) { - // split.reverse(); - // } return split.map(v => new NurbsCurve(v)); + + // return [ + // new NurbsCurve(this.impl, this.uMin, u), + // new NurbsCurve(this.impl, u, this.uMax) + // ]; } - invert() { - return new NurbsCurve(this.verb.reverse()); - } - point(u) { - return pt(this.verb.point(u)); + return pt(this.impl.point(u)); } tessellate(tessTol, scale) { - return impl.curveTessellate(this.data, tessTol, scale).map(p => pt(p)); + return CURVE_CACHING_TESSELLATOR(this.impl, this.uMin, this.uMax, tessTol, scale).map(p => pt(p)); } - intersectCurve(other, tol) { - let isecs = []; - tol = tol || 1e-3; + boundary() { + return [this.uMin, this.uMax]; + } - const eq = (v1, v2) => math.areVectorsEqual3(v1, v2, tol); + intersectCurve(other) { + let isecs = []; + + const eq = (v1, v2) => math.areVectorsEqual3(v1, v2, TOLERANCE); function add(i0) { for (let i1 of isecs) { if (eq(i0.p0, i1.p0)) { - return; - } - } + return; + } + } isecs.push(i0); } function isecOn(c0, c1, u0) { - const p0 = c0.verb.point(u0); - const u1 = c1.verb.closestParam(p0); - const p1 = c1.verb.point(u1); + const p0 = c0.impl.point(u0); + const u1 = c1.impl.param(p0); + if (!c1.isInside(u1)) { + return; + } + const p1 = c1.impl.point(u1); if (eq(p0, p1)) { if (c0 === other) { add({u0: u1, u1: u0, p0: p1, p1: p0}); } else { add({u0, u1, p0, p1}); } - } } - - isecOn(this, other, 0); - isecOn(this, other, 1); - isecOn(other, this, 0); - isecOn(other, this, 1); - impl.curveIntersect(this.data, other.data, tol).forEach(i => add({ - u0: i.u0, - u1: i.u1, - p0: i.point0, - p1: i.point1 - })); + isecOn(this, other, this.uMin); + isecOn(this, other, this.uMax); + isecOn(other, this, other.uMin); + isecOn(other, this, other.uMax); + + curveIntersect( + this.impl, other.impl, + this.boundary(), other.boundary(), + CURVE_CACHING_TESSELLATOR, CURVE_CACHING_TESSELLATOR + ).forEach(i => add(i)); + isecs.forEach(i => { i.p0 = pt(i.p0); i.p1 = pt(i.p1); }); isecs = isecs.filter(({u0, u1}) => { let collinearFactor = Math.abs(this.tangentAtParam(u0).dot(other.tangentAtParam(u1))); - return !math.areEqual(collinearFactor, 1, tol); + return !math.areEqual(collinearFactor, 1); }); return isecs; -} + } - static createByPoints(points, degeree) { - points = points.map(p => p.data()); - return new NurbsCurve(new verb.geom.NurbsCurve.byPoints(points, degeree)); + isInside(u) { + return u >= this.uMin && u <= this.uMax; + } + + invert() { + return new NurbsCurve(this.impl.invert()); } } +const CURVE_CACHING_TESSELLATOR = function(curve, min, max, tessTol, scale) { + return cache('tess', [min, max, tessTol, scale], curve, () => degree1OptTessellator(curve, min, max, tessTol, scale)); +}; + +function degree1OptTessellator(curve, min, max, tessTol, scale) { + if (curve.degree() === 1) { + return curve.degree1Tess().map(u => curve.point(u)); + } + return curveTess(curve, min, max, tessTol, scale); +} + NurbsCurve.createLinearNurbs = function(a, b) { - return new NurbsCurve(new verb.geom.Line(a.data(), b.data())); + return new NurbsCurve(new NurbsCurveImpl(new verb.geom.Line(a.data(), b.data()))); }; export class NurbsSurface extends Surface { @@ -179,6 +296,7 @@ export class NurbsSurface extends Surface { } normalInMiddle() { + //TODO: use domain! return this.normalUV(0.5, 0.5); } @@ -196,7 +314,7 @@ export class NurbsSurface extends Surface { wp.x *= -1; } wp.__3D = pt3d; - return wp; + return wp; } workingPointTo3D(wp) { @@ -210,7 +328,8 @@ export class NurbsSurface extends Surface { return wp.__3D; } - static isMirrored(surface) { + static isMirrored(surface) { + //TODO: use domain! let a = surface.point(0, 0); let b = surface.point(1, 0); let c = surface.point(1, 1); @@ -218,19 +337,23 @@ export class NurbsSurface extends Surface { } intersectSurfaceForSameClass(other, tol) { - const curves = impl.surfaceIntersect(this.data, other.data, tol); + let curves = ext.surfaceIntersect(this.data, other.data, tol); let inverted = this.inverted !== other.inverted; - return curves.map(curve => new NurbsCurve(inverted ? curve.reverse() : curve)); + if (inverted) { + curves = curves.map(curve => ext.curveInvert(curve)); + } + curves.forEach(curve => ext.normalizeCurveParametrizationIfNeeded(curve)) + return curves.map(curve => new NurbsCurve(new NurbsCurveImpl(newVerbCurve(curve)))); } - + invert() { return new NurbsSurface(this.verb, !this.inverted); } isoCurve(param, useV) { const data = verb.eval.Make.surfaceIsocurve(this.verb._data, param, useV); - const isoCurve = new verb.geom.NurbsCurve(data); - return new NurbsCurve(isoCurve); + const isoCurve = newVerbCurve(data); + return new NurbsCurve(new NurbsCurveImpl(isoCurve)); } isoCurveAlignU(param) { @@ -245,7 +368,28 @@ export class NurbsSurface extends Surface { NurbsSurface.WORKING_POINT_SCALE_FACTOR = 1000; NurbsSurface.WORKING_POINT_UNSCALE_FACTOR = 1 / NurbsSurface.WORKING_POINT_SCALE_FACTOR; +NurbsSurface.loft = function(curve1, curve2) { + return new NurbsSurface(verb.geom.NurbsSurface.byLoftingCurves([curve1.impl.verb, curve2.impl.verb], 1)); +}; + +function newVerbCurve(data) { + return new verb.geom.NurbsCurve(data); +} + function pt(data) { return new Point().set3(data); } +function cache(id, keys, obj, op) { + id = '__cache__:' + id + ':' + keys.join('/'); + if (!obj[id]) { + obj[id] = op(); + } + return obj[id]; +} + +const surTess = verb.eval.Tess.rationalSurfaceAdaptive; +verb.eval.Tess.rationalSurfaceAdaptive = function(surface, opts) { + const keys = [opts ? opts.maxDepth: 'undefined']; + return cache('tess', keys, surface, () => surTess(surface, opts)); +} diff --git a/web/app/brep/geom/tolerance.js b/web/app/brep/geom/tolerance.js index 17055bbd..d69638e8 100644 --- a/web/app/brep/geom/tolerance.js +++ b/web/app/brep/geom/tolerance.js @@ -1,3 +1,20 @@ -export const TOLERANCE = 1e-6; +import {areEqual} from "../../math/math"; + +export const TOLERANCE = 1e-5; export const TOLERANCE_SQ = TOLERANCE * TOLERANCE; -export const TOLERANCE_HALF = TOLERANCE * 0.5; + +export const EPSILON = 1e-12; +export const EPSILON_SQ = EPSILON * EPSILON; + + +//tolerance used for parametric domain which is between 0..1 +export const TOLERANCE_01 = TOLERANCE * 1e-2; +export const TOLERANCE_01_SQ = TOLERANCE * TOLERANCE; + +export function eqTol(a, b) { + return areEqual(a, b, TOLERANCE); +} + +export function eqEps(a, b) { + return areEqual(a, b, EPSILON); +} \ No newline at end of file diff --git a/web/app/brep/operations/boolean.js b/web/app/brep/operations/boolean.js index 687511c5..ad2c1b7e 100644 --- a/web/app/brep/operations/boolean.js +++ b/web/app/brep/operations/boolean.js @@ -13,7 +13,7 @@ import {TOLERANCE} from '../geom/tolerance'; const DEBUG = { OPERANDS_MODE: false, - LOOP_DETECTION: true, + LOOP_DETECTION: false, FACE_FACE_INTERSECTION: false, NOOP: () => {} }; @@ -372,6 +372,7 @@ function intersectFaces(shell1, shell2, operationType) { let curves = face1.surface.intersectSurface(face2.surface, TOLERANCE); for (let curve of curves) { + // __DEBUG__.AddCurve(curve); curve = fixCurveDirection(curve, face1.surface, face2.surface, operationType); const nodes = []; collectNodesOfIntersectionOfFace(curve, face1, nodes); diff --git a/web/app/brep/operations/evolve-face.js b/web/app/brep/operations/evolve-face.js index 167eef02..fe9fa20f 100644 --- a/web/app/brep/operations/evolve-face.js +++ b/web/app/brep/operations/evolve-face.js @@ -17,7 +17,7 @@ export function evolveFace(originFace, loops) { function createFaces(nestedLoop, level) { let surface; - __DEBUG__.AddPointPolygon(nestedLoop.workingPolygon) + // __DEBUG__.AddPointPolygon(nestedLoop.workingPolygon) if (nestedLoop.inverted) { surface = invertSurface(surface); } else { diff --git a/web/app/brep/topo/loop.js b/web/app/brep/topo/loop.js index e04b6558..002a4acd 100644 --- a/web/app/brep/topo/loop.js +++ b/web/app/brep/topo/loop.js @@ -35,14 +35,13 @@ export class Loop extends TopoObject { tess() { let out = []; for (let e of this.halfEdges) { - let curvePoints = e.edge.curve.verb.tessellate(100000); + let curvePoints = e.edge.curve.tessellate(); if (e.inverted) { curvePoints.reverse(); } curvePoints.pop(); for (let point of curvePoints) { - let p = Point.fromData(point); - out.push(p); + out.push(point); } } return out; diff --git a/web/app/math/math.js b/web/app/math/math.js index 676e17f7..d89f711d 100644 --- a/web/app/math/math.js +++ b/web/app/math/math.js @@ -53,10 +53,12 @@ export function areEqual(v1, v2, tolerance) { } export function areVectorsEqual(v1, v2, tolerance) { + //TODO: use tolerance_SQ and length_SQ return areEqual(distanceAB3(v1, v2), 0, tolerance); } export function areVectorsEqual3(v1, v2, tolerance) { + //TODO: use tolerance_SQ and length_SQ return areEqual(distance3(v1[0], v1[1], v1[2], v2[0], v2[1], v2[2]), 0, tolerance); } diff --git a/web/app/math/vec.js b/web/app/math/vec.js index 0ebc9feb..9900ba3b 100644 --- a/web/app/math/vec.js +++ b/web/app/math/vec.js @@ -1,10 +1,21 @@ -export function __mul(v, scalar, out) { - out[0] = scalar * v[0]; - out[1] = scalar * v[1]; - out[2] = scalar * v[2]; +function scalarOperand(v, out, func) { + for (let i = 0; i < v.length; i++) { + out[i] = func(v[i]); + } return out; } +function vectorOperand(v1, v2, out, func) { + for (let i = 0; i < v1.length; i++) { + out[i] = func(v1[i], v2[i]); + } + return out; +} + +export function __mul(v, scalar, out) { + return scalarOperand(v, out, x => x * scalar); +} + export function _mul(v, scalar) { return __mul(v, scalar, v); } @@ -14,10 +25,7 @@ export function mul(v, scalar) { } export function __div(v, scalar, out) { - out[0] = v[0] / scalar; - out[1] = v[1] / scalar; - out[2] = v[2] / scalar; - return out; + return scalarOperand(v, out, x => x / scalar); } export function _div(v, scalar) { @@ -30,10 +38,7 @@ export function div(v, scalar) { export function __add(v1, v2, out) { - out[0] = v1[0] + v2[0]; - out[1] = v1[1] + v2[1]; - out[2] = v1[2] + v2[2]; - return out; + return vectorOperand(v1, v2, out, (x1, x2) => x1 + x2); } export function _add(v1, v2) { @@ -45,10 +50,7 @@ export function add(v1, v2) { } export function __sub(v1, v2, out) { - out[0] = v1[0] - v2[0]; - out[1] = v1[1] - v2[1]; - out[2] = v1[2] - v2[2]; - return out; + return vectorOperand(v1, v2, out, (x1, x2) => x1 - x2); } export function _sub(v1, v2) { @@ -59,12 +61,8 @@ export function sub(v1, v2) { return __sub(v1, v2, []); } - export function __negate(v, out) { - out[0] = - v[0]; - out[1] = - v[1]; - out[2] = - v[2]; - return out; + return scalarOperand(v, out, x => -x); } export function _negate(v) { @@ -77,7 +75,11 @@ export function negate(v) { export function dot(v1, v2) { - return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]; + let sum = 0; + for (let i = 0; i < v1.length; i++) { + sum += v1[i] * v2[i]; + } + return sum; } export function __cross(v1, v2, out) { diff --git a/web/lib/verb.js b/web/lib/verb.js index 24bc42b1..67640744 100644 --- a/web/lib/verb.js +++ b/web/lib/verb.js @@ -4063,6 +4063,12 @@ verb_eval_Eval.rationalSurfaceNormal = function(surface,u,v) { return verb_core_Vec.cross(derivs[1][0],derivs[0][1]); }; verb_eval_Eval.rationalSurfaceDerivatives = function(surface,u,v,numDerivs) { + if (surface.knotsU[surface.knotsU.length - 1] === u) { + u -= 1e-8; + } + if (surface.knotsV[surface.knotsV.length - 1] === v) { + v -= 1e-8; + } if(numDerivs == null) numDerivs = 1; var ders = verb_eval_Eval.surfaceDerivatives(surface,u,v,numDerivs); var Aders = verb_eval_Eval.rational2d(ders);