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;
      bool HQ; // High quality, some background object required
    } _lattice;
    
    _lattice llist[N_SKETCH]; // Arguments to lattice
    
    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.HQ ? 3 : 1)*inp.width && d > Delta*.1) {
    	normalize (&inp.n);
    	if (inp.HQ) {
    	  double b = 0, c = 0;
    	  foreach_dimension() {
    	    b += sq(cam.O.x - a.x);
    	    c += sq(cam.O.x - cam.poi.x);
    	  }
    	  b = sqrt(b);
    	  c = sqrt(c);
    	  double ps = cam.fov*b/(cam.nx*c);
    	  coord np = {0,0,0}, rp = {0, 0, 0};
    	  np.z = 1.;
    	  rp.z = r.dir.z;
    	  if(dist == fabs(dx)) 
    	    rp.y = r.dir.y;
    	  else
    	    rp.x = r.dir.x;
    	  normalize (&rp);
    	  double cc = fabs(dot(rp, np));
    	  double ss = cc > 0. ? ps/cc : ps;
    	  //printf ("%g %g %g\n", ss, ps, cc);
    	  double temp = dist + ss/4. < inp.width ? 1 :
    	    dist + ss/4. < 2*inp.width ?
    	    1 - (dist + ss/4. - inp.width)/(inp.width) : 0;
    	  temp += fabs(dist - ss/4.) < inp.width ? 1 :
    	    fabs(dist - ss/4.) < 2*inp.width ?
    	    1 - (fabs(dist - ss/4.) - inp.width)/(inp.width) : 0;
    	  temp = (2. - temp)/2.;
    	  inp.mat.T = max(min(temp, 1), 0);
    	  foreach_dimension()
    	    a.x -= 0.04*Delta*r.dir.x;
    	}
    	if (get_color_material (r, a, inp.n, inp.mat, px)) 
    	  return d - Delta/20.;
          }
        }
      }
      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] = 1;
      if (!any_comp(inp.n))
        inp.n.z = 1.;
      if (!inp.width)
        inp.width = cam.fov/(40*N);
      if (inp.HQ)
        inp.width /= 1.5;
      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];
          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) - (vals[1] - inp.val));
    	if (w < 0 || w > 1)
    	  continue;
    	double mrdir = 0;
    	normalize (&r.dir);
    	foreach_dimension() {
    	  a.x = w*_a[1].x + (1. - w)*_a[0].x;
    	  if (fabs(r.dir.x) > mrdir) {
    	    mrdir = fabs(r.dir.x);
    	    d = (a.x - r.O.x)/r.dir.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 && d > 0) {
        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++;
    }

    Image in the scene

    Place an image in the scene

    #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;
      unsigned char greenscreen[3]; // Transparent color code (non black green screen)
      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
      coord Ori = {X0, Y0, Z0};
      foreach_dimension() {
        if (inp.n.z) {
          int ix = max(min(inp.res*(a.x - Ori.x)/L0, inp.res - 1), 0);
          int jy = max(min(inp.res*(Ori.y + 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 (any_col (inp.greenscreen)) {
        int count = 0;
        for (int i = 0; i < 3; i++) 
          if (inp.greenscreen[i] == inp.mat.col[i])
    	count++;
        if (count == 3) //its the green-screen val
          return HUGE;
      }
      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 -interpolate Nearest -filter point -resize %dx%d^  \
                     -gravity center -extent %dx%d -depth 8 2not8exist5.ppm",
    	     inp.fname, inp.res, inp.res, inp.res, inp.res);
        assert(!system (cmd));
        FILE * fpi = fopen ("2not8exist5.ppm", "r");
        char line_one[2];
        int height, width, max_color;
        fscanf(fpi, "%s\n%d %d\n%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;
    }

    Display text

    It is possible to add text using convert and image

    typedef struct _sketch_text {
      char * str;
      int fs;              // Font size (compare against .res)
      unsigned char tc[3]; // text color
      char * pos;          // Imagemagick Position (center, north, south west etc..)
      int res;             // Text image resolution
      char * ops;          // Options to pass to convert
      coord n;
      double alpha;
      coord up;
      material mat;
    } _sketch_text;
      
    bool sketch_text (struct _sketch_text inp) {
      char fname[99] = "3not9exist6.ppm";
      char cmd[999];
      if (!inp.fs) 
        inp.fs = 100;
      if (!inp.res)
        inp.res = 1000;
      unsigned char tc[3] = {0,0,0}; //black
      char d_ops[2] = " ";
      if (!inp.ops)
        inp.ops = d_ops;
      if (any_col(inp.tc)) {
        tc[0] = inp.tc[0];
        tc[1] = inp.tc[1];
        tc[2] = inp.tc[2];
      }
      unsigned char gs[3]; //Set Green screen color
      for (int i = 0; i < 3; i++) {
        int tci = tc[i];
        gs[i] = tci - 2 <= 0 ? tc[i] + 2 : tc[i] - 2;
      }
      char tmp[99] = "northwest";
      if (!inp.pos)
        inp.pos = tmp;
      sprintf (cmd, "convert -quiet %s -size %dx%d -background \"rgb(%d, %d, %d)\"\
    	   -fill \"rgb(%d,%d,%d)\" -pointsize %d -gravity %s caption:'%s' -depth 8 %s",
    	   inp.ops, inp.res, inp.res, gs[0], gs[1], gs[2], tc[0], tc[1], tc[2],
    	   inp.fs, inp.pos, inp.str, fname);
      //printf ("%s \n", cmd);
      system (cmd);
      struct _image inpn;
      inpn.fname = fname; 
      inpn.greenscreen[0] = gs[0];
      inpn.greenscreen[1] = gs[1];
      inpn.greenscreen[2] = gs[2];
      inpn.res = inp.res;
      inpn.n = inp.n;
      inpn.alpha = inp.alpha;
      inpn.up = inp.up;
      inpn.mat = inp.mat;
      inpn.imgrgb = NULL;
      return image (inpn);
    }

    A volumetric object

    It is possible to sketch volumetric objects as a smoke concentration field. (0 ~ mval < s[])

    #include "radix.h" //A sorting algorithm
    bool sketch_a_volume = false;
    typedef struct _volume {
      scalar s;
      double sc;            // 1/e scale for s times length
      double mval;          // Minimal value (positive and close to zero)
      unsigned char col[3]; // prescribed Color
      bool cols;            // Color by scalar value (color map)
      vector colorv;        // Color map 3D vector direction xyz -> rgb
      colormap map;         // color Map (default is cool_warm)
      double min;           // min cbar
      double max;           // max cbar
      double shading;       // shading effect for default lights
    } _volume;
    
    _volume vollist; //max 1
    
    double get_attenuation (coord a, struct _volume inp) {
      ray r = {.O = a, .dir = lights[1].dir};
      foreach_dimension()
        r.dir.x *= -1;
      normalize (&r.dir);
      scalar s = inp.s;
      double sl = 0;
      foreach_ray_cell_intersection_volume (r, HUGE, s.possible) {
        coord mean = {0, 0, 0};
        double len = 0;
        foreach_dimension() {
          mean.x = (_a[0].x + _a[1].x)/2.;
          len += sq(_a[0].x - _a[1].x);
        }
        double val = interpolate (inp.s, mean.x, mean.y, mean.z);
        if (val > inp.mval) {
          len = sqrt(len);
          if (len > 0)
    	sl += fabs(val)*len;
        }
      }
      return sl;
    }
    
    trace
    void mod_pixel_volume (struct _volume inp, ray r, double DG, unsigned char * px) {
      scalar s = inp.s;
      int np = 50, incf = 2; // Array size, increase factor;
      double * dvl = NULL;
      int elm = 0;
      if (inp.cols) {
        elm = 3 + (inp.shading > 0) + 3*(inp.colorv.x.i > 0); 
        dvl = malloc (elm*sizeof(double)*np); //Distance, values and length;
      }
      int celld = 0; // Number if cells
      double sl = 0; //integral of s times length;
      foreach_ray_cell_intersection_volume(r, DG, s.possible) {
        coord mean = {0, 0, 0};
        double len = 0;
        foreach_dimension() {
          mean.x = (_a[0].x + _a[1].x)/2.;
          len += sq(_a[0].x - _a[1].x);
        }
        double val = interpolate (inp.s, mean.x, mean.y, mean.z);
        if (val > inp.mval) {
          len = sqrt(len);
          if (len > 0)
    	sl += fabs(val)*len;
          if (inp.cols) { // Record data 
    	celld++;
    	if (celld >= np) { // increase array size
    	  np *= incf;
    	  dvl = realloc (dvl, elm*sizeof(double)*np);
    	}
    	double dirm = 0;
    	foreach_dimension() { //distance
    	  if (fabs(r.dir.x) > dirm)
    	    dvl[elm*(celld - 1)] = (mean.x - r.O.x)/r.dir.x;
    	}
    	assert (dvl[elm*(celld - 1)] > 0);
    	dvl[elm*(celld - 1) + 1] = val;
    	dvl[elm*(celld - 1) + 2] = len;
    	if (inp.shading)
    	  dvl[elm*(celld - 1) + 3] = get_attenuation (mean, inp);
    	if (inp.colorv.x.i) {
    	  double xp = mean.x, yp = mean.y, zp = mean.z;
    	  int dim = 0;
    	  foreach_dimension()
    	    dvl[elm*(celld - 1) + 4 + dim++] = interpolate (inp.colorv.x, xp, yp, zp);
    	}
          }
        }
      }
      if (inp.cols && sl > 0) { // Compute effective absorbed color
        int ind[celld], dista[celld];
        for (int it = 0; it < celld; it++) 
          dista[it] = (int)(1 << (depth() + 1))*(dvl[elm*it]/L0); // distance in integer units
        radixsort (celld, dista, ind);
        double col[3] = {0, 0, 0};
        double cmap[NCMAP][3];
        inp.map (cmap);
        double sl2 = 0;
        double TW = 0;
        for (int it = 0; it < celld; it++) {
          unsigned char coll[3];
          if (inp.colorv.x.i) {
    	int mx = 0;
    	for (int i = 0; i < 3; i++) {
    	  coll[i] = 255*min(1, fabs(dvl[elm*ind[it] + 4 + i])/inp.max); 
    	  if (coll[i] > mx)
    	    mx = coll[i];
    	}
    	for (int i = 0; i < 3; i++) 
    	  coll[i] = mx > 0 ? min (255, (255*coll[i])/mx) : 1; 
          } else 
    	colormap_pigmentation (coll, cmap, dvl[elm*ind[it] + 1], inp.min, inp.max);
          double sle = dvl[elm*ind[it] + 2]*fabs(dvl[elm*ind[it] + 1]);
          double Weight = exp(-sl2/inp.sc) - exp(-(sl2 + sle)/inp.sc);
          double fac = 1;
          if (inp.shading) { 
    	double I = lights[0].I;
    	I += lights[1].I*exp(-dvl[elm*ind[it] + 3]/(smoke.att));
    	fac =  min (max(cam.f(I) - cam.f(cam.min),0)/cam.f(cam.max/cam.min), 1);
          }
          TW += Weight;
          for (int i = 0; i < 3; i++) 
    	col[i] += (double)coll[i]*Weight*fac;
          sl2 += sle;
        }
        if (TW > 0) {
          for (int i = 0; i < 3; i++) 
    	inp.col[i] = (unsigned char)(col[i]/TW);
        }
      }
      if (sl > 0) {
        double W = exp(-sl/inp.sc);
        px[0] = (1 - W)*inp.col[0] + W*px[0];
        px[1] = (1 - W)*inp.col[1] + W*px[1];
        px[2] = (1 - W)*inp.col[2] + W*px[2];
      }
      free (dvl); dvl = NULL;
    }
    
    bool volume (struct _volume inp) {
      if (!inp.sc)
        inp.sc = 1.;
      if (!any_col(inp.col)) {
        inp.col[0] = 255;
        inp.col[1] = 255;
        inp.col[2] = 255;
      }
      if (!inp.max)
        inp.max = 1;
      if (!inp.mval)
        inp.mval = 0;
      if (!inp.min)
        inp.min = inp.mval;
      if (inp.shading) {
        scalar s = inp.s;
        smoke = s;
        s.att = inp.sc/inp.shading;
      }
      if (!inp.map)
        inp.map = cool_warm;
      vollist = inp;
      sketch_a_volume = true;
      return sketch_a_volume;
    }

    Triangles

    Draw a list of triangles, stored as a nodata terminated list of vertex triplets (e.g. as read by inport_stl).

    typedef struct _triangles {
      coord * c;
      coord bb[2];  //Bounding box (automatically computed)
      material mat;
    } _triangles;
    
    _triangles tlist[N_SKETCH];
    
    trace
    double get_pixel_triangles (struct _triangles inp, ray r, double DG, unsigned char * px) {
      coord * p = inp.c;
      double dist = DG;
      coord a, n;
      bool any = false;
      if (ray_box (r, inp.bb) < dist) {
        while (p->x != nodata) {
          double d = ray_triangle (r, p);
          if (d < dist) {
    	dist = d;
    	any = true;
    	foreach_dimension()
    	  a.x = r.O.x + r.dir.x*d;
    	coord a1, a2;
    	foreach_dimension() {
    	  a1.x = p[1].x - p[0].x;
    	  a2.x = p[2].x - p[0].x;
    	}
    	n = cross (a1, a2);
    	normalize (&n);
          }
          p += dimension; //3?
        }
        if (any) 
          if (get_color_material (r, a, n, inp.mat, px))
    	return dist;
      }
      return HUGE;
    }
    
    struct _triangles get_bb (struct _triangles in) {
      struct _triangles inp = in;
      foreach_dimension() {
        inp.bb[0].x = HUGE;
        inp.bb[1].x = -HUGE;
      }
      while (in.c->x < nodata) {
        foreach_dimension() {
          if (in.c->x < inp.bb[0].x)
    	inp.bb[0].x = in.c->x;
          if (in.c->x > inp.bb[1].x)
    	inp.bb[1].x = in.c->x;
        }
        in.c++;
      }  
      return inp; 
    }
    
    bool triangles (struct _triangles inp) {
      sketch_list[functot] = &get_pixel_triangles;
      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]= 200, inp.mat.col[1] =  150, inp.mat.col[2] = 100;
      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;
      }
      inp = get_bb (inp);
      tlist[functot] = inp;
      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;
      sketch_a_volume = false;
    }

    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);
        else if (sketch_list[func] == &get_pixel_triangles) 
          dist = sketch_list[func](tlist[func],  r, Distg, pt);
        if (dist < Distg && dist > 0) {
          Distg = dist;
          for (int i = 0; i < 3; i++)
    	px[i] = pt[i];
        }
      }
      if (sketch_a_volume) // Modify this pixel color 
        mod_pixel_volume (vollist, r, Distg, px);
      return Distg;
    }
    
    static inline double point_value (scalar s, Point point) {
      return s[];
    }

    Watch the sketch process

    You may watch the sketching process live in a (minimalistic) X11 window. You should have the libX11 library installed and compile using something like:

    $ qcc -O2 -DWATCH_ALONG bwatch.c -lm -lX11
    $ ./a.out

    Note that it may slow the rendering, especially with openmp.

    #if WATCH_ALONG
    #include <X11/Xlib.h> 
    GC gc_2;
    Window win = (Window)NULL;
    Display *dsp;
    int screen;
    Colormap cmap;
    
    void XSetup (void) {
      if (win == (Window)NULL) {
    #if _OPENMP
        XInitThreads();
    #pragma omp barrier
    #endif
        dsp    = XOpenDisplay(NULL);
        screen = DefaultScreen(dsp);
        win = XCreateWindow(dsp, DefaultRootWindow(dsp),
    			0, 0, cam.nx, cam.ny,
    			0, 0, 0, 0, 0, 0);
        XGCValues gcvalues_2;
        gcvalues_2.function = GXcopy;
        gcvalues_2.plane_mask = AllPlanes;
        gcvalues_2.foreground = 0x00FF00;
        gcvalues_2.background = 0xFFFFFF;
        gc_2 = XCreateGC(dsp, win,
    		     GCFunction|GCPlaneMask|GCForeground|GCBackground,
    		     &gcvalues_2);
        XStoreName (dsp, win, "Bwatch: Watch Along\0");
        XMapWindow(dsp, win);
        cmap= XDefaultColormap(dsp, screen);
       }
    }
    
    void XSketchPx (int x, int y, unsigned char px[3]) {
      XLockDisplay (dsp);
      XColor   pX11;
      pX11.red = 256*px[0];
      pX11.green = 256*px[1];;
      pX11.blue = 256*px[2];;        
      if (XAllocColor(dsp,cmap,&pX11)==0)
        assert (0);
      XSetForeground(dsp,gc_2,pX11.pixel);
      XDrawPoint(dsp, win, gc_2, x, cam.ny - y );
      if (tid() == 0)
        XFlush (dsp);
      XUnlockDisplay (dsp);
    }
    
    void XStop (void) {
      XDestroyWindow(dsp, win);
      XCloseDisplay(dsp);
      win = (Window)NULL;
    }
    
    event stop_X11 (t = end) {
      XStop();
    }
    #endif

    Preparation

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

    void preparation_steps (void){
    #if WATCH_ALONG
      XSetup();
    #endif
      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 not s
          }
          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});
          foreach() {
    	coord n = {0};
    	foreach_dimension()
    	  n.x = v.x[];
    	if (any_comp(n)) {
    	  normalize (&n);
    	  foreach_dimension()
    	    v.x[] = n.x;
    	}
          }
          boundary ((scalar*){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});
          }
        }
      }
      if (sketch_a_volume) { //Pre compute cells with volume
        scalar s = vollist.s;
        char pname [99];
        sprintf (pname, "p_for_%s_volume", s.name);
        scalar p = new_scalar (pname);
        s.possible = p;
        find_possible_vol (s, vollist.mval);
      }
    }
    
    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});
          }
      }
      if (sketch_a_volume) {
        smoke.i = 0;
        scalar s = vollist.s;
        s = s.possible;
        delete ({s});
      }
    }

    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();
      
      foreach_ray() {
        unsigned char px[3] = {BGR, BGG, BGB};
        get_color_ray (_r, px);
        fwrite (px, 1, 3, fp);
    #if WATCH_ALONG
        XSketchPx (ii, jj, px);
    #endif
      }
      fflush (fp);
      aftercare();
      return true;
    }
    
    #if _OPENMP
    trace
    bool store_OPENMP (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();
    #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 WATCH_ALONG
    	XSketchPx (ii, jj, px);
    #endif
          }
          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
    	}
        }
      }  
      fflush (fp);
      aftercare();
      return true;
    }
    // Overload calls to store
    #define store store_OPENMP
    #endif //_OPENMP

    Store adaptive

    With raycaster (or this one) installed, we can adaptively sample the scene (lossy compresion). It can be effecient to generate high resolution movies.

    #include <sys/wait.h>
    
    #define ParentRead      read_pipe[0]
    #define ParentWrite     write_pipe[1]
    #define ChildRead       write_pipe[0]
    #define ChildWrite      read_pipe[1]
    #define BGRa (0)
    #define BGGa (25.*(b[0] + b[1]) + 120)
    #define BGBa (50.*b[1] + 120)
    
    typedef struct _sa {
      FILE * fp;  
      int ml;     // max level of output image (cam.nx/ny are ignored)
      double tol; // tolerance
    } _sa;
    
    trace
    bool store_adaptive (struct _sa inp) {
      double tol = 25;
      int ml = 10;
      if (inp.tol)
        tol = inp.tol;
      if (inp.ml)
        ml = inp.ml;
      if (lights[0].type == 0) // No lights
        default_lights();
      preparation_steps();
      // Organize pipe and fork
      int read_pipe[2];
      int write_pipe[2];  
      pipe (read_pipe);
      pipe (write_pipe);
      pid_t p;
      p = fork();
      if (p == 0) { //child calls `raycaster`
        close (ParentRead);
        close (ParentWrite);
        if (system("which ./raycaster > /dev/null 2>&1")) {
          fprintf (stderr, "Error: Adaptive raycaster is not installed\n");
          assert(0);
        } else {
          dup2 (ChildRead, STDIN_FILENO);
          dup2 (ChildWrite, STDOUT_FILENO);
          dup2 (fileno(inp.fp), STDERR_FILENO);
          char m1[9], m2[15];
          sprintf (m1, "%d", ml);
          sprintf (m2, "%g", tol);
          execl ("raycaster", "raycaster", m1, m2, NULL);
        }
      } else { // Parent computes color codes
        close (ChildRead);
        close (ChildWrite);
        int a = sq(1 << (ml - 2)) + 1;
        int al = a;
        double * pos = malloc (2*al*sizeof(double));
        unsigned char * pxs = malloc (3*al*sizeof(unsigned char));
        while (a) {
          read (ParentRead, &a, sizeof(int));
          if (a > al) {
    	al = 2*a;
    	pos = realloc (pos, 2*al*sizeof(double));
    	pxs = realloc (pxs, 3*al*sizeof(unsigned char));
          }
          for (int i = 0; i < a; i++)
          	read (ParentRead, &pos[2*i], 2*sizeof(double));
          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);
    #if _OPENMP
    #pragma omp parallel
    #endif
          for (int i = tid(); i < a; i += npe()) {
    	coord apoint;
    	double b[2];
    	b[0] = pos[2*(i)];
    	b[1] = pos[2*(i) + 1];
    	foreach_dimension() 
    	  apoint.x = (cam.poi.x
    		      + cam.fov*(b[0] - 0.5)*hori.x
    		      + cam.fov*(b[1] - 0.5)*vert.x);  
    	// New ray
    	ray _r;
    	_r.depth = 0;
    	_r.O = cam.O;
    	foreach_dimension()
    	  _r.dir.x = apoint.x - cam.O.x;
    	normalize (&_r.dir);
    	// Get and write color
    	pxs[3*i] = BGRa;
    	pxs[3*i + 1] = BGGa;
    	pxs[3*i + 2] = BGBa;
    	get_color_ray (_r, &pxs[3*i]);
          }
          write (ParentWrite, &pxs[0], 3*a*sizeof(unsigned char));
        }
        wait (NULL); //wait until the child is done writing the  image
        free (pos); pos = NULL;
        free (pxs); pxs = NULL;
      }
      aftercare();
      return true;
    }