ray cast and vector perturbation

web/app/brep/operations/boolean.js

@@ -8,11 +8,13 @@ import PIP from '../../3d/tess/pip';
 import * as math from '../../math/math';
 import {eqEps, eqTol, eqSqTol, TOLERANCE, ueq, veq} from '../geom/tolerance';
 import {Ray} from "../utils/ray";
+import pickPointInside2dPolygon from "../utils/pickPointInPolygon";
 
 const DEBUG = {
   OPERANDS_MODE: false,
   LOOP_DETECTION: false,
   FACE_FACE_INTERSECTION: false,
+  RAY_CAST: false,
   NOOP: () => {}
 };
 
@@ -147,6 +149,7 @@ function detectLoops(face) {
       if (!candidates) {
         break;
       }
+      candidates = candidates.filter(c => c.vertexB !== edge.vertexA || !isSameEdge(c, edge));
       edge = findMaxTurningLeft(edge, candidates, surface);
       if (seen.has(edge)) {
         break;
@@ -236,19 +239,20 @@ function filterFacesByInvalidEnclose(faces) {
 
 function isPointInsideSolid(pt, initDir, solid) {
   let ray = new Ray(pt, initDir, 3000);
-  const pertrubStep = 5;
-  for (let i = 0; i < 360; i+=pertrubStep) {
+  for (let i = 0; i < 1; ++i) {
     let res = rayCastSolidImpl(ray, solid);
     if (res !== null) {
       return res;    
     }
-    ray.pertrub(pertrubStep);
+    ray.pertrub();
   }
   return false; 
 }
 
 function rayCastSolidImpl(ray, solid) {
-  __DEBUG__.AddCurve(ray.curve, 0xffffff);
+  if (DEBUG.RAY_CAST) {
+    __DEBUG__.AddCurve(ray.curve, 0xffffff);  
+  }
   let closestDistanceSq = -1;
   let inside = false;
   let hitEdge = false;
@@ -261,10 +265,10 @@ function rayCastSolidImpl(ray, solid) {
     }  
   }
 
-
-
   for (let face of solid.faces) {
-    __DEBUG__.AddFace(face, 0xffff00);
+    if (DEBUG.RAY_CAST) {
+      __DEBUG__.AddFace(face, 0xffff00);
+    }
     let pip = face.data[MY].pip;
     let uvs = face.surface.intersectWithCurve(ray.curve);
     
@@ -303,31 +307,28 @@ function rayCastSolidImpl(ray, solid) {
   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';
+function pickPointOnFace(face) {
+  let wp = pickPointInside2dPolygon(face.createWorkingPolygon());
+  if (wp === null) {
+    return null;
   }
-  return face.surface.workingPointTo3D(res);
+  return face.surface.workingPointTo3D(wp);
 }
 
 function filterByRayCast(faces, a, b, isIntersection) {
   
   let result = [];
   for (let face of faces) {
-    __DEBUG__.Clear();
-    __DEBUG__.AddFace(face, 0x00ff00);
-    let pt = guessPointOnFace(face);
+    if (DEBUG.RAY_CAST) {
+      __DEBUG__.Clear();
+      __DEBUG__.AddFace(face, 0x00ff00);
+    }
+
+    let pt = pickPointOnFace(face);
+    if (pt === null) {
+      continue;
+    }
+    
     let normal = face.surface.normal(pt);
 
     let insideA = face.data.__origin.shell === a || isPointInsideSolid(pt, normal, a);
@@ -601,7 +602,9 @@ function intersectFaces(shell1, shell2, operationType) {
       let curves = face1.surface.intersectSurface(face2.surface);
 
       for (let curve of curves) {
-        // __DEBUG__.AddCurve(curve);
+        if (DEBUG.FACE_FACE_INTERSECTION) {
+          __DEBUG__.AddCurve(curve);
+        }
         curve = fixCurveDirection(curve, face1.surface, face2.surface, operationType);
         const nodes = [];
         collectNodesOfIntersectionOfFace(curve, face1, nodes);
@@ -865,7 +868,7 @@ class SolveData {
 
 
 function createPIPForFace(face) {
-  let workingPolygon = [face.outerLoop, ...face.innerLoops].map(loop => loop.tess().map(pt => face.surface.workingPoint(pt)));
+  let workingPolygon = face.createWorkingPolygon();
   let [inner, ...outers] = workingPolygon;
   return {
     pip: PIP(inner, outers),
@@ -888,13 +891,15 @@ class FaceSolveData {
     for (let he of this.face.edges) {
       this.addToGraph(he);
     }
-    this.removeOppositeEdges();
   }
 
   addToGraph(he) {
     // __DEBUG__.Clear();
     // __DEBUG__.AddFace(he.loop.face);
     // __DEBUG__.AddHalfEdge(he, 0xffffff);
+    if (this.isNewOppositeEdge(he)) {
+      return;
+    }
     let list = this.vertexToEdge.get(he.vertexA);
     if (!list) {
       list = [];
@@ -911,17 +916,24 @@ class FaceSolveData {
     this.graphEdges.push(he);
   }
 
+  isNewOppositeEdge(e1) {
+    if (!isNew(e1)) {
+      return false;
+    }
+    let others = this.vertexToEdge.get(e1.vertexB);
+    if (others) {
+      for (let e2 of others) {
+        if (e1.vertexA === e2.vertexB && isSameEdge(e1, e2)) {
+          return true;
+        }
+      }
+    }
+    return false;    
+  }
+  
   removeOppositeEdges() {
     let toRemove = new Set();
     for (let e1 of this.graphEdges) {
-      let others = this.vertexToEdge.get(e1.vertexB);
-      for (let e2 of others) {
-        if (e1 === e2) continue;
-        if (e1.vertexA === e2.vertexB && isSameEdge(e1, e2)) {
-          toRemove.add(e1);
-          toRemove.add(e2);
-        }
-      }
     }
     for (let e of toRemove) {
       removeFromListInMap(this.vertexToEdge, e.vertexA, e);

web/app/brep/topo/edge.js

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

web/app/brep/topo/face.js

@@ -13,22 +13,10 @@ 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;
+
+  createWorkingPolygon() {
+    return [this.outerLoop, ...this.innerLoops].map(loop => loop.tess().map(pt => this.surface.workingPoint(pt)));
   }
-  
-  
 }
 
 export function* loopsGenerator(face) {

web/app/brep/topo/shell.js

@@ -18,15 +18,6 @@ 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/pickPointInPolygon.js

@@ -0,0 +1,47 @@
+import libtess from 'libtess'
+import {area} from "../../math/math";
+
+export default function pickPointInside2dPolygon(polygon) {
+  function vertexCallback(data, tr) {
+    tr.points[tr.counter] = data;
+    tr.counter ++;
+    if (tr.counter === 3) {
+      let trArea = Math.abs(area(tr.points));
+      if (trArea > tr.bestArea) {
+        tr.bestArea = trArea;
+        tr.bestTr = Array.from(tr.points); 
+      }
+      tr.counter = 0;
+    }
+  }
+
+  const tessy = new libtess.GluTesselator();
+  tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_VERTEX_DATA, vertexCallback);
+  tessy.gluTessNormal(0, 0, 1);
+  const tracker = {
+    points: [],
+    bestTr: null,
+    bestArea: -1,
+    counter: 0
+  };
+  tessy.gluTessBeginPolygon(tracker);
+
+  for (let path of polygon) {
+    tessy.gluTessBeginContour();
+    for (let p of path) {
+      tessy.gluTessVertex([p.x, p.y, 0], p);
+    }
+    tessy.gluTessEndContour();
+  }
+  tessy.gluTessEndPolygon();
+
+  if (tracker.bestTr === null) {
+    return null;
+  }
+  
+  let center = tracker.bestTr[0].copy();
+  center._plus(tracker.bestTr[1]);
+  center._plus(tracker.bestTr[2]);
+  center._divide(3);
+  return center;
+}

web/app/brep/utils/ray.js

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

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

@@ -1,22 +1,27 @@
-export default function pertrub([x, y, z], angleStep) {
+export default function pertrub([x, y, z]) {
+
+  let s = x + y + z;
+  x = pertrubFloat(x + 3 + s);
+  y = pertrubFloat(y + 5 + s);
+  z = pertrubFloat(z + 7 + s);
 
-  //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),
+    x/r,
+    y/r,
+    z/r
   ];
-
 }
 
-const _2PI = 2 * Math.PI;
-
-function anglePertrub(angle, angleStep) {
-	return (angle + 0.75 * Math.PI + angleStep/_2PI) % _2PI;
+function pertrubFloat(x) {
+  return xorshift32(Math.round(x * 1e8)) ;
 }
+
+function xorshift32(x) {
+  x ^= x << 13;
+  x ^= x >> 17;
+  x ^= x << 5;
+  return x;
+}
+

web/app/math/vector.js

@@ -36,6 +36,14 @@ Vector.prototype._multiply = function(scalar) {
   return this.set(this.x * scalar, this.y * scalar, this.z * scalar);
 };
 
+Vector.prototype.divide = function(scalar) {
+  return new Vector(this.x / scalar, this.y / scalar, this.z / scalar);
+};
+
+Vector.prototype._divide = function(scalar) {
+  return this.set(this.x / scalar, this.y / scalar, this.z / scalar);
+};
+
 Vector.prototype.dot = function(vector) {
   return this.x * vector.x + this.y * vector.y + this.z * vector.z;
 };