Cubic Hermite spline interpolation curve

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])
+    ]);
+}
+

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);
 }

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
+  });
+}

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);
+}

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();
+  }
+}

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() {

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();
 }

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);
 

web/app/brep/geom/impl/intersectionCurve.js

@@ -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;'
-  }
-}

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);

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;
+}

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);

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();
 }
 
 

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;
+}