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_v2d3.diderot
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4236 - (download) (annotate)
Wed Jul 20 03:02:00 2016 UTC (3 years, 6 months ago) by cchiw
File size: 10205 byte(s)
added generic cases for trace|det and added test cases
/* fs3d-vec.diderot: function sampler, 3D vector synthetic fields

This is based on fs3d-scl.diderot in this same directory, after much
copy-and-pasting.  As with that program, the `-which` option will
determine which function is sampled; look for `(0 == which)` in the
code to see the start of the function definitions.

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

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

This program's printed output needs to be captured to generate
NRRD header with full orientation info. A vector-valued identity
function is sampled via:

	./fs2d-vec -which 0 | unu save -f nrrd -o ident.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 2)") = 10;
input vec3 parm0 ("parameters that functions may use") = [0,0,0];
input vec3 parm1 ("more parameters that functions may use") = [0,0,0];
input vec3 parm2 ("even more parameters that functions may use") = [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 --
// z=1
input vec2 base1_z0 = [5, 0];    // b*x, c*y
input vec2 xsq1_z0 = [0, 0];     // d*x^2, g*y*x^2
input vec2 ysq1_z0 = [0, 0];     // f*y^2, h*x*y^2
input vec3 diag1_z0 = [0, 7, 0];    // a, e*x*y,  i*y^2*x^2
// z=z
input vec2 base1_z1 = [0, 0];    // b*x, c*y
input vec2 xsq1_z1 = [0, 0];     // d*x^2, g*y*x^2
input vec2 ysq1_z1 = [8, 0];     // f*y^2, h*x*y^2
input vec3 diag1_z1 = [0, 0, 0];    // a, e*x*y,  i*y^2*x^2
// z=z*z
input vec2 base1_z2 = [0, 0];    // b*x, c*y
input vec2 xsq1_z2 = [0, 0];     // d*x^2, g*y*x^2
input vec2 ysq1_z2 = [3, 0];     // f*y^2, h*x*y^2
input vec3 diag1_z2 = [0, 0, 0];    // a, e*x*y,  i*y^2*x^2

//-- field defined by coefficients --
// z=1
input vec2 base2_z0 = [0, 0];    // b*x, c*y
input vec2 xsq2_z0 = [0, 2];     // d*x^2, g*y*x^2
input vec2 ysq2_z0 = [0, 0];     // f*y^2, h*x*y^2
input vec3 diag2_z0 = [0, 0, 5];    // a, e*x*y,  i*y^2*x^2
// z=z
input vec2 base2_z1 = [0, 0];    // b*x, c*y
input vec2 xsq2_z1 = [0, 0];     // d*x^2, g*y*x^2
input vec2 ysq2_z1 = [0, 0];     // f*y^2, h*x*y^2
input vec3 diag2_z1 = [1, 0, 0];    // a, e*x*y,  i*y^2*x^2
// z=z*z
input vec2 base2_z2 = [0, 0];    // b*x, c*y
input vec2 xsq2_z2 = [0, 0];     // d*x^2, g*y*x^2
input vec2 ysq2_z2 = [0, 0];     // f*y^2, h*x*y^2
input vec3 diag2_z2 = [0, 0, 9];    // a, e*x*y,  i*y^2*x^2


function real cvt(real x, real y, real z, vec2 base, vec2 xsq, vec2 ysq, vec3 diag){
    real b = base[0];
    real c = base[1];
    real d = xsq[0];
    real g = xsq[1];
    real f = ysq[0];
    real h = ysq[1];
    real a = diag[0];
    real e = diag[1];
    real i = diag[2];
    real t0 = b*x + c*y;
    real t1 = d*x^2 + g*y*x^2;
    real t2 = f*y^2 + h*x*y^2;
    real t3 = a + e*x*y + i*(y^2)*(x^2);
    return z*(t0+t1+t2+t3);
}

function real grad_x(real x, real y, real z, vec2 base, vec2 xsq, vec2 ysq, vec3 diag){
    real b = base[0];
    real c = base[1];
    real d = xsq[0];
    real g = xsq[1];
    real f = ysq[0];
    real h = ysq[1];
    real a = diag[0];
    real e = diag[1];
    real i = diag[2];
    real t0 = b + d*2*x + g*y*x*2+h*y^2+ e*y + i*(y^2)*(x*2);
    return z*t0;
}
function real grad_y(real x, real y, real z, vec2 base, vec2 xsq, vec2 ysq, vec3 diag){
    real b = base[0];
    real c = base[1];
    real d = xsq[0];
    real g = xsq[1];
    real f = ysq[0];
    real h = ysq[1];
    real a = diag[0];
    real e = diag[1];
    real i = diag[2];
    real t0 = c+ g*x^2+f*y*2 + h*x*y*2+  e*x + i*(y*2)*(x^2);
    return z*t0;
}


function real grad_z(real x, real y, real z, vec2 base, vec2 xsq, vec2 ysq, vec3 diag){
    real b = base[0];
    real c = base[1];
    real d = xsq[0];
    real g = xsq[1];
    real f = ysq[0];
    real h = ysq[1];
    real a = diag[0];
    real e = diag[1];
    real i = diag[2];
    real t0 = b*x + c*y;
    real t1 = d*x^2 + g*y*x^2;
    real t2 = f*y^2 + h*x*y^2;
    real t3 = a + e*x*y + i*(y^2)*(x^2);
    real cat =(t0+t1+t2+t3);
    real ret =  z*cat;
    //print ("\n\n\n b:",b,"c",c,"\nd",d,"g:",g,"f",f,"h",h,"\na:",a,"e",e,"i",i);
    //print ("\nt0: ",t0," t1: ",t1," t2: ",t2," t3: ",t3);
    //print ("\ncat: ",cat, "ret",ret);
    return ret;
}


function vec2 func(vec3 pos) {
   real x = pos[0];
   real y = pos[1];
   real z = pos[2];
   vec2 ret = [0,0];
   if (0 == which) {
       // x = 0;
       // y = 0.1;
       // z = 0.2;
        real axis1_z0 = cvt(x,y,1,base1_z0,xsq1_z0,ysq1_z0,diag1_z0);
        real axis1_z1 = cvt(x,y,z,base1_z1,xsq1_z1,ysq1_z1,diag1_z1);
        real axis1_z2 = cvt(x,y,z*z,base1_z2,xsq1_z2,ysq1_z2,diag1_z2);
        real axis1_val = axis1_z0+axis1_z1+axis1_z2;
        real axis2_z0 = cvt(x,y,1,base2_z0,xsq2_z0,ysq2_z0,diag2_z0);
        real axis2_z1 = cvt(x,y,z,base2_z1,xsq2_z1,ysq2_z1,diag2_z1);
        real axis2_z2 = cvt(x,y,z*z,base2_z2,xsq2_z2,ysq2_z2,diag2_z2);
        real axis2_val = axis2_z0+axis2_z1+axis2_z2;
        ret =[axis1_val, axis2_val];
        print ("\n reg pos:", x,"_",y,"_",z, " ret:",ret);
    }   else if (1 == which) {
        //created to print out gradient of specific position
        x = 0;
        y = 0.1;
        z = 0.2;
       // print ("\n *********************************");
        real axis1_z0 =  grad_x(x,y,0,base1_z0,xsq1_z0,ysq1_z0,diag1_z0);
        real axis1_z1 =  grad_x(x,y,z,base1_z1,xsq1_z1,ysq1_z1,diag1_z1);
        real axis1_z2 = grad_x(x,y,z*z,base1_z2,xsq1_z2,ysq1_z2,diag1_z2);
        real axis1_val = axis1_z0+axis1_z1+axis1_z2;
       // print ("\n axis1_val:",axis1_val);
        real axis2_z0 = grad_y(x,y,0,base1_z0,xsq1_z0,ysq1_z0,diag1_z0);
        real axis2_z1 = grad_y(x,y,z,base1_z1,xsq1_z1,ysq1_z1,diag1_z1);
        real axis2_z2 = grad_y(x,y,z*z,base1_z2,xsq1_z2,ysq1_z2,diag1_z2);
        real axis2_val = axis2_z0+axis2_z1+axis2_z2;
       // print ("\n axis2_val:",axis2_val);
       // print ("\n grad pos:", x,"_",y,"_",z);
        real axis3_z0 = grad_z(x,y,0,base1_z0,xsq1_z0,ysq1_z0,diag1_z0);
       // print ("\n axis3_z0:",axis3_z0);
        real axis3_z1 = grad_z(x,y,1,base1_z1,xsq1_z1,ysq1_z1,diag1_z1);
       // print ("\n axis3_z1:",axis3_z1);
        real axis3_z2 = grad_z(x,y,2*z,base1_z2,xsq1_z2,ysq1_z2,diag1_z2);
       // print ("\n axis3_z2:",axis3_z2);
        real axis3_val = axis3_z0+axis3_z1+axis3_z2;
       // print ("\n axis3_val:",axis3_val);
        ret =[axis1_val, axis2_val];
       // print ("\n grad pos:", x,"_",y,"_",z, " ret:",ret);
        }
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 vec2 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: 4\n");
        print("sizes: 2 ", sz0, " ", sz1, " ", sz2, "\n");
        print("kinds: 2-vector space space space\n");
        print("centers: none 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: none (", 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