/**
# Distance field computation from a 3D model
The goal is to build a distance field representation of a 3D
[CAD](https://en.wikipedia.org/wiki/Computer-aided_design) model. */
#include "grid/octree.h"
#include "utils.h"
#include "distance.h"
#include "fractions.h"
#include "view.h"
int main()
{
/**
We get the 3D model from the [Large Geometric Models
Archive](http://www.cc.gatech.edu/projects/large_models/). This is
in [PLY format](https://en.wikipedia.org/wiki/PLY_%28file_format%29)
which we need to convert to
[STL](https://en.wikipedia.org/wiki/STL_%28file_format%29) using
[meshlab](http://www.meshlab.net/). */
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);
while (adapt_wavelet ({d}, (double[]){5e-4*L0}, 10).nf);
/**
We display an isosurface of the distance function coloured with the
level of refinement. */
view (fov = 15.65, quat = {-0.52,0.31,0.38,-0.7},
tx = -0.045, ty = 0.015, width = 640, height = 480, bg = {1,1,1});
isosurface ("d", 0, color = "level", min = 5, max = 10);
save ("horse.png");
/**
We also compute the 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 phi[];
foreach_vertex()
phi[] = (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.;
boundary ({phi});
scalar f[];
face vector s[];
fractions (phi, f, s);
/**
Finally we display the surface reconstructed from volume fractions. */
clear();
draw_vof ("f", "s");
draw_vof ("f", "s", edges = true, lw = 0.5);
save ("vof.png");
}
/**
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.](distance/horse.png)
![Reconstructed VOF surface.](distance/vof.png)
## See also
* [Computation of a levelset field from a contour](/src/test/basilisk.c)
*/