[Teem-commit] teem/src/ten tendHelix.c,1.14,1.15

SourceForge Headquarters 225 Broadway Suite 1600 San Diego, CA 92101 +1 (858) 454-5900

Update of /cvsroot/teem/teem/src/ten
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv23792

Modified Files:
	tendHelix.c 
Log Message:
now usefully using space directions, space origin, and measurement frame fields; used for Deft debugging

Index: tendHelix.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/tendHelix.c,v
retrieving revision 1.14
retrieving revision 1.15
diff -C2 -d -r1.14 -r1.15
*** tendHelix.c	11 Dec 2005 00:49:43 -0000	1.14
--- tendHelix.c	27 Dec 2005 16:58:27 -0000	1.15
***************
*** 27,69 ****
  char *_tend_helixInfoL =
    (INFO
!    ". The main utility of such a field is to expose the handedness of the "
!    "raster coordinate frame, and to see if it agrees with the coordinate "
!    "frame in which the tensor is expressed.  With the default settings, the "
!    "region of high anisotropy is a right-handed helix (same as DNA), "
!    "provided that the three axis sizes given to the \"-s\" option (below) "
!    "correspond to a right-handed ordered basis (\"slow equals fast cross "
!    "medium\").  In addition, the tensors twist relative to the helix. "
     "If you trace paths guided by the principal eigenvector of the tensors, "
     "along the surface of the helical cylinder, you get another "
     "right-handed helix, as if the the tensor field is modeling the result "
!    "if twisting a set of fibers into single-straded helical bundle.  "
!    "On the other hand, "
!    "if you sample in a straight light along the diameter of the cylinder, "
!    "the principal axis of the tensors twist in a *left*-handed helix. ");
  
  void
! tend_helixDoit(Nrrd *nout, float bnd,
!                float r, float R, float S, float angle, float ev[3]) {
!   int sx, sy, sz, xi, yi, zi, idx;
    double th, t0, t1, t2, t3, v1, v2,
!     wpos[3], vpos[3],
      W2H[9], H2W[9], H2C[9], C2H[9], fv[3], rv[3], uv[3], mA[9], mB[9], inside,
      tmp[3], len;
    float *out;
  
-   out = (float*)nout->data;
    sx = nout->axis[1].size;
    sy = nout->axis[2].size;
    sz = nout->axis[3].size;
!   idx = 0;
    for (zi=0; zi<sz; zi++) {
      fprintf(stderr, "zi = %d/%d\n", zi, sz);
-     vpos[2] = AIR_AFFINE(0, zi, sz-1, nout->axis[3].min, nout->axis[3].max);
      for (yi=0; yi<sy; yi++) {
!       vpos[1] = AIR_AFFINE(0, yi, sy-1, nout->axis[2].min, nout->axis[2].max);
!       for (xi=0; xi<sx; xi++, idx++) {
!         vpos[0] = AIR_AFFINE(0, xi, sx-1, nout->axis[1].min,nout->axis[1].max);
  
! #define WPOS(pos, th) ELL_3V_SET((pos), R*cos(th), R*sin(th), S*(th)/(2*AIR_PI))
  #define VAL(th) (WPOS(wpos, th), ELL_3V_DIST(wpos, vpos))
  #define RR 0.61803399
--- 27,72 ----
  char *_tend_helixInfoL =
    (INFO
!    ". The main utility of such a field is to debug handling of coordinate "
!    "systems in tensor field visualization.  The \"space directions\" and "
!    "\"space origin\" fields of the NRRD header determines the mapping from "
!    "coordinates in the index space of the image to coordinates in the "
!    "world space in which the image is "
!    "sampled.  The \"measurement frame\" field determines the mapping from "
!    "the coordinates of the tensor itself, to coordinates of the world space. "
!    "When these are correctly handled, the "
!    "region of high anisotropy is a right-handed helix (same as DNA). "
!    "Using differing axes sizes (via \"-s\") helps make sure that the "
!    "raster ordering of axes is correct.  In addition, the tensors twist "
!    "relative to the helix, which exposes handling of the measurement frame. "
     "If you trace paths guided by the principal eigenvector of the tensors, "
     "along the surface of the helical cylinder, you get another "
     "right-handed helix, as if the the tensor field is modeling the result "
!    "if twisting a set of fibers into single-stranded helical bundle. ");
  
  void
! tend_helixDoit(Nrrd *nout, double bnd,
!                double orig[3], double i2w[9], double mf[9],
!                double r, double R, double S, double angle, double ev[3]) {
!   int sx, sy, sz, xi, yi, zi;
    double th, t0, t1, t2, t3, v1, v2,
!     wpos[3], vpos[3], mfT[9],
      W2H[9], H2W[9], H2C[9], C2H[9], fv[3], rv[3], uv[3], mA[9], mB[9], inside,
      tmp[3], len;
    float *out;
  
    sx = nout->axis[1].size;
    sy = nout->axis[2].size;
    sz = nout->axis[3].size;
!   out = (float*)nout->data;
!   ELL_3M_TRANSPOSE(mfT, mf);
    for (zi=0; zi<sz; zi++) {
      fprintf(stderr, "zi = %d/%d\n", zi, sz);
      for (yi=0; yi<sy; yi++) {
!       for (xi=0; xi<sx; xi++) {
!         ELL_3V_SET(tmp, xi, yi, zi);
!         ELL_3MV_MUL(vpos, i2w, tmp);
!         ELL_3V_INCR(vpos, orig);
  
! #define WPOS(pos, th) ELL_3V_SET((pos),R*cos(th), R*sin(th), S*(th)/(2*AIR_PI))
  #define VAL(th) (WPOS(wpos, th), ELL_3V_DIST(wpos, vpos))
  #define RR 0.61803399
***************
*** 83,87 ****
          v1 = VAL(t1);
          v2 = VAL(t2);
!         while ( t3-t0 > 0.0001*(AIR_ABS(t1+t2)) ) {
            if (v1 < v2) {
              SHIFT3(t3, t2, t1, CC*t0 + RR*t2);
--- 86,90 ----
          v1 = VAL(t1);
          v2 = VAL(t2);
!         while ( t3-t0 > 0.000001*(AIR_ABS(t1)+AIR_ABS(t2)) ) {
            if (v1 < v2) {
              SHIFT3(t3, t2, t1, CC*t0 + RR*t2);
***************
*** 92,101 ****
            }
          }
! 
!         /* well-written code is self-documenting */
  
          WPOS(wpos, t1);
          ELL_3V_SUB(wpos, vpos, wpos);
!         ELL_3V_SET(fv, -R*sin(t1), R*cos(t1), S/AIR_PI);
          ELL_3V_NORM(fv, fv, len);
          ELL_3V_COPY(rv, wpos);
--- 95,104 ----
            }
          }
!         /* t1 (and t2) are now the th for which the point on the helix
!            (R*cos(th), R*sin(th), S*(th)/(2*AIR_PI)) is closest to vpos */
  
          WPOS(wpos, t1);
          ELL_3V_SUB(wpos, vpos, wpos);
!         ELL_3V_SET(fv, -R*sin(t1), R*cos(t1), S/AIR_PI);  /* helix tangent */
          ELL_3V_NORM(fv, fv, len);
          ELL_3V_COPY(rv, wpos);
***************
*** 104,111 ****
          ELL_3V_SCALE(tmp, -len, fv);
          ELL_3V_ADD2(rv, rv, tmp);
!         ELL_3V_NORM(rv, rv, len);
          ELL_3V_CROSS(uv, rv, fv);
!         ELL_3V_NORM(uv, uv, len);
!         ELL_3MV_ROW0_SET(W2H, uv);
          ELL_3MV_ROW1_SET(W2H, rv);
          ELL_3MV_ROW2_SET(W2H, fv);
--- 107,115 ----
          ELL_3V_SCALE(tmp, -len, fv);
          ELL_3V_ADD2(rv, rv, tmp);
!         ELL_3V_NORM(rv, rv, len);  /* rv now normal to helix, closest to 
!                                       pointing to vpos */
          ELL_3V_CROSS(uv, rv, fv);
!         ELL_3V_NORM(uv, uv, len);  /* (rv,fv,uv) now right-handed frame */
!         ELL_3MV_ROW0_SET(W2H, uv); /* as is (uv,rv,fv) */
          ELL_3MV_ROW1_SET(W2H, rv);
          ELL_3MV_ROW2_SET(W2H, fv);
***************
*** 121,129 ****
          ELL_3M_MUL(mB, mA, H2C);
          ELL_3M_MUL(mA, mB, W2H);
!         ELL_3M_MUL(mB, C2H, mA);
!         ELL_3M_MUL(mA, H2W, mB);
          
!         TEN_M2T_TT(out + 7*idx, float, mA);
!         (out + 7*idx)[0] = 1.0;
        }
      }
--- 125,136 ----
          ELL_3M_MUL(mB, mA, H2C);
          ELL_3M_MUL(mA, mB, W2H);
!         ELL_3M_MUL(mB, mA, mf);
!         ELL_3M_MUL(mA, C2H, mB);
!         ELL_3M_MUL(mB, H2W, mA);
!         ELL_3M_MUL(mA, mfT, mB);
          
!         TEN_M2T_TT(out, float, mA);
!         out[0] = 1.0;
!         out += 7;
        }
      }
***************
*** 141,146 ****
    int size[3];
    Nrrd *nout;
!   float R, r, S, bnd, angle, ev[3];
!   double min[4], max[4];
    char *outS;
  
--- 148,153 ----
    int size[3];
    Nrrd *nout;
!   double R, r, S, bnd, angle, ev[3], ip[3], iq[4], mp[3], mq[4], tmp[9],
!     orig[3], i2w[9], rot[9], mf[9], spd[4][3];
    char *outS;
  
***************
*** 150,168 ****
               "slightly different sizes here, to expose errors in interpreting "
               "axis ordering (e.g. \"-s 39 40 41\")");
!   hestOptAdd(&hopt, "min", "min corner", airTypeDouble, 3, 3, min+1,
!              "-2 -2 -2",
!              "location of low corner of sampled tensor volume");
!   hestOptAdd(&hopt, "max", "max corner", airTypeDouble, 3, 3, max+1, "2 2 2",
!              "location of high corner of sampled tensor volume");
!   hestOptAdd(&hopt, "b", "boundary", airTypeFloat, 1, 1, &bnd, "0.05",
               "parameter governing how fuzzy the boundary between high and "
               "low anisotropy is. Use \"-b 0\" for no fuzziness");
!   hestOptAdd(&hopt, "r", "little radius", airTypeFloat, 1, 1, &r, "0.5",
               "(minor) radius of cylinder tracing helix");
!   hestOptAdd(&hopt, "R", "big radius", airTypeFloat, 1, 1, &R, "1.2",
               "(major) radius of helical turns");
!   hestOptAdd(&hopt, "S", "spacing", airTypeFloat, 1, 1, &S, "2",
               "spacing between turns of helix (along its axis)");
!   hestOptAdd(&hopt, "a", "angle", airTypeFloat, 1, 1, &angle, "1.0",
               "maximal angle of twist of tensors along path.  There is no "
               "twist at helical core of path, and twist increases linearly "
--- 157,176 ----
               "slightly different sizes here, to expose errors in interpreting "
               "axis ordering (e.g. \"-s 39 40 41\")");
!   hestOptAdd(&hopt, "ip", "image orientation", airTypeDouble, 3, 3, ip,
!              "0 0 0",
!              "quaternion quotient space orientation of image");
!   hestOptAdd(&hopt, "mp", "measurement orientation", airTypeDouble, 3, 3, mp,
!              "0 0 0",
!              "quaternion quotient space orientation of measurement frame");
!   hestOptAdd(&hopt, "b", "boundary", airTypeDouble, 1, 1, &bnd, "5",
               "parameter governing how fuzzy the boundary between high and "
               "low anisotropy is. Use \"-b 0\" for no fuzziness");
!   hestOptAdd(&hopt, "r", "little radius", airTypeDouble, 1, 1, &r, "30",
               "(minor) radius of cylinder tracing helix");
!   hestOptAdd(&hopt, "R", "big radius", airTypeDouble, 1, 1, &R, "50",
               "(major) radius of helical turns");
!   hestOptAdd(&hopt, "S", "spacing", airTypeDouble, 1, 1, &S, "100",
               "spacing between turns of helix (along its axis)");
!   hestOptAdd(&hopt, "a", "angle", airTypeDouble, 1, 1, &angle, "1.0",
               "maximal angle of twist of tensors along path.  There is no "
               "twist at helical core of path, and twist increases linearly "
***************
*** 170,174 ****
               "positive spacing resulting in a right-handed twist around a "
               "right-handed helix. ");
!   hestOptAdd(&hopt, "ev", "eigenvalues", airTypeFloat, 3, 3, ev, "0.9 0.4 0.2",
               "eigenvalues of tensors (in order) along direction of coil, "
               "circumferential around coil, and radial around coil. ");
--- 178,183 ----
               "positive spacing resulting in a right-handed twist around a "
               "right-handed helix. ");
!   hestOptAdd(&hopt, "ev", "eigenvalues", airTypeDouble, 3, 3, ev,
!              "0.006 0.002 0.001",
               "eigenvalues of tensors (in order) along direction of coil, "
               "circumferential around coil, and radial around coil. ");
***************
*** 193,204 ****
      airMopError(mop); return 1;
    }
-   min[0] = max[0] = AIR_NAN;
-   nrrdAxisInfoSet_nva(nout, nrrdAxisInfoMin, min);
-   nrrdAxisInfoSet_nva(nout, nrrdAxisInfoMax, max);
  
!   tend_helixDoit(nout, bnd, r, R, S, angle, ev);
!   nrrdAxisInfoSpacingSet(nout, 1);
!   nrrdAxisInfoSpacingSet(nout, 2);
!   nrrdAxisInfoSpacingSet(nout, 3);
  
    if (nrrdSave(outS, nout, NULL)) {
--- 202,245 ----
      airMopError(mop); return 1;
    }
  
!   iq[0] = 1.0;
!   ELL_3V_COPY(iq + 1, ip);
!   ell_q_to_3m_d(rot, iq);
!   ELL_3V_SET(orig,
!              -2*R + 2*R/size[0],
!              -2*R + 2*R/size[1],
!              -2*R + 2*R/size[2]);
!   ELL_3M_ZERO_SET(i2w);
!   ELL_3M_DIAG_SET(i2w, 4*R/size[0], 4*R/size[1], 4*R/size[2]);
!   ELL_3MV_MUL(tmp, rot, orig);
!   ELL_3V_COPY(orig, tmp);
!   ELL_3M_MUL(tmp, rot, i2w);
!   ELL_3M_COPY(i2w, tmp);
!   mq[0] = 1.0;
!   ELL_3V_COPY(mq + 1, mp);
!   ell_q_to_3m_d(mf, mq);
!   tend_helixDoit(nout, bnd,
!                  orig, i2w, mf,
!                  r, R, S, angle, ev);
!   nrrdSpaceSet(nout, nrrdSpaceRightAnteriorSuperior);
!   nrrdSpaceOriginSet(nout, orig);
!   ELL_3V_SET(spd[0], AIR_NAN, AIR_NAN, AIR_NAN);
!   ELL_3MV_COL0_GET(spd[1], i2w);
!   ELL_3MV_COL1_GET(spd[2], i2w);
!   ELL_3MV_COL2_GET(spd[3], i2w);
!   nrrdAxisInfoSet_va(nout, nrrdAxisInfoSpaceDirection,
!                      spd[0], spd[1], spd[2], spd[3]);
!   nrrdAxisInfoSet_va(nout, nrrdAxisInfoCenter,
!                      nrrdCenterUnknown, nrrdCenterCell,
!                      nrrdCenterCell, nrrdCenterCell);
!   nout->measurementFrame[0][0] = mf[0];
!   nout->measurementFrame[1][0] = mf[1];
!   nout->measurementFrame[2][0] = mf[2];
!   nout->measurementFrame[0][1] = mf[3];
!   nout->measurementFrame[1][1] = mf[4];
!   nout->measurementFrame[2][1] = mf[5];
!   nout->measurementFrame[0][2] = mf[6];
!   nout->measurementFrame[1][2] = mf[7];
!   nout->measurementFrame[2][2] = mf[8];
  
    if (nrrdSave(outS, nout, NULL)) {

[Teem-commit] teem/src/ten tendHelix.c,1.14,1.15

Tools to process and visualize scientific data and images

[Teem-commit] teem/src/ten tendHelix.c,1.14,1.15