| File: | Testing/meet/probeSS.c |
| Location: | line 1082, column 9 |
| Description: | Value stored to 'ktmp' is never read |
| 1 | /* |
| 2 | Teem: Tools to process and visualize scientific data and images . |
| 3 | Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago |
| 4 | Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann |
| 5 | Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah |
| 6 | |
| 7 | This library is free software; you can redistribute it and/or |
| 8 | modify it under the terms of the GNU Lesser General Public License |
| 9 | (LGPL) as published by the Free Software Foundation; either |
| 10 | version 2.1 of the License, or (at your option) any later version. |
| 11 | The terms of redistributing and/or modifying this software also |
| 12 | include exceptions to the LGPL that facilitate static linking. |
| 13 | |
| 14 | This library is distributed in the hope that it will be useful, |
| 15 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 16 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| 17 | Lesser General Public License for more details. |
| 18 | |
| 19 | You should have received a copy of the GNU Lesser General Public License |
| 20 | along with this library; if not, write to Free Software Foundation, Inc., |
| 21 | 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA |
| 22 | */ |
| 23 | |
| 24 | #include "teem/meet.h" |
| 25 | |
| 26 | /* |
| 27 | ** Tests: |
| 28 | ** ... lots of gage stuff ... |
| 29 | ** |
| 30 | ** The main point of all this is to make sure of two (or maybe 2.5) separate |
| 31 | ** things about gage, the testing of which requires so much in common that one |
| 32 | ** program might as well do them all. The tasks are: |
| 33 | ** |
| 34 | ** 1) ensure that values and their derivatives (where the gageKind supports |
| 35 | ** it) are correctly handled in the multi-value gageKinds (gageKindVec, |
| 36 | ** tenGageKind, tenDwiGageKind), relative to the gageKindScl ground-truth: the |
| 37 | ** per-component values and derivatives have to match those reconstructed from |
| 38 | ** sets of scalar volumes of the components. |
| 39 | ** ==> Also, that gageContextCopy works even on dynamic kinds (like DWIs) |
| 40 | ** |
| 41 | ** 1.5) that there is expected consistency between the scalar, vector, tensor, |
| 42 | ** and DWI properties of a related set of volumes. So the test starts with a |
| 43 | ** diffusion tensor field and generates DWIs, scalars (norm squared of the |
| 44 | ** tensor), and vectors (gradient of norm squared) from this. |
| 45 | ** |
| 46 | ** 2) that scale-space reconstruction works: that sets of pre-blurred volumes |
| 47 | ** can be generated and saved via the utilities in meet, that the smart |
| 48 | ** hermite-spline based scale-space interpolation is working (for all kinds), |
| 49 | ** and that gageContextCopy works on scale-space contexts |
| 50 | */ |
| 51 | |
| 52 | #define BKEY"probeSS" "probeSS" |
| 53 | #define KIND_NUM4 4 |
| 54 | #define KI_SCL0 0 |
| 55 | #define KI_VEC1 1 |
| 56 | #define KI_TEN2 2 |
| 57 | #define KI_DWI3 3 |
| 58 | |
| 59 | /* |
| 60 | static const char *kindStr[KIND_NUM] = { |
| 61 | "scalar", |
| 62 | "vector", |
| 63 | "tensor", |
| 64 | " DWI " |
| 65 | }; |
| 66 | */ |
| 67 | |
| 68 | #define HALTON_BASE100 100 |
| 69 | |
| 70 | static unsigned int volSize[3] = {45, 46, 47}; |
| 71 | |
| 72 | #define NRRD_NEW(nn, mm) \ |
| 73 | (nn) = nrrdNew(); \ |
| 74 | airMopAdd((mm), (nn), (airMopper)nrrdNuke, airMopAlways) |
| 75 | |
| 76 | /* the weird function names of the various local functions here (created in a |
| 77 | rather adhoc and organic way) should be read in usual Teem order: from |
| 78 | right to left */ |
| 79 | |
| 80 | static int |
| 81 | engageGenTensor(gageContext *gctx, Nrrd *ncten, |
| 82 | /* out/in */ |
| 83 | double noiseStdv, unsigned int seed, |
| 84 | unsigned int sx, unsigned int sy, unsigned int sz) { |
| 85 | static const char me[]="engageGenTensor"; |
| 86 | hestParm *hparm; |
| 87 | airArray *smop; |
| 88 | char tmpStr[4][AIR_STRLEN_SMALL(128+1)]; |
| 89 | Nrrd *nclean; |
| 90 | NrrdIter *narg0, *narg1; |
| 91 | const char *helixArgv[] = |
| 92 | /* 0 1 2 3 4 5 6 7 8 9 */ |
| 93 | {"-s", NULL((void*)0), NULL((void*)0), NULL((void*)0), "-v", "0", "-r", "40", "-o", NULL((void*)0), |
| 94 | "-ev", "0.00086", "0.00043", "0.00021", "-bg", "0.003", "-b", "3", |
| 95 | NULL((void*)0)}; |
| 96 | int helixArgc; |
| 97 | gagePerVolume *pvl; |
| 98 | |
| 99 | smop = airMopNew(); |
| 100 | /* NOTE: this is currently the only place where a unrrduCmd |
| 101 | is called from within C code; it was educational to get working. |
| 102 | Learned: |
| 103 | * hest does NOT tolerate having empty or NULL elements of |
| 104 | its argv[]! More error checking for this in hest is needed. |
| 105 | * the "const char **argv" type is not very convenient to |
| 106 | set up in a dynamic way; the per-element setting done below |
| 107 | is certainly awkward |
| 108 | */ |
| 109 | hparm = hestParmNew(); |
| 110 | airMopAdd(smop, hparm, (airMopper)hestParmFree, airMopAlways); |
| 111 | sprintf(tmpStr[0], "%u", sx)__builtin___sprintf_chk (tmpStr[0], 0, __builtin_object_size ( tmpStr[0], 2 > 1 ? 1 : 0), "%u", sx); helixArgv[1] = tmpStr[0]; |
| 112 | sprintf(tmpStr[1], "%u", sy)__builtin___sprintf_chk (tmpStr[1], 0, __builtin_object_size ( tmpStr[1], 2 > 1 ? 1 : 0), "%u", sy); helixArgv[2] = tmpStr[1]; |
| 113 | sprintf(tmpStr[2], "%u", sz)__builtin___sprintf_chk (tmpStr[2], 0, __builtin_object_size ( tmpStr[2], 2 > 1 ? 1 : 0), "%u", sz); helixArgv[3] = tmpStr[2]; |
| 114 | sprintf(tmpStr[3], "tmp-ten.nrrd")__builtin___sprintf_chk (tmpStr[3], 0, __builtin_object_size ( tmpStr[3], 2 > 1 ? 1 : 0), "tmp-ten.nrrd"); |
| 115 | helixArgv[9] = tmpStr[3]; |
| 116 | helixArgc = AIR_CAST(int, sizeof(helixArgv)/sizeof(char *))((int)(sizeof(helixArgv)/sizeof(char *))) - 1; |
| 117 | if (tend_helixCmd.main(helixArgc, helixArgv, me, hparm)) { |
| 118 | /* error already went to stderr, not to any biff container */ |
| 119 | biffAddf(BKEY"probeSS", "%s: problem running tend %s", me, tend_helixCmd.name); |
| 120 | airMopError(smop); return 1; |
| 121 | } |
| 122 | NRRD_NEW(nclean, smop); |
| 123 | if (nrrdLoad(nclean, tmpStr[3], NULL((void*)0))) { |
| 124 | biffAddf(BKEY"probeSS", "%s: trouble loading from new vol %s", me, tmpStr[3]); |
| 125 | airMopError(smop); return 1; |
| 126 | } |
| 127 | |
| 128 | /* add some noise to tensor value; no, this isn't really physical; |
| 129 | since we're adding noise to the tensor and then simulating DWIs, |
| 130 | rather than adding noise to DWIs and then estimating tensor, |
| 131 | but for the purposes of gage testing its fine */ |
| 132 | narg0 = nrrdIterNew(); |
| 133 | narg1 = nrrdIterNew(); |
| 134 | airMopAdd(smop, narg0, (airMopper)nrrdIterNix, airMopAlways); |
| 135 | airMopAdd(smop, narg1, (airMopper)nrrdIterNix, airMopAlways); |
| 136 | nrrdIterSetNrrd(narg0, nclean); |
| 137 | nrrdIterSetValue(narg1, noiseStdv); |
| 138 | airSrandMT(seed); |
| 139 | if (nrrdArithIterBinaryOp(ncten, nrrdBinaryOpNormalRandScaleAdd, |
| 140 | narg0, narg1)) { |
| 141 | biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble noisying output", me); |
| 142 | airMopError(smop); return 1; |
| 143 | } |
| 144 | |
| 145 | /* wrap in gage context */ |
| 146 | if ( !(pvl = gagePerVolumeNew(gctx, ncten, tenGageKind)) |
| 147 | || gagePerVolumeAttach(gctx, pvl) ) { |
| 148 | biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble engaging tensor", me); |
| 149 | return 1; |
| 150 | } |
| 151 | |
| 152 | airMopOkay(smop); |
| 153 | return 0; |
| 154 | } |
| 155 | |
| 156 | /* |
| 157 | ** takes in ncten, measures "S" to nscl, copies that to nsclCopy, |
| 158 | ** wraps those in gctx and gctxComp |
| 159 | */ |
| 160 | static int |
| 161 | engageGenScalar(gageContext *gctx, Nrrd *nscl, |
| 162 | gageContext *gctxComp, Nrrd *nsclCopy, |
| 163 | /* out/in */ |
| 164 | const Nrrd *ncten) { |
| 165 | static const char me[]="engageGenScalar"; |
| 166 | gagePerVolume *pvl; |
| 167 | |
| 168 | if (tenAnisoVolume(nscl, ncten, tenAniso_S, 0)) { |
| 169 | biffMovef(BKEY"probeSS", TENtenBiffKey, "%s: trouble creating scalar volume", me); |
| 170 | return 1; |
| 171 | } |
| 172 | if (nrrdCopy(nsclCopy, nscl)) { |
| 173 | biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble copying scalar volume", me); |
| 174 | return 1; |
| 175 | } |
| 176 | |
| 177 | /* wrap both in gage context */ |
| 178 | if ( !(pvl = gagePerVolumeNew(gctx, nscl, gageKindScl)) |
| 179 | || gagePerVolumeAttach(gctx, pvl) |
| 180 | || !(pvl = gagePerVolumeNew(gctxComp, nsclCopy, gageKindScl)) |
| 181 | || gagePerVolumeAttach(gctxComp, pvl)) { |
| 182 | biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble engaging scalars", me); |
| 183 | return 1; |
| 184 | } |
| 185 | |
| 186 | return 0; |
| 187 | } |
| 188 | |
| 189 | /* |
| 190 | ** Makes a vector volume by measuring the gradient |
| 191 | ** Being the gradient of the given scalar volume is just to make |
| 192 | ** something vaguely interesting, and to test consistency with |
| 193 | ** the the same gradient as measured in the tensor field. |
| 194 | */ |
| 195 | static int |
| 196 | engageGenVector(gageContext *gctx, Nrrd *nvec, |
| 197 | /* out/in */ |
| 198 | const Nrrd *nscl) { |
| 199 | static const char me[]="engageGenVector"; |
| 200 | ptrdiff_t padMin[4] = {0, 0, 0, 0}, padMax[4]; |
| 201 | Nrrd *ntmp; |
| 202 | airArray *smop; |
| 203 | float *vec, *scl; |
| 204 | size_t sx, sy, sz, xi, yi, zi, Px, Mx, Py, My, Pz, Mz; |
| 205 | double dnmX, dnmY, dnmZ; |
| 206 | gagePerVolume *pvl; |
| 207 | |
| 208 | smop = airMopNew(); |
| 209 | if (nrrdTypeFloat != nscl->type) { |
| 210 | biffAddf(BKEY"probeSS", "%s: expected %s not %s type", me, |
| 211 | airEnumStr(nrrdType, nrrdTypeFloat), |
| 212 | airEnumStr(nrrdType, nscl->type)); |
| 213 | airMopError(smop); return 1; |
| 214 | } |
| 215 | NRRD_NEW(ntmp, smop); |
| 216 | sx = nscl->axis[0].size; |
| 217 | sy = nscl->axis[1].size; |
| 218 | sz = nscl->axis[2].size; |
| 219 | ELL_4V_SET(padMax, 2,((padMax)[0] = (2), (padMax)[1] = (((ptrdiff_t)(sx-1))), (padMax )[2] = (((ptrdiff_t)(sy-1))), (padMax)[3] = (((ptrdiff_t)(sz- 1)))) |
| 220 | AIR_CAST(ptrdiff_t, sx-1),((padMax)[0] = (2), (padMax)[1] = (((ptrdiff_t)(sx-1))), (padMax )[2] = (((ptrdiff_t)(sy-1))), (padMax)[3] = (((ptrdiff_t)(sz- 1)))) |
| 221 | AIR_CAST(ptrdiff_t, sy-1),((padMax)[0] = (2), (padMax)[1] = (((ptrdiff_t)(sx-1))), (padMax )[2] = (((ptrdiff_t)(sy-1))), (padMax)[3] = (((ptrdiff_t)(sz- 1)))) |
| 222 | AIR_CAST(ptrdiff_t, sz-1))((padMax)[0] = (2), (padMax)[1] = (((ptrdiff_t)(sx-1))), (padMax )[2] = (((ptrdiff_t)(sy-1))), (padMax)[3] = (((ptrdiff_t)(sz- 1)))); |
| 223 | /* we do axinsert and pad in order to keep all the per-axis info */ |
| 224 | if (nrrdAxesInsert(ntmp, nscl, 0) |
| 225 | || nrrdPad_nva(nvec, ntmp, padMin, padMax, |
| 226 | nrrdBoundaryPad, 0.0)) { |
| 227 | biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble", me); |
| 228 | airMopError(smop); return 1; |
| 229 | } |
| 230 | dnmX = 0.5/nrrdSpaceVecNorm(nscl->spaceDim, nscl->axis[0].spaceDirection); |
| 231 | dnmY = 0.5/nrrdSpaceVecNorm(nscl->spaceDim, nscl->axis[1].spaceDirection); |
| 232 | dnmZ = 0.5/nrrdSpaceVecNorm(nscl->spaceDim, nscl->axis[2].spaceDirection); |
| 233 | vec = AIR_CAST(float *, nvec->data)((float *)(nvec->data)); |
| 234 | scl = AIR_CAST(float *, nscl->data)((float *)(nscl->data)); |
| 235 | #define INDEX(xj, yj, zj) (xj + sx*(yj + sy*zj)) |
| 236 | for (zi=0; zi<sz; zi++) { |
| 237 | Mz = (zi == 0 ? 0 : zi-1); |
| 238 | Pz = AIR_MIN(zi+1, sz-1)((zi+1) < (sz-1) ? (zi+1) : (sz-1)); |
| 239 | for (yi=0; yi<sy; yi++) { |
| 240 | My = (yi == 0 ? 0 : yi-1); |
| 241 | Py = AIR_MIN(yi+1, sy-1)((yi+1) < (sy-1) ? (yi+1) : (sy-1)); |
| 242 | for (xi=0; xi<sx; xi++) { |
| 243 | Mx = (xi == 0 ? 0 : xi-1); |
| 244 | Px = AIR_MIN(xi+1, sx-1)((xi+1) < (sx-1) ? (xi+1) : (sx-1)); |
| 245 | vec[0 + 3*INDEX(xi, yi, zi)] = |
| 246 | AIR_CAST(float, (scl[INDEX(Px, yi, zi)] - scl[INDEX(Mx, yi, zi)])*dnmX)((float)((scl[INDEX(Px, yi, zi)] - scl[INDEX(Mx, yi, zi)])*dnmX )); |
| 247 | vec[1 + 3*INDEX(xi, yi, zi)] = |
| 248 | AIR_CAST(float, (scl[INDEX(xi, Py, zi)] - scl[INDEX(xi, My, zi)])*dnmY)((float)((scl[INDEX(xi, Py, zi)] - scl[INDEX(xi, My, zi)])*dnmY )); |
| 249 | vec[2 + 3*INDEX(xi, yi, zi)] = |
| 250 | AIR_CAST(float, (scl[INDEX(xi, yi, Pz)] - scl[INDEX(xi, yi, Mz)])*dnmZ)((float)((scl[INDEX(xi, yi, Pz)] - scl[INDEX(xi, yi, Mz)])*dnmZ )); |
| 251 | } |
| 252 | } |
| 253 | } |
| 254 | #undef INDEX |
| 255 | |
| 256 | /* wrap in gage context */ |
| 257 | if ( !(pvl = gagePerVolumeNew(gctx, nvec, gageKindVec)) |
| 258 | || gagePerVolumeAttach(gctx, pvl) ) { |
| 259 | biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble engaging vectors", me); |
| 260 | return 1; |
| 261 | } |
| 262 | airMopError(smop); |
| 263 | return 0; |
| 264 | } |
| 265 | |
| 266 | /* |
| 267 | ** make a DWI volume by simulating DWIs from given tensor |
| 268 | ** this includes generating a new gradient set |
| 269 | */ |
| 270 | static int |
| 271 | genDwi(Nrrd *ndwi, Nrrd *ngrad, |
| 272 | /* out/in */ |
| 273 | unsigned int gradNum, double bval, const Nrrd *ncten) { |
| 274 | static const char me[]="genDwi"; |
| 275 | tenGradientParm *gparm; |
| 276 | tenExperSpec *espec; |
| 277 | Nrrd *nten, *nb0; |
| 278 | NrrdIter *narg0, *narg1; |
| 279 | size_t cropMin[4] = {1, 0, 0, 0}, cropMax[4]; |
| 280 | airArray *smop; |
| 281 | |
| 282 | smop = airMopNew(); |
| 283 | gparm = tenGradientParmNew(); |
| 284 | airMopAdd(smop, gparm, (airMopper)tenGradientParmNix, airMopAlways); |
| 285 | espec = tenExperSpecNew(); |
| 286 | airMopAdd(smop, espec, (airMopper)tenExperSpecNix, airMopAlways); |
| 287 | gparm->verbose = 0; |
| 288 | gparm->minMean = 0.002; |
| 289 | gparm->seed = 4242; |
| 290 | gparm->insertZeroVec = AIR_TRUE1; |
| 291 | if (tenGradientGenerate(ngrad, gradNum, gparm) |
| 292 | || tenExperSpecGradSingleBValSet(espec, AIR_FALSE0 /* insertB0 */, |
| 293 | bval, |
| 294 | AIR_CAST(double *, ngrad->data)((double *)(ngrad->data)), |
| 295 | AIR_CAST(unsigned int,((unsigned int)(ngrad->axis[1].size)) |
| 296 | ngrad->axis[1].size)((unsigned int)(ngrad->axis[1].size)))) { |
| 297 | biffMovef(BKEY"probeSS", TENtenBiffKey, "%s: trouble generating grads or espec", me); |
| 298 | airMopError(smop); return 1; |
| 299 | } |
| 300 | NRRD_NEW(nten, smop); |
| 301 | NRRD_NEW(nb0, smop); |
| 302 | ELL_4V_SET(cropMax,((cropMax)[0] = (ncten->axis[0].size-1), (cropMax)[1] = (ncten ->axis[1].size-1), (cropMax)[2] = (ncten->axis[2].size- 1), (cropMax)[3] = (ncten->axis[3].size-1)) |
| 303 | ncten->axis[0].size-1,((cropMax)[0] = (ncten->axis[0].size-1), (cropMax)[1] = (ncten ->axis[1].size-1), (cropMax)[2] = (ncten->axis[2].size- 1), (cropMax)[3] = (ncten->axis[3].size-1)) |
| 304 | ncten->axis[1].size-1,((cropMax)[0] = (ncten->axis[0].size-1), (cropMax)[1] = (ncten ->axis[1].size-1), (cropMax)[2] = (ncten->axis[2].size- 1), (cropMax)[3] = (ncten->axis[3].size-1)) |
| 305 | ncten->axis[2].size-1,((cropMax)[0] = (ncten->axis[0].size-1), (cropMax)[1] = (ncten ->axis[1].size-1), (cropMax)[2] = (ncten->axis[2].size- 1), (cropMax)[3] = (ncten->axis[3].size-1)) |
| 306 | ncten->axis[3].size-1)((cropMax)[0] = (ncten->axis[0].size-1), (cropMax)[1] = (ncten ->axis[1].size-1), (cropMax)[2] = (ncten->axis[2].size- 1), (cropMax)[3] = (ncten->axis[3].size-1)); |
| 307 | if (nrrdSlice(nb0, ncten, 0, 0) |
| 308 | || nrrdCrop(nten, ncten, cropMin, cropMax)) { |
| 309 | biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble slicing or cropping ten vol", me); |
| 310 | airMopError(smop); return 1; |
| 311 | } |
| 312 | narg0 = nrrdIterNew(); |
| 313 | narg1 = nrrdIterNew(); |
| 314 | airMopAdd(smop, narg0, (airMopper)nrrdIterNix, airMopAlways); |
| 315 | airMopAdd(smop, narg1, (airMopper)nrrdIterNix, airMopAlways); |
| 316 | nrrdIterSetValue(narg1, 50000.0); |
| 317 | nrrdIterSetNrrd(narg0, nb0); |
| 318 | if (nrrdArithIterBinaryOp(nb0, nrrdBinaryOpMultiply, |
| 319 | narg0, narg1)) { |
| 320 | biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble generating b0 vol", me); |
| 321 | airMopError(smop); return 1; |
| 322 | } |
| 323 | if (tenModelSimulate(ndwi, nrrdTypeUShort, espec, |
| 324 | tenModel1Tensor2, |
| 325 | nb0, nten, AIR_TRUE1 /* keyValueSet */)) { |
| 326 | biffMovef(BKEY"probeSS", TENtenBiffKey, "%s: trouble simulating DWI vol", me); |
| 327 | airMopError(smop); return 1; |
| 328 | } |
| 329 | |
| 330 | airMopOkay(smop); |
| 331 | return 0; |
| 332 | } |
| 333 | |
| 334 | int |
| 335 | engageMopDiceVector(gageContext *gctx, Nrrd *nvecComp[3], |
| 336 | /* out/in */ |
| 337 | airArray *mop, const Nrrd* nvec) { |
| 338 | static const char me[]="engageMopDiceVector"; |
| 339 | gagePerVolume *pvl; |
| 340 | unsigned int ci; |
| 341 | char stmp[AIR_STRLEN_SMALL(128+1)]; |
| 342 | |
| 343 | if (!( 4 == nvec->dim && 3 == nvec->axis[0].size )) { |
| 344 | biffAddf(BKEY"probeSS", "%s: expected 4-D 3-by-X nrrd (not %u-D %s-by-X)", |
| 345 | me, nvec->dim, airSprintSize_t(stmp, nvec->axis[0].size)); |
| 346 | return 1; |
| 347 | } |
| 348 | for (ci=0; ci<3; ci++) { |
| 349 | NRRD_NEW(nvecComp[ci], mop); |
| 350 | if (nrrdSlice(nvecComp[ci], nvec, 0, ci)) { |
| 351 | biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble getting component %u", me, ci); |
| 352 | return 1; |
| 353 | } |
| 354 | if ( !(pvl = gagePerVolumeNew(gctx, nvecComp[ci], gageKindScl)) |
| 355 | || gagePerVolumeAttach(gctx, pvl) ) { |
| 356 | biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble engaging component %u", me, ci); |
| 357 | return 1; |
| 358 | } |
| 359 | } |
| 360 | |
| 361 | return 0; |
| 362 | } |
| 363 | |
| 364 | int |
| 365 | engageMopDiceTensor(gageContext *gctx, Nrrd *nctenComp[7], |
| 366 | /* out/in */ |
| 367 | airArray *mop, const Nrrd* ncten) { |
| 368 | static const char me[]="engageMopDiceTensor"; |
| 369 | gagePerVolume *pvl; |
| 370 | unsigned int ci; |
| 371 | |
| 372 | if (tenTensorCheck(ncten, nrrdTypeFloat, AIR_TRUE1 /* want4F */, |
| 373 | AIR_TRUE1 /* useBiff */)) { |
| 374 | biffMovef(BKEY"probeSS", TENtenBiffKey, "%s: didn't get tensor volume", me); |
| 375 | return 1; |
| 376 | } |
| 377 | for (ci=0; ci<7; ci++) { |
| 378 | NRRD_NEW(nctenComp[ci], mop); |
| 379 | if (nrrdSlice(nctenComp[ci], ncten, 0, ci)) { |
| 380 | biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble getting component %u", me, ci); |
| 381 | return 1; |
| 382 | } |
| 383 | if ( !(pvl = gagePerVolumeNew(gctx, nctenComp[ci], gageKindScl)) |
| 384 | || gagePerVolumeAttach(gctx, pvl) ) { |
| 385 | biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble engaging component %u", me, ci); |
| 386 | return 1; |
| 387 | } |
| 388 | } |
| 389 | |
| 390 | return 0; |
| 391 | } |
| 392 | |
| 393 | int |
| 394 | engageMopDiceDwi(gageContext *gctx, Nrrd ***ndwiCompP, |
| 395 | /* out/in */ |
| 396 | airArray *mop, const Nrrd* ndwi) { |
| 397 | static const char me[]="mopDiceDwi"; |
| 398 | Nrrd **ndwiComp; |
| 399 | size_t dwiNum; |
| 400 | char stmp[AIR_STRLEN_SMALL(128+1)]; |
| 401 | gagePerVolume *pvl; |
| 402 | unsigned int ci; |
| 403 | |
| 404 | if (!( 4 == ndwi->dim )) { |
| 405 | biffAddf(BKEY"probeSS", "%s: wanted 4D volume (not %u)", me, ndwi->dim); |
| 406 | return 1; |
| 407 | } |
| 408 | if (!( nrrdKindList == ndwi->axis[0].kind && |
| 409 | nrrdKindSpace == ndwi->axis[1].kind && |
| 410 | nrrdKindSpace == ndwi->axis[2].kind && |
| 411 | nrrdKindSpace == ndwi->axis[3].kind )) { |
| 412 | biffAddf(BKEY"probeSS", "%s: wanted kinds %s,3x%s, not %s,%s,%s,%s", me, |
| 413 | airEnumStr(nrrdKind, nrrdKindList), |
| 414 | airEnumStr(nrrdKind, nrrdKindSpace), |
| 415 | airEnumStr(nrrdKind, ndwi->axis[0].kind), |
| 416 | airEnumStr(nrrdKind, ndwi->axis[1].kind), |
| 417 | airEnumStr(nrrdKind, ndwi->axis[2].kind), |
| 418 | airEnumStr(nrrdKind, ndwi->axis[3].kind)); |
| 419 | return 1; |
| 420 | } |
| 421 | dwiNum = ndwi->axis[0].size; |
| 422 | if (!(ndwiComp = AIR_CALLOC(dwiNum, Nrrd *)(Nrrd **)(calloc((dwiNum), sizeof(Nrrd *))))) { |
| 423 | biffAddf(BKEY"probeSS", "%s: couldn't alloc %s Nrrd*", me, |
| 424 | airSprintSize_t(stmp, dwiNum)); |
| 425 | return 1; |
| 426 | } |
| 427 | airMopAdd(mop, ndwiComp, airFree, airMopAlways); |
| 428 | *ndwiCompP = ndwiComp; |
| 429 | for (ci=0; ci<dwiNum; ci++) { |
| 430 | NRRD_NEW(ndwiComp[ci], mop); |
| 431 | if (nrrdSlice(ndwiComp[ci], ndwi, 0, ci)) { |
| 432 | biffMovef(BKEY"probeSS", NRRDnrrdBiffKey, "%s: trouble getting component %u", me, ci); |
| 433 | return 1; |
| 434 | } |
| 435 | if ( !(pvl = gagePerVolumeNew(gctx, ndwiComp[ci], gageKindScl)) |
| 436 | || gagePerVolumeAttach(gctx, pvl) ) { |
| 437 | biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble engaging component %u", me, ci); |
| 438 | return 1; |
| 439 | } |
| 440 | } |
| 441 | return 0; |
| 442 | } |
| 443 | |
| 444 | typedef struct { |
| 445 | char name[AIR_STRLEN_SMALL(128+1)]; |
| 446 | const double **aptr; /* array of pointers to (const) answers */ |
| 447 | gageItemSpec *ispec; /* array of gageItemSpecs (not pointers to them) */ |
| 448 | unsigned int *alen; /* array of answer lengths */ |
| 449 | unsigned int *sidx; /* array of index (within san) of copied answer */ |
| 450 | unsigned int anum; /* lenth of aptr, ispec, alen, sidx arrays */ |
| 451 | double *san; /* single buffer for copy + concating answers */ |
| 452 | unsigned int slen; /* length of single buffer (sum of all alen[i]) */ |
| 453 | } multiAnswer; |
| 454 | |
| 455 | static void |
| 456 | multiAnswerInit(multiAnswer *man) { |
| 457 | |
| 458 | /* HEY sloppy: not resetting name */ |
| 459 | man->aptr = airFree(AIR_VOIDP(man->aptr)((void *)(man->aptr))); |
| 460 | man->ispec = airFree(man->ispec); |
| 461 | man->alen = airFree(man->alen); |
| 462 | man->sidx = airFree(man->sidx); |
| 463 | man->anum = 0; |
| 464 | man->san = airFree(man->san); |
| 465 | man->slen = 0; |
| 466 | } |
| 467 | |
| 468 | static multiAnswer* |
| 469 | multiAnswerNew(char *name) { |
| 470 | multiAnswer *man; |
| 471 | |
| 472 | man = AIR_CALLOC(1, multiAnswer)(multiAnswer*)(calloc((1), sizeof(multiAnswer))); |
| 473 | airStrcpy(man->name, AIR_STRLEN_SMALL(128+1), name); |
| 474 | multiAnswerInit(man); |
| 475 | return man; |
| 476 | } |
| 477 | |
| 478 | void |
| 479 | multiAnswerLenSet(multiAnswer *man, unsigned int anum) { |
| 480 | /* |
| 481 | static const char me[]="multiAnswerLenSet"; |
| 482 | |
| 483 | fprintf(stderr, "!%s: %s hello -> answer number = %u\n", me, |
| 484 | man->name, anum); |
| 485 | */ |
| 486 | man->aptr = AIR_CALLOC(anum, const double *)(const double **)(calloc((anum), sizeof(const double *))); |
| 487 | man->ispec = AIR_CALLOC(anum, gageItemSpec)(gageItemSpec*)(calloc((anum), sizeof(gageItemSpec))); |
| 488 | man->alen = AIR_CALLOC(anum, unsigned int)(unsigned int*)(calloc((anum), sizeof(unsigned int))); |
| 489 | man->sidx = AIR_CALLOC(anum, unsigned int)(unsigned int*)(calloc((anum), sizeof(unsigned int))); |
| 490 | man->anum = anum; |
| 491 | /* don't know answer lengths yet */ |
| 492 | man->san = airFree(man->san); |
| 493 | man->slen = 0; |
| 494 | return; |
| 495 | } |
| 496 | |
| 497 | static multiAnswer* |
| 498 | multiAnswerNix(multiAnswer *man) { |
| 499 | |
| 500 | airFree(AIR_VOIDP(man->aptr)((void *)(man->aptr))); |
| 501 | airFree(man->ispec); |
| 502 | airFree(man->alen); |
| 503 | airFree(man->sidx); |
| 504 | airFree(man->san); |
| 505 | airFree(man); |
| 506 | return NULL((void*)0); |
| 507 | } |
| 508 | |
| 509 | static void |
| 510 | multiAnswerAdd(multiAnswer *man, unsigned int ansIdx, |
| 511 | const gageContext *gctx, const gagePerVolume *pvl, |
| 512 | unsigned int item) { |
| 513 | /* |
| 514 | static const char me[]="multiAnswerAdd"; |
| 515 | |
| 516 | fprintf(stderr, "!%s: %s hello (aidx = %u)\n", me, man->name, ansIdx); |
| 517 | */ |
| 518 | man->aptr[ansIdx] = gageAnswerPointer(gctx, pvl, item); |
| 519 | man->ispec[ansIdx].kind = pvl->kind; |
| 520 | man->ispec[ansIdx].item = item; |
| 521 | man->alen[ansIdx] = gageAnswerLength(gctx, pvl, item); |
| 522 | /* HEY hack: doing this tally and allocation only on the last time |
| 523 | we're called, presuming calls with increasing ansIdx */ |
| 524 | man->slen = 0; |
| 525 | if (ansIdx == man->anum-1) { |
| 526 | unsigned int ai; |
| 527 | /* fprintf(stderr, "!%s: %s hello aidx reached anum-1 = %u\n", me, man->name, man->anum-1); */ |
| 528 | for (ai=0; ai<man->anum; ai++) { |
| 529 | man->sidx[ai] = man->slen; |
| 530 | man->slen += man->alen[ai]; |
| 531 | } |
| 532 | man->san = AIR_CALLOC(man->slen, double)(double*)(calloc((man->slen), sizeof(double))); |
| 533 | /* |
| 534 | fprintf(stderr, "!%s: (%s) slen = %u\n", "multiAnswerAdd", |
| 535 | pvl->kind->name, man->slen); |
| 536 | */ |
| 537 | } |
| 538 | return; |
| 539 | } |
| 540 | |
| 541 | static void |
| 542 | multiAnswerCollect(multiAnswer *man) { |
| 543 | unsigned int ai; |
| 544 | /* |
| 545 | static const char me[]="multiAnswerCollect"; |
| 546 | |
| 547 | fprintf(stderr, "%s: %s hi\n", me, man->name); |
| 548 | */ |
| 549 | for (ai=0; ai<man->anum; ai++) { |
| 550 | /* |
| 551 | fprintf(stderr, "!%s: (%s/%s) ai=%u/%u to %p+%u=%p for %u doubles\n", "multiAnswerCollect", |
| 552 | man->ispec[ai].kind->name, |
| 553 | airEnumStr(man->ispec[ai].kind->enm, man->ispec[ai].item), |
| 554 | ai, man->anum, man->san, man->sidx[ai], man->san + man->sidx[ai], |
| 555 | man->alen[ai]); |
| 556 | */ |
| 557 | memcpy(man->san + man->sidx[ai],__builtin___memcpy_chk (man->san + man->sidx[ai], man-> aptr[ai], man->alen[ai]*sizeof(double), __builtin_object_size (man->san + man->sidx[ai], 0)) |
| 558 | man->aptr[ai],__builtin___memcpy_chk (man->san + man->sidx[ai], man-> aptr[ai], man->alen[ai]*sizeof(double), __builtin_object_size (man->san + man->sidx[ai], 0)) |
| 559 | man->alen[ai]*sizeof(double))__builtin___memcpy_chk (man->san + man->sidx[ai], man-> aptr[ai], man->alen[ai]*sizeof(double), __builtin_object_size (man->san + man->sidx[ai], 0)); |
| 560 | } |
| 561 | return; |
| 562 | } |
| 563 | |
| 564 | static int |
| 565 | multiAnswerCompare(multiAnswer *man1, multiAnswer *man2) { |
| 566 | static const char me[]="multiAnswerCompare"; |
| 567 | unsigned int si, slen; |
| 568 | |
| 569 | |
| 570 | #if 1 |
| 571 | if (man1->slen != man2->slen) { |
| 572 | biffAddf(BKEY"probeSS", "%s: man1->slen %u != man2->slen %u\n", me, |
| 573 | man1->slen, man2->slen); |
| 574 | return 1; |
| 575 | } |
| 576 | slen = man1->slen; |
| 577 | #else |
| 578 | slen = AIR_MIN(man1->slen, man2->slen)((man1->slen) < (man2->slen) ? (man1->slen) : (man2 ->slen)); |
| 579 | #endif |
| 580 | for (si=0; si<slen; si++) { |
| 581 | if (man1->san[si] != man2->san[si]) { |
| 582 | /* HEY should track down which part of which answer, |
| 583 | in man1 and man2, is different, which was the |
| 584 | purpose of recording ispec and sidx */ |
| 585 | biffAddf(BKEY"probeSS", "%s: man1->san[si] %.17g != man2->san[si] %.17g", me, |
| 586 | man1->san[si], man2->san[si]); |
| 587 | return 1; |
| 588 | } |
| 589 | } |
| 590 | return 0; |
| 591 | } |
| 592 | |
| 593 | /* |
| 594 | ** setting up gageContexts for the first of the two tasks listed above: making |
| 595 | ** sure per-component information is handled correctly. NOTE: The combination |
| 596 | ** of function calls made here is very atypical for a Teem-using program |
| 597 | */ |
| 598 | static int |
| 599 | updateQueryKernelSetTask1(gageContext *gctxComp[KIND_NUM4], |
| 600 | gageContext *gctx[KIND_NUM4], int gSetup, |
| 601 | multiAnswer *manComp[KIND_NUM4], |
| 602 | multiAnswer *man[KIND_NUM4], |
| 603 | NrrdKernel *kpack[3], |
| 604 | double support) { |
| 605 | static const char me[]="updateQueryKernelSetTask1"; |
| 606 | double parm1[NRRD_KERNEL_PARMS_NUM8], parmV[NRRD_KERNEL_PARMS_NUM8]; |
| 607 | unsigned int kindIdx; |
| 608 | |
| 609 | if (4 != KIND_NUM4) { |
| 610 | biffAddf(BKEY"probeSS", "%s: sorry, confused: KIND_NUM %u != 4", |
| 611 | me, KIND_NUM4); |
| 612 | return 1; |
| 613 | } |
| 614 | parm1[0] = 1.0; |
| 615 | parmV[0] = support; |
| 616 | if (gSetup) { |
| 617 | for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) { |
| 618 | if (0 /* (no-op for formatting) */ |
| 619 | || gageKernelSet(gctxComp[kindIdx], gageKernel00, kpack[0], parm1) |
| 620 | || gageKernelSet(gctxComp[kindIdx], gageKernel11, kpack[1], parm1) |
| 621 | || gageKernelSet(gctxComp[kindIdx], gageKernel22, kpack[2], parm1) |
| 622 | || gageKernelSet(gctx[kindIdx], gageKernel00, kpack[0], parmV) |
| 623 | || gageKernelSet(gctx[kindIdx], gageKernel11, kpack[1], parmV) |
| 624 | || gageKernelSet(gctx[kindIdx], gageKernel22, kpack[2], parmV)) { |
| 625 | biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble setting kernel (kind %u)", |
| 626 | me, kindIdx); |
| 627 | return 1; |
| 628 | } |
| 629 | } |
| 630 | /* Note that the contexts for the diced-up volumes of components are always |
| 631 | of kind gageKindScl, and all the items are from the gageScl* enum */ |
| 632 | for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) { |
| 633 | unsigned int pvi; |
| 634 | /* |
| 635 | fprintf(stderr, "!%s: gctx[%u] has %u pvls\n", me, kindIdx, gctx[kindIdx]->pvlNum); |
| 636 | fprintf(stderr, "!%s: gctxComp[%u] has %u pvls\n", me, kindIdx, gctxComp[kindIdx]->pvlNum); |
| 637 | */ |
| 638 | for (pvi=0; pvi<gctxComp[kindIdx]->pvlNum; pvi++) { |
| 639 | if (0 /* (no-op for formatting) */ |
| 640 | || gageQueryItemOn(gctxComp[kindIdx], gctxComp[kindIdx]->pvl[pvi], gageSclValue) |
| 641 | || gageQueryItemOn(gctxComp[kindIdx], gctxComp[kindIdx]->pvl[pvi], gageSclGradVec) |
| 642 | || gageQueryItemOn(gctxComp[kindIdx], gctxComp[kindIdx]->pvl[pvi], gageSclHessian)) { |
| 643 | biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble setting query (kind %u) on pvi %u", |
| 644 | me, kindIdx, pvi); |
| 645 | return 1; |
| 646 | } |
| 647 | } |
| 648 | if (gageUpdate(gctxComp[kindIdx])) { |
| 649 | biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble updating comp gctx %u", |
| 650 | me, kindIdx); |
| 651 | return 1; |
| 652 | } |
| 653 | } |
| 654 | /* For the original contexts, we have to use the kind-specific items that |
| 655 | correspond to the values and derivatives */ |
| 656 | if (0 /* (no-op for formatting) */ |
| 657 | || gageQueryItemOn(gctx[KI_SCL0], gctx[KI_SCL0]->pvl[0], gageSclValue) |
| 658 | || gageQueryItemOn(gctx[KI_SCL0], gctx[KI_SCL0]->pvl[0], gageSclGradVec) |
| 659 | || gageQueryItemOn(gctx[KI_SCL0], gctx[KI_SCL0]->pvl[0], gageSclHessian) |
| 660 | |
| 661 | || gageQueryItemOn(gctx[KI_VEC1], gctx[KI_VEC1]->pvl[0], gageVecVector) |
| 662 | || gageQueryItemOn(gctx[KI_VEC1], gctx[KI_VEC1]->pvl[0], gageVecJacobian) |
| 663 | || gageQueryItemOn(gctx[KI_VEC1], gctx[KI_VEC1]->pvl[0], gageVecHessian) |
| 664 | |
| 665 | || gageQueryItemOn(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageTensor) |
| 666 | || gageQueryItemOn(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageTensorGrad) |
| 667 | || gageQueryItemOn(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageS) |
| 668 | || gageQueryItemOn(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageSGradVec) |
| 669 | || gageQueryItemOn(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageHessian) |
| 670 | || gageQueryItemOn(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageSHessian) |
| 671 | |
| 672 | || gageQueryItemOn(gctx[KI_DWI3], gctx[KI_DWI3]->pvl[0], tenDwiGageAll) |
| 673 | || gageQueryItemOn(gctx[KI_DWI3], gctx[KI_DWI3]->pvl[0], tenDwiGageTensorLLS)) { |
| 674 | biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble setting item", me); |
| 675 | return 1; |
| 676 | } |
| 677 | for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) { |
| 678 | if (gageUpdate(gctx[kindIdx])) { |
| 679 | biffMovef(BKEY"probeSS", GAGEgageBiffKey, "%s: trouble updating single gctx %u", |
| 680 | me, kindIdx); |
| 681 | return 1; |
| 682 | } |
| 683 | } |
| 684 | } /* if (gSetup) */ |
| 685 | |
| 686 | for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) { |
| 687 | multiAnswerInit(man[kindIdx]); |
| 688 | multiAnswerInit(manComp[kindIdx]); |
| 689 | } |
| 690 | for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) { |
| 691 | unsigned int pvi; |
| 692 | multiAnswerLenSet(manComp[kindIdx], (KI_DWI3 != kindIdx ? 3 : 1)*gctxComp[kindIdx]->pvlNum); |
| 693 | for (pvi=0; pvi<gctxComp[kindIdx]->pvlNum; pvi++) { |
| 694 | multiAnswerAdd(manComp[kindIdx], pvi, |
| 695 | gctxComp[kindIdx], gctxComp[kindIdx]->pvl[pvi], |
| 696 | gageSclValue); |
| 697 | } |
| 698 | if (KI_DWI3 != kindIdx) { |
| 699 | for (pvi=0; pvi<gctxComp[kindIdx]->pvlNum; pvi++) { |
| 700 | multiAnswerAdd(manComp[kindIdx], pvi + gctxComp[kindIdx]->pvlNum, |
| 701 | gctxComp[kindIdx], gctxComp[kindIdx]->pvl[pvi], |
| 702 | gageSclGradVec); |
| 703 | } |
| 704 | for (pvi=0; pvi<gctxComp[kindIdx]->pvlNum; pvi++) { |
| 705 | multiAnswerAdd(manComp[kindIdx], pvi + 2*gctxComp[kindIdx]->pvlNum, |
| 706 | gctxComp[kindIdx], gctxComp[kindIdx]->pvl[pvi], |
| 707 | gageSclHessian); |
| 708 | } |
| 709 | } |
| 710 | multiAnswerLenSet(man[kindIdx], KI_DWI3 != kindIdx ? 3 : 1); |
| 711 | switch(kindIdx) { |
| 712 | case KI_SCL0: |
| 713 | multiAnswerAdd(man[kindIdx], 0, gctx[kindIdx], gctx[kindIdx]->pvl[0], gageSclValue); |
| 714 | multiAnswerAdd(man[kindIdx], 1, gctx[kindIdx], gctx[kindIdx]->pvl[0], gageSclGradVec); |
| 715 | multiAnswerAdd(man[kindIdx], 2, gctx[kindIdx], gctx[kindIdx]->pvl[0], gageSclHessian); |
| 716 | break; |
| 717 | case KI_VEC1: |
| 718 | multiAnswerAdd(man[kindIdx], 0, gctx[kindIdx], gctx[kindIdx]->pvl[0], gageVecVector); |
| 719 | multiAnswerAdd(man[kindIdx], 1, gctx[kindIdx], gctx[kindIdx]->pvl[0], gageVecJacobian); |
| 720 | multiAnswerAdd(man[kindIdx], 2, gctx[kindIdx], gctx[kindIdx]->pvl[0], gageVecHessian); |
| 721 | break; |
| 722 | case KI_TEN2: |
| 723 | multiAnswerAdd(man[kindIdx], 0, gctx[kindIdx], gctx[kindIdx]->pvl[0], tenGageTensor); |
| 724 | multiAnswerAdd(man[kindIdx], 1, gctx[kindIdx], gctx[kindIdx]->pvl[0], tenGageTensorGrad); |
| 725 | multiAnswerAdd(man[kindIdx], 2, gctx[kindIdx], gctx[kindIdx]->pvl[0], tenGageHessian); |
| 726 | break; |
| 727 | case KI_DWI3: |
| 728 | multiAnswerAdd(man[kindIdx], 0, gctx[kindIdx], gctx[kindIdx]->pvl[0], tenDwiGageAll); |
| 729 | break; |
| 730 | } |
| 731 | } |
| 732 | |
| 733 | return 0; |
| 734 | } |
| 735 | |
| 736 | static int |
| 737 | probeTask1(gageContext *gctxComp[KIND_NUM4], |
| 738 | gageContext *gctx[KIND_NUM4], |
| 739 | multiAnswer *manComp[KIND_NUM4], |
| 740 | multiAnswer *man[KIND_NUM4], |
| 741 | unsigned int probeNum, NrrdKernel *kpack[3], |
| 742 | double ksupport) { |
| 743 | static const char me[]="probeTask1"; |
| 744 | unsigned int kindIdx, probeIdx, errNumMax, tenErrNum, sclErrNum; |
| 745 | double pos[3], upos[3], minp[3], maxp[3], |
| 746 | tenDiff[7], tenAvg[7], tenErr, tenErrMax, |
| 747 | vecDiff[3], vecAvg[3], vecErr, vecErrMax, vecErrNum, |
| 748 | sclDiff, sclAvg, sclErr, sclErrMax, errNumFrac; |
| 749 | const double *dwiTenEstP, |
| 750 | *tenTenP, *tenTenNormP, *tenTenNormGradP, |
| 751 | *sclSclP, *sclGradP, *vecVecP; |
| 752 | |
| 753 | ELL_3V_SET(minp, ksupport, ksupport, ksupport)((minp)[0] = (ksupport), (minp)[1] = (ksupport), (minp)[2] = ( ksupport)); |
| 754 | ELL_3V_SET(maxp, volSize[0]-1-ksupport, volSize[1]-1-ksupport, volSize[2]-1-ksupport)((maxp)[0] = (volSize[0]-1-ksupport), (maxp)[1] = (volSize[1] -1-ksupport), (maxp)[2] = (volSize[2]-1-ksupport)); |
| 755 | |
| 756 | /* this is all for task 1.5 */ |
| 757 | sclSclP = gageAnswerPointer(gctx[KI_SCL0], gctx[KI_SCL0]->pvl[0], gageSclValue); |
| 758 | sclGradP = gageAnswerPointer(gctx[KI_SCL0], gctx[KI_SCL0]->pvl[0], gageSclGradVec); |
| 759 | vecVecP = gageAnswerPointer(gctx[KI_VEC1], gctx[KI_VEC1]->pvl[0], gageVecVector); |
| 760 | tenTenP = gageAnswerPointer(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageTensor); |
| 761 | tenTenNormP = gageAnswerPointer(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageS); |
| 762 | tenTenNormGradP = gageAnswerPointer(gctx[KI_TEN2], gctx[KI_TEN2]->pvl[0], tenGageSGradVec); |
| 763 | dwiTenEstP = gageAnswerPointer(gctx[KI_DWI3], gctx[KI_DWI3]->pvl[0], tenDwiGageTensorLLS); |
| 764 | tenErrNum = sclErrNum = 0; |
| 765 | vecErrNum = 0.0; |
| 766 | errNumFrac = 0.02; |
| 767 | |
| 768 | if (nrrdKernelBoxSupportDebug == kpack[0]) { |
| 769 | /* not actually here for any derivatives, mainly to check on |
| 770 | tensor esimation (in addition to usual check of multivariate as |
| 771 | set of scalars) */ |
| 772 | tenErrMax = 0.0002; |
| 773 | vecErrMax = AIR_NAN(airFloatQNaN.f); |
| 774 | sclErrMax = AIR_NAN(airFloatQNaN.f); |
| 775 | } else if (nrrdKernelCos4SupportDebug == kpack[0]) { |
| 776 | tenErrMax = AIR_NAN(airFloatQNaN.f); |
| 777 | vecErrMax = AIR_NAN(airFloatQNaN.f); |
| 778 | sclErrMax = AIR_NAN(airFloatQNaN.f); |
| 779 | } else if (nrrdKernelCatmullRomSupportDebug == kpack[0]) { |
| 780 | tenErrMax = 0.14; /* honestly, not sure how meaningful this test |
| 781 | is, given how significant we're allowing the |
| 782 | error to be. might be more meaningful if |
| 783 | this was with comparison to a pre-computed |
| 784 | non-linear least-squares fit, instead of the |
| 785 | (log-based) linear least-squares. */ |
| 786 | vecErrMax = 0.0016; |
| 787 | sclErrMax = 0.0008; |
| 788 | } else { |
| 789 | biffAddf(BKEY"probeSS", "%s: unexpected kpack[0] %s\n", me, kpack[0]->name); |
| 790 | return 1; |
| 791 | } |
| 792 | errNumMax = AIR_CAST(unsigned int, errNumFrac*probeNum)((unsigned int)(errNumFrac*probeNum)); |
| 793 | for (probeIdx=0; probeIdx<probeNum; probeIdx++) { |
| 794 | airHalton(upos, probeIdx+HALTON_BASE100, airPrimeList + 0, 3); |
| 795 | pos[0] = AIR_AFFINE(0.0, upos[0], 1.0, minp[0], maxp[0])( ((double)(maxp[0])-(minp[0]))*((double)(upos[0])-(0.0)) / ( (double)(1.0)-(0.0)) + (minp[0])); |
| 796 | pos[1] = AIR_AFFINE(0.0, upos[1], 1.0, minp[1], maxp[1])( ((double)(maxp[1])-(minp[1]))*((double)(upos[1])-(0.0)) / ( (double)(1.0)-(0.0)) + (minp[1])); |
| 797 | pos[2] = AIR_AFFINE(0.0, upos[2], 1.0, minp[2], maxp[2])( ((double)(maxp[2])-(minp[2]))*((double)(upos[2])-(0.0)) / ( (double)(1.0)-(0.0)) + (minp[2])); |
| 798 | for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) { |
| 799 | #define PROBE(gctx, what, pos) \ |
| 800 | if (gageProbeSpace((gctx), (pos)[0], (pos)[1], (pos)[2], \ |
| 801 | AIR_TRUE1 /* indexSpace */, \ |
| 802 | AIR_FALSE0 /* clamp */)) { \ |
| 803 | biffAddf(BKEY"probeSS", "%s: probe (%s) error(%d): %s", me, \ |
| 804 | what, (gctx)->errNum, (gctx)->errStr); \ |
| 805 | return 1; \ |
| 806 | } |
| 807 | |
| 808 | PROBE(gctxComp[kindIdx], "comp", pos); |
| 809 | PROBE(gctx[kindIdx], "single", pos); |
| 810 | |
| 811 | #undef PROBE |
| 812 | |
| 813 | multiAnswerCollect(man[kindIdx]); |
| 814 | multiAnswerCollect(manComp[kindIdx]); |
| 815 | if (multiAnswerCompare(manComp[kindIdx], man[kindIdx])) { |
| 816 | biffAddf(BKEY"probeSS", "%s: on point %u of kindIdx %u", me, |
| 817 | probeIdx, kindIdx); |
| 818 | return 1; |
| 819 | } |
| 820 | } |
| 821 | |
| 822 | /* see if probed tensors mostly agree with |
| 823 | tensors estimated from probed DWIs */ |
| 824 | if (AIR_EXISTS(tenErrMax)(((int)(!((tenErrMax) - (tenErrMax)))))) { |
| 825 | TEN_T_SUB(tenDiff, dwiTenEstP, tenTenP)( (tenDiff)[0] = ((dwiTenEstP)[0] + (tenTenP)[0])/2.0, (tenDiff )[1] = (dwiTenEstP)[1] - (tenTenP)[1], (tenDiff)[2] = (dwiTenEstP )[2] - (tenTenP)[2], (tenDiff)[3] = (dwiTenEstP)[3] - (tenTenP )[3], (tenDiff)[4] = (dwiTenEstP)[4] - (tenTenP)[4], (tenDiff )[5] = (dwiTenEstP)[5] - (tenTenP)[5], (tenDiff)[6] = (dwiTenEstP )[6] - (tenTenP)[6]); |
| 826 | TEN_T_LERP(tenAvg, 0.5, dwiTenEstP, tenTenP)( (tenAvg)[0] = (((0.5))*(((tenTenP)[0]) - ((dwiTenEstP)[0])) + ((dwiTenEstP)[0])), (tenAvg)[1] = (((0.5))*(((tenTenP)[1]) - ((dwiTenEstP)[1])) + ((dwiTenEstP)[1])), (tenAvg)[2] = ((( 0.5))*(((tenTenP)[2]) - ((dwiTenEstP)[2])) + ((dwiTenEstP)[2] )), (tenAvg)[3] = (((0.5))*(((tenTenP)[3]) - ((dwiTenEstP)[3] )) + ((dwiTenEstP)[3])), (tenAvg)[4] = (((0.5))*(((tenTenP)[4 ]) - ((dwiTenEstP)[4])) + ((dwiTenEstP)[4])), (tenAvg)[5] = ( ((0.5))*(((tenTenP)[5]) - ((dwiTenEstP)[5])) + ((dwiTenEstP)[ 5])), (tenAvg)[6] = (((0.5))*(((tenTenP)[6]) - ((dwiTenEstP)[ 6])) + ((dwiTenEstP)[6]))); |
| 827 | tenErr = TEN_T_NORM(tenDiff)(sqrt(( (tenDiff)[1]*(tenDiff)[1] + 2*(tenDiff)[2]*(tenDiff)[ 2] + 2*(tenDiff)[3]*(tenDiff)[3] + (tenDiff)[4]*(tenDiff)[4] + 2*(tenDiff)[5]*(tenDiff)[5] + (tenDiff)[6]*(tenDiff)[6] )))/TEN_T_NORM(tenAvg)(sqrt(( (tenAvg)[1]*(tenAvg)[1] + 2*(tenAvg)[2]*(tenAvg)[2] + 2*(tenAvg)[3]*(tenAvg)[3] + (tenAvg)[4]*(tenAvg)[4] + 2*(tenAvg )[5]*(tenAvg)[5] + (tenAvg)[6]*(tenAvg)[6] ))); |
| 828 | if (tenErr > tenErrMax) { |
| 829 | tenErrNum++; |
| 830 | if (tenErrNum > errNumMax) { |
| 831 | biffAddf(BKEY"probeSS", "%s: (probe %u) tenErr %g > %g too many times %u > %u", |
| 832 | me, probeIdx, tenErr, tenErrMax, tenErrNum, errNumMax); |
| 833 | return 1; |
| 834 | } |
| 835 | } |
| 836 | } |
| 837 | |
| 838 | /* see if invariant gradient learned from tensor volume agrees with |
| 839 | volume of pre-computed invariant gradients, and with |
| 840 | gradient of pre-computed invariants */ |
| 841 | if (AIR_EXISTS(vecErrMax)(((int)(!((vecErrMax) - (vecErrMax)))))) { |
| 842 | ELL_3V_SUB(vecDiff, sclGradP, vecVecP)((vecDiff)[0] = (sclGradP)[0] - (vecVecP)[0], (vecDiff)[1] = ( sclGradP)[1] - (vecVecP)[1], (vecDiff)[2] = (sclGradP)[2] - ( vecVecP)[2]); |
| 843 | ELL_3V_LERP(vecAvg, 0.5, sclGradP, vecVecP)((vecAvg)[0] = (((0.5))*(((vecVecP)[0]) - ((sclGradP)[0])) + ( (sclGradP)[0])), (vecAvg)[1] = (((0.5))*(((vecVecP)[1]) - ((sclGradP )[1])) + ((sclGradP)[1])), (vecAvg)[2] = (((0.5))*(((vecVecP) [2]) - ((sclGradP)[2])) + ((sclGradP)[2]))); |
| 844 | vecErr = ELL_3V_LEN(vecDiff)(sqrt((((vecDiff))[0]*((vecDiff))[0] + ((vecDiff))[1]*((vecDiff ))[1] + ((vecDiff))[2]*((vecDiff))[2])))/sqrt(ELL_3V_LEN(vecAvg)(sqrt((((vecAvg))[0]*((vecAvg))[0] + ((vecAvg))[1]*((vecAvg)) [1] + ((vecAvg))[2]*((vecAvg))[2])))); |
| 845 | if (vecErr > vecErrMax) { |
| 846 | vecErrNum += 0.5; |
| 847 | if (vecErrNum > errNumMax) { |
| 848 | biffAddf(BKEY"probeSS", "%s: (probe %u) (A) vecErr %g > %g too many times %g > %u", |
| 849 | me, probeIdx, vecErr, vecErrMax, vecErrNum, errNumMax); |
| 850 | return 1; |
| 851 | } |
| 852 | } |
| 853 | ELL_3V_SUB(vecDiff, tenTenNormGradP, vecVecP)((vecDiff)[0] = (tenTenNormGradP)[0] - (vecVecP)[0], (vecDiff )[1] = (tenTenNormGradP)[1] - (vecVecP)[1], (vecDiff)[2] = (tenTenNormGradP )[2] - (vecVecP)[2]); |
| 854 | ELL_3V_LERP(vecAvg, 0.5, tenTenNormGradP, vecVecP)((vecAvg)[0] = (((0.5))*(((vecVecP)[0]) - ((tenTenNormGradP)[ 0])) + ((tenTenNormGradP)[0])), (vecAvg)[1] = (((0.5))*(((vecVecP )[1]) - ((tenTenNormGradP)[1])) + ((tenTenNormGradP)[1])), (vecAvg )[2] = (((0.5))*(((vecVecP)[2]) - ((tenTenNormGradP)[2])) + ( (tenTenNormGradP)[2]))); |
| 855 | vecErr = ELL_3V_LEN(vecDiff)(sqrt((((vecDiff))[0]*((vecDiff))[0] + ((vecDiff))[1]*((vecDiff ))[1] + ((vecDiff))[2]*((vecDiff))[2])))/sqrt(ELL_3V_LEN(vecAvg)(sqrt((((vecAvg))[0]*((vecAvg))[0] + ((vecAvg))[1]*((vecAvg)) [1] + ((vecAvg))[2]*((vecAvg))[2])))); |
| 856 | if (vecErr > vecErrMax) { |
| 857 | vecErrNum += 0.5; |
| 858 | if (vecErrNum > errNumMax) { |
| 859 | biffAddf(BKEY"probeSS", "%s: (probe %u) (B) vecErr %g > %g too many times %g > %u", |
| 860 | me, probeIdx, vecErr, vecErrMax, vecErrNum, errNumMax); |
| 861 | return 1; |
| 862 | } |
| 863 | } |
| 864 | } |
| 865 | |
| 866 | /* see if invariant learned from tensor volume agrees with |
| 867 | volume of precomputed invariants */ |
| 868 | if (AIR_EXISTS(sclErrMax)(((int)(!((sclErrMax) - (sclErrMax)))))) { |
| 869 | sclDiff = sclSclP[0] - tenTenNormP[0]; |
| 870 | sclAvg = (sclSclP[0] + tenTenNormP[0])/2; |
| 871 | sclErr = AIR_ABS(sclDiff)((sclDiff) > 0.0f ? (sclDiff) : -(sclDiff))/sqrt(AIR_ABS(sclAvg)((sclAvg) > 0.0f ? (sclAvg) : -(sclAvg))); |
| 872 | if (sclErr > sclErrMax) { |
| 873 | sclErrNum++; |
| 874 | if (sclErrNum > errNumMax) { |
| 875 | biffAddf(BKEY"probeSS", "%s: (probe %u) (B) sclErr %g > %g too many times %u > %u", |
| 876 | me, probeIdx, sclErr, sclErrMax, sclErrNum, errNumMax); |
| 877 | return 1; |
| 878 | } |
| 879 | } |
| 880 | } |
| 881 | } |
| 882 | |
| 883 | return 0; |
| 884 | } |
| 885 | |
| 886 | static const char *testInfo = "for testing gage"; |
| 887 | |
| 888 | int |
| 889 | main(int argc, const char **argv) { |
| 890 | const char *me; |
| 891 | char *err = NULL((void*)0); |
| 892 | hestOpt *hopt=NULL((void*)0); |
| 893 | hestParm *hparm; |
| 894 | airArray *mop; |
| 895 | |
| 896 | const gageKind *kind[KIND_NUM4]; |
| 897 | char name[KIND_NUM4][AIR_STRLEN_SMALL(128+1)] = { "scl", "vec", "ten", "dwi" }; |
| 898 | char nameComp[KIND_NUM4][AIR_STRLEN_SMALL(128+1)] = { "sclComp", "vecComp", "tenComp", "dwiComp" }; |
| 899 | char *kernS; |
| 900 | gageKind *dwikind = NULL((void*)0); |
| 901 | gageContext *gctxComp[KIND_NUM4], *gctxCompCopy[KIND_NUM4], |
| 902 | *gctx[KIND_NUM4], *gctxCopy[KIND_NUM4]; |
| 903 | multiAnswer *manComp[KIND_NUM4], *man[KIND_NUM4]; |
| 904 | Nrrd *nin[KIND_NUM4], |
| 905 | /* these are the volumes that are used in gctxComp[] */ |
| 906 | *nsclCopy, *nvecComp[3], *nctenComp[7], **ndwiComp, |
| 907 | *ngrad; /* need access to list of gradients used to make DWIs; |
| 908 | (this is not the gradient of a scalar field) */ |
| 909 | double bval = 1000, noiseStdv=0.0001, ksupport; |
| 910 | unsigned int kindIdx, probeNum, |
| 911 | gradNum = 10; /* small number so that missing one will produce |
| 912 | a big reconstruction error */ |
| 913 | NrrdKernel *kpack[3]; |
| 914 | |
| 915 | kind[0] = gageKindScl; |
| 916 | kind[1] = gageKindVec; |
| 917 | kind[2] = tenGageKind; |
| 918 | kind[3] = NULL((void*)0); /* dwi */ |
| 919 | |
| 920 | me = argv[0]; |
| 921 | mop = airMopNew(); |
| 922 | hparm = hestParmNew(); |
| 923 | airMopAdd(mop, hparm, (airMopper)hestParmFree, airMopAlways); |
| 924 | /* learn things from hest */ |
| 925 | hestOptAdd(&hopt, "supp", "r", airTypeDouble, 1, 1, &ksupport, "1.0", |
| 926 | "kernel support"); |
| 927 | hestOptAdd(&hopt, "pnum", "N", airTypeUInt, 1, 1, &probeNum, "100", |
| 928 | "# of probes"); |
| 929 | hestOptAdd(&hopt, "k", "kern", airTypeString, 1, 1, &kernS, NULL((void*)0), |
| 930 | "kernel to use; can be: box, cos, or ctmr"); |
| 931 | hestParseOrDie(hopt, argc-1, argv+1, hparm, me, testInfo, |
| 932 | AIR_TRUE1, AIR_TRUE1, AIR_TRUE1); |
| 933 | airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways); |
| 934 | airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways); |
| 935 | |
| 936 | if (!strcmp(kernS, "box")) { |
| 937 | ELL_3V_SET(kpack,((kpack)[0] = (nrrdKernelBoxSupportDebug), (kpack)[1] = (nrrdKernelZero ), (kpack)[2] = (nrrdKernelZero)) |
| 938 | nrrdKernelBoxSupportDebug,((kpack)[0] = (nrrdKernelBoxSupportDebug), (kpack)[1] = (nrrdKernelZero ), (kpack)[2] = (nrrdKernelZero)) |
| 939 | nrrdKernelZero,((kpack)[0] = (nrrdKernelBoxSupportDebug), (kpack)[1] = (nrrdKernelZero ), (kpack)[2] = (nrrdKernelZero)) |
| 940 | nrrdKernelZero)((kpack)[0] = (nrrdKernelBoxSupportDebug), (kpack)[1] = (nrrdKernelZero ), (kpack)[2] = (nrrdKernelZero)); |
| 941 | } else if (!strcmp(kernS, "cos")) { |
| 942 | ELL_3V_SET(kpack,((kpack)[0] = (nrrdKernelCos4SupportDebug), (kpack)[1] = (nrrdKernelCos4SupportDebugD ), (kpack)[2] = (nrrdKernelCos4SupportDebugDD)) |
| 943 | nrrdKernelCos4SupportDebug,((kpack)[0] = (nrrdKernelCos4SupportDebug), (kpack)[1] = (nrrdKernelCos4SupportDebugD ), (kpack)[2] = (nrrdKernelCos4SupportDebugDD)) |
| 944 | nrrdKernelCos4SupportDebugD,((kpack)[0] = (nrrdKernelCos4SupportDebug), (kpack)[1] = (nrrdKernelCos4SupportDebugD ), (kpack)[2] = (nrrdKernelCos4SupportDebugDD)) |
| 945 | nrrdKernelCos4SupportDebugDD)((kpack)[0] = (nrrdKernelCos4SupportDebug), (kpack)[1] = (nrrdKernelCos4SupportDebugD ), (kpack)[2] = (nrrdKernelCos4SupportDebugDD)); |
| 946 | } else if (!strcmp(kernS, "ctmr")) { |
| 947 | ELL_3V_SET(kpack,((kpack)[0] = (nrrdKernelCatmullRomSupportDebug), (kpack)[1] = (nrrdKernelCatmullRomSupportDebugD), (kpack)[2] = (nrrdKernelCatmullRomSupportDebugDD )) |
| 948 | nrrdKernelCatmullRomSupportDebug,((kpack)[0] = (nrrdKernelCatmullRomSupportDebug), (kpack)[1] = (nrrdKernelCatmullRomSupportDebugD), (kpack)[2] = (nrrdKernelCatmullRomSupportDebugDD )) |
| 949 | nrrdKernelCatmullRomSupportDebugD,((kpack)[0] = (nrrdKernelCatmullRomSupportDebug), (kpack)[1] = (nrrdKernelCatmullRomSupportDebugD), (kpack)[2] = (nrrdKernelCatmullRomSupportDebugDD )) |
| 950 | nrrdKernelCatmullRomSupportDebugDD)((kpack)[0] = (nrrdKernelCatmullRomSupportDebug), (kpack)[1] = (nrrdKernelCatmullRomSupportDebugD), (kpack)[2] = (nrrdKernelCatmullRomSupportDebugDD )); |
| 951 | } else { |
| 952 | fprintf(stderr__stderrp, "%s: kernel \"%s\" not recognized", me, kernS); |
| 953 | airMopError(mop); return 1; |
| 954 | } |
| 955 | fprintf(stderr__stderrp, "kpack = %s, %s, %s\n", kpack[0]->name, kpack[1]->name, kpack[2]->name); |
| 956 | |
| 957 | /* This was a tricky bug: Adding gageContextNix(gctx) to the mop |
| 958 | as soon as a gctx is created makes sense. But, in the first |
| 959 | version of this code, gageContextNix was added to the mop |
| 960 | BEFORE tenDwiGageKindNix was added, which meant that gageContextNix |
| 961 | was being called AFTER tenDwiGageKindNix. However, that meant |
| 962 | that existing pvls had their pvl->kind being freed from under them, |
| 963 | but gagePerVolumeNix certainly needs to look at and possibly call |
| 964 | pvl->kind->pvlDataNix in order clean up pvl->data. *Therefore*, |
| 965 | we have to add tenDwiGageKindNix to the mop prior to gageContextNix, |
| 966 | even though nothing about the (dyanamic) kind has been set yet. |
| 967 | If we had something like smart pointers, then tenDwiGageKindNix |
| 968 | would not have freed something that an extant pvl was using */ |
| 969 | dwikind = tenDwiGageKindNew(); |
| 970 | airMopAdd(mop, dwikind, (airMopper)tenDwiGageKindNix, airMopAlways); |
| 971 | |
| 972 | #define GAGE_CTX_NEW(gg, mm) \ |
| 973 | (gg) = gageContextNew(); \ |
| 974 | airMopAdd((mm), (gg), (airMopper)gageContextNix, airMopAlways); \ |
| 975 | gageParmSet((gg), gageParmRenormalize, AIR_FALSE0); \ |
| 976 | gageParmSet((gg), gageParmCheckIntegrals, AIR_TRUE1) |
| 977 | |
| 978 | for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) { |
| 979 | GAGE_CTX_NEW(gctx[kindIdx], mop); |
| 980 | GAGE_CTX_NEW(gctxComp[kindIdx], mop); |
| 981 | /* |
| 982 | fprintf(stderr, "%s: kind %s (%u): gctx=%p, gctxComp=%p\n", me, |
| 983 | kindStr[kindIdx], kindIdx, |
| 984 | AIR_VOIDP(gctx[kindIdx]), |
| 985 | AIR_VOIDP(gctxComp[kindIdx])); |
| 986 | */ |
| 987 | man[kindIdx] = multiAnswerNew(name[kindIdx]); |
| 988 | manComp[kindIdx] = multiAnswerNew(nameComp[kindIdx]); |
| 989 | airMopAdd(mop, man[kindIdx], (airMopper)multiAnswerNix, airMopAlways); |
| 990 | airMopAdd(mop, manComp[kindIdx], (airMopper)multiAnswerNix, airMopAlways); |
| 991 | } |
| 992 | #undef GAGE_CTX_NEW |
| 993 | |
| 994 | for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) { |
| 995 | NRRD_NEW(nin[kindIdx], mop); |
| 996 | } |
| 997 | NRRD_NEW(ngrad, mop); |
| 998 | NRRD_NEW(nsclCopy, mop); |
| 999 | if (engageGenTensor(gctx[KI_TEN2], nin[KI_TEN2], |
| 1000 | /* out/in */ |
| 1001 | noiseStdv, 42 /* RNG seed, could be a command-line option */, |
| 1002 | volSize[0], volSize[1], volSize[2]) |
| 1003 | || engageGenScalar(gctx[KI_SCL0], nin[KI_SCL0], |
| 1004 | gctxComp[KI_SCL0], nsclCopy, |
| 1005 | /* out/in */ |
| 1006 | nin[KI_TEN2]) |
| 1007 | || engageGenVector(gctx[KI_VEC1], nin[KI_VEC1], |
| 1008 | /* out/in */ |
| 1009 | nin[KI_SCL0]) |
| 1010 | /* engage'ing of nin[KI_DWI] done below */ |
| 1011 | || genDwi(nin[KI_DWI3], ngrad, |
| 1012 | /* out/in */ |
| 1013 | gradNum /* for B0 */, bval, nin[KI_TEN2]) |
| 1014 | || engageMopDiceVector(gctxComp[KI_VEC1], nvecComp, |
| 1015 | /* out/in */ |
| 1016 | mop, nin[KI_VEC1]) |
| 1017 | || engageMopDiceTensor(gctxComp[KI_TEN2], nctenComp, |
| 1018 | /* out/in */ |
| 1019 | mop, nin[KI_TEN2]) |
| 1020 | || engageMopDiceDwi(gctxComp[KI_DWI3], &ndwiComp, |
| 1021 | /* out/in */ |
| 1022 | mop, nin[KI_DWI3])) { |
| 1023 | airMopAdd(mop, err = biffGetDone(BKEY"probeSS"), airFree, airMopAlways); |
| 1024 | fprintf(stderr__stderrp, "%s: trouble creating volumes:\n%s", me, err); |
| 1025 | airMopError(mop); return 1; |
| 1026 | } |
| 1027 | if (tenDwiGageKindSet(dwikind, -1 /* thresh */, 0 /* soft */, |
| 1028 | bval, 1 /* valueMin */, |
| 1029 | ngrad, NULL((void*)0), |
| 1030 | tenEstimate1MethodLLS, |
| 1031 | tenEstimate2MethodPeled, |
| 1032 | 424242)) { |
| 1033 | airMopAdd(mop, err = biffGetDone(TENtenBiffKey), airFree, airMopAlways); |
| 1034 | fprintf(stderr__stderrp, "%s: trouble creating DWI kind:\n%s", me, err); |
| 1035 | airMopError(mop); return 1; |
| 1036 | } |
| 1037 | /* access through kind[] is const, but not through dwikind */ |
| 1038 | kind[KI_DWI3] = dwikind; |
| 1039 | /* engage dwi vol */ |
| 1040 | { |
| 1041 | gagePerVolume *pvl; |
| 1042 | if ( !(pvl = gagePerVolumeNew(gctx[KI_DWI3], nin[KI_DWI3], kind[KI_DWI3])) |
| 1043 | || gagePerVolumeAttach(gctx[KI_DWI3], pvl) ) { |
| 1044 | airMopAdd(mop, err = biffGetDone(GAGEgageBiffKey), airFree, airMopAlways); |
| 1045 | fprintf(stderr__stderrp, "%s: trouble creating DWI context:\n%s", me, err); |
| 1046 | airMopError(mop); return 1; |
| 1047 | } |
| 1048 | /* |
| 1049 | fprintf(stderr, "%s: %s gctx=%p gets pvl %p\n", me, |
| 1050 | kindStr[KI_DWI], AIR_VOIDP(gctx[KI_DWI]), AIR_VOIDP(pvl)); |
| 1051 | */ |
| 1052 | } |
| 1053 | |
| 1054 | /* make sure kinds can parse back to themselves */ |
| 1055 | /* the messiness here is in part because of problems with the gage |
| 1056 | API, and the weirdness of how meetGageKindParse is either allocating |
| 1057 | something, or not, depending on its input. There is no way to |
| 1058 | refer to the "name" of a dwi kind without having allocated something. */ |
| 1059 | for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) { |
| 1060 | if (!(kind[kindIdx]->dynamicAlloc)) { |
| 1061 | if (kind[kindIdx] != meetGageKindParse(kind[kindIdx]->name)) { |
| 1062 | fprintf(stderr__stderrp, "%s: static kind[%u]->name %s wasn't parsed\n", me, |
| 1063 | kindIdx, kind[kindIdx]->name); |
| 1064 | airMopError(mop); return 1; |
| 1065 | } |
| 1066 | } else { |
| 1067 | gageKind *ktmp; |
| 1068 | ktmp = meetGageKindParse(kind[kindIdx]->name); |
| 1069 | if (!ktmp) { |
| 1070 | fprintf(stderr__stderrp, "%s: dynamic kind[%u]->name %s wasn't parsed\n", me, |
| 1071 | kindIdx, kind[kindIdx]->name); |
| 1072 | airMopError(mop); return 1; |
| 1073 | } |
| 1074 | if (airStrcmp(ktmp->name, kind[kindIdx]->name)) { |
| 1075 | fprintf(stderr__stderrp, "%s: parsed dynamic kind[%u]->name %s didn't match %s\n", me, |
| 1076 | kindIdx, ktmp->name, kind[kindIdx]->name); |
| 1077 | airMopError(mop); return 1; |
| 1078 | } |
| 1079 | if (!airStrcmp(TEN_DWI_GAGE_KIND_NAME"dwi", ktmp->name)) { |
| 1080 | /* HEY: there needs to be a generic way of freeing |
| 1081 | a dynamic kind, perhaps a new nixer field in gageKind */ |
| 1082 | ktmp = tenDwiGageKindNix(ktmp); |
Value stored to 'ktmp' is never read | |
| 1083 | } |
| 1084 | } |
| 1085 | } |
| 1086 | /* |
| 1087 | for (kindIdx=0; kindIdx<KIND_NUM; kindIdx++) { |
| 1088 | printf("%s: %s (%s) -> (%u) %u\n", me, kind[kindIdx]->name, |
| 1089 | kind[kindIdx]->dynamicAlloc ? "dynamic" : "static", |
| 1090 | kind[kindIdx]->baseDim, kind[kindIdx]->valLen); |
| 1091 | } |
| 1092 | */ |
| 1093 | |
| 1094 | /* |
| 1095 | nrrdSave("tmp-cten.nrrd", nin[KI_TEN], NULL); |
| 1096 | nrrdSave("tmp-scl.nrrd", nin[KI_SCL], NULL); |
| 1097 | nrrdSave("tmp-vec.nrrd", nin[KI_VEC], NULL); |
| 1098 | nrrdSave("tmp-dwi.nrrd", nin[KI_DWI], NULL); |
| 1099 | */ |
| 1100 | |
| 1101 | /* ========================== TASK 1 */ |
| 1102 | if (updateQueryKernelSetTask1(gctxComp, gctx, AIR_TRUE1, manComp, man, kpack, ksupport) |
| 1103 | || probeTask1(gctxComp, gctx, manComp, man, probeNum, kpack, ksupport)) { |
| 1104 | airMopAdd(mop, err = biffGetDone(BKEY"probeSS"), airFree, airMopAlways); |
| 1105 | fprintf(stderr__stderrp, "%s: trouble (orig, orig):\n%s", me, err); |
| 1106 | airMopError(mop); return 1; |
| 1107 | } |
| 1108 | /* testing gageContextCopy on multi-variate volumes */ |
| 1109 | for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) { |
| 1110 | if (!( gctxCopy[kindIdx] = gageContextCopy(gctx[kindIdx]) )) { |
| 1111 | airMopAdd(mop, err = biffGetDone(GAGEgageBiffKey), airFree, airMopAlways); |
| 1112 | fprintf(stderr__stderrp, "%s: trouble w/ multi-variate contexts:\n%s", me, err); |
| 1113 | airMopError(mop); return 1; |
| 1114 | } |
| 1115 | airMopAdd(mop, gctxCopy[kindIdx], (airMopper)gageContextNix, airMopAlways); |
| 1116 | } |
| 1117 | if (updateQueryKernelSetTask1(gctxComp, gctxCopy, AIR_FALSE0, manComp, man, kpack, ksupport) |
| 1118 | || probeTask1(gctxComp, gctxCopy, manComp, man, probeNum, kpack, ksupport)) { |
| 1119 | airMopAdd(mop, err = biffGetDone(BKEY"probeSS"), airFree, airMopAlways); |
| 1120 | fprintf(stderr__stderrp, "%s: trouble (orig, copy):\n%s", me, err); |
| 1121 | airMopError(mop); return 1; |
| 1122 | } |
| 1123 | for (kindIdx=0; kindIdx<KIND_NUM4; kindIdx++) { |
| 1124 | if (!( gctxCompCopy[kindIdx] = gageContextCopy(gctxComp[kindIdx]) )) { |
| 1125 | airMopAdd(mop, err = biffGetDone(GAGEgageBiffKey), airFree, airMopAlways); |
| 1126 | fprintf(stderr__stderrp, "%s: trouble w/ component contexts:\n%s", me, err); |
| 1127 | airMopError(mop); return 1; |
| 1128 | } |
| 1129 | airMopAdd(mop, gctxCompCopy[kindIdx], (airMopper)gageContextNix, airMopAlways); |
| 1130 | } |
| 1131 | if (updateQueryKernelSetTask1(gctxCompCopy, gctxCopy, AIR_FALSE0, manComp, man, kpack, ksupport) |
| 1132 | || probeTask1(gctxCompCopy, gctxCopy, manComp, man, probeNum, kpack, ksupport)) { |
| 1133 | airMopAdd(mop, err = biffGetDone(BKEY"probeSS"), airFree, airMopAlways); |
| 1134 | fprintf(stderr__stderrp, "%s: trouble (copy, copy):\n%s", me, err); |
| 1135 | airMopError(mop); return 1; |
| 1136 | } |
| 1137 | |
| 1138 | /* ========================== TASK 2 */ |
| 1139 | /* create scale-space stacks with tent, ctmr, and hermite */ |
| 1140 | /* save them, free, and read them back in */ |
| 1141 | /* gageContextCopy on stack */ |
| 1142 | /* pick a scale in-between tau samples */ |
| 1143 | /* for all the tau's half-way between tau samples in scale: |
| 1144 | blur at that tau to get correct values |
| 1145 | check that error with hermite is lower than ctmr is lower than tent */ |
| 1146 | /* for all tau samples: |
| 1147 | blur at that tau to (re-)get correct values |
| 1148 | check that everything agrees there */ |
| 1149 | |
| 1150 | |
| 1151 | /* single probe with high verbosity */ |
| 1152 | |
| 1153 | airMopOkay(mop); |
| 1154 | return 0; |
| 1155 | } |
| 1156 | #undef NRRD_NEW |