Menu
▾
▴
[Teem-commit] teem/src/gage GNUmakefile,1.13,1.14 gage.h,1.58,1.59 miscGage.c,1.9,1.10 print.c,1.14,1.15 pvl.c,1.9,1.10 scl.c,1.29,1.30 sclanswer.c,1.14,1.15 sclfilter.c,1.14,1.15 st.c,1.5,1.6 update.c,1.11,1.12 vecGage.c,1.13,1.14
Update of /cvsroot/teem/teem/src/gage In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv29091 Modified Files: GNUmakefile gage.h miscGage.c print.c pvl.c scl.c sclanswer.c sclfilter.c st.c update.c vecGage.c Log Message: massive internal changes, but hopefully only the followin are visible: API CHANGE: unsigned int queries are gone; now using gageQuery, a bit-vector based on an array of uchars. So, where you used to say: unsigned int query; ... gageQuerySet(ctx, pvl, 1 << gageSclValue | 1 << gageSclGradMag) you'll now say: gageQuery query; ... GAGE_QUERY_RESET(query); GAGE_QUERY_ITEM_ON(query, gageSclValue); GAGE_QUERY_ITEM_ON(query, gageSclGradMag); gageQuerySet(ctx, pvl, query); or, less annoyingly (perhaps): gageQueryReset(ctx, pvl); gageQueryItemOn(ctx, pvl, gageSclValue); gageQueryItemOn(ctx, pvl, gageSclGradMag); API CHANGE: many of the arrays in the gageKind have been replaced by a single array of things called gageItemEntry, which is must easier to maintain by hand. See gage.h for details. Index: GNUmakefile =================================================================== RCS file: /cvsroot/teem/teem/src/gage/GNUmakefile,v retrieving revision 1.13 retrieving revision 1.14 diff -C2 -d -r1.13 -r1.14 *** GNUmakefile 7 Jan 2004 15:34:29 -0000 1.13 --- GNUmakefile 13 Feb 2004 22:56:15 -0000 1.14 *************** *** 47,55 **** $(L).PUBLIC_HEADERS = gage.h $(L).PRIVATE_HEADERS = privateGage.h ! $(L).OBJS = defaultsGage.o miscGage.o print.o filter.o update.o \ ! scl.o sclfilter.o sclprint.o sclanswer.o \ ! vecGage.o vecprint.o ! $(L).OBJS = defaultsGage.o miscGage.o shape.o pvl.o ctx.o update.o filter.o \ ! print.o scl.o sclanswer.o sclprint.o sclfilter.o \ vecGage.o vecprint.o st.o $(L).TESTS = --- 47,53 ---- $(L).PUBLIC_HEADERS = gage.h $(L).PRIVATE_HEADERS = privateGage.h ! $(L).OBJS = defaultsGage.o miscGage.o scl.o kind.o \ ! shape.o pvl.o ctx.o update.o filter.o \ ! print.o sclanswer.o sclprint.o sclfilter.o \ vecGage.o vecprint.o st.o $(L).TESTS = Index: gage.h =================================================================== RCS file: /cvsroot/teem/teem/src/gage/gage.h,v retrieving revision 1.58 retrieving revision 1.59 diff -C2 -d -r1.58 -r1.59 *** gage.h 16 Jan 2004 17:35:12 -0000 1.58 --- gage.h 13 Feb 2004 22:56:15 -0000 1.59 *************** *** 55,58 **** --- 55,65 ---- /* + ** terminology: gage is intended to measure multiple things at one + ** point in a volume. The SET of ALL the things that you want + ** gage to measure is called the "QUERY". Each of the many quantities + ** comprising the query are called "ITEM"s. + */ + + /* ******** gage_t ** *************** *** 181,188 **** #define GAGE_KERNEL_NUM 6 /* ** modifying the enums below (scalar, vector, etc query quantities) ! ** necesitates modifying the associated arrays in arrays.c, the arrays ! ** in enums.c, and obviously the "answer" method itself. */ --- 188,226 ---- #define GAGE_KERNEL_NUM 6 + #define GAGE_ITEM_PREREQ_NUM 5 + /* + ******** gageItemEntry struct + ** + ** necessary information about each item supported by the kind, which + ** is defined at compile-time in a per-kind table (at least it is for + ** the scalar, vector, and tensor kinds) + */ + typedef struct { + int enumVal, /* the item's enum value */ + answerLength, /* how many gage_t's are needed to store the answer */ + needDeriv, /* what kind of derivative info is immediately needed + for this item (not recursively expanded). This is + NO LONGER a bitvector: values are 0, 1, 2, ... */ + prereq[GAGE_ITEM_PREREQ_NUM], + /* what are the other items this item depends on + (not recusively expanded). The move away from + bit-vectors based on primitive types is what + necessitates conveying this with a stupid little + array. */ + parentItem, /* the enum value of an item (answerLength > 1) + containing the smaller value for which we are + merely an alias + OR: -1 if there's no parent */ + parentIndex; /* where our value starts in parents value + OR: -1 if there's no parent */ + } gageItemEntry; + /* ** modifying the enums below (scalar, vector, etc query quantities) ! ** necesitates modifying: ! ** - the central item table ! ** - the associated arrays in arrays.c ! ** - the arrays in enums.c, ! ** - the "answer" method itself. */ *************** *** 190,200 **** ******** gageScl* enum ** ! ** all the things that gage can measure in a scalar volume. The query is ! ** formed by a bitwise-or of left-shifts of 1 by these values: ! ** (1<<gageSclValue)|(1<<gageSclGradMag)|(1<<gageScl2ndDD) ! ** queries for the value, gradient magnitude, and 2nd directional derivative. ! ** The un-bit-shifted values are required for gage to index arrays like ! ** gageSclAnsOffset[], _gageSclNeedDeriv[], _gageSclPrereq[], etc, and ! ** for the gageScl airEnum. ** ** NOTE: although it is currently listed that way, it is not necessary --- 228,232 ---- ******** gageScl* enum ** ! ** all the "items" that gage can measure in a scalar volume. ** ** NOTE: although it is currently listed that way, it is not necessary *************** *** 209,270 **** enum { gageSclUnknown=-1, /* -1: nobody knows */ ! gageSclValue, /* 0: "v", data value: *GT */ gageSclGradVec, /* 1: "grad", gradient vector, un-normalized: GT[3] */ ! gageSclGradMag, /* 2: "gm", gradient magnitude: *GT */ gageSclNormal, /* 3: "n", gradient vector, normalized: GT[3] */ gageSclNPerp, /* 4: "np", projection onto tangent plane: GT[9] */ gageSclHessian, /* 5: "h", Hessian: GT[9] (column-order) */ ! gageSclLaplacian, /* 6: "l", Laplacian: Dxx + Dyy + Dzz: *GT */ gageSclHessEval, /* 7: "heval", Hessian's eigenvalues: GT[3] */ ! gageSclHessEvec, /* 8: "hevec", Hessian's eigenvectors: GT[9] */ ! gageScl2ndDD, /* 9: "2d", 2nd dir.deriv. along gradient: *GT */ ! gageSclGeomTens, /* 10: "gten", sym. matx w/ evals 0,K1,K2 and evecs grad, ! curvature directions: GT[9] */ ! gageSclK1, /* 11: "k1", 1st principle curvature: *GT */ ! gageSclK2, /* 12: "k2", 2nd principle curvature (k2 <= k1): *GT */ ! gageSclCurvedness, /* 13: "cv", L2 norm of K1, K2 (not Koen.'s "C"): *GT */ ! gageSclShapeTrace, /* 14, "st", (K1+K2)/Curvedness: *GT */ ! gageSclShapeIndex, /* 15: "si", Koen.'s shape index, ("S"): *GT */ ! gageSclMeanCurv, /* 16: "mc", mean curvature (K1 + K2)/2: *GT */ ! gageSclGaussCurv, /* 17: "gc", gaussian curvature K1*K2: *GT */ ! gageSclCurvDir, /* 18: "cdir", principle curvature directions: GT[6] */ ! gageSclFlowlineCurv,/* 19: "fc", curvature of normal streamline: *GT */ gageSclLast }; ! #define GAGE_SCL_MAX 19 ! #define GAGE_SCL_TOTAL_ANS_LENGTH 63 ! ! /* ! ******** GAGE_SCL_*_BIT #defines ! ** already bit-shifted for you, so that query: ! ** (1<<gageSclValue)|(1<<gageSclGradMag)|(1<<gageScl2ndDD) ! ** is same as: ! ** GAGE_SCL_VALUE_BIT | GAGE_SCL_GRADMAG_BIT | GAGE_SCL_2NDDD_BIT ! */ ! #define GAGE_SCL_VALUE_BIT (1<<0) ! #define GAGE_SCL_GRADVEC_BIT (1<<1) ! #define GAGE_SCL_GRADMAG_BIT (1<<2) ! #define GAGE_SCL_NORMAL_BIT (1<<3) ! #define GAGE_SCL_NPERP_BIT (1<<4) ! #define GAGE_SCL_HESSIAN_BIT (1<<5) ! #define GAGE_SCL_LAPLACIAN_BIT (1<<6) ! #define GAGE_SCL_HESSEVAL_BIT (1<<7) ! #define GAGE_SCL_HESSEVEC_BIT (1<<8) ! #define GAGE_SCL_2NDDD_BIT (1<<9) ! #define GAGE_SCL_GEOMTENS_BIT (1<<10) ! #define GAGE_SCL_K1_BIT (1<<11) ! #define GAGE_SCL_K2_BIT (1<<12) ! #define GAGE_SCL_CURVEDNESS_BIT (1<<13) ! #define GAGE_SCL_SHAPETRACE_BIT (1<<14) ! #define GAGE_SCL_SHAPEINDEX_BIT (1<<15) ! #define GAGE_SCL_MEANCURV_BIT (1<<16) ! #define GAGE_SCL_GAUSSCURV_BIT (1<<17) ! #define GAGE_SCL_CURVDIR_BIT (1<<18) ! #define GAGE_SCL_FLOWLINECURV_BIT (1<<19) /* ******** gageVec* enum ** ! ** all the things that gage knows how to measure in a 3-vector volume ** ** The strings gives one of the gageVec airEnum identifiers, and GT[x] --- 241,280 ---- enum { gageSclUnknown=-1, /* -1: nobody knows */ ! gageSclValue, /* 0: "v", data value: GT[1] */ gageSclGradVec, /* 1: "grad", gradient vector, un-normalized: GT[3] */ ! gageSclGradMag, /* 2: "gm", gradient magnitude: GT[1] */ gageSclNormal, /* 3: "n", gradient vector, normalized: GT[3] */ gageSclNPerp, /* 4: "np", projection onto tangent plane: GT[9] */ 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] */ ! gageSclCurvedness, /* 19: "cv", 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 /* ******** gageVec* enum ** ! ** all the "items" that gage knows how to measure in a 3-vector volume ** ** The strings gives one of the gageVec airEnum identifiers, and GT[x] *************** *** 274,278 **** gageVecUnknown=-1, /* -1: nobody knows */ gageVecVector, /* 0: "v", component-wise-intrpolated (CWI) vec: GT[3] */ ! gageVecLength, /* 1: "l", length of CWI vector: *GT */ gageVecNormalized, /* 2: "n", normalized CWI vector: GT[3] */ gageVecJacobian, /* 3: "j", component-wise Jacobian: GT[9] --- 284,288 ---- gageVecUnknown=-1, /* -1: nobody knows */ gageVecVector, /* 0: "v", component-wise-intrpolated (CWI) vec: GT[3] */ ! gageVecLength, /* 1: "l", length of CWI vector: GT[1] */ gageVecNormalized, /* 2: "n", normalized CWI vector: GT[3] */ gageVecJacobian, /* 3: "j", component-wise Jacobian: GT[9] *************** *** 280,310 **** 1:dv_y/dx 4:dv_y/dy 7:dv_y/dz 2:dv_z/dx 5:dv_z/dy 8:dv_z/dz */ ! gageVecDivergence, /* 4: "d", divergence (based on Jacobian): *GT */ gageVecCurl, /* 5: "c", curl (based on Jacobian): GT[3] */ ! gageVecGradient0, /* 6: "g1", gradient of 1st component of vector: GT[3] */ ! gageVecGradient1, /* 7: "g2", gradient of 2nd component of vector: GT[3] */ ! gageVecGradient2, /* 8: "g3", gradient of 3rd component of vector: GT[3] */ gageVecMultiGrad, /* 9: "mg", sum of outer products of gradients: GT[9] */ ! gageVecMGFrob, /* 10: "mgfrob", frob norm of multi-gradient: *GT */ gageVecMGEval, /* 11: "mgeval", eigenvalues of multi-gradient: GT[3] */ gageVecMGEvec, /* 12: "mgevec", eigenvectors of multi-gradient: GT[9] */ gageVecLast }; ! #define GAGE_VEC_MAX 12 ! #define GAGE_VEC_TOTAL_ANS_LENGTH 51 ! ! #define GAGE_VEC_VECTOR_BIT (1<<0) ! #define GAGE_VEC_LENGTH_BIT (1<<1) ! #define GAGE_VEC_NORMALIZED_BIT (1<<2) ! #define GAGE_VEC_JACOBIAN_BIT (1<<3) ! #define GAGE_VEC_DIVERGENCE_BIT (1<<4) ! #define GAGE_VEC_CURL_BIT (1<<5) ! #define GAGE_VEC_GRADIENT0_BIT (1<<6) ! #define GAGE_VEC_GRADIENT1_BIT (1<<7) ! #define GAGE_VEC_GRADIENT2_BIT (1<<8) ! #define GAGE_VEC_MULTIGRAD_BIT (1<<9) ! #define GAGE_VEC_MGFROB_BIT (1<<10) ! #define GAGE_VEC_MGEVAL_BIT (1<<11) ! #define GAGE_VEC_MGEVEC_BIT (1<<12) struct gageKind_t; /* dumb forward declaraction, ignore */ --- 290,305 ---- 1:dv_y/dx 4:dv_y/dy 7:dv_y/dz 2:dv_z/dx 5:dv_z/dy 8:dv_z/dz */ ! gageVecDivergence, /* 4: "d", divergence (based on Jacobian): GT[1] */ gageVecCurl, /* 5: "c", curl (based on Jacobian): GT[3] */ ! gageVecGradient0, /* 6: "g0", gradient of 1st component of vector: GT[3] */ ! gageVecGradient1, /* 7: "g1", gradient of 2nd component of vector: GT[3] */ ! gageVecGradient2, /* 8: "g2", gradient of 3rd component of vector: GT[3] */ gageVecMultiGrad, /* 9: "mg", sum of outer products of gradients: GT[9] */ ! gageVecMGFrob, /* 10: "mgfrob", frob norm of multi-gradient: GT[1] */ gageVecMGEval, /* 11: "mgeval", eigenvalues of multi-gradient: GT[3] */ gageVecMGEvec, /* 12: "mgevec", eigenvectors of multi-gradient: GT[9] */ gageVecLast }; ! #define GAGE_VEC_ITEM_MAX 12 struct gageKind_t; /* dumb forward declaraction, ignore */ *************** *** 423,426 **** --- 418,460 ---- /* + ******** gageQuery typedef + ** + ** this short, fixed-length array is used as a bit-vector to store + ** all the items that comprise a query. Its length sets an upper + ** bound on the maximum item value that a gageKind may use. + ** + ** The important thing to keep in mind is that a variable of type + ** gageKind ends up being passed by reference, even though the + ** syntax of its usage doesn't immediately announce that. + ** + ** gageKindCheck currently has the role of verifying that a gageKind's + ** maximum item value is within the bounds set here. Using + ** GAGE_QUERY_BYTES_NUM == 8 gives a max item value of 63, which is + ** far above anything being used now. + */ + #define GAGE_QUERY_BYTES_NUM 8 + #define GAGE_ITEM_MAX ((8*GAGE_QUERY_BYTES_NUM)-1) + typedef unsigned char gageQuery[GAGE_QUERY_BYTES_NUM]; + + /* + ******** GAGE_QUERY_RESET, GAGE_QUERY_TEST + ******** GAGE_QUERY_ON, GAGE_QUERY_OFF + ** + ** Macros for manipulating a gageQuery. + ** + ** airSanity ensures that an unsigned char is in fact 8 bits + */ + #define GAGE_QUERY_RESET(q) (q[0]=q[1]=q[2]=q[3]=q[4]=q[5]=q[6]=q[7]=0) + #define GAGE_QUERY_COPY(p, q) \ + p[0] = q[0]; p[1] = q[1]; p[2] = q[2]; p[3] = q[3]; \ + p[4] = q[4]; p[5] = q[5]; p[6] = q[6]; p[7] = q[7] + #define GAGE_QUERY_EQUAL(p, q) \ + (p[0] == q[0] && p[1] == q[1] && p[2] == q[2] && p[3] == q[3] \ + && p[4] == q[4] && p[5] == q[5] && p[6] == q[6] && p[7] == q[7]) + #define GAGE_QUERY_ITEM_TEST(q, i) (q[i/8] & (1 << (i % 8))) + #define GAGE_QUERY_ITEM_ON(q, i) (q[i/8] |= (1 << (i % 8))) + #define GAGE_QUERY_ITEM_OFF(q, i) (q[i/8] &= ^(1 << (i % 8))) + + /* ******** gageContext struct ** *************** *** 482,486 **** thisIsACopy; /* I'm a copy */ struct gageKind_t *kind; /* what kind of volume is this pervolume for */ ! unsigned int query; /* the query, recursively expanded */ int needD[3]; /* which derivatives need to be calculated for the query (above) in this volume */ --- 516,520 ---- thisIsACopy; /* I'm a copy */ struct gageKind_t *kind; /* what kind of volume is this pervolume for */ ! gageQuery query; /* the query, recursively expanded */ int needD[3]; /* which derivatives need to be calculated for the query (above) in this volume */ *************** *** 508,513 **** /* nrrd{F,D}Lookup[] element, according to npad->type and gage_t */ ! gage_t *ans; /* array of length kind->totalAnsLen, to hold ! all answers */ } gagePerVolume; --- 542,547 ---- /* nrrd{F,D}Lookup[] element, according to npad->type and gage_t */ ! gage_t *answer; /* main buffer to hold all the answers */ ! gage_t **directAnswer; /* array of pointers into answer */ } gagePerVolume; *************** *** 526,536 **** (0 for scalars, 1 for vectors) */ valLen, /* number of scalars per data point */ ! queryMax, /* such as GAGE_SCL_MAX */ ! *ansLength, /* such as gageSclAnsLength */ ! *ansOffset, /* such as gageSclAnsOffset */ ! totalAnsLen, /* such as GAGE_SCL_TOTAL_ANS_LENGTH */ ! *needDeriv; /* such as _gageSclNeedDeriv */ ! unsigned int *queryPrereq; /* such as _gageSclPrereq; */ ! void (*iv3Print)(FILE *, /* such as _gageSclIv3Print() */ gageContext *, --- 560,566 ---- (0 for scalars, 1 for vectors) */ valLen, /* number of scalars per data point */ ! itemMax; /* such as GAGE_SCL_ITEM_MAX */ ! gageItemEntry *table; /* array of gageItemEntry's, indexed ! by the item value */ void (*iv3Print)(FILE *, /* such as _gageSclIv3Print() */ gageContext *, *************** *** 542,556 **** } gageKind; - - /* the "answer structs" are gone- now we're just using the bare array - "ans" in the gagePerVolume, more easily used with help from - gageAnswerPointer() or ... */ - /* ! ******** #define GAGE_ANSWER_POINTER() ** ! ** a non-error-checking version of gageAnswerPointer() */ ! #define GAGE_ANSWER_POINTER(pvl, m) ((pvl)->ans + (pvl)->kind->ansOffset[(m)]) /* defaultsGage.c */ --- 572,584 ---- } gageKind; /* ! ******** gageItemSpec struct ** ! ** dumb little way to store an item/kind pair */ ! typedef struct { ! int item; /* the quantity we care about */ ! gageKind *kind; /* what kind it comes from */ ! } gageItemSpec; /* defaultsGage.c */ *************** *** 580,586 **** extern void gageParmReset(gageParm *parm); extern void gagePointReset(gagePoint *point); /* print.c */ ! extern void gageQueryPrint(FILE *file, gageKind *kind, unsigned int query); /* sclfilter.c */ --- 608,621 ---- extern void gageParmReset(gageParm *parm); extern void gagePointReset(gagePoint *point); + extern gageItemSpec *gageItemSpecNew(void); + extern gageItemSpec *gageItemSpecNix(gageItemSpec *isp); + + /* kind.c */ + extern int gageKindCheck(gageKind *kind); + extern int gageKindTotalAnswerLength(gageKind *kind); + extern int gageKindAnswerOffset(gageKind *kind, int item); /* print.c */ ! extern void gageQueryPrint(FILE *file, gageKind *kind, gageQuery query); /* sclfilter.c */ *************** *** 600,605 **** /* scl.c */ - extern gage_export int gageSclAnsLength[GAGE_SCL_MAX+1]; - extern gage_export int gageSclAnsOffset[GAGE_SCL_MAX+1]; extern gage_export airEnum *gageScl; extern gage_export gageKind *gageKindScl; --- 635,638 ---- *************** *** 608,613 **** implement the "vec" kind, and could be used as examples of what it takes to create a new gageKind) */ - extern gage_export int gageVecAnsLength[GAGE_VEC_MAX+1]; - extern gage_export int gageVecAnsOffset[GAGE_VEC_MAX+1]; extern gage_export airEnum *gageVec; extern gage_export gageKind *gageKindVec; --- 641,644 ---- *************** *** 638,644 **** gagePerVolume *pvl, gageNixer_t *nixer); extern gage_t *gageAnswerPointer(gageContext *ctx, ! gagePerVolume *pvl, int measure); ! extern int gageQuerySet(gageContext *ctx, ! gagePerVolume *pvl, unsigned int query); /* ctx.c */ --- 669,676 ---- gagePerVolume *pvl, gageNixer_t *nixer); extern gage_t *gageAnswerPointer(gageContext *ctx, ! gagePerVolume *pvl, int item); ! extern int gageQueryReset(gageContext *ctx, gagePerVolume *pvl); ! extern int gageQuerySet(gageContext *ctx, gagePerVolume *pvl, gageQuery query); ! extern int gageQueryItemOn(gageContext *ctx, gagePerVolume *pvl, int item); /* ctx.c */ Index: miscGage.c =================================================================== RCS file: /cvsroot/teem/teem/src/gage/miscGage.c,v retrieving revision 1.9 retrieving revision 1.10 diff -C2 -d -r1.9 -r1.10 *** miscGage.c 7 Jan 2004 15:34:29 -0000 1.9 --- miscGage.c 13 Feb 2004 22:56:15 -0000 1.10 *************** *** 214,215 **** --- 214,236 ---- return; } + + gageItemSpec * + gageItemSpecNew(void) { + gageItemSpec *isp; + + isp = (gageItemSpec *)calloc(1, sizeof(gageItemSpec)); + if (isp) { + isp->item = -1; + isp->kind = NULL; + } + return isp; + } + + gageItemSpec * + gageItemSpecNix(gageItemSpec *isp) { + + if (isp) { + AIR_FREE(isp); + } + return NULL; + } Index: print.c =================================================================== RCS file: /cvsroot/teem/teem/src/gage/print.c,v retrieving revision 1.14 retrieving revision 1.15 diff -C2 -d -r1.14 -r1.15 *** print.c 7 Jan 2004 15:34:29 -0000 1.14 --- print.c 13 Feb 2004 22:56:15 -0000 1.15 *************** *** 165,178 **** void ! gageQueryPrint (FILE *file, gageKind *kind, unsigned int query) { ! unsigned int q; ! fprintf(file, "%s query = %u ...\n", kind->name, query); ! q = kind->queryMax+1; do { ! q--; ! if ((1<<q) & query) { ! fprintf(file, " %3d: %s\n", q, airEnumStr(kind->enm, q)); } ! } while (q); } --- 165,178 ---- void ! gageQueryPrint (FILE *file, gageKind *kind, gageQuery query) { ! int ii; ! fprintf(file, "%s query = ...\n", kind->name); ! ii = kind->itemMax+1; do { ! ii--; ! if (GAGE_QUERY_ITEM_TEST(query, ii)) { ! fprintf(file, " %3d: %s\n", ii, airEnumStr(kind->enm, ii)); } ! } while (ii); } Index: pvl.c =================================================================== RCS file: /cvsroot/teem/teem/src/gage/pvl.c,v retrieving revision 1.9 retrieving revision 1.10 diff -C2 -d -r1.9 -r1.10 *** pvl.c 7 Jan 2004 15:34:29 -0000 1.9 --- pvl.c 13 Feb 2004 22:56:15 -0000 1.10 *************** *** 52,56 **** char me[]="gagePerVolumeNew", err[AIR_STRLEN_MED]; gagePerVolume *pvl; ! int i; if (!( nin && kind )) { --- 52,56 ---- char me[]="gagePerVolumeNew", err[AIR_STRLEN_MED]; gagePerVolume *pvl; ! int ii; if (!( nin && kind )) { *************** *** 70,74 **** pvl->verbose = gageDefVerbose; pvl->kind = kind; ! pvl->query = 0; pvl->needD[0] = pvl->needD[1] = pvl->needD[2] = AIR_FALSE; pvl->nin = nin; --- 70,74 ---- pvl->verbose = gageDefVerbose; pvl->kind = kind; ! GAGE_QUERY_RESET(pvl->query); pvl->needD[0] = pvl->needD[1] = pvl->needD[2] = AIR_FALSE; pvl->nin = nin; *************** *** 77,90 **** pvl->padInfo = NULL; pvl->npad = NULL; ! for (i=0; i<GAGE_PVL_FLAG_NUM; i++) { ! pvl->flag[i] = AIR_FALSE; } pvl->iv3 = pvl->iv2 = pvl->iv1 = NULL; pvl->lup = nrrdLOOKUP[nin->type]; ! pvl->ans = (gage_t *)calloc(kind->totalAnsLen, sizeof(gage_t)); ! if (!pvl->ans) { ! sprintf(err, "%s: couldn't alloc answer array", me); biffAdd(GAGE, err); return NULL; } pvl->flag[gagePvlFlagVolume] = AIR_TRUE; --- 77,95 ---- pvl->padInfo = NULL; pvl->npad = NULL; ! for (ii=0; ii<GAGE_PVL_FLAG_NUM; ii++) { ! pvl->flag[ii] = AIR_FALSE; } pvl->iv3 = pvl->iv2 = pvl->iv1 = NULL; pvl->lup = nrrdLOOKUP[nin->type]; ! pvl->answer = (gage_t *)calloc(gageKindTotalAnswerLength(kind), ! sizeof(gage_t)); ! pvl->directAnswer = (gage_t **)calloc(kind->itemMax+1, sizeof(gage_t*)); ! if (!(pvl->answer && pvl->directAnswer)) { ! sprintf(err, "%s: couldn't alloc answer and directAnswer arrays", me); biffAdd(GAGE, err); return NULL; } + for (ii=0; ii<=kind->itemMax; ii++) { + pvl->directAnswer[ii] = pvl->answer + gageKindAnswerOffset(kind, ii); + } pvl->flag[gagePvlFlagVolume] = AIR_TRUE; *************** *** 102,105 **** --- 107,111 ---- char me[]="gagePerVolumeCopy", err[AIR_STRLEN_MED]; gagePerVolume *nvl; + int ii; nvl = (gagePerVolume *)calloc(1, sizeof(gagePerVolume)); *************** *** 122,130 **** nvl->iv2 = (gage_t *)calloc(fd*fd*nvl->kind->valLen, sizeof(gage_t)); nvl->iv1 = (gage_t *)calloc(fd*nvl->kind->valLen, sizeof(gage_t)); ! nvl->ans = (gage_t *)calloc(nvl->kind->totalAnsLen, sizeof(gage_t)); ! if (!( nvl->iv3 && nvl->iv2 && nvl->iv1 && nvl->ans )) { sprintf(err, "%s: couldn't allocate all caches", me); biffAdd(GAGE, err); return NULL; } return nvl; --- 128,143 ---- nvl->iv2 = (gage_t *)calloc(fd*fd*nvl->kind->valLen, sizeof(gage_t)); nvl->iv1 = (gage_t *)calloc(fd*nvl->kind->valLen, sizeof(gage_t)); ! nvl->answer = (gage_t *)calloc(gageKindTotalAnswerLength(nvl->kind), ! sizeof(gage_t)); ! nvl->directAnswer = (gage_t **)calloc(nvl->kind->itemMax+1, ! sizeof(gage_t*)); ! if (!( nvl->iv3 && nvl->iv2 && nvl->iv1 ! && nvl->answer && nvl->directAnswer )) { sprintf(err, "%s: couldn't allocate all caches", me); biffAdd(GAGE, err); return NULL; } + for (ii=0; ii<=pvl->kind->itemMax; ii++) { + nvl->directAnswer[ii] = nvl->answer + gageKindAnswerOffset(pvl->kind, ii); + } return nvl; *************** *** 148,152 **** pvl->nixer(pvl->npad, pvl->kind, pvl); } ! AIR_FREE(pvl->ans); AIR_FREE(pvl); return NULL; --- 161,166 ---- pvl->nixer(pvl->npad, pvl->kind, pvl); } ! AIR_FREE(pvl->answer); ! AIR_FREE(pvl->directAnswer); AIR_FREE(pvl); return NULL; *************** *** 183,194 **** ** way of getting a pointer to a specific answer in a pervolume's ans array ** - ** Basically just a wrapper around GAGE_ANSWER_POINTER with error checking */ gage_t * ! gageAnswerPointer (gageContext *ctx, gagePerVolume *pvl, int measure) { gage_t *ret; ! if (pvl && !airEnumValCheck(pvl->kind->enm, measure)) { ! ret = GAGE_ANSWER_POINTER(pvl, measure); } else { ret = NULL; --- 197,207 ---- ** way of getting a pointer to a specific answer in a pervolume's ans array ** */ gage_t * ! gageAnswerPointer (gageContext *ctx, gagePerVolume *pvl, int item) { gage_t *ret; ! if (pvl && !airEnumValCheck(pvl->kind->enm, item)) { ! ret = pvl->answer + gageKindAnswerOffset(pvl->kind, item); } else { ret = NULL; *************** *** 197,200 **** --- 210,232 ---- } + int + gageQueryReset(gageContext *ctx, gagePerVolume *pvl) { + char me[]="gageQueryReset", err[AIR_STRLEN_MED]; + + if (!( pvl )) { + sprintf(err, "%s: got NULL pointer", me); + biffAdd(GAGE, err); return 1; + } + if (pvl->thisIsACopy) { + sprintf(err, "%s: can't operate on a pervolume copy", me); + biffAdd(GAGE, err); return 1; + } + + GAGE_QUERY_RESET(pvl->query); + + return 0; + } + + /* ******** gageQuerySet() *************** *** 204,212 **** ** ** Sets: pvl->query */ int ! gageQuerySet (gageContext *ctx, gagePerVolume *pvl, unsigned int query) { char me[]="gageQuerySet", err[AIR_STRLEN_MED]; ! unsigned int mask, lastq, q; if (!( pvl )) { --- 236,248 ---- ** ** Sets: pvl->query + ** + ** the gageContext is not actually used here, but I'm cautiously + ** including it in case its used in the future. */ int ! gageQuerySet (gageContext *ctx, gagePerVolume *pvl, gageQuery query) { char me[]="gageQuerySet", err[AIR_STRLEN_MED]; ! gageQuery lastQuery; ! int pi, ii; if (!( pvl )) { *************** *** 218,227 **** biffAdd(GAGE, err); return 1; } ! mask = (1U << (pvl->kind->queryMax+1)) - 1; ! if (query != (query & mask)) { ! sprintf(err, "%s: invalid bits set in query", me); ! biffAdd(GAGE, err); return 1; ! } ! pvl->query = query; if (pvl->verbose) { fprintf(stderr, "%s: original ", me); --- 254,258 ---- biffAdd(GAGE, err); return 1; } ! GAGE_QUERY_COPY(pvl->query, query); if (pvl->verbose) { fprintf(stderr, "%s: original ", me); *************** *** 230,241 **** /* recursive expansion of prerequisites */ do { ! lastq = pvl->query; ! q = pvl->kind->queryMax+1; do { ! q--; ! if ((1<<q) & pvl->query) ! pvl->query |= pvl->kind->queryPrereq[q]; ! } while (q); ! } while (pvl->query != lastq); if (pvl->verbose) { fprintf(stderr, "%s: expanded ", me); --- 261,277 ---- /* recursive expansion of prerequisites */ do { ! GAGE_QUERY_COPY(lastQuery, pvl->query); ! ii = pvl->kind->itemMax+1; do { ! ii--; ! if (GAGE_QUERY_ITEM_TEST(pvl->query, ii)) { ! for (pi=0; pi<GAGE_ITEM_PREREQ_NUM; pi++) { ! if (-1 != pvl->kind->table[ii].prereq[pi]) { ! GAGE_QUERY_ITEM_ON(pvl->query, pvl->kind->table[ii].prereq[pi]); ! } ! } ! } ! } while (ii); ! } while (!GAGE_QUERY_EQUAL(pvl->query, lastQuery)); if (pvl->verbose) { fprintf(stderr, "%s: expanded ", me); *************** *** 246,247 **** --- 282,306 ---- return 0; } + + int + gageQueryItemOn(gageContext *ctx, gagePerVolume *pvl, int item) { + char me[]="gageQueryItemOn", err[AIR_STRLEN_MED]; + + if (!( pvl )) { + sprintf(err, "%s: got NULL pointer", me); + biffAdd(GAGE, err); return 1; + } + if (pvl->thisIsACopy) { + sprintf(err, "%s: can't operate on a pervolume copy", me); + biffAdd(GAGE, err); return 1; + } + + GAGE_QUERY_ITEM_ON(pvl->query, item); + if (gageQuerySet(ctx, pvl, pvl->query)) { + sprintf(err, "%s: trouble", me); + biffAdd(GAGE, err); return 1; + } + + return 0; + } + Index: scl.c =================================================================== RCS file: /cvsroot/teem/teem/src/gage/scl.c,v retrieving revision 1.29 retrieving revision 1.30 diff -C2 -d -r1.29 -r1.30 *** scl.c 7 Jan 2004 15:34:29 -0000 1.29 --- scl.c 13 Feb 2004 22:56:15 -0000 1.30 *************** *** 22,156 **** /* ! gageSclUnknown=-1, -1: nobody knows ! gageSclValue, * 0: "v", data value: *GT * ! gageSclGradVec, * 1: "grad", gradient vector, un-normalized: GT[3] * ! gageSclGradMag, * 2: "gm", gradient magnitude: *GT * ! gageSclNormal, * 3: "n", gradient vector, normalized: GT[3] * ! gageSclNPerp, * 4: "np", projection onto tangent plane: GT[9] * ! gageSclHessian, * 5: "h", Hessian: GT[9] (column-order) * ! gageSclLaplacian, * 6: "l", Laplacian: Dxx + Dyy + Dzz: *GT * ! gageSclHessEval, * 7: "heval", Hessian's eigenvalues: GT[3] * ! gageSclHessEvec, * 8: "hevec", Hessian's eigenvectors: GT[9] * ! gageScl2ndDD, * 9: "2d", 2nd dir.deriv. along gradient: *GT * ! gageSclGeomTens, * 10: "gten", sym. matx w/ evals 0,K1,K2 and evecs grad, ! curvature directions: GT[9] * ! gageSclK1, * 11: "k1", 1st principle curvature: *GT * ! gageSclK2, * 12: "k2", 2nd principle curvature (k2 <= k1): *GT * ! gageSclCurvedness, * 13: "cv", L2 norm of K1, K2 (not Koen.'s "C"): *GT * ! gageSclShapeTrace, * 14, "st", (K1+K2)/Curvedness: *GT * ! gageSclShapeIndex, * 15: "si", Koen.'s shape index, ("S"): *GT * ! gageSclMeanCurv, * 16: "mc", mean curvature (K1 + K2)/2: *GT * ! gageSclGaussCurv, * 17: "gc", gaussian curvature K1*K2: *GT * ! gageSclCurvDir, * 18: "cdir", principle curvature directions: GT[6] * ! gageSclFlowlineCurv,* 19: "nc", curvature of normal streamline: *GT * ! 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ! */ ! ! /* ! ******** gageSclAnsLength[] ! ** ! ** the number of gage_t used for each answer ! */ ! int ! gageSclAnsLength[GAGE_SCL_MAX+1] = { ! 1, 3, 1, 3, 9, 9, 1, 3, 9, 1, 9, 1, 1, 1, 1, 1, 1, 1, 6, 1 ! }; ! ! /* ! ******** gageSclAnsOffset[] ! ** ! ** the index into the answer array of the first element of the answer ! */ ! int ! gageSclAnsOffset[GAGE_SCL_MAX+1] = { ! 0, 1, 4, 5, 8, 17, 26, 27, 30, 39, 40, 49, 50, 51, 52, 53, 54, 55, 56, 62 ! /* --> 63 == GAGE_SCL_TOTAL_ANS_LENGTH */ ! }; ! ! /* ! ** _gageSclNeedDeriv[] ** ! ** each value is a BIT FLAG representing the different value/derivatives ! ** that are needed to calculate the quantity. ** ! ** 1: need value interpolation reconstruction (as with k00) ! ** 2: need first derivatives (as with k11) ! ** 4: need second derivatives (as with k22) ! */ ! int ! _gageSclNeedDeriv[GAGE_SCL_MAX+1] = { ! 1, 2, 2, 2, 2, 4, 4, 4, 4, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 ! }; ! ! /* ! ** _gageSclPrereq[] ! ** ! ** this records the measurements which are needed as ingredients for any ! ** given measurement, but it is not necessarily the recursive expansion of ! ** that requirement (that role is performed by gageQuerySet()) */ ! unsigned int ! _gageSclPrereq[GAGE_SCL_MAX+1] = { ! /* 0: gageSclValue */ ! 0, ! ! /* 1: gageSclGradVec */ ! 0, ! ! /* 2: gageSclGradMag */ ! (1<<gageSclGradVec), ! ! /* 3: gageSclNormal */ ! (1<<gageSclGradVec) | (1<<gageSclGradMag), ! ! /* 3: gageSclNPerp */ ! (1<<gageSclNormal), ! ! /* 5: gageSclHessian */ ! 0, ! ! /* 6: gageSclLaplacian */ ! (1<<gageSclHessian), /* not really true, but this is simpler */ ! ! /* 7: gageSclHessEval */ ! (1<<gageSclHessian), ! ! /* 8: gageSclHessEvec */ ! (1<<gageSclHessian) | (1<<gageSclHessEval), ! ! /* 9: gageScl2ndDD */ ! (1<<gageSclHessian) | (1<<gageSclNormal), ! ! /* 10: gageSclGeomTens */ ! (1<<gageSclHessian) | (1<<gageSclNPerp) | (1<<gageSclGradMag), ! ! /* 11: gageSclK1 */ ! (1<<gageSclCurvedness) | (1<<gageSclShapeTrace), ! ! /* 12: gageSclK2 */ ! (1<<gageSclCurvedness) | (1<<gageSclShapeTrace), ! ! /* 13: gageSclCurvedness */ ! (1<<gageSclGeomTens), ! ! /* 14: gageSclShapeTrace */ ! (1<<gageSclGeomTens), ! ! /* 15: gageSclShapeIndex */ ! (1<<gageSclK1) | (1<<gageSclK2), ! ! /* 16: gageSclMeanCurv */ ! (1<<gageSclK1) | (1<<gageSclK2), ! ! /* 17: gageSclGaussCurv */ ! (1<<gageSclK1) | (1<<gageSclK2), ! ! /* 18: gageSclCurvDir */ ! (1<<gageSclGeomTens) | (1<<gageSclK1) | (1<<gageSclK2), ! ! /* 19: gageSclFlowlineCurv */ ! /* this is because of how answer code uses sHess, nPerp, nProj */ ! (1<<gageSclGeomTens) ! }; --- 22,65 ---- /* ! ** _gageSclTable ** ! ** the static array of item information for the scalar kind. ** ! ** General notes about setting these things up: ! ** - Yes, you'll need more than 80 columns of display ! ** - You do need to explicitily initialize all the prerequisite elements, ! ** because the compiler will initialize them to zero (which is a valid ! ** item value) */ ! gageItemEntry ! _gageSclTable[GAGE_SCL_ITEM_MAX+1] = { ! /* enum value len,deriv, prereqs, parent item, index*/ ! {gageSclValue, 1, 0, {-1, -1, -1, -1, -1}, -1, -1}, ! {gageSclGradVec, 3, 1, {-1, -1, -1, -1, -1}, -1, -1}, ! {gageSclGradMag, 1, 1, {gageSclGradVec, -1, -1, -1, -1}, -1, -1}, ! {gageSclNormal, 3, 1, {gageSclGradVec, gageSclGradMag, -1, -1, -1}, -1, -1}, ! {gageSclNPerp, 9, 1, {gageSclNormal, -1, -1, -1, -1}, -1, -1}, ! {gageSclHessian, 9, 2, {gageSclHessian, -1, -1, -1, -1}, -1, -1}, ! {gageSclLaplacian, 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}, ! {gageSclHessEval1, 1, 2, {gageSclHessEval, -1, -1, -1, -1}, gageSclHessEval, 1}, ! {gageSclHessEval2, 1, 2, {gageSclHessEval, -1, -1, -1, -1}, gageSclHessEval, 2}, ! {gageSclHessEvec, 9, 2, {gageSclHessian, gageSclHessEval, -1, -1, -1}, -1, -1}, ! {gageSclHessEvec0, 3, 2, {gageSclHessEvec, -1, -1, -1, -1}, gageSclHessEvec, 0}, ! {gageSclHessEvec1, 3, 2, {gageSclHessEvec, -1, -1, -1, -1}, gageSclHessEvec, 3}, ! {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, {gageSclCurvedness, gageSclShapeTrace, -1, -1, -1}, -1, -1}, ! {gageSclK2, 1, 2, {gageSclCurvedness, gageSclShapeTrace, -1, -1, -1}, -1, -1}, ! {gageSclCurvedness, 1, 2, {gageSclGeomTens, -1, -1, -1, -1}, -1, -1}, ! {gageSclShapeTrace, 1, 2, {gageSclGeomTens, -1, -1, -1, -1}, -1, -1}, ! {gageSclShapeIndex, 1, 2, {gageSclK1, gageSclK2, -1, -1, -1}, -1, -1}, ! {gageSclMeanCurv, 1, 2, {gageSclK1, gageSclK2, -1, -1, -1}, -1, -1}, ! {gageSclGaussCurv, 1, 2, {gageSclK1, gageSclK2, -1, -1, -1}, -1, -1}, ! {gageSclCurvDir1, 3, 2, {gageSclGeomTens, gageSclK1, gageSclK2, -1, -1}, -1, -1}, ! {gageSclCurvDir2, 3, 2, {gageSclGeomTens, gageSclK1, gageSclK2, -1, -1}, -1, -1}, ! {gageSclFlowlineCurv, 1, 2, {gageSclGeomTens, -1, -1, -1, -1}, -1, -1} }; *************** *** 166,170 **** --- 75,85 ---- "Laplacian", "Hessian eigenvalues", + "Hessian eigenvalue[0]", + "Hessian eigenvalue[1]", + "Hessian eigenvalue[2]", "Hessian eigenvectors", + "Hessian eigenvector[0]", + "Hessian eigenvector[1]", + "Hessian eigenvector[2]", "2nd DD along gradient", "geometry tensor", *************** *** 176,180 **** "mean curvature", "Gaussian curvature", ! "curvature directions", "flowline curvature" }; --- 91,96 ---- "mean curvature", "Gaussian curvature", ! "1st curvature direction", ! "2nd curvature direction", "flowline curvature" }; *************** *** 191,195 **** --- 107,117 ---- "Laplacian", "Hessian's eigenvalues", + "Hessian's 1st eigenvalue", + "Hessian's 2nd eigenvalue", + "Hessian's 3rd eigenvalue", "Hessian's eigenvectors", + "Hessian's 1st eigenvector", + "Hessian's 2nd eigenvector", + "Hessian's 3rd eigenvector", "2nd directional derivative along gradient", "geometry tensor", *************** *** 201,205 **** "mean curvature = (K1+K2)/2", "gaussian curvature = K1*K2", ! "curvature directions", "curvature of normal streamline" }; --- 123,128 ---- "mean curvature = (K1+K2)/2", "gaussian curvature = K1*K2", ! "1st principal curvature direction", ! "2nd principal curvature direction", "curvature of normal streamline" }; *************** *** 216,220 **** --- 139,149 ---- gageSclLaplacian, gageSclHessEval, + gageSclHessEval0, + gageSclHessEval1, + gageSclHessEval2, gageSclHessEvec, + gageSclHessEvec0, + gageSclHessEvec1, + gageSclHessEvec2, gageScl2ndDD, gageSclGeomTens, *************** *** 226,253 **** gageSclMeanCurv, gageSclGaussCurv, ! gageSclCurvDir, gageSclFlowlineCurv }; ! #define GS_V gageSclValue ! #define GS_GV gageSclGradVec ! #define GS_GM gageSclGradMag ! #define GS_N gageSclNormal ! #define GS_NP gageSclNPerp ! #define GS_H gageSclHessian ! #define GS_L gageSclLaplacian ! #define GS_HA gageSclHessEval ! #define GS_HE gageSclHessEvec ! #define GS_2D gageScl2ndDD ! #define GS_GT gageSclGeomTens ! #define GS_K1 gageSclK1 ! #define GS_K2 gageSclK2 ! #define GS_CV gageSclCurvedness ! #define GS_ST gageSclShapeTrace ! #define GS_SI gageSclShapeIndex ! #define GS_MC gageSclMeanCurv ! #define GS_GC gageSclGaussCurv ! #define GS_CD gageSclCurvDir ! #define GS_FC gageSclFlowlineCurv char --- 155,190 ---- gageSclMeanCurv, gageSclGaussCurv, ! gageSclCurvDir1, ! gageSclCurvDir2, gageSclFlowlineCurv }; ! #define GS_V gageSclValue ! #define GS_GV gageSclGradVec ! #define GS_GM gageSclGradMag ! #define GS_N gageSclNormal ! #define GS_NP gageSclNPerp ! #define GS_H gageSclHessian ! #define GS_L gageSclLaplacian ! #define GS_HA gageSclHessEval ! #define GS_HA0 gageSclHessEval0 ! #define GS_HA1 gageSclHessEval1 ! #define GS_HA2 gageSclHessEval2 ! #define GS_HE gageSclHessEvec ! #define GS_HE0 gageSclHessEvec0 ! #define GS_HE1 gageSclHessEvec1 ! #define GS_HE2 gageSclHessEvec2 ! #define GS_2D gageScl2ndDD ! #define GS_GT gageSclGeomTens ! #define GS_K1 gageSclK1 ! #define GS_K2 gageSclK2 ! #define GS_CV gageSclCurvedness ! #define GS_ST gageSclShapeTrace ! #define GS_SI gageSclShapeIndex ! #define GS_MC gageSclMeanCurv ! #define GS_GC gageSclGaussCurv ! #define GS_C1 gageSclCurvDir1 ! #define GS_C2 gageSclCurvDir2 ! #define GS_FC gageSclFlowlineCurv char *************** *** 262,266 **** --- 199,209 ---- "l", "lapl", "laplacian", "heval", "h eval", "hessian eval", "hessian eigenvalues", + "heval0", + "heval1", + "heval2", "hevec", "h evec", "hessian evec", "hessian eigenvectors", + "hevec0", + "hevec1", + "hevec2", "2d", "2dd", "2nddd", "2nd", "2nd dd", "2nd dd along gradient", "gten", "geoten", "geomten", "geometry tensor", *************** *** 272,276 **** "mc", "mcurv", "meancurv", "mean curvature", "gc", "gcurv", "gausscurv", "gaussian curvature", ! "cdir", "c dir", "curvdir", "curv dir", "curvature directions", "fc", "flowlinecurv", "flowline curv", "flowline curvature", "" --- 215,220 ---- "mc", "mcurv", "meancurv", "mean curvature", "gc", "gcurv", "gausscurv", "gaussian curvature", ! "cdir1", "c dir1", "curvdir1", "curv dir1", "curvature direction 1", ! "cdir2", "c dir2", "curvdir2", "curv dir2", "curvature direction 2", "fc", "flowlinecurv", "flowline curv", "flowline curvature", "" *************** *** 287,291 **** --- 231,241 ---- GS_L, GS_L, GS_L, GS_HA, GS_HA, GS_HA, GS_HA, + GS_HA0, + GS_HA1, + GS_HA2, GS_HE, GS_HE, GS_HE, GS_HE, + GS_HE0, + GS_HE1, + GS_HE2, GS_2D, GS_2D, GS_2D, GS_2D, GS_2D, GS_2D, GS_GT, GS_GT, GS_GT, GS_GT, *************** *** 297,301 **** GS_MC, GS_MC, GS_MC, GS_MC, GS_GC, GS_GC, GS_GC, GS_GC, ! GS_CD, GS_CD, GS_CD, GS_CD, GS_CD, GS_FC, GS_FC, GS_FC, GS_FC }; --- 247,252 ---- GS_MC, GS_MC, GS_MC, GS_MC, GS_GC, GS_GC, GS_GC, GS_GC, ! GS_C1, GS_C1, GS_C1, GS_C1, GS_C1, ! GS_C2, GS_C2, GS_C2, GS_C2, GS_C2, GS_FC, GS_FC, GS_FC, GS_FC }; *************** *** 304,308 **** _gageScl = { "gageScl", ! GAGE_SCL_MAX+1, _gageSclStr, _gageSclVal, _gageSclDesc, --- 255,259 ---- _gageScl = { "gageScl", ! GAGE_SCL_ITEM_MAX+1, _gageSclStr, _gageSclVal, _gageSclDesc, *************** *** 319,328 **** 0, 1, ! GAGE_SCL_MAX, ! gageSclAnsLength, ! gageSclAnsOffset, ! GAGE_SCL_TOTAL_ANS_LENGTH, ! _gageSclNeedDeriv, ! _gageSclPrereq, _gageSclIv3Print, _gageSclFilter, --- 270,275 ---- 0, 1, ! GAGE_SCL_ITEM_MAX, ! _gageSclTable, _gageSclIv3Print, _gageSclFilter, *************** *** 331,333 **** gageKind * gageKindScl = &_gageKindScl; - --- 278,279 ---- Index: sclanswer.c =================================================================== RCS file: /cvsroot/teem/teem/src/gage/sclanswer.c,v retrieving revision 1.14 retrieving revision 1.15 diff -C2 -d -r1.14 -r1.15 *** sclanswer.c 7 Jan 2004 15:34:29 -0000 1.14 --- sclanswer.c 13 Feb 2004 22:56:15 -0000 1.15 *************** *** 24,55 **** _gageSclAnswer (gageContext *ctx, gagePerVolume *pvl) { char me[]="_gageSclAnswer"; - unsigned int query; gage_t *ans, gmag=0, *hess, *norm, *gvec, *gten, *k1, *k2, sHess[9], curv=0; double tmpMat[9], tmpVec[3], hevec[9], heval[3]; - int *offset; gage_t len, gp1[3], gp2[3], *nPerp, nProj[9], ncTen[9]; double T, N, D; ! query = pvl->query; ! ans = pvl->ans; /* convenience pointers for work below */ ! offset = gageKindScl->ansOffset; ! hess = ans + offset[gageSclHessian]; ! gvec = ans + offset[gageSclGradVec]; ! norm = ans + offset[gageSclNormal]; ! nPerp = ans + offset[gageSclNPerp]; ! gten = ans + offset[gageSclGeomTens]; ! k1 = ans + offset[gageSclK1]; ! k2 = ans + offset[gageSclK2]; ! if (1 & (query >> gageSclValue)) { /* done if doV */ if (ctx->verbose) { fprintf(stderr, "%s: val = % 15.7f\n", me, ! (double)(ans[offset[gageSclValue]])); } } ! if (1 & (query >> gageSclGradVec)) { /* done if doD1 */ if (ctx->verbose) { --- 24,51 ---- _gageSclAnswer (gageContext *ctx, gagePerVolume *pvl) { char me[]="_gageSclAnswer"; gage_t *ans, gmag=0, *hess, *norm, *gvec, *gten, *k1, *k2, sHess[9], curv=0; double tmpMat[9], tmpVec[3], hevec[9], heval[3]; gage_t len, gp1[3], gp2[3], *nPerp, nProj[9], ncTen[9]; double T, N, D; ! ans = pvl->answer; /* convenience pointers for work below */ ! hess = pvl->directAnswer[gageSclHessian]; ! gvec = pvl->directAnswer[gageSclGradVec]; ! norm = pvl->directAnswer[gageSclNormal]; ! nPerp = pvl->directAnswer[gageSclNPerp]; ! gten = pvl->directAnswer[gageSclGeomTens]; ! k1 = pvl->directAnswer[gageSclK1]; ! k2 = pvl->directAnswer[gageSclK2]; ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclValue)) { /* done if doV */ if (ctx->verbose) { fprintf(stderr, "%s: val = % 15.7f\n", me, ! (double)(pvl->directAnswer[gageSclValue][0])); } } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclGradVec)) { /* done if doD1 */ if (ctx->verbose) { *************** *** 58,69 **** } } ! if (1 & (query >> gageSclGradMag)) { /* this is the true value of gradient magnitude */ ! gmag = ans[offset[gageSclGradMag]] = sqrt(ELL_3V_DOT(gvec, gvec)); } /* NB: it would seem that gageParmGradMagMin is completely ignored ... */ ! if (1 & (query >> gageSclNormal)) { if (gmag) { ELL_3V_SCALE(norm, 1.0/gmag, gvec); --- 54,65 ---- } } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclGradMag)) { /* this is the true value of gradient magnitude */ ! gmag = pvl->directAnswer[gageSclGradMag][0] = sqrt(ELL_3V_DOT(gvec, gvec)); } /* NB: it would seem that gageParmGradMagMin is completely ignored ... */ ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclNormal)) { if (gmag) { ELL_3V_SCALE(norm, 1.0/gmag, gvec); *************** *** 76,80 **** } } ! if (1 & (query >> gageSclNPerp)) { /* nPerp = I - outer(norm, norm) */ /* NB: this sets both nPerp and nProj */ --- 72,76 ---- } } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclNPerp)) { /* nPerp = I - outer(norm, norm) */ /* NB: this sets both nPerp and nProj */ *************** *** 85,89 **** nPerp[8] += 1; } ! if (1 & (query >> gageSclHessian)) { /* done if doD2 */ if (ctx->verbose) { --- 81,85 ---- nPerp[8] += 1; } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessian)) { /* done if doD2 */ if (ctx->verbose) { *************** *** 92,103 **** } } ! if (1 & (query >> gageSclLaplacian)) { ! ans[offset[gageSclLaplacian]] = hess[0] + hess[4] + hess[8]; if (ctx->verbose) { fprintf(stderr, "%s: lapl = %g + %g + %g = %g\n", me, ! hess[0], hess[4], hess[8], ans[offset[gageSclLaplacian]]); } } ! if (1 & (query >> gageSclHessEval)) { ELL_3M_COPY(tmpMat, hess); /* HEY: look at the return value for root multiplicity? */ --- 88,100 ---- } } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclLaplacian)) { ! pvl->directAnswer[gageSclLaplacian][0] = hess[0] + hess[4] + hess[8]; if (ctx->verbose) { fprintf(stderr, "%s: lapl = %g + %g + %g = %g\n", me, ! hess[0], hess[4], hess[8], ! pvl->directAnswer[gageSclLaplacian][0]); } } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessEval)) { ELL_3M_COPY(tmpMat, hess); /* HEY: look at the return value for root multiplicity? */ *************** *** 105,118 **** mismatch between double and gage_t */ ell_3m_eigensolve_d(heval, hevec, tmpMat, AIR_TRUE); ! ELL_3V_COPY(ans+offset[gageSclHessEval], heval); } ! if (1 & (query >> gageSclHessEvec)) { ! ELL_3M_COPY(ans+offset[gageSclHessEvec], hevec); } ! if (1 & (query >> gageScl2ndDD)) { ELL_3MV_MUL(tmpVec, hess, norm); ! ans[offset[gageScl2ndDD]] = ELL_3V_DOT(norm, tmpVec); } ! if (1 & (query >> gageSclGeomTens)) { if (gmag > ctx->parm.gradMagCurvMin) { /* parm.curvNormalSide applied here to determine the sense of the --- 102,115 ---- mismatch between double and gage_t */ ell_3m_eigensolve_d(heval, hevec, tmpMat, AIR_TRUE); ! ELL_3V_COPY(pvl->directAnswer[gageSclHessEval], heval); } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessEvec)) { ! ELL_3M_COPY(pvl->directAnswer[gageSclHessEvec], hevec); } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageScl2ndDD)) { ELL_3MV_MUL(tmpVec, hess, norm); ! pvl->directAnswer[gageScl2ndDD][0] = ELL_3V_DOT(norm, tmpVec); } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclGeomTens)) { if (gmag > ctx->parm.gradMagCurvMin) { /* parm.curvNormalSide applied here to determine the sense of the *************** *** 144,157 **** } } ! if (1 && (query >> gageSclCurvedness)) { ! curv = ans[offset[gageSclCurvedness]] = ELL_3M_FROB(gten); } ! if (1 && (query >> gageSclShapeTrace)) { ! ans[offset[gageSclShapeTrace]] = (curv ! ? ELL_3M_TRACE(gten)/curv ! : 0); } ! if ( (1 && (query >> gageSclK1)) || ! (1 && (query >> gageSclK2)) ){ T = ELL_3M_TRACE(gten); N = curv; --- 141,154 ---- } } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclCurvedness)) { ! curv = pvl->directAnswer[gageSclCurvedness][0] = ELL_3M_FROB(gten); } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclShapeTrace)) { ! pvl->directAnswer[gageSclShapeTrace][0] = (curv ! ? ELL_3M_TRACE(gten)/curv ! : 0); } ! if ( (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclK1)) || ! (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclK2)) ){ T = ELL_3M_TRACE(gten); N = curv; *************** *** 170,183 **** k2[0] = 0.5*(T - D); } ! if (1 && (query >> gageSclMeanCurv)) { ! ans[offset[gageSclMeanCurv]] = (*k1 + *k2)/2; } ! if (1 && (query >> gageSclGaussCurv)) { ! ans[offset[gageSclGaussCurv]] = (*k1)*(*k2); } ! if (1 && (query >> gageSclShapeIndex)) { ! ans[offset[gageSclShapeIndex]] = -(2/M_PI)*atan2(*k1 + *k2, *k1 - *k2); } ! if (1 & (query >> gageSclCurvDir)) { /* HEY: this only works when K1, K2, 0 are all well mutually distinct, since these are the eigenvalues of the geometry tensor, and this --- 167,181 ---- k2[0] = 0.5*(T - D); } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclMeanCurv)) { ! pvl->directAnswer[gageSclMeanCurv][0] = (*k1 + *k2)/2; } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclGaussCurv)) { ! pvl->directAnswer[gageSclGaussCurv][0] = (*k1)*(*k2); } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclShapeIndex)) { ! pvl->directAnswer[gageSclShapeIndex][0] = ! -(2/M_PI)*atan2(*k1 + *k2, *k1 - *k2); } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclCurvDir1)) { /* HEY: this only works when K1, K2, 0 are all well mutually distinct, since these are the eigenvalues of the geometry tensor, and this *************** *** 186,195 **** ELL_3M_DIAG_SET(tmpMat, gten[0] - *k1, gten[4]- *k1, gten[8] - *k1); ell_3m_1d_nullspace_d(tmpVec, tmpMat); ! ELL_3V_COPY(ans+offset[gageSclCurvDir]+0, tmpVec); ELL_3M_DIAG_SET(tmpMat, gten[0] - *k2, gten[4] - *k2, gten[8] - *k2); ell_3m_1d_nullspace_d(tmpVec, tmpMat); ! ELL_3V_COPY(ans+offset[gageSclCurvDir]+3, tmpVec); } ! if (1 & (query >> gageSclFlowlineCurv)) { if (gmag >= ctx->parm.gradMagCurvMin) { /* because of the gageSclGeomTens prerequisite, sHess, nPerp, and --- 184,199 ---- ELL_3M_DIAG_SET(tmpMat, gten[0] - *k1, gten[4]- *k1, gten[8] - *k1); ell_3m_1d_nullspace_d(tmpVec, tmpMat); ! ELL_3V_COPY(pvl->directAnswer[gageSclCurvDir1], tmpVec); ! } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclCurvDir2)) { ! /* HEY: this only works when K1, K2, 0 are all well mutually distinct, ! since these are the eigenvalues of the geometry tensor, and this ! code assumes that the eigenspaces are all one-dimensional */ ! ELL_3M_COPY(tmpMat, gten); ELL_3M_DIAG_SET(tmpMat, gten[0] - *k2, gten[4] - *k2, gten[8] - *k2); ell_3m_1d_nullspace_d(tmpVec, tmpMat); ! ELL_3V_COPY(pvl->directAnswer[gageSclCurvDir2], tmpVec); } ! if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclFlowlineCurv)) { if (gmag >= ctx->parm.gradMagCurvMin) { /* because of the gageSclGeomTens prerequisite, sHess, nPerp, and *************** *** 201,205 **** ELL_3M_ZERO_SET(ncTen); } ! ans[offset[gageSclFlowlineCurv]] = sqrt(ELL_3M_FROB(ncTen)); } return; --- 205,209 ---- ELL_3M_ZERO_SET(ncTen); } ! pvl->directAnswer[gageSclFlowlineCurv][0] = sqrt(ELL_3M_FROB(ncTen)); } return; Index: sclfilter.c =================================================================== RCS file: /cvsroot/teem/teem/src/gage/sclfilter.c,v retrieving revision 1.14 retrieving revision 1.15 diff -C2 -d -r1.14 -r1.15 *** sclfilter.c 7 Jan 2004 15:34:29 -0000 1.14 --- sclfilter.c 13 Feb 2004 22:56:15 -0000 1.15 *************** *** 397,406 **** _gageSclFilter (gageContext *ctx, gagePerVolume *pvl) { char me[]="_gageSclFilter"; ! int fd, *offset; ! gage_t *fw00, *fw11, *fw22, *ans; fd = GAGE_FD(ctx); - offset = gageKindScl->ansOffset; - ans = pvl->ans; if (!ctx->parm.k3pack) { fprintf(stderr, "!%s: sorry, 6-pack filtering not implemented\n", me); --- 397,404 ---- _gageSclFilter (gageContext *ctx, gagePerVolume *pvl) { char me[]="_gageSclFilter"; ! int fd; ! gage_t *fw00, *fw11, *fw22; fd = GAGE_FD(ctx); if (!ctx->parm.k3pack) { fprintf(stderr, "!%s: sorry, 6-pack filtering not implemented\n", me); *************** *** 415,421 **** gageScl3PFilter2(pvl->iv3, pvl->iv2, pvl->iv1, fw00, fw11, fw22, ! ans + offset[gageSclValue], ! ans + offset[gageSclGradVec], ! ans + offset[gageSclHessian], pvl->needD[0], pvl->needD[1], pvl->needD[2]); break; --- 413,419 ---- gageScl3PFilter2(pvl->iv3, pvl->iv2, pvl->iv1, fw00, fw11, fw22, ! pvl->directAnswer[gageSclValue], ! pvl->directAnswer[gageSclGradVec], ! pvl->directAnswer[gageSclHessian], pvl->needD[0], pvl->needD[1], pvl->needD[2]); break; *************** *** 423,429 **** gageScl3PFilter4(pvl->iv3, pvl->iv2, pvl->iv1, fw00, fw11, fw22, ! ans + offset[gageSclValue], ! ans + offset[gageSclGradVec], ! ans + offset[gageSclHessian], pvl->needD[0], pvl->needD[1], pvl->needD[2]); break; --- 421,427 ---- gageScl3PFilter4(pvl->iv3, pvl->iv2, pvl->iv1, fw00, fw11, fw22, ! pvl->directAnswer[gageSclValue], ! pvl->directAnswer[gageSclGradVec], ! pvl->directAnswer[gageSclHessian], pvl->needD[0], pvl->needD[1], pvl->needD[2]); break; *************** *** 432,438 **** pvl->iv3, pvl->iv2, pvl->iv1, fw00, fw11, fw22, ! ans + offset[gageSclValue], ! ans + offset[gageSclGradVec], ! ans + offset[gageSclHessian], pvl->needD[0], pvl->needD[1], pvl->needD[2]); break; --- 430,436 ---- pvl->iv3, pvl->iv2, pvl->iv1, fw00, fw11, fw22, ! pvl->directAnswer[gageSclValue], ! pvl->directAnswer[gageSclGradVec], ! pvl->directAnswer[gageSclHessian], pvl->needD[0], pvl->needD[1], pvl->needD[2]); break; Index: st.c =================================================================== RCS file: /cvsroot/teem/teem/src/gage/st.c,v retrieving revision 1.5 retrieving revision 1.6 diff -C2 -d -r1.5 -r1.6 *** st.c 7 Jan 2004 15:34:29 -0000 1.5 --- st.c 13 Feb 2004 22:56:15 -0000 1.6 *************** *** 89,92 **** --- 89,93 ---- _ixi, _iyi, _izi, ixi, iyi, izi, wi, *coordCache; gageContext *ctx; + gageQuery query; gagePerVolume *pvl; airArray *mop; *************** *** 141,145 **** if (!E) E |= gageKernelSet(ctx, gageKernel00, gk0->kernel, gk0->parm); if (!E) E |= gageKernelSet(ctx, gageKernel11, gk1->kernel, gk1->parm); ! if (!E) E |= gageQuerySet(ctx, pvl, 1 << gageSclGradVec); if (!E) E |= gageUpdate(ctx); if (E) { --- 142,148 ---- if (!E) E |= gageKernelSet(ctx, gageKernel00, gk0->kernel, gk0->parm); if (!E) E |= gageKernelSet(ctx, gageKernel11, gk1->kernel, gk1->parm); ! if (!E) GAGE_QUERY_RESET(query); ! if (!E) GAGE_QUERY_ITEM_ON(query, gageSclGradVec); ! if (!E) E |= gageQuerySet(ctx, pvl, query); if (!E) E |= gageUpdate(ctx); if (E) { Index: update.c =================================================================== RCS file: /cvsroot/teem/teem/src/gage/update.c,v retrieving revision 1.11 retrieving revision 1.12 diff -C2 -d -r1.11 -r1.12 *** update.c 7 Jan 2004 15:34:29 -0000 1.11 --- update.c 13 Feb 2004 22:56:15 -0000 1.12 *************** *** 55,59 **** char me[]="_gagePvlNeedDUpdate"; gagePerVolume *pvl; ! int i, q, d, needD[3]; if (ctx->verbose) fprintf(stderr, "%s: hello\n", me); --- 55,59 ---- char me[]="_gagePvlNeedDUpdate"; gagePerVolume *pvl; ! int i, q, needD[3]; if (ctx->verbose) fprintf(stderr, "%s: hello\n", me); *************** *** 62,72 **** if (pvl->flag[gagePvlFlagQuery]) { ELL_3V_SET(needD, 0, 0, 0); ! q = pvl->kind->queryMax+1; do { q--; ! if (pvl->query & (1 << q)) { ! for (d=0; d<=2; d++) { ! needD[d] |= (pvl->kind->needDeriv[q] & (1 << d)); ! } } } while (q); --- 62,70 ---- if (pvl->flag[gagePvlFlagQuery]) { ELL_3V_SET(needD, 0, 0, 0); ! q = pvl->kind->itemMax+1; do { q--; ! if (GAGE_QUERY_ITEM_TEST(pvl->query, q)) { ! needD[pvl->kind->table[q].needDeriv] = 1; } } while (q); Index: vecGage.c =================================================================== RCS file: /cvsroot/teem/teem/src/gage/vecGage.c,v retrieving revision 1.13 retrieving revision 1.14 diff -C2 -d -r1.13 -r1.14 *** vecGage.c 7 Jan 2004 15:34:29 -0000 1.13 --- vecGage.c 13 Feb 2004 22:56:15 -0000 1.14 *************** *** 21,124 **** #include "privateGage.h" ! /* ! gageVecVector, 0: "v", component-wise-intrpolated (CWI) vec: GT[3] ! gageVecLength, 1: "l", length of CWI vector: *GT ! gageVecNormalized, 2: "n", normalized CWI vector: GT[3] ! gageVecJacobian, 3: "j", component-wise Jacobian: GT[9] ! gageVecDivergence, 4: "d", divergence (based on Jacobian): *GT ! gageVecCurl, 5: "c", curl (based on Jacobian): GT[3] ! gageVecGradient0, 6: "g1", gradient of 1st component of vector: GT[3] ! gageVecGradient1, 7: "g2", gradient of 2nd component of vector: GT[3] ! gageVecGradient2, 8: "g3", gradient of 3rd component of vector: GT[3] ! gageVecMultiGrad, 9: "mg", sum of outer products of gradients: GT[9] ! gageVecMGFrob, 10: "mgfrob", frob norm of multi-gradient: *GT ! gageVecMGEval, 11: "mgeval", eigenvalues of multi-gradient: GT[3] ! gageVecMGEvec, 12: "mgevec", ... [truncated message content] |
