sandbox/Antoonvh/bwatch-iterators.h

    Iterators bwatch

    Concept

    Iterate over

    • All camera rays
    • Each cell that the ray passes
    • Each cell where a ray intersects a facet
    // A helper function
    bool interfacial2 (scalar f, Point point) {
      if (f[] > 1e-6 && f[] < 1 - 1e-6)
        return true;
      return false;
    }
    
    @def foreach_ray() 
    coord cr, proj, vert, hori;
    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;
    }
    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;
        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);  
        }
        ray _r;
        _r.depth = 0;
        _r.O = cam.O;
        foreach_dimension()
          _r.dir.x = apoint.x - cam.O.x;
        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();
    @

    Workhorses

    sketch should use optimized iterators,

    #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;
      long 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();
    @

    Iterator for volumetric rendering

    long int find_possible_vol (scalar s, double mval) {
      scalar posi = s.possible;
      long 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();
    @