src/distance.h

Computation of a distance field from a discretised curve or surface
- Usage
  - Examples
  - Tests

Computation of a distance field from a discretised curve or surface

The goal is to compute the (signed) minimal distance from a set of oriented segments (2d) or triangles (3D) to any point on the grid. This signed distance function is also dynamically updated when the mesh is refined or coarsened.

The scheme is robust and will give results even for inconsistent surface representations (i.e. surfaces with holes, manifold egdes, incompatible orientations etc…). Faces do not need to be properly connected i.e. the scheme works also for “triangle soups”.

We need a few utility function such as computation of the minimal distance between a point and a segment or a point and a triangle.

#include <stdint.h>
#include "PointTriangle.h"

The 3D triangulated surfaces are defined using the (binary) STL format which can be exported from most CAD modelling programs. The input_stl() function reads such a file a returns an array of triplets of vertex coordinates defining the triangles.

trace
coord * input_stl (FILE * fp)
{
  Array * a = array_new();
  char tag[6];

  if (fgets (tag, 6, fp) != tag) {
    fprintf (stderr, "input_stl(): error reading tag\n");
    exit(1);
  }
  rewind (fp);
  if (!strcmp (tag, "solid")) { /* ASCII file */
    fprintf (stderr, "input_stl(): ASCII STL not implemented yet "
	     "(use binary instead)\n");
    exit(1);
  }
  else { /* binary file */
    uint32_t nf;
    char header[80];
    unsigned i;

    if (fread (header, sizeof (char), 80, fp) != 80) {
      fprintf (stderr, "Input file is not a valid STL file\n"
	       "stdin: incomplete header\n");
      exit (1);
    }
    if (fread (&nf, sizeof (uint32_t), 1, fp) != 1) {
      fprintf (stderr, "Input file is not a valid STL file\n"
	       "stdin: missing number of facets\n");
      exit (1);
    }
    i = nf;
    while (i > 0) {
      float x = 0, y = 0, z = 0;
      dimensional ((x == Delta, y == Delta, z == Delta));
      unsigned j;
      uint16_t attbytecount;

      if (fread (&x, sizeof (float), 1, fp) != 1) {
	fprintf (stderr, "Input file is not a valid STL file\n"
	       "stdin: missing normal x-coordinate\n");
	exit (1);
      }
      if (fread (&y, sizeof (float), 1, fp) != 1) {
	fprintf (stderr, "Input file is not a valid STL file\n"
	       "stdin: missing normal y-coordinate\n");
	exit (1);
      }
      if (fread (&z, sizeof (float), 1, fp) != 1) {
	fprintf (stderr, "Input file is not a valid STL file\n"
	       "stdin: missing normal z-coordinate\n");
	exit (1);
      }

      for (j = 0; j < 3; j++) {
	if (fread (&x, sizeof (float), 1, fp) != 1) {
	  fprintf (stderr, "Input file is not a valid STL file\n"
		   "stdin: missing vertex x-coordinate\n");
	  exit (1);
	}
	if (fread (&y, sizeof (float), 1, fp) != 1) {
	  fprintf (stderr, "Input file is not a valid STL file\n"
		   "stdin: missing vertex y-coordinate\n");
	  exit (1);
	}
	if (fread (&z, sizeof (float), 1, fp) != 1) {
	  fprintf (stderr, "Input file is not a valid STL file\n"
		   "stdin: missing vertex z-coordinate\n");
	  exit (1);
	}
	coord p = {x,y,z};
	array_append (a, &p, sizeof(coord));
      }

      if (fread (&attbytecount, sizeof (uint16_t), 1, fp) != 1) {
	fprintf (stderr, "Input file is not a valid STL file\n"
	       "stdin: missing attribute byte count\n");
	exit (1);
      }
      i--;
    }
  }
  coord p = {nodata};
  array_append (a, &p, sizeof(coord));
  return (coord *) array_shrink (a);
}

In 2D dimensions, the file format is that used by gnuplot i.e. pairs of 2D vertex coordinates separated by blank lines. An easy way to create these files is to use a vector graphics editing program (e.g. inkscape) and save the curve as an .eps file, then convert it to gnuplot format using

pstoedit -f gnuplot -flat 0.1 file.eps file.gnu

The function below reads such a file and returns an array of pairs of coordinates.

trace
coord * input_xy (FILE * fp)
{
  Array * a = array_new();
  coord p = {0}, last, * la = NULL;
  while (!feof(fp)) {
    if (fscanf (fp, "%lf %lf", &p.x, &p.y) == 2) {
      if (la) {
	array_append (a, la, sizeof(coord));
	array_append (a, &p, sizeof(coord));
      }
      last = p, la = &last;
    }
    else {
      int c;
      while ((c = fgetc(fp)) != EOF && c != '\n');
      la = NULL;
    }
  }
  p.x = nodata;
  array_append (a, &p, sizeof(coord));
  return (coord *) array_shrink (a);
}

This function computes the coordinates of the bounding box of a set of segments or triangles.

void bounding_box (coord * p, coord * min, coord * max)
{
  foreach_dimension()
    (*min).x = HUGE, (*max).x = - HUGE;
  while (p->x != nodata) {
    foreach_dimension() {
      if ((*p).x < (*min).x)
	(*min).x = (*p).x;
      if ((*p).x > (*max).x)
	(*max).x = (*p).x;
    }
    p++;
  }
}

An extra field, holding a pointer to the elements (segments or triangles) intersecting the neighborhood of the cell, is associated with the distance function. The neighborhood is a sphere centered on the cell center and with a diameter 3\Delta.

attribute {
  scalar surface;
}

#define double_to_pointer(d) (*((void **) &(d)))
#define BSIZE 3. // if larger than 1, cells overlap

To compute the minimal distance and its sign, we need to store information on the closest (2 in 2D, 12 in 3D) elements.

typedef struct {
  double d2; // minimal distance
  coord * v; // the element
  int type;  // 0,1,2: vertex, 3: edge, 4: face
} closest_t;

#if dimension == 2
#  define ND 2
#else // dimension == 3
#  define ND 12
static double vertex_angle (coord * p, int type)
{
  if (type >= 3) // edge or face
    return pi;
  coord * v = p + type, * v1 = p + (type + 1)%3, * v2 = p + (type + 2)%3;
  coord vv1 = vecdiff (*v1, *v), vv2 = vecdiff(*v2,*v);
  return acos(vecdot(vv1,vv2)/sqrt(vecdot(vv1,vv1)*vecdot(vv2,vv2)));
}

static coord face_normal (coord * q, int type)
{
  coord ab = vecdiff(*(q+1),*q), ac = vecdiff(*(q+2),*q);
  coord n = vecdotproduct(ab,ac);
  double nn = sqrt(vecdot(n,n));
  assert (nn > 0.);
  nn = vertex_angle(q, type)/nn;
  foreach_dimension()
    n.x *= nn;
  return n;
}
#endif // dimension == 3

static void update_distance (Point point, coord ** i, scalar d)
{
  scalar surface = d.surface;
  Array * a = array_new();
  coord c = {x,y,z}, closest = {0};
  closest_t q[ND];
  for (int i = 0; i < ND; i++)
    q[i].d2 = HUGE;
  int nd = 0;
  double r2 = sq(BSIZE*Delta/2.);
  bool first = (level == 0);
  while (*i) {
    coord * p = *i;
#if dimension == 2
    coord r;
    double s, d2 = PointSegmentDistance (&c, p, p + 1, &r, &s);
#elif dimension == 3
    double s, t, d2 = PointTriangleDistance (&c, p, p + 1, p + 2, &s, &t);
#endif
    // keep pointers/distances/types of up to ND closest elements
    for (int i = 0; i < ND; i++)
      if (d2 < q[i].d2) {
	for (int j = ND - 1; j > i; j--)
	  q[j] = q[j-1];
	q[i].d2 = d2, q[i].v = p;
#if dimension == 2
	// vertices
	if (s == 0.)
	  q[i].type = 0;
	else if (s == 1.)
	  q[i].type = 1;
	else
	  // edge
	  q[i].type = 3;
	if (i == 0)
	  closest = r;
#elif dimension == 3
	// vertices
	if (s == 0. && t == 0.)
	  q[i].type = 0;
	else if (s == 1. && t == 0.)
	  q[i].type = 1;
	else if (s == 0. && t == 1.)
	  q[i].type = 2;
	else if (s == 0. || t == 0. || s + t == 1.)
	  // edge
	  q[i].type = 3;
	else
	  // face
	  q[i].type = 4;
	if (i == 0)
	  foreach_dimension()
	    closest.x = ((*q[0].v).x*(1. - s - t) + s*(*(q[0].v+1)).x +
			 t*(*(q[0].v+2)).x);
#endif // dimension == 3
	if (i >= nd)
	  nd = i + 1;
	break;
      }
    // add elements which are close enough to the local list
    if (d2 < r2 || first)
      array_append (a, &p, sizeof(coord *));
    first = false, i++;
  }
  if (a->len) {
    // set surface[] to list, ended with NULL
    coord * p = NULL;
    array_append (a, &p, sizeof(coord *));
    p = (coord *) array_shrink (a);
    assert (sizeof(real) >= sizeof(void *));
    memcpy (&surface[], &p, sizeof(void *));

    int orient;
#if dimension == 2
    if (q[0].type == 3)
      // edge
      orient = PointSegmentOrientation (&c, q[0].v, q[0].v + 1);
    else {
      // vertex
      if (nd == 1) // a single vertex, cannot find sign
	// get sign from parent
	orient = sign(bilinear (point, d));
      else { // two vertices
	orient = PointSegmentOrientation (&c, q[0].v, q[0].v + 1);
	if (orient != PointSegmentOrientation (&c, q[1].v, q[1].v + 1)) {
	  coord n = {0};
	  for (int i = 0; i < 2; i++) {
#if DEBUG	    
	    fprintf (stderr, "q %g %g %g %g\n",
		     x, y, q[i].v->x - x, q[i].v->y - y);
#endif
	    coord ab = vecdiff(*(q[i].v + 1),*q[i].v);
	    double nn = sqrt(vecdot(ab,ab));
	    assert (nn > 0.);
	    n.x -= ab.y/nn, n.y += ab.x/nn;
	  }
#if DEBUG
	  fprintf (stderr, "vertex %g %g %g %g %g %g %g %g\n",
		   x, y, closest.x - x, closest.y - y,
		   closest.x, closest.y, n.x, n.y);
#endif
	  coord diff = vecdiff(closest,c);
	  orient = sign(vecdot(n,diff));
	}
      }
    }
#elif dimension == 3
    if (q[0].type < 3) {
      coord n = face_normal (q[0].v, q[0].type), diff = vecdiff(closest,c);
      int nv = 1;
      orient = sign(vecdot(n,diff));
      for (int i = 1; i < nd && q[i].type < 3; i++)
	if (vecdist2 (closest, *(q[i].v + q[i].type)) < sq(1e-6)) {
	  coord n1 = face_normal (q[i].v, q[i].type);
	  foreach_dimension()
	    n.x += n1.x;
	  nv++;
	  if (orient > -2 && sign(vecdot(n1,diff)) != orient)
	    orient = -2;
	}
      if (nv < 3) // less than 3 vertices, cannot find sign
	// get sign from parent
	orient = sign(bilinear (point, d));
      else if (orient == -2) {
	// the vertices do not have the same orientation
	// get the proper orientation from the pseudo-normal n
	orient = sign(vecdot(n,diff));
#if DEBUG
	fprintf (stderr, "vertex %g %g %g %g %g %g %d %g %g %g %g %g %g %d\n",
		 x, y, z, closest.x - x, closest.y - y, closest.z - z, orient,
		 closest.x, closest.y, closest.z, n.x, n.y, n.z, nv);
#endif
      }
    }
    else if (q[0].type == 3) {
      // edge
      if (nd == 1 || q[1].type != 3) // a single edge, cannot find sign
	// get sign from parent
	orient = sign(bilinear (point, d));
      else { // two edges
	orient = PointTriangleOrientation (&c, q[0].v, q[0].v+1, q[0].v+2);
	if (orient !=
	    PointTriangleOrientation (&c, q[1].v, q[1].v+1, q[1].v+2)) {
	  coord n1 = face_normal (q[0].v, 3), n2 = face_normal (q[1].v, 3), n;
	  foreach_dimension()
	    n.x = n1.x + n2.x;
	  coord diff = vecdiff(closest,c);
	  orient = sign(vecdot(n,diff));
#if DEBUG
	  fprintf (stderr, "edge %g %g %g %g %g %g %d %g %g %g %g %g %g\n",
		   x, y, z, closest.x - x, closest.y - y, closest.z - z, orient,
		   closest.x, closest.y, closest.z, n.x, n.y, n.z);
#endif
	}
#if DEBUG
	else
	  fprintf (stderr, "edge %g %g %g %g %g %g %d\n",
		   x, y, z, closest.x - x, closest.y - y, closest.z - z, orient);
#endif
      }
    }
    else { // face
#if DEBUG
      fprintf (stderr, "face %g %g %g %g %g %g\n",
	       x, y, z, closest.x - x, closest.y - y, closest.z - z);
#endif
      orient = PointTriangleOrientation (&c, q[0].v, q[0].v+1, q[0].v+2);
    }
#endif // dimension == 3
    d[] = sqrt (q[0].d2)*orient; 
  }
  else { // !a->len
    free (a);
    surface[] = 0.;
    if (level > 0)
      d[] = bilinear (point, d);
    else
      d[] = 0.;
  }
}

#undef ND

static void refine_distance (Point point, scalar d)
{
  scalar surface = d.surface;
  if (surface[] == 0.)
    foreach_child() {
      surface[] = 0.;
      d[] = bilinear (point, d);
    }
  else {
    coord ** ap = (coord **) double_to_pointer (surface[]);
    int s = 0;
    foreach_child() {
      update_distance (point, ap, d);
      s += sign(d[]);
    }

To increase robustness to inconsistent input, we check whether all children are included within the minimum distance sphere. If this is the case then the children and parent must have the same orientation. We enforce this, using the “average” orientation.

    if (fabs(d[]) > sqrt(dimension)/4.*Delta) {
      if (abs(s) != 1 << dimension) {
	s = sign(s);
	foreach_child()
	  d[] = s*fabs(d[]);
      }
      if (sign(d[]) != sign(s))
	d[] = - d[];
    }
  }
}

static void restriction_distance (Point point, scalar d) {}

static void coarsen_distance (Point point, scalar d) {
  scalar surface = d.surface;
  foreach_child()
    free (double_to_pointer (surface[]));
}

static void delete_distance (scalar d) {
  scalar surface = d.surface;
  foreach_level (0)
    free (*((void **)double_to_pointer (surface[])));
  for (int l = 0; l <= depth(); l++)
    foreach_level (l)
      free (double_to_pointer (surface[]));
  delete ({surface});
}

trace
void distance (scalar d, coord * p)
{
  scalar surface = d.surface;
  if (surface.i)
    delete_distance (d);
  surface = new scalar;
  surface.restriction = no_restriction;
#if TREE
  surface.prolongation = no_restriction;
  surface.refine = no_restriction; // handled by refine_distance()
  d.prolongation = refine_bilinear;
  d.refine = refine_distance;  
  d.coarsen = coarsen_distance;
  d.dirty = true;
#endif
  d.surface = surface;
  d.delete = delete_distance;
  d.restriction = restriction_distance;

  Array * a = array_new();
  while (p->x != nodata) {
#if dimension == 3
    // filter degenerate triangles
    coord ab = vecdiff(*(p+1),*p), ac = vecdiff(*(p+2),*p);
    coord n = vecdotproduct(ab,ac);
    if (vecdot(n,n) > 0.)
#endif
      array_append (a, &p, sizeof (coord *));
    p += dimension;
  }
  p = NULL;
  array_append (a, &p, sizeof (coord *));
  p = (coord *) array_shrink (a);

  foreach_level (0, noauto)
    update_distance (point, (coord **) p, d);
  free (p);
  
  boundary_level ({d}, 0);
  for (int l = 0; l < depth(); l++) {
    foreach_coarse_level (l)
      refine_distance (point, d);
    boundary_level ({d}, l + 1);
  }
}

Usage

Examples

Tests

Computation of a levelset field from a contour