diff --git a/web/app/brep/geom/curves/basicCurve.js b/web/app/brep/geom/curves/basicCurve.js deleted file mode 100644 index e69de29b..00000000 diff --git a/web/app/brep/geom/curves/bezierCubic.js b/web/app/brep/geom/curves/bezierCubic.js new file mode 100644 index 00000000..79bd78a2 --- /dev/null +++ b/web/app/brep/geom/curves/bezierCubic.js @@ -0,0 +1,36 @@ +import * as vec from '../../../math/vec'; + +export function cubicBezierPoint(p0, p1, p2, p3, t) { + const mt = 1 - t; + const mt2 = mt * mt; + const mt3 = mt2 * mt; + const t2 = t * t; + const t3 = t2 * t; + + return vec.polynomial( + [mt3 , 3 * mt2 * t, 3 * mt * t2, t3], + [p0 , p1 , p2 , p3]); +} + +export function cubicBezierDer1(p0, p1, p2, p3, t) { + const mt = 1 - t; + const mt2 = mt * mt; + const t2 = t * t; + + return vec.polynomial( + [3 * mt2 , 6 * mt * t , 3 * t2], + [vec.sub(p1, p0), vec.sub(p2, p1), vec.sub(p3, p2)]); +} + +export function cubicBezierDer2(p0, p1, p2, p3, t) { + return vec.polynomial( + [ + 6 * (1 - t), + 6 * t + ], + [ + vec.polynomial([1, -2, 1], [p2, p1, p0]), + vec.polynomial([1, -2, 1], [p3, p2, p1]) + ]); +} + diff --git a/web/app/brep/geom/curves/brepCurve.js b/web/app/brep/geom/curves/brepCurve.js index ee42c996..9b18f48b 100644 --- a/web/app/brep/geom/curves/brepCurve.js +++ b/web/app/brep/geom/curves/brepCurve.js @@ -12,18 +12,10 @@ import cache from "../impl/cache"; export default class 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; + [uMin, uMax] = this.impl.domain(); + this.uMin = uMin; + this.uMax = uMax; } translate(vector) { @@ -176,7 +168,7 @@ const CURVE_CACHING_TESSELLATOR = function(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 curve.knots().map(u => curve.point(u)); } return curveTess(curve, min, max, tessTol, scale); } diff --git a/web/app/brep/geom/curves/cubicHermiteIntepolation.js b/web/app/brep/geom/curves/cubicHermiteIntepolation.js new file mode 100644 index 00000000..a8119e2f --- /dev/null +++ b/web/app/brep/geom/curves/cubicHermiteIntepolation.js @@ -0,0 +1,74 @@ +import * as vec from '../../../math/vec'; +import {cubicBezierDer1, cubicBezierDer2, cubicBezierPoint} from './bezierCubic'; + +export default function CubicHermiteInterpolation(points, tangents) { + let n = points.length; + let knots = new Array(n).fill().map((e,i) => i); + let beziers = []; + for (let i = 1; i < n; i++) { + let p0 = points[i - 1]; + let p3 = points[i]; + let tangent1 = tangents[i - 1]; + let tangent2 = tangents[i]; + let length = vec.length(vec.sub(p3, p0)) * 0.5; + let p1 = vec.add(p0, vec.mul(tangent1, length)); + let p2 = vec.sub(p3, vec.mul(tangent2, length)); + beziers.push({p0, p1, p2, p3}); + } + + function evalPatch(p0, p1, p2, p3, u, num) { + switch (num) { + case 0: return cubicBezierPoint(p0, p1, p2, p3, u); + case 1: return cubicBezierDer1(p0, p1, p2, p3, u); + case 2: return cubicBezierDer2(p0, p1, p2, p3, u); + default: throw 'illegal derivative order for cubic bezier curve'; + } + } + + function localizeParam(u) { + let pieceIndex; + if (u >= n - 1) { + pieceIndex = beziers.length - 1; + u = 1; + } else { + pieceIndex = Math.floor(u); + u = u % 1; + } + return [pieceIndex, u]; + } + + function evaluate(u, num) { + let [pieceIndex, uL] = localizeParam(u); + let {p0, p1, p2, p3} = beziers[pieceIndex]; + let out = []; + for (let i = 0; i <= num; ++i) { + out.push(evalPatch(p0, p1, p2, p3, uL, i)); + } + return out; + } + + function point(u, num) { + let [pieceIndex, uL] = localizeParam(u); + let {p0, p1, p2, p3} = beziers[pieceIndex]; + return cubicBezierPoint(p0, p1, p2, p3, uL); + } + + function param(point) { + throw 'unsupported'; + } + + function transform(point) { + throw 'unsupported'; + } + + Object.assign(this, { + domain: () => [0, n - 1], + degree: () => 3, + eval: evaluate, + point, + param, + transform, + knots: () => knots, + // invert + }); +} diff --git a/web/app/brep/geom/curves/intersectionCurve.js b/web/app/brep/geom/curves/intersectionCurve.js new file mode 100644 index 00000000..794c46c6 --- /dev/null +++ b/web/app/brep/geom/curves/intersectionCurve.js @@ -0,0 +1,90 @@ +import {TOLERANCE} from '../tolerance'; +import {surfaceClosestParam} from '../impl/nurbs-ext'; +import * as vec from '../../../math/vec'; +import CubicHermiteInterpolation from './cubicHermiteIntepolation'; + +export class IntersectionCurve { + + constructor(exactPoints, surfaceA, surfaceB) { + + let tangents = []; + for (let pt of exactPoints) { + let tangent = curveSSTangent(pt, surfaceA, surfaceB); + tangents.push(vec._normalize(tangent)); + } + + this.surfaceA = surfaceA; + this.surfaceB = surfaceB; + this.approx = new CubicHermiteInterpolation(exactPoints, tangents); + + this.exactify = (pt) => { + let uvA = surfaceClosestParam(surfaceA, pt); + let uvB = surfaceClosestParam(surfaceB, pt); + return verb.eval.Intersect.surfacesAtPointWithEstimate(surfaceA,surfaceB,uvA,uvB,TOLERANCE).point; + } + } + + domain() { + return this.approx.domain(); + } + + eval(u, num) { + let pt = this.point(u); + + // let [uA, vA] = surfaceClosestParam(this.surfaceA, pt); + // let [uB, vB] = surfaceClosestParam(this.surfaceB, pt); + + let out = []; + let approxEval = this.approx.eval(u, num); + for (let i = 0; i < num + 1; ++i) { + if (i === 0) { + out.push(pt); + } else { + // let nA = vec.cross(dA[i][0], dA[0][i]); + // let nB = vec.cross(dB[i][0], dB[0][i]); + // out.push(vec.cross(nA, nB)); + + out[i] = approxEval[i]; + } + } + return out; + } + + point(u) { + let pt = this.approx.point(u); + return this.exactify(pt); + } + + param(point) { + let pointOnCurve = this.exactify(point); + return this.approx.param(pointOnCurve); + } + + knots() { + return this.approx.knots(); + } + + degree() { + return undefined; + } + + transform(tr) { + throw 'unsupported;' + } + + invert() { + throw 'unsupported;' + } +} + +function curveSSTangent(pt, surfaceA, surfaceB) { + let [uA, vA] = surfaceClosestParam(surfaceA, pt); + let [uB, vB] = surfaceClosestParam(surfaceB, pt); + + let dA = verb.eval.Eval.rationalSurfaceDerivatives(surfaceA, uA, vA, 1); + let dB = verb.eval.Eval.rationalSurfaceDerivatives(surfaceB, uB, vB, 1); + + let nA = vec.cross(dA[1][0], dA[0][1]); + let nB = vec.cross(dB[1][0], dB[0][1]); + return vec.cross(nA, nB); +} \ No newline at end of file diff --git a/web/app/brep/geom/curves/invertedCurve.js b/web/app/brep/geom/curves/invertedCurve.js new file mode 100644 index 00000000..0fef9f13 --- /dev/null +++ b/web/app/brep/geom/curves/invertedCurve.js @@ -0,0 +1,56 @@ +import * as vec from '../../../math/vec'; + +export default class InvertedCurve { + + constructor(curve) { + this.curve = curve; + let [uMin, uMax] = this.curve.domain(); + this.uMin = uMin; + this.uMax = uMax; + } + + wrapParam(u) { + return this.uMax - (u - this.uMin); + } + + domain() { + return this.curve.domain(); + } + + degree() { + return this.curve.degree(); + } + + transform(tr) { + return new InvertedCurve(this.curve.transform(tr)); + } + + point(u) { + return this.curve.point(this.wrapParam(u)); + } + + param(point) { + return this.wrapParam(this.curve.param(point)); + } + + eval(u, num) { + let res = this.curve.eval(this.wrapParam(u), num); + if (res.length > 1) { + vec._negate(res[1]) + } + return eval; + + } + + knots() { + return this.curve.knots(); + } + + invert() { + return this.curve; + } + + split(u) { + return this.curve.split(this.wrapParam(u)).map(c => new InvertedCurve(c)).reverse(); + } +} diff --git a/web/app/brep/geom/curves/nurbsCurve.js b/web/app/brep/geom/curves/nurbsCurve.js index 6eb61d05..4e97ca1e 100644 --- a/web/app/brep/geom/curves/nurbsCurve.js +++ b/web/app/brep/geom/curves/nurbsCurve.js @@ -1,6 +1,6 @@ import * as ext from '../impl/nurbs-ext'; import {newVerbCurve} from "../impl/nurbs"; - +import {distinctKnots} from '../impl/nurbs-ext'; export default class NurbsCurve { @@ -13,10 +13,6 @@ export default class NurbsCurve { return ext.curveDomain(this.data); } - degree1Tess() { - return ext.distinctKnots(this.data); - } - degree() { return this.data.degree; } @@ -37,8 +33,8 @@ export default class NurbsCurve { return verb.eval.Eval.rationalCurveDerivatives( this.data, u, num ); } - optimalSplits() { - return this.data.knots.length - 1; + knots() { + return distinctKnots(this.data); } invert() { diff --git a/web/app/brep/geom/curves/parametricCurve.js b/web/app/brep/geom/curves/parametricCurve.js index 3c86fea9..4bbdb306 100644 --- a/web/app/brep/geom/curves/parametricCurve.js +++ b/web/app/brep/geom/curves/parametricCurve.js @@ -5,19 +5,15 @@ interface ParametricCurve { degree(): number; - degree1Tess(): number[][]; - eval(u: number, num: number): number[]; point(param: number): number[]; param(point: number[]): number; - transform(tr); + transform(tr); - optimalSplits(): number; - - normalizeParametrization(); + knots(): number[]; invert(); } diff --git a/web/app/brep/geom/impl/curve/curve-tess.js b/web/app/brep/geom/impl/curve/curve-tess.js index be6a19cf..235684d7 100644 --- a/web/app/brep/geom/impl/curve/curve-tess.js +++ b/web/app/brep/geom/impl/curve/curve-tess.js @@ -6,19 +6,15 @@ export default function curveTess(curve, min, max, tessTol, scale) { 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 knots = curve.knots(); let splits = [min]; - for (let i = 1; i < nSplits; ++i) { - splits.push(min + i * splitStep); + for (let split of knots) { + if (split > min && split < max) { + splits.push(split); + } } splits.push(max); diff --git a/web/app/brep/geom/impl/intersectionCurve.js b/web/app/brep/geom/impl/intersectionCurve.js deleted file mode 100644 index d3fa8064..00000000 --- a/web/app/brep/geom/impl/intersectionCurve.js +++ /dev/null @@ -1,69 +0,0 @@ -import {TOLERANCE} from '../tolerance'; -import {curveClosestParam, curveDomain, curvePoint, surfaceClosestParam} from './nurbs-ext'; -import * as vec from "../../../math/vec"; - -export class IntersectionCurve { - - constructor(approxPolyline, surfaceA, surfaceB) { - this.surfaceA = surfaceA; - this.surfaceB = surfaceB; - this.approxPolyline = approxPolyline; - this.exactify = (pt) => { - let uvA = surfaceClosestParam(surfaceA, pt); - let uvB = surfaceClosestParam(surfaceB, pt); - return verb.eval.Intersect.surfacesAtPointWithEstimate(surfaceA,surfaceB,uvA,uvB,TOLERANCE).point; - } - } - - domain() { - return curveDomain(this.approxPolyline); - } - - eval(u, num) { - let pt = this.point(u); - - let [uA, vA] = surfaceClosestParam(this.surfaceA, pt); - let [uB, vB] = surfaceClosestParam(this.surfaceB, pt); - - let dA = verb.eval.Eval.rationalSurfaceDerivatives(this.surfaceA, uA, vA, num); - let dB = verb.eval.Eval.rationalSurfaceDerivatives(this.surfaceB, uB, vB, num); - - let out = []; - for (let i = 0; i < num + 1; ++i) { - if (i === 0) { - out.push(pt); - } else { - let nA = vec.cross(dA[i][0], dA[0][i]); - let nB = vec.cross(dB[i][0], dB[0][i]); - out.push(vec.cross(nA, nB)); - } - } - return out; - } - - point(u) { - let pt = curvePoint(this.approxPolyline, u); - return this.exactify(pt); - } - - param(point) { - let pointOnCurve = this.exactify(point); - return curveClosestParam(this.approxPolyline, pointOnCurve); - } - - optimalSplits() { - return this.approxPolyline.knots.length - 3; - } - - degree() { - return undefined; - } - - transform(tr) { - throw 'unsupported;' - } - - invert() { - throw 'unsupported;' - } -} \ No newline at end of file diff --git a/web/app/brep/geom/impl/nurbs-ext.js b/web/app/brep/geom/impl/nurbs-ext.js index fab6ad99..0dfd5a80 100644 --- a/web/app/brep/geom/impl/nurbs-ext.js +++ b/web/app/brep/geom/impl/nurbs-ext.js @@ -83,6 +83,8 @@ export function curveClosestParam(curve, point) { return verb.eval.Analyze.rationalCurveClosestParam(curve, point); } +export const surfaceClosestParam = verb.eval.Analyze.rationalSurfaceClosestParam; + export function surfaceIntersect(surface0, surface1) { const tess0 = verb.eval.Tess.rationalSurfaceAdaptive(surface0); const tess1 = verb.eval.Tess.rationalSurfaceAdaptive(surface1); @@ -135,7 +137,7 @@ export function surfaceIntersect(surface0, surface1) { return nurbses; } -function meshesIntersect(mesh0,mesh1, TOLERANCE, TOLERANCE_SQ, TOLERANCE_01) { +export function meshesIntersect(mesh0,mesh1, TOLERANCE, TOLERANCE_SQ, TOLERANCE_01) { let bbtree0 = new verb.core.LazyMeshBoundingBoxTree(mesh0); let bbtree1 = new verb.core.LazyMeshBoundingBoxTree(mesh1); let bbints = verb.eval.Intersect.boundingBoxTrees(bbtree0,bbtree1,TOLERANCE); diff --git a/web/app/brep/geom/intersection/surfaceSurface.js b/web/app/brep/geom/intersection/surfaceSurface.js new file mode 100644 index 00000000..4bbb3961 --- /dev/null +++ b/web/app/brep/geom/intersection/surfaceSurface.js @@ -0,0 +1,78 @@ +import {TOLERANCE, TOLERANCE_01, TOLERANCE_SQ} from '../tolerance'; +import {curveDomain, curvePoint, meshesIntersect, surfaceMaxDegree} from '../impl/nurbs-ext'; +import {IntersectionCurve} from '../curves/intersectionCurve'; +import * as vec from '../../../math/vec'; + +export function surfaceIntersect(surfaceA, surfaceB) { + const tessA = verb.eval.Tess.rationalSurfaceAdaptive(surfaceA); + const tessB = verb.eval.Tess.rationalSurfaceAdaptive(surfaceB); + + 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(surfaceA, tessA); + fixTessNaNPoitns(surfaceB, tessB); + + const resApprox = meshesIntersect(tessA, tessB, TOLERANCE, TOLERANCE_SQ, TOLERANCE_01); + const exactPls = resApprox.map(pl => pl.map(inter => + verb.eval.Intersect.surfacesAtPointWithEstimate(surfaceA, surfaceB, inter.uv0, inter.uv1, TOLERANCE) + )); + + let curves = []; + for (let pl of exactPls) { + let points = pl.map(ip => ip.point); + + points = removeSuperfluousPoints(points, 30*30); //5*5 + // points.forEach(__DEBUG__.AddPoint3); + curves.push(new IntersectionCurve(points, surfaceA, surfaceB)) + } + return curves; +} + +//it's Douglas–Peucker and it's not well suited here +function removeSuperfluousPoints(points, tolSq) { + + let out = []; + let stack = [[0, points.length - 1]]; + out.push(points[0]); + while (stack.length !== 0) { + let stackItem = stack.pop(); + if (!Array.isArray(stackItem)) { + out.push(points[stackItem]); + continue; + } + let [from, to] = stackItem; + let maxDistSq = tolSq; + let hero = -1; + let v = vec._normalize(vec.sub(points[to], points[from])); + + for (let i = from + 1; i < to; i ++) { + let proj = vec.dot(v, vec.sub(points[i], points[from])); + let vA = vec.add(points[from], vec.mul(v, proj)); + let vX = vec.sub(points[i], vA); + let perpDistSq = vec.lengthSq(vX); + if (perpDistSq > maxDistSq) { + hero = i; + maxDistSq = perpDistSq; + } + } + if (hero !== -1) { + if (to - hero > 1) { + stack.push([hero, to]); + } + stack.push(hero); + if (hero - from > 1) { + stack.push([from, hero]); + } + } + } + out.push(points[points.length - 1]); + return out; +} \ No newline at end of file diff --git a/web/app/cad/debugPlugin.js b/web/app/cad/debugPlugin.js index 3db3ace9..4b380b33 100644 --- a/web/app/cad/debugPlugin.js +++ b/web/app/cad/debugPlugin.js @@ -142,8 +142,8 @@ function addGlobalDebugActions({viewer, cadScene, cadRegistry}) { debugVolumeGroup.add(mesh); viewer.render(); }, - AddCurve: (curve, color) => { - __DEBUG__.AddPolyLine( curve.tessellate(), color); + AddCurve: (curve, color, scale) => { + __DEBUG__.AddPolyLine( curve.tessellate(undefined, scale), color); }, AddVerbCurve: (curve, color) => { __DEBUG__.AddPolyLine(curve.tessellate().map(p => new Vector().set3(p)), color); diff --git a/web/app/cad/sandbox.js b/web/app/cad/sandbox.js index 31d7e3b5..60ed3eb4 100644 --- a/web/app/cad/sandbox.js +++ b/web/app/cad/sandbox.js @@ -4,9 +4,8 @@ import BrepBuilder from '../brep/brep-builder' import * as BREPPrimitives from '../brep/brep-primitives' import {NurbsSurface} from "../brep/geom/impl/nurbs"; import BrepCurve from '../brep/geom/curves/brepCurve'; -import {surfaceIntersect} from '../brep/geom/impl/nurbs-ext'; import NurbsCurve from "../brep/geom/curves/nurbsCurve"; - +import {surfaceIntersect} from '../brep/geom/intersection/surfaceSurface'; export function runSandbox({bus, services: { viewer, cadScene, cadRegistry, tpi, tpi: {addShellOnScene} }}) { @@ -187,16 +186,22 @@ export function runSandbox({bus, services: { viewer, cadScene, cadRegistry, tpi, let surfaceB = cadRegistry.findFace('1:4').brepFace.surface; - let curve = surfaceIntersect(surfaceA.data, surfaceB.data)[0]; + let curves = surfaceIntersect(surfaceA.data, surfaceB.data); // curve.approxPolyline. - - console.dir(curve); - __DEBUG__.AddCurve(curve); + + for (let ic of curves) { + + let curve = new BrepCurve(ic); + + console.dir(curve); + __DEBUG__.AddCurve(curve, 0xffffff, 10); + __DEBUG__.HideSolids(); + } } cylinderAndPlaneIntersect(); - + // curvesIntersect(); } diff --git a/web/app/math/vec.js b/web/app/math/vec.js index c07e2414..f3e63edf 100644 --- a/web/app/math/vec.js +++ b/web/app/math/vec.js @@ -118,3 +118,15 @@ export function lengthSq(v) { export function length(v) { return Math.sqrt(lengthSq(v)); } + +export function polynomial(coefs, vectors) { + let out = []; + out.length = vectors[0].length; + out.fill(0); + for (let i = 0; i < vectors.length; i++) { + for (let j = 0; j < out.length; j++) { + out[j] += vectors[i][j] * coefs[i]; + } + } + return out; +}