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 && "
	  "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;
    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);
  while (adapt_wavelet ({d}, (double[]){5e-4*L0}, 10).nf);

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 φ[];
    φ[] = (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);

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.

Isosurface of the distance function coloured with level of refinement.

Isosurface of the distance function coloured with level of refinement.

Reconstructed VOF surface.

Reconstructed VOF surface.

See also