From: <kin...@us...> - 2008-02-19 13:58:42
|
Revision: 3742 http://teem.svn.sourceforge.net/teem/?rev=3742&view=rev Author: kindlmann Date: 2008-02-19 05:58:47 -0800 (Tue, 19 Feb 2008) Log Message: ----------- made some of the triple conversion code more robust, added tenTripleCalc, tenTripleCalcSingle_f, tenTripleCalcSingle_d, tenTripleConvert Modified Paths: -------------- teem/trunk/src/ten/GNUmakefile teem/trunk/src/ten/miscTen.c teem/trunk/src/ten/ten.h teem/trunk/src/ten/triple.c Added Paths: ----------- teem/trunk/src/ten/tendTconv.c teem/trunk/src/ten/tendTriple.c Modified: teem/trunk/src/ten/GNUmakefile =================================================================== --- teem/trunk/src/ten/GNUmakefile 2008-02-18 07:55:08 UTC (rev 3741) +++ teem/trunk/src/ten/GNUmakefile 2008-02-19 13:58:47 UTC (rev 3742) @@ -53,6 +53,7 @@ triple.o \ tendFlotsam.o tendGrads.o tendAnplot.o tendAnvol.o tendEval.o \ tendEvec.o tendSten.o tendExpand.o tendEvq.o tendPoint.o \ + tendTriple.o tendTconv.o \ tendAnhist.o tendMake.o tendSatin.o tendShrink.o tendGlyph.o \ tendFiber.o tendEpireg.o tendBmat.o tendEstim.o tendSim.o \ tendSlice.o tendEllipse.o tendEvecrgb.o tendNorm.o tendAnscale.o \ Modified: teem/trunk/src/ten/miscTen.c =================================================================== --- teem/trunk/src/ten/miscTen.c 2008-02-18 07:55:08 UTC (rev 3741) +++ teem/trunk/src/ten/miscTen.c 2008-02-19 13:58:47 UTC (rev 3742) @@ -106,7 +106,7 @@ if (nrrdBasicInfoCopy(nout, nin, NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) { sprintf(err, "%s:", me); - biffAdd(NRRD, err); return 1; + biffAdd(TEN, err); return 1; } return 0; @@ -206,7 +206,7 @@ if (nrrdBasicInfoCopy(nout, nin, NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) { sprintf(err, "%s:", me); - biffAdd(NRRD, err); return 1; + biffAdd(TEN, err); return 1; } return 0; Modified: teem/trunk/src/ten/ten.h =================================================================== --- teem/trunk/src/ten/ten.h 2008-02-18 07:55:08 UTC (rev 3741) +++ teem/trunk/src/ten/ten.h 2008-02-19 13:58:47 UTC (rev 3742) @@ -1049,6 +1049,11 @@ int dstType, const float src[3], const int srcType); +TEN_EXPORT void tenTripleCalcSingle_d(double dst[3], + int ttype, double ten[7]); +TEN_EXPORT void tenTripleCalcSingle_f(float dst[3], + int ttype, float ten[7]); +TEN_EXPORT int tenTripleCalc(Nrrd *nout, int ttype, const Nrrd *nten); TEN_EXPORT int tenTripleConvert(Nrrd *nout, int dstType, const Nrrd *nin, int srcType); @@ -1445,6 +1450,8 @@ F(anvol) \ F(anscale) \ F(anhist) \ +F(triple) \ +F(tconv) \ F(point) \ F(slice) \ F(norm) \ Added: teem/trunk/src/ten/tendTconv.c =================================================================== --- teem/trunk/src/ten/tendTconv.c (rev 0) +++ teem/trunk/src/ten/tendTconv.c 2008-02-19 13:58:47 UTC (rev 3742) @@ -0,0 +1,74 @@ +/* + Teem: Tools to process and visualize scientific data and images + Copyright (C) 2006, 2005 Gordon Kindlmann + Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public License + (LGPL) as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + The terms of redistributing and/or modifying this software also + include exceptions to the LGPL that facilitate static linking. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + +#include "ten.h" +#include "privateTen.h" + +#define INFO "Convert between different shape triples" +char *_tend_tconvInfoL = + (INFO + ". The triples can be eignvalues, invariants (J, K, R), " + "and lots of other things."); + +int +tend_tconvMain(int argc, char **argv, char *me, hestParm *hparm) { + int pret; + hestOpt *hopt = NULL; + char *perr, *err; + airArray *mop; + + int ttype[2]; + Nrrd *nin, *nout; + char *outS; + + hestOptAdd(&hopt, "t", "inType outType", airTypeEnum, 2, 2, ttype, NULL, + "given input and desired output type of triples", + NULL, tenTripleType); + hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &nin, "-", + "input array of triples", NULL, NULL, nrrdHestNrrd); + hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-", + "output array"); + + mop = airMopNew(); + airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways); + USAGE(_tend_tconvInfoL); + PARSE(); + airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways); + + nout = nrrdNew(); + airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways); + if (tenTripleConvert(nout, ttype[1], nin, ttype[0])) { + airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways); + fprintf(stderr, "%s: trouble converting:\n%s\n", me, err); + airMopError(mop); return 1; + } + if (nrrdSave(outS, nout, NULL)) { + airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways); + fprintf(stderr, "%s: trouble writing:\n%s\n", me, err); + airMopError(mop); return 1; + } + + airMopOkay(mop); + return 0; +} +/* TEND_CMD(anvol, INFO); */ +unrrduCmd tend_tconvCmd = { "tconv", INFO, tend_tconvMain }; Added: teem/trunk/src/ten/tendTriple.c =================================================================== --- teem/trunk/src/ten/tendTriple.c (rev 0) +++ teem/trunk/src/ten/tendTriple.c 2008-02-19 13:58:47 UTC (rev 3742) @@ -0,0 +1,73 @@ +/* + Teem: Tools to process and visualize scientific data and images + Copyright (C) 2006, 2005 Gordon Kindlmann + Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public License + (LGPL) as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + The terms of redistributing and/or modifying this software also + include exceptions to the LGPL that facilitate static linking. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + +#include "ten.h" +#include "privateTen.h" + +#define INFO "Compute volume of shape triples" +char *_tend_tripleInfoL = + (INFO + ". The triple can be eignvalues, invariants (J, K, R), " + "and lots of other things."); + +int +tend_tripleMain(int argc, char **argv, char *me, hestParm *hparm) { + int pret; + hestOpt *hopt = NULL; + char *perr, *err; + airArray *mop; + + int ttype; + Nrrd *nin, *nout; + char *outS; + + hestOptAdd(&hopt, "t", "type", airTypeEnum, 1, 1, &ttype, NULL, + "desired output triple type", NULL, tenTripleType); + hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &nin, "-", + "input tensor volume", NULL, NULL, nrrdHestNrrd); + hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-", + "output triple volume"); + + mop = airMopNew(); + airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways); + USAGE(_tend_tripleInfoL); + PARSE(); + airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways); + + nout = nrrdNew(); + airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways); + if (tenTripleCalc(nout, ttype, nin)) { + airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways); + fprintf(stderr, "%s: trouble:\n%s\n", me, err); + airMopError(mop); return 1; + } + if (nrrdSave(outS, nout, NULL)) { + airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways); + fprintf(stderr, "%s: trouble writing:\n%s\n", me, err); + airMopError(mop); return 1; + } + + airMopOkay(mop); + return 0; +} +/* TEND_CMD(anvol, INFO); */ +unrrduCmd tend_tripleCmd = { "triple", INFO, tend_tripleMain }; Modified: teem/trunk/src/ten/triple.c =================================================================== --- teem/trunk/src/ten/triple.c 2008-02-18 07:55:08 UTC (rev 3741) +++ teem/trunk/src/ten/triple.c 2008-02-19 13:58:47 UTC (rev 3742) @@ -101,7 +101,7 @@ K1 = 3*MU1; stdv = sqrt(MU2); K2 = SQRT3*stdv; - K3 = SQRT2*MU3/(stdv*stdv*stdv); + K3 = stdv ? SQRT2*MU3/(stdv*stdv*stdv) : 0; } static void @@ -111,8 +111,8 @@ _mu_ev(mu, ev); R1 = sqrt(ev[0]*ev[0] + ev[1]*ev[1] + ev[2]*ev[2]); stdv = sqrt(MU2); - R2 = (3/SQRT2)*stdv/R1; - R3 = SQRT2*MU3/(stdv*stdv*stdv); + R2 = R1 ? (3/SQRT2)*stdv/R1 : 0; + R3 = stdv ? SQRT2*MU3/(stdv*stdv*stdv) : 0; } static void @@ -121,8 +121,8 @@ R1 = sqrt(3*(MU1*MU1 + MU2)); stdv = sqrt(MU2); - R2 = (3/SQRT2)*stdv/R1; - R3 = SQRT2*MU3/(stdv*stdv*stdv); + R2 = R1 ? (3/SQRT2)*stdv/R1 : 0; + R3 = stdv ? SQRT2*MU3/(stdv*stdv*stdv) : 0; } static void @@ -140,7 +140,7 @@ wp[0] = MU1; stdv = sqrt(MU2); wp[1] = SQRT2*stdv; - mode = SQRT2*MU3/(stdv*stdv*stdv); + mode = stdv ? SQRT2*MU3/(stdv*stdv*stdv) : 0; mode = AIR_CLAMP(-1, mode, 1); wp[2] = acos(AIR_CLAMP(-1, mode, 1))/3; } @@ -161,7 +161,7 @@ R2 = sqrt(J1*J1 - 3*J2)/R1; _mu_j(mu, j); stdv = sqrt(MU2); - R3 = SQRT2*MU3/(stdv*stdv*stdv); + R3 = stdv ? SQRT2*MU3/(stdv*stdv*stdv) : 0; } static void @@ -375,15 +375,154 @@ ELL_3V_COPY(_dst, dst); } +void +tenTripleCalcSingle_d(double dst[3], int ttype, double ten[7]) { + double eval[3]; + + /* in time this can become more efficient ... */ + switch (ttype) { + case tenTripleTypeEigenvalue: + tenEigensolve_d(dst, NULL, ten); + break; + case tenTripleTypeMoment: + case tenTripleTypeXYZ: + case tenTripleTypeRThetaZ: + case tenTripleTypeRThetaPhi: + case tenTripleTypeK: + case tenTripleTypeJ: + case tenTripleTypeWheelParm: + tenEigensolve_d(eval, NULL, ten); + tenTripleConvertSingle_d(dst, ttype, eval, tenTripleTypeEigenvalue); + break; + case tenTripleTypeR: + dst[0] = sqrt(_tenAnisoTen_d[tenAniso_S](ten)); + dst[1] = _tenAnisoTen_d[tenAniso_FA](ten); + dst[2] = _tenAnisoTen_d[tenAniso_Mode](ten); + break; + default: + /* what on earth? */ + ELL_3V_SET(dst, AIR_NAN, AIR_NAN, AIR_NAN); + } + return; +} + +void +tenTripleCalcSingle_f(float dst[3], int ttype, float ten[7]) { + double dst_d[3], ten_d[7]; + + TEN_T_COPY(ten_d, ten); + tenTripleCalcSingle_d(dst_d, ttype, ten_d); + ELL_3V_COPY_TT(dst, float, dst_d); + return; +} + int +tenTripleCalc(Nrrd *nout, int ttype, const Nrrd *nten) { + char me[]="tenTripleCalc", err[BIFF_STRLEN]; + size_t II, NN, size[NRRD_DIM_MAX]; + double (*ins)(void *, size_t, double), (*lup)(const void *, size_t); + + if (!( nout && nten )) { + sprintf(err, "%s: got NULL pointer", me); + biffAdd(TEN, err); return 1; + } + if (airEnumValCheck(tenTripleType, ttype)) { + sprintf(err, "%s: got invalid %s (%d)", me, + tenTripleType->name, ttype); + biffAdd(TEN, err); return 1; + } + if (tenTensorCheck(nten, nrrdTypeDefault, AIR_FALSE, AIR_TRUE)) { + sprintf(err, "%s: didn't get a valid DT array", me); + biffAdd(TEN, err); return 1; + } + if (!( nrrdTypeFloat == nten->type || + nrrdTypeDouble == nten->type )) { + sprintf(err, "%s: need input type %s or %s, not %s\n", me, + airEnumStr(nrrdType, nrrdTypeFloat), + airEnumStr(nrrdType, nrrdTypeFloat), + airEnumStr(nrrdType, nten->type)); + } + + nrrdAxisInfoGet_nva(nten, nrrdAxisInfoSize, size); + size[0] = 3; + if (nrrdMaybeAlloc_nva(nout, nten->type, nten->dim, size)) { + sprintf(err, "%s: couldn't alloc output", me); + biffMove(TEN, err, NRRD); return 1; + } + + NN = nrrdElementNumber(nten)/7; + lup = nrrdDLookup[nten->type]; + ins = nrrdDInsert[nten->type]; + for (II=0; II<NN; II++) { + double ten[7], trip[3]; + unsigned int vv; + for (vv=0; vv<7; vv++) { + ten[vv] = lup(nten->data, vv + 7*II); + } + tenTripleCalcSingle_d(trip, ttype, ten); + for (vv=0; vv<3; vv++) { + ins(nout->data, vv + 3*II, trip[vv]); + } + } + if (nrrdAxisInfoCopy(nout, nten, NULL, (NRRD_AXIS_INFO_SIZE_BIT))) { + sprintf(err, "%s: couldn't copy axis info", me); + biffMove(TEN, err, NRRD); return 1; + } + nout->axis[0].kind = nrrdKindUnknown; + if (nrrdBasicInfoCopy(nout, nten, + NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) { + sprintf(err, "%s:", me); + biffAdd(NRRD, err); return 1; + } + + return 0; +} + +int tenTripleConvert(Nrrd *nout, int dstType, const Nrrd *nin, int srcType) { char me[]="tenTripleConvert", err[BIFF_STRLEN]; + size_t II, NN; + double (*ins)(void *, size_t, double), (*lup)(const void *, size_t); - AIR_UNUSED(nout); - AIR_UNUSED(dstType); - AIR_UNUSED(nin); - AIR_UNUSED(srcType); + if (!( nout && nin )) { + sprintf(err, "%s: got NULL pointer", me); + biffAdd(TEN, err); return 1; + } + if ( airEnumValCheck(tenTripleType, dstType) || + airEnumValCheck(tenTripleType, srcType) ) { + sprintf(err, "%s: got invalid %s dst (%d) or src (%d)", me, + tenTripleType->name, dstType, srcType); + biffAdd(TEN, err); return 1; + } + if (3 != nin->axis[0].size) { + sprintf(err, "%s: need axis[0].size 3, not " _AIR_SIZE_T_CNV, me, + nin->axis[0].size); + biffAdd(TEN, err); return 1; + } + if (nrrdTypeBlock == nin->type) { + sprintf(err, "%s: input has non-scalar %s type", + me, airEnumStr(nrrdType, nrrdTypeBlock)); + biffAdd(TEN, err); return 1; + } + if (nrrdCopy(nout, nin)) { + sprintf(err, "%s: couldn't initialize output", me); + biffMove(TEN, err, NRRD); return 1; + } + lup = nrrdDLookup[nin->type]; + ins = nrrdDInsert[nout->type]; + NN = nrrdElementNumber(nin)/3; + for (II=0; II<NN; II++) { + double src[3], dst[3]; + src[0] = lup(nin->data, 0 + 3*II); + src[1] = lup(nin->data, 1 + 3*II); + src[2] = lup(nin->data, 2 + 3*II); + tenTripleConvertSingle_d(dst, dstType, src, srcType); + ins(nout->data, 0 + 3*II, dst[0]); + ins(nout->data, 1 + 3*II, dst[1]); + ins(nout->data, 2 + 3*II, dst[2]); + } + return 0; } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |