sandbox/Antoonvh/bwatch.h

    Bwatch

    Bwatch is a ray-casting toolbox written in Basilisk C. It should work with Octrees or multigrid3D.

    #include "fractions.h"

    A list of sketch functions is stored. Rendering takes place when an image is store()d.

    #define N_SKETCH 256          // Maximum sketch functions
    typedef double (*sketch)();   // a Sketch function prototype   
    sketch sketch_list[N_SKETCH]; // List sketching functions
    int functot;                  // number of sketching function calls
    scalar * shading = NULL;      // VOF objects that cast shade
    scalar * Smoke   = NULL;      // Concentration fields that cast shade
    //Prim * Obj  = NULL;         // Primitives should cast a shade as well
    int max_ray_depth = 5;        // Number of ray recasts
    attribute {                   // attenuation coefficient of smoke fields  
      double att;                 
    }
    
    // A camera class and "the" cam
    typedef struct Camera {
      coord O;         // camera pos
      coord up;        // upward pointing vector (related to image vertical)
      coord rhs;       // In case up-vector is parallel to the viewing direction
      coord poi;       // Pointing towards
      double fov;      // horizontal width fustrum through poi.
      int nx, ny;      // Resolution
      double min, max; // Sensitivities
      double (*f)();   // Transfer function
    } Camera;
    
    Camera cam = {.O   = {0., 0., 3},   // view from back to front (-z)
    	      .up  = {0., 1., 0.},
    	      .rhs = {1., 0., 0.},  
    	      .poi = {0., 0., 0.},  // view from back to front (-z)
    	      .fov = 1.1,
    	      .nx = 450,
    	      .ny = 400,
    	      .min = 0.1,
    	      .max = 100,
    	      .f = log};

    User function for the Camera setup

    bool watch (struct Camera inp) {
      foreach_dimension() {
        if (inp.O.x)
          cam.O.x = inp.O.x;
        if (inp.up.x)
          cam.up.x = inp.up.x;
        if (inp.rhs.x)
          cam.rhs.x = inp.rhs.x;
        if (inp.poi.x)
          cam.poi.x = inp.poi.x;
      }
      if (inp.fov)
        cam.fov = inp.fov;
      if (inp.nx)
        cam.nx = inp.nx;
      if (inp.ny)
        cam.ny = inp.ny;
      if (inp.min)
        cam.min = inp.min;
      if (inp.max)
        cam.max = inp.max;
      if (inp.f)
        cam.f = inp.f;
      return true;
    }
    
    // A ray class
    typedef struct ray {
      coord O;
      coord dir;
      int   depth;   // for recasted rays
    } ray;
    
    // Lights
    typedef struct light {
      int type;             // ambient (1), Sun (2), point (3), ...
      unsigned char col[3]; // Color
      double I;             // Intensity (at cam.poi)
      coord O;              // Location of point light
      coord dir;            // Direction of sun light
    } light;
    
    light lights[N_SKETCH]; //Reuse N_SKETCH

    Material properties

    Properties for matierial objects

    #include "utils.h"
    typedef struct material {
      unsigned char col[3];   // prescribed color
      scalar s;               // color from field data
      colormap map;           // Corresponding colorbar ...
      vector v;               // RGB from vector
      double min, max;        // ... cbar range
      bool linear;            // Cbar interpolation;
      double SPR;             // Specular albedo   (Blinn)
      double SPexp;           // Specular exponent (Blinn)
      double ind;             // refraction index (normal points outwards)
      double T;               // Transparancy  (0 - 1 intended)
      double R;               // Reflectivity  (0 - 1 intended)
      bool dull;              // Dull object: no light interaction just `col`
      unsigned char col2[3];  // Or a gradient to col2 (not so dull!)
      coord n1;               // normal for `col`, col2 at n = -n1;
    } material;

    A few (vector) helper functions

    void bias (ray * r) {
      foreach_dimension()
        r->O.x += r->dir.x*L0*1e-4;
    }
    
    bool any_col (unsigned char px[3]) {
      if (px[0] || px[1] || px[2])
        return true;
      return false;
    }
    
    bool any_comp (coord a) {
      if (a.x || a.y || a.z)
        return true;
      return false;
    }
    
    coord cross (coord a, coord b) {
      coord new;
      foreach_dimension()
        new.x = a.y*b.z - a.z*b.y;
      return new;
    }
    
    double dot (coord a, coord b) {
      double d = 0;
      foreach_dimension()
        d += a.x*b.x;
      return d;
    }
    
    coord reflect (coord in, coord n) {
      coord out;
      foreach_dimension()
        out.x = in.x - 2*dot (in, n)*n.x;
      return out;
    }
    
    // Thanks to scratch pixel:
    // https://www.scratchapixel.com/lessons/3d-basic-rendering/
    // introduction-to-shading/reflection-refraction-fresnel
    coord refract(coord in, coord normal, double ind) { 
      normalize (&in);
      normalize (&normal);
      double cosi = dot (in, normal); 
      double etai = 1, etat = ind; 
      coord n;
      foreach_dimension()
        n.x = normal.x; 
      if (cosi < 0) 
        cosi = -cosi;
      else {
        etai = etat;
        etat = 1.;
        foreach_dimension()
          n.x = -normal.x;
      }
      double eta = etai/etat; 
      double k = 1. - sq(eta) * (1. - sq(cosi));
      assert (k >= 0);
      coord out = {0.};
      foreach_dimension()
        out.x = k < 0 ? 0 : eta*in.x + (eta*cosi - sqrt(k))*n.x;
      return out; 
    }
    
    double vec_length (coord v) {
      double l = 0;
      foreach_dimension()
        l += sq(v.x);
      return l > 0 ? sqrt(l): 0;
    }
    
    // Thanks you scratchpixel:
    // https://www.scratchapixel.com/lessons/3d-basic-rendering/
    // introduction-to-shading/reflection-refraction-fresnel
    double fresnel(coord in, coord normal, double ind) { 
      normalize (&in);
      normalize (&normal);
      double cosi = dot(in, normal); 
      double etai = 1, etat = ind;
      if (cosi > 0) {
        etai = ind;
        etat = 1.;
      } 
      // Compute sini using Snell's law
      double sint = etai/etat * sqrt(max(0., 1 - sq(cosi))); 
      // Total internal reflection
      if (sint >= 1)
        return 1.;
      else { 
        double cost = sqrt(max(0., 1 - sq(sint))); 
        cosi = fabs(cosi); 
        double Rs = ((etat*cosi) - (etai*cost)) / ((etat*cosi) + (etai*cost)); 
        double Rp = ((etai*cosi) - (etat*cost)) / ((etai*cosi) + (etat*cost)); 
        return (sq(Rs) + sq(Rp))/2; 
      } 
    }

    Intersections

    One may compute

    • ray-cell intersection
    • ray-plane intersection
    • segment-facet intersection
    static inline double ray_cell_intersect (ray r, Point point, coord a[2]) {
      coord cc = {x,y,z};
      double in = -HUGE, out = HUGE;
      foreach_dimension() {
        in  = max(in,  r.dir.x != 0 ? ((cc.x - (sign(r.dir.x)*Delta/2.)) - r.O.x)/r.dir.x : -HUGE);
        out = min(out, r.dir.x != 0 ? ((cc.x + (sign(r.dir.x)*Delta/2.)) - r.O.x)/r.dir.x :  HUGE);
      }
      if (in >= out || out < 0) 
        return HUGE;
      in = in < 0 ? 0 : in; //The origin is in the cell.
      foreach_dimension() {
        a[0].x = r.O.x +  in*r.dir.x;
        a[1].x = r.O.x + out*r.dir.x;
      }
      return in;
    }
    
    static inline double ray_plane_intersect (ray r, coord n, double alpha, coord * a) {
      double ldotn = 0;
      foreach_dimension()
        ldotn += n.x*r.dir.x;
      if (ldotn == 0) //parallel
        return HUGE;
      double d = alpha;
      foreach_dimension()
        d -= n.x*r.O.x;
      d /= ldotn;
      if (d <= 0) //plane certainly not in view.
        return -1;
      d = d > 1e4*L0 ? 1e4*L0 : d;
      foreach_dimension() {
        a[0].x = r.O.x + d*r.dir.x;
      }
      return d;
    }
    
    static inline double ray_sphere_intersect (ray r, coord C, double R, coord * a, coord * n) {
      normalize (&r.dir);
      coord oc;
      foreach_dimension()
        oc.x = r.O.x - C.x;
      double det = (sq(dot(r.dir, oc)) - (sq(vec_length (oc)) - sq(R)));
      if (det < 0)
        return HUGE;
      det = sqrt(det);
      double dist  = -dot(r.dir, oc) - det;
      double dist2 = -dot(r.dir, oc) + det;
      if (dist < 0 && dist2 < 0)
        return HUGE;
      else if (dist < 0)
        dist = dist2;
      else if (dist > dist2)
        dist = dist2;
      foreach_dimension() {
        a->x = r.O.x + r.dir.x*dist;
        n->x = a->x - C.x;
      }
      normalize (n);
      return dist;
    }
    
    // use ray_plane_intersect?
    static inline bool segment_facet_intersect (coord a[2], scalar f, Point point, coord * n) {
      coord cc = {x, y, z};
      n[0]  = mycs (point, f); //pointing outwards.
      double alpha = plane_alpha (f[], n[0]);
      double ALP = 0, ALP2 = 0;
      foreach_dimension() {
        ALP  += n[0].x*(a[0].x - cc.x)/Delta;
        ALP2 += n[0].x*(a[1].x - cc.x)/Delta;
      }
      if ((ALP2 - alpha)/(ALP - alpha) > 0.05) // 3% 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;
    }

    Iterators

    see,

    attribute {
      scalar possible;
      vector center;
      vector normals;
      scalar vofd; 
    }
    
    #include "bwatch-iterators.h"

    Color from colorbar

    Quantative data can be mapped to an RGB-color code via a colorbar. The function below is a version of colormap_color() in output.h.

    #include "utils.h"
    
    bool colormap_pigmentation (unsigned char c[3], double cmap[NCMAP][3],
    			    double val, double min, double max) {
      if (val == nodata) {
        c[0] = c[1] = c[2] = 0; // nodata is black
        return false;
      }
      int i;
      double coef;
      if (max != min)
        val = (val - min)/(max - min);
      else
        val = 0.;
      if (val <= 0.) i = 0, coef = 0.;
      else if (val >= 1.) i = NCMAP - 2, coef = 1.;
      else {
        i = val*(NCMAP - 1);
        coef = val*(NCMAP - 1) - i;
      }
      assert (i < NCMAP - 1);
      for (int j = 0; j < 3; j++)
        c[j] = 255*(cmap[i][j]*(1. - coef) + cmap[i + 1][j]*coef);
      return true;
     }

    Lightning

    Lights can be shaded by vof facets (in shading) and should be reflected by reflectors (in rlist)

    bool shaded (coord start, light source) {
      ray l;
      l.O = start;
      foreach_dimension()
        l.dir.x = -source.dir.x;
      for (scalar c in shading) 
        foreach_ray_facet_intersection(l, c)
          return true;
      return false;
    }
    
    double Intensity_by_reflection (coord start, light source) {
      return 0;
    }
    
    void default_lights (void) {
      // Ambient and sun;
      lights[0].type = 1;
      lights[0].I = 2;
      lights[0].col[0] = 255;
      lights[0].col[1] = 255;
      lights[0].col[2] = 255;
      
      lights[1].type = 2;
      lights[1].I = 80;
      lights[1].dir = (coord){-1, -1.8, -2};
      lights[1].col[0] = 255;
      lights[1].col[1] = 255;
      lights[1].col[2] = 250;
    }
    
    double diffuse_Lambert (coord n, coord a) {
      double I = 0;
      int l = 0; //light
      while (lights[l].type > 0) { //Fixme: No mirrored light sources yet
        if (lights[l].type > 1) {
          double Il = 0;
          coord ldir = {0,0,0};
          if (lights[l].type == 2)
    	ldir = lights[l].dir;
          if (!shaded (a, lights[l])) {
    	normalize (&ldir);
    	foreach_dimension()
    	  Il -= n.x * ldir.x;
    	I += max (0, lights[l].I*Il);
          }
        }
        l++;
      }
      return I;
    }
    
    bool specular (ray r, coord n, coord a, double gloss, double alR, unsigned char c[3]) {
      int l = 0;
      normalize (&n);
      while (lights[l].type > 0) { //Fixme: No mirrored light sources yet
        if (lights[l].type > 1) {
          if (lights[l].type == 2) { //Fixme: Only sun
    	if (shaded (a, lights[l]))
    	  return false;
    	coord R = reflect (lights[l].dir, n);
    	normalize (&R);
    	double RdN = -dot(r.dir, R);
    	if (RdN < 0)
    	  return false;
    	double I = max(alR*lights[l].I*pow(RdN, gloss), 1e-6);
    	for (int i = 0; i < 3; i++) {
    	  c[i] = lights[l].col[i]* min(max(log(I) - log(cam.min),0)/log(cam.max/cam.min), 1);
    	  c[i] = min (c[i], 255);
    	}
          }
        }
        l++;
      }
      return true;
    }
    
    double light_intensity (coord n, coord a) {
      double I = diffuse_Lambert(n, a);
      int l = 0;
      //ambients
      while (lights[l].type > 0) {
        if (lights[l].type == 1) 
          I += lights[l].I;
        l++;
      }
      I = min (max(cam.f(I) - cam.f(cam.min),0)/cam.f(cam.max/cam.min), 1);
      return I;  
    }
    
    struct _get_pixel{
      unsigned char c[3];
      coord n;
      coord a;
      ray r;
      double SPexp;
      double SPR;
    };
      
    void get_pixel (struct _get_pixel * inp) {
      double I = light_intensity (inp->n, inp->a);
      for (int i = 0; i < 3; i++)
        inp->c[i] *= I;
      if (inp->SPexp && inp->SPR) {
        unsigned char tmp[3];
        if (specular (inp->r, inp->n , inp->a, inp->SPexp, inp->SPR, tmp))
          for (int i = 0; i < 3; i++) 
    	inp->c[i] = min(inp->c[i] + tmp[i], 255);
      }
    }

    Pixel color from a material hit (maybe recursive)

    Mind the presidence certain properties take over others:

    1. Dull
    2. Refractive
    3. Fully reflective (R = 1)
    4. Prescribed color
    5. Color from scalar field data and color bar
    6. RGB Color code from vector field (not implemented atm)

    4, 5 and 6 maybe supplemented with either 7 or 8,

    1. Reflection (R < 1)
    2. Transmission (transparancy, without refraction)
    double get_color_ray(); //Prototype
    
    trace
    bool get_color_material (ray r, coord a, coord n, material mat, unsigned char px[3]) {
      // Dull object?
      if (mat.dull) {
        if (any_comp(mat.n1) && any_comp(n)) {
          normalize (&n); normalize (&mat.n1);
          double l = (dot (n, mat.n1) + 1.)/2.;
          for (int i = 0; i < 3; i++)
    	px[i] = l*mat.col[i] + (1. - l)*mat.col2[i];
          return true;
        }
        for (int i = 0; i < 3; i++)
          px[i] = mat.col[i];
        return true;
      }
      // Refractive object?
      else if (mat.ind) {
        if (r.depth <= max_ray_depth) {
          normalize (&n);  normalize (&r.dir);
          ray refl, refr;
          refl.depth = refr.depth = r.depth + 1;
          refl.O = refr.O = a;
          refl.dir = reflect (r.dir, n);
          normalize (&refl.dir);
          double w = fresnel (r.dir, n, mat.ind);
          if (w < 1) {
    	refr.dir = refract (r.dir, n, mat.ind);
    	normalize (&refr.dir);
          }
          bias (&refl); bias (&refr);
          unsigned char ptrefl[3], ptrefr[3];
          get_color_ray (refl, ptrefl);
          if (w < 1) 
    	get_color_ray (refr, ptrefr);
          else
    	w = 1;
          for (int i = 0; i < 3; i++)
    	px[i] = w*ptrefl[i] + (1 - w)*ptrefr[i];
          if (mat.SPexp && mat.SPR) {
    	unsigned char tmp[3];
    	if (specular (r, n , a, mat.SPexp, mat.SPR, tmp))
    	  for (int i = 0; i < 3; i++) 
    	    px[i] = min(px[i] + tmp[i], 255);
          }
          return true;
        }
        return false;
      }
      // Full Reflector?
      if (mat.R >= 1) {
        if (r.depth <= max_ray_depth) {
          ray new = {.O = a, .dir = reflect (r.dir, n), .depth = r.depth + 1};
          bias (&new);
          if ((get_color_ray (new, px) == HUGE))
    	for (int i = 0; i < 3; i++)
    	  px[i] = mat.col[i];
          return true;
        }
        return false;
      }
      // It must be a "tradional" object:
      assert (any_col (mat.col) || mat.s.i || mat.v.y.i || mat.R >= 1 );
      struct _get_pixel out;
      out.r = r;
      out.SPexp = mat.SPexp;
      out.SPR   = mat.SPR;
      out.n = n;
      out.a = a;
      if (any_col (mat.col))
        for (int i = 0; i < 3; i++)
          out.c[i] = mat.col[i];
      else if (mat.s.i) {
        double cmap[NCMAP][3];
        scalar s = mat.s;
        mat.map (cmap);
        foreach_dimension()
          a.x += L0/(1<<depth())*1e-5*r.dir.x;
        Point point = locate (a.x, a.y, a.z);
        if (point.level < 0)
          return false;
        double val = mat.linear ?
          interpolate (s, a.x, a.y, a.z) : s[];
        colormap_pigmentation (out.c, cmap, val, mat.min, mat.max);
      }
      else if (mat.v.y.i) {
        vector v = mat.v;
        foreach_dimension()
          a.x += L0/(1<<depth())*1e-5*r.dir.x;
        Point point = locate (a.x, a.y, a.z);
        if (point.level < 0)
          return false;
        double val = mat.linear ?
          interpolate (v.x, a.x, a.y, a.z) : v.x[];
        if (val == nodata)
          return false;
        out.c[0] =  255*min(max(val - mat.min, 0.)/(mat.max - mat.min), 1.);
        val = mat.linear ?
          interpolate (v.y, a.x, a.y, a.z) : v.y[];
        if (val == nodata)
          return false;
        out.c[1] =  255*min(max(val - mat.min, 0.)/(mat.max - mat.min), 1.);
        val = mat.linear ?
          interpolate (v.z, a.x, a.y, a.z) : v.z[];
        if (val == nodata)
          return false;
        out.c[2] =  255*min(max(val - mat.min, 0.)/(mat.max - mat.min), 1.);
      }
      get_pixel (&out);
      for (int i = 0; i < 3; i++)
        px[i] = out.c[i];
      if ((!mat.R && !mat.T) || r.depth > max_ray_depth)
        return true;
      if (mat.R) {
        unsigned char temp[3];
        ray new = {.O = a, .dir = reflect (r.dir, n), .depth = r.depth + 1};
        bias (&new);
        if ((get_color_ray (new, temp) != HUGE))
          for (int i = 0; i < 3; i++)
    	px[i] = (1 - mat.R)*px[i] + mat.R*temp[i];
        return true;
      }
      if (mat.T) {
        unsigned char temp[3];
        ray new = {.O = a, .dir = r.dir, .depth = r.depth + 1};
        bias (&new);
        if ((get_color_ray (new, temp) != HUGE))
          for (int i = 0; i < 3; i++)
    	px[i] = (1 - mat.T)*px[i] + mat.T*temp[i];
        return true;
      }
      
      //this should not happen:  
      assert (0);
      for (int i = 0; i < 3; i++)
        px[i] = 0; // Darkness
      return false;
    }

    User interface

    The sketch and store functions are elsewhere:

    #include "sketch.h"

    Tests

    Usage

    All pages using bwatch

    to do

    1. A flexible camera
    2. MG-accelerated foreach_segment_3D() ray-casting
    3. Do something ray-ish
    4. implement and combine more than one sketch function
    5. Implement a fun sketch function
    6. Consider lights and colours
    7. Reflection of light sources
    8. Volumetric rendering
    9. Use bwatch for something
    10. Think about optimizations to mitigate the long rendering times
    11. Fill gaps in vof facets
    12. Refractive objects
    13. Implement point-light sources
    14. Primitives cast shades
    15. Partially reflective and transparant (alpha buffering) objects
    16. Gouraud shading for facets
    17. Isosurfaces
    18. RGB-code from vector field
    19. Parallel acceleration with OpenMP