Menu
▾
▴
[Teem-commit] SF.net SVN: teem:[7455] teem/trunk/src/ten
|
From: <kin...@us...> - 2025-09-18 15:34:20
|
Revision: 7455
http://sourceforge.net/p/teem/code/7455
Author: kindlmann
Date: 2025-09-18 15:34:17 +0000 (Thu, 18 Sep 2025)
Log Message:
-----------
ABI CHANGE: as part of shifting some ten/test/*.c programs to the new typed hestOptAdd_ functions, discovered that some fields in the tenEMBimodalParm struct were doubles when they should have been (boolean) ints: twoStage and verbose. They were set by the non-type-checked hestOptAdd() as ints, which still worked, but the type error has now been fixed
Modified Paths:
--------------
teem/trunk/src/ten/ten.h
teem/trunk/src/ten/test/cntr.c
teem/trunk/src/ten/test/igrt.c
teem/trunk/src/ten/test/odf-hist.c
teem/trunk/src/ten/test/teigen.c
teem/trunk/src/ten/test/tem.c
teem/trunk/src/ten/test/tg.c
teem/trunk/src/ten/test/tqgl.c
teem/trunk/src/ten/test/tt.c
Modified: teem/trunk/src/ten/ten.h
===================================================================
--- teem/trunk/src/ten/ten.h 2025-09-18 08:28:52 UTC (rev 7454)
+++ teem/trunk/src/ten/ten.h 2025-09-18 15:34:17 UTC (rev 7455)
@@ -911,11 +911,11 @@
minDelta, /* convergence test for maximization */
minFraction, /* smallest fraction (in 0.0 to 1.0) that material
1 or 2 can legitimately have */
- minConfidence, /* smallest confidence value that the model fitting
- is allowed to have */
- twoStage, /* wacky two-stage fitting */
+ minConfidence; /* smallest confidence value that the model fitting
+ is allowed to have */
+ unsigned int maxIteration; /* cap on # of non-convergent iters allowed */
+ int twoStage, /* wacky two-stage fitting */
verbose; /* output messages and/or progress images */
- unsigned int maxIteration; /* cap on # of non-convergent iters allowed */
/* ----- internal ----- */
double *histo, /* double version of histogram */
*pp1, *pp2, /* pre-computed posterior probabilities for the
Modified: teem/trunk/src/ten/test/cntr.c
===================================================================
--- teem/trunk/src/ten/test/cntr.c 2025-09-18 08:28:52 UTC (rev 7454)
+++ teem/trunk/src/ten/test/cntr.c 2025-09-18 15:34:17 UTC (rev 7455)
@@ -17,7 +17,6 @@
along with this library; if not, see <https://www.gnu.org/licenses/>.
*/
-
#include "../ten.h"
const char *info = ("does contraction between 2 2nd-order "
@@ -27,24 +26,25 @@
main(int argc, const char *argv[]) {
const char *me;
char *err;
- hestOpt *hopt=NULL;
+ hestOpt *hopt = NULL;
airArray *mop;
- char *outS;
- Nrrd *_ncov, *ncov, *_nten[2], *nten[2], *nout;
+ Nrrd *ncov, *nten[2], *nout;
double *cc, *t0, *t1, *out, ww[21];
size_t nn, ii;
mop = airMopNew();
me = argv[0];
- hestOptAdd(&hopt, "i4", "volume", airTypeOther, 1, 1, &_ncov, NULL,
- "4th-order tensor volume", NULL, NULL, nrrdHestNrrd);
- hestOptAdd(&hopt, "i2", "v0 v1", airTypeOther, 2, 2, _nten, NULL,
- "two 2nd-order tensor volumes", NULL, NULL, nrrdHestNrrd);
- hestOptAdd(&hopt, "o", "filename", airTypeString, 1, 1, &outS, "-",
- "file to write output nrrd to");
- hestParseOrDie(hopt, argc-1, argv+1, NULL,
- me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
+ Nrrd *_ncov;
+ hestOptAdd_1_Other(&hopt, "i4", "volume", &_ncov, NULL, "4th-order tensor volume",
+ nrrdHestNrrd);
+ Nrrd *_nten[2];
+ hestOptAdd_2_Other(&hopt, "i2", "v0 v1", _nten, NULL, "two 2nd-order tensor volumes",
+ nrrdHestNrrd);
+ char *outS;
+ hestOptAdd_1_String(&hopt, "o", "filename", &outS, "-",
+ "file to write output nrrd to");
+ hestParseOrDie(hopt, argc - 1, argv + 1, NULL, me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
@@ -56,17 +56,17 @@
return 1;
}
if (!(4 == _ncov->dim && 21 == _ncov->axis[0].size)) {
- fprintf(stderr, "%s: didn't get a 4-D 21-by-X volume (got %u-D %u-by-X)\n",
- me, _ncov->dim, AIR_UINT(_ncov->axis[0].size));
+ fprintf(stderr, "%s: didn't get a 4-D 21-by-X volume (got %u-D %u-by-X)\n", me,
+ _ncov->dim, AIR_UINT(_ncov->axis[0].size));
airMopError(mop);
return 1;
}
- if (!(nrrdElementNumber(_ncov)/21 == nrrdElementNumber(_nten[0])/7
- && nrrdElementNumber(_nten[0])/7 == nrrdElementNumber(_nten[1])/7)) {
+ if (!(nrrdElementNumber(_ncov) / 21 == nrrdElementNumber(_nten[0]) / 7
+ && nrrdElementNumber(_nten[0]) / 7 == nrrdElementNumber(_nten[1]) / 7)) {
fprintf(stderr, "%s: number voxels %u %u %u don't all match\n", me,
- AIR_UINT(nrrdElementNumber(_ncov)/21),
- AIR_UINT(nrrdElementNumber(_nten[0])/7),
- AIR_UINT(nrrdElementNumber(_nten[1])/7));
+ AIR_UINT(nrrdElementNumber(_ncov) / 21),
+ AIR_UINT(nrrdElementNumber(_nten[0]) / 7),
+ AIR_UINT(nrrdElementNumber(_nten[1]) / 7));
airMopError(mop);
return 1;
}
@@ -83,8 +83,8 @@
|| nrrdConvert(nten[0], _nten[0], nrrdTypeDouble)
|| nrrdConvert(nten[1], _nten[1], nrrdTypeDouble)) {
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
- fprintf(stderr, "%s: trouble converting to %s:\n%s\n",
- me, airEnumStr(nrrdType, nrrdTypeDouble), err);
+ fprintf(stderr, "%s: trouble converting to %s:\n%s\n", me,
+ airEnumStr(nrrdType, nrrdTypeDouble), err);
airMopError(mop);
return 1;
}
@@ -101,12 +101,27 @@
out = AIR_CAST(double *, nout->data);
nn = nrrdElementNumber(nout);
- ww[ 0] = 1*1; ww[ 1] = 2*1; ww[ 2] = 2*1; ww[ 3] = 1*1; ww[ 4] = 2*1; ww[ 5] = 1*1;
- /* */ ww[ 6] = 2*2; ww[ 7] = 2*2; ww[ 8] = 1*2; ww[ 9] = 2*2; ww[10] = 1*2;
- /* */ ww[11] = 2*2; ww[12] = 1*2; ww[13] = 2*2; ww[14] = 1*2;
- /* */ ww[15] = 1*1; ww[16] = 2*1; ww[17] = 1*1;
- /* */ ww[18] = 2*2; ww[19] = 1*2;
- /* */ ww[20] = 1*1;
+ ww[0] = 1 * 1;
+ ww[1] = 2 * 1;
+ ww[2] = 2 * 1;
+ ww[3] = 1 * 1;
+ ww[4] = 2 * 1;
+ ww[5] = 1 * 1;
+ /* */ ww[6] = 2 * 2;
+ ww[7] = 2 * 2;
+ ww[8] = 1 * 2;
+ ww[9] = 2 * 2;
+ ww[10] = 1 * 2;
+ /* */ ww[11] = 2 * 2;
+ ww[12] = 1 * 2;
+ ww[13] = 2 * 2;
+ ww[14] = 1 * 2;
+ /* */ ww[15] = 1 * 1;
+ ww[16] = 2 * 1;
+ ww[17] = 1 * 1;
+ /* */ ww[18] = 2 * 2;
+ ww[19] = 1 * 2;
+ /* */ ww[20] = 1 * 1;
/*
for (ii=0; ii<21; ii++) {
@@ -114,14 +129,26 @@
}
*/
- for (ii=0; ii<nn; ii++) {
+ for (ii = 0; ii < nn; ii++) {
- out[ii] = (+ cc[ 0]*ww[ 0]*t0[1]*t1[1] + cc[ 1]*ww[ 1]*t0[2]*t1[1] + cc[ 2]*ww[ 2]*t0[3]*t1[1] + cc[ 3]*ww[ 3]*t0[4]*t1[1] + cc[ 4]*ww[ 4]*t0[5]*t1[1] + cc[ 5]*ww[ 5]*t0[6]*t1[1] +
- + cc[ 1]*ww[ 1]*t0[1]*t1[2] + cc[ 6]*ww[ 6]*t0[2]*t1[2] + cc[ 7]*ww[ 7]*t0[3]*t1[2] + cc[ 8]*ww[ 8]*t0[4]*t1[2] + cc[ 9]*ww[ 9]*t0[5]*t1[2] + cc[10]*ww[10]*t0[6]*t1[2] +
- + cc[ 2]*ww[ 2]*t0[1]*t1[3] + cc[ 7]*ww[ 7]*t0[2]*t1[3] + cc[11]*ww[11]*t0[3]*t1[3] + cc[12]*ww[12]*t0[4]*t1[3] + cc[13]*ww[13]*t0[5]*t1[3] + cc[14]*ww[14]*t0[6]*t1[3] +
- + cc[ 3]*ww[ 3]*t0[1]*t1[4] + cc[ 8]*ww[ 8]*t0[2]*t1[4] + cc[12]*ww[12]*t0[3]*t1[4] + cc[15]*ww[15]*t0[4]*t1[4] + cc[16]*ww[16]*t0[5]*t1[4] + cc[17]*ww[17]*t0[6]*t1[4] +
- + cc[ 4]*ww[ 4]*t0[1]*t1[5] + cc[ 9]*ww[ 9]*t0[2]*t1[5] + cc[13]*ww[13]*t0[3]*t1[5] + cc[16]*ww[16]*t0[4]*t1[5] + cc[18]*ww[18]*t0[5]*t1[5] + cc[19]*ww[19]*t0[6]*t1[5] +
- + cc[ 5]*ww[ 5]*t0[1]*t1[6] + cc[10]*ww[10]*t0[2]*t1[6] + cc[14]*ww[14]*t0[3]*t1[6] + cc[17]*ww[17]*t0[4]*t1[6] + cc[19]*ww[19]*t0[5]*t1[6] + cc[20]*ww[20]*t0[6]*t1[6]);
+ out[ii] = (+cc[0] * ww[0] * t0[1] * t1[1] + cc[1] * ww[1] * t0[2] * t1[1]
+ + cc[2] * ww[2] * t0[3] * t1[1] + cc[3] * ww[3] * t0[4] * t1[1]
+ + cc[4] * ww[4] * t0[5] * t1[1] + cc[5] * ww[5] * t0[6] * t1[1]
+ + +cc[1] * ww[1] * t0[1] * t1[2] + cc[6] * ww[6] * t0[2] * t1[2]
+ + cc[7] * ww[7] * t0[3] * t1[2] + cc[8] * ww[8] * t0[4] * t1[2]
+ + cc[9] * ww[9] * t0[5] * t1[2] + cc[10] * ww[10] * t0[6] * t1[2]
+ + +cc[2] * ww[2] * t0[1] * t1[3] + cc[7] * ww[7] * t0[2] * t1[3]
+ + cc[11] * ww[11] * t0[3] * t1[3] + cc[12] * ww[12] * t0[4] * t1[3]
+ + cc[13] * ww[13] * t0[5] * t1[3] + cc[14] * ww[14] * t0[6] * t1[3]
+ + +cc[3] * ww[3] * t0[1] * t1[4] + cc[8] * ww[8] * t0[2] * t1[4]
+ + cc[12] * ww[12] * t0[3] * t1[4] + cc[15] * ww[15] * t0[4] * t1[4]
+ + cc[16] * ww[16] * t0[5] * t1[4] + cc[17] * ww[17] * t0[6] * t1[4]
+ + +cc[4] * ww[4] * t0[1] * t1[5] + cc[9] * ww[9] * t0[2] * t1[5]
+ + cc[13] * ww[13] * t0[3] * t1[5] + cc[16] * ww[16] * t0[4] * t1[5]
+ + cc[18] * ww[18] * t0[5] * t1[5] + cc[19] * ww[19] * t0[6] * t1[5]
+ + +cc[5] * ww[5] * t0[1] * t1[6] + cc[10] * ww[10] * t0[2] * t1[6]
+ + cc[14] * ww[14] * t0[3] * t1[6] + cc[17] * ww[17] * t0[4] * t1[6]
+ + cc[19] * ww[19] * t0[5] * t1[6] + cc[20] * ww[20] * t0[6] * t1[6]);
/* 0:xxxx 1:xxxy 2:xxxz 3:xxyy 4:xxyz 5:xxzz
* 6:xyxy 7:xyxz 8:xyyy 9:xyyz 10:xyzz
Modified: teem/trunk/src/ten/test/igrt.c
===================================================================
--- teem/trunk/src/ten/test/igrt.c 2025-09-18 08:28:52 UTC (rev 7454)
+++ teem/trunk/src/ten/test/igrt.c 2025-09-18 15:34:17 UTC (rev 7455)
@@ -17,7 +17,6 @@
along with this library; if not, see <https://www.gnu.org/licenses/>.
*/
-
#include "../ten.h"
const char *info = ("tests invariant grads and rotation tangents.");
@@ -24,36 +23,32 @@
int
main(int argc, const char *argv[]) {
- const char *me;
- hestOpt *hopt=NULL;
- airArray *mop;
+ const char *me = argv[0];
+ hestOpt *hopt = NULL;
+ airArray *mop = airMopNew();
- double _ten[6], ten[7], minnorm, igrt[6][7], eval[3], evec[9],
- pp[3], qq[4], rot[9], matA[9], matB[9], tmp;
- int doK, ret, ii, jj;
- mop = airMopNew();
-
- me = argv[0];
- hestOptAdd(&hopt, NULL, "tensor", airTypeDouble, 6, 6, _ten, NULL,
- "tensor value");
- hestOptAdd(&hopt, "mn", "minnorm", airTypeDouble, 1, 1, &minnorm,
- "0.00001",
+ double _ten[6];
+ hestOptAdd_N_Double(&hopt, NULL, "tensor", 6, _ten, NULL, "tensor value");
+ double minnorm;
+ hestOptAdd(&hopt, "mn", "minnorm", airTypeDouble, 1, 1, &minnorm, "0.00001",
"minimum norm before special handling");
- hestOptAdd(&hopt, "k", NULL, airTypeInt, 0, 0, &doK, NULL,
- "Use K invariants, instead of R (the default)");
- hestOptAdd(&hopt, "p", "x y z", airTypeDouble, 3, 3, pp, "0 0 0",
- "location in quaternion quotient space");
- hestParseOrDie(hopt, argc-1, argv+1, NULL,
- me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
+ int doK;
+ hestOptAdd_Flag(&hopt, "k", &doK, "Use K invariants, instead of R (the default)");
+ double pp[3];
+ hestOptAdd_3_Double(&hopt, "p", "x y z", pp, "0 0 0",
+ "location in quaternion quotient space");
+ hestParseOrDie(hopt, argc - 1, argv + 1, NULL, me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
- ELL_6V_COPY(ten+1, _ten);
+ double ten[7];
+ ELL_6V_COPY(ten + 1, _ten);
ten[0] = 1.0;
- fprintf(stderr, "input tensor = %f %f %f %f %f %f\n",
- ten[1], ten[2], ten[3], ten[4], ten[5], ten[6]);
+ fprintf(stderr, "input tensor = %f %f %f %f %f %f\n", ten[1], ten[2], ten[3],
+ ten[4], ten[5], ten[6]);
+ double qq[4], rot[9], matA[9], matB[9], tmp;
ELL_4V_SET(qq, 1, pp[0], pp[1], pp[2]);
ELL_4V_NORM(qq, qq, tmp);
ell_q_to_3m_d(rot, qq);
@@ -63,14 +58,13 @@
ELL_3M_MUL(matA, matB, rot);
TEN_M2T(ten, matA);
- fprintf(stderr, "rotated tensor = %f %f %f %f %f %f\n",
- ten[1], ten[2], ten[3], ten[4], ten[5], ten[6]);
+ fprintf(stderr, "rotated tensor = %f %f %f %f %f %f\n", ten[1], ten[2], ten[3],
+ ten[4], ten[5], ten[6]);
- ret = tenEigensolve_d(eval, evec, ten);
- fprintf(stderr, "eigensystem: %s: %g %g %g\n",
- airEnumDesc(ell_cubic_root, ret),
+ double igrt[6][7], eval[3], evec[9];
+ int ret = tenEigensolve_d(eval, evec, ten);
+ fprintf(stderr, "eigensystem: %s: %g %g %g\n", airEnumDesc(ell_cubic_root, ret),
eval[0], eval[1], eval[2]);
-
if (doK) {
tenInvariantGradientsK_d(igrt[0], igrt[1], igrt[2], ten, minnorm);
} else {
@@ -79,32 +73,25 @@
tenRotationTangents_d(igrt[3], igrt[4], igrt[5], evec);
fprintf(stderr, "invariant gradients and rotation tangents:\n");
- for (ii=0; ii<=2; ii++) {
- fprintf(stderr, " %s_%d: (norm=%g) %f %f %f %f %f %f\n",
- doK ? "K" : "R", ii+1,
- TEN_T_NORM(igrt[ii]),
- igrt[ii][1], igrt[ii][2], igrt[ii][3],
- igrt[ii][4], igrt[ii][5],
- igrt[ii][6]);
+ for (int ii = 0; ii <= 2; ii++) {
+ fprintf(stderr, " %s_%d: (norm=%g) %f %f %f %f %f %f\n", doK ? "K" : "R",
+ ii + 1, TEN_T_NORM(igrt[ii]), igrt[ii][1], igrt[ii][2], igrt[ii][3],
+ igrt[ii][4], igrt[ii][5], igrt[ii][6]);
}
- for (ii=3; ii<=5; ii++) {
- fprintf(stderr, "phi_%d: (norm=%g) %f %f %f %f %f %f\n",
- ii-2,
- TEN_T_NORM(igrt[ii]),
- igrt[ii][1], igrt[ii][2], igrt[ii][3],
- igrt[ii][4], igrt[ii][5],
- igrt[ii][6]);
+ for (int ii = 3; ii <= 5; ii++) {
+ fprintf(stderr, "phi_%d: (norm=%g) %f %f %f %f %f %f\n", ii - 2,
+ TEN_T_NORM(igrt[ii]), igrt[ii][1], igrt[ii][2], igrt[ii][3], igrt[ii][4],
+ igrt[ii][5], igrt[ii][6]);
}
fprintf(stderr, "dot products:\n");
- for (ii=0; ii<=5; ii++) {
- for (jj=ii+1; jj<=5; jj++) {
+ for (int ii = 0; ii <= 5; ii++) {
+ for (int jj = ii + 1; jj <= 5; jj++) {
fprintf(stderr, "%d,%d==%f ", ii, jj, TEN_T_DOT(igrt[ii], igrt[jj]));
}
fprintf(stderr, "\n");
}
-
airMopOkay(mop);
return 0;
}
Modified: teem/trunk/src/ten/test/odf-hist.c
===================================================================
--- teem/trunk/src/ten/test/odf-hist.c 2025-09-18 08:28:52 UTC (rev 7454)
+++ teem/trunk/src/ten/test/odf-hist.c 2025-09-18 15:34:17 UTC (rev 7455)
@@ -27,28 +27,31 @@
hestOpt *hopt = NULL;
airArray *mop;
- char *errS, *outS, *covarS;
- Nrrd *_nodf, *nvec, *nhist, *ncovar;
- unsigned int bins;
+ char *errS;
+ Nrrd *nhist, *ncovar;
size_t size[NRRD_DIM_MAX];
- float min;
mop = airMopNew();
me = argv[0];
- hestOptAdd(&hopt, "i", "odf", airTypeOther, 1, 1, &_nodf, NULL,
- "ODF volume to analyze", NULL, NULL, nrrdHestNrrd);
- hestOptAdd(&hopt, "v", "odf", airTypeOther, 1, 1, &nvec, NULL,
- "list of vectors by which odf is sampled", NULL, NULL, nrrdHestNrrd);
- hestOptAdd(&hopt, "min", "min", airTypeFloat, 1, 1, &min, "0.0",
- "ODF values below this are ignored, and per-voxel ODF is "
- "normalized to have sum 1.0. Use \"nan\" to subtract out "
- "the per-voxel min.");
- hestOptAdd(&hopt, "b", "bins", airTypeUInt, 1, 1, &bins, "128",
- "number of bins in histograms");
- hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-", "output file");
- hestOptAdd(&hopt, "co", "covariance out", airTypeString, 1, 1, &covarS, "covar.nrrd",
- "covariance output file");
+ Nrrd *_nodf;
+ hestOptAdd_1_Other(&hopt, "i", "odf", &_nodf, NULL, "ODF volume to analyze",
+ nrrdHestNrrd);
+ Nrrd *nvec;
+ hestOptAdd_1_Other(&hopt, "v", "odf", &nvec, NULL,
+ "list of vectors by which odf is sampled", nrrdHestNrrd);
+ float min;
+ hestOptAdd_1_Float(&hopt, "min", "min", &min, "0.0",
+ "ODF values below this are ignored, and per-voxel ODF is "
+ "normalized to have sum 1.0. Use \"nan\" to subtract out "
+ "the per-voxel min.");
+ unsigned int bins;
+ hestOptAdd_1_UInt(&hopt, "b", "bins", &bins, "128", "number of bins in histograms");
+ char *outS;
+ hestOptAdd_1_String(&hopt, "o", "nout", &outS, "-", "output file");
+ char *covarS;
+ hestOptAdd_1_String(&hopt, "co", "covariance out", &covarS, "covar.nrrd",
+ "covariance output file");
hestParseOrDie(hopt, argc - 1, argv + 1, NULL, me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
Modified: teem/trunk/src/ten/test/teigen.c
===================================================================
--- teem/trunk/src/ten/test/teigen.c 2025-09-18 08:28:52 UTC (rev 7454)
+++ teem/trunk/src/ten/test/teigen.c 2025-09-18 15:34:17 UTC (rev 7455)
@@ -17,46 +17,29 @@
along with this library; if not, see <https://www.gnu.org/licenses/>.
*/
-
#include "../ten.h"
const char *info = ("tests tenEigensolve_d and new stand-alone function.");
-#define ROOT_TRIPLE 2 /* ell_cubic_root_triple */
-#define ROOT_SINGLE_DOUBLE 3 /* ell_cubic_root_single_double */
-#define ROOT_THREE 4 /* ell_cubic_root_three */
+#define ROOT_TRIPLE 2 /* ell_cubic_root_triple */
+#define ROOT_SINGLE_DOUBLE 3 /* ell_cubic_root_single_double */
+#define ROOT_THREE 4 /* ell_cubic_root_three */
-#define ABS(a) (((a) > 0.0f ? (a) : -(a)))
-#define VEC_SET(v, a, b, c) \
- ((v)[0] = (a), (v)[1] = (b), (v)[2] = (c))
-#define VEC_DOT(v1, v2) \
- ((v1)[0]*(v2)[0] + (v1)[1]*(v2)[1] + (v1)[2]*(v2)[2])
-#define VEC_CROSS(v3, v1, v2) \
- ((v3)[0] = (v1)[1]*(v2)[2] - (v1)[2]*(v2)[1], \
- (v3)[1] = (v1)[2]*(v2)[0] - (v1)[0]*(v2)[2], \
- (v3)[2] = (v1)[0]*(v2)[1] - (v1)[1]*(v2)[0])
-#define VEC_ADD(v1, v2) \
- ((v1)[0] += (v2)[0], \
- (v1)[1] += (v2)[1], \
- (v1)[2] += (v2)[2])
-#define VEC_SUB(v1, v2) \
- ((v1)[0] -= (v2)[0], \
- (v1)[1] -= (v2)[1], \
- (v1)[2] -= (v2)[2])
-#define VEC_SCL(v1, s) \
- ((v1)[0] *= (s), \
- (v1)[1] *= (s), \
- (v1)[2] *= (s))
-#define VEC_LEN(v) (sqrt(VEC_DOT(v,v)))
-#define VEC_NORM(v, len) ((len) = VEC_LEN(v), VEC_SCL(v, 1.0/len))
-#define VEC_SCL_SUB(v1, s, v2) \
- ((v1)[0] -= (s)*(v2)[0], \
- (v1)[1] -= (s)*(v2)[1], \
- (v1)[2] -= (s)*(v2)[2])
-#define VEC_COPY(v1, v2) \
- ((v1)[0] = (v2)[0], \
- (v1)[1] = (v2)[1], \
- (v1)[2] = (v2)[2])
+#define ABS(a) (((a) > 0.0f ? (a) : -(a)))
+#define VEC_SET(v, a, b, c) ((v)[0] = (a), (v)[1] = (b), (v)[2] = (c))
+#define VEC_DOT(v1, v2) ((v1)[0] * (v2)[0] + (v1)[1] * (v2)[1] + (v1)[2] * (v2)[2])
+#define VEC_CROSS(v3, v1, v2) \
+ ((v3)[0] = (v1)[1] * (v2)[2] - (v1)[2] * (v2)[1], \
+ (v3)[1] = (v1)[2] * (v2)[0] - (v1)[0] * (v2)[2], \
+ (v3)[2] = (v1)[0] * (v2)[1] - (v1)[1] * (v2)[0])
+#define VEC_ADD(v1, v2) ((v1)[0] += (v2)[0], (v1)[1] += (v2)[1], (v1)[2] += (v2)[2])
+#define VEC_SUB(v1, v2) ((v1)[0] -= (v2)[0], (v1)[1] -= (v2)[1], (v1)[2] -= (v2)[2])
+#define VEC_SCL(v1, s) ((v1)[0] *= (s), (v1)[1] *= (s), (v1)[2] *= (s))
+#define VEC_LEN(v) (sqrt(VEC_DOT(v, v)))
+#define VEC_NORM(v, len) ((len) = VEC_LEN(v), VEC_SCL(v, 1.0 / len))
+#define VEC_SCL_SUB(v1, s, v2) \
+ ((v1)[0] -= (s) * (v2)[0], (v1)[1] -= (s) * (v2)[1], (v1)[2] -= (s) * (v2)[2])
+#define VEC_COPY(v1, v2) ((v1)[0] = (v2)[0], (v1)[1] = (v2)[1], (v1)[2] = (v2)[2])
/*
** All the three given vectors span only a 2D space, and this finds
@@ -65,8 +48,7 @@
** products to line up before summing.
*/
void
-nullspace1(double ret[3],
- const double r0[3], const double r1[3], const double r2[3]) {
+nullspace1(double ret[3], const double r0[3], const double r1[3], const double r2[3]) {
double crs[3];
/* ret = r0 x r1 */
@@ -96,8 +78,8 @@
** mutually vectors perpendicular to that span
*/
void
-nullspace2(double reta[3], double retb[3],
- const double r0[3], const double r1[3], const double r2[3]) {
+nullspace2(double reta[3], double retb[3], const double r0[3], const double r1[3],
+ const double r2[3]) {
double sqr[3], sum[3];
int idx;
@@ -114,14 +96,12 @@
}
/* find largest component, to get most stable expression for a
perpendicular vector */
- sqr[0] = sum[0]*sum[0];
- sqr[1] = sum[1]*sum[1];
- sqr[2] = sum[2]*sum[2];
+ sqr[0] = sum[0] * sum[0];
+ sqr[1] = sum[1] * sum[1];
+ sqr[2] = sum[2] * sum[2];
idx = 0;
- if (sqr[0] < sqr[1])
- idx = 1;
- if (sqr[idx] < sqr[2])
- idx = 2;
+ if (sqr[0] < sqr[1]) idx = 1;
+ if (sqr[idx] < sqr[2]) idx = 2;
/* reta will be perpendicular to sum */
if (0 == idx) {
VEC_SET(reta, sum[1] - sum[2], -sum[0], sum[0]);
@@ -170,10 +150,8 @@
**
*/
int
-evals(double eval[3],
- const double _M00, const double _M01, const double _M02,
- const double _M11, const double _M12,
- const double _M22) {
+evals(double eval[3], const double _M00, const double _M01, const double _M02,
+ const double _M11, const double _M12, const double _M22) {
#include "teigen-evals-A.c"
@@ -183,14 +161,11 @@
}
int
-evals_evecs(double eval[3], double evec[9],
- const double _M00, const double _M01, const double _M02,
- const double _M11, const double _M12,
- const double _M22) {
+evals_evecs(double eval[3], double evec[9], const double _M00, const double _M01,
+ const double _M02, const double _M11, const double _M12, const double _M22) {
double r0[3], r1[3], r2[3], crs[3], len, dot;
- double mean, norm, rnorm, Q, R, QQQ, D, theta,
- M00, M01, M02, M11, M12, M22;
+ double mean, norm, rnorm, Q, R, QQQ, D, theta, M00, M01, M02, M11, M12, M22;
double epsilon = 1.0E-12;
int roots;
@@ -206,7 +181,7 @@
** subtract out the eigenvalue mean (will add back to evals later);
** helps with numerical stability
*/
- mean = (M00 + M11 + M22)/3.0;
+ mean = (M00 + M11 + M22) / 3.0;
M00 -= mean;
M11 -= mean;
M22 -= mean;
@@ -215,10 +190,9 @@
** divide out L2 norm of eigenvalues (will multiply back later);
** this too seems to help with stability
*/
- norm = sqrt(M00*M00 + 2*M01*M01 + 2*M02*M02 +
- M11*M11 + 2*M12*M12 +
- M22*M22);
- rnorm = norm ? 1.0/norm : 1.0;
+ norm = sqrt(M00 * M00 + 2 * M01 * M01 + 2 * M02 * M02 + M11 * M11 + 2 * M12 * M12
+ + M22 * M22);
+ rnorm = norm ? 1.0 / norm : 1.0;
M00 *= rnorm;
M01 *= rnorm;
M02 *= rnorm;
@@ -228,21 +202,22 @@
/* this code is a mix of prior Teem code and ideas from Eberly's
"Eigensystems for 3 x 3 Symmetric Matrices (Revisited)" */
- Q = (M01*M01 + M02*M02 + M12*M12 - M00*M11 - M00*M22 - M11*M22)/3.0;
- QQQ = Q*Q*Q;
- R = (M00*M11*M22 + M02*(2*M01*M12 - M02*M11)
- - M00*M12*M12 - M01*M01*M22)/2.0;
- D = QQQ - R*R;
+ Q = (M01 * M01 + M02 * M02 + M12 * M12 - M00 * M11 - M00 * M22 - M11 * M22) / 3.0;
+ QQQ = Q * Q * Q;
+ R = (M00 * M11 * M22 + M02 * (2 * M01 * M12 - M02 * M11) - M00 * M12 * M12
+ - M01 * M01 * M22)
+ / 2.0;
+ D = QQQ - R * R;
if (D > epsilon) {
/* three distinct roots- this is the most common case */
double mm, ss, cc;
- theta = atan2(sqrt(D), R)/3.0;
+ theta = atan2(sqrt(D), R) / 3.0;
mm = sqrt(Q);
ss = sin(theta);
cc = cos(theta);
- eval[0] = 2*mm*cc;
- eval[1] = mm*(-cc + sqrt(3.0)*ss);
- eval[2] = mm*(-cc - sqrt(3.0)*ss);
+ eval[0] = 2 * mm * cc;
+ eval[1] = mm * (-cc + sqrt(3.0) * ss);
+ eval[2] = mm * (-cc - sqrt(3.0) * ss);
roots = ROOT_THREE;
/* else D is near enough to zero */
} else if (R < -epsilon || epsilon < R) {
@@ -250,13 +225,13 @@
/* one double root and one single root */
U = airCbrt(R); /* cube root function */
if (U > 0) {
- eval[0] = 2*U;
+ eval[0] = 2 * U;
eval[1] = -U;
eval[2] = -U;
} else {
eval[0] = -U;
eval[1] = -U;
- eval[2] = 2*U;
+ eval[2] = 2 * U;
}
roots = ROOT_SINGLE_DOUBLE;
} else {
@@ -271,57 +246,76 @@
VEC_SET(r1, M01, 0.0, M12);
VEC_SET(r2, M02, M12, 0.0);
if (ROOT_THREE == roots) {
- r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
- nullspace1(evec+0, r0, r1, r2);
- r0[0] = M00 - eval[1]; r1[1] = M11 - eval[1]; r2[2] = M22 - eval[1];
- nullspace1(evec+3, r0, r1, r2);
- r0[0] = M00 - eval[2]; r1[1] = M11 - eval[2]; r2[2] = M22 - eval[2];
- nullspace1(evec+6, r0, r1, r2);
+ r0[0] = M00 - eval[0];
+ r1[1] = M11 - eval[0];
+ r2[2] = M22 - eval[0];
+ nullspace1(evec + 0, r0, r1, r2);
+ r0[0] = M00 - eval[1];
+ r1[1] = M11 - eval[1];
+ r2[2] = M22 - eval[1];
+ nullspace1(evec + 3, r0, r1, r2);
+ r0[0] = M00 - eval[2];
+ r1[1] = M11 - eval[2];
+ r2[2] = M22 - eval[2];
+ nullspace1(evec + 6, r0, r1, r2);
} else if (ROOT_SINGLE_DOUBLE == roots) {
if (eval[1] == eval[2]) {
/* one big (eval[0]) , two small (eval[1,2]) */
- r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
- nullspace1(evec+0, r0, r1, r2);
- r0[0] = M00 - eval[1]; r1[1] = M11 - eval[1]; r2[2] = M22 - eval[1];
- nullspace2(evec+3, evec+6, r0, r1, r2);
- }
- else {
+ r0[0] = M00 - eval[0];
+ r1[1] = M11 - eval[0];
+ r2[2] = M22 - eval[0];
+ nullspace1(evec + 0, r0, r1, r2);
+ r0[0] = M00 - eval[1];
+ r1[1] = M11 - eval[1];
+ r2[2] = M22 - eval[1];
+ nullspace2(evec + 3, evec + 6, r0, r1, r2);
+ } else {
/* two big (eval[0,1]), one small (eval[2]) */
- r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
- nullspace2(evec+0, evec+3, r0, r1, r2);
- r0[0] = M00 - eval[2]; r1[1] = M11 - eval[2]; r2[2] = M22 - eval[2];
- nullspace1(evec+6, r0, r1, r2);
+ r0[0] = M00 - eval[0];
+ r1[1] = M11 - eval[0];
+ r2[2] = M22 - eval[0];
+ nullspace2(evec + 0, evec + 3, r0, r1, r2);
+ r0[0] = M00 - eval[2];
+ r1[1] = M11 - eval[2];
+ r2[2] = M22 - eval[2];
+ nullspace1(evec + 6, r0, r1, r2);
}
} else {
/* ROOT_TRIPLE == roots; use any basis for eigenvectors */
- VEC_SET(evec+0, 1, 0, 0);
- VEC_SET(evec+3, 0, 1, 0);
- VEC_SET(evec+6, 0, 0, 1);
+ VEC_SET(evec + 0, 1, 0, 0);
+ VEC_SET(evec + 3, 0, 1, 0);
+ VEC_SET(evec + 6, 0, 0, 1);
}
/* we always make sure its really orthonormal; keeping fixed the
eigenvector associated with the largest-magnitude eigenvalue */
if (ABS(eval[0]) > ABS(eval[2])) {
/* normalize evec+0 but don't move it */
- VEC_NORM(evec+0, len);
- dot = VEC_DOT(evec+0, evec+3); VEC_SCL_SUB(evec+3, dot, evec+0);
- VEC_NORM(evec+3, len);
- dot = VEC_DOT(evec+0, evec+6); VEC_SCL_SUB(evec+6, dot, evec+0);
- dot = VEC_DOT(evec+3, evec+6); VEC_SCL_SUB(evec+6, dot, evec+3);
- VEC_NORM(evec+6, len);
+ VEC_NORM(evec + 0, len);
+ dot = VEC_DOT(evec + 0, evec + 3);
+ VEC_SCL_SUB(evec + 3, dot, evec + 0);
+ VEC_NORM(evec + 3, len);
+ dot = VEC_DOT(evec + 0, evec + 6);
+ VEC_SCL_SUB(evec + 6, dot, evec + 0);
+ dot = VEC_DOT(evec + 3, evec + 6);
+ VEC_SCL_SUB(evec + 6, dot, evec + 3);
+ VEC_NORM(evec + 6, len);
} else {
/* normalize evec+6 but don't move it */
- VEC_NORM(evec+6, len);
- dot = VEC_DOT(evec+6, evec+3); VEC_SCL_SUB(evec+3, dot, evec+6);
- VEC_NORM(evec+3, len);
- dot = VEC_DOT(evec+3, evec+0); VEC_SCL_SUB(evec+0, dot, evec+3);
- dot = VEC_DOT(evec+6, evec+0); VEC_SCL_SUB(evec+0, dot, evec+6);
- VEC_NORM(evec+0, len);
+ VEC_NORM(evec + 6, len);
+ dot = VEC_DOT(evec + 6, evec + 3);
+ VEC_SCL_SUB(evec + 3, dot, evec + 6);
+ VEC_NORM(evec + 3, len);
+ dot = VEC_DOT(evec + 3, evec + 0);
+ VEC_SCL_SUB(evec + 0, dot, evec + 3);
+ dot = VEC_DOT(evec + 6, evec + 0);
+ VEC_SCL_SUB(evec + 0, dot, evec + 6);
+ VEC_NORM(evec + 0, len);
}
/* to be nice, make it right-handed */
- VEC_CROSS(crs, evec+0, evec+3);
- if (0 > VEC_DOT(crs, evec+6)) {
- VEC_SCL(evec+6, -1);
+ VEC_CROSS(crs, evec + 0, evec + 3);
+ if (0 > VEC_DOT(crs, evec + 6)) {
+ VEC_SCL(evec + 6, -1);
}
/* multiply back by eigenvalue L2 norm */
@@ -337,7 +331,6 @@
return roots;
}
-
void
testeigen(double tt[7], double eval[3], double evec[9]) {
double mat[9], dot[3], cross[3];
@@ -345,17 +338,14 @@
TEN_T2M(mat, tt);
printf("evals %g %g %g\n", eval[0], eval[1], eval[2]);
- printf("evec0 (%g) %g %g %g\n",
- ELL_3V_LEN(evec + 0), evec[0], evec[1], evec[2]);
- printf("evec1 (%g) %g %g %g\n",
- ELL_3V_LEN(evec + 3), evec[3], evec[4], evec[5]);
- printf("evec2 (%g) %g %g %g\n",
- ELL_3V_LEN(evec + 6), evec[6], evec[7], evec[8]);
+ printf("evec0 (%g) %g %g %g\n", ELL_3V_LEN(evec + 0), evec[0], evec[1], evec[2]);
+ printf("evec1 (%g) %g %g %g\n", ELL_3V_LEN(evec + 3), evec[3], evec[4], evec[5]);
+ printf("evec2 (%g) %g %g %g\n", ELL_3V_LEN(evec + 6), evec[6], evec[7], evec[8]);
printf("Mv - lv: (len) X Y Z (should be ~zeros)\n");
- for (ii=0; ii<3; ii++) {
+ for (ii = 0; ii < 3; ii++) {
double uu[3], vv[3], dd[3];
- ELL_3MV_MUL(uu, mat, evec + 3*ii);
- ELL_3V_SCALE(vv, eval[ii], evec + 3*ii);
+ ELL_3MV_MUL(uu, mat, evec + 3 * ii);
+ ELL_3V_SCALE(vv, eval[ii], evec + 3 * ii);
ELL_3V_SUB(dd, uu, vv);
printf("%d: (%g) %g %g %g\n", ii, ELL_3V_LEN(dd), dd[0], dd[1], dd[2]);
}
@@ -362,10 +352,9 @@
dot[0] = ELL_3V_DOT(evec + 0, evec + 3);
dot[1] = ELL_3V_DOT(evec + 0, evec + 6);
dot[2] = ELL_3V_DOT(evec + 3, evec + 6);
- printf("pairwise dots: (%g) %g %g %g\n",
- ELL_3V_LEN(dot), dot[0], dot[1], dot[2]);
- ELL_3V_CROSS(cross, evec+0, evec+3);
- printf("right-handed: %g\n", ELL_3V_DOT(evec+6, cross));
+ printf("pairwise dots: (%g) %g %g %g\n", ELL_3V_LEN(dot), dot[0], dot[1], dot[2]);
+ ELL_3V_CROSS(cross, evec + 0, evec + 3);
+ printf("right-handed: %g\n", ELL_3V_DOT(evec + 6, cross));
return;
}
@@ -372,30 +361,29 @@
int
main(int argc, const char *argv[]) {
const char *me;
- hestOpt *hopt=NULL;
+ hestOpt *hopt = NULL;
airArray *mop;
- double _tt[6], tt[7], ss, pp[3], qq[4], rot[9], mat1[9], mat2[9], tmp,
- evalA[3], evecA[9], evalB[3], evecB[9];
- int roots;
-
mop = airMopNew();
me = argv[0];
- hestOptAdd(&hopt, NULL, "m00 m01 m02 m11 m12 m22",
- airTypeDouble, 6, 6, _tt, NULL, "symmtric matrix coeffs");
- hestOptAdd(&hopt, "p", "vec", airTypeDouble, 3, 3, pp, "0 0 0",
- "rotation as P vector");
- hestOptAdd(&hopt, "s", "scl", airTypeDouble, 1, 1, &ss, "1.0",
- "scaling");
- hestParseOrDie(hopt, argc-1, argv+1, NULL,
- me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
+ double _tt[6];
+ hestOptAdd_N_Double(&hopt, NULL, "m00 m01 m02 m11 m12 m22", 6, _tt, NULL,
+ "symmtric matrix coeffs");
+ double pp[3];
+ hestOptAdd_3_Double(&hopt, "p", "vec", pp, "0 0 0", "rotation as P vector");
+ double ss;
+ hestOptAdd_1_Double(&hopt, "s", "scl", &ss, "1.0", "scaling");
+ hestParseOrDie(hopt, argc - 1, argv + 1, NULL, me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
+ double tt[7];
ELL_6V_COPY(tt + 1, _tt);
tt[0] = 1.0;
TEN_T_SCALE(tt, ss, tt);
+ double qq[4], rot[9], mat1[9], mat2[9], tmp, evalA[3], evecA[9], evalB[3], evecB[9];
+ int roots;
ELL_4V_SET(qq, 1, pp[0], pp[1], pp[2]);
ELL_4V_NORM(qq, qq, tmp);
ell_q_to_3m_d(rot, qq);
@@ -410,8 +398,8 @@
ell_3m_mul_d(mat1, mat2, rot);
TEN_M2T(tt, mat1);
- printf("input matrix = \n %g %g %g\n %g %g\n %g\n",
- tt[1], tt[2], tt[3], tt[4], tt[5], tt[6]);
+ printf("input matrix = \n %g %g %g\n %g %g\n %g\n", tt[1], tt[2], tt[3], tt[4], tt[5],
+ tt[6]);
printf("================== tenEigensolve_d ==================\n");
roots = tenEigensolve_d(evalA, evecA, tt);
@@ -430,10 +418,9 @@
printf("================== new eigensolve ==================\n");
roots = evals(evalB, tt[1], tt[2], tt[3], tt[4], tt[5], tt[6]);
- printf("%s roots: %g %g %g\n", airEnumStr(ell_cubic_root, roots),
- evalB[0], evalB[1], evalB[2]);
- roots = evals_evecs(evalB, evecB,
- tt[1], tt[2], tt[3], tt[4], tt[5], tt[6]);
+ printf("%s roots: %g %g %g\n", airEnumStr(ell_cubic_root, roots), evalB[0], evalB[1],
+ evalB[2]);
+ roots = evals_evecs(evalB, evecB, tt[1], tt[2], tt[3], tt[4], tt[5], tt[6]);
printf("%s roots\n", airEnumStr(ell_cubic_root, roots));
testeigen(tt, evalB, evecB);
Modified: teem/trunk/src/ten/test/tem.c
===================================================================
--- teem/trunk/src/ten/test/tem.c 2025-09-18 08:28:52 UTC (rev 7454)
+++ teem/trunk/src/ten/test/tem.c 2025-09-18 15:34:17 UTC (rev 7455)
@@ -17,7 +17,6 @@
along with this library; if not, see <https://www.gnu.org/licenses/>.
*/
-
#include "../ten.h"
const char *info = ("Test EM bimodal histogram fitting.");
@@ -26,11 +25,9 @@
main(int argc, const char *argv[]) {
const char *me;
char *err;
- hestOpt *hopt=NULL;
+ hestOpt *hopt = NULL;
airArray *mop;
- Nrrd *nhisto;
tenEMBimodalParm *biparm;
- double minprob[2];
mop = airMopNew();
me = argv[0];
@@ -38,18 +35,16 @@
biparm = tenEMBimodalParmNew();
airMopAdd(mop, biparm, (airMopper)tenEMBimodalParmNix, airMopAlways);
- hestOptAdd(&hopt, NULL, "histogram", airTypeOther, 1, 1, &nhisto, NULL,
- "The 1-D histogram to analyize", NULL, NULL, nrrdHestNrrd);
- hestOptAdd(&hopt, "ts", "two stage", airTypeInt, 0, 0,
- &(biparm->twoStage), NULL,
- "use two-stage processing");
- hestOptAdd(&hopt, "v", "verbose", airTypeInt, 1, 1, &(biparm->verbose), "1",
- "verbosity level");
- hestOptAdd(&hopt, "mp", "minprob 1,2", airTypeDouble, 2, 2, minprob, "0 0",
- "minimum significant posterior probabilies, for first and "
- "second stages");
- hestParseOrDie(hopt, argc-1, argv+1, NULL,
- me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
+ Nrrd *nhisto;
+ hestOptAdd_1_Other(&hopt, NULL, "histogram", &nhisto, NULL,
+ "The 1-D histogram to analyize", nrrdHestNrrd);
+ hestOptAdd_Flag(&hopt, "ts", &(biparm->twoStage), "use two-stage processing");
+ hestOptAdd_1_Int(&hopt, "v", "verbose", &(biparm->verbose), "1", "verbosity level");
+ double minprob[2];
+ hestOptAdd_2_Double(&hopt, "mp", "minprob 1,2", minprob, "0 0",
+ "minimum significant posterior probabilies, for first and "
+ "second stages");
+ hestParseOrDie(hopt, argc - 1, argv + 1, NULL, me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
@@ -63,14 +58,13 @@
return 1;
}
fprintf(stderr, "%s: bimodal histogram fit\n", me);
- fprintf(stderr, "material 1 (%g%%): mean = %g, stdv = %g\n",
- 100*(biparm->fraction1), biparm->mean1, biparm->stdv1);
+ fprintf(stderr, "material 1 (%g%%): mean = %g, stdv = %g\n", 100 * (biparm->fraction1),
+ biparm->mean1, biparm->stdv1);
fprintf(stderr, "material 2 (%g%%): mean = %g, stdv = %g\n",
- 100*(1 - biparm->fraction1), biparm->mean2, biparm->stdv2);
- fprintf(stderr, " ---> optimal threshold = %g (confidence = %g)\n",
- biparm->threshold, biparm->confidence);
+ 100 * (1 - biparm->fraction1), biparm->mean2, biparm->stdv2);
+ fprintf(stderr, " ---> optimal threshold = %g (confidence = %g)\n", biparm->threshold,
+ biparm->confidence);
airMopOkay(mop);
return 0;
}
-
Modified: teem/trunk/src/ten/test/tg.c
===================================================================
--- teem/trunk/src/ten/test/tg.c 2025-09-18 08:28:52 UTC (rev 7454)
+++ teem/trunk/src/ten/test/tg.c 2025-09-18 15:34:17 UTC (rev 7455)
@@ -54,32 +54,33 @@
int
main(int argc, const char *argv[]) {
const char *me;
- char *err, *outS;
+ char *err;
hestOpt *hopt = NULL;
airArray *mop;
- int xi, yi, zi, samp;
- float *tdata;
- double clp[2], xyz[3], q[4], len;
- double mD[9], mRF[9], mRI[9], mT[9];
- Nrrd *nten;
mop = airMopNew();
me = argv[0];
- hestOptAdd(&hopt, "n", "# samples", airTypeInt, 1, 1, &samp, "4",
- "number of samples along each edge of cube");
- hestOptAdd(&hopt, "c", "cl cp", airTypeDouble, 2, 2, clp, NULL,
- "shape of tensor to use; \"cl\" and \"cp\" are cl1 "
- "and cp1 values, both in [0.0,1.0]");
- hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-",
- "output file to save tensors into");
+ int samp;
+ hestOptAdd_1_Int(&hopt, "n", "# samples", &samp, "4",
+ "number of samples along each edge of cube");
+ double clp[2];
+ hestOptAdd_2_Double(&hopt, "c", "cl cp", clp, NULL,
+ "shape of tensor to use; \"cl\" and \"cp\" are cl1 "
+ "and cp1 values, both in [0.0,1.0]");
+ char *outS;
+ hestOptAdd_1_String(&hopt, "o", "nout", &outS, "-",
+ "output file to save tensors into");
hestParseOrDie(hopt, argc - 1, argv + 1, NULL, me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
- nten = nrrdNew();
+ Nrrd *nten = nrrdNew();
airMopAdd(mop, nten, (airMopper)nrrdNuke, airMopAlways);
+ float *tdata;
+ double xyz[3], q[4], len;
+ double mD[9], mRF[9], mRI[9], mT[9];
_clp2xyz(xyz, clp);
fprintf(stderr, "%s: want evals = %g %g %g\n", me, xyz[0], xyz[1], xyz[2]);
@@ -93,9 +94,9 @@
ELL_3M_IDENTITY_SET(mD);
ELL_3M_DIAG_SET(mD, xyz[0], xyz[1], xyz[2]);
tdata = (float *)nten->data;
- for (zi = 0; zi < samp; zi++) {
- for (yi = 0; yi < samp; yi++) {
- for (xi = 0; xi < samp; xi++) {
+ for (int zi = 0; zi < samp; zi++) {
+ for (int yi = 0; yi < samp; yi++) {
+ for (int xi = 0; xi < samp; xi++) {
q[0] = 1.0;
q[1] = AIR_AFFINE(-0.5, (float)xi, samp - 0.5, -1, 1);
q[2] = AIR_AFFINE(-0.5, (float)yi, samp - 0.5, -1, 1);
Modified: teem/trunk/src/ten/test/tqgl.c
===================================================================
--- teem/trunk/src/ten/test/tqgl.c 2025-09-18 08:28:52 UTC (rev 7454)
+++ teem/trunk/src/ten/test/tqgl.c 2025-09-18 15:34:17 UTC (rev 7455)
@@ -17,22 +17,17 @@
along with this library; if not, see <https://www.gnu.org/licenses/>.
*/
-
#include "../ten.h"
#include "../privateTen.h"
/* BAD gordon */
extern double _tenQGL_Kdist(const double RThZA[3], const double RThZB[3]);
-extern void _tenQGL_Klog(double klog[3],
- const double RThZA[3], const double RThZB[3]);
-extern void _tenQGL_Kexp(double RThZB[3],
- const double RThZA[3], const double klog[3]);
+extern void _tenQGL_Klog(double klog[3], const double RThZA[3], const double RThZB[3]);
+extern void _tenQGL_Kexp(double RThZB[3], const double RThZA[3], const double klog[3]);
extern double _tenQGL_Rdist(const double RThPhA[3], const double RThPhB[3]);
-extern void _tenQGL_Rlog(double rlog[3],
- const double RThPhA[3], const double RThPhB[3]);
-extern void _tenQGL_Rexp(double RThPhB[3],
- const double RThPhA[3], const double rlog[3]);
+extern void _tenQGL_Rlog(double rlog[3], const double RThPhA[3], const double RThPhB[3]);
+extern void _tenQGL_Rexp(double RThPhB[3], const double RThPhA[3], const double rlog[3]);
/* normalized gradients of k or r invariants, in XYZ space,
to help determine if path is really a loxodrome */
@@ -40,8 +35,7 @@
kgrads(double grad[3][3], const double eval[3]) {
double rtz[3];
- tenTripleConvertSingle_d(rtz, tenTripleTypeRThetaZ,
- eval, tenTripleTypeEigenvalue);
+ tenTripleConvertSingle_d(rtz, tenTripleTypeRThetaZ, eval, tenTripleTypeEigenvalue);
ELL_3V_SET(grad[0], cos(rtz[1]), sin(rtz[1]), 0);
ELL_3V_SET(grad[1], -sin(rtz[1]), cos(rtz[1]), 0);
@@ -52,17 +46,13 @@
rgrads(double grad[3][3], const double eval[3]) {
double rtp[3];
- tenTripleConvertSingle_d(rtp, tenTripleTypeRThetaPhi,
- eval, tenTripleTypeEigenvalue);
+ tenTripleConvertSingle_d(rtp, tenTripleTypeRThetaPhi, eval, tenTripleTypeEigenvalue);
- ELL_3V_SET(grad[0],
- cos(rtp[1])*sin(rtp[2]),
- sin(rtp[1])*sin(rtp[2]),
- cos(rtp[2]));
+ ELL_3V_SET(grad[0], cos(rtp[1]) * sin(rtp[2]), sin(rtp[1]) * sin(rtp[2]), cos(rtp[2]));
ELL_3V_SET(grad[1], -sin(rtp[1]), cos(rtp[1]), 0);
ELL_3V_SET(grad[2],
- cos(rtp[1])*cos(rtp[2]),
- sin(rtp[1])*cos(rtp[2]),
+ cos(rtp[1]) * cos(rtp[2]),
+ sin(rtp[1]) * cos(rtp[2]),
-sin(rtp[2]));
}
@@ -73,67 +63,47 @@
int
main(int argc, const char *argv[]) {
const char *me;
- hestOpt *hopt=NULL;
- airArray *mop;
+ hestOpt *hopt = NULL;
- double tripA[3], tripB[3], evalA[3], evalB[3],
- rt_A[3], rt_B[3], trip[3], eval[3], lasteval[3], lastxyz[3],
- logAB[3], ndist;
- int ittype, ottype, ptype, rttype;
- unsigned int NN, ii;
- tenInterpParm *tip;
-
- void (*interp)(double oeval[3], const double evalA[3],
- const double evalB[3], const double tt);
- double (*qdist)(const double RTh_A[3], const double RTh_B[3]);
- void (*qlog)(double klog[3],
- const double RThZA[3], const double RThZB[3]);
- void (*qexp)(double RThZB[3],
- const double RThZA[3], const double klog[3]);
- void (*grads)(double grad[3][3], const double eval[3]);
-
me = argv[0];
- mop = airMopNew();
- tip = tenInterpParmNew();
+ airArray *mop = airMopNew();
+ tenInterpParm *tip = tenInterpParmNew();
airMopAdd(mop, tip, (airMopper)tenInterpParmNix, airMopAlways);
- hestOptAdd(&hopt, "a", "start", airTypeDouble, 3, 3, tripA, NULL,
- "start triple of values");
- hestOptAdd(&hopt, "b", "end", airTypeDouble, 3, 3, tripB, NULL,
- "end triple of values");
- hestOptAdd(&hopt, "it", "type", airTypeEnum, 1, 1, &ittype, NULL,
- "type of given start and end triples", NULL, tenTripleType);
- hestOptAdd(&hopt, "ot", "type", airTypeEnum, 1, 1, &ottype, NULL,
- "type of triples for output", NULL, tenTripleType);
- hestOptAdd(&hopt, "p", "type", airTypeEnum, 1, 1, &ptype, NULL,
- "type of path interpolation", NULL, tenInterpType);
- hestOptAdd(&hopt, "n", "# steps", airTypeUInt, 1, 1, &NN, "100",
- "number of steps along path");
+ double tripA[3];
+ hestOptAdd_3_Double(&hopt, "a", "start", tripA, NULL, "start triple of values");
+ double tripB[3];
+ hestOptAdd_3_Double(&hopt, "b", "end", tripB, NULL, "end triple of values");
+ int ittype;
+ hestOptAdd_1_Enum(&hopt, "it", "type", &ittype, NULL,
+ "type of given start and end triples", tenTripleType);
+ int ottype;
+ hestOptAdd_1_Enum(&hopt, "ot", "type", &ottype, NULL, "type of triples for output",
+ tenTripleType);
+ int ptype;
+ hestOptAdd_1_Enum(&hopt, "p", "type", &ptype, NULL, "type of path interpolation",
+ tenInterpType);
+ unsigned int NN;
+ hestOptAdd_1_UInt(&hopt, "n", "# steps", &NN, "100", "number of steps along path");
- hestOptAdd(&hopt, "v", "verbosity", airTypeInt, 1, 1,
- &(tip->verbose), "0", "verbosity");
- hestOptAdd(&hopt, "s", "stepsize", airTypeDouble, 1, 1,
- &(tip->convStep), "1", "step size in update");
- hestOptAdd(&hopt, "r", "recurse", airTypeInt, 0, 0,
- &(tip->enableRecurse), NULL,
+ hestOptAdd(&hopt, "v", "verbosity", airTypeInt, 1, 1, &(tip->verbose), "0",
+ "verbosity");
+ hestOptAdd(&hopt, "s", "stepsize", airTypeDouble, 1, 1, &(tip->convStep), "1",
+ "step size in update");
+ hestOptAdd(&hopt, "r", "recurse", airTypeInt, 0, 0, &(tip->enableRecurse), NULL,
"enable recursive solution, when useful");
- hestOptAdd(&hopt, "mn", "minnorm", airTypeDouble, 1, 1,
- &(tip->minNorm), "0.000001",
+ hestOptAdd(&hopt, "mn", "minnorm", airTypeDouble, 1, 1, &(tip->minNorm), "0.000001",
"minnorm of something");
- hestOptAdd(&hopt, "mi", "maxiter", airTypeUInt, 1, 1,
- &(tip->maxIter), "0",
+ hestOptAdd(&hopt, "mi", "maxiter", airTypeUInt, 1, 1, &(tip->maxIter), "0",
"if non-zero, max # iterations for computation");
- hestOptAdd(&hopt, "c", "conv", airTypeDouble, 1, 1,
- &(tip->convEps), "0.0001",
+ hestOptAdd(&hopt, "c", "conv", airTypeDouble, 1, 1, &(tip->convEps), "0.0001",
"convergence threshold of length fraction");
- hestParseOrDie(hopt, argc-1, argv+1, NULL,
- me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
+ hestParseOrDie(hopt, argc - 1, argv + 1, NULL, me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
- if (!( tenInterpTypeQuatGeoLoxK == ptype
- || tenInterpTypeQuatGeoLoxR == ptype )) {
+ if (!(tenInterpTypeQuatGeoLoxK == ptype || tenInterpTypeQuatGeoLoxR == ptype)) {
fprintf(stderr, "%s: need type %s or %s, not %s\n", me,
airEnumStr(tenInterpType, tenInterpTypeQuatGeoLoxK),
airEnumStr(tenInterpType, tenInterpTypeQuatGeoLoxR),
@@ -142,6 +112,13 @@
return 1;
}
+ void (*interp)(double oeval[3], const double evalA[3], const double evalB[3],
+ const double tt);
+ double (*qdist)(const double RTh_A[3], const double RTh_B[3]);
+ void (*qlog)(double klog[3], const double RThZA[3], const double RThZB[3]);
+ void (*qexp)(double RThZB[3], const double RThZA[3], const double klog[3]);
+ void (*grads)(double grad[3][3], const double eval[3]);
+ int rttype;
if (tenInterpTypeQuatGeoLoxK == ptype) {
interp = tenQGLInterpTwoEvalK;
qdist = _tenQGL_Kdist;
@@ -159,30 +136,28 @@
}
fprintf(stderr, "%s: (%s) %f %f %f \n--%s--> %f %f %f\n", me,
- airEnumStr(tenTripleType, ittype),
- tripA[0], tripA[1], tripA[2],
- airEnumStr(tenInterpType, ptype),
- tripB[0], tripB[1], tripB[2]);
+ airEnumStr(tenTripleType, ittype), tripA[0], tripA[1], tripA[2],
+ airEnumStr(tenInterpType, ptype), tripB[0], tripB[1], tripB[2]);
+ double evalA[3], evalB[3], rt_A[3], rt_B[3], trip[3], eval[3], lasteval[3], lastxyz[3],
+ logAB[3];
tenTripleConvertSingle_d(evalA, tenTripleTypeEigenvalue, tripA, ittype);
tenTripleConvertSingle_d(evalB, tenTripleTypeEigenvalue, tripB, ittype);
tenTripleConvertSingle_d(rt_A, rttype, tripA, ittype);
tenTripleConvertSingle_d(rt_B, rttype, tripB, ittype);
- ndist = 0;
+ double ndist = 0;
ELL_3V_SET(lasteval, AIR_NAN, AIR_NAN, AIR_NAN);
ELL_3V_SET(lastxyz, AIR_NAN, AIR_NAN, AIR_NAN);
qlog(logAB, rt_A, rt_B);
- fprintf(stderr, "%s: log = %g %g %g (%g)\n", me,
- logAB[0], logAB[1], logAB[2], ELL_3V_LEN(logAB));
- for (ii=0; ii<NN; ii++) {
+ fprintf(stderr, "%s: log = %g %g %g (%g)\n", me, logAB[0], logAB[1], logAB[2],
+ ELL_3V_LEN(logAB));
+ for (unsigned int ii = 0; ii < NN; ii++) {
double tt, xyz[3], dot[3], ll[3], prayRT[3], prayO[3];
- tt = AIR_AFFINE(0, ii, NN-1, 0.0, 1.0);
+ tt = AIR_AFFINE(0, ii, NN - 1, 0.0, 1.0);
interp(eval, evalA, evalB, tt);
- tenTripleConvertSingle_d(trip, ottype,
- eval, tenTripleTypeEigenvalue);
- tenTripleConvertSingle_d(xyz, tenTripleTypeXYZ,
- eval, tenTripleTypeEigenvalue);
+ tenTripleConvertSingle_d(trip, ottype, eval, tenTripleTypeEigenvalue);
+ tenTripleConvertSingle_d(xyz, tenTripleTypeXYZ, eval, tenTripleTypeEigenvalue);
ELL_3V_SCALE(ll, tt, logAB);
qexp(prayRT, rt_A, ll);
tenTripleConvertSingle_d(prayO, ottype, prayRT, rttype);
@@ -198,17 +173,13 @@
} else {
ELL_3V_SET(dot, 0, 0, 0);
}
- printf("%03u %g %g %g %g %g %g 00 %g %g %g\n", ii,
- trip[0], prayO[0],
- trip[1], prayO[1],
- trip[2], prayO[2],
- dot[0], dot[1], dot[2]);
+ printf("%03u %g %g %g %g %g %g 00 %g %g %g\n", ii, trip[0], prayO[0],
+ trip[1], prayO[1], trip[2], prayO[2], dot[0], dot[1], dot[2]);
ELL_3V_COPY(lasteval, eval);
ELL_3V_COPY(lastxyz, xyz);
}
- fprintf(stderr, "%s: dist %g =?= %g\n", me,
- qdist(rt_A, rt_B), ndist);
+ fprintf(stderr, "%s: dist %g =?= %g\n", me, qdist(rt_A, rt_B), ndist);
airMopOkay(mop);
return 0;
Modified: teem/trunk/src/ten/test/tt.c
===================================================================
--- teem/trunk/src/ten/test/tt.c 2025-09-18 08:28:52 UTC (rev 7454)
+++ teem/trunk/src/ten/test/tt.c 2025-09-18 15:34:17 UTC (rev 7455)
@@ -103,45 +103,44 @@
int
main(int argc, const char *argv[]) {
- const char *me;
- char *err, *outS;
+ const char *me = argv[0];
+ char *err;
hestOpt *hopt = NULL;
- airArray *mop;
- int sx, sy, xi, yi, samp, version, whole, right;
- float *tdata;
- double p[3], xyz[3], q[4], len, hackcp = 0, maxca;
- double ca, cp, mD[9], mRF[9], mRI[9], mT[9], hack;
- Nrrd *nten;
- mop = airMopNew();
-
- me = argv[0];
- hestOptAdd(&hopt, "n", "# samples", airTypeInt, 1, 1, &samp, "4",
- "number of glyphs along each edge of triangle");
- hestOptAdd(&hopt, "p", "x y z", airTypeDouble, 3, 3, p, NULL,
- "location in quaternion quotient space");
- hestOptAdd(&hopt, "ca", "max ca", airTypeDouble, 1, 1, &maxca, "0.8",
- "maximum ca to use at bottom edge of triangle");
- hestOptAdd(&hopt, "r", NULL, airTypeInt, 0, 0, &right, NULL,
- "sample a right-triangle-shaped region, instead of "
- "a roughly equilateral triangle. ");
- hestOptAdd(&hopt, "w", NULL, airTypeInt, 0, 0, &whole, NULL,
- "sample the whole triangle of constant trace, "
- "instead of just the "
- "sixth of it in which the eigenvalues have the "
- "traditional sorted order. ");
- hestOptAdd(&hopt, "hack", "hack", airTypeDouble, 1, 1, &hack, "0.04",
- "this is a hack");
- hestOptAdd(&hopt, "v", "version", airTypeInt, 1, 1, &version, "1",
- "which version of the Westin metrics to use to parameterize "
- "triangle; \"1\" for ISMRM 97, \"2\" for MICCAI 99");
- hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-",
- "output file to save tensors into");
+ airArray *mop = airMopNew();
+ int samp;
+ hestOptAdd_1_Int(&hopt, "n", "# samples", &samp, "4",
+ "number of glyphs along each edge of triangle");
+ double p[3];
+ hestOptAdd_3_Double(&hopt, "p", "x y z", p, NULL,
+ "location in quaternion quotient space");
+ double maxca;
+ hestOptAdd_1_Double(&hopt, "ca", "max ca", &maxca, "0.8",
+ "maximum ca to use at bottom edge of triangle");
+ int right;
+ hestOptAdd_Flag(&hopt, "r", &right,
+ "sample a right-triangle-shaped region, instead of "
+ "a roughly equilateral triangle. ");
+ int whole;
+ hestOptAdd_Flag(&hopt, "w", &whole,
+ "sample the whole triangle of constant trace, "
+ "instead of just the "
+ "sixth of it in which the eigenvalues have the "
+ "traditional sorted order. ");
+ double hack;
+ hestOptAdd_1_Double(&hopt, "hack", "hack", &hack, "0.04", "this is a hack");
+ int version;
+ hestOptAdd_1_Int(&hopt, "v", "version", &version, "1",
+ "which version of the Westin metrics to use to parameterize "
+ "triangle; \"1\" for ISMRM 97, \"2\" for MICCAI 99");
+ char *outS;
+ hestOptAdd_1_String(&hopt, "o", "nout", &outS, "-",
+ "output file to save tensors into");
hestParseOrDie(hopt, argc - 1, argv + 1, NULL, me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
- nten = nrrdNew();
+ Nrrd *nten = nrrdNew();
airMopAdd(mop, nten, (airMopper)nrrdNuke, airMopAlways);
if (!(1 == version || 2 == version)) {
@@ -149,6 +148,7 @@
airMopError(mop);
return 1;
}
+ int sx, sy;
if (right) {
sx = samp;
sy = (int)(1.0 * samp / sqrt(3.0));
@@ -163,6 +163,9 @@
airMopError(mop);
return 1;
}
+ float *tdata;
+ double xyz[3], q[4], len, hackcp = 0;
+ double ca, cp, mD[9], mRF[9], mRI[9], mT[9];
q[0] = 1.0;
q[1] = p[0];
q[2] = p[1];
@@ -252,7 +255,7 @@
nten->axis[2].spacing = (sx - 1) / (sqrt(3.0) * (sy - 1));
nten->axis[3].spacing = 1;
} else {
- for (yi = 0; yi < samp; yi++) {
+ for (int yi = 0; yi < samp; yi++) {
if (whole) {
ca = AIR_AFFINE(0, yi, samp - 1, 0.0, 1.0);
} else {
@@ -259,7 +262,7 @@
ca = AIR_AFFINE(0, yi, samp - 1, hack, maxca);
hackcp = AIR_AFFINE(0, yi, samp - 1, hack, 0);
}
- for (xi = 0; xi <= yi; xi++) {
+ for (int xi = 0; xi <= yi; xi++) {
if (whole) {
cp = AIR_AFFINE(0, xi, samp - 1, 0.0, 1.0);
} else {
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
