teem-commit — gets email notification of all Teem svn commits

[Teem-commit] teem/src/ell eigen.c,1.17,1.18 ell.h,1.27,1.28 ellMacros.h,1.39,1.40 genmat.c,1.4,1.5 miscEll.c,1.5,1.6 quat.c,1.3,1.4

[<kin...@us...>]

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

Modified Files:
	eigen.c ell.h ellMacros.h genmat.c miscEll.c quat.c 
Log Message:
ell's storage of matrices is now row-major instead of column major.  Debugging in progress

Index: eigen.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ell/eigen.c,v
retrieving revision 1.17
retrieving revision 1.18
diff -C2 -d -r1.17 -r1.18
*** eigen.c	19 Feb 2004 11:10:34 -0000	1.17
--- eigen.c	1 Mar 2004 10:19:44 -0000	1.18
***************
*** 103,108 ****
  **
  ** the given matrix is assumed to have a nullspace of dimension one.
! ** The COLUMNS of the matrix are n+0, n+3, n+6.  A normalized vector
! ** which spans the nullspace is put into ans.
  **
  ** The given nullspace matrix is NOT modified.  
--- 103,107 ----
  **
  ** the given matrix is assumed to have a nullspace of dimension one.
! ** A normalized vector which spans the nullspace is put into ans.
  **
  ** The given nullspace matrix is NOT modified.  
***************
*** 111,117 ****
  */
  void
! ell_3m_1d_nullspace_d(double ans[3], double n[9]) {
!   double t[9], norm;
    
    /* find the three cross-products of pairs of column vectors of n */
    ELL_3V_CROSS(t+0, n+0, n+3);
--- 110,117 ----
  */
  void
! ell_3m_1d_nullspace_d(double ans[3], double _n[9]) {
!   double t[9], n[9], norm;
    
+   ELL_3M_TRANSPOSE(n, _n);
    /* find the three cross-products of pairs of column vectors of n */
    ELL_3V_CROSS(t+0, n+0, n+3);
***************
*** 134,138 ****
  **
  ** the given matrix is assumed to have a nullspace of dimension two.
- ** The COLUMNS of the matrix are n+0, n+3, n+6
  **
  ** The given nullspace matrix is NOT modified 
--- 134,137 ----
***************
*** 144,148 ****
    double n[9], tmp[3], norm;
  
!   ELL_3M_COPY(n, _n);
    _ell_align3_d(n);
    ELL_3V_ADD3(tmp, n+0, n+3, n+6);
--- 143,147 ----
    double n[9], tmp[3], norm;
  
!   ELL_3M_TRANSPOSE(n, _n);
    _ell_align3_d(n);
    ELL_3V_ADD3(tmp, n+0, n+3, n+6);
***************
*** 163,168 ****
  ** finds eigenvalues of given matrix.
  **
- ** m+0, m+3, m+6, are the COLUMNS of the matrix
- **
  ** returns information about the roots according to ellCubeRoot enum,
  ** see header for ellCubic for details.
--- 162,165 ----
***************
*** 190,199 ****
    */
    A = -m[0] - m[4] - m[8];
!   B = m[0]*m[4] - m[1]*m[3] 
!     + m[0]*m[8] - m[2]*m[6] 
!     + m[4]*m[8] - m[5]*m[7];
!   C = (m[2]*m[4] - m[1]*m[5])*m[6]
!     + (m[0]*m[5] - m[2]*m[3])*m[7]
!     + (m[1]*m[3] - m[0]*m[4])*m[8];
    roots = ell_cubic(eval, A, B, C, newton);
    if (ell_cubic_root_three == roots 
--- 187,196 ----
    */
    A = -m[0] - m[4] - m[8];
!   B = m[0]*m[4] - m[3]*m[1] 
!     + m[0]*m[8] - m[6]*m[2] 
!     + m[4]*m[8] - m[7]*m[5];
!   C = (m[6]*m[4] - m[3]*m[7])*m[2]
!     + (m[0]*m[7] - m[6]*m[1])*m[5]
!     + (m[3]*m[1] - m[0]*m[4])*m[8];
    roots = ell_cubic(eval, A, B, C, newton);
    if (ell_cubic_root_three == roots 
***************
*** 208,213 ****
  ******** ell_3m_eigensolve_d()
  **
! ** finds eigenvalues and eigenvectors of given matrix m,
! ** m+0, m+3, m+6, are the COLUMNS of the matrix.
  **
  ** returns information about the roots according to ellCubeRoot enum,
--- 205,209 ----
  ******** ell_3m_eigensolve_d()
  **
! ** finds eigenvalues and eigenvectors of given matrix m
  **
  ** returns information about the roots according to ellCubeRoot enum,
***************
*** 228,234 ****
    /* if (ell_debug) {
      printf("ell_3m_eigensolve: input matrix:\n");
!     printf("{{%20.15f,\t%20.15f,\t%20.15f},\n", m[0], m[3], m[6]);
!     printf(" {%20.15f,\t%20.15f,\t%20.15f},\n", m[1], m[4], m[7]);
!     printf(" {%20.15f,\t%20.15f,\t%20.15f}};\n",m[2], m[5], m[8]);
      } */
    
--- 224,230 ----
    /* if (ell_debug) {
      printf("ell_3m_eigensolve: input matrix:\n");
!     printf("{{%20.15f,\t%20.15f,\t%20.15f},\n", m[0], m[1], m[2]);
!     printf(" {%20.15f,\t%20.15f,\t%20.15f},\n", m[3], m[4], m[5]);
!     printf(" {%20.15f,\t%20.15f,\t%20.15f}};\n",m[6], m[7], m[8]);
      } */
    

Index: ell.h
===================================================================
RCS file: /cvsroot/teem/teem/src/ell/ell.h,v
retrieving revision 1.27
retrieving revision 1.28
diff -C2 -d -r1.27 -r1.28
*** ell.h	19 Feb 2004 03:52:47 -0000	1.27
--- ell.h	1 Mar 2004 10:19:44 -0000	1.28
***************
*** 51,55 ****
  
  /* 
! ** Note: matrix element ordering is:
  **
  **     0   3   6
--- 51,55 ----
  
  /* 
! ** Note: *** PRE-TEEM 1.7 *** matrix element ordering was:
  **
  **     0   3   6
***************
*** 62,66 ****
  **   3   7  11  15
  **
! ** and all vectors are COLUMN vectors
  */
  
--- 62,77 ----
  **   3   7  11  15
  **
! ** as of TEEM 1.7, matrix element ordering is:
! **
! **     0   1   2
! **     3   4   5      for 3x3
! **     6   7   8
! **
! **   0   1   2   3
! **   4   5   6   7    for 4x4
! **   8   9  10  11
! **  12  13  14  15
! **
! ** all vectors are still, logically, COLUMN vectors
  */
  
***************
*** 164,168 ****
  
  /* genmat.c */
! TEEM_API int ell_Nm_valid(Nrrd *mat);
  TEEM_API int ell_Nm_tran(Nrrd *dest, Nrrd *src);
  TEEM_API int ell_Nm_mul(Nrrd *dest, Nrrd *A, Nrrd *B);
--- 175,179 ----
  
  /* genmat.c */
! TEEM_API int ell_Nm_check(Nrrd *mat);
  TEEM_API int ell_Nm_tran(Nrrd *dest, Nrrd *src);
  TEEM_API int ell_Nm_mul(Nrrd *dest, Nrrd *A, Nrrd *B);
***************
*** 171,181 ****
  
  /* cubic.c */
! TEEM_API int ell_cubic(double root[3], double A, double B, double C, int newton);
  
  /* eigen.c */
  TEEM_API void ell_3m_1d_nullspace_d(double ans[3], double n[9]);
! TEEM_API void ell_3m_2d_nullspace_d(double ans0[3], double ans1[3], double n[9]);
! TEEM_API int ell_3m_eigenvalues_d(double eval[3], double m[9], 
! 				  int newton);
  TEEM_API int ell_3m_eigensolve_d(double eval[3], double evec[9], double m[9],
  				 int newton);
--- 182,194 ----
  
  /* cubic.c */
! TEEM_API int ell_cubic(double root[3],
! 		       double A, double B, double C, int newton);
  
  /* eigen.c */
  TEEM_API void ell_3m_1d_nullspace_d(double ans[3], double n[9]);
! TEEM_API void ell_3m_2d_nullspace_d(double ans0[3],
! 				    double ans1[3], double n[9]);
! TEEM_API int ell_3m_eigenvalues_d(double eval[3], 
! 				  double m[9], int newton);
  TEEM_API int ell_3m_eigensolve_d(double eval[3], double evec[9], double m[9],
  				 int newton);

Index: ellMacros.h
===================================================================
RCS file: /cvsroot/teem/teem/src/ell/ellMacros.h,v
retrieving revision 1.39
retrieving revision 1.40
diff -C2 -d -r1.39 -r1.40
*** ellMacros.h	23 Feb 2004 02:54:22 -0000	1.39
--- ellMacros.h	1 Mar 2004 10:19:44 -0000	1.40
***************
*** 87,91 ****
        : 1))
  
! #define _ELL_2M_DET(a,b,c,d) ((a)*(d) - (c)*(b))
  
  #define ELL_2M_DET(m) _ELL_2M_DET((m)[0],(m)[1],(m)[2],(m)[3])
--- 87,91 ----
        : 1))
  
! #define _ELL_2M_DET(a,b,c,d) ((a)*(d) - (b)*(c))
  
  #define ELL_2M_DET(m) _ELL_2M_DET((m)[0],(m)[1],(m)[2],(m)[3])
***************
*** 93,99 ****
  /*
  ** the 3x3 matrix-related macros assume that the matrix indexing is:
! ** 0  3  6
! ** 1  4  7
! ** 2  5  8
  */
  
--- 93,99 ----
  /*
  ** the 3x3 matrix-related macros assume that the matrix indexing is:
! ** 0  1  2
! ** 3  4  5
! ** 6  7  8
  */
  
***************
*** 298,368 ****
  
  #define ELL_3MV_COL0_GET(v, m) \
!   (ELL_3V_SET((v), (m)[0], (m)[1], (m)[2]))
  
  #define ELL_3MV_COL1_GET(v, m) \
!   (ELL_3V_SET((v), (m)[3], (m)[4], (m)[5]))
  
  #define ELL_3MV_COL2_GET(v, m) \
!   (ELL_3V_SET((v), (m)[6], (m)[7], (m)[8]))
  
  #define ELL_3MV_ROW0_GET(v, m) \
!   (ELL_3V_SET((v), (m)[0], (m)[3], (m)[6]))
  
  #define ELL_3MV_ROW1_GET(v, m) \
!   (ELL_3V_SET((v), (m)[1], (m)[4], (m)[7]))
  
  #define ELL_3MV_ROW2_GET(v, m) \
!   (ELL_3V_SET((v), (m)[2], (m)[5], (m)[8]))
  
  #define ELL_3MV_COL0_SET(m, v) \
!   (ELL_3V_GET((m)[0], (m)[1], (m)[2], (v)))
  
  #define ELL_3MV_COL1_SET(m, v) \
!   (ELL_3V_GET((m)[3], (m)[4], (m)[5], (v)))
  
  #define ELL_3MV_COL2_SET(m, v) \
!   (ELL_3V_GET((m)[6], (m)[7], (m)[8], (v)))
  
  #define ELL_3MV_ROW0_SET(m, v) \
!   (ELL_3V_GET((m)[0], (m)[3], (m)[6], (v)))
  
  #define ELL_3MV_ROW1_SET(m, v) \
!   (ELL_3V_GET((m)[1], (m)[4], (m)[7], (v)))
  
  #define ELL_3MV_ROW2_SET(m, v) \
!   (ELL_3V_GET((m)[2], (m)[5], (m)[8], (v)))
  
  #define ELL_3MV_OUTER(m, v1, v2) \
!   (ELL_3V_SCALE((m)+0, (v2)[0], (v1)), \
!    ELL_3V_SCALE((m)+3, (v2)[1], (v1)), \
!    ELL_3V_SCALE((m)+6, (v2)[2], (v1)))
  
  #define ELL_3MV_OUTER_ADD(m, v1, v2) \
!   (ELL_3V_SCALE_ADD2((m)+0, 1, (m)+0, (v2)[0], (v1)), \
!    ELL_3V_SCALE_ADD2((m)+3, 1, (m)+3, (v2)[1], (v1)), \
!    ELL_3V_SCALE_ADD2((m)+6, 1, (m)+6, (v2)[2], (v1)))
  
  #define ELL_3MV_MUL(v2, m, v1) \
-   ((v2)[0] = (m)[0]*(v1)[0] + (m)[3]*(v1)[1] + (m)[6]*(v1)[2], \
-    (v2)[1] = (m)[1]*(v1)[0] + (m)[4]*(v1)[1] + (m)[7]*(v1)[2], \
-    (v2)[2] = (m)[2]*(v1)[0] + (m)[5]*(v1)[1] + (m)[8]*(v1)[2])
- 
- #define ELL_3MV_TMUL(v2, m, v1) \
    ((v2)[0] = (m)[0]*(v1)[0] + (m)[1]*(v1)[1] + (m)[2]*(v1)[2], \
     (v2)[1] = (m)[3]*(v1)[0] + (m)[4]*(v1)[1] + (m)[5]*(v1)[2], \
     (v2)[2] = (m)[6]*(v1)[0] + (m)[7]*(v1)[1] + (m)[8]*(v1)[2])
  
  #define ELL_3M_MUL(m3, m1, m2)                                    \
!   ((m3)[0] = (m1)[0]*(m2)[0] + (m1)[3]*(m2)[1] + (m1)[6]*(m2)[2], \
!    (m3)[1] = (m1)[1]*(m2)[0] + (m1)[4]*(m2)[1] + (m1)[7]*(m2)[2], \
!    (m3)[2] = (m1)[2]*(m2)[0] + (m1)[5]*(m2)[1] + (m1)[8]*(m2)[2], \
                                                                    \
!    (m3)[3] = (m1)[0]*(m2)[3] + (m1)[3]*(m2)[4] + (m1)[6]*(m2)[5], \
!    (m3)[4] = (m1)[1]*(m2)[3] + (m1)[4]*(m2)[4] + (m1)[7]*(m2)[5], \
!    (m3)[5] = (m1)[2]*(m2)[3] + (m1)[5]*(m2)[4] + (m1)[8]*(m2)[5], \
                                                                    \
!    (m3)[6] = (m1)[0]*(m2)[6] + (m1)[3]*(m2)[7] + (m1)[6]*(m2)[8], \
!    (m3)[7] = (m1)[1]*(m2)[6] + (m1)[4]*(m2)[7] + (m1)[7]*(m2)[8], \
!    (m3)[8] = (m1)[2]*(m2)[6] + (m1)[5]*(m2)[7] + (m1)[8]*(m2)[8])
  
  #define ELL_3M_SCALE_SET(m, x, y, z)  \
--- 298,368 ----
  
  #define ELL_3MV_COL0_GET(v, m) \
!   (ELL_3V_SET((v), (m)[0], (m)[3], (m)[6]))
  
  #define ELL_3MV_COL1_GET(v, m) \
!   (ELL_3V_SET((v), (m)[1], (m)[4], (m)[7]))
  
  #define ELL_3MV_COL2_GET(v, m) \
!   (ELL_3V_SET((v), (m)[2], (m)[5], (m)[8]))
  
  #define ELL_3MV_ROW0_GET(v, m) \
!   (ELL_3V_SET((v), (m)[0], (m)[1], (m)[2]))
  
  #define ELL_3MV_ROW1_GET(v, m) \
!   (ELL_3V_SET((v), (m)[3], (m)[4], (m)[5]))
  
  #define ELL_3MV_ROW2_GET(v, m) \
!   (ELL_3V_SET((v), (m)[6], (m)[7], (m)[8]))
  
  #define ELL_3MV_COL0_SET(m, v) \
!   (ELL_3V_GET((m)[0], (m)[3], (m)[6], (v)))
  
  #define ELL_3MV_COL1_SET(m, v) \
!   (ELL_3V_GET((m)[1], (m)[4], (m)[7], (v)))
  
  #define ELL_3MV_COL2_SET(m, v) \
!   (ELL_3V_GET((m)[2], (m)[5], (m)[8], (v)))
  
  #define ELL_3MV_ROW0_SET(m, v) \
!   (ELL_3V_GET((m)[0], (m)[1], (m)[2], (v)))
  
  #define ELL_3MV_ROW1_SET(m, v) \
!   (ELL_3V_GET((m)[3], (m)[4], (m)[5], (v)))
  
  #define ELL_3MV_ROW2_SET(m, v) \
!   (ELL_3V_GET((m)[6], (m)[7], (m)[8], (v)))
  
  #define ELL_3MV_OUTER(m, v1, v2) \
!   (ELL_3V_SCALE((m)+0, (v1)[0], (v2)), \
!    ELL_3V_SCALE((m)+3, (v1)[1], (v2)), \
!    ELL_3V_SCALE((m)+6, (v1)[2], (v2)))
  
  #define ELL_3MV_OUTER_ADD(m, v1, v2) \
!   (ELL_3V_SCALE_INCR((m)+0, (v1)[0], (v2)), \
!    ELL_3V_SCALE_INCR((m)+3, (v1)[1], (v2)), \
!    ELL_3V_SCALE_INCR((m)+6, (v1)[2], (v2)))
  
  #define ELL_3MV_MUL(v2, m, v1) \
    ((v2)[0] = (m)[0]*(v1)[0] + (m)[1]*(v1)[1] + (m)[2]*(v1)[2], \
     (v2)[1] = (m)[3]*(v1)[0] + (m)[4]*(v1)[1] + (m)[5]*(v1)[2], \
     (v2)[2] = (m)[6]*(v1)[0] + (m)[7]*(v1)[1] + (m)[8]*(v1)[2])
  
+ #define ELL_3MV_TMUL(v2, m, v1) \
+   ((v2)[0] = (m)[0]*(v1)[0] + (m)[3]*(v1)[1] + (m)[6]*(v1)[2], \
+    (v2)[1] = (m)[1]*(v1)[0] + (m)[4]*(v1)[1] + (m)[7]*(v1)[2], \
+    (v2)[2] = (m)[2]*(v1)[0] + (m)[5]*(v1)[1] + (m)[8]*(v1)[2])
+ 
  #define ELL_3M_MUL(m3, m1, m2)                                    \
!   ((m3)[0] = (m1)[0]*(m2)[0] + (m1)[1]*(m2)[3] + (m1)[2]*(m2)[6], \
!    (m3)[1] = (m1)[0]*(m2)[1] + (m1)[1]*(m2)[4] + (m1)[2]*(m2)[7], \
!    (m3)[2] = (m1)[0]*(m2)[2] + (m1)[1]*(m2)[5] + (m1)[2]*(m2)[8], \
                                                                    \
!    (m3)[3] = (m1)[3]*(m2)[0] + (m1)[4]*(m2)[3] + (m1)[5]*(m2)[6], \
!    (m3)[4] = (m1)[3]*(m2)[1] + (m1)[4]*(m2)[4] + (m1)[5]*(m2)[7], \
!    (m3)[5] = (m1)[3]*(m2)[2] + (m1)[4]*(m2)[5] + (m1)[5]*(m2)[8], \
                                                                    \
!    (m3)[6] = (m1)[6]*(m2)[0] + (m1)[7]*(m2)[3] + (m1)[8]*(m2)[6], \
!    (m3)[7] = (m1)[6]*(m2)[1] + (m1)[7]*(m2)[4] + (m1)[8]*(m2)[7], \
!    (m3)[8] = (m1)[6]*(m2)[2] + (m1)[7]*(m2)[5] + (m1)[8]*(m2)[8])
  
  #define ELL_3M_SCALE_SET(m, x, y, z)  \
***************
*** 373,387 ****
  #define ELL_3M_ROTATE_X_SET(m, th)               \
    (ELL_3V_SET((m)+ 0,  1 ,     0    ,     0   ), \
!    ELL_3V_SET((m)+ 3,  0 ,  cos(th) , +sin(th)), \
!    ELL_3V_SET((m)+ 6,  0 , -sin(th) ,  cos(th)))
  
  #define ELL_3M_ROTATE_Y_SET(m, th)               \
!   (ELL_3V_SET((m)+ 0,  cos(th) ,  0 , -sin(th)), \
     ELL_3V_SET((m)+ 3,     0    ,  1 ,     0   ), \
!    ELL_3V_SET((m)+ 6, +sin(th) ,  0 ,  cos(th)))
  
  #define ELL_3M_ROTATE_Z_SET(m, th)               \
!   (ELL_3V_SET((m)+ 0,  cos(th) , +sin(th) ,  0), \
!    ELL_3V_SET((m)+ 3, -sin(th) ,  cos(th) ,  0), \
     ELL_3V_SET((m)+ 6,     0    ,     0    ,  1))
  
--- 373,387 ----
  #define ELL_3M_ROTATE_X_SET(m, th)               \
    (ELL_3V_SET((m)+ 0,  1 ,     0    ,     0   ), \
!    ELL_3V_SET((m)+ 3,  0 ,  cos(th) , -sin(th)), \
!    ELL_3V_SET((m)+ 6,  0 , +sin(th) ,  cos(th)))
  
  #define ELL_3M_ROTATE_Y_SET(m, th)               \
!   (ELL_3V_SET((m)+ 0,  cos(th) ,  0 , +sin(th)), \
     ELL_3V_SET((m)+ 3,     0    ,  1 ,     0   ), \
!    ELL_3V_SET((m)+ 6, -sin(th) ,  0 ,  cos(th)))
  
  #define ELL_3M_ROTATE_Z_SET(m, th)               \
!   (ELL_3V_SET((m)+ 0,  cos(th) , -sin(th) ,  0), \
!    ELL_3V_SET((m)+ 3, +sin(th) ,  cos(th) ,  0), \
     ELL_3V_SET((m)+ 6,     0    ,     0    ,  1))
  
***************
*** 389,396 ****
  ** the 4x4 matrix-related macros assume that the matrix indexing is:
  **
! ** 0   4   8  12
! ** 1   5   9  13
! ** 2   6  10  14
! ** 3   7  11  15
  */
  
--- 389,396 ----
  ** the 4x4 matrix-related macros assume that the matrix indexing is:
  **
! **  0   1   2   3
! **  4   5   6   7
! **  8   9  10  11
! ** 12  13  14  15
  */
  
***************
*** 505,563 ****
     ELL_SWAP2((m)[11],(m)[14],(t)))
  
! #define ELL_4MV_COL0_GET(v, m) \
    (ELL_4V_SET((v), (m)[ 0], (m)[ 1], (m)[ 2], (m)[ 3]))
  
! #define ELL_4MV_COL1_GET(v, m) \
    (ELL_4V_SET((v), (m)[ 4], (m)[ 5], (m)[ 6], (m)[ 7]))
  
! #define ELL_4MV_COL2_GET(v, m) \
!   (ELL_4V_SET((v), (m)[ 8], (m)[ 9], (m)[10], (m)[11])
  
! #define ELL_4MV_COL3_GET(v, m) \
    (ELL_4V_SET((v), (m)[12], (m)[13], (m)[14], (m)[15]))
  
! #define ELL_4MV_ROW0_GET(v, m) \
    (ELL_4V_SET((v), (m)[ 0], (m)[ 4], (m)[ 8], (m)[12]))
  
! #define ELL_4MV_ROW1_GET(v, m) \
    (ELL_4V_SET((v), (m)[ 1], (m)[ 5], (m)[ 9], (m)[13]))
  
! #define ELL_4MV_ROW2_GET(v, m) \
    (ELL_4V_SET((v), (m)[ 2], (m)[ 6], (m)[10], (m)[14]))
  
! #define ELL_4MV_ROW3_GET(v, m) \
    (ELL_4V_SET((v), (m)[ 3], (m)[ 7], (m)[11], (m)[15]))
  
! #define ELL_4MV_COL0_SET(m, v) \
    (ELL_4V_GET((m)[ 0], (m)[ 1], (m)[ 2], (m)[ 3], (v)))
  
! #define ELL_4MV_COL1_SET(m, v) \
    (ELL_4V_GET((m)[ 4], (m)[ 5], (m)[ 6], (m)[ 7], (v)))
  
! #define ELL_4MV_COL2_SET(m, v) \
    (ELL_4V_GET((m)[ 8], (m)[ 9], (m)[10], (m)[11], (v)))
  
! #define ELL_4MV_COL3_SET(m, v) \
    (ELL_4V_GET((m)[12], (m)[13], (m)[14], (m)[15], (v)))
  
! #define ELL_4MV_ROW0_SET(m, v) \
    (ELL_4V_GET((m)[ 0], (m)[ 4], (m)[ 8], (m)[12], (v)))
  
! #define ELL_4MV_ROW1_SET(m, v) \
    (ELL_4V_GET((m)[ 1], (m)[ 5], (m)[ 9], (m)[13], (v)))
  
! #define ELL_4MV_ROW2_SET(m, v) \
    (ELL_4V_GET((m)[ 2], (m)[ 6], (m)[10], (m)[14], (v)))
  
! #define ELL_4MV_ROW3_SET(m, v) \
    (ELL_4V_GET((m)[ 3], (m)[ 7], (m)[11], (m)[15], (v)))
  
  #define ELL_4MV_MUL(v2, m, v1)                                              \
-   ((v2)[0]=(m)[ 0]*(v1)[0]+(m)[ 4]*(v1)[1]+(m)[ 8]*(v1)[2]+(m)[12]*(v1)[3], \
-    (v2)[1]=(m)[ 1]*(v1)[0]+(m)[ 5]*(v1)[1]+(m)[ 9]*(v1)[2]+(m)[13]*(v1)[3], \
-    (v2)[2]=(m)[ 2]*(v1)[0]+(m)[ 6]*(v1)[1]+(m)[10]*(v1)[2]+(m)[14]*(v1)[3], \
-    (v2)[3]=(m)[ 3]*(v1)[0]+(m)[ 7]*(v1)[1]+(m)[11]*(v1)[2]+(m)[15]*(v1)[3])
- 
- #define ELL_4MV_TMUL(v2, m, v1)                                             \
    ((v2)[0]=(m)[ 0]*(v1)[0]+(m)[ 1]*(v1)[1]+(m)[ 2]*(v1)[2]+(m)[ 3]*(v1)[3], \
     (v2)[1]=(m)[ 4]*(v1)[0]+(m)[ 5]*(v1)[1]+(m)[ 6]*(v1)[2]+(m)[ 7]*(v1)[3], \
--- 505,557 ----
     ELL_SWAP2((m)[11],(m)[14],(t)))
  
! #define ELL_4MV_ROW0_GET(v, m) \
    (ELL_4V_SET((v), (m)[ 0], (m)[ 1], (m)[ 2], (m)[ 3]))
  
! #define ELL_4MV_ROW1_GET(v, m) \
    (ELL_4V_SET((v), (m)[ 4], (m)[ 5], (m)[ 6], (m)[ 7]))
  
! #define ELL_4MV_ROW2_GET(v, m) \
!   (ELL_4V_SET((v), (m)[ 8], (m)[ 9], (m)[10], (m)[11]))
  
! #define ELL_4MV_ROW3_GET(v, m) \
    (ELL_4V_SET((v), (m)[12], (m)[13], (m)[14], (m)[15]))
  
! #define ELL_4MV_COL0_GET(v, m) \
    (ELL_4V_SET((v), (m)[ 0], (m)[ 4], (m)[ 8], (m)[12]))
  
! #define ELL_4MV_COL1_GET(v, m) \
    (ELL_4V_SET((v), (m)[ 1], (m)[ 5], (m)[ 9], (m)[13]))
  
! #define ELL_4MV_COL2_GET(v, m) \
    (ELL_4V_SET((v), (m)[ 2], (m)[ 6], (m)[10], (m)[14]))
  
! #define ELL_4MV_COL3_GET(v, m) \
    (ELL_4V_SET((v), (m)[ 3], (m)[ 7], (m)[11], (m)[15]))
  
! #define ELL_4MV_ROW0_SET(m, v) \
    (ELL_4V_GET((m)[ 0], (m)[ 1], (m)[ 2], (m)[ 3], (v)))
  
! #define ELL_4MV_ROW1_SET(m, v) \
    (ELL_4V_GET((m)[ 4], (m)[ 5], (m)[ 6], (m)[ 7], (v)))
  
! #define ELL_4MV_ROW2_SET(m, v) \
    (ELL_4V_GET((m)[ 8], (m)[ 9], (m)[10], (m)[11], (v)))
  
! #define ELL_4MV_ROW3_SET(m, v) \
    (ELL_4V_GET((m)[12], (m)[13], (m)[14], (m)[15], (v)))
  
! #define ELL_4MV_COL0_SET(m, v) \
    (ELL_4V_GET((m)[ 0], (m)[ 4], (m)[ 8], (m)[12], (v)))
  
! #define ELL_4MV_COL1_SET(m, v) \
    (ELL_4V_GET((m)[ 1], (m)[ 5], (m)[ 9], (m)[13], (v)))
  
! #define ELL_4MV_COL2_SET(m, v) \
    (ELL_4V_GET((m)[ 2], (m)[ 6], (m)[10], (m)[14], (v)))
  
! #define ELL_4MV_COL3_SET(m, v) \
    (ELL_4V_GET((m)[ 3], (m)[ 7], (m)[11], (m)[15], (v)))
  
  #define ELL_4MV_MUL(v2, m, v1)                                              \
    ((v2)[0]=(m)[ 0]*(v1)[0]+(m)[ 1]*(v1)[1]+(m)[ 2]*(v1)[2]+(m)[ 3]*(v1)[3], \
     (v2)[1]=(m)[ 4]*(v1)[0]+(m)[ 5]*(v1)[1]+(m)[ 6]*(v1)[2]+(m)[ 7]*(v1)[3], \
***************
*** 565,568 ****
--- 559,568 ----
     (v2)[3]=(m)[12]*(v1)[0]+(m)[13]*(v1)[1]+(m)[14]*(v1)[2]+(m)[15]*(v1)[3])
  
+ #define ELL_4MV_TMUL(v2, m, v1)                                             \
+   ((v2)[0]=(m)[ 0]*(v1)[0]+(m)[ 4]*(v1)[1]+(m)[ 8]*(v1)[2]+(m)[12]*(v1)[3], \
+    (v2)[1]=(m)[ 1]*(v1)[0]+(m)[ 5]*(v1)[1]+(m)[ 9]*(v1)[2]+(m)[13]*(v1)[3], \
+    (v2)[2]=(m)[ 2]*(v1)[0]+(m)[ 6]*(v1)[1]+(m)[10]*(v1)[2]+(m)[14]*(v1)[3], \
+    (v2)[3]=(m)[ 3]*(v1)[0]+(m)[ 7]*(v1)[1]+(m)[11]*(v1)[2]+(m)[15]*(v1)[3])
+ 
  #define ELL_34V_HOMOG(v2, v1) \
    ((v2)[0] = (v1)[0]/(v1)[3], \
***************
*** 577,593 ****
  
  /*
- ** the ELL_4M_SET... macros are setting the matrix one _column_
- ** at a time- so the matrix components appear below in transpose
- **
  ** These macros are intended to be used as aids with homogeneous transforms
  */
  
! #define ELL_4M_COLS_SET(m, a, b, c, d)  \
!   (ELL_4V_COPY((m)+ 0, a),              \
!    ELL_4V_COPY((m)+ 4, b),              \
!    ELL_4V_COPY((m)+ 8, c),              \
!    ELL_4V_COPY((m)+12, d))
! 
! #define ELL_4M_ROWS_SET(m, a, b, c, d)                 \
    (ELL_4V_SET((m)+ 0, (a)[0], (b)[0], (c)[0], (d)[0]), \
     ELL_4V_SET((m)+ 4, (a)[1], (b)[1], (c)[1], (d)[1]), \
--- 577,584 ----
  
  /*
  ** These macros are intended to be used as aids with homogeneous transforms
  */
  
! #define ELL_4M_COLS_SET(m, a, b, c, d)                 \
    (ELL_4V_SET((m)+ 0, (a)[0], (b)[0], (c)[0], (d)[0]), \
     ELL_4V_SET((m)+ 4, (a)[1], (b)[1], (c)[1], (d)[1]), \
***************
*** 595,598 ****
--- 586,595 ----
     ELL_4V_SET((m)+12, (a)[3], (b)[3], (c)[3], (d)[3]))
  
+ #define ELL_4M_ROWS_SET(m, a, b, c, d)  \
+   (ELL_4V_COPY((m)+ 0, a),              \
+    ELL_4V_COPY((m)+ 4, b),              \
+    ELL_4V_COPY((m)+ 8, c),              \
+    ELL_4V_COPY((m)+12, d))
+ 
  #define ELL_4M_IDENTITY_SET(m) \
    (ELL_4V_SET((m)+ 0,  1 ,  0 ,  0 , 0), \
***************
*** 614,660 ****
  
  #define ELL_4M_TRANSLATE_SET(m, x, y, z) \
!   (ELL_4V_SET((m)+ 0,  1 ,  0 ,  0 , 0), \
!    ELL_4V_SET((m)+ 4,  0 ,  1 ,  0 , 0), \
!    ELL_4V_SET((m)+ 8,  0 ,  0 ,  1 , 0), \
!    ELL_4V_SET((m)+12, (x), (y), (z), 1))
  
  #define ELL_4M_ROTATE_X_SET(m, th)                   \
    (ELL_4V_SET((m)+ 0,  1 ,     0    ,     0    , 0), \
!    ELL_4V_SET((m)+ 4,  0 ,  cos(th) , +sin(th) , 0), \
!    ELL_4V_SET((m)+ 8,  0 , -sin(th) ,  cos(th) , 0), \
     ELL_4V_SET((m)+12,  0 ,     0    ,     0    , 1))
  
  #define ELL_4M_ROTATE_Y_SET(m, th)                   \
!   (ELL_4V_SET((m)+ 0,  cos(th) ,  0 , -sin(th) , 0), \
     ELL_4V_SET((m)+ 4,     0    ,  1 ,     0    , 0), \
!    ELL_4V_SET((m)+ 8, +sin(th) ,  0 ,  cos(th) , 0), \
     ELL_4V_SET((m)+12,     0    ,  0 ,     0    , 1))
  
  #define ELL_4M_ROTATE_Z_SET(m, th)                   \
!   (ELL_4V_SET((m)+ 0,  cos(th) , +sin(th) ,  0 , 0), \
!    ELL_4V_SET((m)+ 4, -sin(th) ,  cos(th) ,  0 , 0), \
     ELL_4V_SET((m)+ 8,     0    ,     0    ,  1 , 0), \
     ELL_4V_SET((m)+12,     0    ,     0    ,  0 , 1))
  
  #define ELL_4M_MUL(n, l, m)                                                 \
!   ((n)[ 0]=(l)[ 0]*(m)[ 0]+(l)[ 4]*(m)[ 1]+(l)[ 8]*(m)[ 2]+(l)[12]*(m)[ 3], \
!    (n)[ 1]=(l)[ 1]*(m)[ 0]+(l)[ 5]*(m)[ 1]+(l)[ 9]*(m)[ 2]+(l)[13]*(m)[ 3], \
!    (n)[ 2]=(l)[ 2]*(m)[ 0]+(l)[ 6]*(m)[ 1]+(l)[10]*(m)[ 2]+(l)[14]*(m)[ 3], \
!    (n)[ 3]=(l)[ 3]*(m)[ 0]+(l)[ 7]*(m)[ 1]+(l)[11]*(m)[ 2]+(l)[15]*(m)[ 3], \
                                                                              \
!    (n)[ 4]=(l)[ 0]*(m)[ 4]+(l)[ 4]*(m)[ 5]+(l)[ 8]*(m)[ 6]+(l)[12]*(m)[ 7], \
!    (n)[ 5]=(l)[ 1]*(m)[ 4]+(l)[ 5]*(m)[ 5]+(l)[ 9]*(m)[ 6]+(l)[13]*(m)[ 7], \
!    (n)[ 6]=(l)[ 2]*(m)[ 4]+(l)[ 6]*(m)[ 5]+(l)[10]*(m)[ 6]+(l)[14]*(m)[ 7], \
!    (n)[ 7]=(l)[ 3]*(m)[ 4]+(l)[ 7]*(m)[ 5]+(l)[11]*(m)[ 6]+(l)[15]*(m)[ 7], \
                                                                              \
!    (n)[ 8]=(l)[ 0]*(m)[ 8]+(l)[ 4]*(m)[ 9]+(l)[ 8]*(m)[10]+(l)[12]*(m)[11], \
!    (n)[ 9]=(l)[ 1]*(m)[ 8]+(l)[ 5]*(m)[ 9]+(l)[ 9]*(m)[10]+(l)[13]*(m)[11], \
!    (n)[10]=(l)[ 2]*(m)[ 8]+(l)[ 6]*(m)[ 9]+(l)[10]*(m)[10]+(l)[14]*(m)[11], \
!    (n)[11]=(l)[ 3]*(m)[ 8]+(l)[ 7]*(m)[ 9]+(l)[11]*(m)[10]+(l)[15]*(m)[11], \
                                                                              \
!    (n)[12]=(l)[ 0]*(m)[12]+(l)[ 4]*(m)[13]+(l)[ 8]*(m)[14]+(l)[12]*(m)[15], \
!    (n)[13]=(l)[ 1]*(m)[12]+(l)[ 5]*(m)[13]+(l)[ 9]*(m)[14]+(l)[13]*(m)[15], \
!    (n)[14]=(l)[ 2]*(m)[12]+(l)[ 6]*(m)[13]+(l)[10]*(m)[14]+(l)[14]*(m)[15], \
!    (n)[15]=(l)[ 3]*(m)[12]+(l)[ 7]*(m)[13]+(l)[11]*(m)[14]+(l)[15]*(m)[15])
  
  #define ELL_34M_EXTRACT(m, l) \
--- 611,657 ----
  
  #define ELL_4M_TRANSLATE_SET(m, x, y, z) \
!   (ELL_4V_SET((m)+ 0,  1 ,  0 ,  0 , (x)), \
!    ELL_4V_SET((m)+ 4,  0 ,  1 ,  0 , (y)), \
!    ELL_4V_SET((m)+ 8,  0 ,  0 ,  1 , (z)), \
!    ELL_4V_SET((m)+12,  0 ,  0 ,  0 ,  1))
  
  #define ELL_4M_ROTATE_X_SET(m, th)                   \
    (ELL_4V_SET((m)+ 0,  1 ,     0    ,     0    , 0), \
!    ELL_4V_SET((m)+ 4,  0 ,  cos(th) , -sin(th) , 0), \
!    ELL_4V_SET((m)+ 8,  0 , +sin(th) ,  cos(th) , 0), \
     ELL_4V_SET((m)+12,  0 ,     0    ,     0    , 1))
  
  #define ELL_4M_ROTATE_Y_SET(m, th)                   \
!   (ELL_4V_SET((m)+ 0,  cos(th) ,  0 , +sin(th) , 0), \
     ELL_4V_SET((m)+ 4,     0    ,  1 ,     0    , 0), \
!    ELL_4V_SET((m)+ 8, -sin(th) ,  0 ,  cos(th) , 0), \
     ELL_4V_SET((m)+12,     0    ,  0 ,     0    , 1))
  
  #define ELL_4M_ROTATE_Z_SET(m, th)                   \
!   (ELL_4V_SET((m)+ 0,  cos(th) , -sin(th) ,  0 , 0), \
!    ELL_4V_SET((m)+ 4, +sin(th) ,  cos(th) ,  0 , 0), \
     ELL_4V_SET((m)+ 8,     0    ,     0    ,  1 , 0), \
     ELL_4V_SET((m)+12,     0    ,     0    ,  0 , 1))
  
  #define ELL_4M_MUL(n, l, m)                                                 \
!   ((n)[ 0]=(l)[ 0]*(m)[ 0]+(l)[ 1]*(m)[ 4]+(l)[ 2]*(m)[ 8]+(l)[ 3]*(m)[12], \
!    (n)[ 1]=(l)[ 0]*(m)[ 1]+(l)[ 1]*(m)[ 5]+(l)[ 2]*(m)[ 9]+(l)[ 3]*(m)[13], \
!    (n)[ 2]=(l)[ 0]*(m)[ 2]+(l)[ 1]*(m)[ 6]+(l)[ 2]*(m)[10]+(l)[ 3]*(m)[14], \
!    (n)[ 3]=(l)[ 0]*(m)[ 3]+(l)[ 1]*(m)[ 7]+(l)[ 2]*(m)[11]+(l)[ 3]*(m)[15], \
                                                                              \
!    (n)[ 4]=(l)[ 4]*(m)[ 0]+(l)[ 5]*(m)[ 4]+(l)[ 6]*(m)[ 8]+(l)[ 7]*(m)[12], \
!    (n)[ 5]=(l)[ 4]*(m)[ 1]+(l)[ 5]*(m)[ 5]+(l)[ 6]*(m)[ 9]+(l)[ 7]*(m)[13], \
!    (n)[ 6]=(l)[ 4]*(m)[ 2]+(l)[ 5]*(m)[ 6]+(l)[ 6]*(m)[10]+(l)[ 7]*(m)[14], \
!    (n)[ 7]=(l)[ 4]*(m)[ 3]+(l)[ 5]*(m)[ 7]+(l)[ 6]*(m)[11]+(l)[ 7]*(m)[15], \
                                                                              \
!    (n)[ 8]=(l)[ 8]*(m)[ 0]+(l)[ 9]*(m)[ 4]+(l)[10]*(m)[ 8]+(l)[11]*(m)[12], \
!    (n)[ 9]=(l)[ 8]*(m)[ 1]+(l)[ 9]*(m)[ 5]+(l)[10]*(m)[ 9]+(l)[11]*(m)[13], \
!    (n)[10]=(l)[ 8]*(m)[ 2]+(l)[ 9]*(m)[ 6]+(l)[10]*(m)[10]+(l)[11]*(m)[14], \
!    (n)[11]=(l)[ 8]*(m)[ 3]+(l)[ 9]*(m)[ 7]+(l)[10]*(m)[11]+(l)[11]*(m)[15], \
                                                                              \
!    (n)[12]=(l)[12]*(m)[ 0]+(l)[13]*(m)[ 4]+(l)[14]*(m)[ 8]+(l)[15]*(m)[12], \
!    (n)[13]=(l)[12]*(m)[ 1]+(l)[13]*(m)[ 5]+(l)[14]*(m)[ 9]+(l)[15]*(m)[13], \
!    (n)[14]=(l)[12]*(m)[ 2]+(l)[13]*(m)[ 6]+(l)[14]*(m)[10]+(l)[15]*(m)[14], \
!    (n)[15]=(l)[12]*(m)[ 3]+(l)[13]*(m)[ 7]+(l)[14]*(m)[11]+(l)[15]*(m)[15])
  
  #define ELL_34M_EXTRACT(m, l) \
***************
*** 679,691 ****
                             (m)[ 9], (m)[10], (m)[11], \
                             (m)[13], (m)[14], (m)[15]) \
!    - (m)[ 4] * _ELL_3M_DET((m)[ 1], (m)[ 2], (m)[ 3], \
!                            (m)[ 9], (m)[10], (m)[11], \
!                            (m)[13], (m)[14], (m)[15]) \
!    + (m)[ 8] * _ELL_3M_DET((m)[ 1], (m)[ 2], (m)[ 3], \
!                            (m)[ 5], (m)[ 6], (m)[ 7], \
!                            (m)[13], (m)[14], (m)[15]) \
!    - (m)[12] * _ELL_3M_DET((m)[ 1], (m)[ 2], (m)[ 3], \
!                            (m)[ 5], (m)[ 6], (m)[ 7], \
!                            (m)[ 9], (m)[10], (m)[11]))
  
  #define ELL_Q_MUL(q3, q1, q2)                                            \
--- 676,688 ----
                             (m)[ 9], (m)[10], (m)[11], \
                             (m)[13], (m)[14], (m)[15]) \
!    - (m)[ 1] * _ELL_3M_DET((m)[ 4], (m)[ 6], (m)[ 7], \
!                            (m)[ 8], (m)[10], (m)[11], \
!                            (m)[12], (m)[14], (m)[15]) \
!    + (m)[ 2] * _ELL_3M_DET((m)[ 4], (m)[ 5], (m)[ 7], \
!                            (m)[ 8], (m)[ 9], (m)[11], \
!                            (m)[12], (m)[13], (m)[15]) \
!    - (m)[ 3] * _ELL_3M_DET((m)[ 4], (m)[ 5], (m)[ 6], \
!                            (m)[ 8], (m)[ 9], (m)[10], \
!                            (m)[12], (m)[13], (m)[14]))
  
  #define ELL_Q_MUL(q3, q1, q2)                                            \
***************
*** 704,734 ****
  
  #define ELL_Q_TO_3M(m, q)                                                   \
!   (ELL_3V_SET((m)+0,                                                        \
               (q)[0]*(q)[0] + (q)[1]*(q)[1] - (q)[2]*(q)[2] - (q)[3]*(q)[3], \
-              2*((q)[1]*(q)[2] + (q)[0]*(q)[3]),                             \
-              2*((q)[1]*(q)[3] - (q)[0]*(q)[2])),                            \
-   ELL_3V_SET((m)+3,                                                         \
               2*((q)[1]*(q)[2] - (q)[0]*(q)[3]),                             \
               (q)[0]*(q)[0] - (q)[1]*(q)[1] + (q)[2]*(q)[2] - (q)[3]*(q)[3], \
!              2*((q)[2]*(q)[3] + (q)[0]*(q)[1])),                            \
    ELL_3V_SET((m)+6,                                                         \
!              2*((q)[1]*(q)[3] + (q)[0]*(q)[2]),                             \
!              2*((q)[2]*(q)[3] - (q)[0]*(q)[1]),                             \
               (q)[0]*(q)[0] - (q)[1]*(q)[1] - (q)[2]*(q)[2] + (q)[3]*(q)[3]))
  
  #define ELL_Q_TO_4M(m, q)                                                   \
!   (ELL_4V_SET((m)+0,                                                        \
               (q)[0]*(q)[0] + (q)[1]*(q)[1] - (q)[2]*(q)[2] - (q)[3]*(q)[3], \
!              2*((q)[1]*(q)[2] + (q)[0]*(q)[3]),                             \
!              2*((q)[1]*(q)[3] - (q)[0]*(q)[2]),                             \
               0),                                                            \
    ELL_4V_SET((m)+4,                                                         \
!              2*((q)[1]*(q)[2] - (q)[0]*(q)[3]),                             \
               (q)[0]*(q)[0] - (q)[1]*(q)[1] + (q)[2]*(q)[2] - (q)[3]*(q)[3], \
!              2*((q)[2]*(q)[3] + (q)[0]*(q)[1]),                             \
               0),                                                            \
    ELL_4V_SET((m)+8,                                                         \
!              2*((q)[1]*(q)[3] + (q)[0]*(q)[2]),                             \
!              2*((q)[2]*(q)[3] - (q)[0]*(q)[1]),                             \
               (q)[0]*(q)[0] - (q)[1]*(q)[1] - (q)[2]*(q)[2] + (q)[3]*(q)[3], \
               0),                                                            \
--- 701,731 ----
  
  #define ELL_Q_TO_3M(m, q)                                                   \
!  (ELL_3V_SET((m)+0,                                                         \
               (q)[0]*(q)[0] + (q)[1]*(q)[1] - (q)[2]*(q)[2] - (q)[3]*(q)[3], \
               2*((q)[1]*(q)[2] - (q)[0]*(q)[3]),                             \
+              2*((q)[1]*(q)[3] + (q)[0]*(q)[2])),                            \
+   ELL_3V_SET((m)+3,                                                         \
+              2*((q)[1]*(q)[2] + (q)[0]*(q)[3]),                             \
               (q)[0]*(q)[0] - (q)[1]*(q)[1] + (q)[2]*(q)[2] - (q)[3]*(q)[3], \
!              2*((q)[2]*(q)[3] - (q)[0]*(q)[1])),                            \
    ELL_3V_SET((m)+6,                                                         \
!              2*((q)[1]*(q)[3] - (q)[0]*(q)[2]),                             \
!              2*((q)[2]*(q)[3] + (q)[0]*(q)[1]),                             \
               (q)[0]*(q)[0] - (q)[1]*(q)[1] - (q)[2]*(q)[2] + (q)[3]*(q)[3]))
  
  #define ELL_Q_TO_4M(m, q)                                                   \
!  (ELL_4V_SET((m)+0,                                                         \
               (q)[0]*(q)[0] + (q)[1]*(q)[1] - (q)[2]*(q)[2] - (q)[3]*(q)[3], \
!              2*((q)[1]*(q)[2] - (q)[0]*(q)[3]),                             \
!              2*((q)[1]*(q)[3] + (q)[0]*(q)[2]),                             \
               0),                                                            \
    ELL_4V_SET((m)+4,                                                         \
!              2*((q)[1]*(q)[2] + (q)[0]*(q)[3]),                             \
               (q)[0]*(q)[0] - (q)[1]*(q)[1] + (q)[2]*(q)[2] - (q)[3]*(q)[3], \
!              2*((q)[2]*(q)[3] - (q)[0]*(q)[1]),                             \
               0),                                                            \
    ELL_4V_SET((m)+8,                                                         \
!              2*((q)[1]*(q)[3] - (q)[0]*(q)[2]),                             \
!              2*((q)[2]*(q)[3] + (q)[0]*(q)[1]),                             \
               (q)[0]*(q)[0] - (q)[1]*(q)[1] - (q)[2]*(q)[2] + (q)[3]*(q)[3], \
               0),                                                            \

Index: genmat.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ell/genmat.c,v
retrieving revision 1.4
retrieving revision 1.5
diff -C2 -d -r1.4 -r1.5
*** genmat.c	7 Jan 2004 15:34:29 -0000	1.4
--- genmat.c	1 Mar 2004 10:19:44 -0000	1.5
***************
*** 21,56 ****
  #include "ell.h"
  
- /*
- ** ========================
- ** NOTE NOTE NOTE NOTE NOTE
- ** ========================
- ** All the functions in this file use a matrix ordering which
- ** is opposite if the (current) ordering in the rest of ell.
- **
- ** Rest of ell:
- **  0 . . .
- **  1 . . .
- **  2 . . .
- **  . . . .
- **
- ** This file:
- **  0 1 2 .
- **  . . . .
- **  . . . .
- **  . . . .
- */
- 
  int
! ell_Nm_valid(Nrrd *mat) {
!   char me[]="ell_Nm_valid", err[AIR_STRLEN_MED];
  
!   if (!mat) {
!     sprintf(err, "%s: got NULL pointer", me);
!     biffAdd(ELL, err); return 0;
    }
-   /* nrrdCheck(mat)? */
    if (!( 2 == mat->dim )) {
      sprintf(err, "%s: nrrd must be 2-D (not %d-D)", me, mat->dim);
!     biffAdd(ELL, err); return 0;
    }
    if (!( nrrdTypeDouble == mat->type )) {
--- 21,35 ----
  #include "ell.h"
  
  int
! ell_Nm_check(Nrrd *mat) {
!   char me[]="ell_Nm_check", err[AIR_STRLEN_MED];
  
!   if (nrrdCheck(mat)) {
!     sprintf(err, "%s: basic nrrd validity check failed", me);
!     biffMove(ELL, err, NRRD); return 1;
    }
    if (!( 2 == mat->dim )) {
      sprintf(err, "%s: nrrd must be 2-D (not %d-D)", me, mat->dim);
!     biffAdd(ELL, err); return 1;
    }
    if (!( nrrdTypeDouble == mat->type )) {
***************
*** 58,65 ****
  	    airEnumStr(nrrdType, nrrdTypeDouble),
  	    airEnumStr(nrrdType, mat->type));
!     biffAdd(ELL, err); return 0;
    }
  
!   return 1;
  }
  
--- 37,44 ----
  	    airEnumStr(nrrdType, nrrdTypeDouble),
  	    airEnumStr(nrrdType, mat->type));
!     biffAdd(ELL, err); return 1;
    }
  
!   return 0;
  }
  
***************
*** 68,72 ****
    char me[]="ell_Nm_tran", err[AIR_STRLEN_MED];
  
!   if (!( ntrn && ell_Nm_valid(nmat) )) {
      sprintf(err, "%s: NULL or invalid args", me);
      biffAdd(ELL, err); return 1;
--- 47,51 ----
    char me[]="ell_Nm_tran", err[AIR_STRLEN_MED];
  
!   if (!( ntrn && !ell_Nm_check(nmat) )) {
      sprintf(err, "%s: NULL or invalid args", me);
      biffAdd(ELL, err); return 1;
***************
*** 98,102 ****
    int LL, MM, NN, ll, mm, nn;
    
!   if (!( nAB && ell_Nm_valid(nA) && ell_Nm_valid(nB) )) {
      sprintf(err, "%s: NULL or invalid args", me);
      biffAdd(ELL, err); return 1;
--- 77,81 ----
    int LL, MM, NN, ll, mm, nn;
    
!   if (!( nAB && !ell_Nm_check(nA) && !ell_Nm_check(nB) )) {
      sprintf(err, "%s: NULL or invalid args", me);
      biffAdd(ELL, err); return 1;
***************
*** 313,317 ****
    int NN;
  
!   if (!( ninv && ell_Nm_valid(nmat) )) {
      sprintf(err, "%s: NULL or invalid args", me);
      biffAdd(ELL, err); return 1;
--- 292,296 ----
    int NN;
  
!   if (!( ninv && !ell_Nm_check(nmat) )) {
      sprintf(err, "%s: NULL or invalid args", me);
      biffAdd(ELL, err); return 1;
***************
*** 353,357 ****
    int ret=0;
    
!   if (!( ninv && ell_Nm_valid(nA) )) {
      sprintf(err, "%s: NULL or invalid args", me);
      biffAdd(ELL, err); return 1;
--- 332,336 ----
    int ret=0;
    
!   if (!( ninv && !ell_Nm_check(nA) )) {
      sprintf(err, "%s: NULL or invalid args", me);
      biffAdd(ELL, err); return 1;

Index: miscEll.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ell/miscEll.c,v
retrieving revision 1.5
retrieving revision 1.6
diff -C2 -d -r1.5 -r1.6
*** miscEll.c	7 Jan 2004 15:34:29 -0000	1.5
--- miscEll.c	1 Mar 2004 10:19:44 -0000	1.6
***************
*** 37,45 ****
  
    fprintf(f, "% 15.7f % 15.7f % 15.7f\n", 
! 	  s[0], s[3], s[6]);
    fprintf(f, "% 15.7f % 15.7f % 15.7f\n", 
! 	  s[1], s[4], s[7]);
    fprintf(f, "% 15.7f % 15.7f % 15.7f\n", 
! 	  s[2], s[5], s[8]);
  }
  
--- 37,45 ----
  
    fprintf(f, "% 15.7f % 15.7f % 15.7f\n", 
! 	  s[0], s[1], s[2]);
    fprintf(f, "% 15.7f % 15.7f % 15.7f\n", 
! 	  s[3], s[4], s[5]);
    fprintf(f, "% 15.7f % 15.7f % 15.7f\n", 
! 	  s[6], s[7], s[8]);
  }
  
***************
*** 55,63 ****
  
    fprintf(f, "% 31.15f % 31.15f % 31.15f\n", 
! 	  s[0], s[3], s[6]);
    fprintf(f, "% 31.15f % 31.15f % 31.15f\n", 
! 	  s[1], s[4], s[7]);
    fprintf(f, "% 31.15f % 31.15f % 31.15f\n", 
! 	  s[2], s[5], s[8]);
  }
  
--- 55,63 ----
  
    fprintf(f, "% 31.15f % 31.15f % 31.15f\n", 
! 	  s[0], s[1], s[2]);
    fprintf(f, "% 31.15f % 31.15f % 31.15f\n", 
! 	  s[3], s[4], s[5]);
    fprintf(f, "% 31.15f % 31.15f % 31.15f\n", 
! 	  s[6], s[7], s[8]);
  }
  
***************
*** 73,83 ****
  
    fprintf(f, "% 15.7f % 15.7f % 15.7f % 15.7f\n", 
! 	  s[0], s[4], s[8], s[12]);
    fprintf(f, "% 15.7f % 15.7f % 15.7f % 15.7f\n", 
! 	  s[1], s[5], s[9], s[13]);
    fprintf(f, "% 15.7f % 15.7f % 15.7f % 15.7f\n", 
! 	  s[2], s[6], s[10], s[14]);
    fprintf(f, "% 15.7f % 15.7f % 15.7f % 15.7f\n", 
! 	  s[3], s[7], s[11], s[15]);
  }
  
--- 73,83 ----
  
    fprintf(f, "% 15.7f % 15.7f % 15.7f % 15.7f\n", 
! 	  s[ 0], s[ 1], s[ 2], s[ 3]);
    fprintf(f, "% 15.7f % 15.7f % 15.7f % 15.7f\n", 
! 	  s[ 4], s[ 5], s[ 6], s[ 7]);
    fprintf(f, "% 15.7f % 15.7f % 15.7f % 15.7f\n", 
! 	  s[ 8], s[ 9], s[10], s[11]);
    fprintf(f, "% 15.7f % 15.7f % 15.7f % 15.7f\n", 
! 	  s[12], s[13], s[14], s[15]);
  }
  
***************
*** 93,103 ****
  
    fprintf(f, "% 31.15f % 31.15f % 31.15f % 31.15f\n", 
! 	  s[0], s[4], s[8], s[12]);
    fprintf(f, "% 31.15f % 31.15f % 31.15f % 31.15f\n", 
! 	  s[1], s[5], s[9], s[13]);
    fprintf(f, "% 31.15f % 31.15f % 31.15f % 31.15f\n", 
! 	  s[2], s[6], s[10], s[14]);
    fprintf(f, "% 31.15f % 31.15f % 31.15f % 31.15f\n", 
! 	  s[3], s[7], s[11], s[15]);
  }
  
--- 93,103 ----
  
    fprintf(f, "% 31.15f % 31.15f % 31.15f % 31.15f\n", 
! 	  s[ 0], s[ 1], s[ 2], s[ 3]);
    fprintf(f, "% 31.15f % 31.15f % 31.15f % 31.15f\n", 
! 	  s[ 4], s[ 5], s[ 6], s[ 7]);
    fprintf(f, "% 31.15f % 31.15f % 31.15f % 31.15f\n", 
! 	  s[ 8], s[ 9], s[10], s[11]);
    fprintf(f, "% 31.15f % 31.15f % 31.15f % 31.15f\n", 
! 	  s[12], s[13], s[14], s[15]);
  }
  

Index: quat.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ell/quat.c,v
retrieving revision 1.3
retrieving revision 1.4
diff -C2 -d -r1.3 -r1.4
*** quat.c	7 Jan 2004 15:34:29 -0000	1.3
--- quat.c	1 Mar 2004 10:19:44 -0000	1.4
***************
*** 26,37 ****
  
  /*
!  0  3  6
!  1  4  7
!  2  5  8 
  
!  0   4   8  12
!  1   5   9  13
!  2   6  10  14
!  3   7  11  15
  */
  
--- 26,37 ----
  
  /*
!  0  1  2
!  3  4  5
!  6  7  8
  
!  0   1   2   3
!  4   5   6   7
!  8   9  10  11
! 12  13  14  15
  */
  
***************
*** 51,60 ****
    s[Y] = 1 - m[i0] + m[i4] - m[i8];       \
    s[Z] = 1 - m[i0] - m[i4] + m[i8];       \
!   wx = m[i5] - m[i7];                     \
!   wy = m[i6] - m[i2];                     \
!   wz = m[i1] - m[i3];                     \
!   xy = m[i1] + m[i3];                     \
!   xz = m[i2] + m[i6];                     \
!   yz = m[i5] + m[i7];                     \
    mi =  s[W] > s[X] ?  W : X;             \
    mi = s[mi] > s[Y] ? mi : Y;             \
--- 51,60 ----
    s[Y] = 1 - m[i0] + m[i4] - m[i8];       \
    s[Z] = 1 - m[i0] - m[i4] + m[i8];       \
!   wx = m[i7] - m[i5];                     \
!   wy = m[i2] - m[i6];                     \
!   wz = m[i3] - m[i1];                     \
!   xy = m[i3] + m[i1];                     \
!   xz = m[i6] + m[i2];                     \
!   yz = m[i7] + m[i5];                     \
    mi =  s[W] > s[X] ?  W : X;             \
    mi = s[mi] > s[Y] ? mi : Y;             \
***************
*** 109,121 ****
    ELL_3V_SET(m+0,                    \
               w*w + x*x - y*y - z*z,  \
-              2*(x*y + w*z),          \
-              2*(x*z - w*y));         \
-   ELL_3V_SET(m+3,                    \
               2*(x*y - w*z),          \
               w*w - x*x + y*y - z*z,  \
!              2*(y*z + w*x));         \
    ELL_3V_SET(m+6,                    \
!              2*(x*z + w*y),          \
!              2*(y*z - w*x),          \
               w*w - x*x - y*y + z*z)
  
--- 109,121 ----
    ELL_3V_SET(m+0,                    \
               w*w + x*x - y*y - z*z,  \
               2*(x*y - w*z),          \
+              2*(x*z + w*y));         \
+   ELL_3V_SET(m+3,                    \
+              2*(x*y + w*z),          \
               w*w - x*x + y*y - z*z,  \
!              2*(y*z - w*x));         \
    ELL_3V_SET(m+6,                    \
!              2*(x*z - w*y),          \
!              2*(y*z + w*x),          \
               w*w - x*x - y*y + z*z)
  
***************
*** 137,151 ****
    ELL_4V_SET(m+0,                    \
               w*w + x*x - y*y - z*z,  \
!              2*(x*y + w*z),          \
!              2*(x*z - w*y),          \
               0);                     \
    ELL_4V_SET(m+4,                    \
!              2*(x*y - w*z),          \
               w*w - x*x + y*y - z*z,  \
!              2*(y*z + w*x),          \
               0);                     \
    ELL_4V_SET(m+8,                    \
!              2*(x*z + w*y),          \
!              2*(y*z - w*x),          \
               w*w - x*x - y*y + z*z,  \
               0);                     \
--- 137,151 ----
    ELL_4V_SET(m+0,                    \
               w*w + x*x - y*y - z*z,  \
!              2*(x*y - w*z),          \
!              2*(x*z + w*y),          \
               0);                     \
    ELL_4V_SET(m+4,                    \
!              2*(x*y + w*z),          \
               w*w - x*x + y*y - z*z,  \
!              2*(y*z - w*x),          \
               0);                     \
    ELL_4V_SET(m+8,                    \
!              2*(x*z - w*y),          \
!              2*(y*z + w*x),          \
               w*w - x*x - y*y + z*z,  \
               0);                     \

[Teem-commit] teem/src/gage sclanswer.c,1.19,1.20

[<kin...@us...>]

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

Modified Files:
	sclanswer.c 
Log Message:
ell's storage of matrices is now row-major instead of column major.  Debugging in progress

Index: sclanswer.c
===================================================================
RCS file: /cvsroot/teem/teem/src/gage/sclanswer.c,v
retrieving revision 1.19
retrieving revision 1.20
diff -C2 -d -r1.19 -r1.20
*** sclanswer.c	28 Feb 2004 17:01:53 -0000	1.19
--- sclanswer.c	1 Mar 2004 10:19:45 -0000	1.20
***************
*** 119,123 ****
        ELL_3M_SCALE(sHess, -(ctx->parm.curvNormalSide)/gmag, hess);
        
!       /* san->gten = nPerp * sHess * nPerp */
        ELL_3M_MUL(tmpMat, sHess, nPerp);
        ELL_3M_MUL(gten, nPerp, tmpMat);
--- 119,123 ----
        ELL_3M_SCALE(sHess, -(ctx->parm.curvNormalSide)/gmag, hess);
        
!       /* gten = nPerp * sHess * nPerp */
        ELL_3M_MUL(tmpMat, sHess, nPerp);
        ELL_3M_MUL(gten, nPerp, tmpMat);

[Teem-commit] teem/src/ell/test invert.c,1.4,1.5

[<kin...@us...>]

Update of /cvsroot/teem/teem/src/ell/test
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv1922/ell/test

Modified Files:
	invert.c 
Log Message:
ell's storage of matrices is now row-major instead of column major.  Debugging in progress

Index: invert.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ell/test/invert.c,v
retrieving revision 1.4
retrieving revision 1.5
diff -C2 -d -r1.4 -r1.5
*** invert.c	7 Jan 2004 15:34:29 -0000	1.4
--- invert.c	1 Mar 2004 10:19:45 -0000	1.5
***************
*** 20,24 ****
  #include "../ell.h"
  
! char *invInfo = "Tests ell_Nm_inv and ell_Nm_pseudo_inv";
  
  int
--- 20,25 ----
  #include "../ell.h"
  
! char *invInfo = ("Tests ell_Nm_inv and ell_Nm_pseudo_inv, "
! 		 "and ell_{3,4}inv_d where possible ");
  
  int
***************
*** 31,34 ****
--- 32,36 ----
    int (*func)(Nrrd *, Nrrd *);
  
+   double m3[9], m4[16];
    me = argv[0];
    mop = airMopNew();
***************
*** 54,57 ****
--- 56,73 ----
    
    nrrdConvert(nmat, nin, nrrdTypeDouble);
+   if (3 == nmat->axis[0].size && 3 == nmat->axis[1].size) {
+     ell_3m_inv_d(m3, nmat->data);
+     fprintf(stderr, "%s: input:\n", me);
+     ell_3m_print_d(stderr, nmat->data);
+     fprintf(stderr, "%s: inverse:\n", me);
+     ell_3m_print_d(stderr, m3);
+   }
+   if (4 == nmat->axis[0].size && 4 == nmat->axis[1].size) {
+     ell_4m_inv_d(m4, nmat->data);
+     fprintf(stderr, "%s: input:\n", me);
+     ell_4m_print_d(stderr, nmat->data);
+     fprintf(stderr, "%s: inverse:\n", me);
+     ell_4m_print_d(stderr, m4);
+   }
    func = (nmat->axis[0].size == nmat->axis[1].size
  	  ? ell_Nm_inv

[Teem-commit] teem/src/bin emap.c,1.9,1.10

[<kin...@us...>]

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

Modified Files:
	emap.c 
Log Message:
added -err option for looking at error images

Index: emap.c
===================================================================
RCS file: /cvsroot/teem/teem/src/bin/emap.c,v
retrieving revision 1.9
retrieving revision 1.10
diff -C2 -d -r1.9 -r1.10
*** emap.c	27 Feb 2004 18:16:50 -0000	1.9
--- emap.c	1 Mar 2004 10:17:43 -0000	1.10
***************
*** 41,46 ****
    char *me, *outS, *errS, *debugS;
    airArray *mop;
!   float amb[3], *linfo, *debug, *map, W[3], V[3] /* , diff */;
!   int li, ui, vi, qn, bits, method;
    limnLight *light;
    limnCamera *cam;
--- 41,46 ----
    char *me, *outS, *errS, *debugS;
    airArray *mop;
!   float amb[3], *linfo, *debug, *map, W[3], V[3], diff;
!   int li, ui, vi, qn, bits, method, doerr;
    limnLight *light;
    limnCamera *cam;
***************
*** 78,81 ****
--- 78,84 ----
  	     "right-ward, V increases downward.  The back sphere half is "
  	     "rendered as though the front half was removed");
+   hestOptAdd(&hopt, "err", NULL, airTypeInt, 0, 0, &doerr, NULL,
+ 	     "If using \"-d\", make the image represent the error between the "
+ 	     "real and quantized vector");
    hestParseOrDie(hopt, argc-1, argv+1, hparm, me, emapInfo,
  		 AIR_TRUE, AIR_TRUE, AIR_TRUE);
***************
*** 162,186 ****
  	ELL_3MV_MUL(W, V2W, V);
  	qn = limnVtoQN_f[method](W);
! 	/* 
! 	limnQNtoV_f[method](V, qn);
! 	ELL_3V_SUB(W, W, V);
! 	diff = ELL_3V_LEN(W);
! 	ELL_3V_SET(debug + 3*(ui + 1024*vi), diff, diff, diff);
! 	*/
! 	ELL_3V_COPY(debug + 3*(ui + 1024*vi), map + 3*qn);
  
- 	
  	/* second, the far side of the sphere */
  	ELL_3V_SET(V, u, v, w);
  	ELL_3MV_MUL(W, V2W, V);
  	qn = limnVtoQN_f[method](W);
! 	/*
! 	limnQNtoV_f[method](V, qn);
! 	ELL_3V_SUB(W, W, V);
! 	diff = ELL_3V_LEN(W);
! 	ELL_3V_SET(debug + 3*(ui + 512 + 1024*vi), diff, diff, diff);
! 	*/
! 	ELL_3V_COPY(debug + 3*(ui + 512 + 1024*vi), map + 3*qn);
! 
        }
      }
--- 165,189 ----
  	ELL_3MV_MUL(W, V2W, V);
  	qn = limnVtoQN_f[method](W);
! 	if (doerr) {
! 	  limnQNtoV_f[method](V, qn);
! 	  ELL_3V_SUB(W, W, V);
! 	  diff = ELL_3V_LEN(W);
! 	  ELL_3V_SET(debug + 3*(ui + 1024*vi), diff, diff, diff);
! 	} else {
! 	  ELL_3V_COPY(debug + 3*(ui + 1024*vi), map + 3*qn);
! 	}
  
  	/* second, the far side of the sphere */
  	ELL_3V_SET(V, u, v, w);
  	ELL_3MV_MUL(W, V2W, V);
  	qn = limnVtoQN_f[method](W);
! 	if (doerr) {
! 	  limnQNtoV_f[method](V, qn);
! 	  ELL_3V_SUB(W, W, V);
! 	  diff = ELL_3V_LEN(W);
! 	  ELL_3V_SET(debug + 3*(ui + 512 + 1024*vi), diff, diff, diff);
! 	} else {
! 	  ELL_3V_COPY(debug + 3*(ui + 512 + 1024*vi), map + 3*qn);
! 	}
        }
      }

[Teem-commit] teem/src/mite thread.c,1.11,1.12

[<kin...@us...>]

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

Modified Files:
	thread.c 
Log Message:
fixed memory leak in miteThreadNix

Index: thread.c
===================================================================
RCS file: /cvsroot/teem/teem/src/mite/thread.c,v
retrieving revision 1.11
retrieving revision 1.12
diff -C2 -d -r1.11 -r1.12
*** thread.c	29 Feb 2004 13:42:41 -0000	1.11
--- thread.c	29 Feb 2004 15:49:10 -0000	1.12
***************
*** 69,72 ****
--- 69,73 ----
  
    AIR_FREE(mtt->ansMiteVal);
+   AIR_FREE(mtt->directAnsMiteVal);
    AIR_FREE(mtt);
    return NULL;

[Teem-commit] teem/src/gage kind.c,1.1,1.2

[<kin...@us...>]

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

Modified Files:
	kind.c 
Log Message:
fixed big bug in _gageKindAnswerOffset

Index: kind.c
===================================================================
RCS file: /cvsroot/teem/teem/src/gage/kind.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -C2 -d -r1.1 -r1.2
*** kind.c	13 Feb 2004 23:07:15 -0000	1.1
--- kind.c	29 Feb 2004 15:48:38 -0000	1.2
***************
*** 62,65 ****
--- 62,77 ----
      pindex = kind->table[ii].parentIndex;
      if (-1 != pitem) {
+       if (0 == ii) {
+ 	sprintf(err, "%s: first item (index 0) of kind \"%s\" can't "
+ 		"be a sub-item (wanted parent index %d)", 
+ 		me, kind->name, pitem);
+ 	biffAdd(GAGE, err); return 1;
+       }
+       if (!(AIR_IN_CL(0, pitem, kind->itemMax))) {
+ 	sprintf(err, "%s: item %d of kind \"%s\" wants parent item %d "
+ 		"outside valid range [0..%d]",
+ 		me, ii, kind->name, pitem, kind->itemMax);
+ 	biffAdd(GAGE, err); return 1;
+       }
        if (-1 != kind->table[pitem].parentItem) {
  	sprintf(err, "%s: item %d of kind \"%s\" has parent %d which "
***************
*** 101,111 ****
  }
  
  int
  _gageKindAnswerOffset(gageKind *kind, int item) {
  
!   return (item
! 	  ? (kind->table[item-1].answerLength 
! 	     + _gageKindAnswerOffset(kind, item-1))
! 	  : 0);
  }
  
--- 113,148 ----
  }
  
+ /*
+ ** _gageKindAnswerOffset
+ **
+ ** return the location of the item in the master answer array 
+ **
+ ** I don't think this will work if there are sub-sub-items
+ */
  int
  _gageKindAnswerOffset(gageKind *kind, int item) {
+   int parent, ii;
  
!   if (!item) {
!     /* the first item always has zero offset */
!     return 0;
!   }
! 
!   /* else we're not the first */
!   parent = kind->table[item].parentItem;
!   if (-1 != parent) {
!     /* we're a sub-item */
!     return (kind->table[item].parentIndex 
! 	    + _gageKindAnswerOffset(kind, parent));
!   }
! 
!   /* else we're not a sub-item: find the first previous non-sub-item */
!   ii = item-1;
!   while (-1 != kind->table[ii].parentItem) {
!     /* gageKindCheck ensures that item 0 is not a sub-item */
!     ii--;
!   }
!   return (kind->table[ii].answerLength
! 	  + _gageKindAnswerOffset(kind, ii));
  }
  
***************
*** 120,128 ****
    }
  
!   /* here's why its important that there are no sub-sub-items */
!   return (-1 == kind->table[item].parentItem
! 	  ? _gageKindAnswerOffset(kind, item)
! 	  : (kind->table[item].parentIndex
! 	     + _gageKindAnswerOffset(kind, kind->table[item].parentItem)));
  }
  
--- 157,161 ----
    }
  
!   return _gageKindAnswerOffset(kind, item);
  }
  

[Teem-commit] teem/src/mite kindnot.c,1.1,1.2 mite.h,1.27,1.28 ray.c,1.18,1.19 renderMite.c,1.11,1.12 shade.c,1.1,1.2 thread.c,1.10,1.11 txf.c,1.17,1.18 user.c,1.15,1.16

[<kin...@us...>]

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

Modified Files:
	kindnot.c mite.h ray.c renderMite.c shade.c thread.c txf.c 
	user.c 
Log Message:
attempt to fix and optimize computation of the miteVals per-thread

Index: kindnot.c
===================================================================
RCS file: /cvsroot/teem/teem/src/mite/kindnot.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -C2 -d -r1.1 -r1.2
*** kindnot.c	26 Feb 2004 06:41:18 -0000	1.1
--- kindnot.c	29 Feb 2004 13:42:41 -0000	1.2
***************
*** 109,115 ****
  
  /*
! ** again, this is not a true gageKind- mainly because these items
  ** depend on items in different gageKinds (scalar and vector).  So,
! ** the prerequisites here are all blank.  Go look in _miteQueryAdd()
  ** to see these items' true prereqs
  */
--- 109,115 ----
  
  /*
! ** again, this is not a true gageKind- mainly because these items may
  ** depend on items in different gageKinds (scalar and vector).  So,
! ** the prerequisites here are all blank.  Go look in miteQueryAdd()
  ** to see these items' true prereqs
  */

Index: mite.h
===================================================================
RCS file: /cvsroot/teem/teem/src/mite/mite.h,v
retrieving revision 1.27
retrieving revision 1.28
diff -C2 -d -r1.27 -r1.28
*** mite.h	26 Feb 2004 06:26:10 -0000	1.27
--- mite.h	29 Feb 2004 13:42:41 -0000	1.28
***************
*** 220,223 ****
--- 220,226 ----
    miteShadeSpec *shpec;       /* information based on muu->shadeStr */
    double time0;               /* rendering start time */
+   gageQuery queryMite;        /* record of the miteVal quantities which 
+ 				 we'll need to compute per-sample */
+   int queryMiteNonzero;       /* shortcut miteVal computation if possible */
  
    /* as long as there's no mutex around how the miteThreads are
***************
*** 318,321 ****
--- 321,325 ----
    gageContext *gctx;            /* per-thread context */
    gage_t *ansScl,               /* pointer to gageKindScl answer vector */
+     *nPerp, *geomTens,          /* convenience pointers into ansScl */
      *ansVec,                    /* pointer to gageKindVec answer vector */
      *ansTen,                    /* pointer to tenGageKind answer vector */
***************
*** 323,326 ****
--- 327,331 ----
  				   unlike ans{Scl,Vec,Ten} is allocated by
  				   mite instead of by gage */
+     **directAnsMiteVal,         /* pointers into ansMiteVal */
      *shadeVec0, *shadeVec1, 
      *shadeScl0, *shadeScl1;     /* pointers into the ans* arrays above,
***************
*** 359,365 ****
  TEEM_API void miteVariablePrint(char *buff, const gageItemSpec *isp);
  TEEM_API int miteNtxfCheck(const Nrrd *ntxf);
! TEEM_API void miteQueryAdd(gageQuery queryScl, 
! 			   gageQuery queryVec, 
! 			   gageQuery queryTen, gageItemSpec *isp);
  
  
--- 364,370 ----
  TEEM_API void miteVariablePrint(char *buff, const gageItemSpec *isp);
  TEEM_API int miteNtxfCheck(const Nrrd *ntxf);
! TEEM_API void miteQueryAdd(gageQuery queryScl, gageQuery queryVec, 
! 			   gageQuery queryTen, gageQuery queryMite,
! 			   gageItemSpec *isp);
  
  
***************
*** 373,379 ****
  TEEM_API int miteShadeSpecParse(miteShadeSpec *shpec, char *shadeStr);
  TEEM_API void miteShadeSpecPrint(char *buff, const miteShadeSpec *shpec);
! TEEM_API void miteShadeSpecQueryAdd(gageQuery queryScl, 
! 				    gageQuery queryVec, 
! 				    gageQuery queryTen, miteShadeSpec *shpec);
  
  /* renderMite.c */
--- 378,384 ----
  TEEM_API int miteShadeSpecParse(miteShadeSpec *shpec, char *shadeStr);
  TEEM_API void miteShadeSpecPrint(char *buff, const miteShadeSpec *shpec);
! TEEM_API void miteShadeSpecQueryAdd(gageQuery queryScl, gageQuery queryVec, 
! 				    gageQuery queryTen, gageQuery queryMite,
! 				    miteShadeSpec *shpec);
  
  /* renderMite.c */

Index: ray.c
===================================================================
RCS file: /cvsroot/teem/teem/src/mite/ray.c,v
retrieving revision 1.18
retrieving revision 1.19
diff -C2 -d -r1.18 -r1.19
*** ray.c	26 Feb 2004 06:26:10 -0000	1.18
--- ray.c	29 Feb 2004 13:42:41 -0000	1.19
***************
*** 138,146 ****
    }
  
    if (1-mtt->TT >= muu->opacNear1) {
-     /* early ray termination */
      mtt->TT = 0.0;
      return 0.0;
    }
    mtt->samples += 1;
    if (gageProbe(mtt->gctx,
--- 138,149 ----
    }
  
+   /* early ray termination */
    if (1-mtt->TT >= muu->opacNear1) {
      mtt->TT = 0.0;
      return 0.0;
    }
+ 
+   /* do probing at this location to determine values of everything
+      that might appear in the txf domain */
    mtt->samples += 1;
    if (gageProbe(mtt->gctx,
***************
*** 149,186 ****
      biffAdd(MITE, err); return AIR_NAN;
    }
!   /* HEY: NONE of this should be done if the txfs don't need any miteVal */
!   mtt->ansMiteVal[miteValXw] = samplePosWorld[0];
!   mtt->ansMiteVal[miteValXi] = samplePosIndex[0];
!   mtt->ansMiteVal[miteValYw] = samplePosWorld[1];
!   mtt->ansMiteVal[miteValYi] = samplePosIndex[1];
!   mtt->ansMiteVal[miteValZw] = samplePosWorld[2];
!   mtt->ansMiteVal[miteValZi] = samplePosIndex[2];
!   mtt->ansMiteVal[miteValTw] = rayT;
!   mtt->ansMiteVal[miteValTi] = num;
! #if 0
  
!   what if shading is none?
  
!   mtt->ansMiteVal[miteValNdotV] =
!     -muu->normalSide*ELL_3V_DOT(mtt->V, mtt->shadeVec0);
!   mtt->ansMiteVal[miteValNdotL] =
!     -muu->normalSide*ELL_3V_DOT(mtt->shadeVec0, muu->lit->dir[0]);
! #endif
!   if (!muu->normalSide) {
!     mtt->ansMiteVal[miteValNdotV] = AIR_ABS(mtt->ansMiteVal[miteValNdotV]);
!     mtt->ansMiteVal[miteValNdotL] = AIR_ABS(mtt->ansMiteVal[miteValNdotL]);
    }
- 
- #if 0
-   this is bullshit
- 
-   ELL_3MV_MUL(kn, mtt->nPerp, mtt->V);
-   ELL_3V_NORM(kn, kn, len);
-   ELL_3MV_MUL(knd, mtt->gten, kn);
-   mtt->ansMiteVal[miteValGTdotV] = ELL_3V_DOT(knd, kn);
- #endif
    
    memcpy(mtt->range, muu->rangeInit, MITE_RANGE_NUM*sizeof(mite_t));
    _miteStageRun(mtt);
    if (mtt->range[miteRangeAlpha]) {
      /* fprintf(stderr, "%s: mtt->TT = %g\n", me, mtt->TT); */
--- 152,197 ----
      biffAdd(MITE, err); return AIR_NAN;
    }
!   if (mrr->queryMiteNonzero) {
!     /* There is some optimal trade-off between slowing things down
!        with too many branches on all possible checks of queryMite,
!        and slowing things down with doing the work of setting them all.
!        This code has not been profiled whatsoever */
!     mtt->directAnsMiteVal[miteValXw][0] = samplePosWorld[0];
!     mtt->directAnsMiteVal[miteValXi][0] = samplePosIndex[0];
!     mtt->directAnsMiteVal[miteValYw][0] = samplePosWorld[1];
!     mtt->directAnsMiteVal[miteValYi][0] = samplePosIndex[1];
!     mtt->directAnsMiteVal[miteValZw][0] = samplePosWorld[2];
!     mtt->directAnsMiteVal[miteValZi][0] = samplePosIndex[2];
!     mtt->directAnsMiteVal[miteValTw][0] = rayT;
!     mtt->directAnsMiteVal[miteValTi][0] = num;
  
!     if (mtt->shadeVec0 && 
! 	(GAGE_QUERY_ITEM_TEST(mrr->queryMite, miteValNdotV)
! 	 || GAGE_QUERY_ITEM_TEST(mrr->queryMite, miteValNdotL))) {
!       mtt->directAnsMiteVal[miteValNdotV][0] =
! 	-muu->normalSide*ELL_3V_DOT(mtt->shadeVec0, mtt->V);
!       mtt->directAnsMiteVal[miteValNdotL][0] =
! 	-muu->normalSide*ELL_3V_DOT(mtt->shadeVec0, muu->lit->dir[0]);
!       if (!muu->normalSide) {
! 	mtt->directAnsMiteVal[miteValNdotV][0] = 
! 	  AIR_ABS(mtt->directAnsMiteVal[miteValNdotV][0]);
! 	mtt->directAnsMiteVal[miteValNdotL][0] = 
! 	  AIR_ABS(mtt->directAnsMiteVal[miteValNdotL][0]);
!       }
!     }
  
!     if (GAGE_QUERY_ITEM_TEST(mrr->queryMite, miteValGTdotV)) {
!       ELL_3MV_MUL(kn, mtt->nPerp, mtt->V);
!       ELL_3V_NORM(kn, kn, len);
!       ELL_3MV_MUL(knd, mtt->geomTens, kn);
!       mtt->ansMiteVal[miteValGTdotV] = ELL_3V_DOT(knd, kn);
!     }
    }
    
+   /* initialize txf range quantities, and apply all txfs */
    memcpy(mtt->range, muu->rangeInit, MITE_RANGE_NUM*sizeof(mite_t));
    _miteStageRun(mtt);
+ 
+   /* if there's opacity, do shading and compositing */
    if (mtt->range[miteRangeAlpha]) {
      /* fprintf(stderr, "%s: mtt->TT = %g\n", me, mtt->TT); */

Index: renderMite.c
===================================================================
RCS file: /cvsroot/teem/teem/src/mite/renderMite.c,v
retrieving revision 1.11
retrieving revision 1.12
diff -C2 -d -r1.11 -r1.12
*** renderMite.c	26 Feb 2004 06:26:10 -0000	1.11
--- renderMite.c	29 Feb 2004 13:42:41 -0000	1.12
***************
*** 34,43 ****
      mrr->ntxf = NULL;
      mrr->ntxfNum = 0;
-     mrr->shpec = miteShadeSpecNew();
      mrr->sclPvlIdx = -1;
      mrr->vecPvlIdx = -1;
      mrr->tenPvlIdx = -1;
      airMopAdd(mrr->rmop, mrr->shpec,
  	      (airMopper)miteShadeSpecNix, airMopAlways);
    }
    return mrr;
--- 34,46 ----
      mrr->ntxf = NULL;
      mrr->ntxfNum = 0;
      mrr->sclPvlIdx = -1;
      mrr->vecPvlIdx = -1;
      mrr->tenPvlIdx = -1;
+     mrr->shpec = miteShadeSpecNew();
      airMopAdd(mrr->rmop, mrr->shpec,
  	      (airMopper)miteShadeSpecNix, airMopAlways);
+     mrr->time0 = AIR_NAN;
+     GAGE_QUERY_RESET(mrr->queryMite);
+     mrr->queryMiteNonzero = AIR_FALSE;
    }
    return mrr;
***************
*** 82,93 ****
    GAGE_QUERY_RESET(queryVec);
    GAGE_QUERY_RESET(queryTen);
    for (T=0; T<muu->ntxfNum; T++) {
      for (axi=1; axi<muu->ntxf[T]->dim; axi++) {
        miteVariableParse(&isp, muu->ntxf[T]->axis[axi].label);
!       miteQueryAdd(queryScl, queryVec, queryTen, &isp);
      }
    }
    miteShadeSpecParse((*mrrP)->shpec, muu->shadeStr);
!   miteShadeSpecQueryAdd(queryScl, queryVec, queryTen, (*mrrP)->shpec);
  
    E = 0;
--- 85,99 ----
    GAGE_QUERY_RESET(queryVec);
    GAGE_QUERY_RESET(queryTen);
+   GAGE_QUERY_RESET((*mrrP)->queryMite);
    for (T=0; T<muu->ntxfNum; T++) {
      for (axi=1; axi<muu->ntxf[T]->dim; axi++) {
        miteVariableParse(&isp, muu->ntxf[T]->axis[axi].label);
!       miteQueryAdd(queryScl, queryVec, queryTen, (*mrrP)->queryMite, &isp);
      }
    }
    miteShadeSpecParse((*mrrP)->shpec, muu->shadeStr);
!   miteShadeSpecQueryAdd(queryScl, queryVec, queryTen, (*mrrP)->queryMite,
! 			(*mrrP)->shpec);
!   (*mrrP)->queryMiteNonzero = GAGE_QUERY_NONZERO((*mrrP)->queryMite);
  
    E = 0;

Index: shade.c
===================================================================
RCS file: /cvsroot/teem/teem/src/mite/shade.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -C2 -d -r1.1 -r1.2
*** shade.c	26 Feb 2004 06:41:18 -0000	1.1
--- shade.c	29 Feb 2004 13:42:41 -0000	1.2
***************
*** 193,199 ****
  
  void
! miteShadeSpecQueryAdd(gageQuery queryScl, 
! 		      gageQuery queryVec, 
! 		      gageQuery queryTen, miteShadeSpec *shpec) {
    if (shpec) {
      switch(shpec->method) {
--- 193,199 ----
  
  void
! miteShadeSpecQueryAdd(gageQuery queryScl, gageQuery queryVec, 
! 		      gageQuery queryTen, gageQuery queryMite,
! 		      miteShadeSpec *shpec) {
    if (shpec) {
      switch(shpec->method) {
***************
*** 202,212 ****
        break;
      case miteShadeMethodPhong:
!       miteQueryAdd(queryScl, queryVec, queryTen, shpec->vec0);
        break;
      case miteShadeMethodLitTen:
!       miteQueryAdd(queryScl, queryVec, queryTen, shpec->vec0);
!       miteQueryAdd(queryScl, queryVec, queryTen, shpec->vec1);
!       miteQueryAdd(queryScl, queryVec, queryTen, shpec->scl0);
!       miteQueryAdd(queryScl, queryVec, queryTen, shpec->scl1);
        break;
      default:
--- 202,212 ----
        break;
      case miteShadeMethodPhong:
!       miteQueryAdd(queryScl, queryVec, queryTen, queryMite, shpec->vec0);
        break;
      case miteShadeMethodLitTen:
!       miteQueryAdd(queryScl, queryVec, queryTen, queryMite, shpec->vec0);
!       miteQueryAdd(queryScl, queryVec, queryTen, queryMite, shpec->vec1);
!       miteQueryAdd(queryScl, queryVec, queryTen, queryMite, shpec->scl0);
!       miteQueryAdd(queryScl, queryVec, queryTen, queryMite, shpec->scl1);
        break;
      default:

Index: thread.c
===================================================================
RCS file: /cvsroot/teem/teem/src/mite/thread.c,v
retrieving revision 1.10
retrieving revision 1.11
diff -C2 -d -r1.10 -r1.11
*** thread.c	26 Feb 2004 06:26:10 -0000	1.10
--- thread.c	29 Feb 2004 13:42:41 -0000	1.11
***************
*** 23,40 ****
  miteThread *
  miteThreadNew() {
    miteThread *mtt;
    
    mtt = (miteThread *)calloc(1, sizeof(miteThread));
    if (!mtt) {
!     return NULL;
    }
  
    mtt->gctx = NULL;
    mtt->ansScl = mtt->ansVec = mtt->ansTen = NULL;
!   mtt->shadeVec0 = mtt->shadeVec1 = NULL;
!   mtt->shadeScl0 = mtt->shadeScl1 = NULL;
    mtt->ansMiteVal = 
      (gage_t *)calloc(gageKindTotalAnswerLength(miteValGageKind), 
  		     sizeof(gage_t));
    mtt->verbose = 0;
    mtt->thrid = -1;
--- 23,57 ----
  miteThread *
  miteThreadNew() {
+   char me[]="miteThreadNew", err[AIR_STRLEN_MED];
    miteThread *mtt;
+   int ii;
    
    mtt = (miteThread *)calloc(1, sizeof(miteThread));
    if (!mtt) {
!     sprintf(err, "%s: couldn't calloc miteThread", me);
!     biffAdd(MITE, err); return NULL;
    }
  
    mtt->gctx = NULL;
    mtt->ansScl = mtt->ansVec = mtt->ansTen = NULL;
!   mtt->shadeVec0 = NULL;
!   mtt->shadeVec1 = NULL;
!   mtt->shadeScl0 = NULL;
!   mtt->shadeScl1 = NULL;
!   /* were miteVal a full-fledged gageKind, the following would
!      be done by gagePerVolumeNew */
    mtt->ansMiteVal = 
      (gage_t *)calloc(gageKindTotalAnswerLength(miteValGageKind), 
  		     sizeof(gage_t));
+   mtt->directAnsMiteVal = 
+     (gage_t **)calloc(miteValGageKind->itemMax+1, sizeof(gage_t*));
+   if (!(mtt->ansMiteVal && mtt->directAnsMiteVal)) {
+     sprintf(err, "%s: couldn't calloc miteVal answer arrays", me);
+     biffAdd(MITE, err); return NULL;
+   }
+   for (ii=0; ii<=miteValGageKind->itemMax; ii++) {
+     mtt->directAnsMiteVal[ii] = mtt->ansMiteVal 
+       + gageKindAnswerOffset(miteValGageKind, ii);
+   }
    mtt->verbose = 0;
    mtt->thrid = -1;
***************
*** 81,87 ****
    }
  
!   (*mttP)->ansScl = (-1 != mrr->sclPvlIdx
! 		     ? (*mttP)->gctx->pvl[mrr->sclPvlIdx]->answer
! 		     : NULL);
    (*mttP)->ansVec = (-1 != mrr->vecPvlIdx
  		     ? (*mttP)->gctx->pvl[mrr->vecPvlIdx]->answer
--- 98,112 ----
    }
  
!   if (-1 != mrr->sclPvlIdx) {
!     (*mttP)->ansScl = (*mttP)->gctx->pvl[mrr->sclPvlIdx]->answer;
!     (*mttP)->nPerp = ((*mttP)->ansScl 
! 		      + gageKindAnswerOffset(gageKindScl, gageSclNPerp));
!     (*mttP)->geomTens = ((*mttP)->ansScl
! 			 + gageKindAnswerOffset(gageKindScl, gageSclGeomTens));
!   } else {
!     (*mttP)->ansScl = NULL;
!     (*mttP)->nPerp = NULL;
!     (*mttP)->geomTens = NULL;
!   }
    (*mttP)->ansVec = (-1 != mrr->vecPvlIdx
  		     ? (*mttP)->gctx->pvl[mrr->vecPvlIdx]->answer
***************
*** 90,99 ****
  		     ? (*mttP)->gctx->pvl[mrr->tenPvlIdx]->answer
  		     : NULL);
- 
- #if 0
-   (*mttP)->nPerp = (*mttP)->ans + gageKindScl->ansOffset[gageSclNPerp];
-   (*mttP)->gten = (*mttP)->ans + gageKindScl->ansOffset[gageSclGeomTens];
- #endif
- 
    (*mttP)->thrid = whichThread;
    (*mttP)->samples = 0;
--- 115,118 ----

Index: txf.c
===================================================================
RCS file: /cvsroot/teem/teem/src/mite/txf.c,v
retrieving revision 1.17
retrieving revision 1.18
diff -C2 -d -r1.17 -r1.18
*** txf.c	28 Feb 2004 15:11:05 -0000	1.17
--- txf.c	29 Feb 2004 13:42:42 -0000	1.18
***************
*** 270,278 ****
  ** up explicitly depending on the libraries supplying those gageKinds
  ** (like how mite now must depend on ten)
  */
  void
! miteQueryAdd(gageQuery queryScl, 
! 	     gageQuery queryVec, 
! 	     gageQuery queryTen, gageItemSpec *isp) {
    char me[]="miteQueryAdd";
    
--- 270,282 ----
  ** up explicitly depending on the libraries supplying those gageKinds
  ** (like how mite now must depend on ten)
+ **
+ ** The queryMite argument is a little odd- its not a real gage kind,
+ ** but we use it to organize which of the miteVal quantities we take
+ ** the time to compute in miteSample().
  */
  void
! miteQueryAdd(gageQuery queryScl, gageQuery queryVec, 
! 	     gageQuery queryTen, gageQuery queryMite,
! 	     gageItemSpec *isp) {
    char me[]="miteQueryAdd";
    
***************
*** 284,287 ****
--- 288,295 ----
      GAGE_QUERY_ITEM_ON(queryTen, isp->item);
    } else if (miteValGageKind == isp->kind) {
+     /* regardless of whether the mite query requires scl, vec, or ten
+        queries, we add it to the quantites that have to be computed
+        per-sample */
+     GAGE_QUERY_ITEM_ON(queryMite, isp->item);
      /* HEY: some these have useful analogs for tensor data, but I
         won't be able to express them.  This means that while Phong

Index: user.c
===================================================================
RCS file: /cvsroot/teem/teem/src/mite/user.c,v
retrieving revision 1.15
retrieving revision 1.16
diff -C2 -d -r1.15 -r1.16
*** user.c	26 Feb 2004 06:26:10 -0000	1.15
--- user.c	29 Feb 2004 13:42:42 -0000	1.16
***************
*** 78,82 ****
    int T, axi, gotOpac;
    gageItemSpec isp;
!   gageQuery queryScl, queryVec, queryTen;
    miteShadeSpec *shpec;
    airArray *mop;
--- 78,82 ----
    int T, axi, gotOpac;
    gageItemSpec isp;
!   gageQuery queryScl, queryVec, queryTen, queryMite;
    miteShadeSpec *shpec;
    airArray *mop;
***************
*** 95,98 ****
--- 95,99 ----
    GAGE_QUERY_RESET(queryVec);
    GAGE_QUERY_RESET(queryTen);
+   GAGE_QUERY_RESET(queryMite);  /* not actually used here */
    for (T=0; T<muu->ntxfNum; T++) {
      if (miteNtxfCheck(muu->ntxf[T])) {
***************
*** 104,108 ****
      for (axi=1; axi<muu->ntxf[T]->dim; axi++) {
        miteVariableParse(&isp, muu->ntxf[T]->axis[axi].label);
!       miteQueryAdd(queryScl, queryVec, queryTen, &isp);
      }
      gotOpac |= !!strchr(muu->ntxf[T]->axis[0].label, 'A');
--- 105,109 ----
      for (axi=1; axi<muu->ntxf[T]->dim; axi++) {
        miteVariableParse(&isp, muu->ntxf[T]->axis[axi].label);
!       miteQueryAdd(queryScl, queryVec, queryTen, queryMite, &isp);
      }
      gotOpac |= !!strchr(muu->ntxf[T]->axis[0].label, 'A');
***************
*** 123,127 ****
      biffAdd(MITE, err); airMopError(mop); return 1;
    }
!   miteShadeSpecQueryAdd(queryScl, queryVec, queryTen, shpec);
    if (GAGE_QUERY_NONZERO(queryScl) && !(muu->nsin)) {
      sprintf(err, "%s: txf or shading require %s volume, but don't have one",
--- 124,128 ----
      biffAdd(MITE, err); airMopError(mop); return 1;
    }
!   miteShadeSpecQueryAdd(queryScl, queryVec, queryTen, queryMite, shpec);
    if (GAGE_QUERY_NONZERO(queryScl) && !(muu->nsin)) {
      sprintf(err, "%s: txf or shading require %s volume, but don't have one",

[Teem-commit] teem/src/gage gage.h,1.61,1.62 scl.c,1.35,1.36 sclanswer.c,1.18,1.19

[<kin...@us...>]

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

Modified Files:
	gage.h scl.c sclanswer.c 
Log Message:
added frobenious norm of hessian ("hf") to gageSclKind

Index: gage.h
===================================================================
RCS file: /cvsroot/teem/teem/src/gage/gage.h,v
retrieving revision 1.61
retrieving revision 1.62
diff -C2 -d -r1.61 -r1.62
*** gage.h	26 Feb 2004 06:19:28 -0000	1.61
--- gage.h	28 Feb 2004 17:01:53 -0000	1.62
***************
*** 242,269 ****
    gageSclHessian,     /*  5: "h", Hessian: GT[9] (column-order) */
    gageSclLaplacian,   /*  6: "l", Laplacian: Dxx + Dyy + Dzz: GT[1] */
!   gageSclHessEval,    /*  7: "heval", Hessian's eigenvalues: GT[3] */
!   gageSclHessEval0,   /*  8: "heval0", Hessian's 1st eigenvalue: GT[1] */
!   gageSclHessEval1,   /*  9: "heval1", Hessian's 2nd eigenvalue: GT[1] */
!   gageSclHessEval2,   /* 10: "heval2", Hessian's 3rd eigenvalue: GT[1] */
!   gageSclHessEvec,    /* 11: "hevec", Hessian's eigenvectors: GT[9] */
!   gageSclHessEvec0,   /* 12: "hevec0", Hessian's 1st eigenvector: GT[3] */
!   gageSclHessEvec1,   /* 13: "hevec1", Hessian's 2nd eigenvector: GT[3] */
!   gageSclHessEvec2,   /* 14: "hevec2", Hessian's 3rd eigenvector: GT[3] */
!   gageScl2ndDD,       /* 15: "2d", 2nd dir.deriv. along gradient: GT[1] */
!   gageSclGeomTens,    /* 16: "gten", sym. matx w/ evals {0, K1, K2} and
  			     evecs {grad, cdir0, cdir1}: GT[9] */
!   gageSclK1,          /* 17: "k1", 1st principle curvature: GT[1] */
!   gageSclK2,          /* 18: "k2", 2nd principle curvature (k2 <= k1): GT[1] */
!   gageSclTotalCurv,   /* 19: "tc", L2 norm(K1,K2) (not Koen.'s "C"): GT[1] */
!   gageSclShapeTrace,  /* 20, "st", (K1+K2)/Curvedness: GT[1] */
!   gageSclShapeIndex,  /* 21: "si", Koen.'s shape index, ("S"): GT[1] */
!   gageSclMeanCurv,    /* 22: "mc", mean curvature (K1 + K2)/2: GT[1] */
!   gageSclGaussCurv,   /* 23: "gc", gaussian curvature K1*K2: GT[1] */
!   gageSclCurvDir1,    /* 24: "cdir1", 1st principle curv direction: GT[3] */
!   gageSclCurvDir2,    /* 25: "cdir2", 2nd principle curv direction: GT[3] */
!   gageSclFlowlineCurv,/* 26: "fc", curvature of normal streamline: GT[1] */
    gageSclLast
  };
! #define GAGE_SCL_ITEM_MAX     26
  
  /*
--- 242,270 ----
    gageSclHessian,     /*  5: "h", Hessian: GT[9] (column-order) */
    gageSclLaplacian,   /*  6: "l", Laplacian: Dxx + Dyy + Dzz: GT[1] */
!   gageSclHessFrob,    /*  7: "hf", Frobenius normal of Hessian: GT[1] */
!   gageSclHessEval,    /*  8: "heval", Hessian's eigenvalues: GT[3] */
!   gageSclHessEval0,   /*  9: "heval0", Hessian's 1st eigenvalue: GT[1] */
!   gageSclHessEval1,   /* 10: "heval1", Hessian's 2nd eigenvalue: GT[1] */
!   gageSclHessEval2,   /* 11: "heval2", Hessian's 3rd eigenvalue: GT[1] */
!   gageSclHessEvec,    /* 12: "hevec", Hessian's eigenvectors: GT[9] */
!   gageSclHessEvec0,   /* 13: "hevec0", Hessian's 1st eigenvector: GT[3] */
!   gageSclHessEvec1,   /* 14: "hevec1", Hessian's 2nd eigenvector: GT[3] */
!   gageSclHessEvec2,   /* 15: "hevec2", Hessian's 3rd eigenvector: GT[3] */
!   gageScl2ndDD,       /* 16: "2d", 2nd dir.deriv. along gradient: GT[1] */
!   gageSclGeomTens,    /* 17: "gten", sym. matx w/ evals {0, K1, K2} and
  			     evecs {grad, cdir0, cdir1}: GT[9] */
!   gageSclK1,          /* 18: "k1", 1st principle curvature: GT[1] */
!   gageSclK2,          /* 19: "k2", 2nd principle curvature (k2 <= k1): GT[1] */
!   gageSclTotalCurv,   /* 20: "tc", L2 norm(K1,K2) (not Koen.'s "C"): GT[1] */
!   gageSclShapeTrace,  /* 21, "st", (K1+K2)/Curvedness: GT[1] */
!   gageSclShapeIndex,  /* 22: "si", Koen.'s shape index, ("S"): GT[1] */
!   gageSclMeanCurv,    /* 23: "mc", mean curvature (K1 + K2)/2: GT[1] */
!   gageSclGaussCurv,   /* 24: "gc", gaussian curvature K1*K2: GT[1] */
!   gageSclCurvDir1,    /* 25: "cdir1", 1st principle curv direction: GT[3] */
!   gageSclCurvDir2,    /* 26: "cdir2", 2nd principle curv direction: GT[3] */
!   gageSclFlowlineCurv,/* 27: "fc", curvature of normal streamline: GT[1] */
    gageSclLast
  };
! #define GAGE_SCL_ITEM_MAX     27
  
  /*

Index: scl.c
===================================================================
RCS file: /cvsroot/teem/teem/src/gage/scl.c,v
retrieving revision 1.35
retrieving revision 1.36
diff -C2 -d -r1.35 -r1.36
*** scl.c	28 Feb 2004 15:24:01 -0000	1.35
--- scl.c	28 Feb 2004 17:01:53 -0000	1.36
***************
*** 42,45 ****
--- 42,46 ----
    {gageSclHessian,       9,  2,  {gageSclHessian, -1, -1, -1, -1},                        -1,  -1},
    {gageSclLaplacian,     1,  2,  {gageSclHessian, -1, -1, -1, -1},                        -1,  -1},
+   {gageSclHessFrob,      1,  2,  {gageSclHessian, -1, -1, -1, -1},                        -1,  -1},
    {gageSclHessEval,      3,  2,  {gageSclHessian, -1, -1, -1, -1},                        -1,  -1},
    {gageSclHessEval0,     1,  2,  {gageSclHessEval, -1, -1, -1, -1},                       gageSclHessEval,  0},
***************
*** 74,77 ****
--- 75,79 ----
    "Hessian",
    "Laplacian",
+   "Frob(Hessian)",
    "Hessian eigenvalues",
    "Hessian eigenvalue[0]",
***************
*** 106,109 ****
--- 108,112 ----
    "3x3 Hessian matrix",
    "Laplacian",
+   "Frobenius norm of Hessian",
    "Hessian's eigenvalues",
    "Hessian's 1st eigenvalue",
***************
*** 138,141 ****
--- 141,145 ----
    gageSclHessian,
    gageSclLaplacian,
+   gageSclHessFrob,
    gageSclHessEval,
    gageSclHessEval0,
***************
*** 167,170 ****
--- 171,175 ----
  #define GS_H   gageSclHessian
  #define GS_L   gageSclLaplacian
+ #define GS_HF  gageSclHessFrob
  #define GS_HA  gageSclHessEval
  #define GS_HA0 gageSclHessEval0
***************
*** 198,201 ****
--- 203,207 ----
    "h", "hess", "hessian",
    "l", "lapl", "laplacian",
+   "hf",
    "heval", "h eval", "hessian eval", "hessian eigenvalues",
    "heval0",
***************
*** 230,233 ****
--- 236,240 ----
    GS_H, GS_H, GS_H, 
    GS_L, GS_L, GS_L, 
+   GS_HF,
    GS_HA, GS_HA, GS_HA, GS_HA, 
    GS_HA0,

Index: sclanswer.c
===================================================================
RCS file: /cvsroot/teem/teem/src/gage/sclanswer.c,v
retrieving revision 1.18
retrieving revision 1.19
diff -C2 -d -r1.18 -r1.19
*** sclanswer.c	19 Feb 2004 06:45:12 -0000	1.18
--- sclanswer.c	28 Feb 2004 17:01:53 -0000	1.19
***************
*** 95,98 ****
--- 95,101 ----
      }
    }
+   if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessFrob)) {
+     pvl->directAnswer[gageSclHessFrob][0] = ELL_3M_FROB(hess);
+   }
    if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessEval)) {
      ELL_3M_COPY(tmpMat, hess);

[Teem-commit] teem/src/gage scl.c,1.34,1.35

[<kin...@us...>]

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

Modified Files:
	scl.c 
Log Message:
fixed HUGE bug in computation of geometry tensor

Index: scl.c
===================================================================
RCS file: /cvsroot/teem/teem/src/gage/scl.c,v
retrieving revision 1.34
retrieving revision 1.35
diff -C2 -d -r1.34 -r1.35
*** scl.c	28 Feb 2004 15:09:38 -0000	1.34
--- scl.c	28 Feb 2004 15:24:01 -0000	1.35
***************
*** 51,55 ****
    {gageSclHessEvec2,     3,  2,  {gageSclHessEvec, -1, -1, -1, -1},                       gageSclHessEvec,  6},
    {gageScl2ndDD,         1,  2,  {gageSclHessian, gageSclNormal, -1, -1, -1},             -1,  -1},
!   {gageSclGeomTens,      1,  2,  {gageSclHessian, gageSclNPerp, gageSclGradMag, -1, -1},  -1,  -1},
    {gageSclK1,            1,  2,  {gageSclTotalCurv, gageSclShapeTrace, -1, -1, -1},      -1,  -1},
    {gageSclK2,            1,  2,  {gageSclTotalCurv, gageSclShapeTrace, -1, -1, -1},      -1,  -1},
--- 51,55 ----
    {gageSclHessEvec2,     3,  2,  {gageSclHessEvec, -1, -1, -1, -1},                       gageSclHessEvec,  6},
    {gageScl2ndDD,         1,  2,  {gageSclHessian, gageSclNormal, -1, -1, -1},             -1,  -1},
!   {gageSclGeomTens,      9,  2,  {gageSclHessian, gageSclNPerp, gageSclGradMag, -1, -1},  -1,  -1},
    {gageSclK1,            1,  2,  {gageSclTotalCurv, gageSclShapeTrace, -1, -1, -1},      -1,  -1},
    {gageSclK2,            1,  2,  {gageSclTotalCurv, gageSclShapeTrace, -1, -1, -1},      -1,  -1},
***************
*** 207,211 ****
    "hevec2",
    "2d", "2dd", "2nddd", "2nd", "2nd dd", "2nd dd along gradient",
!   "gten", "geoten", "geomten", "geometry tensor",
    "k1", "kap1", "kappa1",
    "k2", "kap2", "kappa2",
--- 207,211 ----
    "hevec2",
    "2d", "2dd", "2nddd", "2nd", "2nd dd", "2nd dd along gradient",
!   "gt", "gten", "geoten", "geomten", "geometry tensor",
    "k1", "kap1", "kappa1",
    "k2", "kap2", "kappa2",
***************
*** 239,243 ****
    GS_HE2,
    GS_2D, GS_2D, GS_2D, GS_2D, GS_2D, GS_2D,
!   GS_GT, GS_GT, GS_GT, GS_GT, 
    GS_K1, GS_K1, GS_K1,
    GS_K2, GS_K2, GS_K2,
--- 239,243 ----
    GS_HE2,
    GS_2D, GS_2D, GS_2D, GS_2D, GS_2D, GS_2D,
!   GS_GT, GS_GT, GS_GT, GS_GT, GS_GT, 
    GS_K1, GS_K1, GS_K1,
    GS_K2, GS_K2, GS_K2,

[Teem-commit] teem/src/mite txf.c,1.16,1.17

[<kin...@us...>]

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

Modified Files:
	txf.c 
Log Message:
can now use all the new quantized normal methods

Index: txf.c
===================================================================
RCS file: /cvsroot/teem/teem/src/mite/txf.c,v
retrieving revision 1.16
retrieving revision 1.17
diff -C2 -d -r1.16 -r1.17
*** txf.c	27 Feb 2004 16:43:04 -0000	1.16
--- txf.c	28 Feb 2004 15:11:05 -0000	1.17
***************
*** 163,167 ****
    char me[]="miteNtxfCheck", err[AIR_STRLEN_MED], *rangeStr, *domStr;
    gageItemSpec isp;
!   int rii, axi;
  
    if (nrrdCheck(ntxf)) {
--- 163,167 ----
    char me[]="miteNtxfCheck", err[AIR_STRLEN_MED], *rangeStr, *domStr;
    gageItemSpec isp;
!   int rii, axi, log2;
  
    if (nrrdCheck(ntxf)) {
***************
*** 226,237 ****
      if (3 == isp.kind->table[isp.item].answerLength) {
        /* has to be right length for one of the quantization schemes */
!       if (!( limnQNBins[limnQN16checker] == ntxf->axis[axi].size ||
! 	     limnQNBins[limnQN14checker] == ntxf->axis[axi].size ||
! 	     limnQNBins[limnQN12checker] == ntxf->axis[axi].size )) {
! 	sprintf(err, "%s: vector %s can be quantized into %d, %d, or %d bins "
! 		"but not %d", me, domStr, limnQNBins[limnQN16checker],
! 		limnQNBins[limnQN14checker], limnQNBins[limnQN12checker],
! 		ntxf->axis[axi].size);
  	biffAdd(MITE, err); return 1;
        }
      } else {
--- 226,240 ----
      if (3 == isp.kind->table[isp.item].answerLength) {
        /* has to be right length for one of the quantization schemes */
!       log2 = airLog2(ntxf->axis[axi].size);
!       if (-1 == log2) {
! 	sprintf(err, "%s: txf axis size for %s must be power of 2 (not %d)",
! 		me, domStr, ntxf->axis[axi].size);
  	biffAdd(MITE, err); return 1;
+       } else {
+ 	if (!( AIR_IN_CL(8, log2, 16) )) {
+ 	  sprintf(err, "%s: log_2 of txf axis size for %s should be in "
+ 		  "range [8,16] (not %d)", me, domStr, log2);
+ 	  biffAdd(MITE, err); return 1;
+ 	}
        }
      } else {
***************
*** 447,451 ****
  _miteStageSet(miteThread *mtt, miteRender *mrr) {
    char me[]="_miteStageSet", err[AIR_STRLEN_MED];
!   int ni, di, stageIdx, rii, stageNum;
    Nrrd *ntxf;
    miteStage *stage;
--- 450,454 ----
  _miteStageSet(miteThread *mtt, miteRender *mrr) {
    char me[]="_miteStageSet", err[AIR_STRLEN_MED];
!   int ni, di, stageIdx, rii, stageNum, log2;
    Nrrd *ntxf;
    miteStage *stage;
***************
*** 485,499 ****
  	  stage->qn = NULL;
  	} else if (3 == isp.kind->table[isp.item].answerLength) {
! 	  if (limnQNBins[limnQN16checker] == ntxf->axis[di].size) {
! 	    stage->qn = limnVtoQN_GT[limnQN16checker];
! 	  } else if (limnQNBins[limnQN14checker] == ntxf->axis[di].size) {
! 	    stage->qn = limnVtoQN_GT[limnQN14checker];
! 	  } else if (limnQNBins[limnQN12checker] == ntxf->axis[di].size) {
! 	    stage->qn = limnVtoQN_GT[limnQN12checker];
! 	  } else {
  	    sprintf(err, "%s: txf axis %d size %d not usable for vector "
  		    "txf domain variable %s", me,
  		    di, ntxf->axis[di].size, ntxf->axis[di].label);
  	    biffAdd(MITE, err); return 1;
  	  }
  	} else {
--- 488,508 ----
  	  stage->qn = NULL;
  	} else if (3 == isp.kind->table[isp.item].answerLength) {
! 	  log2 = airLog2(ntxf->axis[di].size);
! 	  switch(log2) {
! 	  case 8:  stage->qn = limnVtoQN_GT[ limnQN8checker]; break;
! 	  case 9:  stage->qn = limnVtoQN_GT[ limnQN9checker]; break;
! 	  case 10: stage->qn = limnVtoQN_GT[limnQN10checker]; break;
! 	  case 11: stage->qn = limnVtoQN_GT[limnQN11checker]; break;
! 	  case 12: stage->qn = limnVtoQN_GT[limnQN12checker]; break;
! 	  case 13: stage->qn = limnVtoQN_GT[limnQN13checker]; break;
! 	  case 14: stage->qn = limnVtoQN_GT[limnQN14checker]; break;
! 	  case 15: stage->qn = limnVtoQN_GT[limnQN15checker]; break;
! 	  case 16: stage->qn = limnVtoQN_GT[limnQN16checker]; break;
! 	  default:
  	    sprintf(err, "%s: txf axis %d size %d not usable for vector "
  		    "txf domain variable %s", me,
  		    di, ntxf->axis[di].size, ntxf->axis[di].label);
  	    biffAdd(MITE, err); return 1;
+ 	    break;
  	  }
  	} else {

[Teem-commit] teem/src/hest GNUmakefile,1.10,1.11

[<kin...@us...>]

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

Modified Files:
	GNUmakefile 
Log Message:
better working template for hest usage

Index: GNUmakefile
===================================================================
RCS file: /cvsroot/teem/teem/src/hest/GNUmakefile,v
retrieving revision 1.10
retrieving revision 1.11
diff -C2 -d -r1.10 -r1.11
*** GNUmakefile	7 Jan 2004 15:34:29 -0000	1.10
--- GNUmakefile	28 Feb 2004 15:10:22 -0000	1.11
***************
*** 49,53 ****
  $(L).OBJS = defaultsHest.o methodsHest.o usage.o parseHest.o
  $(L).TESTS = test/ex0 test/ex1 test/ex2 test/ex3 test/ex4 test/ex5 \
! 	test/strings test/bday
  ####
  ####
--- 49,53 ----
  $(L).OBJS = defaultsHest.o methodsHest.o usage.o parseHest.o
  $(L).TESTS = test/ex0 test/ex1 test/ex2 test/ex3 test/ex4 test/ex5 \
! 	test/strings test/bday test/tmpl
  ####
  ####

[Teem-commit] teem/src/hest/test tmpl.c,1.3,1.4

[<kin...@us...>]

Update of /cvsroot/teem/teem/src/hest/test
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv23780/hest/test

Modified Files:
	tmpl.c 
Log Message:
better working template for hest usage

Index: tmpl.c
===================================================================
RCS file: /cvsroot/teem/teem/src/hest/test/tmpl.c,v
retrieving revision 1.3
retrieving revision 1.4
diff -C2 -d -r1.3 -r1.4
*** tmpl.c	7 Jan 2004 15:34:29 -0000	1.3
--- tmpl.c	28 Feb 2004 15:10:22 -0000	1.4
***************
*** 20,23 ****
--- 20,25 ----
  #include "../hest.h"
  
+ char *tmplInfo = ("shows minimal usage of hest with mops");
+ 
  int
  main(int argc, char *argv[]) {
***************
*** 27,30 ****
--- 29,34 ----
    airArray *mop;
  
+   int size[3];
+ 
    me = argv[0];
    mop = airMopNew();
***************
*** 32,41 ****
    hopt = NULL;
    airMopAdd(mop, hparm, (airMopper)hestParmFree, airMopAlways);
!   /*
    hestOptAdd(&hopt, "s", "sx sy sz", airTypeInt, 3, 3, size, "128 128 128",
  	     "dimensions of output volume");
-   */
    hestParseOrDie(hopt, argc-1, argv+1, hparm,
! 		 me, quadInfo, AIR_TRUE, AIR_TRUE, AIR_TRUE);
    airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
    airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
--- 36,44 ----
    hopt = NULL;
    airMopAdd(mop, hparm, (airMopper)hestParmFree, airMopAlways);
! 
    hestOptAdd(&hopt, "s", "sx sy sz", airTypeInt, 3, 3, size, "128 128 128",
  	     "dimensions of output volume");
    hestParseOrDie(hopt, argc-1, argv+1, hparm,
! 		 me, tmplInfo, AIR_TRUE, AIR_TRUE, AIR_TRUE);
    airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
    airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);

[Teem-commit] teem/src/gage scl.c,1.33,1.34

[<kin...@us...>]

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

Modified Files:
	scl.c 
Log Message:
reinstated "n" as synonym for scalar gradient normalized

Index: scl.c
===================================================================
RCS file: /cvsroot/teem/teem/src/gage/scl.c,v
retrieving revision 1.33
retrieving revision 1.34
diff -C2 -d -r1.33 -r1.34
*** scl.c	26 Feb 2004 06:20:25 -0000	1.33
--- scl.c	28 Feb 2004 15:09:38 -0000	1.34
***************
*** 193,197 ****
    "gv", "gvec", "gradvec", "grad vec", "gradient vector",
    "gm", "gmag", "gradmag", "grad mag", "gradient magnitude",
!   "gn", /* pray "n", */ "normal", "gnorm", "normg", "norm", "normgrad", \
         "norm grad", "normalized gradient",
    "np", "nperp", 
--- 193,197 ----
    "gv", "gvec", "gradvec", "grad vec", "gradient vector",
    "gm", "gmag", "gradmag", "grad mag", "gradient magnitude",
!   "gn", "n", "normal", "gnorm", "normg", "norm", "normgrad", \
         "norm grad", "normalized gradient",
    "np", "nperp", 
***************
*** 226,230 ****
    GS_GV, GS_GV, GS_GV, GS_GV, GS_GV, 
    GS_GM, GS_GM, GS_GM, GS_GM, GS_GM,
!   GS_N, /* pray GS_N, */ GS_N, GS_N, GS_N, GS_N, GS_N, GS_N, GS_N,
    GS_NP, GS_NP,
    GS_H, GS_H, GS_H, 
--- 226,230 ----
    GS_GV, GS_GV, GS_GV, GS_GV, GS_GV, 
    GS_GM, GS_GM, GS_GM, GS_GM, GS_GM,
!   GS_N,  GS_N, GS_N, GS_N, GS_N, GS_N, GS_N, GS_N, GS_N,
    GS_NP, GS_NP,
    GS_H, GS_H, GS_H, 

[Teem-commit] teem/src/air miscAir.c,1.22,1.23

[<kin...@us...>]

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

Modified Files:
	miscAir.c 
Log Message:
airLog2 was broken

Index: miscAir.c
===================================================================
RCS file: /cvsroot/teem/teem/src/air/miscAir.c,v
retrieving revision 1.22
retrieving revision 1.23
diff -C2 -d -r1.22 -r1.23
*** miscAir.c	19 Feb 2004 06:45:11 -0000	1.22
--- miscAir.c	28 Feb 2004 15:08:23 -0000	1.23
***************
*** 419,423 ****
    if (n < 2)
      return -1;
!   return airLog2(n/2.0);
  }
  
--- 419,423 ----
    if (n < 2)
      return -1;
!   return 1 + airLog2(n/2.0);
  }
  

[Teem-commit] teem/src/nrrd apply1D.c,1.23,1.24

[<kin...@us...>]

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

Modified Files:
	apply1D.c 
Log Message:
added airMop usage, renamed some variables for good measure, TEAT approved

Index: apply1D.c
===================================================================
RCS file: /cvsroot/teem/teem/src/nrrd/apply1D.c,v
retrieving revision 1.23
retrieving revision 1.24
diff -C2 -d -r1.23 -r1.24
*** apply1D.c	27 Feb 2004 18:09:57 -0000	1.23
--- apply1D.c	27 Feb 2004 19:10:31 -0000	1.24
***************
*** 873,900 ****
  */
  int
! nrrdApply1DSubstitution(Nrrd *nout, const Nrrd *nin, const Nrrd *nsubst) {
    char me[]="nrrdApply1DSubstitution", err[AIR_STRLEN_MED];
    double (*lup)(const void *, size_t);
    double (*ins)(void *, size_t, double);
!   Nrrd *dsubst;
!   double val, sval, *d;
    size_t ii, num;
    int jj, asize0, asize1, changed;
  
!   if (!(nout && nsubst && nin)) {
      sprintf(err, "%s: got NULL pointer", me);
      biffAdd(NRRD, err); return 1;
    }
!   if (nrrdTypeBlock == nin->type || nrrdTypeBlock == nsubst->type) {
      sprintf(err, "%s: input or substitution type is %s, need scalar",
  	    me, airEnumStr(nrrdType, nrrdTypeBlock));
      biffAdd(NRRD, err); return 1;
    }
!   if (2 != nsubst->dim) {
      sprintf(err, "%s: substitution table has to be 2-D, not %d-D",
! 	    me, nsubst->dim);
      biffAdd(NRRD, err); return 1;
    }
!   nrrdAxisInfoGet(nsubst, nrrdAxisInfoSize, &asize0, &asize1);
    if (2 != asize0) {
      sprintf(err, "%s: substitution table has to be 2xN, not %dxN",
--- 873,901 ----
  */
  int
! nrrdApply1DSubstitution(Nrrd *nout, const Nrrd *nin, const Nrrd *_nsubst) {
    char me[]="nrrdApply1DSubstitution", err[AIR_STRLEN_MED];
    double (*lup)(const void *, size_t);
    double (*ins)(void *, size_t, double);
!   Nrrd *nsubst;
!   double val, *subst;
    size_t ii, num;
    int jj, asize0, asize1, changed;
+   airArray *mop;
  
!   if (!(nout && _nsubst && nin)) {
      sprintf(err, "%s: got NULL pointer", me);
      biffAdd(NRRD, err); return 1;
    }
!   if (nrrdTypeBlock == nin->type || nrrdTypeBlock == _nsubst->type) {
      sprintf(err, "%s: input or substitution type is %s, need scalar",
  	    me, airEnumStr(nrrdType, nrrdTypeBlock));
      biffAdd(NRRD, err); return 1;
    }
!   if (2 != _nsubst->dim) {
      sprintf(err, "%s: substitution table has to be 2-D, not %d-D",
! 	    me, _nsubst->dim);
      biffAdd(NRRD, err); return 1;
    }
!   nrrdAxisInfoGet(_nsubst, nrrdAxisInfoSize, &asize0, &asize1);
    if (2 != asize0) {
      sprintf(err, "%s: substitution table has to be 2xN, not %dxN",
***************
*** 908,926 ****
      }
    }
!   if (nrrdConvert(dsubst = nrrdNew(), nsubst, nrrdTypeDouble)) {
      sprintf(err, "%s: couldn't create double copy of substitution table", me);
!     biffAdd(NRRD, err); dsubst = nrrdNuke(dsubst); return 1;
    }
    lup = nrrdDLookup[nout->type];
    ins = nrrdDInsert[nout->type];
    num = nrrdElementNumber(nout);
-   d = (double *)dsubst->data;
    for (ii=0; ii<num; ii++) {
      val = lup(nout->data, ii);
      changed = AIR_FALSE;
      for (jj=0; jj<asize1; jj++) {
!       sval = d[jj*2+0];
!       if (val == sval) {
! 	val = d[jj*2+1];
  	changed = AIR_TRUE;
        }
--- 909,930 ----
      }
    }
! 
!   mop = airMopNew();
!   nsubst = nrrdNew();
!   airMopAdd(mop, nsubst, (airMopper)nrrdNuke, airMopAlways);
!   if (nrrdConvert(nsubst, _nsubst, nrrdTypeDouble)) {
      sprintf(err, "%s: couldn't create double copy of substitution table", me);
!     biffAdd(NRRD, err); airMopError(mop); return 1;
    }
    lup = nrrdDLookup[nout->type];
    ins = nrrdDInsert[nout->type];
+   subst = (double *)nsubst->data;
    num = nrrdElementNumber(nout);
    for (ii=0; ii<num; ii++) {
      val = lup(nout->data, ii);
      changed = AIR_FALSE;
      for (jj=0; jj<asize1; jj++) {
!       if (val == subst[jj*2+0]) {
! 	val = subst[jj*2+1];
  	changed = AIR_TRUE;
        }
***************
*** 930,934 ****
      }
    }
!   dsubst = nrrdNuke(dsubst);
    return 0;
  }
--- 934,939 ----
      }
    }
! 
!   airMopOkay(mop);
    return 0;
  }

[Teem-commit] teem/src/limn env.c,1.6,1.7 qn.c,1.18,1.19

[<kin...@us...>]

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

Modified Files:
	env.c qn.c 
Log Message:
added a whole bunch more normal quantization methods

Index: env.c
===================================================================
RCS file: /cvsroot/teem/teem/src/limn/env.c,v
retrieving revision 1.6
retrieving revision 1.7
diff -C2 -d -r1.6 -r1.7
*** env.c	27 Feb 2004 16:42:53 -0000	1.6
--- env.c	27 Feb 2004 18:16:50 -0000	1.7
***************
*** 38,98 ****
    case limnQN16checker:
      sx = sy = 256;
-     if (nrrdMaybeAlloc(map, nrrdTypeFloat, 3, 3, sx, sy)) {
-       sprintf(err, "%s: couldn't alloc output", me);
-       biffMove(LIMN, err, NRRD); return 1;
-     }
-     mapData = map->data;
-     for (qn=0; qn<=sx*sy-1; qn++) {
-       limnQNtoV_f[qnMethod](vec, qn);
-       cb(mapData + 3*qn, vec, data);
-     }
      break;
    case limnQN14checker:
      sx = sy = 128;
-     if (nrrdMaybeAlloc(map, nrrdTypeFloat, 3, 3, sx, sy)) {
-       sprintf(err, "%s: couldn't alloc output", me);
-       biffMove(LIMN, err, NRRD); return 1;
-     }
-     mapData = map->data;
-     for (qn=0; qn<=sx*sy-1; qn++) {
-       limnQNtoV_f[qnMethod](vec, qn);
-       cb(mapData + 3*qn, vec, data);
-     }
      break;
    case limnQN12checker:
      sx = sy = 64;
-     if (nrrdMaybeAlloc(map, nrrdTypeFloat, 3, 3, sx, sy)) {
-       sprintf(err, "%s: couldn't alloc output", me);
-       biffMove(LIMN, err, NRRD); return 1;
-     }
-     mapData = map->data;
-     for (qn=0; qn<=sx*sy-1; qn++) {
-       limnQNtoV_f[qnMethod](vec, qn);
-       cb(mapData + 3*qn, vec, data);
-     }
      break;
    case limnQN10checker:
      sx = sy = 32;
-     if (nrrdMaybeAlloc(map, nrrdTypeFloat, 3, 3, sx, sy)) {
-       sprintf(err, "%s: couldn't alloc output", me);
-       biffMove(LIMN, err, NRRD); return 1;
-     }
-     mapData = map->data;
-     for (qn=0; qn<=sx*sy-1; qn++) {
-       limnQNtoV_f[qnMethod](vec, qn);
-       cb(mapData + 3*qn, vec, data);
-     }
      break;
    case limnQN8checker:
      sx = sy = 16;
!     if (nrrdMaybeAlloc(map, nrrdTypeFloat, 3, 3, sx, sy)) {
!       sprintf(err, "%s: couldn't alloc output", me);
!       biffMove(LIMN, err, NRRD); return 1;
!     }
!     mapData = map->data;
!     for (qn=0; qn<=sx*sy-1; qn++) {
!       limnQNtoV_f[qnMethod](vec, qn);
!       cb(mapData + 3*qn, vec, data);
!     }
      break;
    default:
--- 38,69 ----
    case limnQN16checker:
      sx = sy = 256;
      break;
    case limnQN14checker:
      sx = sy = 128;
      break;
    case limnQN12checker:
      sx = sy = 64;
      break;
    case limnQN10checker:
      sx = sy = 32;
      break;
    case limnQN8checker:
      sx = sy = 16;
!     break;
!   case limnQN15checker:
!     sx = 128;
!     sy = 256;
!     break;
!   case limnQN13checker:
!     sx = 64;
!     sy = 128;
!     break;
!   case limnQN11checker:
!     sx = 32;
!     sy = 64;
!     break;
!   case limnQN9checker:
!     sx = 16;
!     sy = 32;
      break;
    default:
***************
*** 100,103 ****
--- 71,83 ----
      biffAdd(LIMN, err); return 1;
    }
+   if (nrrdMaybeAlloc(map, nrrdTypeFloat, 3, 3, sx, sy)) {
+     sprintf(err, "%s: couldn't alloc output", me);
+     biffMove(LIMN, err, NRRD); return 1;
+   }
+   mapData = map->data;
+   for (qn=0; qn<=sx*sy-1; qn++) {
+     limnQNtoV_f[qnMethod](vec, qn);
+     cb(mapData + 3*qn, vec, data);
+   }
  
    return 0;

Index: qn.c
===================================================================
RCS file: /cvsroot/teem/teem/src/limn/qn.c,v
retrieving revision 1.18
retrieving revision 1.19
diff -C2 -d -r1.18 -r1.19
*** qn.c	27 Feb 2004 16:43:04 -0000	1.18
--- qn.c	27 Feb 2004 18:16:50 -0000	1.19
***************
*** 219,240 ****
    (vec)[2] = z*n
  
- void
- _limnQN16checker_QNtoV_f(float *vec, int qn) {
-   int ui, vi;
-   double u, v, x, y, z, n;
- 
-   /* _16_QtoV(vec, qn); */
-   _EVEN_QtoV(8, vec, qn);
- }
- 
- void
- _limnQN16checker_QNtoV_d(double *vec, int qn) {
-   int ui, vi;
-   double u, v, x, y, z, n;
- 
-   /* _16_QtoV(vec, qn); */
-   _EVEN_QtoV(8, vec, qn);
- }
- 
  #define _16_VtoQ(vec) \
    x = (vec)[0]; \
--- 219,222 ----
***************
*** 260,263 ****
--- 242,261 ----
    }
  
+ void
+ _limnQN16checker_QNtoV_f(float *vec, int qn) {
+   int ui, vi;
+   double u, v, x, y, z, n;
+ 
+   _EVEN_QtoV(8, vec, qn);
+ }
+ 
+ void
+ _limnQN16checker_QNtoV_d(double *vec, int qn) {
+   int ui, vi;
+   double u, v, x, y, z, n;
+ 
+   _EVEN_QtoV(8, vec, qn);
+ }
+ 
  int
  _limnQN16checker_VtoQN_f(float *vec) {
***************
*** 265,272 ****
    int xi, yi, ui, vi;
    
-   /*
-   _16_VtoQ(vec);
-   return (vi << 8) | ui;
-   */
    _EVEN_VtoQ(8, vec);
  }
--- 263,266 ----
***************
*** 277,284 ****
    int xi, yi, ui, vi;
  
-   /*
-   _16_VtoQ(vec);
-   return (vi << 8) | ui;
-   */
    _EVEN_VtoQ(8, vec);
  }
--- 271,274 ----
***************
*** 286,289 ****
--- 276,314 ----
  /* ----------------------------------------------------------------  */
  
+ /* 15 bit --> HNB == 7 */
+ 
+ #define _ODD_QtoV(HNB, vec, qn) \
+   ui = qn & ((1<<HNB)-1); \
+   vi = (qn >> HNB) & ((1<<HNB)-1); \
+   zi = (qn >> (2*HNB)) & 0x01; \
+   u = AIR_AFFINE(-0.5, ui, ((1<<HNB)-1)+0.5, -0.5, 0.5); \
+   v = AIR_AFFINE(-0.5, vi, ((1<<HNB)-1)+0.5, -0.5, 0.5); \
+   x =  u + v; \
+   y =  u - v; \
+   z = 1 - AIR_ABS(x) - AIR_ABS(y); \
+   z *= (zi << 1) - 1; \
+   n = 1.0/sqrt(x*x + y*y + z*z); \
+   vec[0] = x*n; \
+   vec[1] = y*n; \
+   vec[2] = z*n
+ 
+ #define _ODD_VtoQ(HNB, vec) \
+   x = vec[0]; \
+   y = vec[1]; \
+   z = vec[2]; \
+   L = AIR_ABS(x) + AIR_ABS(y) + AIR_ABS(z); \
+   if (!L) { \
+     return 0; \
+   } \
+   x /= L; \
+   y /= L; \
+   u = x + y; \
+   v = x - y; \
+   AIR_INDEX(-1, u, 1, (1<<HNB), ui); \
+   AIR_INDEX(-1, v, 1, (1<<HNB), vi); \
+   zi = (z > 0); \
+   return (zi << (2*HNB)) | (vi << HNB) | ui
+ 
+ 
  void
  _limnQN15checker_QNtoV_f(float *vec, int qn) {
***************
*** 291,307 ****
    float u, v, x, y, z, n;
    
!   ui = qn & 0x7F;
!   vi = (qn >> 7) & 0x7F;
!   zi = (qn >> 14) & 0x01;
!   u = AIR_AFFINE(-0.5, ui, 127.5, -0.5, 0.5);
!   v = AIR_AFFINE(-0.5, vi, 127.5, -0.5, 0.5);
!   x =  u + v;
!   y =  u - v;
!   z = 1 - AIR_ABS(x) - AIR_ABS(y);
!   z *= (zi << 1) - 1;
!   n = 1.0/sqrt(x*x + y*y + z*z);
!   vec[0] = x*n; 
!   vec[1] = y*n; 
!   vec[2] = z*n;
  }
  
--- 316,328 ----
    float u, v, x, y, z, n;
    
!   _ODD_QtoV(7, vec, qn);
! }
! 
! void
! _limnQN15checker_QNtoV_d(double *vec, int qn) {
!   int ui, vi, zi;
!   double u, v, x, y, z, n;
!   
!   _ODD_QtoV(7, vec, qn);
  }
  
***************
*** 311,329 ****
    int ui, vi, zi;
    
!   x = vec[0];
!   y = vec[1];
!   z = vec[2];
!   L = AIR_ABS(x) + AIR_ABS(y) + AIR_ABS(z);
!   if (!L) {
!     return 0;
!   }
!   x /= L;
!   y /= L;
!   u = x + y;
!   v = x - y;
!   AIR_INDEX(-1, u, 1, 128, ui);
!   AIR_INDEX(-1, v, 1, 128, vi);
!   zi = (z > 0);
!   return (zi << 14) | (vi << 7) | ui;
  }
  
--- 332,344 ----
    int ui, vi, zi;
    
!   _ODD_VtoQ(7, vec);
! }
! 
! int
! _limnQN15checker_VtoQN_d(double *vec) {
!   double L, u, v, x, y, z;
!   int ui, vi, zi;
!   
!   _ODD_VtoQ(7, vec);
  }
  
***************
*** 365,368 ****
--- 380,417 ----
  
  void
+ _limnQN13checker_QNtoV_f(float *vec, int qn) {
+   int ui, vi, zi;
+   float u, v, x, y, z, n;
+   
+   _ODD_QtoV(6, vec, qn);
+ }
+ 
+ void
+ _limnQN13checker_QNtoV_d(double *vec, int qn) {
+   int ui, vi, zi;
+   double u, v, x, y, z, n;
+   
+   _ODD_QtoV(6, vec, qn);
+ }
+ 
+ int
+ _limnQN13checker_VtoQN_f(float *vec) {
+   float L, u, v, x, y, z;
+   int ui, vi, zi;
+   
+   _ODD_VtoQ(6, vec);
+ }
+ 
+ int
+ _limnQN13checker_VtoQN_d(double *vec) {
+   double L, u, v, x, y, z;
+   int ui, vi, zi;
+   
+   _ODD_VtoQ(6, vec);
+ }
+ 
+ /* ----------------------------------------------------------- */
+ 
+ void
  _limnQN12checker_QNtoV_f(float *vec, int qn) {
    int ui, vi;
***************
*** 398,401 ****
--- 447,585 ----
  /* ----------------------------------------------------------- */
  
+ void
+ _limnQN11checker_QNtoV_f(float *vec, int qn) {
+   int ui, vi, zi;
+   float u, v, x, y, z, n;
+   
+   _ODD_QtoV(5, vec, qn);
+ }
+ 
+ void
+ _limnQN11checker_QNtoV_d(double *vec, int qn) {
+   int ui, vi, zi;
+   double u, v, x, y, z, n;
+   
+   _ODD_QtoV(5, vec, qn);
+ }
+ 
+ int
+ _limnQN11checker_VtoQN_f(float *vec) {
+   float L, u, v, x, y, z;
+   int ui, vi, zi;
+   
+   _ODD_VtoQ(5, vec);
+ }
+ 
+ int
+ _limnQN11checker_VtoQN_d(double *vec) {
+   double L, u, v, x, y, z;
+   int ui, vi, zi;
+   
+   _ODD_VtoQ(5, vec);
+ }
+ 
+ /* ----------------------------------------------------------- */
+ 
+ void
+ _limnQN10checker_QNtoV_f(float *vec, int qn) {
+   int ui, vi;
+   double u, v, x, y, z, n;
+ 
+   _EVEN_QtoV(5, vec, qn);
+ }
+ 
+ void
+ _limnQN10checker_QNtoV_d(double *vec, int qn) {
+   int ui, vi;
+   double u, v, x, y, z, n;
+ 
+   _EVEN_QtoV(5, vec, qn);
+ }
+ 
+ int
+ _limnQN10checker_VtoQN_f(float *vec) {
+   double L, x, y, z;
+   int xi, yi, ui, vi;
+   
+   _EVEN_VtoQ(5, vec);
+ }
+ 
+ int
+ _limnQN10checker_VtoQN_d(double *vec) {
+   double L, x, y, z;
+   int xi, yi, ui, vi;
+   
+   _EVEN_VtoQ(5, vec);
+ }
+ 
+ /* ----------------------------------------------------------- */
+ void
+ _limnQN9checker_QNtoV_f(float *vec, int qn) {
+   int ui, vi, zi;
+   float u, v, x, y, z, n;
+   
+   _ODD_QtoV(4, vec, qn);
+ }
+ 
+ void
+ _limnQN9checker_QNtoV_d(double *vec, int qn) {
+   int ui, vi, zi;
+   double u, v, x, y, z, n;
+   
+   _ODD_QtoV(4, vec, qn);
+ }
+ 
+ int
+ _limnQN9checker_VtoQN_f(float *vec) {
+   float L, u, v, x, y, z;
+   int ui, vi, zi;
+   
+   _ODD_VtoQ(4, vec);
+ }
+ 
+ int
+ _limnQN9checker_VtoQN_d(double *vec) {
+   double L, u, v, x, y, z;
+   int ui, vi, zi;
+   
+   _ODD_VtoQ(4, vec);
+ }
+ 
+ /* ----------------------------------------------------------- */
+ 
+ void
+ _limnQN8checker_QNtoV_f(float *vec, int qn) {
+   int ui, vi;
+   double u, v, x, y, z, n;
+ 
+   _EVEN_QtoV(4, vec, qn);
+ }
+ 
+ void
+ _limnQN8checker_QNtoV_d(double *vec, int qn) {
+   int ui, vi;
+   double u, v, x, y, z, n;
+ 
+   _EVEN_QtoV(4, vec, qn);
+ }
+ 
+ int
+ _limnQN8checker_VtoQN_f(float *vec) {
+   double L, x, y, z;
+   int xi, yi, ui, vi;
+   
+   _EVEN_VtoQ(4, vec);
+ }
+ 
+ int
+ _limnQN8checker_VtoQN_d(double *vec) {
+   double L, x, y, z;
+   int xi, yi, ui, vi;
+   
+   _EVEN_VtoQ(4, vec);
+ }
+ 
+ /* ----------------------------------------------------------- */
+ 
  void (*
  limnQNtoV_f[LIMN_QN_MAX+1])(float *, int) = {
***************
*** 406,414 ****
    _limnQN15checker_QNtoV_f,
    _limnQN14checker_QNtoV_f,
    _limnQN12checker_QNtoV_f,
!   NULL,
!   NULL,
!   NULL,
!   NULL,
  };
    
--- 590,599 ----
    _limnQN15checker_QNtoV_f,
    _limnQN14checker_QNtoV_f,
+   _limnQN13checker_QNtoV_f,
    _limnQN12checker_QNtoV_f,
!   _limnQN11checker_QNtoV_f,
!   _limnQN10checker_QNtoV_f,
!   _limnQN9checker_QNtoV_f,
!   _limnQN8checker_QNtoV_f
  };
    
***************
*** 419,429 ****
    NULL,
    _limnQN16checker_QNtoV_d,
!   NULL,
    _limnQN14checker_QNtoV_d,
    _limnQN12checker_QNtoV_d,
!   NULL,
!   NULL,
!   NULL,
!   NULL,
  };
    
--- 604,615 ----
    NULL,
    _limnQN16checker_QNtoV_d,
!   _limnQN15checker_QNtoV_d,
    _limnQN14checker_QNtoV_d,
+   _limnQN13checker_QNtoV_d,
    _limnQN12checker_QNtoV_d,
!   _limnQN11checker_QNtoV_d,
!   _limnQN10checker_QNtoV_d,
!   _limnQN9checker_QNtoV_d,
!   _limnQN8checker_QNtoV_d
  };
    
***************
*** 436,445 ****
    _limnQN15checker_VtoQN_f,
    _limnQN14checker_VtoQN_f,
!   NULL,
    _limnQN12checker_VtoQN_f,
!   NULL,
!   NULL,
!   NULL,
!   NULL
  };
  
--- 622,631 ----
    _limnQN15checker_VtoQN_f,
    _limnQN14checker_VtoQN_f,
!   _limnQN13checker_VtoQN_f,
    _limnQN12checker_VtoQN_f,
!   _limnQN11checker_VtoQN_f,
!   _limnQN10checker_VtoQN_f,
!   _limnQN9checker_VtoQN_f,
!   _limnQN8checker_VtoQN_f
  };
  
***************
*** 450,461 ****
    NULL,
    _limnQN16checker_VtoQN_d,
!   NULL,
    _limnQN14checker_VtoQN_d,
!   NULL,
    _limnQN12checker_VtoQN_d,
!   NULL,
!   NULL,
!   NULL,
!   NULL
  };
  
--- 636,647 ----
    NULL,
    _limnQN16checker_VtoQN_d,
!   _limnQN15checker_VtoQN_d,
    _limnQN14checker_VtoQN_d,
!   _limnQN13checker_VtoQN_d,
    _limnQN12checker_VtoQN_d,
!   _limnQN11checker_VtoQN_d,
!   _limnQN10checker_VtoQN_d,
!   _limnQN9checker_VtoQN_d,
!   _limnQN8checker_VtoQN_d
  };
  
***************
*** 468,471 ****
--- 654,658 ----
    (1 << 15),
    (1 << 14),
+   (1 << 13),
    (1 << 12),
    (1 << 11),

[Teem-commit] teem/src/bin emap.c,1.8,1.9

[<kin...@us...>]

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

Modified Files:
	emap.c 
Log Message:
added a whole bunch more normal quantization methods

Index: emap.c
===================================================================
RCS file: /cvsroot/teem/teem/src/bin/emap.c,v
retrieving revision 1.8
retrieving revision 1.9
diff -C2 -d -r1.8 -r1.9
*** emap.c	27 Feb 2004 16:42:52 -0000	1.8
--- emap.c	27 Feb 2004 18:16:50 -0000	1.9
***************
*** 41,45 ****
    char *me, *outS, *errS, *debugS;
    airArray *mop;
!   float amb[3], *linfo, *debug, *map, W[3], V[3], diff;
    int li, ui, vi, qn, bits, method;
    limnLight *light;
--- 41,45 ----
    char *me, *outS, *errS, *debugS;
    airArray *mop;
!   float amb[3], *linfo, *debug, *map, W[3], V[3] /* , diff */;
    int li, ui, vi, qn, bits, method;
    limnLight *light;
***************
*** 146,149 ****
--- 146,150 ----
      ndebug = nrrdNew();
      nrrdMaybeAlloc(ndebug, nrrdTypeFloat, 3, 3, 1024, 512);
+     airMopAdd(mop, ndebug, (airMopper)nrrdNuke, airMopAlways);
      debug = ndebug->data;
      for (vi=0; vi<=511; vi++) {

[Teem-commit] teem/src release.txt,1.27,1.28

[<kin...@us...>]

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

Modified Files:
	release.txt 
Log Message:
added a whole bunch more normal quantization methods

Index: release.txt
===================================================================
RCS file: /cvsroot/teem/teem/src/release.txt,v
retrieving revision 1.27
retrieving revision 1.28
diff -C2 -d -r1.27 -r1.28
*** release.txt	29 Oct 2003 18:49:29 -0000	1.27
--- release.txt	27 Feb 2004 18:16:50 -0000	1.28
***************
*** 96,99 ****
--- 96,100 ----
  setenv TEEM_PNG_LPATH -L/home/sci/gk/usr/irix6.n32/lib
  setenv TEEM_PNG_IPATH -I/home/sci/gk/usr/irix6.n32/include
+ setenv PATH ${PATH}:/home/sci/gk/usr/local/irix6.n32/bin
  rapture/ray/muse:
  setenv TEEM_ARCH irix6.64

[Teem-commit] teem/src/nrrd apply1D.c,1.22,1.23

[<ik...@us...>]

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

Modified Files:
	apply1D.c 
Log Message:
added nrrdConvert to subst; could probably use airMop; waiting for TEAC blessing

Index: apply1D.c
===================================================================
RCS file: /cvsroot/teem/teem/src/nrrd/apply1D.c,v
retrieving revision 1.22
retrieving revision 1.23
diff -C2 -d -r1.22 -r1.23
*** apply1D.c	27 Feb 2004 06:05:16 -0000	1.22
--- apply1D.c	27 Feb 2004 18:09:57 -0000	1.23
***************
*** 876,885 ****
    char me[]="nrrdApply1DSubstitution", err[AIR_STRLEN_MED];
    double (*lup)(const void *, size_t);
-   double (*slup)(const void *, size_t);
    double (*ins)(void *, size_t, double);
!   double val, sval;
    size_t ii, num;
!   int jj, asize0, asize1;
!   int changed;
  
    if (!(nout && nsubst && nin)) {
--- 876,884 ----
    char me[]="nrrdApply1DSubstitution", err[AIR_STRLEN_MED];
    double (*lup)(const void *, size_t);
    double (*ins)(void *, size_t, double);
!   Nrrd *dsubst;
!   double val, sval, *d;
    size_t ii, num;
!   int jj, asize0, asize1, changed;
  
    if (!(nout && nsubst && nin)) {
***************
*** 909,923 ****
      }
    }
    lup = nrrdDLookup[nout->type];
-   slup = nrrdDLookup[nsubst->type];
    ins = nrrdDInsert[nout->type];
    num = nrrdElementNumber(nout);
    for (ii=0; ii<num; ii++) {
      val = lup(nout->data, ii);
      changed = AIR_FALSE;
      for (jj=0; jj<asize1; jj++) {
!       sval = slup(nsubst->data, jj*2+0);
        if (val == sval) {
! 	val = slup(nsubst->data, jj*2+1);
  	changed = AIR_TRUE;
        }
--- 908,926 ----
      }
    }
+   if (nrrdConvert(dsubst = nrrdNew(), nsubst, nrrdTypeDouble)) {
+     sprintf(err, "%s: couldn't create double copy of substitution table", me);
+     biffAdd(NRRD, err); dsubst = nrrdNuke(dsubst); return 1;
+   }
    lup = nrrdDLookup[nout->type];
    ins = nrrdDInsert[nout->type];
    num = nrrdElementNumber(nout);
+   d = (double *)dsubst->data;
    for (ii=0; ii<num; ii++) {
      val = lup(nout->data, ii);
      changed = AIR_FALSE;
      for (jj=0; jj<asize1; jj++) {
!       sval = d[jj*2+0];
        if (val == sval) {
! 	val = d[jj*2+1];
  	changed = AIR_TRUE;
        }
***************
*** 927,930 ****
--- 930,934 ----
      }
    }
+   dsubst = nrrdNuke(dsubst);
    return 0;
  }

[Teem-commit] teem/src/bin emap.c,1.7,1.8

[<kin...@us...>]

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

Modified Files:
	emap.c 
Log Message:
starting to expand number of quantized normal schemes

Index: emap.c
===================================================================
RCS file: /cvsroot/teem/teem/src/bin/emap.c,v
retrieving revision 1.7
retrieving revision 1.8
diff -C2 -d -r1.7 -r1.8
*** emap.c	9 Feb 2004 07:20:46 -0000	1.7
--- emap.c	27 Feb 2004 16:42:52 -0000	1.8
***************
*** 24,29 ****
  #include <teem/limn.h>
  
! char *emapInfo = ("Creates environment maps based on limn's 16checker "
! 		  "normal quantization method.  By taking into account "
  		  "camera parameters, this allows for both lights in "
  		  "both world and view space.  Solely out of laziness, "
--- 24,29 ----
  #include <teem/limn.h>
  
! char *emapInfo = ("Creates environment maps based on limn's \"checker\" "
! 		  "normal quantization methods.  By taking into account "
  		  "camera parameters, this allows for both lights in "
  		  "both world and view space.  Solely out of laziness, "
***************
*** 41,46 ****
    char *me, *outS, *errS, *debugS;
    airArray *mop;
!   float amb[3], *linfo, *debug, *map, W[3], V[3];
!   int li, ui, vi, qn;
    limnLight *light;
    limnCamera *cam;
--- 41,46 ----
    char *me, *outS, *errS, *debugS;
    airArray *mop;
!   float amb[3], *linfo, *debug, *map, W[3], V[3], diff;
!   int li, ui, vi, qn, bits, method;
    limnLight *light;
    limnCamera *cam;
***************
*** 57,60 ****
--- 57,63 ----
  	     "input nrrd containing light information",
  	     NULL, NULL, nrrdHestNrrd);
+   hestOptAdd(&hopt, "b", "# bits", airTypeInt, 1, 1, &bits, "16",
+ 	     "number of bits to use for normal quantization, "
+ 	     "between 8 and 16 inclusive. ");
    hestOptAdd(&hopt, "amb", "ambient RGB", airTypeFloat, 3, 3, amb, "0 0 0",
  	     "ambient light color");
***************
*** 80,83 ****
--- 83,102 ----
    airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
  
+   switch(bits) {
+   case 16: method = limnQN16checker; break;
+   case 15: method = limnQN15checker; break;
+   case 14: method = limnQN14checker; break;
+   case 13: method = limnQN13checker; break;
+   case 12: method = limnQN12checker; break;
+   case 11: method = limnQN11checker; break;
+   case 10: method = limnQN10checker; break;
+   case 9: method = limnQN9checker; break;
+   case 8: method = limnQN8checker; break;
+   default:
+     fprintf(stderr, "%s: requested #bits (%d) not in valid range [8,16]\n",
+ 	    me, bits);
+     airMopError(mop); return 1;
+   }
+ 
    if (!(nrrdTypeFloat == nlight->type &&
  	2 == nlight->dim && 
***************
*** 111,115 ****
    nmap=nrrdNew();
    airMopAdd(mop, nmap, (airMopper)nrrdNuke, airMopAlways);
!   if (limnEnvMapFill(nmap, limnLightDiffuseCB, limnQN16checker, light)) {
      airMopAdd(mop, errS = biffGetDone(LIMN), airFree, airMopAlways);
      fprintf(stderr, "%s: problem making environment map:\n%s\n", me, errS);
--- 130,134 ----
    nmap=nrrdNew();
    airMopAdd(mop, nmap, (airMopper)nrrdNuke, airMopAlways);
!   if (limnEnvMapFill(nmap, limnLightDiffuseCB, method, light)) {
      airMopAdd(mop, errS = biffGetDone(LIMN), airFree, airMopAlways);
      fprintf(stderr, "%s: problem making environment map:\n%s\n", me, errS);
***************
*** 141,152 ****
  	ELL_3V_SET(V, u, v, -w);
  	ELL_3MV_MUL(W, V2W, V);
! 	qn = limnVtoQN_f[limnQN16checker](W);
  	ELL_3V_COPY(debug + 3*(ui + 1024*vi), map + 3*qn);
  	
  	/* second, the far side of the sphere */
  	ELL_3V_SET(V, u, v, w);
  	ELL_3MV_MUL(W, V2W, V);
! 	qn = limnVtoQN_f[limnQN16checker](W);
  	ELL_3V_COPY(debug + 3*(ui + 512 + 1024*vi), map + 3*qn);
        }
      }
--- 160,185 ----
  	ELL_3V_SET(V, u, v, -w);
  	ELL_3MV_MUL(W, V2W, V);
! 	qn = limnVtoQN_f[method](W);
! 	/* 
! 	limnQNtoV_f[method](V, qn);
! 	ELL_3V_SUB(W, W, V);
! 	diff = ELL_3V_LEN(W);
! 	ELL_3V_SET(debug + 3*(ui + 1024*vi), diff, diff, diff);
! 	*/
  	ELL_3V_COPY(debug + 3*(ui + 1024*vi), map + 3*qn);
+ 
  	
  	/* second, the far side of the sphere */
  	ELL_3V_SET(V, u, v, w);
  	ELL_3MV_MUL(W, V2W, V);
! 	qn = limnVtoQN_f[method](W);
! 	/*
! 	limnQNtoV_f[method](V, qn);
! 	ELL_3V_SUB(W, W, V);
! 	diff = ELL_3V_LEN(W);
! 	ELL_3V_SET(debug + 3*(ui + 512 + 1024*vi), diff, diff, diff);
! 	*/
  	ELL_3V_COPY(debug + 3*(ui + 512 + 1024*vi), map + 3*qn);
+ 
        }
      }

[Teem-commit] teem/src/limn env.c,1.5,1.6 limn.h,1.49,1.50 qn.c,1.17,1.18

[<kin...@us...>]

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

Modified Files:
	env.c limn.h qn.c 
Log Message:
starting to expand number of quantized normal schemes

Index: env.c
===================================================================
RCS file: /cvsroot/teem/teem/src/limn/env.c,v
retrieving revision 1.5
retrieving revision 1.6
diff -C2 -d -r1.5 -r1.6
*** env.c	9 Feb 2004 07:06:18 -0000	1.5
--- env.c	27 Feb 2004 16:42:53 -0000	1.6
***************
*** 44,48 ****
      mapData = map->data;
      for (qn=0; qn<=sx*sy-1; qn++) {
!       limnQNtoV_f[limnQN16checker](vec, qn);
        cb(mapData + 3*qn, vec, data);
      }
--- 44,96 ----
      mapData = map->data;
      for (qn=0; qn<=sx*sy-1; qn++) {
!       limnQNtoV_f[qnMethod](vec, qn);
!       cb(mapData + 3*qn, vec, data);
!     }
!     break;
!   case limnQN14checker:
!     sx = sy = 128;
!     if (nrrdMaybeAlloc(map, nrrdTypeFloat, 3, 3, sx, sy)) {
!       sprintf(err, "%s: couldn't alloc output", me);
!       biffMove(LIMN, err, NRRD); return 1;
!     }
!     mapData = map->data;
!     for (qn=0; qn<=sx*sy-1; qn++) {
!       limnQNtoV_f[qnMethod](vec, qn);
!       cb(mapData + 3*qn, vec, data);
!     }
!     break;
!   case limnQN12checker:
!     sx = sy = 64;
!     if (nrrdMaybeAlloc(map, nrrdTypeFloat, 3, 3, sx, sy)) {
!       sprintf(err, "%s: couldn't alloc output", me);
!       biffMove(LIMN, err, NRRD); return 1;
!     }
!     mapData = map->data;
!     for (qn=0; qn<=sx*sy-1; qn++) {
!       limnQNtoV_f[qnMethod](vec, qn);
!       cb(mapData + 3*qn, vec, data);
!     }
!     break;
!   case limnQN10checker:
!     sx = sy = 32;
!     if (nrrdMaybeAlloc(map, nrrdTypeFloat, 3, 3, sx, sy)) {
!       sprintf(err, "%s: couldn't alloc output", me);
!       biffMove(LIMN, err, NRRD); return 1;
!     }
!     mapData = map->data;
!     for (qn=0; qn<=sx*sy-1; qn++) {
!       limnQNtoV_f[qnMethod](vec, qn);
!       cb(mapData + 3*qn, vec, data);
!     }
!     break;
!   case limnQN8checker:
!     sx = sy = 16;
!     if (nrrdMaybeAlloc(map, nrrdTypeFloat, 3, 3, sx, sy)) {
!       sprintf(err, "%s: couldn't alloc output", me);
!       biffMove(LIMN, err, NRRD); return 1;
!     }
!     mapData = map->data;
!     for (qn=0; qn<=sx*sy-1; qn++) {
!       limnQNtoV_f[qnMethod](vec, qn);
        cb(mapData + 3*qn, vec, data);
      }

Index: limn.h
===================================================================
RCS file: /cvsroot/teem/teem/src/limn/limn.h,v
retrieving revision 1.49
retrieving revision 1.50
diff -C2 -d -r1.49 -r1.50
*** limn.h	19 Feb 2004 03:44:07 -0000	1.49
--- limn.h	27 Feb 2004 16:42:53 -0000	1.50
***************
*** 270,283 ****
  */
  enum {
!   limnQNUnknown,     /* 0 */
!   limnQN16checker,   /* 1 */
!   limnQN16simple,    /* 2 */
!   limnQN16border1,   /* 3 */
!   limnQN15checker,   /* 4 */
!   limnQN14checker,   /* 5 */
!   limnQN12checker,   /* 6 */
    limnQNLast
  };
! #define LIMN_QN_MAX      6
  
  enum {
--- 270,288 ----
  */
  enum {
!   limnQNUnknown,     /*  0 */
!   limnQN16simple,    /*  1 */
!   limnQN16border1,   /*  2 */
!   limnQN16checker,   /*  3 */
!   limnQN15checker,   /*  4 */
!   limnQN14checker,   /*  5 */
!   limnQN13checker,   /*  6 */
!   limnQN12checker,   /*  7 */
!   limnQN11checker,   /*  8 */
!   limnQN10checker,   /*  9 */
!   limnQN9checker,    /* 10 */
!   limnQN8checker,    /* 11 */
    limnQNLast
  };
! #define LIMN_QN_MAX     11
  
  enum {
***************
*** 371,375 ****
  typedef void (*limnEnvMapCB)(float rgb[3], float vec[3], void *data);
  TEEM_API int limnEnvMapFill(Nrrd *envMap, limnEnvMapCB cb, 
! 			  int qnMethod, void *data);
  TEEM_API void limnLightDiffuseCB(float rgb[3], float vec[3], void *_lit);
  TEEM_API int limnEnvMapCheck(Nrrd *envMap);
--- 376,380 ----
  typedef void (*limnEnvMapCB)(float rgb[3], float vec[3], void *data);
  TEEM_API int limnEnvMapFill(Nrrd *envMap, limnEnvMapCB cb, 
! 			    int qnMethod, void *data);
  TEEM_API void limnLightDiffuseCB(float rgb[3], float vec[3], void *_lit);
  TEEM_API int limnEnvMapCheck(Nrrd *envMap);
***************
*** 408,414 ****
  				 float x, float y, float z);
  TEEM_API int limnObjectEdgeAdd(limnObject *obj, int partIdx, int lookIdx,
! 			       int faceIdxIdx, int vertIdxIdx0, int vertIdxIdx1);
! TEEM_API int limnObjectFaceAdd(limnObject *obj, int partIdx,
! 			       int lookIdx, int sideNum, int *vertIdxIdx);
  
  /* io.c */
--- 413,420 ----
  				 float x, float y, float z);
  TEEM_API int limnObjectEdgeAdd(limnObject *obj, int partIdx, int lookIdx,
! 			       int faceIdxIdx, int vertIdxIdx0,
! 			       int vertIdxIdx1);
! TEEM_API int limnObjectFaceAdd(limnObject *obj, int partIdx, int lookIdx,
! 			       int sideNum, int *vertIdxIdx);
  
  /* io.c */
***************
*** 425,431 ****
  TEEM_API int limnObjectPolarSphereAdd(limnObject *obj, int lookIdx, int axis,
  				      int thetaRes, int phiRes);
! TEEM_API int limnObjectConeAdd(limnObject *obj, int lookIdx, int axis, int res);
! TEEM_API int limnObjectPolarSuperquadAdd(limnObject *obj, int lookIdx, int axis,
! 					 float A, float B,
  					 int thetaRes, int phiRes);
  
--- 431,438 ----
  TEEM_API int limnObjectPolarSphereAdd(limnObject *obj, int lookIdx, int axis,
  				      int thetaRes, int phiRes);
! TEEM_API int limnObjectConeAdd(limnObject *obj, int lookIdx,
! 			       int axis, int res);
! TEEM_API int limnObjectPolarSuperquadAdd(limnObject *obj, int lookIdx,
! 					 int axis, float A, float B,
  					 int thetaRes, int phiRes);
  
***************
*** 435,439 ****
  TEEM_API int limnObjectSpaceTransform(limnObject *obj, limnCamera *cam,
  				      limnWindow *win, int space);
! TEEM_API int limnObjectPartTransform(limnObject *obj, int partIdx, float tx[16]);
  TEEM_API int limnObjectDepthSortParts(limnObject *obj);
  TEEM_API int limnObjectDepthSortFaces(limnObject *obj);
--- 442,447 ----
  TEEM_API int limnObjectSpaceTransform(limnObject *obj, limnCamera *cam,
  				      limnWindow *win, int space);
! TEEM_API int limnObjectPartTransform(limnObject *obj, int partIdx,
! 				     float tx[16]);
  TEEM_API int limnObjectDepthSortParts(limnObject *obj);
  TEEM_API int limnObjectDepthSortFaces(limnObject *obj);
***************
*** 441,445 ****
  
  /* renderLimn.c */
! TEEM_API int limnObjectRender(limnObject *obj, limnCamera *cam, limnWindow *win);
  TEEM_API int limnObjectPSDraw(limnObject *obj, limnCamera *cam,
  			      Nrrd *envMap, limnWindow *win);
--- 449,454 ----
  
  /* renderLimn.c */
! TEEM_API int limnObjectRender(limnObject *obj, limnCamera *cam,
! 			      limnWindow *win);
  TEEM_API int limnObjectPSDraw(limnObject *obj, limnCamera *cam,
  			      Nrrd *envMap, limnWindow *win);
***************
*** 453,457 ****
  				   limnSplineTypeSpec *spec);
  TEEM_API limnSpline *limnSplineNix(limnSpline *spline);
! TEEM_API int limnSplineNrrdCleverFix(Nrrd *nout, Nrrd *nin, int info, int type);
  TEEM_API limnSpline *limnSplineCleverNew(Nrrd *ncpt, int info,
  					 limnSplineTypeSpec *spec);
--- 462,467 ----
  				   limnSplineTypeSpec *spec);
  TEEM_API limnSpline *limnSplineNix(limnSpline *spline);
! TEEM_API int limnSplineNrrdCleverFix(Nrrd *nout, Nrrd *nin,
! 				     int info, int type);
  TEEM_API limnSpline *limnSplineCleverNew(Nrrd *ncpt, int info,
  					 limnSplineTypeSpec *spec);

Index: qn.c
===================================================================
RCS file: /cvsroot/teem/teem/src/limn/qn.c,v
retrieving revision 1.17
retrieving revision 1.18
diff -C2 -d -r1.17 -r1.18
*** qn.c	9 Feb 2004 07:06:18 -0000	1.17
--- qn.c	27 Feb 2004 16:43:04 -0000	1.18
***************
*** 91,94 ****
--- 91,151 ----
  /* ----------------------------------------------------------------  */
  
+ void
+ _limnQN16border1_QNtoV_f(float *vec, int qn) {
+   int ui, vi;
+   float u, v, x, y, z, n;
+   
+   ui = qn & 0xFF;
+   vi = qn >> 8;
+   u = AIR_AFFINE(0.5, ui, 254.5, -0.5, 0.5);
+   v = AIR_AFFINE(0.5, vi, 254.5, -0.5, 0.5);
+   x =  u + v;
+   y =  u - v;
+   z = 1 - AIR_ABS(x) - AIR_ABS(y);
+   z *= (((ui ^ vi) & 0x01) << 1) - 1;
+   n = 1.0/sqrt(x*x + y*y + z*z);
+   vec[0] = x*n; 
+   vec[1] = y*n; 
+   vec[2] = z*n;
+ }
+ 
+ int
+ _limnQN16border1_VtoQN_f(float *vec) {
+   float L, u, v, x, y, z;
+   int ui, vi, zi;
+   char me[]="limnQNVto16PB1";
+   
+   x = vec[0];
+   y = vec[1];
+   z = vec[2];
+   L = AIR_ABS(x) + AIR_ABS(y) + AIR_ABS(z);
+   if (!L) {
+     return 0;
+   }
+   x /= L;
+   y /= L;
+   u = x + y;
+   v = x - y;
+   AIR_INDEX(-1, u, 1, 254, ui); ui++;
+   AIR_INDEX(-1, v, 1, 254, vi); vi++;
+   zi = (ui ^ vi) & 0x01;
+   if (!zi && z > 1.0/128.0) {
+     ui -= (((ui >> 7) & 0x01) << 1) - 1;
+   } 
+   else if (zi && z < -1.0/128.0) {
+     vi -= (((vi >> 7) & 0x01) << 1) - 1;
+   }
+   zi = (ui ^ vi) & 0x01;
+   if (!zi && z > 1.0/127.0) {
+     printf("%s: panic01\n", me);
+   } 
+   else if (zi && z < -1.0/127.0) {
+     printf("%s: panic02\n", me);
+   }
+   return (vi << 8) | ui;
+ }
+ 
+ /* ----------------------------------------------------------------  */
+ 
    /*
      x =  [  1.0   1.0 ] u
***************
*** 103,109 ****
    /* xor of low bits == 0 --> z<0 ; xor == 1 --> z>0 */
  
! /* May 11 2003 GK verified that except at equator seam,
     the error due to quantization is unbiased */
  
  #define _16_QtoV(vec, qn) \
    ui = (qn) & 0xFF; \
--- 160,207 ----
    /* xor of low bits == 0 --> z<0 ; xor == 1 --> z>0 */
  
! /* May 11 2003 GK (visually) verified that except at equator seam,
     the error due to quantization is unbiased */
  
+ #define _EVEN_QtoV(HNB, vec, qn) \
+   ui = (qn) & ((1<<HNB)-1); \
+   vi = ((qn) >> HNB) & ((1<<HNB)-1); \
+   u = AIR_AFFINE(0, ui, ((1<<HNB)-1), -0.5, 0.5); \
+   v = AIR_AFFINE(0, vi, ((1<<HNB)-1), -0.5, 0.5); \
+   x =  u + v; \
+   y =  u - v; \
+   z = 1 - AIR_ABS(x) - AIR_ABS(y); \
+   /* would this be better served by a branch? */ \
+   z *= (((ui ^ vi) & 0x01) << 1) - 1; \
+   n = 1.0/sqrt(x*x + y*y + z*z); \
+   (vec)[0] = x*n; \
+   (vec)[1] = y*n; \
+   (vec)[2] = z*n
+ 
+ #define _EVEN_VtoQ(HNB, vec) \
+   x = (vec)[0]; \
+   y = (vec)[1]; \
+   z = (vec)[2]; \
+   L = AIR_ABS(x) + AIR_ABS(y) + AIR_ABS(z); \
+   if (!L) { \
+     return 0; \
+   } \
+   x /= L; \
+   y /= L; \
+   if (z > 0) { \
+     AIR_INDEX(-1.0, x, 1.0, ((1<<HNB)-1), xi); \
+     AIR_INDEX(-1.0 - 1.0/((1<<HNB)-1), y, 1.0 + 1.0/((1<<HNB)-1), \
+               (1<<HNB), yi); \
+     ui = xi + yi - ((1<<(HNB-1))-1); \
+     vi = xi - yi + (1<<(HNB-1)); \
+   } \
+   else { \
+     AIR_INDEX(-1.0 - 1.0/((1<<HNB)-1), x, 1.0 + 1.0/((1<<HNB)-1), \
+               (1<<HNB), xi); \
+     AIR_INDEX(-1, y, 1, ((1<<HNB)-1), yi); \
+     ui = xi + yi - ((1<<(HNB-1))-1); \
+     vi = xi - yi + ((1<<(HNB-1))-1); \
+   } \
+   return (vi << HNB) | ui
+ 
  #define _16_QtoV(vec, qn) \
    ui = (qn) & 0xFF; \
***************
*** 126,130 ****
    double u, v, x, y, z, n;
  
!   _16_QtoV(vec, qn);
  }
  
--- 224,229 ----
    double u, v, x, y, z, n;
  
!   /* _16_QtoV(vec, qn); */
!   _EVEN_QtoV(8, vec, qn);
  }
  
***************
*** 134,138 ****
    double u, v, x, y, z, n;
  
!   _16_QtoV(vec, qn);
  }
  
--- 233,238 ----
    double u, v, x, y, z, n;
  
!   /* _16_QtoV(vec, qn); */
!   _EVEN_QtoV(8, vec, qn);
  }
  
***************
*** 164,170 ****
    double L, x, y, z;
    int xi, yi, ui, vi;
! 
    _16_VtoQ(vec);
    return (vi << 8) | ui;
  }
  
--- 264,273 ----
    double L, x, y, z;
    int xi, yi, ui, vi;
!   
!   /*
    _16_VtoQ(vec);
    return (vi << 8) | ui;
+   */
+   _EVEN_VtoQ(8, vec);
  }
  
***************
*** 174,236 ****
    int xi, yi, ui, vi;
  
    _16_VtoQ(vec);
    return (vi << 8) | ui;
! }
! 
! /* ----------------------------------------------------------------  */
! 
! void
! _limnQN16border1_QNtoV_f(float *vec, int qn) {
!   int ui, vi;
!   float u, v, x, y, z, n;
!   
!   ui = qn & 0xFF;
!   vi = qn >> 8;
!   u = AIR_AFFINE(0.5, ui, 254.5, -0.5, 0.5);
!   v = AIR_AFFINE(0.5, vi, 254.5, -0.5, 0.5);
!   x =  u + v;
!   y =  u - v;
!   z = 1 - AIR_ABS(x) - AIR_ABS(y);
!   z *= (((ui ^ vi) & 0x01) << 1) - 1;
!   n = 1.0/sqrt(x*x + y*y + z*z);
!   vec[0] = x*n; 
!   vec[1] = y*n; 
!   vec[2] = z*n;
! }
! 
! int
! _limnQN16border1_VtoQN_f(float *vec) {
!   float L, u, v, x, y, z;
!   int ui, vi, zi;
!   char me[]="limnQNVto16PB1";
!   
!   x = vec[0];
!   y = vec[1];
!   z = vec[2];
!   L = AIR_ABS(x) + AIR_ABS(y) + AIR_ABS(z);
!   if (!L) {
!     return 0;
!   }
!   x /= L;
!   y /= L;
!   u = x + y;
!   v = x - y;
!   AIR_INDEX(-1, u, 1, 254, ui); ui++;
!   AIR_INDEX(-1, v, 1, 254, vi); vi++;
!   zi = (ui ^ vi) & 0x01;
!   if (!zi && z > 1.0/128.0) {
!     ui -= (((ui >> 7) & 0x01) << 1) - 1;
!   } 
!   else if (zi && z < -1.0/128.0) {
!     vi -= (((vi >> 7) & 0x01) << 1) - 1;
!   }
!   zi = (ui ^ vi) & 0x01;
!   if (!zi && z > 1.0/127.0) {
!     printf("%s: panic01\n", me);
!   } 
!   else if (zi && z < -1.0/127.0) {
!     printf("%s: panic02\n", me);
!   }
!   return (vi << 8) | ui;
  }
  
--- 277,285 ----
    int xi, yi, ui, vi;
  
+   /*
    _16_VtoQ(vec);
    return (vi << 8) | ui;
!   */
!   _EVEN_VtoQ(8, vec);
  }
  
***************
*** 281,299 ****
  /* ----------------------------------------------------------- */
  
- #define _14_QNtoV(vec, qn) \
-   ui = (qn) & 0x7F; \
-   vi = ((qn) >> 7) & 0x7F; \
-   u = AIR_AFFINE(0, ui, 127, -0.5, 0.5); \
-   v = AIR_AFFINE(0, vi, 127, -0.5, 0.5); \
-   x =  u + v; \
-   y =  u - v; \
-   z = 1 - AIR_ABS(x) - AIR_ABS(y); \
-   /* would this be better served by a branch? */ \
-   z *= (((ui ^ vi) & 0x01) << 1) - 1; \
-   n = 1.0/sqrt(x*x + y*y + z*z); \
-   (vec)[0] = x*n; \
-   (vec)[1] = y*n; \
-   (vec)[2] = z*n
- 
  void
  _limnQN14checker_QNtoV_f(float *vec, int qn) {
--- 330,333 ----
***************
*** 301,305 ****
    double u, v, x, y, z, n;
  
!   _14_QNtoV(vec, qn);
  }
  
--- 335,339 ----
    double u, v, x, y, z, n;
  
!   _EVEN_QtoV(7, vec, qn);
  }
  
***************
*** 309,338 ****
    double u, v, x, y, z, n;
  
!   _14_QNtoV(vec, qn);
  }
  
- #define _14_VtoQN(vec) \
-   x = (vec)[0]; \
-   y = (vec)[1]; \
-   z = (vec)[2]; \
-   L = AIR_ABS(x) + AIR_ABS(y) + AIR_ABS(z); \
-   if (!L) { \
-     return 0; \
-   } \
-   x /= L; \
-   y /= L; \
-   if (z > 0) { \
-     AIR_INDEX(-1.0, x, 1.0, 127, xi); \
-     AIR_INDEX(-1.0 - 1.0/127, y, 1.0 + 1.0/127, 128, yi); \
-     ui = xi + yi - 63; \
-     vi = xi - yi + 64; \
-   } \
-   else { \
-     AIR_INDEX(-1.0 - 1.0/127, x, 1.0 + 1.0/127, 128, xi); \
-     AIR_INDEX(-1, y, 1, 127, yi); \
-     ui = xi + yi - 63; \
-     vi = xi - yi + 63; \
-   }
- 
  int
  _limnQN14checker_VtoQN_f(float *vec) {
--- 343,349 ----
    double u, v, x, y, z, n;
  
!   _EVEN_QtoV(7, vec, qn);
  }
  
  int
  _limnQN14checker_VtoQN_f(float *vec) {
***************
*** 340,345 ****
    int xi, yi, ui, vi;
  
!   _14_VtoQN(vec);
!   return (vi << 7) | ui;
  }
  
--- 351,355 ----
    int xi, yi, ui, vi;
  
!   _EVEN_VtoQ(7, vec);
  }
  
***************
*** 349,373 ****
    int xi, yi, ui, vi;
    
!   _14_VtoQN(vec);
!   return (vi << 7) | ui;
  }
  
  /* ----------------------------------------------------------- */
  
- #define _12_QNtoV(vec, qn) \
-   ui = (qn) & 0x3F; \
-   vi = ((qn) >> 6) & 0x3F; \
-   u = AIR_AFFINE(0, ui, 63, -0.5, 0.5); \
-   v = AIR_AFFINE(0, vi, 63, -0.5, 0.5); \
-   x =  u + v; \
-   y =  u - v; \
-   z = 1 - AIR_ABS(x) - AIR_ABS(y); \
-   /* would this be better served by a branch? */ \
-   z *= (((ui ^ vi) & 0x01) << 1) - 1; \
-   n = 1.0/sqrt(x*x + y*y + z*z); \
-   (vec)[0] = x*n; \
-   (vec)[1] = y*n; \
-   (vec)[2] = z*n
- 
  void
  _limnQN12checker_QNtoV_f(float *vec, int qn) {
--- 359,367 ----
    int xi, yi, ui, vi;
    
!   _EVEN_VtoQ(7, vec);
  }
  
  /* ----------------------------------------------------------- */
  
  void
  _limnQN12checker_QNtoV_f(float *vec, int qn) {
***************
*** 375,379 ****
    double u, v, x, y, z, n;
  
!   _12_QNtoV(vec, qn);
  }
  
--- 369,373 ----
    double u, v, x, y, z, n;
  
!   _EVEN_QtoV(6, vec, qn);
  }
  
***************
*** 383,412 ****
    double u, v, x, y, z, n;
  
!   _12_QNtoV(vec, qn);
  }
  
- #define _12_VtoQN(vec) \
-   x = vec[0]; \
-   y = vec[1]; \
-   z = vec[2]; \
-   L = AIR_ABS(x) + AIR_ABS(y) + AIR_ABS(z); \
-   if (!L) { \
-     return 0; \
-   } \
-   x /= L; \
-   y /= L; \
-   if (z > 0) { \
-     AIR_INDEX(-1.0, x, 1.0, 63, xi); \
-     AIR_INDEX(-1.0 - 1.0/63, y, 1.0 + 1.0/63, 64, yi); \
-     ui = xi + yi - 31; \
-     vi = xi - yi + 32; \
-   } \
-   else { \
-     AIR_INDEX(-1.0 - 1.0/63, x, 1.0 + 1.0/63, 64, xi); \
-     AIR_INDEX(-1, y, 1, 63, yi); \
-     ui = xi + yi - 31; \
-     vi = xi - yi + 31; \
-   }
- 
  int
  _limnQN12checker_VtoQN_f(float *vec) {
--- 377,383 ----
    double u, v, x, y, z, n;
  
!   _EVEN_QtoV(6, vec, qn);
  }
  
  int
  _limnQN12checker_VtoQN_f(float *vec) {
***************
*** 414,419 ****
    int xi, yi, ui, vi;
    
!   _12_VtoQN(vec);
!   return (vi << 6) | ui;
  }
  
--- 385,389 ----
    int xi, yi, ui, vi;
    
!   _EVEN_VtoQ(6, vec);
  }
  
***************
*** 423,428 ****
    int xi, yi, ui, vi;
    
!   _12_VtoQN(vec);
!   return (vi << 6) | ui;
  }
  
--- 393,397 ----
    int xi, yi, ui, vi;
    
!   _EVEN_VtoQ(6, vec);
  }
  
***************
*** 432,441 ****
  limnQNtoV_f[LIMN_QN_MAX+1])(float *, int) = {
    NULL,
-   _limnQN16checker_QNtoV_f,
    _limnQN16simple_QNtoV_f,
    _limnQN16border1_QNtoV_f,
    _limnQN15checker_QNtoV_f,
    _limnQN14checker_QNtoV_f,
!   _limnQN12checker_QNtoV_f
  };
    
--- 401,414 ----
  limnQNtoV_f[LIMN_QN_MAX+1])(float *, int) = {
    NULL,
    _limnQN16simple_QNtoV_f,
    _limnQN16border1_QNtoV_f,
+   _limnQN16checker_QNtoV_f,
    _limnQN15checker_QNtoV_f,
    _limnQN14checker_QNtoV_f,
!   _limnQN12checker_QNtoV_f,
!   NULL,
!   NULL,
!   NULL,
!   NULL,
  };
    
***************
*** 443,452 ****
  limnQNtoV_d[LIMN_QN_MAX+1])(double *, int) = {
    NULL,
-   _limnQN16checker_QNtoV_d,
    NULL,
    NULL,
    NULL,
    _limnQN14checker_QNtoV_d,
!   _limnQN12checker_QNtoV_d
  };
    
--- 416,429 ----
  limnQNtoV_d[LIMN_QN_MAX+1])(double *, int) = {
    NULL,
    NULL,
    NULL,
+   _limnQN16checker_QNtoV_d,
    NULL,
    _limnQN14checker_QNtoV_d,
!   _limnQN12checker_QNtoV_d,
!   NULL,
!   NULL,
!   NULL,
!   NULL,
  };
    
***************
*** 454,463 ****
  limnVtoQN_f[LIMN_QN_MAX+1])(float *vec) = {
    NULL,
-   _limnQN16checker_VtoQN_f,
    _limnQN16simple_VtoQN_f,
    _limnQN16border1_VtoQN_f,
    _limnQN15checker_VtoQN_f,
    _limnQN14checker_VtoQN_f,
!   _limnQN12checker_VtoQN_f
  };
  
--- 431,445 ----
  limnVtoQN_f[LIMN_QN_MAX+1])(float *vec) = {
    NULL,
    _limnQN16simple_VtoQN_f,
    _limnQN16border1_VtoQN_f,
+   _limnQN16checker_VtoQN_f,
    _limnQN15checker_VtoQN_f,
    _limnQN14checker_VtoQN_f,
!   NULL,
!   _limnQN12checker_VtoQN_f,
!   NULL,
!   NULL,
!   NULL,
!   NULL
  };
  
***************
*** 465,474 ****
  limnVtoQN_d[LIMN_QN_MAX+1])(double *vec) = {
    NULL,
-   _limnQN16checker_VtoQN_d,
    NULL,
    NULL,
    NULL,
    _limnQN14checker_VtoQN_d,
!   _limnQN12checker_VtoQN_d
  };
  
--- 447,461 ----
  limnVtoQN_d[LIMN_QN_MAX+1])(double *vec) = {
    NULL,
    NULL,
    NULL,
+   _limnQN16checker_VtoQN_d,
    NULL,
    _limnQN14checker_VtoQN_d,
!   NULL,
!   _limnQN12checker_VtoQN_d,
!   NULL,
!   NULL,
!   NULL,
!   NULL
  };
  
***************
*** 478,484 ****
    (1 << 16),
    (1 << 16),
    (1 << 15),
    (1 << 14),
!   (1 << 12)
  };
  
--- 465,476 ----
    (1 << 16),
    (1 << 16),
+   (1 << 16),
    (1 << 15),
    (1 << 14),
!   (1 << 12),
!   (1 << 11),
!   (1 << 10),
!   (1 << 9),
!   (1 << 8)
  };
  

[Teem-commit] teem/src/ten tendGrads.c,1.5,1.6

[<kin...@us...>]

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

Modified Files:
	tendGrads.c 
Log Message:
starting to expand number of quantized normal schemes

Index: tendGrads.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/tendGrads.c,v
retrieving revision 1.5
retrieving revision 1.6
diff -C2 -d -r1.5 -r1.6
*** tendGrads.c	27 Feb 2004 13:02:51 -0000	1.5
--- tendGrads.c	27 Feb 2004 16:43:05 -0000	1.6
***************
*** 31,35 ****
     "faster.  In the mean time, you can run a second pass on the output of "
     "the first pass, using lower drag. A second phase of the algorithm "
!    "tries random signs in gradient directions in trying to find an optimally "
     "balanced set of directions.  This uses a randomized search, so if it "
     "doesn't seem to be finishing in a reasonable amount of time, try "
--- 31,35 ----
     "faster.  In the mean time, you can run a second pass on the output of "
     "the first pass, using lower drag. A second phase of the algorithm "
!    "tries sign changes in gradient directions in trying to find an optimally "
     "balanced set of directions.  This uses a randomized search, so if it "
     "doesn't seem to be finishing in a reasonable amount of time, try "

[Teem-commit] teem/src/mite txf.c,1.15,1.16

[<kin...@us...>]

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

Modified Files:
	txf.c 
Log Message:
starting to expand number of quantized normal schemes

Index: txf.c
===================================================================
RCS file: /cvsroot/teem/teem/src/mite/txf.c,v
retrieving revision 1.15
retrieving revision 1.16
diff -C2 -d -r1.15 -r1.16
*** txf.c	26 Feb 2004 18:41:14 -0000	1.15
--- txf.c	27 Feb 2004 16:43:04 -0000	1.16
***************
*** 203,216 ****
    }
    for (axi=1; axi<ntxf->dim; axi++) {
-     if (!( AIR_EXISTS(ntxf->axis[axi].min) && 
- 	   AIR_EXISTS(ntxf->axis[axi].max) )) {
-       sprintf(err, "%s: min and max of axis %d aren't both set", me, axi);
-       biffAdd(MITE, err); return 1;
-     }
-     if (!( ntxf->axis[axi].min < ntxf->axis[axi].max )) {
-       sprintf(err, "%s: min (%g) not less than max (%g) on axis %d", 
- 	      me, ntxf->axis[axi].min, ntxf->axis[axi].max, axi);
-       biffAdd(MITE, err); return 1;
-     }
      if (1 == ntxf->axis[axi].size) {
        sprintf(err, "%s: # samples on axis %d must be > 1", me, axi);
--- 203,206 ----
***************
*** 245,248 ****
--- 235,249 ----
  	biffAdd(MITE, err); return 1;
        }
+     } else {
+       if (!( AIR_EXISTS(ntxf->axis[axi].min) && 
+ 	     AIR_EXISTS(ntxf->axis[axi].max) )) {
+ 	sprintf(err, "%s: min and max of axis %d aren't both set", me, axi);
+ 	biffAdd(MITE, err); return 1;
+       }
+       if (!( ntxf->axis[axi].min < ntxf->axis[axi].max )) {
+ 	sprintf(err, "%s: min (%g) not less than max (%g) on axis %d", 
+ 		me, ntxf->axis[axi].min, ntxf->axis[axi].max, axi);
+ 	biffAdd(MITE, err); return 1;
+       }
      }
    }

[Teem-commit] teem/src/ten tendGrads.c,1.4,1.5

[<kin...@us...>]

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

Modified Files:
	tendGrads.c 
Log Message:
improved usage info

Index: tendGrads.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/tendGrads.c,v
retrieving revision 1.4
retrieving revision 1.5
diff -C2 -d -r1.4 -r1.5
*** tendGrads.c	27 Feb 2004 12:50:37 -0000	1.4
--- tendGrads.c	27 Feb 2004 13:02:51 -0000	1.5
***************
*** 32,36 ****
     "the first pass, using lower drag. A second phase of the algorithm "
     "tries random signs in gradient directions in trying to find an optimally "
!    "balanced set of directions.");
  
  int
--- 32,38 ----
     "the first pass, using lower drag. A second phase of the algorithm "
     "tries random signs in gradient directions in trying to find an optimally "
!    "balanced set of directions.  This uses a randomized search, so if it "
!    "doesn't seem to be finishing in a reasonable amount of time, try "
!    "restarting.");
  
  int

[Teem-commit] teem/src/ten grads.c,1.1,1.2 ten.h,1.77,1.78 tendGrads.c,1.3,1.4

[<kin...@us...>]

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

Modified Files:
	grads.c ten.h tendGrads.c 
Log Message:
added second phase of tenGradientDistribute in which randomized signs of gradients are tried to minimize the length of the mean gradient vector

Index: grads.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/grads.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -C2 -d -r1.1 -r1.2
*** grads.c	27 Feb 2004 10:28:07 -0000	1.1
--- grads.c	27 Feb 2004 12:50:36 -0000	1.2
***************
*** 31,36 ****
      ret->drag = 0.01;
      ret->dt = 0.05;
-     ret->minVelocity = 0.00001;
      ret->jitter = 0.05;
      ret->srand = AIR_TRUE;
      ret->snap = 0;
--- 31,38 ----
      ret->drag = 0.01;
      ret->dt = 0.05;
      ret->jitter = 0.05;
+     ret->minVelocity = 0.000001;
+     ret->minMean = 0.0001;
+     ret->minMeanImprovement = 0.00001;
      ret->srand = AIR_TRUE;
      ret->snap = 0;
***************
*** 241,250 ****
  }
  
  /*
  ******** tenGradientDistribute
  **
! ** takes the given list of gradients, normalizes their lengths, optionally jitters
! ** their positions, does point repulsion, and then selects a combination of
! ** directions with minimum vector sum.
  */
  int
--- 243,300 ----
  }
  
+ double
+ _tenGradientParty(double *grad, int num) {
+   double mean[3];
+   int ii;
+   
+   ELL_3V_SET(mean, 0, 0, 0);
+   for (ii=0; ii<num; ii++) {
+     if (airRandInt(2)) {
+       ELL_3V_SCALE(grad + 3*ii, -1, grad + 3*ii);
+     }
+     ELL_3V_INCR(mean, grad + 3*ii);
+   }
+   ELL_3V_SCALE(mean, 1.0/num, mean);
+   return ELL_3V_LEN(mean);
+ }
+ 
+ int
+ tenGradientMeanMinimize(Nrrd *nout, Nrrd *nin, tenGradientParm *tgparm) {
+   char me[]="tenGradientMeanMinimize", err[AIR_STRLEN_MED];
+   int num;
+   double *grad, len, lastLen, improv;
+ 
+   if (!nout || tenGradientCheck(nin, nrrdTypeUnknown, 2)) {
+     sprintf(err, "%s: got NULL pointer or invalid input", me);
+     biffAdd(TEN, err); return 1;
+   }
+   if (nrrdConvert(nout, nin, nrrdTypeDouble)) {
+     sprintf(err, "%s: can't initialize output with input", me);
+     biffMove(TEN, err, NRRD); return 1;
+   }
+   num = nout->axis[1].size;
+   grad = (double*)(nout->data);
+    
+   lastLen = _tenGradientParty(grad, num);
+   do {
+     do {
+       len = _tenGradientParty(grad, num);
+     } while (len > lastLen);
+     improv = lastLen - len;
+     lastLen = len;
+     fprintf(stderr, "%s: improvement: %g  (mean length = %g)\n",
+ 	    me, improv, len);
+   } while (improv > tgparm->minMeanImprovement
+ 	   && len > tgparm->minMean);
+   
+   return 0;
+ }
+ 
  /*
  ******** tenGradientDistribute
  **
! ** takes the given list of gradients, normalizes their lengths,
! ** optionally jitters their positions, does point repulsion, and then
! ** selects a combination of directions with minimum vector sum.
  */
  int
***************
*** 344,359 ****
        }
      } else {
!       if (!(iter % 25)) {
! 	fprintf(stderr, "%s: %d: meanVelocity = %g\n", me, iter, meanVelocity);
        }
      }
    }
  
!   /* select output gradients
!      - either do the trivial thing
!      - or run through 2^(#grads) to find one with maximum covariance
!      print covariance matrix eigenvalues
!   */
!   nrrdCopy(nout, npos);
  
    airMopOkay(mop); 
--- 394,409 ----
        }
      } else {
!       if (!(iter % 1000)) {
! 	fprintf(stderr, "%s: iteration = %d: meanVelocity = %g\n",
! 		me, iter, meanVelocity);
        }
      }
    }
  
!   fprintf(stderr, "%s: optimizing balance ... \n", me);
!   if (tenGradientMeanMinimize(nout, npos, tgparm)) {
!     sprintf(err, "%s: failed to minimize vector sum of gradients", me);
!     biffAdd(TEN, err); return 1;
!   }
  
    airMopOkay(mop); 
***************
*** 372,376 ****
    }
    if (!( num >= 3 )) {
!     sprintf(err, "%s: can generate at minimum 3 gradient directions (not %d)", me, num);
      biffAdd(TEN, err); return 1;
    }
--- 422,427 ----
    }
    if (!( num >= 3 )) {
!     sprintf(err, "%s: can generate minimum of 3 gradient directions "
! 	    "(not %d)", me, num);
      biffAdd(TEN, err); return 1;
    }
***************
*** 388,392 ****
    return 0;
  }
- 
- 
- 
--- 439,440 ----

Index: ten.h
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/ten.h,v
retrieving revision 1.77
retrieving revision 1.78
diff -C2 -d -r1.77 -r1.78
*** ten.h	27 Feb 2004 10:23:56 -0000	1.77
--- ten.h	27 Feb 2004 12:50:37 -0000	1.78
***************
*** 378,383 ****
      drag, 
      dt,
      minVelocity,
!     jitter; 
    int srand, snap, single, minIteration, maxIteration;
  } tenGradientParm;
--- 378,385 ----
      drag, 
      dt,
+     jitter,
      minVelocity,
!     minMean,
!     minMeanImprovement;
    int srand, snap, single, minIteration, maxIteration;
  } tenGradientParm;
***************
*** 400,405 ****
  TEEM_API int tenGradientRandom(Nrrd *ngrad, int num, int srand);
  TEEM_API int tenGradientJitter(Nrrd *nout, Nrrd *nin, double dist);
! TEEM_API int tenGradientDistribute(Nrrd *nout, Nrrd *nin, tenGradientParm *tgparm);
! TEEM_API int tenGradientGenerate(Nrrd *nout, int num, tenGradientParm *tgparm);
  
  /* enumsTen.c */
--- 402,411 ----
  TEEM_API int tenGradientRandom(Nrrd *ngrad, int num, int srand);
  TEEM_API int tenGradientJitter(Nrrd *nout, Nrrd *nin, double dist);
! TEEM_API int tenGradientMeanMinimize(Nrrd *nout, Nrrd *nin,
! 				     tenGradientParm *tgparm);
! TEEM_API int tenGradientDistribute(Nrrd *nout, Nrrd *nin,
! 				   tenGradientParm *tgparm);
! TEEM_API int tenGradientGenerate(Nrrd *nout, int num,
! 				 tenGradientParm *tgparm);
  
  /* enumsTen.c */

Index: tendGrads.c
===================================================================
RCS file: /cvsroot/teem/teem/src/ten/tendGrads.c,v
retrieving revision 1.3
retrieving revision 1.4
diff -C2 -d -r1.3 -r1.4
*** tendGrads.c	27 Feb 2004 10:41:03 -0000	1.3
--- tendGrads.c	27 Feb 2004 12:50:37 -0000	1.4
***************
*** 21,33 ****
  #include "privateTen.h"
  
! #define INFO "Calculate gradient directions for DWI acquisition"
  char *_tend_gradsInfoL =
    (INFO
     ", based on a simulation of anti-podal point pairs repelling each other "
!    "on the unit sphere surface. This can either distribute more uniformly a given "
!    "set of gradients, or it can make a new distribution from scratch.  A more "
!    "clever implementation could decrease drag with time, as the solution converges, "
!    "to get closer to the minimum energy configuration faster.  In the mean time, you "
!    "can run a second pass on the output of the first pass, using lower drag. ");
  
  int
--- 21,36 ----
  #include "privateTen.h"
  
! #define INFO "Calculate balanced gradient directions for DWI acquisition"
  char *_tend_gradsInfoL =
    (INFO
     ", based on a simulation of anti-podal point pairs repelling each other "
!    "on the unit sphere surface. This can either distribute more uniformly "
!    "a given set of gradients, or it can make a new distribution from scratch. "
!    "A more clever implementation could decrease drag with time, as the "
!    "solution converges, to get closer to the minimum energy configuration "
!    "faster.  In the mean time, you can run a second pass on the output of "
!    "the first pass, using lower drag. A second phase of the algorithm "
!    "tries random signs in gradient directions in trying to find an optimally "
!    "balanced set of directions.");
  
  int
***************
*** 57,63 ****
    hestOptAdd(&hopt, "dt", "dt", airTypeDouble, 1, 1, &(tgparm->dt), "0.05",
  	     "time increment in solver");
!   hestOptAdd(&hopt, "drag", "drag", airTypeDouble, 1, 1, &(tgparm->drag), "0.01",
! 	     "viscous drag, to keep things stable");
!   hestOptAdd(&hopt, "charge", "charge", airTypeDouble, 1, 1, &(tgparm->charge), "0.2",
  	     "amount of charge on each particle");
    hestOptAdd(&hopt, "single", NULL, airTypeInt, 0, 0, &(tgparm->single), NULL,
--- 60,67 ----
    hestOptAdd(&hopt, "dt", "dt", airTypeDouble, 1, 1, &(tgparm->dt), "0.05",
  	     "time increment in solver");
!   hestOptAdd(&hopt, "drag", "drag", airTypeDouble, 1, 1, &(tgparm->drag),
! 	     "0.0005", "viscous drag, to keep things stable");
!   hestOptAdd(&hopt, "charge", "charge", airTypeDouble, 1, 1,
! 	     &(tgparm->charge), "0.2",
  	     "amount of charge on each particle");
    hestOptAdd(&hopt, "single", NULL, airTypeInt, 0, 0, &(tgparm->single), NULL,
***************
*** 69,77 ****
  	     "positions should be saved out.  By default (not using this "
  	     "option), there is no such snapshot behavior");
!   hestOptAdd(&hopt, "jitter", "jitter", airTypeDouble, 1, 1, &(tgparm->jitter), "0.05",
  	     "amount by which to perturb points when given an input nrrd");
    hestOptAdd(&hopt, "maxiter", "# iters", airTypeInt, 1, 1,
  	     &(tgparm->maxIteration), "1000000",
  	     "max number of iterations for which to run the simulation");
    hestOptAdd(&hopt, "o", "filename", airTypeString, 1, 1, &outS, "-",
  	     "file to write output nrrd to");
--- 73,97 ----
  	     "positions should be saved out.  By default (not using this "
  	     "option), there is no such snapshot behavior");
!   hestOptAdd(&hopt, "jitter", "jitter", airTypeDouble, 1, 1,
! 	     &(tgparm->jitter), "0.05",
  	     "amount by which to perturb points when given an input nrrd");
+   hestOptAdd(&hopt, "minvelo", "vel", airTypeDouble, 1, 1, 
+ 	     &(tgparm->minVelocity), "0.00001",
+ 	     "low threshold on mean velocity of repelling points, "
+ 	     "at which point repulsion phase of algorithm terminates. ");
    hestOptAdd(&hopt, "maxiter", "# iters", airTypeInt, 1, 1,
  	     &(tgparm->maxIteration), "1000000",
  	     "max number of iterations for which to run the simulation");
+   hestOptAdd(&hopt, "minimprov", "delta", airTypeDouble, 1, 1, 
+ 	     &(tgparm->minMeanImprovement), "0.00001",
+ 	     "in the second phase of the algorithm, "
+ 	     "when stochastically optimizing the balance of gradients, "
+ 	     "the (small) improvement in length of mean gradient "
+ 	     "which triggers termination (as further improvements "
+ 	     "are unlikely. ");
+   hestOptAdd(&hopt, "minmean", "len", airTypeDouble, 1, 1,
+ 	     &(tgparm->minMean), "0.0001",
+ 	     "if length of mean gradient falls below this, finish "
+ 	     "the balancing phase");
    hestOptAdd(&hopt, "o", "filename", airTypeString, 1, 1, &outS, "-",
  	     "file to write output nrrd to");