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

SCM Repository

[diderot] View of /branches/ein16/synth/d2/symb/symb_fsd3.diderot
ViewVC logotype

View of /branches/ein16/synth/d2/symb/symb_fsd3.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4456 - (download) (annotate)
Thu Aug 25 04:03:27 2016 UTC (2 years, 9 months ago) by cchiw
File size: 5002 byte(s)
added missing cubic
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 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;

//-- field defined by coefficients --
input tensor[4] setA_z0 =[0,0,0,0];
input tensor[4] setB_z0 =[0,0,0,0];
input tensor[4] setC_z0 =[0,0,0,0];
input tensor[4] setD_z0 =[0,0,0,0];

input tensor[4] setA_z1 =[0,0,0,0];
input tensor[4] setB_z1 =[0,0,0,0];
input tensor[4] setC_z1 =[0,0,0,0];
input tensor[4] setD_z1 =[0,0,0,0];

input tensor[4] setA_z2 =[0,0,0,0];
input tensor[4] setB_z2 =[0,0,0,0];
input tensor[4] setC_z2 =[0,0,0,0];
input tensor[4] setD_z2 =[0,0,0,0];


function real cvt(real x, real y, real z, tensor[4] setA, tensor[4] setB, tensor[4] setC, tensor[4] setD) {
    real a = setA[0];
    real b = setA[1];
    real c = setA[2];
    real d = setA[3];
    real e = setB[0];
    real f = setB[1];
    real g = setB[2];
    real h = setB[3];
    real i = setC[0];
    real j = setC[1];
    real k = setC[2];
    real l = setC[3];
    real m = setD[0];
    real n = setD[1];
    real o = setD[2];
    real p = setD[3];

    // as intended to be represented
    real tA = a + b*y + c*x*y+ d*x;
    real tB = ((e+f*x+g*(x*x))*y*y) + (h*(x*x)*y);
    real tC = i*(x*x) + (j+k*x+l*(x*x))*y*y*y;
    real tD = (x*x*x)*((m*y*y*y)+(n*y*y)+(o*y)+p);
    return  z*(tA+tB+tC+tD);
}

function real func(vec3 pos) {
    real x = pos[0];
    real y = pos[1];
    real z = pos[2];
    real z0 = cvt(x,y,1, setA_z0, setB_z0, setC_z0, setD_z0);
    real z1 = cvt(x,y,z, setA_z1, setB_z1, setC_z1, setD_z1);
    real z2 = cvt(x,y,z*z,setA_z2, setB_z2, setC_z2, setD_z2);
    return z0+z1+z2;


}

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