Menu
▾
▴
[Teem-commit] SF.net SVN: teem:[5252] teem/trunk/src/ten
|
From: <kin...@us...> - 2012-06-11 01:31:52
|
Revision: 5252
http://teem.svn.sourceforge.net/teem/?rev=5252&view=rev
Author: kindlmann
Date: 2012-06-11 01:31:44 +0000 (Mon, 11 Jun 2012)
Log Message:
-----------
tenModel hacking
Modified Paths:
--------------
teem/trunk/src/ten/GNUmakefile
teem/trunk/src/ten/model1Cylinder.c
teem/trunk/src/ten/model1Stick.c
teem/trunk/src/ten/model1Tensor2.c
teem/trunk/src/ten/model1Vector2D.c
teem/trunk/src/ten/modelB0.c
teem/trunk/src/ten/modelBall.c
teem/trunk/src/ten/modelBall1Cylinder.c
teem/trunk/src/ten/modelBall1Stick.c
teem/trunk/src/ten/modelBall1StickEMD.c
teem/trunk/src/ten/modelZero.c
teem/trunk/src/ten/privateTen.h
teem/trunk/src/ten/sources.cmake
teem/trunk/src/ten/ten.h
teem/trunk/src/ten/tenModel.c
teem/trunk/src/ten/tendMfit.c
Added Paths:
-----------
teem/trunk/src/ten/model1Unit2D.c
teem/trunk/src/ten/model2Unit2D.c
Modified: teem/trunk/src/ten/GNUmakefile
===================================================================
--- teem/trunk/src/ten/GNUmakefile 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/GNUmakefile 2012-06-11 01:31:44 UTC (rev 5252)
@@ -51,7 +51,8 @@
$(L).OBJS = tensor.o chan.o aniso.o glyph.o enumsTen.o grads.o miscTen.o \
mod.o estimate.o tenGage.o tenDwiGage.o qseg.o path.o qglox.o \
fiberMethods.o fiber.o epireg.o defaultsTen.o bimod.o bvec.o \
- triple.o experSpec.o tenModel.o modelBall.o model1Stick.o model1Vector2D.o \
+ triple.o experSpec.o tenModel.o modelBall.o model1Stick.o \
+ model1Vector2D.o model1Unit2D.o model2Unit2D.o \
modelBall1Stick.o modelBall1StickEMD.o modelBall1Cylinder.o \
model1Cylinder.o model1Tensor2.o modelZero.o modelB0.o \
tendAbout.o \
Modified: teem/trunk/src/ten/model1Cylinder.c
===================================================================
--- teem/trunk/src/ten/model1Cylinder.c 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/model1Cylinder.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -1,5 +1,6 @@
/*
Teem: Tools to process and visualize scientific data and images
+ Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
@@ -23,16 +24,15 @@
#include "ten.h"
#include "privateTen.h"
-#define DOF_NUM 5
#define PARM_NUM 6
static const tenModelParmDesc
parmDesc[] = {
- /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, 0},
- /* 1 */ {"length", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0},
- /* 2 */ {"radius", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0},
- /* 3 */ {"x", -1.0, 1.0, AIR_TRUE, 0},
- /* 4 */ {"y", -1.0, 1.0, AIR_TRUE, 1},
- /* 5 */ {"z", -1.0, 1.0, AIR_TRUE, 2}
+ /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 1 */ {"length", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 2 */ {"radius", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 3 */ {"x", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 0},
+ /* 4 */ {"y", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 1},
+ /* 5 */ {"z", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 2}
};
static void
@@ -58,12 +58,6 @@
return;
}
-_TEN_PARM_ALLOC
-_TEN_PARM_RAND
-_TEN_PARM_STEP
-_TEN_PARM_DIST
-_TEN_PARM_COPY
-
static char *
parmSprint(char str[AIR_STRLEN_MED], const double *parm) {
sprintf(str, "(%g) %gX%g (%g,%g,%g)", parm[0],
@@ -71,33 +65,13 @@
return str;
}
-static int
-parmConvert(double *parmDst, const double *parmSrc,
- const tenModel *modelSrc) {
- int ret;
+_TEN_PARM_ALLOC
+_TEN_PARM_RAND
+_TEN_PARM_STEP
+_TEN_PARM_DIST
+_TEN_PARM_COPY
+_TEN_PARM_CONVERT_NOOP
- ret = 0;
- parmDst[0] = parmSrc[0];
- if (modelSrc == tenModelBall) {
-
- } else if (modelSrc == tenModel1Stick) {
-
- } else if (modelSrc == tenModelBall1Stick) {
-
- } else if (modelSrc == tenModel1Cylinder) {
-
- } else if (modelSrc == tenModel1Tensor2) {
-
- } else {
- unsigned int ii;
- for (ii=0; ii<PARM_NUM; ii++) {
- parmDst[ii] = AIR_NAN;
- }
- ret = 2;
- }
- return ret;
-}
-
_TEN_SQE
_TEN_SQE_GRAD_CENTDIFF
_TEN_SQE_FIT(tenModel1Cylinder)
@@ -108,7 +82,7 @@
tenModel
_tenModel1Cylinder = {
- TEN_MODEL_STR_CYLINDER,
+ TEN_MODEL_STR_1CYLINDER,
_TEN_MODEL_FIELDS
};
const tenModel *const tenModel1Cylinder = &_tenModel1Cylinder;
Modified: teem/trunk/src/ten/model1Stick.c
===================================================================
--- teem/trunk/src/ten/model1Stick.c 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/model1Stick.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -1,5 +1,6 @@
/*
Teem: Tools to process and visualize scientific data and images
+ Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
@@ -23,15 +24,14 @@
#include "ten.h"
#include "privateTen.h"
-#define DOF_NUM 4
#define PARM_NUM 5
static const tenModelParmDesc
parmDesc[] = {
- /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, 0},
- /* 1 */ {"diffusivity", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0},
- /* 2 */ {"x", -1.0, 1.0, AIR_TRUE, 0},
- /* 3 */ {"y", -1.0, 1.0, AIR_TRUE, 1},
- /* 4 */ {"z", -1.0, 1.0, AIR_TRUE, 2}
+ /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 1 */ {"diffusivity", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 2 */ {"x", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 0},
+ /* 3 */ {"y", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 1},
+ /* 4 */ {"z", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 2}
};
static void
@@ -64,34 +64,8 @@
_TEN_PARM_STEP
_TEN_PARM_DIST
_TEN_PARM_COPY
+_TEN_PARM_CONVERT_NOOP
-static int
-parmConvert(double *parmDst, const double *parmSrc,
- const tenModel *modelSrc) {
- int ret;
-
- ret = 0;
- parmDst[0] = parmSrc[0];
- if (modelSrc == tenModelBall) {
-
- } else if (modelSrc == tenModel1Stick) {
-
- } else if (modelSrc == tenModelBall1Stick) {
-
- } else if (modelSrc == tenModel1Cylinder) {
-
- } else if (modelSrc == tenModel1Tensor2) {
-
- } else {
- unsigned int ii;
- for (ii=0; ii<PARM_NUM; ii++) {
- parmDst[ii] = AIR_NAN;
- }
- ret = 2;
- }
- return ret;
-}
-
_TEN_SQE
_TEN_SQE_GRAD_CENTDIFF
_TEN_SQE_FIT(tenModel1Stick)
Modified: teem/trunk/src/ten/model1Tensor2.c
===================================================================
--- teem/trunk/src/ten/model1Tensor2.c 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/model1Tensor2.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -1,5 +1,6 @@
/*
Teem: Tools to process and visualize scientific data and images
+ Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
@@ -25,19 +26,18 @@
/* NOTE: this model is a single 2nd-order tensor, not a two-tensor model */
-#define DOF_NUM 6
#define PARM_NUM 7
/* 1/sqrt(2) */
#define OST 0.70710678118654752440
static const tenModelParmDesc
parmDesc[] = {
- {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, 0},
- {"Dxx", -TEN_MODEL_DIFF_MAX, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0},
- {"Dxy", -TEN_MODEL_DIFF_MAX*OST, TEN_MODEL_DIFF_MAX*OST, AIR_FALSE, 0},
- {"Dxz", -TEN_MODEL_DIFF_MAX*OST, TEN_MODEL_DIFF_MAX*OST, AIR_FALSE, 0},
- {"Dyy", -TEN_MODEL_DIFF_MAX, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0},
- {"Dyz", -TEN_MODEL_DIFF_MAX*OST, TEN_MODEL_DIFF_MAX*OST, AIR_FALSE, 0},
- {"Dzz", -TEN_MODEL_DIFF_MAX, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0}
+ /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 1 */ {"Dxx", -TEN_MODEL_DIFF_MAX, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 2 */ {"Dxy", -TEN_MODEL_DIFF_MAX*OST, TEN_MODEL_DIFF_MAX*OST, AIR_FALSE, AIR_FALSE, 0},
+ /* 3 */ {"Dxz", -TEN_MODEL_DIFF_MAX*OST, TEN_MODEL_DIFF_MAX*OST, AIR_FALSE, AIR_FALSE, 0},
+ /* 4 */ {"Dyy", -TEN_MODEL_DIFF_MAX, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 5 */ {"Dyz", -TEN_MODEL_DIFF_MAX*OST, TEN_MODEL_DIFF_MAX*OST, AIR_FALSE, AIR_FALSE, 0},
+ /* 6 */ {"Dzz", -TEN_MODEL_DIFF_MAX, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0}
};
static void
@@ -129,10 +129,9 @@
_TEN_NLL_GRAD_STUB
_TEN_NLL_FIT_STUB
-
tenModel
_tenModel1Tensor2 = {
- TEN_MODEL_STR_TENSOR2,
+ TEN_MODEL_STR_1TENSOR2,
_TEN_MODEL_FIELDS
};
const tenModel *const tenModel1Tensor2 = &_tenModel1Tensor2;
Added: teem/trunk/src/ten/model1Unit2D.c
===================================================================
--- teem/trunk/src/ten/model1Unit2D.c (rev 0)
+++ teem/trunk/src/ten/model1Unit2D.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -0,0 +1,74 @@
+/*
+ Teem: Tools to process and visualize scientific data and images
+ Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
+ Copyright (C) 2008, 2007, 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 Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+#include "ten.h"
+#include "privateTen.h"
+
+#define PARM_NUM 2
+static const tenModelParmDesc
+parmDesc[] = {
+ /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 1 */ {"th", 0, 2*AIR_PI, AIR_TRUE, AIR_FALSE, 0}
+};
+
+static void
+simulate(double *dwiSim, const double *parm, const tenExperSpec *espec) {
+ unsigned int ii;
+ double b0, theta, vec[3];
+
+ b0 = parm[0]; /* not used */
+ theta = parm[1];
+ ELL_3V_SET(vec, cos(theta), sin(theta), 0.0);
+ for (ii=0; ii<espec->imgNum; ii++) {
+ dwiSim[ii] = ELL_3V_DOT(vec, espec->grad + 3*ii);
+ }
+ return;
+}
+
+static char *
+parmSprint(char str[AIR_STRLEN_MED], const double *parm) {
+ sprintf(str, "(%g) th=%g", parm[0], parm[1]);
+ return str;
+}
+
+_TEN_PARM_ALLOC
+_TEN_PARM_RAND
+_TEN_PARM_STEP
+_TEN_PARM_DIST
+_TEN_PARM_COPY
+_TEN_PARM_CONVERT_NOOP
+
+_TEN_SQE
+_TEN_SQE_GRAD_CENTDIFF
+_TEN_SQE_FIT(tenModel1Unit2D)
+
+_TEN_NLL
+_TEN_NLL_GRAD_STUB
+_TEN_NLL_FIT_STUB
+
+tenModel
+_tenModel1Unit2D = {
+ TEN_MODEL_STR_1UNIT2D,
+ _TEN_MODEL_FIELDS
+};
+const tenModel *const tenModel1Unit2D = &_tenModel1Unit2D;
Modified: teem/trunk/src/ten/model1Vector2D.c
===================================================================
--- teem/trunk/src/ten/model1Vector2D.c 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/model1Vector2D.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -1,5 +1,6 @@
/*
Teem: Tools to process and visualize scientific data and images
+ Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
@@ -23,13 +24,12 @@
#include "ten.h"
#include "privateTen.h"
-#define DOF_NUM 2
#define PARM_NUM 3
static const tenModelParmDesc
parmDesc[] = {
- /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, 0},
- /* 1 */ {"r", 0.0, 5.0, AIR_FALSE, 0},
- /* 2 */ {"th", 0, 2*AIR_PI, AIR_FALSE, 0}
+ /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 1 */ {"r", 0.0, 5.0, AIR_FALSE, AIR_FALSE, 0},
+ /* 2 */ {"th", 0, 2*AIR_PI, AIR_TRUE, AIR_FALSE, 0}
};
static void
@@ -58,34 +58,8 @@
_TEN_PARM_STEP
_TEN_PARM_DIST
_TEN_PARM_COPY
+_TEN_PARM_CONVERT_NOOP
-static int
-parmConvert(double *parmDst, const double *parmSrc,
- const tenModel *modelSrc) {
- int ret;
-
- ret = 0;
- parmDst[0] = parmSrc[0];
- if (modelSrc == tenModelBall) {
-
- } else if (modelSrc == tenModel1Stick) {
-
- } else if (modelSrc == tenModelBall1Stick) {
-
- } else if (modelSrc == tenModel1Cylinder) {
-
- } else if (modelSrc == tenModel1Tensor2) {
-
- } else {
- unsigned int ii;
- for (ii=0; ii<PARM_NUM; ii++) {
- parmDst[ii] = AIR_NAN;
- }
- ret = 2;
- }
- return ret;
-}
-
_TEN_SQE
_TEN_SQE_GRAD_CENTDIFF
_TEN_SQE_FIT(tenModel1Vector2D)
Added: teem/trunk/src/ten/model2Unit2D.c
===================================================================
--- teem/trunk/src/ten/model2Unit2D.c (rev 0)
+++ teem/trunk/src/ten/model2Unit2D.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -0,0 +1,83 @@
+/*
+ Teem: Tools to process and visualize scientific data and images
+ Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
+ Copyright (C) 2008, 2007, 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 Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+*/
+
+#include "ten.h"
+#include "privateTen.h"
+
+#define PARM_NUM 4
+static const tenModelParmDesc
+parmDesc[] = {
+ /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 1 */ {"th0", 0, 2*AIR_PI, AIR_TRUE, AIR_FALSE, 0},
+ /* 2 */ {"frac", 0, 1, AIR_FALSE, AIR_FALSE, 0},
+ /* 3 */ {"th1", 0, 2*AIR_PI, AIR_TRUE, AIR_FALSE, 0}
+};
+
+static void
+simulate(double *dwiSim, const double *parm, const tenExperSpec *espec) {
+ unsigned int ii;
+ double b0, th0, frac, th1, vec0[3], vec1[3];
+
+ b0 = parm[0]; /* not used */
+ th0 = parm[1];
+ frac = parm[2];
+ th1 = parm[3];
+ ELL_3V_SET(vec0, cos(th0), sin(th0), 0.0);
+ ELL_3V_SET(vec1, cos(th1), sin(th1), 0.0);
+ for (ii=0; ii<espec->imgNum; ii++) {
+ double dot0, dot1;
+ dot0 = ELL_3V_DOT(vec0, espec->grad + 3*ii);
+ dot1 = ELL_3V_DOT(vec1, espec->grad + 3*ii);
+ dwiSim[ii] = AIR_LERP(frac, dot0, dot1);
+ }
+ return;
+}
+
+static char *
+parmSprint(char str[AIR_STRLEN_MED], const double *parm) {
+ sprintf(str, "(%g) (1-f)*th0=%g + (f=%g)*th1=%g",
+ parm[0], parm[1], parm[2], parm[3]);
+ return str;
+}
+
+_TEN_PARM_ALLOC
+_TEN_PARM_RAND
+_TEN_PARM_STEP
+_TEN_PARM_DIST
+_TEN_PARM_COPY
+_TEN_PARM_CONVERT_NOOP
+
+_TEN_SQE
+_TEN_SQE_GRAD_CENTDIFF
+_TEN_SQE_FIT(tenModel2Unit2D)
+
+_TEN_NLL
+_TEN_NLL_GRAD_STUB
+_TEN_NLL_FIT_STUB
+
+tenModel
+_tenModel2Unit2D = {
+ TEN_MODEL_STR_2UNIT2D,
+ _TEN_MODEL_FIELDS
+};
+const tenModel *const tenModel2Unit2D = &_tenModel2Unit2D;
Modified: teem/trunk/src/ten/modelB0.c
===================================================================
--- teem/trunk/src/ten/modelB0.c 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/modelB0.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -1,5 +1,6 @@
/*
Teem: Tools to process and visualize scientific data and images
+ Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
@@ -23,11 +24,10 @@
#include "ten.h"
#include "privateTen.h"
-#define DOF_NUM 1
#define PARM_NUM 1
static const tenModelParmDesc
parmDesc[] = {
- /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, 0},
+ /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, AIR_FALSE, 0},
};
static void
@@ -53,73 +53,16 @@
_TEN_PARM_STEP
_TEN_PARM_DIST
_TEN_PARM_COPY
+_TEN_PARM_CONVERT_NOOP
-static int
-parmConvert(double *parmDst, const double *parmSrc,
- const tenModel *modelSrc) {
- int ret;
-
- ret = 0;
- parmDst[0] = parmSrc[0];
- if (modelSrc == tenModelBall) {
-
- } else if (modelSrc == tenModel1Stick) {
-
- } else if (modelSrc == tenModelBall1Stick) {
-
- } else if (modelSrc == tenModel1Cylinder) {
-
- } else if (modelSrc == tenModel1Tensor2) {
-
- } else {
- unsigned int ii;
- for (ii=0; ii<PARM_NUM; ii++) {
- parmDst[ii] = AIR_NAN;
- }
- ret = 2;
- }
- return ret;
-}
-
_TEN_SQE
_TEN_SQE_GRAD_CENTDIFF
_TEN_SQE_FIT(tenModelB0)
_TEN_NLL
+_TEN_NLL_GRAD_STUB
+_TEN_NLL_FIT_STUB
-static void
-nllGrad(double *grad, const double *parm,
- const tenExperSpec *espec,
- double *dwiBuff, const double *dwiMeas,
- int rician, double sigma) {
-
- AIR_UNUSED(grad);
- AIR_UNUSED(parm);
- AIR_UNUSED(espec);
- AIR_UNUSED(dwiBuff);
- AIR_UNUSED(dwiMeas);
- AIR_UNUSED(rician);
- AIR_UNUSED(sigma);
- return;
-}
-
-static double
-nllFit(double *parm, const tenExperSpec *espec,
- const double *dwiMeas, const double *parmInit,
- int rician, double sigma, int knownB0) {
- unsigned int pp;
-
- AIR_UNUSED(espec);
- AIR_UNUSED(dwiMeas);
- AIR_UNUSED(rician);
- AIR_UNUSED(sigma);
- AIR_UNUSED(knownB0);
- for (pp=0; pp<PARM_NUM; pp++) {
- parm[pp] = parmInit[pp];
- }
- return 0;
-}
-
tenModel
_tenModelB0 = {
TEN_MODEL_STR_B0,
Modified: teem/trunk/src/ten/modelBall.c
===================================================================
--- teem/trunk/src/ten/modelBall.c 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/modelBall.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -1,5 +1,6 @@
/*
Teem: Tools to process and visualize scientific data and images
+ Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
@@ -23,12 +24,11 @@
#include "ten.h"
#include "privateTen.h"
-#define DOF_NUM 2
#define PARM_NUM 2
static const tenModelParmDesc
parmDesc[] = {
- /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, 0},
- /* 1 */ {"diffusivity", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0}
+ /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 1 */ {"diffusivity", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0}
};
static void
@@ -55,76 +55,38 @@
_TEN_PARM_STEP
_TEN_PARM_DIST
_TEN_PARM_COPY
+_TEN_PARM_CONVERT_NOOP
-static int
-parmConvert(double *parmDst, const double *parmSrc,
- const tenModel *modelSrc) {
- int ret;
-
- ret = 0;
- parmDst[0] = parmSrc[0];
- if (modelSrc == tenModelBall) {
-
- } else if (modelSrc == tenModel1Stick) {
-
- } else if (modelSrc == tenModelBall1Stick) {
-
- } else if (modelSrc == tenModel1Cylinder) {
-
- } else if (modelSrc == tenModel1Tensor2) {
-
- } else {
- unsigned int ii;
- for (ii=0; ii<PARM_NUM; ii++) {
- parmDst[ii] = AIR_NAN;
- }
- ret = 2;
- }
- return ret;
-}
-
_TEN_SQE
_TEN_SQE_GRAD_CENTDIFF
_TEN_SQE_FIT(tenModelBall)
_TEN_NLL
+_TEN_NLL_GRAD_STUB
+_TEN_NLL_FIT_STUB
-static void
-nllGrad(double *grad, const double *parm,
- const tenExperSpec *espec,
- double *dwiBuff, const double *dwiMeas,
- int rician, double sigma) {
-
- AIR_UNUSED(grad);
- AIR_UNUSED(parm);
- AIR_UNUSED(espec);
- AIR_UNUSED(dwiBuff);
- AIR_UNUSED(dwiMeas);
- AIR_UNUSED(rician);
- AIR_UNUSED(sigma);
- return;
-}
-
-static double
-nllFit(double *parm, const tenExperSpec *espec,
- const double *dwiMeas, const double *parmInit,
- int rician, double sigma, int knownB0) {
- unsigned int pp;
-
- AIR_UNUSED(espec);
- AIR_UNUSED(dwiMeas);
- AIR_UNUSED(rician);
- AIR_UNUSED(sigma);
- AIR_UNUSED(knownB0);
- for (pp=0; pp<PARM_NUM; pp++) {
- parm[pp] = parmInit[pp];
- }
- return 0;
-}
-
tenModel
_tenModelBall = {
+ /*
TEN_MODEL_STR_BALL,
_TEN_MODEL_FIELDS
+ */
+ "ball",
+ 2,
+ parmDesc,
+ simulate,
+ parmSprint,
+ parmAlloc,
+ parmRand,
+ parmStep,
+ parmDist,
+ parmCopy,
+ parmConvert,
+ sqe,
+ sqeGrad,
+ sqeFit,
+ nll,
+ nllGrad,
+ nllFit
};
const tenModel *const tenModelBall = &_tenModelBall;
Modified: teem/trunk/src/ten/modelBall1Cylinder.c
===================================================================
--- teem/trunk/src/ten/modelBall1Cylinder.c 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/modelBall1Cylinder.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -1,5 +1,6 @@
/*
Teem: Tools to process and visualize scientific data and images
+ Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
@@ -23,18 +24,17 @@
#include "ten.h"
#include "privateTen.h"
-#define DOF_NUM 7
#define PARM_NUM 8
static const tenModelParmDesc
parmDesc[] = {
- /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, 0},
- /* 1 */ {"diffusivity", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0},
- /* 2 */ {"fraction", 0.0, 1.0, AIR_FALSE, 0},
- /* 3 */ {"length", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0},
- /* 4 */ {"radius", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0},
- /* 5 */ {"x", -1.0, 1.0, AIR_TRUE, 0},
- /* 6 */ {"y", -1.0, 1.0, AIR_TRUE, 1},
- /* 7 */ {"z", -1.0, 1.0, AIR_TRUE, 2}
+ /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 1 */ {"diffusivity", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 2 */ {"fraction", 0.0, 1.0, AIR_FALSE, AIR_FALSE, 0},
+ /* 3 */ {"length", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 4 */ {"radius", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 5 */ {"x", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 0},
+ /* 6 */ {"y", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 1},
+ /* 7 */ {"z", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 2}
};
static void
@@ -64,12 +64,6 @@
return;
}
-_TEN_PARM_ALLOC
-_TEN_PARM_RAND
-_TEN_PARM_STEP
-_TEN_PARM_DIST
-_TEN_PARM_COPY
-
static char *
parmSprint(char str[AIR_STRLEN_MED], const double *parm) {
sprintf(str, "(%g) [(1-f) %g + (f=%g) %gX%g (%g,%g,%g)]", parm[0],
@@ -78,33 +72,13 @@
return str;
}
-static int
-parmConvert(double *parmDst, const double *parmSrc,
- const tenModel *modelSrc) {
- int ret;
+_TEN_PARM_ALLOC
+_TEN_PARM_RAND
+_TEN_PARM_STEP
+_TEN_PARM_DIST
+_TEN_PARM_COPY
+_TEN_PARM_CONVERT_NOOP
- ret = 0;
- parmDst[0] = parmSrc[0];
- if (modelSrc == tenModelBall) {
-
- } else if (modelSrc == tenModel1Stick) {
-
- } else if (modelSrc == tenModelBall1Stick) {
-
- } else if (modelSrc == tenModel1Cylinder) {
-
- } else if (modelSrc == tenModel1Tensor2) {
-
- } else {
- unsigned int ii;
- for (ii=0; ii<PARM_NUM; ii++) {
- parmDst[ii] = AIR_NAN;
- }
- ret = 2;
- }
- return ret;
-}
-
_TEN_SQE
_TEN_SQE_GRAD_CENTDIFF
_TEN_SQE_FIT(tenModelBall1Cylinder)
Modified: teem/trunk/src/ten/modelBall1Stick.c
===================================================================
--- teem/trunk/src/ten/modelBall1Stick.c 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/modelBall1Stick.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -1,5 +1,6 @@
/*
Teem: Tools to process and visualize scientific data and images
+ Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
@@ -23,17 +24,16 @@
#include "ten.h"
#include "privateTen.h"
-#define DOF_NUM 6
#define PARM_NUM 7
static const tenModelParmDesc
parmDesc[] = {
- /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, 0},
- /* 1 */ {"diff_ball", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0},
- /* 2 */ {"fraction", 0.0, 1.0, AIR_FALSE, 0},
- /* 3 */ {"diff_stick", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0},
- /* 4 */ {"x", -1.0, 1.0, AIR_TRUE, 0},
- /* 5 */ {"y", -1.0, 1.0, AIR_TRUE, 1},
- /* 6 */ {"z", -1.0, 1.0, AIR_TRUE, 2}
+ /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 1 */ {"diff_ball", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 2 */ {"fraction", 0.0, 1.0, AIR_FALSE, AIR_FALSE, 0},
+ /* 3 */ {"diff_stick", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 4 */ {"x", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 0},
+ /* 5 */ {"y", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 1},
+ /* 6 */ {"z", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 2}
};
static void
@@ -71,34 +71,8 @@
_TEN_PARM_STEP
_TEN_PARM_DIST
_TEN_PARM_COPY
+_TEN_PARM_CONVERT_NOOP
-static int
-parmConvert(double *parmDst, const double *parmSrc,
- const tenModel *modelSrc) {
- int ret;
-
- ret = 0;
- parmDst[0] = parmSrc[0];
- if (modelSrc == tenModelBall) {
-
- } else if (modelSrc == tenModel1Stick) {
-
- } else if (modelSrc == tenModelBall1Stick) {
-
- } else if (modelSrc == tenModel1Cylinder) {
-
- } else if (modelSrc == tenModel1Tensor2) {
-
- } else {
- unsigned int ii;
- for (ii=0; ii<PARM_NUM; ii++) {
- parmDst[ii] = AIR_NAN;
- }
- ret = 2;
- }
- return ret;
-}
-
_TEN_SQE
_TEN_SQE_GRAD_CENTDIFF
_TEN_SQE_FIT(tenModelBall1Stick)
Modified: teem/trunk/src/ten/modelBall1StickEMD.c
===================================================================
--- teem/trunk/src/ten/modelBall1StickEMD.c 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/modelBall1StickEMD.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -1,5 +1,6 @@
/*
Teem: Tools to process and visualize scientific data and images
+ Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
@@ -23,16 +24,15 @@
#include "ten.h"
#include "privateTen.h"
-#define DOF_NUM 5
#define PARM_NUM 6
static const tenModelParmDesc
parmDesc[] = {
- /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, 0},
- /* 1 */ {"diffusivity", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, 0},
- /* 2 */ {"fraction", 0.0, 1.0, AIR_FALSE, 0},
- /* 3 */ {"x", -1.0, 1.0, AIR_TRUE, 0},
- /* 4 */ {"y", -1.0, 1.0, AIR_TRUE, 1},
- /* 5 */ {"z", -1.0, 1.0, AIR_TRUE, 2}
+ /* 0 */ {"B0", 0.0, TEN_MODEL_B0_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 1 */ {"diffusivity", 0.0, TEN_MODEL_DIFF_MAX, AIR_FALSE, AIR_FALSE, 0},
+ /* 2 */ {"fraction", 0.0, 1.0, AIR_FALSE, AIR_FALSE, 0},
+ /* 3 */ {"x", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 0},
+ /* 4 */ {"y", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 1},
+ /* 5 */ {"z", -1.0, 1.0, AIR_FALSE, AIR_TRUE, 2}
};
static void
@@ -68,34 +68,8 @@
_TEN_PARM_STEP
_TEN_PARM_DIST
_TEN_PARM_COPY
+_TEN_PARM_CONVERT_NOOP
-static int
-parmConvert(double *parmDst, const double *parmSrc,
- const tenModel *modelSrc) {
- int ret;
-
- ret = 0;
- parmDst[0] = parmSrc[0];
- if (modelSrc == tenModelBall) {
-
- } else if (modelSrc == tenModel1Stick) {
-
- } else if (modelSrc == tenModelBall1Stick) {
-
- } else if (modelSrc == tenModel1Cylinder) {
-
- } else if (modelSrc == tenModel1Tensor2) {
-
- } else {
- unsigned int ii;
- for (ii=0; ii<PARM_NUM; ii++) {
- parmDst[ii] = AIR_NAN;
- }
- ret = 2;
- }
- return ret;
-}
-
_TEN_SQE
_TEN_SQE_GRAD_CENTDIFF
_TEN_SQE_FIT(tenModelBall1StickEMD)
Modified: teem/trunk/src/ten/modelZero.c
===================================================================
--- teem/trunk/src/ten/modelZero.c 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/modelZero.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -1,5 +1,6 @@
/*
Teem: Tools to process and visualize scientific data and images
+ Copyright (C) 2012, 2011, 2010, 2009 University of Chicago
Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
@@ -23,12 +24,11 @@
#include "ten.h"
#include "privateTen.h"
-#define DOF_NUM 0
#define PARM_NUM 0
static const tenModelParmDesc
parmDesc[] = {
/* dummy to avoid compiler error */
- {"dummy", 0.0, 0.0, AIR_FALSE, 0},
+ {"dummy", 0.0, 0.0, AIR_FALSE, AIR_FALSE, 0},
};
static void
@@ -59,123 +59,74 @@
static void
parmRand(double *parm, airRandMTState *rng, int knownB0) {
-
AIR_UNUSED(parm);
AIR_UNUSED(rng);
AIR_UNUSED(knownB0);
- return;
}
static void
parmStep(double *parm1, const double scl,
const double *grad, const double *parm0) {
-
AIR_UNUSED(parm1);
AIR_UNUSED(scl);
AIR_UNUSED(grad);
AIR_UNUSED(parm0);
- return;
}
static double
parmDist(const double *parmA, const double *parmB) {
-
AIR_UNUSED(parmA);
AIR_UNUSED(parmB);
- return 0;
+ return 0.0;
}
static void
parmCopy(double *parmA, const double *parmB) {
-
AIR_UNUSED(parmA);
AIR_UNUSED(parmB);
- return;
}
static int
parmConvert(double *parmDst, const double *parmSrc,
const tenModel *modelSrc) {
-
AIR_UNUSED(parmDst);
AIR_UNUSED(parmSrc);
AIR_UNUSED(modelSrc);
- return 2;
+ return 0;
}
-static double
-sqe(const double *parm, const tenExperSpec *espec,
- double *dwiBuff, const double *dwiMeas, int knownB0) {
+_TEN_SQE
- simulate(dwiBuff, parm, espec);
- return _tenExperSpec_sqe(dwiMeas, dwiBuff, espec, knownB0);
-}
-
static void
sqeGrad(double *grad, const double *parm0,
const tenExperSpec *espec,
double *dwiBuff, const double *dwiMeas,
int knownB0) {
-
AIR_UNUSED(grad);
AIR_UNUSED(parm0);
AIR_UNUSED(espec);
AIR_UNUSED(dwiBuff);
AIR_UNUSED(dwiMeas);
AIR_UNUSED(knownB0);
- return;
}
-static double
-sqeFit(double *parm, double *convFrac, const tenExperSpec *espec,
- double *dwiBuff, const double *dwiMeas,
- const double *parmInit, int knownB0,
- unsigned int minIter, unsigned int maxIter, double convEps) {
+_TEN_SQE_FIT(tenModelZero)
- AIR_UNUSED(parm);
- AIR_UNUSED(convFrac);
- AIR_UNUSED(espec);
- AIR_UNUSED(dwiBuff);
- AIR_UNUSED(dwiMeas);
- AIR_UNUSED(parmInit);
- AIR_UNUSED(knownB0);
- AIR_UNUSED(minIter);
- AIR_UNUSED(maxIter);
- AIR_UNUSED(convEps);
- return 0;
-}
-
_TEN_NLL
+_TEN_NLL_GRAD_STUB
-static void
-nllGrad(double *grad, const double *parm,
- const tenExperSpec *espec,
- double *dwiBuff, const double *dwiMeas,
- int rician, double sigma) {
-
- AIR_UNUSED(grad);
- AIR_UNUSED(parm);
- AIR_UNUSED(espec);
- AIR_UNUSED(dwiBuff);
- AIR_UNUSED(dwiMeas);
- AIR_UNUSED(rician);
- AIR_UNUSED(sigma);
- return;
-}
-
static double
nllFit(double *parm, const tenExperSpec *espec,
const double *dwiMeas, const double *parmInit,
int rician, double sigma, int knownB0) {
-
AIR_UNUSED(parm);
- AIR_UNUSED(parmInit);
AIR_UNUSED(espec);
AIR_UNUSED(dwiMeas);
+ AIR_UNUSED(parmInit);
AIR_UNUSED(rician);
AIR_UNUSED(sigma);
AIR_UNUSED(knownB0);
- return 0;
+ return AIR_NAN;
}
tenModel
Modified: teem/trunk/src/ten/privateTen.h
===================================================================
--- teem/trunk/src/ten/privateTen.h 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/privateTen.h 2012-06-11 01:31:44 UTC (rev 5252)
@@ -161,6 +161,7 @@
TEN_EXPORT double _tenModelSqeFitSingle(const tenModel *model,
double *testParm, double *grad,
double *parm, double *convFrac,
+ unsigned int *itersTaken,
const tenExperSpec *espec,
double *dwiBuff,
const double *dwiMeas,
@@ -222,7 +223,17 @@
\
for (ii=0; ii<PARM_NUM; ii++) { \
parm1[ii] = scl*grad[ii] + parm0[ii]; \
- parm1[ii] = AIR_CLAMP(parmDesc[ii].min, parm1[ii], parmDesc[ii].max); \
+ if (parmDesc[ii].cyclic) { \
+ double delta = parmDesc[ii].max - parmDesc[ii].min; \
+ while (parm1[ii] > parmDesc[ii].max) { \
+ parm1[ii] -= delta; \
+ } \
+ while (parm1[ii] < parmDesc[ii].min) { \
+ parm1[ii] += delta; \
+ } \
+ } else { \
+ parm1[ii] = AIR_CLAMP(parmDesc[ii].min, parm1[ii], parmDesc[ii].max); \
+ } \
if (parmDesc[ii].vec3 && 2 == parmDesc[ii].vecIdx) { \
double len; \
ELL_3V_NORM(parm1 + ii - 2, parm1 + ii - 2, len); \
@@ -234,12 +245,26 @@
static double \
parmDist(const double *parmA, const double *parmB) { \
unsigned int ii; \
- double dist, dd; \
+ double dist; \
\
dist = 0; \
for (ii=0; ii<PARM_NUM; ii++) { \
- dd = (parmA[ii] - parmB[ii])/(parmDesc[ii].max - parmDesc[ii].min); \
- dist += dd*dd; \
+ double dp1, delta; \
+ delta = parmDesc[ii].max - parmDesc[ii].min; \
+ dp1 = parmA[ii] - parmB[ii]; \
+ if (parmDesc[ii].cyclic) { \
+ double dp2, dp3; \
+ dp2 = dp1 + delta; \
+ dp3 = dp1 - delta; \
+ dp1 *= dp1; \
+ dp2 *= dp2; \
+ dp3 *= dp3; \
+ dp1 = AIR_MIN(dp1, dp2); \
+ dp1 = AIR_MIN(dp1, dp3); \
+ dist += dp1/(delta*delta); \
+ } else { \
+ dist += dp1*dp1/(delta*delta); \
+ } \
} \
return sqrt(dist); \
}
@@ -255,6 +280,20 @@
return; \
}
+#define _TEN_PARM_CONVERT_NOOP \
+static int \
+parmConvert(double *parmDst, const double *parmSrc, \
+ const tenModel *modelSrc) { \
+ unsigned int ii; \
+ \
+ AIR_UNUSED(modelSrc); \
+ parmDst[0] = parmSrc[0]; \
+ for (ii=1; ii<PARM_NUM; ii++) { \
+ parmDst[ii] = AIR_NAN; \
+ } \
+ return 1; \
+}
+
#define _TEN_SQE \
static double \
sqe(const double *parm, const tenExperSpec *espec, \
@@ -305,7 +344,8 @@
#define _TEN_SQE_FIT(model) \
static double \
-sqeFit(double *parm, double *convFrac, const tenExperSpec *espec, \
+sqeFit(double *parm, double *convFrac, unsigned int *itersTaken, \
+ const tenExperSpec *espec, \
double *dwiBuff, const double *dwiMeas, \
const double *parmInit, int knownB0, \
unsigned int minIter, unsigned int maxIter, double convEps) { \
@@ -313,7 +353,7 @@
\
return _tenModelSqeFitSingle((model), \
testparm, grad, \
- parm, convFrac, \
+ parm, convFrac, itersTaken, \
espec, dwiBuff, dwiMeas, \
parmInit, knownB0, \
minIter, maxIter, convEps); \
@@ -321,7 +361,8 @@
#define _TEN_SQE_FIT_STUB \
static double \
-sqeFit(double *parm, double *convFrac, const tenExperSpec *espec, \
+sqeFit(double *parm, double *convFrac, unsigned int *itersTaken, \
+ const tenExperSpec *espec, \
double *dwiBuff, const double *dwiMeas, \
const double *parmInit, int knownB0, \
unsigned int minIter, unsigned int maxIter, double convEps) { \
@@ -388,7 +429,7 @@
}
#define _TEN_MODEL_FIELDS \
- DOF_NUM, PARM_NUM, parmDesc, \
+ PARM_NUM, parmDesc, \
simulate, \
parmSprint, parmAlloc, parmRand, \
parmStep, parmDist, parmCopy, parmConvert, \
Modified: teem/trunk/src/ten/sources.cmake
===================================================================
--- teem/trunk/src/ten/sources.cmake 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/sources.cmake 2012-06-11 01:31:44 UTC (rev 5252)
@@ -73,6 +73,8 @@
modelBall.c
model1Stick.c
model1Vector2D.c
+ model1Unit2D.c
+ model2Unit2D.c
modelBall1StickEMD.c
modelBall1Stick.c
modelBall1Cylinder.c
Modified: teem/trunk/src/ten/ten.h
===================================================================
--- teem/trunk/src/ten/ten.h 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/ten.h 2012-06-11 01:31:44 UTC (rev 5252)
@@ -1152,6 +1152,7 @@
typedef struct {
char name[AIR_STRLEN_SMALL]; /* name */
double min, max; /* bounds */
+ int cyclic; /* non-zero if effectively min == max */
int vec3; /* non-zero if this is a coefficient
of a UNIT-LENGTH 3-vector */
unsigned int vecIdx; /* *if* this param is part of vector,
@@ -1174,7 +1175,7 @@
**
** NOTE: the current convention is to *always* have the non-DW T2
** image value ("B0") be the first value in the parameter vector that
-** a tenModel is used to do describe. The B0 value will be found either
+** a tenModel is used to do describe. The B0 value will be found either
** by trivial means (i.e copied from the image data), or by a different
** method than the rest of the model parameters, but it is stored along with
** the parameters because
@@ -1185,10 +1186,11 @@
** to be able to support this.
** On the other hand, the B0 value may not be stored along with the rest of
** the parm vec in the case of saving out whole nrrd of parm vecs.
+** The "parmNum" field below therefore always includes the B0 value
*/
typedef struct tenModel_t {
char name[AIR_STRLEN_SMALL];
- unsigned int dofNum, parmNum;
+ unsigned int parmNum;
const tenModelParmDesc *parmDesc;
/* noise free simulation */
void (*simulate)(double *dwiSim, const double *parm,
@@ -1211,7 +1213,8 @@
const tenExperSpec *espec,
double *dwiBuff, const double *dwiMeas,
int knownB0);
- double (*sqeFit)(double *parm, double *convFrac, const tenExperSpec *espec,
+ double (*sqeFit)(double *parm, double *convFrac, unsigned int *itersTaken,
+ const tenExperSpec *espec,
double *dwiBuff, const double *dwiMeas,
const double *parmInit, int knownB0,
unsigned int minIter, unsigned int maxIter,
@@ -1645,7 +1648,8 @@
const Nrrd *nB0, /* maybe NULL */
const Nrrd *nparm,
int keyValueSet);
-TEN_EXPORT int tenModelSqeFit(Nrrd *nparm, Nrrd **nsqeP,
+TEN_EXPORT int tenModelSqeFit(Nrrd *nparm,
+ Nrrd **nsqeP, Nrrd **nconvP, Nrrd **niterP,
const tenModel *model,
const tenExperSpec *espec, const Nrrd *ndwi,
int knownB0, int saveB0, int typeOut,
@@ -1683,6 +1687,12 @@
/* model1Vector2D.c */
#define TEN_MODEL_STR_1VECTOR2D "1vector2d"
TEN_EXPORT const tenModel *const tenModel1Vector2D;
+/* model1Unit2D.c */
+#define TEN_MODEL_STR_1UNIT2D "1unit2d"
+TEN_EXPORT const tenModel *const tenModel1Unit2D;
+/* model2Unit2D.c */
+#define TEN_MODEL_STR_2UNIT2D "2unit2d"
+TEN_EXPORT const tenModel *const tenModel2Unit2D;
/* model1Stick.c */
#define TEN_MODEL_STR_1STICK "1stick"
TEN_EXPORT const tenModel *const tenModel1Stick;
@@ -1696,10 +1706,10 @@
#define TEN_MODEL_STR_BALL1CYLINDER "ball1cylinder"
TEN_EXPORT const tenModel *const tenModelBall1Cylinder;
/* model1Cylinder.c */
-#define TEN_MODEL_STR_CYLINDER "1cylinder"
+#define TEN_MODEL_STR_1CYLINDER "1cylinder"
TEN_EXPORT const tenModel *const tenModel1Cylinder;
/* model1Tensor2.c: 2nd-order tensor (one of them), not two-tensor */
-#define TEN_MODEL_STR_TENSOR2 "1tensor2"
+#define TEN_MODEL_STR_1TENSOR2 "1tensor2"
TEN_EXPORT const tenModel *const tenModel1Tensor2;
/* mod.c */
Modified: teem/trunk/src/ten/tenModel.c
===================================================================
--- teem/trunk/src/ten/tenModel.c 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/tenModel.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -41,15 +41,19 @@
ret = tenModel1Stick;
} else if (!strcmp(str, TEN_MODEL_STR_1VECTOR2D)) {
ret = tenModel1Vector2D;
+ } else if (!strcmp(str, TEN_MODEL_STR_1UNIT2D)) {
+ ret = tenModel1Unit2D;
+ } else if (!strcmp(str, TEN_MODEL_STR_2UNIT2D)) {
+ ret = tenModel2Unit2D;
} else if (!strcmp(str, TEN_MODEL_STR_BALL1STICKEMD)) {
ret = tenModelBall1StickEMD;
} else if (!strcmp(str, TEN_MODEL_STR_BALL1STICK)) {
ret = tenModelBall1Stick;
} else if (!strcmp(str, TEN_MODEL_STR_BALL1CYLINDER)) {
ret = tenModelBall1Cylinder;
- } else if (!strcmp(str, TEN_MODEL_STR_CYLINDER)) {
+ } else if (!strcmp(str, TEN_MODEL_STR_1CYLINDER)) {
ret = tenModel1Cylinder;
- } else if (!strcmp(str, TEN_MODEL_STR_TENSOR2)) {
+ } else if (!strcmp(str, TEN_MODEL_STR_1TENSOR2)) {
ret = tenModel1Tensor2;
} else {
/* we don't currently have a tenModelUnknown */
@@ -335,12 +339,13 @@
** _tenModelSqeFitSingle
**
** callable function (as opposed to tenModel method) for doing
-** sqe fitting. Requires PARM_NUM length buffers testParm and grad
+** sqe fitting. Returns the sqe at the converged fit location
+** Requires PARM_NUM length buffers testParm and grad
*/
double
_tenModelSqeFitSingle(const tenModel *model,
double *testParm, double *grad,
- double *parm, double *convFrac,
+ double *parm, double *convFrac, unsigned int *itersTaken,
const tenExperSpec *espec,
double *dwiBuff, const double *dwiMeas,
const double *parmInit, int knownB0,
@@ -356,8 +361,8 @@
val = model->sqe(parm, espec, dwiBuff, dwiMeas, knownB0);
model->sqeGrad(grad, parm, espec, dwiBuff, dwiMeas, knownB0);
- opp = 2;
- bak = 0.1;
+ opp = 1.2; /* opportunistic step size increase */
+ bak = 0.2; /* scaling back because of bad step */
iter = 0;
dist = convEps*8;
do {
@@ -380,31 +385,33 @@
: (iter > maxIter) || dist < convEps);
} while (!done);
*convFrac = dist/convEps;
+ *itersTaken = iter;
return val;
}
int
-tenModelSqeFit(Nrrd *nparm, Nrrd **nsqeP,
+tenModelSqeFit(Nrrd *nparm,
+ Nrrd **nsqeP, Nrrd **nconvP, Nrrd **niterP,
const tenModel *model,
const tenExperSpec *espec, const Nrrd *ndwi,
int knownB0, int saveB0, int typeOut,
unsigned int minIter, unsigned int maxIter,
unsigned int starts, double convEps, airRandMTState *_rng) {
static const char me[]="tenModelSqeFit";
- double *ddwi, *dwibuff, sqe, sqeBest, convFrac,
+ double *ddwi, *dwibuff, sqe, sqeBest,
*dparm, *dparmBest,
(*ins)(void *v, size_t I, double d),
(*lup)(const void *v, size_t I);
airArray *mop;
- unsigned int saveParmNum, dwiNum, ii, lablen;
+ unsigned int saveParmNum, dwiNum, ii, lablen, itersTaken;
size_t szOut[NRRD_DIM_MAX], II, numSamp;
int axmap[NRRD_DIM_MAX], erraxmap[NRRD_DIM_MAX];
const char *dwi;
char *parm;
airRandMTState *rng;
- Nrrd *nsqe;
-
- /* nsqeP can be NULL */
+ Nrrd *nsqe, *nconv, *niter;
+
+ /* nsqeP, nconvP, niterP can be NULL */
if (!( nparm && model && espec && ndwi )) {
biffAddf(TEN, "%s: got NULL pointer", me);
return 1;
@@ -458,15 +465,41 @@
nsqe = *nsqeP;
if (!nsqe) {
nsqe = nrrdNew();
- if (nrrdMaybeAlloc_nva(nsqe, typeOut, ndwi->dim-1, szOut+1)) {
- biffMovef(TEN, NRRD, "%s: couldn't allocate error output", me);
- airMopError(mop); return 1;
- }
*nsqeP = nsqe;
}
+ if (nrrdMaybeAlloc_nva(nsqe, typeOut, ndwi->dim-1, szOut+1)) {
+ biffMovef(TEN, NRRD, "%s: couldn't allocate error output", me);
+ airMopError(mop); return 1;
+ }
} else {
nsqe = NULL;
}
+ if (nconvP) {
+ nconv = *nconvP;
+ if (!nconv) {
+ nconv = nrrdNew();
+ *nconvP = nconv;
+ }
+ if (nrrdMaybeAlloc_nva(nconv, nrrdTypeDouble, ndwi->dim-1, szOut+1)) {
+ biffMovef(TEN, NRRD, "%s: couldn't allocate conv output", me);
+ airMopError(mop); return 1;
+ }
+ } else {
+ nconv = NULL;
+ }
+ if (niterP) {
+ niter = *niterP;
+ if (!niter) {
+ niter = nrrdNew();
+ *niterP = niter;
+ }
+ if (nrrdMaybeAlloc_nva(niter, nrrdTypeUInt, ndwi->dim-1, szOut+1)) {
+ biffMovef(TEN, NRRD, "%s: couldn't allocate iter output", me);
+ airMopError(mop); return 1;
+ }
+ } else {
+ niter = NULL;
+ }
ddwi = AIR_CALLOC(espec->imgNum, double);
dwibuff = AIR_CALLOC(espec->imgNum, double);
if (!(ddwi && dwibuff)) {
@@ -488,8 +521,10 @@
ins = nrrdDInsert[typeOut];
parm = AIR_CAST(char *, nparm->data);
dwi = AIR_CAST(char *, ndwi->data);
+ itersTaken = 0;
for (II=0; II<numSamp; II++) {
- unsigned int ss;
+ double cvf, convFrac;
+ unsigned int ss, itak;
for (ii=0; ii<dwiNum; ii++) {
ddwi[ii] = lup(dwi, ii);
}
@@ -499,20 +534,29 @@
dparm[0] = tenExperSpecKnownB0Get(espec, ddwi);
}
model->rand(dparm, rng, knownB0);
- sqe = model->sqeFit(dparm, &convFrac, espec, dwibuff, ddwi,
+ sqe = model->sqeFit(dparm, &cvf, &itak,
+ espec, dwibuff, ddwi,
dparm, knownB0, minIter, maxIter, convEps);
if (sqe < sqeBest) {
sqeBest = sqe;
model->copy(dparmBest, dparm);
+ itersTaken = itak;
+ convFrac = cvf;
}
}
for (ii=0; ii<saveParmNum; ii++) {
ins(parm, ii, saveB0 ? dparmBest[ii] : dparmBest[ii+1]);
}
+ /* save things about fitting into nrrds */
if (nsqeP) {
ins(nsqe->data, II, sqeBest);
}
- /* possibly save sqeBest */
+ if (nconvP) {
+ nrrdDInsert[nrrdTypeDouble](nconv->data, II, convFrac);
+ }
+ if (niterP) {
+ nrrdDInsert[nrrdTypeUInt](niter->data, II, itersTaken);
+ }
parm += saveParmNum*nrrdTypeSize[typeOut];
dwi += espec->imgNum*nrrdTypeSize[ndwi->type];
}
@@ -548,6 +592,40 @@
airMopError(mop); return 1;
}
}
+ if (nconvP) {
+ if (nrrdAxisInfoCopy(nconv, ndwi, erraxmap, NRRD_AXIS_INFO_SIZE_BIT)
+ || nrrdBasicInfoCopy(nconv, ndwi,
+ NRRD_BASIC_INFO_DATA_BIT
+ | NRRD_BASIC_INFO_TYPE_BIT
+ | NRRD_BASIC_INFO_BLOCKSIZE_BIT
+ | NRRD_BASIC_INFO_DIMENSION_BIT
+ | NRRD_BASIC_INFO_CONTENT_BIT
+ | NRRD_BASIC_INFO_COMMENTS_BIT
+ | (nrrdStateKeyValuePairsPropagate
+ ? 0
+ : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {
+ biffMovef(TEN, NRRD,
+ "%s: couldn't copy axis or basic info to conv out", me);
+ airMopError(mop); return 1;
+ }
+ }
+ if (niterP) {
+ if (nrrdAxisInfoCopy(niter, ndwi, erraxmap, NRRD_AXIS_INFO_SIZE_BIT)
+ || nrrdBasicInfoCopy(niter, ndwi,
+ NRRD_BASIC_INFO_DATA_BIT
+ | NRRD_BASIC_INFO_TYPE_BIT
+ | NRRD_BASIC_INFO_BLOCKSIZE_BIT
+ | NRRD_BASIC_INFO_DIMENSION_BIT
+ | NRRD_BASIC_INFO_CONTENT_BIT
+ | NRRD_BASIC_INFO_COMMENTS_BIT
+ | (nrrdStateKeyValuePairsPropagate
+ ? 0
+ : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {
+ biffMovef(TEN, NRRD,
+ "%s: couldn't copy axis or basic info to iter out", me);
+ airMopError(mop); return 1;
+ }
+ }
lablen = (strlen(tenModelPrefixStr)
+ (saveB0 ? strlen("B0+") : 0)
+ strlen(model->name)
Modified: teem/trunk/src/ten/tendMfit.c
===================================================================
--- teem/trunk/src/ten/tendMfit.c 2012-06-08 20:18:54 UTC (rev 5251)
+++ teem/trunk/src/ten/tendMfit.c 2012-06-11 01:31:44 UTC (rev 5252)
@@ -35,8 +35,8 @@
char *perr, *err;
airArray *mop;
- Nrrd *nin, *nout, *nterr;
- char *outS, *terrS, *modS;
+ Nrrd *nin, *nout, *nterr, *nconv, *niter;
+ char *outS, *terrS, *convS, *iterS, *modS;
int knownB0, saveB0, verbose, mlfit, typeOut;
unsigned int maxIter, minIter, starts;
double sigma, eps;
@@ -86,6 +86,13 @@
hestOptAdd(&hopt, "eo", "filename", airTypeString, 1, 1, &terrS, "",
"Giving a filename here allows you to save out the per-sample "
"fitting error. By default, no such error is saved.");
+ hestOptAdd(&hopt, "co", "filename", airTypeString, 1, 1, &convS, "",
+ "Giving a filename here allows you to save out the per-sample "
+ "convergence fraction. By default, no such error is saved.");
+ hestOptAdd(&hopt, "io", "filename", airTypeString, 1, 1, &iterS, "",
+ "Giving a filename here allows you to save out the per-sample "
+ "number of iterations needed for fitting. "
+ "By default, no such error is saved.");
mop = airMopNew();
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
@@ -94,6 +101,8 @@
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
nterr = NULL;
+ nconv = NULL;
+ niter = NULL;
espec = tenExperSpecNew();
airMopAdd(mop, espec, (airMopper)tenExperSpecNix, airMopAlways);
nout = nrrdNew();
@@ -110,6 +119,8 @@
}
if (tenModelSqeFit(nout,
airStrlen(terrS) ? &nterr : NULL,
+ airStrlen(convS) ? &nconv : NULL,
+ airStrlen(iterS) ? &niter : NULL,
model, espec, nin,
knownB0, saveB0, typeOut,
minIter, maxIter, starts, eps, NULL)) {
@@ -119,7 +130,9 @@
}
if (nrrdSave(outS, nout, NULL)
- || (airStrlen(terrS) && nrrdSave(terrS, nterr, NULL))) {
+ || (airStrlen(terrS) && nrrdSave(terrS, nterr, NULL))
+ || (airStrlen(convS) && nrrdSave(convS, nconv, NULL))
+ || (airStrlen(iterS) && nrrdSave(iterS, niter, NULL))) {
airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
fprintf(stderr, "%s: trouble writing output:\n%s\n", me, err);
airMopError(mop); return 1;
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|