3d raycast for face filtering

web/app/3d/debug.js

@@ -3,6 +3,7 @@ import {nurbsToThreeGeom, triangulateToThree} from './scene/brep-scene-object'
 import {createSolidMaterial} from './scene/scene-object'
 import DPR from '../utils/dpr'
 import Vector from "../math/vector";
+import {NurbsCurve} from "../brep/geom/impl/nurbs";
 
 export const DEBUG = true;
 
@@ -135,6 +136,9 @@ function addGlobalDebugActions(app) {
     AddCurve: (curve, color) => {
       __DEBUG__.AddPolyLine( curve.tessellate(), color);
     },
+    AddVerbCurve: (curve, color) => {
+      __DEBUG__.AddPolyLine(curve.tessellate().map(p => new Vector().set3(p)), color);
+    },
     AddNurbsCorners: (nurbs) => {
       __DEBUG__.AddPoint(nurbs.point(0, 0), 0xff0000);
       __DEBUG__.AddPoint(nurbs.point(1, 0), 0x00ff00);

web/app/brep/geom/impl/nurbs-ext.js

@@ -108,7 +108,7 @@ export function surfaceIntersect(surface0, surface1) {
   });
 
 
-  //temporary workaround
+  //TODO: temporary workaround
   return exactPls.map(pl => verb.eval.Make.polyline(pl.map(ip => ip.point)));
   
   return exactPls.map(function(x) {

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

@@ -363,6 +363,10 @@ export class NurbsSurface extends Surface {
   isoCurveAlignV(param) {
     return this.isoCurve(param, false);
   }
+  
+  intersectWithCurve(curve) {
+    return verb.geom.Intersect.curveAndSurface(curve.impl.verb, this.verb, TOLERANCE).map(({uv}) => uv);
+  }
 }
 
 NurbsSurface.WORKING_POINT_SCALE_FACTOR = 1000;

web/app/brep/geom/tolerance.js

@@ -15,6 +15,10 @@ export function eqTol(a, b) {
   return areEqual(a, b, TOLERANCE);
 }
 
+export function eqSqTol(a, b) {
+  return areEqual(a, b, TOLERANCE_SQ);
+}
+
 export function eqEps(a, b) {
   return areEqual(a, b, EPSILON);
 }

web/app/brep/operations/boolean.js

@@ -4,12 +4,14 @@ import {Loop} from '../topo/loop';
 import {Shell} from '../topo/shell';
 import {Vertex} from '../topo/vertex';
 import {evolveFace} from './evolve-face'
+import PIP from '../../3d/tess/pip';
 import * as math from '../../math/math';
-import {eqEps, eqTol, TOLERANCE, ueq, veq} from '../geom/tolerance';
+import {eqEps, eqTol, eqSqTol, TOLERANCE, ueq, veq} from '../geom/tolerance';
+import {Ray} from "../utils/ray";
 
 const DEBUG = {
-  OPERANDS_MODE: true,
-  LOOP_DETECTION: true,
+  OPERANDS_MODE: false,
+  LOOP_DETECTION: false,
   FACE_FACE_INTERSECTION: false,
   NOOP: () => {}
 };
@@ -49,6 +51,7 @@ export function invert( shell ) {
       loop.link();
     }
   }
+  shell.data.inverted = !shell.data.inverted;
   BREPValidator.validateToConsole(shell);
 }
 
@@ -90,7 +93,7 @@ export function BooleanAlgorithm( shell1, shell2, type ) {
     loopsToFaces(faceData.face, faceData.detectedLoops, faces);
   }
 
-  faces = filterFaces(faces);
+  faces = filterFaces(faces, shell1, shell2, type !== TYPE.UNION);
   
   
   const result = new Shell();
@@ -177,7 +180,6 @@ export function mergeVertices(shell1, shell2) {
   }
 }
 
-
 function filterFacesByInvalidEnclose(faces) {
 
   function encloses(f1, f2, testee, overEdge) {
@@ -229,10 +231,123 @@ function filterFacesByInvalidEnclose(faces) {
       }
     }
   }
-  return faces.filter(f => !invalidFaces.has(f));
+  return Array.from(faces).filter(f => !invalidFaces.has(f));
 }
 
-function filterFaces(faces) {
+function isPointInsideSolid(pt, initDir, solid) {
+  let ray = new Ray(pt, initDir, 3000);
+  const pertrubStep = 5;
+  for (let i = 0; i < 360; i+=pertrubStep) {
+    let res = rayCastSolidImpl(ray, solid);
+    if (res !== null) {
+      return res;    
+    }
+    ray.pertrub(pertrubStep);
+  }
+  return false; 
+}
+
+function rayCastSolidImpl(ray, solid) {
+  __DEBUG__.AddCurve(ray.curve, 0xffffff);
+  let closestDistanceSq = -1;
+  let inside = false;
+  let hitEdge = false;
+
+  let edgeDistancesSq = [];
+  for (let e of solid.edges) {
+    let points = e.curve.intersectCurve(ray.curve, TOLERANCE);
+    for (let {p0} of points) {
+      edgeDistancesSq.push(ray.pt.distanceToSquared(p0));
+    }  
+  }
+
+
+
+  for (let face of solid.faces) {
+    __DEBUG__.AddFace(face, 0xffff00);
+    let pip = face.data[MY].pip;
+    let uvs = face.surface.intersectWithCurve(ray.curve);
+    
+    for (let uv of uvs) {
+      let normal = face.surface.normalUV(uv[0], uv[1]);
+      let dotPr = normal.dot(ray.dir);
+      if (eqTol(dotPr, 0)) {
+        continue;
+      }
+      let pt = face.surface.point(uv[0], uv[1]);
+      let wpt = face.surface.createWorkingPoint(uv, pt); 
+      let pipClass = pip(wpt);
+      if (pipClass.inside) {
+        let distSq = ray.pt.distanceToSquared(pt);
+         if (closestDistanceSq === -1 || distSq < closestDistanceSq) {
+          hitEdge = false; 
+          for (let edgeDistSq of edgeDistancesSq) {
+            if (eqSqTol(edgeDistSq, distSq)) {
+              hitEdge = true;
+            }    
+          }
+          closestDistanceSq = distSq;
+          inside = dotPr > 0;
+        }
+      }
+    }    
+  }
+
+  if (hitEdge) {
+    return null;
+  }
+
+  if (solid.data.inverted) {
+    inside = !inside;
+  }
+  return inside;
+}
+
+function guessPointOnFace(face) {
+  //TODO:
+  let {pip, workingPolygon: [poly]} = createPIPForFace(face);
+  let [a, b] = poly;
+  let ab = b.minus(a); 
+  let len = ab.length();
+  let unit = ab.normalize();
+  let res = a.plus(unit.multiply(len * 0.5));
+  let offVec = unit.multiply(100); //???
+  res.x += - offVec.y;
+  res.y +=   offVec.x;
+
+  if (!pip(res).inside) {
+    throw 'bummer';
+  }
+  return face.surface.workingPointTo3D(res);
+}
+
+function filterByRayCast(faces, a, b, isIntersection) {
+  
+  let result = [];
+  for (let face of faces) {
+    __DEBUG__.Clear();
+    __DEBUG__.AddFace(face, 0x00ff00);
+    let pt = guessPointOnFace(face);
+    let normal = face.surface.normal(pt);
+
+    let insideA = face.data.__origin.shell === a || isPointInsideSolid(pt, normal, a);
+    let insideB = face.data.__origin.shell === b || isPointInsideSolid(pt, normal, b);
+    if (isIntersection) {
+      if (insideA && insideB) {
+        result.push(face);
+      }
+    } else {
+      if (insideA || insideB) {
+        result.push(face);
+      }
+    }
+  }
+  return result;
+}
+
+function filterFaces(faces, a, b, isIntersection) {
+
+  return filterByRayCast(faces, a, b, isIntersection);
 
   function isFaceContainNewEdge(face) {
     for (let e of face.edges) {
@@ -246,8 +361,8 @@ function filterFaces(faces) {
   const validFaces = new Set(faces);
   const result = new Set();
   for (let face of faces) {
-    // __DEBUG__.Clear();
-    // __DEBUG__.AddFace(face);
+    __DEBUG__.Clear();
+    __DEBUG__.AddFace(face);
     traverseFaces(face, validFaces, (it) => {
       if (result.has(it) || isFaceContainNewEdge(it)) {
         result.add(face);
@@ -255,7 +370,7 @@ function filterFaces(faces) {
       }
     });
   }
-  return filterFacesByInvalidEnclose(result);
+  return result;//filterFacesByInvalidEnclose(result);
 }
 
 function traverseFaces(face, validFaces, callback) {
@@ -573,7 +688,7 @@ function intersectCurveWithEdge(curve, edge, result) {
       vertex = vertexFactory.create(point.p0);
     }
 
-    __DEBUG__.AddVertex(vertex);
+    // __DEBUG__.AddVertex(vertex);
 
     result.push(new Node(vertex, edge, curve, u1));
   }
@@ -748,6 +863,16 @@ class SolveData {
   }
 }
 
+
+function createPIPForFace(face) {
+  let workingPolygon = [face.outerLoop, ...face.innerLoops].map(loop => loop.tess().map(pt => face.surface.workingPoint(pt)));
+  let [inner, ...outers] = workingPolygon;
+  return {
+    pip: PIP(inner, outers),
+    workingPolygon
+  }
+}
+
 class FaceSolveData {
   constructor(face) {
     this.face = face;
@@ -755,6 +880,7 @@ class FaceSolveData {
     face.innerLoops.push(this.loopOfNew);
     this.vertexToEdge = new Map();
     this.graphEdges = [];
+    Object.assign(this, createPIPForFace(face))
   }
 
   initGraph() {

web/app/brep/operations/evolve-face.js

@@ -27,6 +27,7 @@ export function evolveFace(originFace, loops) {
     const loop = nestedLoop.loop;
     const newFace = new Face(surface);
     Object.assign(newFace.data, originFace.data);
+    newFace.data.__origin = originFace;
     newFace.outerLoop = loop;
     loop.face = newFace;
     out.push(newFace);
@@ -59,6 +60,7 @@ export function evolveFace(originFace, loops) {
 function getNestedLoops(face, brepLoops) {
   function NestedLoop(loop) {
     this.loop = loop;
+    //FIXME
     this.workingPolygon = loop.asPolygon().map(p => face.surface.workingPoint(p));
     this.inverted = !isCCW(this.workingPolygon);
     this.pip = PIP(this.workingPolygon);

web/app/brep/topo/edge.js

@@ -79,4 +79,9 @@ class HalfEdge extends TopoObject {
     }
     he.manifoldHolder = this;
   }
+  
+  clone() {
+    let clone
+    Object.assign(clone.data, this.data);
+  }
 }

web/app/brep/topo/face.js

@@ -13,6 +13,22 @@ export class Face extends TopoObject {
     this.defineIterable('loops', () => loopsGenerator(this));
     this.defineIterable('edges', () => halfEdgesGenerator(this))
   }
+  
+  clone() {
+    function cloneLoop(source, dest) {
+      source.halfEdges.forEach(edge => dest.halfEdges.push(edge.clone()));
+      Object.assign(dest.data, source.data);
+      return dest;
+    }
+    
+    let clone = new Face(this.surface);
+    cloneLoop(this.outerLoop, clone.outerLoop);
+    this.innerLoops.forEach(loop => clone.innerLoops.push(cloneLoop(loop, new Loop(clone))));
+    Object.assign(clone.data, this.data);
+    return clone;
+  }
+  
+  
 }
 
 export function* loopsGenerator(face) {

web/app/brep/topo/shell.js

@@ -18,6 +18,15 @@ export class Shell extends TopoObject {
       e.halfEdge2.vertexA.edges.add(e.halfEdge2);
     }
   }
+  
+  clone() {
+    let clone = new Shell();
+    this.faces.forEach(face => clone.faces.push(face.clone));
+    clone.faces.forEach(face => face.shell = clone);
+    Object.assign(clone.data, this.data);
+    return clone;
+  }
+  
 }
 
 export function* verticesGenerator(shell) {

web/app/brep/utils/ray.js

@@ -0,0 +1,21 @@
+import pertrub from './vector-petrub';
+import {NurbsCurve, NurbsCurveImpl} from '../geom/impl/nurbs';
+
+
+export class Ray {
+  constructor(pt, dir, reachableDistance) {
+    this.pt = pt;
+    this.dir = dir;
+    this.reachableDistance = reachableDistance;
+    this.updateCurve();
+  }
+
+  updateCurve() {
+    this.curve = NurbsCurve.createLinearNurbs(this.pt, this.pt.plus(this.dir.multiply(this.reachableDistance)));
+  }
+
+  pertrub(angleStep) {
+    this.dir.set3(pertrub(this.dir.data(), angleStep));
+    this.updateCurve();
+  }
+} 

web/app/brep/utils/vector-petrub.js

@@ -0,0 +1,22 @@
+export default function pertrub([x, y, z], angleStep) {
+
+  //convert to spherical coordinate system
+  let r = Math.sqrt(x*x + y*y + z*z);
+  let teta = Math.acos(z / r);
+  let phi = Math.atan2(y, x);
+
+  phi = anglePertrub(phi, angleStep);
+  teta = anglePertrub(teta, angleStep);
+  return [
+    r * Math.sin(teta) * Math.cos(phi),
+    r * Math.sin(teta) * Math.sin(phi),
+    r * Math.cos(teta),
+  ];
+
+}
+
+const _2PI = 2 * Math.PI;
+
+function anglePertrub(angle, angleStep) {
+	return (angle + 0.75 * Math.PI + angleStep/_2PI) % _2PI;
+}