Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

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

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

Wed Jul 20 03:02:00 2016 UTC (3 years, 6 months ago) by cchiw
File size: 5194 byte(s)
`added generic cases for trace|det and added test cases`
```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;

//-- field defined by coefficients --
// z=1
input vec2 base_z0 = [0, 0];    // b*x, c*y
input vec2 xsq_z0 = [0, 0];     // d*x^2, g*y*x^2
input vec2 ysq_z0 = [0, 0];     // f*y^2, h*x*y^2
input vec3 diag_z0 = [0, 0, 0];    // a, e*x*y,  i*y^2*x^2
// z=z
input vec2 base_z1 = [0, 0];    // b*x, c*y
input vec2 xsq_z1 = [0, 0];     // d*x^2, g*y*x^2
input vec2 ysq_z1 = [0, 0];     // f*y^2, h*x*y^2
input vec3 diag_z1 = [0, 0, 0];    // a, e*x*y,  i*y^2*x^2
// z=z*z
input vec2 base_z2 = [0, 0];    // b*x, c*y
input vec2 xsq_z2 = [0, 0];     // d*x^2, g*y*x^2
input vec2 ysq_z2 = [0, 0];     // f*y^2, h*x*y^2
input vec3 diag_z2 = [0, 0, 0];    // 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 func(vec3 pos) {
real x = pos[0];
real y = pos[1];
real z = pos[2];

real ret = 0;
if (0 == which) {
real z0 = cvt(x,y,1,base_z0,xsq_z0,ysq_z0,diag_z0);
real z1 = cvt(x,y,z,base_z1,xsq_z1,ysq_z1,diag_z1);
real z2 = cvt(x,y,z*z,base_z2,xsq_z2,ysq_z2,diag_z2);
ret = z0+z1+z2;

} 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
```