CRUST algorithm

Curve reconstruction algorithm by Nina Amenta, Marshall Bern, and David Eppstein (1998). Reference: section 5.7 in Devadoss and O’Rourke text.

Example 1

Start with some points. Suppose that these points are sampled from an (unknown) curve.

input=;

Extract the list of Point objects from the input graphic:

pointList=InputForm[input][[1,1]]

Extract the coordinates from the Point objects:

points=Part[#,1] &/@pointList

Compute the Voronoi diagram of the points. We will set a bounding box of [0,1]×[0,1] to make it easier to discard the boundary later.

vMesh=VoronoiMesh[points,{{0,1},{0,1}}];
Show[vMesh,Graphics[pointList]]

Get the coordinates of the Voronoi vertices. Unfortunately, Mathematica also returns coordinates along the boundary of the bounding box.

vVerts=MeshCoordinates[vMesh]

Select only the vertices that are not part of the bounding box. In other words, select only the vertices that don’t have a coordinate of 0 or 1.

vVerts=Select[vVerts,#[[1]]≠0&&#[[1]]≠1&&#[[2]]≠0&&#[[2]]≠1&]

Draw the Voronoi diagram, the Voronoi vertices, and the original points:

Show[vMesh,Graphics[{Black,Point[vVerts],Red,Point[points]}]]

Compute the Delaunay triangulation of the set containing the Voronoi vertices and the original points.

dMesh=DelaunayMesh[Join[points,vVerts]];
Show[dMesh,Graphics[{Black,Point[vVerts],Red,Point[points]}]]

Extract the edges from the Delaunay triangulation. Note these are indexed according to the order in which the points occurred in the combined list.

dEdges=MeshCells[dMesh,1]

Select only those edges with both endpoints in the original point set. These points have the smallest indexes, from 1 to Length[points].

curve=Select[dEdges,#[[1,1]]≤Length[points]&&#[[1,2]]≤Length[points]&]

Now draw the reconstructed curve. The simplest way to do this is to use MeshRegion.

MeshRegion[points,curve]

Example 2

input=;

Extract the list of Point objects from the input graphic:

pointList=InputForm[input][[1,1]]

Extract the coordinates from the Point objects:

points=Part[#,1] &/@pointList

Compute the Voronoi diagram of the points. We will set a bounding box of [-1,2]×[-1,2] to make it easier to discard the boundary later.

vMesh=VoronoiMesh[points,{{-1,2},{-1,2}}];
Show[vMesh,Graphics[pointList]]

Get the coordinates of the Voronoi vertices. Unfortunately, Mathematica also returns coordinates along the boundary of the bounding box.

vVerts=MeshCoordinates[vMesh]

Select only the vertices that are not part of the bounding box. In other words, select only the vertices that don’t have a coordinate of -1 or 2.

vVerts=Select[vVerts,#[[1]]≠-1&&#[[1]]≠2&&#[[2]]≠-1&&#[[2]]≠2&]

Draw the Voronoi diagram, the Voronoi vertices, and the original points:

Show[vMesh,Graphics[{Black,Point[vVerts],Red,Point[points]}]]

Compute the Delaunay triangulation of the set containing the Voronoi vertices and the original points.

dMesh=DelaunayMesh[Join[points,vVerts]];
Show[dMesh,Graphics[{Black,Point[vVerts],Red,Point[points]}]]

Extract the edges from the Delaunay triangulation. Note these are indexed according to the order in which the points occurred in the combined list.

dEdges=MeshCells[dMesh,1]

Select only those edges with both endpoints in the original point set. These points have the smallest indexes, from 1 to Length[points].

curve=Select[dEdges,#[[1,1]]≤Length[points]&&#[[1,2]]≤Length[points]&]

Now draw the reconstructed curve. The simplest way to do this is to use MeshRegion.

MeshRegion[points,curve]