# Distance field computation from a 3D model

The goal is to build a distance field representation of a 3D CAD model.

``````#include "grid/octree.h"
#include "utils.h"
#include "distance.h"
#include "fractions.h"

int main()
{``````

We get the 3D model from the Large Geometric Models Archive. This is in PLY format which we need to convert to STL using meshlab.

``````  system ("test -f distance.stl || "
"(wget http://www-static.cc.gatech.edu/data_files/large_models/horse.ply.gz && "
"gunzip -f horse.ply.gz && "
"meshlabserver -i horse.ply -o distance.stl)");``````

We read the STL file, compute the bounding box of the model and set the domain center and size using this bounding box.

``````  coord * p = input_stl (fopen ("distance.stl", "r"));
coord min, max;
bounding_box (p, &min, &max);
double maxl = -HUGE;
foreach_dimension()
if (max.x - min.x > maxl)
maxl = max.x - min.x;

init_grid (8);
size (1.2*maxl);
origin ((max.x + min.x)/2. - L0/2,
(max.y + min.y)/2. - L0/2,
(max.z + min.z)/2. - L0/2);``````

We initialize the distance field on the coarse initial mesh and refine it adaptively until the threshold error (on distance) is reached.

``````  scalar d[];
distance (d, p);

We display an isosurface of the distance function coloured with the level of refinement.

``````  FILE * fp =
popen ("gfsview-batch3D distance.gfv | convert ppm:- horse.png", "w");
output_gfs (fp);
fprintf (fp, "Save stdout { format = PPM }\n");``````

We also compute volume and surface fractions from the distance field. We first construct a vertex field interpolated from the centered field and then call the appropriate VOF functions.

``````  vertex scalar φ[];
foreach_vertex()
φ[] = (d[] + d[-1] + d[0,-1] + d[-1,-1] +
d[0,0,-1] + d[-1,0,-1] + d[0,-1,-1] + d[-1,-1,-1])/8.;
scalar f[];
face vector s[];
fractions (φ, f, s);``````

Finally we display the surface reconstructed from volume fractions.

``````  FILE * fp1 =
popen ("gfsview-batch3D distance.vof.gfv | "
"convert ppm:- -resize 640x480 vof.png", "w");
output_gfs (fp1);
fprintf (fp1, "Save stdout { format = PPM width = 1280 height = 960 }\n");

#if 0
FILE * fp2 = popen ("gfsview3D -s distance.gfv", "w");
output_gfs (fp2);
#endif
}``````

Note that the “tail” of the horse is an artefact due to an inconsistency in the surface mesh, which is self-intersecting near this point.