I'm trying to understand a Monte-Carlo Laplace/Poisson PDE solver: http://www.cs.cmu.edu/~kmcrane/Projects/MonteCarloGeometryProcessing/paper.pdf

This method inspired by random-walks and ray-tracing has an advantage over FEM mesh methods in that it requires little to no mesh pre-processing. I do not expect great performance from Mathematica but that's okay - I'm just trying to play with concepts in the paper.

I want to perform a walk on spheres in Mathematica shown in this diagram in the paper:

walk on sphere


  1. Given a mesh such as mesh = DiscretizeGraphics@ExampleData[{"Geometry3D", "Cow"}] how do I get a random point $x_0$ on the interior of the mesh $\Omega$ ? I've tried RandomPoint but it only gives points on the surface.
  2. How do I find the closest point on the mesh surface $\partial\Omega$ to $x_0$ so I can then determine the radius $r_0$ of the first sphere?

Once I've got $x_0, r_0$ I can do RandomPoint on the sphere to generate $x_1$ and Nest this process until the sphere radius is below some threshold.

  • $\begingroup$ Replace DiscretizeGraphics by BoundaryDiscretizeGraphics and use RegionNearest. RegionNearest[mesh] generates a RegionNearestFunction that can (and should) be reused. $\endgroup$ – Henrik Schumacher May 26 '20 at 20:31
  • $\begingroup$ If I do that I get "The first argument BoundaryDiscretizeGraphics[] is expected to be parameter-free" when calling RandomPoint $\endgroup$ – flinty May 26 '20 at 20:34
  • 2
    $\begingroup$ Ah, the cow is not watertight. You may try ExampleData[{"Geometry3D", "Triceratops"}, "BoundaryMeshRegion"] instead (it is already a BoundaryMeshRegion has `RegionDimension equal to 3). $\endgroup$ – Henrik Schumacher May 26 '20 at 20:39
  • $\begingroup$ With the above and rnf=RegionNearest[mesh], the point x0=RandomPoint[mesh] is correctly on the interior, but the nearest point I expected to lie on the surface rnf[x0] is actually the same point x0 $\endgroup$ – flinty May 26 '20 at 21:02
  • 2
    $\begingroup$ @flinty If Henrik doesn't write it up as an answer, then perhaps you could. As you probably know, self-answers are encouraged! $\endgroup$ – MarcoB May 26 '20 at 21:45

Thanks to Henrik Schumacher I got this to work with NestWhileList and visualized the walk on spheres. The other outputs of nextPoint besides RandomPoint are for visualization and the termination of the NestWhileList when the radius is small enough:

mesh = ExampleData[{"Geometry3D", "Triceratops"}, "BoundaryMeshRegion"];
rnf = RegionNearest@RegionBoundary@mesh;

nextPoint[p_] := Block[{r = EuclideanDistance[rnf[p], p]},
 {RandomPoint[Sphere[p, r]], p, r}]

walk = NestWhileList[
   nextPoint[#[[1]]] &, {RandomPoint[mesh], {}, ∞}, #[[3]] > 10^-3 &];

Graphics3D[{Opacity[.0], mesh,
  Opacity[1], Red, Thick, Line[walk[[All, 1]]],
  Blue, Opacity[0.04], 
  Sphere[#[[2]], #[[3]]] & /@ Rest[walk]}]

walk on spheres triceratops mesh

And with simple modifications it works for 2D, here demonstrated with a random polygon:

reg = RandomPolygon[12];
rnf = RegionNearest@RegionBoundary@reg;

nextPoint[p_] := 
 Block[{r = EuclideanDistance[rnf[p], p]}, {RandomPoint[Circle[p, r]], p, r}]

walk = NestWhileList[
   nextPoint[#[[1]]] &, {RandomPoint[reg], {}, ∞}, #[[3]] > 10^-5 &];

Graphics[{Opacity[.1], reg, Opacity[1], Red, Thick, 
  Line[walk[[All, 1]]], Blue, Opacity[0.04], 
  Disk[#[[2]], #[[3]]] & /@ Rest[walk]}]

walk on circles random polygon 2d

  • 3
    $\begingroup$ Love the paper, thanks for bringing it up in a Mathematica context! Would be awesome to see the whole poisson/laplace-solver working in Mathematica in case you are implementing it. $\endgroup$ – Thies Heidecke May 26 '20 at 22:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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