// A helper function
bool interfacial2 (scalar f, Point point) {
  if (f[] > 1e-6 && f[] < 1 - 1e-6)
    return true;
  return false;
}

void set_cr (coord * cr, coord * proj, coord * vert) {
  foreach_dimension() 
    cr->x = cam.poi.x - cam.O.x;
  normalize (cr);
  foreach_dimension() {
    proj->x = cam.up.x*cr->x;
    vert->x = cam.up.x - proj->x;
  }
}

void set_apoint (coord * apoint, coord hori, coord vert, int ii, int jj) {
  foreach_dimension() {
    apoint->x = (cam.poi.x
		 + cam.fov*((ii + 0.5)/cam.nx - 0.5)*hori.x
		 + ((float)cam.ny/cam.nx)*cam.fov*((jj + 0.5)/cam.ny - 0.5)*vert.x);  
  }
}

@def foreach_ray() 
coord cr, proj, vert, hori;
set_cr (&cr, &proj, &vert);
hori = cross (cr, vert);
normalize (&vert);
normalize (&hori);
for (int jj = cam.ny; jj > 0; jj--) {    
  for (int ii = 0; ii < cam.nx; ii++) {
    coord apoint;
    set_apoint (&apoint, hori, vert, ii, jj);
    ray _r;
    _r.depth = 0;
    _r.O = cam.O;
    _r.dir.x = apoint.x - cam.O.x;
    _r.dir.y = apoint.y - cam.O.y;
    _r.dir.z = apoint.z - cam.O.z;
    
    normalize (&_r.dir); 
    @
      @def end_foreach_ray()
      }
 }
@

@def foreach_ray_cell_intersection(r) 
  foreach_cell() { //MG acceleration
  coord _a[2] = {0};
    double _dist = ray_cell_intersect (r, point, _a);
    if (_dist < HUGE) {
      if (is_leaf(cell)) {
	@
	  @def end_foreach_ray_cell_intersection()
	  }
    } else // Ray missed, skip children
  continue;
  } end_foreach_cell();
@

@def foreach_ray_facet_intersection(r, f)
double _distg = HUGE;
foreach_cell() { //MG acceleration
  coord _a[2] = {{0,0,0}, {0,0,0}}; //Pedantic...
  double _dist = ray_cell_intersect (r, point, _a);
  if (_dist < 0 || _dist < _distg) {
    if (interfacial2 (f, point)) {
      if (is_leaf(cell)) {
	coord _n;
	if (segment_facet_intersect (_a, f, point, &_n)) {
	  normalize (&_n);
	  _distg = _dist;
	  @
	    @def end_foreach_ray_facet_intersection()
	    }
      }
    } else  // No fragment, skip children
      continue;
  } else   // Ray missed or blocked, skip children
    continue; 
}end_foreach_cell();
@

#include "PointTriangle.h"

void find_nearby_distances (scalar s) {
  assert (s.vofd.i);
  boundary ({s});
  scalar d = s.vofd;
  foreach() {
    coord cc = {x, y, z};
    double dist = 3*Delta*sign (s[] - 0.5);
    foreach_neighbor (1) {
      if (s[] > 1e-6 && s[] < 1. - 1e-6) {
	coord n = mycs (point, s);
	double alpha = plane_alpha (s[], n);
	coord v[12];
	int m = facets (n, alpha, v, 1.);
	// Using Alexis Berny's triangulation
	for (int j = 0; j < m - 2; j++) {
	  coord t1 = {x + v[0].x*Delta  , y + v[0].y*Delta  , z + v[0].z*Delta}; 
	  coord t2 = {x + v[j+1].x*Delta, y + v[j+1].y*Delta, z + v[j+1].z*Delta};
	  coord t3 = {x + v[j+2].x*Delta, y + v[j+2].y*Delta, z + v[j+2].z*Delta};
	  double s, t, dtr = PointTriangleDistance (&cc, &t1, &t2, &t3, &s, &t);
	  if (dtr > 0)
	    if (fabs(sqrt(dtr)) < fabs(dist))
	      dist = sqrt(dtr)*PointTriangleOrientation (&cc, &t1, &t2, &t3);
	}
      }
    }
    d[] = dist;
  }
  boundary ({d});
}

long int find_possible (scalar s, double val) {
  scalar posi = s.possible;
  int n = 0;
  foreach(reduction (+:n)) {//leaf
    posi[] = 0;
    bool pos = false, neg = false; 
    foreach_neighbor(1) {
      if (s[] < val)
	neg = true;
      if (s[] > val)
	pos = true;
    }
    if (pos && neg) {
      n++;
      posi[] = 1;
    }
  }
  for (int l = depth() - 1; l >=0 ; l--) {
    foreach_coarse_level(l) {
      double p = 0;
      foreach_child()
	if (posi[])
	  p = 1;
      posi[] = p;
    }
  }
  return n;
}

bool maybe (scalar posi, Point point) {
  if (posi[])
    return true;
  return false;
}

@def foreach_possible_ray_equi_intersection(r, s, posi, DG)
foreach_cell() { //MG acceleration
  if (maybe(posi, point)) {
    coord _a[2] = {0};
    double _dist = ray_cell_intersect (r, point, _a);
    if (_dist < DG) {
      if (is_leaf(cell)) {
	@
	  @def end_foreach_possible_ray_equi_intersection()
	  }
    } else // Ray missed or blocked, skip children
      continue;
  } else // Does not contain the isosurface
    continue;
} end_foreach_cell();
@

static inline bool precomp_segment_facet_intersect (coord a[2], double alpha, Point point, coord n) {
  coord cc = {x, y, z};
  double ALP = 0, ALP2 = 0;
  foreach_dimension() {
    ALP  += n.x*(a[0].x - cc.x)/Delta;
    ALP2 += n.x*(a[1].x - cc.x)/Delta;
  }
  if ((ALP2 - alpha)/(ALP - alpha) > 0.05) //~/wiki/sandbox/Antoonvh 5% gap filling
    return false;
  double w = fabs((ALP2 - alpha)) / (fabs(ALP - alpha) + fabs(ALP2 - alpha));
  foreach_dimension() {
    a[1].x = w*a[0].x + (1 - w)*a[1].x;
  }
  {;}
  return true;
}

static inline coord get_normal (scalar f, Point point) {
  coord n1 = {0};
  vector v = f.normals;
  foreach_dimension()
    n1.x = v.x[];
  return n1;
}

static inline double get_alpha (scalar f, Point point) {
  scalar s = f.possible;
  return s[];
}
  
@def foreach_precomp_ray_facet_intersection(r, f, DG)
double _distg = HUGE;
foreach_cell() { //MG acceleration
  if (interfacial2 (f, point)) {
    coord _a[2] = {{0,0,0}, {0,0,0}}; //Pedantic...
    double _dist = ray_cell_intersect (r, point, _a);
    if (_dist >= 0 && _dist < _distg && _dist < DG) {
      if (is_leaf(cell)) {
	coord _n = get_normal (f, point);
	double alpha = get_alpha (f, point);
	if (precomp_segment_facet_intersect (_a, alpha, point, _n)) {
	  normalize (&_n);
	  _distg = _dist;
	  @
	    @def end_foreach_precomp_ray_facet_intersection()
	    }
      }
    } else  // Ray missed or blocked, skip children
      continue;
  } else   // No fragment, skip children
    continue; 
}end_foreach_cell();
@

long int find_possible_vol (scalar s, double mval) {
  scalar posi = s.possible;
  int n = 0;
  foreach(reduction(+:n)) {//leaf
    posi[] = 0;
    bool nv = false;
    foreach_neighbor(1)
      if (s[] > mval)
	nv = true;
    if (nv)  {
      posi[] = 1;
      n++;
    } 
  }
  // restriction
  for (int l = depth() - 1; l >= 0; l--) {
    foreach_coarse_level(l) {
      double p = 0;
      foreach_child()
	if (posi[])
	  p = 1;
      posi[] = p;
    }
  }
  return n;
}

@def foreach_ray_cell_intersection_volume(r, dm, posi) 
foreach_cell() { //MG acceleration
  if (maybe(posi, point)) {
    coord _a[2] = {0};
    double _dist = ray_cell_intersect (r, point, _a);
    if (_dist < dm) {
      if (is_leaf(cell)) { //do something
	@
	  @def end_foreach_ray_cell_intersection_volume()
	  }
    } else // Ray missed or blocked, skip children
      continue;
	} else // Nothing to see, skip children 
    continue;
} end_foreach_cell();
@