Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] View of /tests/examples/fs3d/fs3d-scl.diderot
ViewVC logotype

View of /tests/examples/fs3d/fs3d-scl.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4640 - (download) (annotate)
Tue Sep 27 20:54:47 2016 UTC (2 years, 9 months ago) by glk
File size: 14102 byte(s)
initial result of svn export --username anonsvn --password=anonsvn https://svn.smlnj-gforge.cs.uchicago.edu/svn/diderot/branches/vis15/src/tests/
/* ==========================================
## fs3d-scl.diderot: function sampler, 3D scalar synthetic fields

`fs3d-scl.diderot` is based on [`fs2d-scl.diderot`](../fs2d); some
copy-and-pasting was involved.  Like that program, the value of this
is not so much as a model of what a good or typical Diderot program
looks like, but as a utility for generating datasets for other Diderot
programs to work on.  This same directory has a `fs3d-vec.diderot` for
generating vector fields, and possibly other programs; the
documentation below is for `fs3d-scl.diderot` though the other
programs will have very similar structure and usage.

This programs generates synthetic scalar 3D data on regular grids that
are located and oriented in world space.  The input arguments here
make it easy to sample the same underlying function on different
grids.  Grid orientation is specified by the angle and axis of
rotation.

Assuming the directions at https://github.com/Diderot-Language/examples
this program can be compiled with:

	../../vis12/bin/diderotc --exec fs3d-scl.diderot

The `-which` option will determine which function is sampled; look for
`(0 == which)` in the code to see the start of the function
definitions, this includes 3D functions used to make datasets for
[some](http://people.cs.uchicago.edu/~glk/pubs/#VIS-2003)
[previous](http://people.cs.uchicago.edu/~glk/pubs/#VIS-2009)
[papers](http://people.cs.uchicago.edu/~glk/pubs/#VIS-2014).

This program is unusual in that its printed output needs to be captured
in order to have a NRRD header that records the orientation of the
sampling grid, so using the program involves redirection.  To
get a self-contained parab.nrrd containing a parabola function

	./fs2d-scl -which 3 | unu save -f nrrd -o parab.nrrd
	rm out.nrrd

Note that like in [`fs2d-scl.diderot`](../fs2d), the NRRD header
generated by this program assumes:

1. This program was not compiled with `--double`
2. The program is running on a little-endian machine.

========================================== */
input int sz0 ("# samples on fastest axis") = 52;
input int sz1 ("# samples on medium axis") = 51;
input int sz2 ("# samples on slowest axis") = 50;
input real width ("width of edge of cube region sampled") = 4;
input vec3 axis ("axis (non-normalized) of rotation of sampling grid") = [1,1,1];
input real angle ("angle (in degrees) of rotation of sampling grid") = 0;
input vec3 off ("translation offset, in index space, from origin-centered grid") = [0,0,0];
input int which ("which function to sample (from 0 to 13)") = 0;
input vec4 parm ("parameters that functions may use") = [0,0,0,0];

real theta = angle*π/180;
// unit-length axis of rotation scaled by sin(theta/2)
vec3 snax = sin(theta/2)*normalize(axis);
// quaternion representing rotation
vec4 qq = [cos(theta/2), snax[0], snax[1], snax[2]];
// create rotation matrix by converstion from quaternion
// (which aren't currently natively supported in Diderot)
tensor[3,3] rot = [[qq[0]*qq[0] + qq[1]*qq[1] - qq[2]*qq[2] - qq[3]*qq[3],
                    2*(qq[1]*qq[2] - qq[0]*qq[3]),
                    2*(qq[1]*qq[3] + qq[0]*qq[2])],
                   [2*(qq[1]*qq[2] + qq[0]*qq[3]),
                    qq[0]*qq[0] - qq[1]*qq[1] + qq[2]*qq[2] - qq[3]*qq[3],
                    2*(qq[2]*qq[3] - qq[0]*qq[1])],
                   [2*(qq[1]*qq[3] - qq[0]*qq[2]),
                    2*(qq[2]*qq[3] + qq[0]*qq[1]),
                    qq[0]*qq[0] - qq[1]*qq[1] - qq[2]*qq[2] + qq[3]*qq[3]]];
// per axis inter-sample edge vectors
real space0 = width/(sz0-1);
real space1 = width/(sz1-1);
real space2 = width/(sz2-1);
vec3 edge0 = rot•[space0, 0, 0];
vec3 edge1 = rot•[0, space1, 0];
vec3 edge2 = rot•[0, 0, space2];
// offset in world space
vec3 offws = [off[0]*space0, off[1]*space1, off[2]*space2];
// location of first sample point
vec3 orig = -(edge0*(sz0-1) + edge1*(sz1-1) + edge2*(sz2-1))/2 + offws;

/* Note how a function (like gauss and ball below) may defined by
   equality with a single expression */

// univariate normal distribution
function real gauss(real stdv, real x) =
   exp(-(x*x)/(2*stdv^2))/(stdv*sqrt(2*π)); // have some π

// distance inside sphere of radius r
function real ball(real r, vec3 pos) = max(0, r - |pos|);

// radius R code pointing along positive x
function real cone(real R, vec3 pos) {
   real x = pos[0];
   real y = pos[1];
   real z = pos[2];
   real ret = max(0, R - (sqrt(y^2 + z^2) + x/2));
   ret *= 1 if (x > 0) else 0;
   return ret;
}

// trapezoid function at b in [a,c], with ramp R
function real trpz(real R, real a, real b, real c, real scl) {
   real ww = (b-a)/(c-a);
   real rr = R + 0.000001;
   /* Note the use of Python syntax for conditional expressions;
      chaining them like this may look a little weird at first
      but note that it can sort of be read like a decision tree,
      with the two sides branching off the central "if" part */
   return (0.0
           if ww < -rr
           else (scl*(1 + ww/rr)
                 if ww < 0
                 else (lerp(scl, 1, ww)
                       if ww < 1
                       else (1 - (ww-1)/rr
                             if ww < 1+rr
                             else 0.0))));
}

// rounds v to nearest integer
function real round(real v) = (v+0.5 - fmod(v+0.5, 1)
                               if (v > 0) else
                               -(-v+0.5 - fmod(-v+0.5, 1)));

function real func(vec3 pos) {
   real x = pos[0];
   real y = pos[1];
   real z = pos[2];
   real ret = 0;
   if (0 == which) {        // 0: x ramp
      ret = x;
   } else if (1 == which) { // 1: y ramp
      ret = y;
   } else if (2 == which) { // 2: z ramp
      ret = z;
   } else if (3 == which) { // 3: upwards parabola
      ret = pos•pos;
   } else if (4 == which) { // 4: downwards parabola
      ret = 1 - pos•pos;
   } else if (5 == which) { // 5: scaled downwards parabola
      real sx = 1 + parm[0];
      real sy = 1 + parm[1];
      real sz = 1 + parm[2];
      ret = 1 - (sx*x^2 + sy*y^2 + sz*z^2);
   } else if (6 == which) { // 6: tunable torus
      /* attains max value on circle (minor radius=0) parm[0] sets
         major radius (minus 1); linear ramp parm[1]*x/4 is added on;
         the torus amplitude (scaling) is controlled by parm[2] and
         parm[3]. Note that the intensity cross-section across the
         core is more like a parabola x^2 than |x|; this is to make
         the core circle into a nicer ridge feature */
      real rmaj = parm[0] + 1;
      real rr = sqrt(x^2 + y^2);
      real win = parm[3] + 1;
      real ang = clamp(0, π, lerp(0, π, -rmaj*win, y, rmaj*win));
      real scl = 1 + parm[2]*lerp(0, 1, 1, cos(ang), -1);
      ret = 1 + parm[1]*x/4 - scl*(z^2 + (rr-rmaj)^2);
   } else if (7 == which) { // 7: rounded-cube quartic
      /* Fig 1 of http://people.cs.uchicago.edu/~glk/pubs/#VIS-2003
         used data generated by this function (in a C program); that
         was wholly based on http://dl.acm.org/citation.cfm?id=134082
         This is also http://mathworld.wolfram.com/GoursatsSurface.html
         with (a,b,c) = (0,-1,0), though negated so higher values are
         closer to the origin. The usual isosurface (at isovalue 0)
         fits comfortably inside width 3 */
      ret = x^2 + y^2 + z^2 - x^4 - y^4 - z^4;
   } else if (8 == which) { // 8: Möbius strip
      /* Fig 4 of http://people.cs.uchicago.edu/~glk/pubs/#VIS-2009
         used data generated by this function (in a C program) */
      real dd = parm[2] + 0.85;
      real Rbig = sqrt(x*x + y*y);
      real Rlit = sqrt(z*z + (Rbig-dd)^2);
      real phi = atan2(Rbig-dd, z) - atan2(x, y)/2;
      real a = Rlit*cos(phi);
      real b = Rlit*sin(phi);
      ret = gauss(0.35+parm[0], a)*gauss(0.1+parm[1], b);
   } else if (9 == which) { // 9: three twisty lobes
      /* some methods in http://people.cs.uchicago.edu/~glk/pubs/#VIS-2013
         were debugged with the aid of this dataset (via particles sampling
         its ridge surfaces), though it didn't appear in the paper. The
         isosurface at isovalue 0.5 fits comfortably inside width 2. */
      real mask = lerp(1, 0, -1, erf((|pos| - 0.75)*15), 1);
      real R2 = sqrt(x^2 + y^2);
      real twist = parm[0] + 1.3;
      real phi = atan2(y,x) + z*twist;
      // sorcery here. lobe sharpness 60 matches mask sharpness 15
      ret = mask*pow((1+cos(3*phi))/2, max(0, 60*R2^2 - 0.03));
   } else if (10 == which) { // 10: Cayley cubic surface
      /* see http://mathworld.wolfram.com/CayleyCubic.html
         Fig 2(c) of http://people.cs.uchicago.edu/~glk/pubs/#VIS-2014
         used data generated by this function (in a Diderot program!);
         it assumes width ~= 6. In that paper opacity was limited to
         a sphere of radius 2.7 */
      ret = -(x^2 + y^2 - z*x^2 + z*y^2 + z^2 - 1);
   } else if (11 == which) { // 11: sphere with ramp
      /* Talk for http://people.cs.uchicago.edu/~glk/pubs/#VIS-2015
         used this dataset to show off Canny edge rendering */
      real rad = 1 + parm[0];
      real sharp = 7 + parm[1];
      real ramp = 2 + parm[2];
      ret = ramp*pos[0] + erf(sharp*(rad-|pos|));
   } else if (12 == which) { // 12: little coordinate frame
      /* this is useful for debugging cameras in volume rendering;
         it assumes width >= 3, though much of it is constant zero
         (it doesn't set much in the negative X, Y, or Z half-spaces) */
      // big ball at origin, so local maxima at [0,0,0]
      ret = max(ret, ball(0.3, pos));
      // cylinder along X axis
      ret = max(ret, ball(0.2,  [0,y,z])*trpz(0.1,0,x,1,1));
      // cone at end of X axis w/ local maxima at [1,0,0]
      ret = max(ret, cone(0.25, pos - [1,0,0]));
      // conical-ish thing along Y axis w/ local maxima at [0,1,0]
      ret = max(ret, ball(0.25, [x,0,z])*trpz(0.1,0,y,1,0.5));
      // short rod along Z axis
      ret = max(ret, ball(0.2,  [x,y,0])*trpz(0.1,0,z,0.7,1));
      // ball centered at [0,0,1]
      ret = max(ret, ball(0.25, pos - [0,0,1]));
      // so that isovalue at 1 looks as intended
      ret *= 10;
   } else if (13 == which) { // 13: cube frame
      /* has local maxima at (+-1,+-1,+-1), should use at least width 3,
         but may need a larger width in case of small sampling grid, and
         larger rendering kernels */
      ret = max(ret, ball(0.3, pos-[ 1, 1, 1]));
      ret = max(ret, ball(0.3, pos-[-1, 1, 1]));
      ret = max(ret, ball(0.3, pos-[ 1,-1, 1]));
      ret = max(ret, ball(0.3, pos-[-1,-1, 1]));
      ret = max(ret, ball(0.3, pos-[ 1, 1,-1]));
      ret = max(ret, ball(0.3, pos-[-1, 1,-1]));
      ret = max(ret, ball(0.3, pos-[ 1,-1,-1]));
      ret = max(ret, ball(0.3, pos-[-1,-1,-1]));
      ret = max(ret, ball(0.15, [0,y,z]-[ 0, 1, 1])*trpz(0,-1,x,1,1));
      ret = max(ret, ball(0.15, [0,y,z]-[ 0,-1, 1])*trpz(0,-1,x,1,1));
      ret = max(ret, ball(0.15, [0,y,z]-[ 0, 1,-1])*trpz(0,-1,x,1,1));
      ret = max(ret, ball(0.15, [0,y,z]-[ 0,-1,-1])*trpz(0,-1,x,1,1));
      ret = max(ret, ball(0.15, [x,0,z]-[ 1, 0, 1])*trpz(0,-1,y,1,1));
      ret = max(ret, ball(0.15, [x,0,z]-[-1, 0, 1])*trpz(0,-1,y,1,1));
      ret = max(ret, ball(0.15, [x,0,z]-[ 1, 0,-1])*trpz(0,-1,y,1,1));
      ret = max(ret, ball(0.15, [x,0,z]-[-1, 0,-1])*trpz(0,-1,y,1,1));
      ret = max(ret, ball(0.15, [x,y,0]-[ 1, 1, 0])*trpz(0,-1,z,1,1));
      ret = max(ret, ball(0.15, [x,y,0]-[-1, 1, 0])*trpz(0,-1,z,1,1));
      ret = max(ret, ball(0.15, [x,y,0]-[ 1,-1, 0])*trpz(0,-1,z,1,1));
      ret = max(ret, ball(0.15, [x,y,0]-[-1,-1, 0])*trpz(0,-1,z,1,1));
   } else if (14 == which) { // 14: quadratic patches
      /* a similar function (with width 5) was used for Fig 8(a)
         http://people.cs.uchicago.edu/~glk/pubs/#VIS-2003 */
      real xi = round(x);
      real yi = round(y);
      real xf = x - xi;
      real yf = y - yi;
      real bord = parm[0]+1.0001;
      real marg = parm[1]*0.75+0.0001;
      // trapezoid function at b in [a,c], with ramp R
      real scl = sqrt(trpz(marg, -0.5+marg, bord*xf, 0.5-marg, 1) *
                      trpz(marg, -0.5+marg, bord*yf, 0.5-marg, 1));
      if (yi > xi) {
         scl = 0;
      }
      ret = z + lerp(-width/4, xi*xf^2 + yi*yf^2 - (xi + yi)*0.1, scl);
   } else {
      // update "input int which" annotation above as cases are added
      print("Sorry, no function defined for which = ", which, "\n");
   }
   return ret;
}

strand sample(int idx0, int idx1, int idx2) {
   output real out = 0.0;
   update {
      /* see comment in ../fs2d/fs2d-scl.diderot about the need for these
         conditionals around the print statements */
      if (0 == idx0 && 0 == idx1 && 0 == idx2) {
         print("NRRD0004\n");
         print("# Complete NRRD file format specification at:\n");
         print("# http://teem.sourceforge.net/nrrd/format.html\n");
         /* NOTE: this assumes we haven't been compiled with --double,
            and there isn't currently a way for the program to learn this
            (which in our experience has not been a problem) */
         print("type: float\n");
         print("dimension: 3\n");
         print("sizes: ", sz0, " ", sz1, " ", sz2, "\n");
         print("kinds: space space space\n");
         print("centers: cell cell cell\n");
         // NOTE: this assumes machine endianness
         print("endian: little\n");
         print("encoding: raw\n");
         print("space dimension: 3\n");
         print("space directions: (", edge0[0], ",", edge0[1], ",", edge0[2],
                               ") (", edge1[0], ",", edge1[1], ",", edge1[2],
                               ") (", edge2[0], ",", edge2[1], ",", edge2[2], ")\n");
         print("space origin: (", orig[0], ",", orig[1], ",", orig[2], ")\n");
         print("data file: out.nrrd\n");
         print("byte skip: -1\n");
      }
      out = func(orig + idx0*edge0 + idx1*edge1 + idx2*edge2);
      stabilize;
   }
}
initially [ sample(idx0, idx1, idx2)
            | idx2 in 0..(sz2-1),  // slowest axis
              idx1 in 0..(sz1-1),  // medium axis
              idx0 in 0..(sz0-1)]; // fastest axis

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0