sandbox/Antoonvh/sketch.h

    Sketching functions for Bwatch

    Sketching utilities consist of a few functions and definitions.

    Sketch vof facets

    sketch_vof() mimics bview’s draw_vof() a little

    typedef struct _sketch_vof {
      scalar f;
      bool precomp; //precompute normals
      material mat;
    } _sketch_vof;
    
    _sketch_vof svlist[N_SKETCH];   // Arguments to sketch vof
    
    trace
    double get_pixel_sketch_vof (_sketch_vof inp, ray r, double DG, unsigned char * px) {
      double distance = HUGE;
      coord n, a;
      if (inp.precomp) {
        foreach_precomp_ray_facet_intersection(r, inp.f, DG) {
          distance = _distg;
          n = _n;
          a = _a[0];
        }
      }
      else {
        foreach_ray_facet_intersection(r, inp.f) {
          distance = _distg;
          n = _n;
          a = _a[0];
        }
      }
      if (distance < HUGE && distance < DG)
        if (get_color_material (r, a, n, inp.mat, px))
          return distance;
      return HUGE;
    }
    
    // User function:
    bool sketch_vof (struct _sketch_vof inp) {
      sketch_list[functot] = &get_pixel_sketch_vof;
      if (!inp.mat.dull && !inp.mat.ind && !any_col (inp.mat.col) &&
          !inp.mat.s.i && !inp.mat.v.x.i && !(inp.mat.R >= 1.))
        inp.mat.col[0] = inp.mat.col[1] = inp.mat.col[2] = 250;
      if (!inp.mat.map)
        inp.mat.map = jet;
      if (!inp.mat.min && !inp.mat.max) {
        inp.mat.min = -1;
        inp.mat.max =  1;
      }
      if (!inp.mat.SPexp && !inp.mat.SPR) {
        inp.mat.SPexp = 40;
        inp.mat.SPR = 0.1;
      }
      svlist[functot] = inp;
      functot++;
      shading = list_add (shading, inp.f);
      return true;
    }

    A Slice

    quadriangles may mimic bview’s squares() a bit.

    typedef struct _quadriangles {
      scalar s;
      double alpha;
      coord n;
      material mat;
    } _quadriangles;
    
    _quadriangles qlist[N_SKETCH]; // Arguments to quadriangles
    
    trace
    double get_pixel_quadriangles (_quadriangles inp, ray r, double DG, unsigned char * px) {
      coord a;
      double d = ray_plane_intersect (r, inp.n, inp.alpha, &a);
      if (d <= 0 || d >= DG)
        return HUGE;
      Point point = locate (a.x, a.y, a.z);
      if (point.level < 0) // Only slices in the domain
        return HUGE;
      normalize (&inp.n);
      if (get_color_material (r, a, inp.n, inp.mat, px))
        return d;
      return HUGE;
    }
    
    // User interface
    void quadriangles (struct _quadriangles inp) {
      sketch_list[functot] = &get_pixel_quadriangles;
      if (!any_comp (inp.n))
        inp.n.z = 1;
      if (!inp.mat.dull && !inp.mat.ind && !any_col (inp.mat.col) &&
          !inp.mat.s.i && !inp.mat.v.y.i && !(inp.mat.R >= 1)) {
        if (inp.s.i)
          inp.mat.s = inp.s;
        else 
          inp.mat.col[0] = 250, inp.mat.col[1] = 150, inp.mat.col[2] = 50;
      }
      if (!inp.mat.map)
        inp.mat.map = jet;
      if (!inp.mat.min && !inp.mat.max) {
        inp.mat.min = -L0;
        inp.mat.max =  L0;
      }
      qlist[functot] = inp;
      functot++;
    }

    The grid

    lattice() may mimic bview’s cells() a bit.

    typedef struct _lattice {
      double alpha;
      coord n;
      double width;
      material mat;
    } _lattice;
    
    _lattice llist[N_SKETCH];    // Arguments to quadriangles
    
    double get_pixel_lattice (_lattice inp, ray r, double DG, unsigned char * px) {
      coord a;
      double d = ray_plane_intersect (r, inp.n, inp.alpha, &a);
      if (d < 0 || d > 1e4*L0 || d > DG)
        return HUGE;
      Point point = locate (a.x, a.y, a.z);
      if (point.level < 0) {
        return HUGE;
      }
      coord cc = {x, y, z};
      //grid-alligned slices only
      foreach_dimension() {
        if (inp.n.z) {
          int sx = sign (a.x - cc.x);
          int sy = sign (a.y - cc.y);
          double dx = cc.x + sx*Delta/2 - a.x;
          double dy = cc.y + sy*Delta/2 - a.y;
          double dist = min (fabs(dx), fabs(dy));
          if (dist < inp.width) {
    	normalize (&inp.n);
    	if (get_color_material (r, a, inp.n, inp.mat, px)) {
    	  return d - 1e-4*L0;
    	}
          }
        }
      }
      return HUGE;
    }
    // User interface
    void lattice (struct _lattice inp) {
      sketch_list[functot] = &get_pixel_lattice;
      if (!inp.mat.dull && !inp.mat.ind && !any_col (inp.mat.col) &&
          !inp.mat.s.i && !inp.mat.v.x.i && !(inp.mat.R >= 1)) 
        for (int i = 0; i < 3; i++)
          inp.mat.col[i] = 50;
      if (!any_comp(inp.n))
        inp.n.z = 1.;
      if (!inp.width)
        inp.width = cam.fov/300.;
      llist[functot] = inp;
      functot++;
    }

    A disk

    A round planar disk

    typedef struct _disk {
      double R;
      coord P;
      coord n;
      material mat;
    } _disk;
    
    _disk dlist[N_SKETCH]; // Arguments to disk
    
    double get_pixel_disk (_disk inp, ray r,  double DG, unsigned char * px) {
      coord a;
      double alpha = 0;
      foreach_dimension()
        alpha += inp.n.x*inp.P.x;
      double d = ray_plane_intersect (r, inp.n, alpha, &a);
      if (d >= DG || d < 0.0001) //No further computations
        return HUGE;
      double D = 0;
      foreach_dimension()
        D += sq(a.x - inp.P.x);
      if (D > sq(inp.R))
        return HUGE;
      normalize (&inp.n);
      normalize (&r.dir);
      if (get_color_material (r, a, inp.n, inp.mat, px))
        return d;
      return HUGE;
    }
    
    bool disk (struct _disk inp) {
      sketch_list[functot] = &get_pixel_disk;
      if (!inp.R)
        inp.R = L0/2.;
      if (!any_comp (inp.n))
        inp.n.z = 1;
      if (!inp.mat.dull || !inp.mat.ind || !any_col (inp.mat.col) ||
          !inp.mat.s.i || !inp.mat.v.x.i || !(inp.mat.R >= 1.)) 
        inp.mat.R = 1;
      dlist[functot] = inp;
      functot++;
      return true;
    }

    Sketch a sphere

    typedef struct _sphere {
      double R;      //Radius
      coord C;       //centre
      material mat;
    } _sphere;
    
    _sphere slist[N_SKETCH];    // Arguments to sphere
    
    trace
    double get_pixel_sphere (_sphere inp, ray r, double DG, unsigned char * px) {
      coord a, n;
      double d = ray_sphere_intersect (r, inp.C, inp.R, &a, &n); 
      if (d > 0 && d < DG) { // possibly in sight
        if (get_color_material (r, a, n, inp.mat, px))
          return d;
      }
      return HUGE;
    }
    
    // User interface
    void sphere (struct _sphere inp) {
      sketch_list[functot] = &get_pixel_sphere;
      if (!inp.R)
        inp.R = 0.5;
      if (!inp.mat.dull && !inp.mat.ind && !any_col (inp.mat.col) &&
          !inp.mat.s.i && !inp.mat.v.y.i && !(inp.mat.R >= 1.)) 
       inp.mat.col[0]= 200, inp.mat.col[1] =  100, inp.mat.col[2] = 150;
     if (inp.mat.dull && !any_col(inp.mat.col)) {
       inp.mat.col[0]  = 50, inp.mat.col[1]  = 100, inp.mat.col[2]  = 150;
       inp.mat.col2[0] = 50, inp.mat.col2[1] = 150, inp.mat.col2[2] = 100;
       inp.mat.n1 = (coord){0.1, 0.9, 0};
     }
     if (!inp.mat.SPexp && !inp.mat.SPR) {
        inp.mat.SPexp = 40;
        inp.mat.SPR = 0.1;
      }
     if (!inp.mat.map)
        inp.mat.map = jet;
      if (!inp.mat.min && !inp.mat.max) {
        inp.mat.min = -L0;
        inp.mat.max = L0;
      }
      slist[functot] = inp;
      functot++;
    }

    Equiplane

    An (Curved) Plane where a scalar field s has the value val. It can also be used to sketch a smooth vof-interface reconstruction.

    typedef struct _equiplane {
      scalar s;
      double val;
      bool insideout;
      bool vof;
      material mat;
    } _equiplane;
    
    _equiplane eqlist[N_SKETCH];
    
    trace
    double get_pixel_equiplane (struct _equiplane inp, ray r, double DG, unsigned char * px) {
      scalar s = inp.s;
      // Draw the distance equiplane for vof
      if (inp.vof) {
        s = s.vofd;
        if (inp.insideout) //Invert 
          inp.insideout = false;
        else
          inp.insideout = true;
      }
      scalar posi = s.possible;
      vector v = s.normals;
      double d = HUGE;
      coord n, a;
      foreach_possible_ray_equi_intersection(r, s, posi, DG) {
        if (_dist < d) {
          double vals[2];
          foreach_dimension() {
    	_a[0].x += r.dir.x*L0*1e-7;
    	_a[1].x -= r.dir.x*L0*1e-7;
          }
          for (int i = 0; i < 2; i++) 
    	vals[i] = interpolate(s, _a[i].x, _a[i].y, _a[i].z);
          if ((vals[0] - inp.val) * (vals[1] - inp.val) < 0.) {
    	double w = fabs(vals[0] - inp.val)/(fabs(vals[0] - inp.val) + fabs(vals[1] - inp.val));
    	if (w < 1e-6 || w > 1 - 1e-6)
    	  return HUGE;
    	d = _dist;
    	foreach_dimension() {
    	  a.x = w*_a[1].x + (1. - w)*_a[0].x;
    	}
    	n.x = interpolate (v.x, a.x, a.y, a.z); 
    	n.y = interpolate (v.y, a.x, a.y, a.z); 
    	n.z = interpolate (v.z, a.x, a.y, a.z);
          }
        }
      }
      if (any_comp(n) && d < DG) {
        if (inp.insideout) {
          foreach_dimension()
    	n.x *= -1;
        }
        normalize (&n);
        if (get_color_material (r, a, n, inp.mat, px))
          return d;
      }
      return HUGE;
    }
    
    void equiplane (struct _equiplane inp) {
      sketch_list[functot] = &get_pixel_equiplane;
      if (!inp.mat.dull && !inp.mat.ind && !any_col (inp.mat.col) &&
           !inp.mat.s.i && !inp.mat.v.x.i && !(inp.mat.R >= 1.)) 
         inp.mat.col[0]= 100, inp.mat.col[1] =  200, inp.mat.col[2] = 150;
      if (!inp.mat.SPexp && !inp.mat.SPR) {
        inp.mat.SPexp = 40;
        inp.mat.SPR = 0.1;
      }
      if (!inp.mat.map)
        inp.mat.map = jet;
      if (!inp.mat.min && !inp.mat.max) {
        inp.mat.min = -L0;
        inp.mat.max = L0;
      }
      eqlist[functot] = inp;
      functot++;
    }
    
    #define PIXEL_INDEX (((ix + inp.res*jy)*3))
    
    typedef struct _image {
      char * fname;
      int res;
      unsigned char * imgrgb; //[res*res*3];
      coord n;
      double alpha;
      coord up;
      material mat;
    } _image;
    
    _image ilist[N_SKETCH]; 
    
    double get_pixel_image (struct _image inp, ray r, double DG, unsigned char * px) {
      coord a;
      double d = ray_plane_intersect (r, inp.n, inp.alpha, &a);
      if (d <= 0 || d >= DG)
        return HUGE;
      Point point = locate (a.x, a.y, a.z);
      if (point.level < 0) // Only slices in the domain
        return HUGE;
      normalize (&inp.n);
      {;}
      //Only grid aligned projections
      foreach_dimension() {
        if (inp.n.z) {
          int ix = max(min(inp.res*(a.x - X0)/L0, inp.res - 1), 0);
          int jy = max(min(inp.res*(Y0 + L0 - a.y)/L0, inp.res - 1), 0);
          for (int i = 0; i < 3; i++) 
    	inp.mat.col[i] = max(inp.imgrgb[PIXEL_INDEX + i], 1);
        }
      }
      if (get_color_material (r, a, inp.n, inp.mat, px))
        return d;
      return HUGE;
    }
    
    bool image (struct _image inp) {
      if (!any_comp(inp.n))
        inp.n.z = 1;
      if (!inp.alpha)
        inp.alpha = 0;
      if (!inp.res)
        inp.res = (1 << grid->maxdepth);
      if (inp.fname && inp.imgrgb == NULL) {
        inp.imgrgb = malloc (3*sizeof(unsigned char)*sq(inp.res));
        char cmd[999];
        sprintf (cmd, "convert -quiet %s -resize %dx%d^ \
                     -gravity center -extent %dx%d 2not8exist5.ppm",
    	     inp.fname, inp.res, inp.res, inp.res, inp.res);
        assert(!system (cmd));
        FILE * fpi = fopen ("2not8exist5.ppm", "rb");
        char line_one[2];
        int height, width, max_color;
        fscanf(fpi, "%s %d %d %d\n", line_one, &height, &width, &max_color);
        // printf ("%s %d %d %d\n", line_one, height, width, max_color);
        
        fread(inp.imgrgb, 3*sizeof(unsigned char), sq(inp.res), fpi);
        fclose (fpi);
      } else
        assert (0);
       ilist[functot] = inp;
       sketch_list[functot] = &get_pixel_image;
       functot++;
      return true;
    }

    A plain sheet for sketching…

    plain() mimics bview’s clear() a bit.

    void plain (void) {
      for (int func = 0; func < functot; func++) 
        if (sketch_list[func] == &get_pixel_image) { 
          free (ilist[func].imgrgb);
          ilist[func].imgrgb = NULL;
        }
      functot = 0;
      free (shading);
      shading = NULL;
     
      
    }

    The ray-caster function

    The implementation of a helper function that gets the rgb values of a ray by cycling though all sketch function calls…

    double get_color_ray (ray r, unsigned char * px) {
      double Distg = HUGE;
      for (int func = 0; func < functot; func++) {
        double dist = HUGE;
        unsigned char pt[3];
        if (sketch_list[func] == &get_pixel_quadriangles) 
          dist = sketch_list[func](qlist[func], r, Distg, pt);
        else if (sketch_list[func] == &get_pixel_sketch_vof) 
          dist = sketch_list[func](svlist[func], r, Distg, pt);
        else if (sketch_list[func] == &get_pixel_disk) 
          dist = sketch_list[func](dlist[func], r, Distg, pt);
        else if (sketch_list[func] == &get_pixel_lattice) 
          dist = sketch_list[func](llist[func], r, Distg, pt);
        else if (sketch_list[func] == &get_pixel_sphere) 
          dist = sketch_list[func](slist[func], r, Distg, pt);
        else if (sketch_list[func] == &get_pixel_equiplane) 
          dist = sketch_list[func](eqlist[func], r, Distg, pt);
        else if (sketch_list[func] == &get_pixel_image) 
          dist = sketch_list[func](ilist[func], r, Distg, pt);
        if (dist < Distg) {
          Distg = dist;
          for (int i = 0; i < 3; i++)
    	px[i] = pt[i];
        }
      }
      return Distg;
    }

    Preparation

    Before rendering, it maybe usefull to compute some helper fields to store data for rendering or help guide optimized grid iterators.

    static inline double point_value (scalar s, Point point) {
      return s[];
    }
    
    void preparation_steps (void){
      for (int func = 0; func < functot; func++) {
        //Gradients and possibles for isosurfaces,
        // maybe compute the nearby distances to vof fields.
        if (sketch_list[func] == &get_pixel_equiplane) {
          scalar s = eqlist[func].s;
          if (eqlist[func].vof) {
    	char dname[99];
          	sprintf (dname, "dfor%s%d", s.name, func); 
          	scalar d = new_scalar (dname);
    	s.vofd = d;
          	find_nearby_distances (s);
    	s = d; // we will draw the isosurface of d
          }
          char vname[99], pname[99];
          sprintf (vname, "vfor%s%d", s.name, func);
          sprintf (pname, "pfor%s%d", s.name, func);
          vector v = new_vector (vname);
          scalar p = new_scalar (pname);
          s.normals = v;
          s.possible = p;
          find_possible (s, eqlist[func].val);
          gradients ({s}, {v});
        }
        // Precompute normals and plane alphas for vof facets
        else if (sketch_list[func] == &get_pixel_sketch_vof) {
          if (svlist[func].precomp) {
    	scalar f = svlist[func].f;
    	char vname[99], aname[99];
    	sprintf (vname, "vfor%s%d", f.name, func);
    	sprintf (aname, "afor%s%d", f.name, func);
    	vector v = new_vector (vname);
    	scalar alpha = new_scalar (aname);
    	f.normals = v;
    	f.possible = alpha; //alpha
    	foreach() {
    	  alpha[] = nodata;
    	  foreach_dimension()
    	    v.x[] = nodata;
    	  if (interfacial2 (f, point)) {
    	    coord n = mycs (point, f); //pointing outwards.
    	    double l = vec_length(n);
    	    alpha[] = plane_alpha (point_value(f, point), n)/l;
    	    normalize (&n);
    	    foreach_dimension()
    	      v.x[] = n.x;
    	  }
    	}
    	boundary ({alpha, v});
          }
        } 
      }
    }
    
    void aftercare (void) {
      for (int func = 0; func < functot; func++) {
        if (sketch_list[func] == &get_pixel_equiplane) {
          scalar s = eqlist[func].s;
          if (eqlist[func].vof)
    	s = s.vofd;
          vector v = s.normals;
          scalar posi = s.possible;
          delete ({posi, v});
          if (eqlist[func].vof) 
    	delete ({s});
        }
        else if (sketch_list[func] == &get_pixel_sketch_vof) 
          if (svlist[func].precomp) {
    	scalar s = svlist[func].f;
    	vector v = s.normals;
    	scalar alpha = s.possible;
    	delete ((scalar*){alpha, v});
          }
      }
    }

    Store

    Write a .ppm to the FILE * fp.

    // Default background
    #define BGR (0)
    #define BGG (25.*(ii + jj)/cam.nx + 120)
    #define BGB (50.*jj/cam.nx + 120)
    
    trace
    bool store (FILE * fp) {
      if (lights[0].type == 0) // No lights
        default_lights();
      fprintf (fp, "P6\n%d %d\n%d\n", cam.nx, cam.ny, 255); // .ppm header
      preparation_steps();
    #if _OPENMP
    #pragma omp parallel shared(fp)
      {
        assert (!(cam.nx*cam.ny % npe()));
        unsigned char px[3];
        foreach_ray() {
          int jjj = cam.ny - jj;
          if ((jjj*cam.nx + ii) % npe() == tid()) {
    	px[0]= BGR, px[1] = BGG, px[2] = BGB;
    	get_color_ray (_r, px);
          }
          if ((jjj*cam.nx + ii) % npe() == (npe() - 1)) 
    	for (int i = 0; i < npe(); i++) {
    	  if (tid() == i)
    	    fwrite (px, 1, 3, fp);
            #pragma omp barrier
    	}
        }
      }  
    #else
      foreach_ray() {
        unsigned char px[3] = {BGR, BGG, BGB};
        get_color_ray (_r, px);
        fwrite (px, 1, 3, fp);
      }
    #endif
      fflush (fp);
      aftercare();
      return true;
    }