sgbaird / latticeplane Goto Github PK

PlaneIntersection[\[ScriptCapitalP]_, \[ScriptCapitalU]_] := 
 Module[{\[ScriptCapitalP]1, \[ScriptCapitalB], \[ScriptCapitalM], \
\[ScriptCapitalP]crop, \[ScriptCapitalP]int, shapes, bmesh, d, order, 
   coords, TF},
  If[Head@\[ScriptCapitalP] === 
    InfinitePlane, \[ScriptCapitalP]1 = {\[ScriptCapitalP]}, \
\[ScriptCapitalP]1 = \[ScriptCapitalP]];
  \[ScriptCapitalB] = 
   RegionResize[BoundingRegion[\[ScriptCapitalU]], 
    Scaled@Rationalize[
      1.1(*100$MachineEpsilon*)]];(*for cropping the plane*)
  \[ScriptCapitalP]crop = 
   RegionIntersection[#, \[ScriptCapitalB]] & /@ \[ScriptCapitalP]1;
  \[ScriptCapitalP]crop = 
   CanonicalizePolygon /@ N[\[ScriptCapitalP]crop];
  \[ScriptCapitalP]int = 
   RegionIntersection[#, N@\[ScriptCapitalU]] & /@ 
    N[\[ScriptCapitalP]crop];
  shapes = OpenCascadeShape[#] & /@ \[ScriptCapitalP]int;
  bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[#] & /@ shapes;
  coords = #["Coordinates"] & /@ bmesh;
  TF = CoplanarPoints /@ coords;
  coords = 
   MapThread[
    If[#1, #2, 
      DimensionReduction[#2, 2][#2, "ReconstructedData"]] &, {TF, 
     coords}];
  {d, order} = (FindShortestTour[#] & /@ coords)\[Transpose];
  \[ScriptCapitalP]int = 
   MapThread[Polygon[#1[[#2[[;; -2]]]]] &, {coords, order}]
  ]

Message[CoplanarPoints::pts, $Failed["Coordinates"]] "\!\($Failed[\"Coordinates\"]\) should be a non-empty list of points."