8
$\begingroup$

I want to use geometric shapes in Mathematica to build complex shapes and use my raytracing algorithm on it. I have a working example where we can get the intersections from a combination of a Cone[] and Cuboid[], e.g 

shape1 = Cone[];
shape2 = Cuboid[];
(* add shapes in this list to make a more complicated shape *) 
shapes = {shape1, shape2};

(* this constains the shapes so the shape is considered as a whole *) 
constraints[shapes__] := 
 And[## & @@ (Not /@ 
      Through[(RegionMember[RegionIntersection@##] & @@@ 
          Subsets[{shapes}, {2}])@#]), 
   RegionMember[RegionUnion @@ (RegionBoundary /@ {shapes})]@#] &

direction = {-0.2, -0.2, -1};
point = {0.5, 0.5, 1.5};

line = HalfLine[{point, point + direction}];

intersections[l_, s__] := 
 NSolve[# ∈ l && constraints[s][#], #] &@({x, y, z}[[;; RegionEmbeddingDimension[l]]])

(* find intersection *) 
intersection = intersections[line, ##] & @@ shapes;

points = Point[{x, y, z}] /. intersection;

Graphics3D[{{Opacity[0.2], shapes}, line, {Red, points}}, 
 PlotRange -> {{-1, 1}, {-1, 1}, {-2, 2}}, Axes -> True]

This works well, and we get the external intersections as we expect.

enter image description here

Now, let us try to take the difference between two shapes, modelling something like 

square = Cuboid[];
ball = Ball[{0, 0, 1}, 1];
Region[RegionDifference[square, ball]]

enter image description here

shapes = {RegionDifference[square, ball]};
direction = {0, 0, -1};
point = {0.5, 0.5, 5};
line = HalfLine[{point, point + direction}];
intersection = intersections[line, ##] & @@ shapes

Does not work, with an error that the constraints are "not a quantified system of equations and inequalities"...even though the constraints look fine 

constraints[shapes]
(* (##1 &) @@ 
   Not /@ Through[
     Apply[RegionMember[RegionIntersection[##1]] &, 
       Subsets[{{BooleanRegion[#1 && ! #2 &, {Cuboid[{0, 0, 0}], 
            Ball[{0, 0, 1}, 1]}]}}, {2}], {1}][#1]] && 
  RegionMember[
    RegionUnion @@ 
     RegionBoundary /@ {{BooleanRegion[#1 && ! #2 &, {Cuboid[{0, 0, 
            0}], Ball[{0, 0, 1}, 1]}]}}][#1] & *)
Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ Just to give you starting ideas: have you seen Quilez's treatment? $\endgroup$
    – J. M.'s missing motivation ♦
    Commented May 18, 2020 at 11:55
  • $\begingroup$ No, but this looks very helpful indeed. $\endgroup$
    – Tomi
    Commented May 18, 2020 at 11:59
  • 2
    $\begingroup$ Tim Laska's answer shows some great 3D raytracing using complicated shapes. There might be something there that can help if you haven't already seen it. $\endgroup$
    – MassDefect
    Commented May 18, 2020 at 18:16

3 Answers 3

Reset to default
9
$\begingroup$

This is not a direct answer to your question, but an alternate approach. You could create a list of primitives and a build function that contains the Computational Solid Geometry (CSG).

square = Cuboid[];
ball = Ball[{0, 0, 1}, 1];
buildList = {square, ball};
(* Constraints *)
buildFn = ¬ #2 ∧ #1 &;
reg = Region[
   Style[BooleanRegion[buildFn, buildList], Opacity[0.5], Green]];
direction = {0, 0, -1};
point = {0.5, 0.5, 5};
line = HalfLine[{point, point + direction}];
rint = Region[RegionIntersection[reg, line], 
  BaseStyle -> {Blue, Thick}]; 
intpoints = Point[Transpose@RegionBounds@rint];
Show[reg, rint, Graphics3D[{PointSize[Large], Red, intpoints}]]

Difference

Here is how it would look for the initial case:

shape1 = Cone[];
shape2 = Cuboid[];
buildList = {shape1, shape2};
(* Constraints *)
buildFn = #2 || #1 &;
reg = Region[
   Style[BooleanRegion[buildFn, buildList], Opacity[0.5], Green]];
direction = {-0.2, -0.2, -1};
point = {0.5, 0.5, 1.5};
line = HalfLine[{point, point + direction}];
rint = Region[RegionIntersection[reg, line], 
  BaseStyle -> {Blue, Thick}]; intpoints = 
 Point[Transpose@RegionBounds@rint];
Show[reg, rint, Graphics3D[{PointSize[Large], Red, intpoints}], 
 PlotRange -> All]

Initial Case

Update to Increase Speed

@Tomi mentioned in the comments that speed is a concern. As addressed in my answer to the MSE question Why is Ray Tracing Slow? I created a solver that that used the fast region functions RegionDistanceand RegionNormal to solve a 1000 multiple bounce ray traces in 3D geometry including geometry produced by a commercial CAD package. I will adapt that approach to look at the bouncing of single ray.

Set up the Geometry

The OpenCascadeLink does a pretty good job at constructing geometry that snaps to features while keeping the triangle count down. The following workflow will create the the initial Box-Cone geometry.

Needs["OpenCascadeLink`"]
Needs["NDSolve`FEM`"]
pp = Polygon[{{0, 0, 0}, {0, 0, 1}, {1, 0, 1}}];
shape = OpenCascadeShape[pp];
axis = {{0, 0, 0}, {0, 0, 1}};
sweep = OpenCascadeShapeRotationalSweep[shape, axis];
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[sweep];
Show[Graphics3D[{{Red, pp}, {Blue, Thick, Arrow[axis]}}], 
 bmesh["Wireframe"], Boxed -> False]
cu = OpenCascadeShape[Cuboid[{0, 0, 0}, {1, 1, 1}]];
union = OpenCascadeShapeUnion[cu, sweep];
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[union];
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]
mrd = MeshRegion[bmesh, PlotTheme -> "Lines"]

Solve a Single Ray Trace

The following workflow solves for a single ray trace. Each bounce will cause the ray to attenuate the representative sphere size by 10%. This solves and plots quickly.

(* Set up Region Operators on Differenced Geometry *)
rdf = RegionDistance[mrd];
rnf = RegionNearest[mrd];
(* Setup and run simulation *)
(* Time Increment *)
dt = 0.01;
(* Collision Margin *)
margin = (1 + dt) dt;
(* Conditional Particle Advancer *)
advance[r_, x_, v_, c_] := 
 Block[{xnew = x + dt v}, {rdf[xnew], xnew, v, c}] /; r > margin
advance[r_, x_, v_, c_] := 
 Block[{xnew = x , vnew = v, normal = Normalize[x - rnf[x]]},
   vnew = Normalize[v - 2 v.normal normal];
   xnew += dt vnew;
   {rdf[xnew], xnew, vnew, c + 1}] /; r <= margin
(* Starting Point for Emission *)
sp = {0, 0, 0.25};
nparticles = 1;
ntimesteps = 800;
tabres = Table[
   NestList[
    advance @@ # &, {rdf[sp], 
     sp, { Cos[2 Pi #[[1]]] Sin[Pi #[[2]]], 
        Sin[ Pi #[[2]]] Sin[2 Pi #[[1]]], Cos[ Pi #[[2]]]} &@
      First@RandomReal[1, {1, 2}], 0}, ntimesteps], {i, 1, 
    nparticles}];
epilog[i_] := {ColorData["Rainbow", (#4 - 1)/10], 
    Sphere[#2, 0.04 0.9^#4]} & @@@ tabres[[i]]
Graphics3D[{White, EdgeForm[Thin], Opacity[0.25], mrd, Opacity[1]}~
  Join~epilog[1], Boxed -> False, PlotRange -> RegionBounds[mrd], 
 ViewPoint -> {-1.7742436871276688`, 1.5459832360779067`, 
   2.431459473742817`}, 
 ViewVertical -> {0.052110700162003136`, -0.06948693625348555`, 
   0.9962208794332359`}]

Cone Box Single Ray Trace

A More Complex Case

The following produces a shape with concavity that could find rays that intersect but would be blocked by an intervening surface. Because the solver use a fine time increment, these intersections are not found because the collision of the intervening surface is detected.

pp = Polygon[{{0, 0, 0}, {0, 0, 1}, {1, 0, 1}}];
shape = OpenCascadeShape[pp];
axis = {{0, 0, 0}, {0, 0, 1}};
sweep = OpenCascadeShapeRotationalSweep[shape, axis];
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[sweep];
Show[Graphics3D[{{Red, pp}, {Blue, Thick, Arrow[axis]}}], 
 bmesh["Wireframe"], Boxed -> False]
cu = OpenCascadeShape[Cuboid[{0, 0, 0}, {1, 1, 1}]];
ball = OpenCascadeShape[Ball[{1/2, 1/2, 2.4}, 1.5]];
union = OpenCascadeShapeUnion[cu, sweep, ball];
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[union];
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]
mrd = MeshRegion[bmesh, PlotTheme -> "Lines"]
(* Set up Region Operators on Differenced Geometry *)
rdf = RegionDistance[mrd];
rnf = RegionNearest[mrd];
(* Setup and run simulation *)
(* Time Increment *)
dt = 0.01;
(* Collision Margin *)
margin = (1 + dt) dt;
(* Conditional Particle Advancer *)
advance[r_, x_, v_, c_] := 
 Block[{xnew = x + dt v}, {rdf[xnew], xnew, v, c}] /; r > margin
advance[r_, x_, v_, c_] := 
 Block[{xnew = x , vnew = v, normal = Normalize[x - rnf[x]]},
   vnew = Normalize[v - 2 v.normal normal];
   xnew += dt vnew;
   {rdf[xnew], xnew, vnew, c + 1}] /; r <= margin
(* Starting Point for Emission *)
sp = {0, 0, 0.5};
nparticles = 1;
ntimesteps = 1600;
(*tabres= Table[NestList[advance@@#&,{rdf[sp],sp,{ Cos[2 Pi #[[1]]] \
Sin[Pi #[[2]]],Sin[ Pi #[[2]]] Sin[2 Pi #[[1]]], Cos[ Pi \
#[[2]]]}&@First@RandomReal[1,{1,2}],0},ntimesteps],{i,1,nparticles}];*)


tabres = Table[
   NestList[
    advance @@ # &, {rdf[sp], 
     sp, { Cos[2 Pi #[[1]]] Sin[Pi #[[2]]], 
        Sin[ Pi #[[2]]] Sin[2 Pi #[[1]]], Cos[ Pi #[[2]]]} &@
      First@{{0.3788624698388783`, 0.8749177935911279`}}, 0}, 
    ntimesteps], {i, 1, nparticles}];
epilog[i_] := {ColorData["Rainbow", (#4 - 1)/12], 
    Sphere[#2, 0.04 0.9^#4]} & @@@ tabres[[i]]
Graphics3D[{White, EdgeForm[Thin], Opacity[0.25], mrd, Opacity[1]}~
  Join~epilog[1], Boxed -> False, PlotRange -> RegionBounds[mrd], 
 ViewPoint -> {-3.102894731729034`, -1.0062787100553268`, 
   0.8996929706836663`}, 
 ViewVertical -> {-0.34334064946409365`, -0.07403103185215265`, 
   0.93628874005217`}]

Concave Surface

Improve this answer
$\endgroup$
3
  • $\begingroup$ The disadvantage is that RegionIntersection[] is much slower than NSolve[]. Hence your particle advancer solution, which works very well - but I was hoping to be able to use the capabilities of NSolve[]. $\endgroup$
    – Tomi
    Commented May 19, 2020 at 15:21
  • $\begingroup$ I think I have identified the problem: posted a new question here: mathematica.stackexchange.com/questions/222199/… $\endgroup$
    – Tomi
    Commented May 19, 2020 at 16:22
  • $\begingroup$ @Tomi Thank you for the accept and your kind words. Just a word of caution about OpenCascade. Since it is a CAD package, you probably need to avoid non-manifold geometry (i.e., zero volume sections). Also, I may have discovered bug with the Cone primitive and that is why I revolved a triangle around the axis. I will contact support when I have time. $\endgroup$
    – Tim Laska
    Commented May 20, 2020 at 1:18
5
$\begingroup$

Tim Laska's solution is excellent. It is fast and accurate. However, for completeness, I have a solution for the NDSolve solution, where we can find the intersections instead of the (excellent) particle advancer (i.e. just jump between the intersections instead of advance).

By using the solution from here 

line = HalfLine[{0.5, 0.5, 2}, {0, 0, -1}]


intersection = 
 NSolve[{x, y, z} \[Element] line && 
   RegionMember[
     regionBoundary[RegionDifference[Cuboid[], Ball[]]]][{x, y, 
     z}], {x, y, z}]

regionBoundary[reg_?RegionQ] := 
 Module[{x, y, z}, 
  ImplicitRegion[
   CylindricalDecomposition[RegionMember[reg, {x, y, z}], {x, y, z}, 
    "Boundary"], {x, y, z}]]

Show[{Region[RegionDifference[Cuboid[], Ball[]]], 
  Region[Style[Point[{x, y, z}] /. intersection[[1]], Red]], 
  Region[Style[Point[{x, y, z}] /. intersection[[2]], Red]]}]

enter image description here

Intersections highlighted in red.

Improve this answer
$\endgroup$
2
$\begingroup$

This issue could be solved by "sphere tracing", it's similar to Tim's particle advancer" and reasonable fast which will converge in a very small amount of steps.

RayIntersect[ray_, region_, maxIteration_, maxRadius_, 
  radiusThreshold_] := Module[{ rnf, RegionMarcher, result},
  rnf = RegionNearest[region];
  RegionMarcher[orig_, radius_, step_] := 
   Block[{next, closestPoint, advOp, radiusCurrent},
    closestPoint = rnf[orig];
    advOp = 
     If[Dot[closestPoint - orig, dir /. ray] > 0, 
      Subtract, {a, b} |-> b - a];
    radiusCurrent =  Max[Abs[advOp @@ {closestPoint , orig}]] // N;
    next = orig + radiusCurrent * (dir /. ray);
    Sow[{orig, next, Abs[radiusCurrent], advOp, step + 1}];
    {next, Abs[radiusCurrent], step + 1}
    ] ;
  
  result = 
   NestWhile[RegionMarcher @@ # &, {origin /. ray, 1, 0}, 
    e |-> (e[[2]] < maxRadius && e[[2]] > radiusThreshold && 
       e[[3]] < maxIteration)];
  
  If[result[[2]] >= maxRadius || result[[3]] >= maxIteration, None, 
   result[[1]]]
  ]
RayIntersect[ray_, region_] := 
 RayIntersect[ray, region, 10, 256, 0.01]

And the result:

Graphics3D[{Opacity[.5], Blue, 
  DiscretizeRegion[RegionDifference[square, ball]], 
  Arrow[{point, point + direction}], Red, 
  Map[(e |-> { Sphere[e[[1]], e[[3]]]}), 
   Reap[RayIntersect[{origin -> point, dir -> direction}, 
        RegionDifference[square, ball]]] // Rest // First // First]}, 
 Axes -> True]

sphere tarcing

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.