| File: | src/bin/gprobe.c |
| Location: | line 849, column 5 |
| Description: | Value stored to 'sox' 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 | |
| 25 | #include <stdio.h> |
| 26 | |
| 27 | #include <teem/biff.h> |
| 28 | #include <teem/hest.h> |
| 29 | #include <teem/nrrd.h> |
| 30 | #include <teem/gage.h> |
| 31 | #include <teem/ten.h> |
| 32 | #include <teem/meet.h> |
| 33 | |
| 34 | static void |
| 35 | printans(FILE *file, const double *ans, unsigned int len) { |
| 36 | unsigned int ai; |
| 37 | |
| 38 | AIR_UNUSED(file)(void)(file); |
| 39 | for (ai=0; ai<len; ai++) { |
| 40 | if (ai) { |
| 41 | printf(", "); |
| 42 | } |
| 43 | printf("%g", ans[ai]); |
| 44 | } |
| 45 | } |
| 46 | |
| 47 | static int |
| 48 | gridProbe(gageContext *ctx, gagePerVolume *pvl, int what, |
| 49 | Nrrd *nout, int typeOut, Nrrd *_ngrid, |
| 50 | int indexSpace, int verbose, int clamp, int scaleIsTau, |
| 51 | double eft, double eftVal) { |
| 52 | char me[]="gridProbe"; |
| 53 | Nrrd *ngrid; |
| 54 | airArray *mop; |
| 55 | double *grid, pos[4]; |
| 56 | const double *answer; |
| 57 | unsigned int ansLen, dim, aidx, baseDim, gridDim; |
| 58 | size_t sizeOut[NRRD_DIM_MAX16], coordOut[NRRD_DIM_MAX16], II, NN; |
| 59 | double (*ins)(void *v, size_t I, double d); |
| 60 | char stmp[2][AIR_STRLEN_SMALL(128+1)]; |
| 61 | |
| 62 | if (!(ctx && pvl && nout && _ngrid)) { |
| 63 | biffAddf(GAGEgageBiffKey, "%s: got NULL pointer", me); |
| 64 | return 1; |
| 65 | } |
| 66 | if (airEnumValCheck(nrrdType, typeOut)) { |
| 67 | biffAddf(GAGEgageBiffKey, "%s: type %d not valid", me, typeOut); |
| 68 | return 1; |
| 69 | } |
| 70 | if (!gagePerVolumeIsAttached(ctx, pvl)) { |
| 71 | biffAddf(GAGEgageBiffKey, "%s: given pvl not attached to context", me); |
| 72 | return 1; |
| 73 | } |
| 74 | if (!(2 == _ngrid->dim)) { |
| 75 | biffAddf(GAGEgageBiffKey, "%s: ngrid must be 2 (not %u)", me, _ngrid->dim); |
| 76 | return 1; |
| 77 | } |
| 78 | if ((ctx->stackPos && _ngrid->axis[0].size != 5) |
| 79 | || (!ctx->stackPos && _ngrid->axis[0].size != 4)) { |
| 80 | biffAddf(GAGEgageBiffKey, "%s: if %susing stack, need " |
| 81 | "ngrid->axis[0].size = %u = 1 + %u (not %u)", me, |
| 82 | (ctx->stackPos ? "" : "not "), |
| 83 | (ctx->stackPos ? 4 : 3) + 1, |
| 84 | (ctx->stackPos ? 4 : 3), |
| 85 | AIR_UINT(_ngrid->axis[0].size)((unsigned int)(_ngrid->axis[0].size))); |
| 86 | return 1; |
| 87 | } |
| 88 | |
| 89 | mop = airMopNew(); |
| 90 | ngrid = nrrdNew(); |
| 91 | airMopAdd(mop, ngrid, (airMopper)nrrdNuke, airMopAlways); |
| 92 | if (ctx->stackPos) { |
| 93 | if (nrrdConvert(ngrid, _ngrid, nrrdTypeDouble)) { |
| 94 | biffMovef(GAGEgageBiffKey, NRRDnrrdBiffKey, "%s: trouble converting ngrid", me); |
| 95 | airMopError(mop); return 1; |
| 96 | } |
| 97 | } else { |
| 98 | Nrrd *ntmp; |
| 99 | ptrdiff_t minIdx[2], maxIdx[2]; |
| 100 | minIdx[0] = minIdx[1] = 0; |
| 101 | maxIdx[0] = 4; /* pad by one sample */ |
| 102 | maxIdx[1] = AIR_CAST(ptrdiff_t, _ngrid->axis[1].size-1)((ptrdiff_t)(_ngrid->axis[1].size-1)); /* no padding */ |
| 103 | ntmp = nrrdNew(); |
| 104 | airMopAdd(mop, ntmp, (airMopper)nrrdNuke, airMopAlways); |
| 105 | if (nrrdConvert(ntmp, _ngrid, nrrdTypeDouble) |
| 106 | || nrrdPad_nva(ngrid, ntmp, minIdx, maxIdx, nrrdBoundaryPad, 0.0)) { |
| 107 | biffMovef(GAGEgageBiffKey, NRRDnrrdBiffKey, "%s: trouble converting/padding ngrid", me); |
| 108 | airMopError(mop); return 1; |
| 109 | } |
| 110 | } |
| 111 | grid = AIR_CAST(double *, ngrid->data)((double *)(ngrid->data)); |
| 112 | gridDim = AIR_ROUNDUP_UI(grid[0])((unsigned int)(floor((grid[0])+0.5))); |
| 113 | if (gridDim + 1 != ngrid->axis[1].size) { |
| 114 | biffAddf(GAGEgageBiffKey, "%s: ngrid->axis[1].size = %u but expected %u = 1 + %u", |
| 115 | me, AIR_UINT(ngrid->axis[1].size)((unsigned int)(ngrid->axis[1].size)), |
| 116 | 1 + gridDim, gridDim); |
| 117 | airMopError(mop); return 1; |
| 118 | } |
| 119 | answer = gageAnswerPointer(ctx, pvl, what); |
| 120 | ansLen = pvl->kind->table[what].answerLength; |
| 121 | baseDim = 1 == ansLen ? 0 : 1; |
| 122 | dim = baseDim + gridDim; |
| 123 | if (dim > NRRD_DIM_MAX16) { |
| 124 | biffAddf(GAGEgageBiffKey, "%s: output dimension %u unreasonable", me, dim); |
| 125 | airMopError(mop); return 1; |
| 126 | } |
| 127 | if (ansLen > 1) { |
| 128 | sizeOut[0] = ansLen; |
| 129 | coordOut[0] = 0; |
| 130 | } |
| 131 | NN = 1; |
| 132 | for (aidx=0; aidx<gridDim; aidx++) { |
| 133 | sizeOut[aidx + baseDim] = AIR_ROUNDUP_UI(grid[0 + 5*(aidx+1)])((unsigned int)(floor((grid[0 + 5*(aidx+1)])+0.5))); |
| 134 | NN *= sizeOut[aidx + baseDim]; |
| 135 | coordOut[aidx + baseDim] = 0; |
| 136 | } |
| 137 | if (nrrdMaybeAlloc_nva(nout, typeOut, dim, sizeOut)) { |
| 138 | biffMovef(GAGEgageBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate output", me); |
| 139 | airMopError(mop); return 1; |
| 140 | } |
| 141 | ins = nrrdDInsert[nout->type]; |
| 142 | for (II=0; II<NN; II++) { |
| 143 | int E; |
| 144 | if (verbose && 3 == gridDim |
| 145 | && !coordOut[0+baseDim] && !coordOut[1+baseDim]) { |
| 146 | if (verbose > 1) { |
| 147 | fprintf(stderr__stderrp, "z = "); |
| 148 | } |
| 149 | fprintf(stderr__stderrp, " %s/%s", |
| 150 | airSprintSize_t(stmp[0], coordOut[2+baseDim]), |
| 151 | airSprintSize_t(stmp[1], sizeOut[2+baseDim])); |
| 152 | fflush(stderr__stderrp); |
| 153 | if (verbose > 1) { |
| 154 | fprintf(stderr__stderrp, "\n"); |
| 155 | } |
| 156 | } |
| 157 | ELL_4V_COPY(pos, grid + 1 + 5*0)((pos)[0] = (grid + 1 + 5*0)[0], (pos)[1] = (grid + 1 + 5*0)[ 1], (pos)[2] = (grid + 1 + 5*0)[2], (pos)[3] = (grid + 1 + 5* 0)[3]); |
| 158 | for (aidx=0; aidx<gridDim; aidx++) { |
| 159 | ELL_4V_SCALE_ADD2(pos, 1, pos,((pos)[0] = (1)*(pos)[0] + (((double)(coordOut[aidx + baseDim ])))*(grid + 1 + 5*(1+aidx))[0], (pos)[1] = (1)*(pos)[1] + (( (double)(coordOut[aidx + baseDim])))*(grid + 1 + 5*(1+aidx))[ 1], (pos)[2] = (1)*(pos)[2] + (((double)(coordOut[aidx + baseDim ])))*(grid + 1 + 5*(1+aidx))[2], (pos)[3] = (1)*(pos)[3] + (( (double)(coordOut[aidx + baseDim])))*(grid + 1 + 5*(1+aidx))[ 3]) |
| 160 | AIR_CAST(double, coordOut[aidx + baseDim]),((pos)[0] = (1)*(pos)[0] + (((double)(coordOut[aidx + baseDim ])))*(grid + 1 + 5*(1+aidx))[0], (pos)[1] = (1)*(pos)[1] + (( (double)(coordOut[aidx + baseDim])))*(grid + 1 + 5*(1+aidx))[ 1], (pos)[2] = (1)*(pos)[2] + (((double)(coordOut[aidx + baseDim ])))*(grid + 1 + 5*(1+aidx))[2], (pos)[3] = (1)*(pos)[3] + (( (double)(coordOut[aidx + baseDim])))*(grid + 1 + 5*(1+aidx))[ 3]) |
| 161 | grid + 1 + 5*(1+aidx))((pos)[0] = (1)*(pos)[0] + (((double)(coordOut[aidx + baseDim ])))*(grid + 1 + 5*(1+aidx))[0], (pos)[1] = (1)*(pos)[1] + (( (double)(coordOut[aidx + baseDim])))*(grid + 1 + 5*(1+aidx))[ 1], (pos)[2] = (1)*(pos)[2] + (((double)(coordOut[aidx + baseDim ])))*(grid + 1 + 5*(1+aidx))[2], (pos)[3] = (1)*(pos)[3] + (( (double)(coordOut[aidx + baseDim])))*(grid + 1 + 5*(1+aidx))[ 3]); |
| 162 | } |
| 163 | if (scaleIsTau && ctx->stackPos) { |
| 164 | /* have to convert given tau values to sigma */ |
| 165 | pos[3] = airSigmaOfTau(pos[3]); |
| 166 | } |
| 167 | /* |
| 168 | printf("%s: %u -> (%u %u) -> %g %g %g %g (%s)\n", me, |
| 169 | AIR_UINT(II), |
| 170 | AIR_UINT(coordOut[0+baseDim]), |
| 171 | AIR_UINT(coordOut[1+baseDim]), |
| 172 | pos[0], pos[1], pos[2], pos[3], |
| 173 | indexSpace ? "index" : "world"); |
| 174 | */ |
| 175 | E = (ctx->stackPos |
| 176 | ? gageStackProbeSpace(ctx, pos[0], pos[1], pos[2], pos[3], |
| 177 | indexSpace, clamp) |
| 178 | : gageProbeSpace(ctx, pos[0], pos[1], pos[2], |
| 179 | indexSpace, clamp)); |
| 180 | if (E) { |
| 181 | biffAddf(GAGEgageBiffKey, "%s: trouble at II=%s =(%g,%g,%g,%g):\n%s\n(%d)\n", me, |
| 182 | airSprintSize_t(stmp[0], II), |
| 183 | pos[0], pos[1], pos[2], pos[3], |
| 184 | ctx->errStr, ctx->errNum); |
| 185 | airMopError(mop); return 1; |
| 186 | } |
| 187 | if (1 == ansLen) { |
| 188 | ins(nout->data, II, (ctx->edgeFrac > eft |
| 189 | ? eftVal |
| 190 | : *answer)); |
| 191 | } else { |
| 192 | for (aidx=0; aidx<ansLen; aidx++) { |
| 193 | ins(nout->data, aidx + ansLen*II, (ctx->edgeFrac > eft |
| 194 | ? eftVal |
| 195 | : answer[aidx])); |
| 196 | } |
| 197 | } |
| 198 | NRRD_COORD_INCR(coordOut, sizeOut, dim, baseDim)if ((baseDim) < (dim)) { (coordOut)[(baseDim)]++; { unsigned int ddd; for (ddd=0; ddd+1 < ((dim)-(baseDim)) && ((coordOut)+(baseDim))[ddd] >= ((sizeOut)+(baseDim))[ddd] ; ddd++) { ((coordOut)+(baseDim))[ddd] = 0; ((coordOut)+(baseDim ))[ddd+1]++; } if ((dim)-(baseDim)) { ((coordOut)+(baseDim))[ ((dim)-(baseDim))-1] = ((((coordOut)+(baseDim))[((dim)-(baseDim ))-1]) < (((sizeOut)+(baseDim))[((dim)-(baseDim))-1]-1) ? ( ((coordOut)+(baseDim))[((dim)-(baseDim))-1]) : (((sizeOut)+(baseDim ))[((dim)-(baseDim))-1]-1)); } }; }; |
| 199 | } |
| 200 | if (verbose && verbose <= 1) { |
| 201 | fprintf(stderr__stderrp, "\n"); |
| 202 | } |
| 203 | |
| 204 | if (!indexSpace) { |
| 205 | /* set the output space directions, but (being conservative/cautious) |
| 206 | only do so when grid had specified world-space positions */ |
| 207 | /* HEY: untested! whipped up out of frustration for GLK Bonn talk */ |
| 208 | nout->spaceDim = 3; |
| 209 | if (baseDim) { |
| 210 | nrrdSpaceVecSetNaN(nout->axis[0].spaceDirection); |
| 211 | } |
| 212 | nrrdSpaceVecCopy(nout->spaceOrigin, grid + 1 + 5*0); |
| 213 | for (aidx=0; aidx<gridDim; aidx++) { |
| 214 | nrrdSpaceVecCopy(nout->axis[baseDim+aidx].spaceDirection, |
| 215 | grid + 1 + 5*(1+aidx)); |
| 216 | } |
| 217 | } |
| 218 | airMopOkay(mop); |
| 219 | return 0; |
| 220 | } |
| 221 | |
| 222 | static const char *probeInfo = |
| 223 | ("Shows off the functionality of the gage library. " |
| 224 | "Uses gageProbe() to query various kinds of volumes " |
| 225 | "to learn various measured or derived quantities."); |
| 226 | |
| 227 | int |
| 228 | main(int argc, const char *argv[]) { |
| 229 | gageKind *kind; |
| 230 | const char *me; |
| 231 | char *whatS, *err, *outS, *stackFnameFormat; |
| 232 | hestParm *hparm; |
| 233 | hestOpt *hopt = NULL((void*)0); |
| 234 | NrrdKernelSpec *k00, *k11, *k22, *kSS, *kSSblur; |
| 235 | int what, E=0, renorm, uniformSS, optimSS, verbose, zeroZ, |
| 236 | orientationFromSpacing, probeSpaceIndex, normdSS; |
| 237 | unsigned int iBaseDim, oBaseDim, axi, numSS, seed; |
| 238 | const double *answer; |
| 239 | Nrrd *nin, *_npos, *npos, *_ngrid, *ngrid, *nout, **ninSS=NULL((void*)0); |
| 240 | Nrrd *ngrad=NULL((void*)0), *nbmat=NULL((void*)0); |
| 241 | size_t six, siy, siz, sox, soy, soz; |
| 242 | double bval=0, eps, gmc, rangeSS[2], *pntPos, scale[3], posSS, biasSS, |
| 243 | dsix, dsiy, dsiz, dsox, dsoy, dsoz, edgeFracInfo[2]; |
| 244 | gageContext *ctx; |
| 245 | gagePerVolume *pvl=NULL((void*)0); |
| 246 | double t0, t1, rscl[3], min[3], maxOut[3], maxIn[3]; |
| 247 | airArray *mop; |
| 248 | #define NON_SBP_OPT_NUM5 5 |
| 249 | unsigned int ansLen, *skip, skipNum, pntPosNum, |
| 250 | nonSbpOpi[NON_SBP_OPT_NUM5], nsi; |
| 251 | gageStackBlurParm *sbpIN, *sbpCL, *sbp; |
| 252 | int otype, clamp, scaleIsTau; |
| 253 | char stmp[4][AIR_STRLEN_SMALL(128+1)]; |
| 254 | |
| 255 | me = argv[0]; |
| 256 | /* parse environment variables first, in case they break nrrdDefault* |
| 257 | or nrrdState* variables in a way that nrrdSanity() should see */ |
| 258 | nrrdDefaultGetenv(); |
| 259 | nrrdStateGetenv(); |
| 260 | /* no harm done in making sure we're sane */ |
| 261 | if (!nrrdSanity()) { |
| 262 | fprintf(stderr__stderrp, "******************************************\n"); |
| 263 | fprintf(stderr__stderrp, "******************************************\n"); |
| 264 | fprintf(stderr__stderrp, "\n"); |
| 265 | fprintf(stderr__stderrp, " %s: nrrd sanity check FAILED.\n", me); |
| 266 | fprintf(stderr__stderrp, "\n"); |
| 267 | fprintf(stderr__stderrp, " This means that either nrrd can't work on this " |
| 268 | "platform, or (more likely)\n"); |
| 269 | fprintf(stderr__stderrp, " there was an error in the compilation options " |
| 270 | "and variable definitions\n"); |
| 271 | fprintf(stderr__stderrp, " for how Teem was built here.\n"); |
| 272 | fprintf(stderr__stderrp, "\n"); |
| 273 | fprintf(stderr__stderrp, " %s\n", err = biffGetDone(NRRDnrrdBiffKey)); |
| 274 | fprintf(stderr__stderrp, "\n"); |
| 275 | fprintf(stderr__stderrp, "******************************************\n"); |
| 276 | fprintf(stderr__stderrp, "******************************************\n"); |
| 277 | free(err); |
| 278 | return 1; |
| 279 | } |
| 280 | |
| 281 | mop = airMopNew(); |
| 282 | hparm = hestParmNew(); |
| 283 | airMopAdd(mop, hparm, AIR_CAST(airMopper, hestParmFree)((airMopper)(hestParmFree)), airMopAlways); |
| 284 | hparm->elideSingleOtherType = AIR_TRUE1; |
| 285 | hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &nin, NULL((void*)0), |
| 286 | "input volume", NULL((void*)0), NULL((void*)0), nrrdHestNrrd); |
| 287 | hestOptAdd(&hopt, "k", "kind", airTypeOther, 1, 1, &kind, NULL((void*)0), |
| 288 | "\"kind\" of volume (\"scalar\", \"vector\", " |
| 289 | "\"tensor\", or \"dwi\")", |
| 290 | NULL((void*)0), NULL((void*)0), meetHestGageKind); |
| 291 | hestOptAdd(&hopt, "v", "verbosity", airTypeInt, 1, 1, &verbose, "1", |
| 292 | "verbosity level"); |
| 293 | hestOptAdd(&hopt, "q", "query", airTypeString, 1, 1, &whatS, NULL((void*)0), |
| 294 | "the quantity (scalar, vector, or matrix) to learn by probing"); |
| 295 | hestOptAdd(&hopt, "gmc", "min gradmag", airTypeDouble, 1, 1, &gmc, |
| 296 | "0.0", "For curvature-based queries, use zero when gradient " |
| 297 | "magnitude is below this"); |
| 298 | hestOptAdd(&hopt, "ofs", "ofs", airTypeInt, 0, 0, &orientationFromSpacing, |
| 299 | NULL((void*)0), "If only per-axis spacing is available, use that to " |
| 300 | "contrive full orientation info"); |
| 301 | hestOptAdd(&hopt, "seed", "N", airTypeUInt, 1, 1, &seed, "42", |
| 302 | "RNG seed; mostly for debugging"); |
| 303 | hestOptAdd(&hopt, "c", "bool", airTypeBool, 1, 1, &clamp, "false", |
| 304 | "clamp positions as part of probing"); |
| 305 | hestOptAdd(&hopt, "ev", "thresh val", airTypeDouble, 2, 2, |
| 306 | edgeFracInfo, "1 0", |
| 307 | "if using position clamping (with \"-c true\"), the fraction " |
| 308 | "of values invented for the kernel support, or \"edge frac\" " |
| 309 | "is saved per probe (0 means kernel support was entirely within " |
| 310 | "data; 1 means everything was invented). " |
| 311 | "If this frac exceeds the first \"thresh\" " |
| 312 | "value given here, then the saved value for the probe will be " |
| 313 | "the second value \"val\" given here"); |
| 314 | hestOptAdd(&hopt, "zz", "bool", airTypeBool, 1, 1, &zeroZ, "false", |
| 315 | "enable \"zeroZ\" behavior in gage that partially " |
| 316 | "implements working with 3D images as if they are 2D"); |
| 317 | |
| 318 | hestOptAdd(&hopt, "k00", "kern00", airTypeOther, 1, 1, &k00, |
| 319 | "tent", "kernel for gageKernel00", |
| 320 | NULL((void*)0), NULL((void*)0), nrrdHestKernelSpec); |
| 321 | hestOptAdd(&hopt, "k11", "kern11", airTypeOther, 1, 1, &k11, |
| 322 | "cubicd:1,0", "kernel for gageKernel11", |
| 323 | NULL((void*)0), NULL((void*)0), nrrdHestKernelSpec); |
| 324 | hestOptAdd(&hopt, "k22", "kern22", airTypeOther, 1, 1, &k22, |
| 325 | "cubicdd:1,0", "kernel for gageKernel22", |
| 326 | NULL((void*)0), NULL((void*)0), nrrdHestKernelSpec); |
| 327 | hestOptAdd(&hopt, "rn", NULL((void*)0), airTypeInt, 0, 0, &renorm, NULL((void*)0), |
| 328 | "renormalize kernel weights at each new sample location. " |
| 329 | "\"Accurate\" kernels don't need this; doing it always " |
| 330 | "makes things go slower"); |
| 331 | |
| 332 | nsi = 0; |
| 333 | nonSbpOpi[nsi++] = |
| 334 | hestOptAdd(&hopt, "ssn", "SS #", airTypeUInt, 1, 1, &numSS, |
| 335 | "0", "how many scale-space samples to evaluate, or, " |
| 336 | "0 to turn-off all scale-space behavior"); |
| 337 | nonSbpOpi[nsi++] = |
| 338 | hestOptAdd(&hopt, "ssr", "scale range", airTypeDouble, 2, 2, rangeSS, |
| 339 | "nan nan", "range of scales in scale-space"); |
| 340 | nonSbpOpi[nsi++] = |
| 341 | hestOptAdd(&hopt, "ssu", NULL((void*)0), airTypeInt, 0, 0, &uniformSS, NULL((void*)0), |
| 342 | "do uniform samples along sigma, and not (by default) " |
| 343 | "samples according to the effective diffusion scale"); |
| 344 | nonSbpOpi[nsi++] = |
| 345 | hestOptAdd(&hopt, "sso", NULL((void*)0), airTypeInt, 0, 0, &optimSS, NULL((void*)0), |
| 346 | "if not using \"-ssu\", use pre-computed optimal " |
| 347 | "sigmas when possible"); |
| 348 | nonSbpOpi[nsi++] = |
| 349 | hestOptAdd(&hopt, "kssb", "kernel", airTypeOther, 1, 1, &kSSblur, |
| 350 | "dgauss:1,5", "blurring kernel, to sample scale space", |
| 351 | NULL((void*)0), NULL((void*)0), nrrdHestKernelSpec); |
| 352 | if (nsi != NON_SBP_OPT_NUM5) { |
| 353 | fprintf(stderr__stderrp, "%s: PANIC nsi %u != %u", me, nsi, NON_SBP_OPT_NUM5); |
| 354 | exit(1); |
| 355 | } |
| 356 | hestOptAdd(&hopt, "sbp", "blur spec", airTypeOther, 1, 1, &sbpCL, "", |
| 357 | "complete specification of stack blur parms; " |
| 358 | "over-rides all previous \"ss\" options", |
| 359 | NULL((void*)0), NULL((void*)0), gageHestStackBlurParm); |
| 360 | /* These two options are needed even if sbp is used, because they are *not* |
| 361 | part of the gageStackBlurParm. In meet, this info is handled by the |
| 362 | extraFlag/extraParm construct, which is not available here */ |
| 363 | hestOptAdd(&hopt, "ssnd", NULL((void*)0), airTypeInt, 0, 0, &normdSS, NULL((void*)0), |
| 364 | "normalize derivatives by scale"); |
| 365 | hestOptAdd(&hopt, "ssnb", "bias", airTypeDouble, 1, 1, &biasSS, "0.0", |
| 366 | "bias on scale-based derivative normalization"); |
| 367 | |
| 368 | hestOptAdd(&hopt, "ssf", "SS read/save format", airTypeString, 1, 1, |
| 369 | &stackFnameFormat, "", |
| 370 | "printf-style format (including a \"%u\") for the " |
| 371 | "filenames from which to read " |
| 372 | "pre-blurred volumes computed for the stack, if they " |
| 373 | "exist and match the stack parameters, and where to save " |
| 374 | "them if they had to be re-computed. Leave this as empty " |
| 375 | "string to disable this."); |
| 376 | |
| 377 | hestOptAdd(&hopt, "kssr", "kernel", airTypeOther, 1, 1, &kSS, |
| 378 | "hermite", "kernel for reconstructing from scale space samples", |
| 379 | NULL((void*)0), NULL((void*)0), nrrdHestKernelSpec); |
| 380 | hestOptAdd(&hopt, "sit", "sit", airTypeBool, 1, 1, &scaleIsTau, "false", |
| 381 | "in some places, scale should be interpreted as tau, not " |
| 382 | "sigma. Currently limited to grid probing (via \"-pg\")"); |
| 383 | |
| 384 | hestOptAdd(&hopt, "s", "sclX sclY sxlZ", airTypeDouble, 3, 3, scale, |
| 385 | "1 1 1", |
| 386 | "scaling factor for resampling on each axis " |
| 387 | "(>1.0 : supersampling); use \"-ssp\" (and \"-psi\")" |
| 388 | "to specify scale position of sampling"); |
| 389 | hestOptAdd(&hopt, "ssp", "pos", airTypeDouble, 1, 1, &posSS, "0", |
| 390 | "when using scale-space, scale-position at which to probe"); |
| 391 | hestOptAdd(&hopt, "pg", "nrrd", airTypeOther, 1, 1, &_ngrid, "", |
| 392 | "overrides \"-s\": " |
| 393 | "2-D nrrd which specifies origin and direction vectors " |
| 394 | "for sampling grid", NULL((void*)0), NULL((void*)0), nrrdHestNrrd); |
| 395 | hestOptAdd(&hopt, "pi", "nrrd", airTypeOther, 1, 1, &_npos, "", |
| 396 | "overrides \"-pv\": probes at this list of 3-vec or " |
| 397 | "4-vec positions", NULL((void*)0), NULL((void*)0), nrrdHestNrrd); |
| 398 | hestOptAdd(&hopt, "pp", "pos", airTypeDouble, 3, 4, &pntPos, |
| 399 | "nan nan nan", |
| 400 | "overrides \"-pi\": only sample at this specified point", |
| 401 | &pntPosNum); |
| 402 | hestOptAdd(&hopt, "eps", "epsilon", airTypeDouble, 1, 1, &eps, "0", |
| 403 | "if non-zero, and if query is a scalar, and if using \"pp\" " |
| 404 | "to query at a single point, then do epsilon offset probes " |
| 405 | "to calculate discrete differences, to find the numerical " |
| 406 | "gradient and hessian (for debugging)"); |
| 407 | hestOptAdd(&hopt, "psi", "p", airTypeBool, 1, 1, &probeSpaceIndex, "false", |
| 408 | "whether the probe location specification (by any of " |
| 409 | "the four previous flags) are in index space"); |
| 410 | |
| 411 | hestOptAdd(&hopt, "t", "type", airTypeEnum, 1, 1, &otype, "float", |
| 412 | "type of output volume", NULL((void*)0), nrrdType); |
| 413 | hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-", |
| 414 | "output volume"); |
| 415 | hestParseOrDie(hopt, argc-1, argv+1, hparm, |
| 416 | me, probeInfo, AIR_TRUE1, AIR_TRUE1, AIR_TRUE1); |
| 417 | airMopAdd(mop, hopt, AIR_CAST(airMopper, hestOptFree)((airMopper)(hestOptFree)), airMopAlways); |
| 418 | airMopAdd(mop, hopt, AIR_CAST(airMopper, hestParseFree)((airMopper)(hestParseFree)), airMopAlways); |
| 419 | |
| 420 | what = airEnumVal(kind->enm, whatS); |
| 421 | if (!what) { |
| 422 | /* 0 indeed always means "unknown" for any gageKind */ |
| 423 | fprintf(stderr__stderrp, "%s: couldn't parse \"%s\" as measure of \"%s\" volume\n", |
| 424 | me, whatS, kind->name); |
| 425 | hestUsage(stderr__stderrp, hopt, me, hparm); |
| 426 | hestGlossary(stderr__stderrp, hopt, hparm); |
| 427 | airMopError(mop); return 1; |
| 428 | } |
| 429 | |
| 430 | /* special set-up required for DWI kind */ |
| 431 | if (!strcmp(TEN_DWI_GAGE_KIND_NAME"dwi", kind->name)) { |
| 432 | if (tenDWMRIKeyValueParse(&ngrad, &nbmat, &bval, &skip, &skipNum, nin)) { |
| 433 | airMopAdd(mop, err = biffGetDone(TENtenBiffKey), airFree, airMopAlways); |
| 434 | fprintf(stderr__stderrp, "%s: trouble parsing DWI info:\n%s\n", me, err); |
| 435 | airMopError(mop); return 1; |
| 436 | } |
| 437 | if (skipNum) { |
| 438 | fprintf(stderr__stderrp, "%s: sorry, can't do DWI skipping in tenDwiGage", me); |
| 439 | airMopError(mop); return 1; |
| 440 | } |
| 441 | /* this could stand to use some more command-line arguments */ |
| 442 | if (tenDwiGageKindSet(kind, 50, 1, bval, 0.001, ngrad, nbmat, |
| 443 | tenEstimate1MethodLLS, |
| 444 | tenEstimate2MethodQSegLLS, seed)) { |
| 445 | airMopAdd(mop, err = biffGetDone(TENtenBiffKey), airFree, airMopAlways); |
| 446 | fprintf(stderr__stderrp, "%s: trouble parsing DWI info:\n%s\n", me, err); |
| 447 | airMopError(mop); return 1; |
| 448 | } |
| 449 | } |
| 450 | |
| 451 | /* for setting up pre-blurred scale-space samples */ |
| 452 | if (numSS || sbpCL) { |
| 453 | unsigned int vi; |
| 454 | int recompute, gotOld; |
| 455 | |
| 456 | if (sbpCL) { |
| 457 | /* we got the whole stack blar parm here */ |
| 458 | gotOld = AIR_FALSE0; |
| 459 | for (nsi=0; nsi<NON_SBP_OPT_NUM5; nsi++) { |
| 460 | gotOld |= (hestSourceUser == hopt[nonSbpOpi[nsi]].source); |
| 461 | } |
| 462 | if (gotOld) { |
| 463 | fprintf(stderr__stderrp, "%s: with new -sbp option; can't also use older " |
| 464 | "scale-space options (used", me); |
| 465 | for (nsi=0; nsi<NON_SBP_OPT_NUM5; nsi++) { |
| 466 | if (hestSourceUser == hopt[nonSbpOpi[nsi]].source) { |
| 467 | fprintf(stderr__stderrp, " -%s", hopt[nonSbpOpi[nsi]].flag); |
| 468 | } |
| 469 | } |
| 470 | fprintf(stderr__stderrp, ")\n"); |
| 471 | airMopError(mop); return 1; |
| 472 | } |
| 473 | if (gageStackBlurManage(&ninSS, &recompute, sbpCL, |
| 474 | stackFnameFormat, AIR_TRUE1, NULL((void*)0), |
| 475 | nin, kind)) { |
| 476 | airMopAdd(mop, err = biffGetDone(GAGEgageBiffKey), airFree, airMopAlways); |
| 477 | fprintf(stderr__stderrp, "%s: trouble getting volume stack:\n%s\n", me, err); |
| 478 | airMopError(mop); return 1; |
| 479 | } |
| 480 | if (sbpCL->verbose > 2) { |
| 481 | fprintf(stderr__stderrp, "%s: sampling scale range %g--%g via %s:\n", me, |
| 482 | sbpCL->sigmaRange[0], sbpCL->sigmaRange[1], |
| 483 | airEnumStr(gageSigmaSampling, sbpCL->sigmaSampling)); |
| 484 | for (vi=0; vi<sbpCL->num; vi++) { |
| 485 | fprintf(stderr__stderrp, " sigma[%u] = %g\n", vi, sbpCL->sigma[vi]); |
| 486 | } |
| 487 | } |
| 488 | sbp = sbpCL; |
| 489 | } else { |
| 490 | /* old way of doing things; depending on many separate options |
| 491 | to set numSS, rangeSS, uniformSS, optimSS, etc */ |
| 492 | sbpIN = gageStackBlurParmNew(); |
| 493 | airMopAdd(mop, sbpIN, (airMopper)gageStackBlurParmNix, airMopAlways); |
| 494 | if (gageStackBlurParmVerboseSet(sbpIN, verbose) |
| 495 | || gageStackBlurParmScaleSet(sbpIN, numSS, rangeSS[0], rangeSS[1], |
| 496 | uniformSS, optimSS) |
| 497 | || gageStackBlurParmKernelSet(sbpIN, kSSblur) |
| 498 | || gageStackBlurParmRenormalizeSet(sbpIN, AIR_TRUE1) |
| 499 | || gageStackBlurParmBoundarySet(sbpIN, nrrdBoundaryBleed, AIR_NAN(airFloatQNaN.f)) |
| 500 | || gageStackBlurManage(&ninSS, &recompute, sbpIN, |
| 501 | stackFnameFormat, AIR_TRUE1, NULL((void*)0), |
| 502 | nin, kind)) { |
| 503 | airMopAdd(mop, err = biffGetDone(GAGEgageBiffKey), airFree, airMopAlways); |
| 504 | fprintf(stderr__stderrp, "%s: trouble getting volume stack:\n%s\n", me, err); |
| 505 | airMopError(mop); return 1; |
| 506 | } |
| 507 | if (verbose > 2) { |
| 508 | fprintf(stderr__stderrp, "%s: sampling scale range %g--%g %suniformly:\n", me, |
| 509 | rangeSS[0], rangeSS[1], uniformSS ? "" : "non-"); |
| 510 | for (vi=0; vi<numSS; vi++) { |
| 511 | fprintf(stderr__stderrp, " scalePos[%u] = %g\n", vi, sbpIN->sigma[vi]); |
| 512 | } |
| 513 | } |
| 514 | sbp = sbpIN; |
| 515 | } |
| 516 | airMopAdd(mop, ninSS, airFree, airMopAlways); |
| 517 | } else { |
| 518 | ninSS = NULL((void*)0); |
| 519 | sbpIN = NULL((void*)0); |
| 520 | sbp = NULL((void*)0); |
| 521 | } |
| 522 | |
| 523 | /*** |
| 524 | **** Except for the gageProbe() call in the inner loop below, |
| 525 | **** and the gageContextNix() call at the very end, all the gage |
| 526 | **** calls which set up (and take down) the context and state are here. |
| 527 | ***/ |
| 528 | ctx = gageContextNew(); |
| 529 | airMopAdd(mop, ctx, AIR_CAST(airMopper, gageContextNix)((airMopper)(gageContextNix)), airMopAlways); |
| 530 | gageParmSet(ctx, gageParmGradMagCurvMin, gmc); |
| 531 | gageParmSet(ctx, gageParmVerbose, verbose); |
| 532 | gageParmSet(ctx, gageParmTwoDimZeroZ, zeroZ); |
| 533 | gageParmSet(ctx, gageParmRenormalize, renorm ? AIR_TRUE1 : AIR_FALSE0); |
| 534 | gageParmSet(ctx, gageParmCheckIntegrals, AIR_TRUE1); |
| 535 | gageParmSet(ctx, gageParmOrientationFromSpacing, orientationFromSpacing); |
| 536 | E = 0; |
| 537 | if (!E) E |= !(pvl = gagePerVolumeNew(ctx, nin, kind)); |
| 538 | if (!E) E |= gageKernelSet(ctx, gageKernel00, k00->kernel, k00->parm); |
| 539 | if (!E) E |= gageKernelSet(ctx, gageKernel11, k11->kernel, k11->parm); |
| 540 | if (!E) E |= gageKernelSet(ctx, gageKernel22, k22->kernel, k22->parm); |
| 541 | if (sbp) { |
| 542 | gagePerVolume **pvlSS; |
| 543 | gageParmSet(ctx, gageParmStackUse, AIR_TRUE1); |
| 544 | gageParmSet(ctx, gageParmStackNormalizeDeriv, normdSS); |
| 545 | gageParmSet(ctx, gageParmStackNormalizeDerivBias, biasSS); |
| 546 | if (!E) E |= !(pvlSS = AIR_CAST(gagePerVolume **,((gagePerVolume **)(calloc(sbp->num, sizeof(gagePerVolume * )))) |
| 547 | calloc(sbp->num, sizeof(gagePerVolume *)))((gagePerVolume **)(calloc(sbp->num, sizeof(gagePerVolume * ))))); |
| 548 | if (!E) airMopAdd(mop, pvlSS, (airMopper)airFree, airMopAlways); |
| 549 | if (!E) E |= gageStackPerVolumeNew(ctx, pvlSS, |
| 550 | AIR_CAST(const Nrrd*const*, ninSS)((const Nrrd*const*)(ninSS)), |
| 551 | sbp->num, kind); |
| 552 | if (!E) E |= gageStackPerVolumeAttach(ctx, pvl, pvlSS, |
| 553 | sbp->sigma, sbp->num); |
| 554 | if (!E) E |= gageKernelSet(ctx, gageKernelStack, kSS->kernel, kSS->parm); |
| 555 | } else { |
| 556 | if (!E) E |= gagePerVolumeAttach(ctx, pvl); |
| 557 | } |
| 558 | if (!E) E |= gageQueryItemOn(ctx, pvl, what); |
| 559 | if (!E) E |= gageUpdate(ctx); |
| 560 | if (E) { |
| 561 | airMopAdd(mop, err = biffGetDone(GAGEgageBiffKey), airFree, airMopAlways); |
| 562 | fprintf(stderr__stderrp, "%s: trouble:\n%s\n", me, err); |
| 563 | airMopError(mop); return 1; |
| 564 | } |
| 565 | answer = gageAnswerPointer(ctx, pvl, what); |
| 566 | ansLen = kind->table[what].answerLength; |
| 567 | /*** |
| 568 | **** end gage setup. |
| 569 | ***/ |
| 570 | if (verbose) { |
| 571 | fprintf(stderr__stderrp, "%s: kernel support = %d^3 samples\n", me, |
| 572 | 2*ctx->radius); |
| 573 | } |
| 574 | |
| 575 | if (ELL_3V_EXISTS(pntPos)((((int)(!(((pntPos)[0]) - ((pntPos)[0]))))) && (((int )(!(((pntPos)[1]) - ((pntPos)[1]))))) && (((int)(!((( pntPos)[2]) - ((pntPos)[2]))))))) { |
| 576 | /* only interested in a single point, make sure we have the right |
| 577 | info about the point WRT scale stuff */ |
| 578 | if (sbp) { |
| 579 | if (!(4 == pntPosNum && ELL_4V_EXISTS(pntPos)((((int)(!(((pntPos)[0]) - ((pntPos)[0]))))) && (((int )(!(((pntPos)[1]) - ((pntPos)[1]))))) && (((int)(!((( pntPos)[2]) - ((pntPos)[2]))))) && (((int)(!(((pntPos )[3]) - ((pntPos)[3])))))))) { |
| 580 | fprintf(stderr__stderrp, "%s: need a 4-vec position with scale-space", me); |
| 581 | airMopError(mop); return 1; |
| 582 | } |
| 583 | } else { |
| 584 | if (!(3 == pntPosNum)) { |
| 585 | fprintf(stderr__stderrp, "%s: need a 3-vec position (w/out scale-space)", me); |
| 586 | airMopError(mop); return 1; |
| 587 | } |
| 588 | } |
| 589 | if (sbp |
| 590 | ? gageStackProbeSpace(ctx, |
| 591 | pntPos[0], pntPos[1], pntPos[2], pntPos[3], |
| 592 | probeSpaceIndex, clamp) |
| 593 | : gageProbeSpace(ctx, |
| 594 | pntPos[0], pntPos[1], pntPos[2], |
| 595 | probeSpaceIndex, clamp)) { |
| 596 | fprintf(stderr__stderrp, "%s: trouble probing: (errNum %d) %s\n", me, |
| 597 | ctx->errNum, ctx->errStr); |
| 598 | airMopError(mop); return 1; |
| 599 | } |
| 600 | if (sbp) { |
| 601 | printf("%s: %s(%s:%g,%g,%g,%g) = ", me, airEnumStr(kind->enm, what), |
| 602 | probeSpaceIndex ? "index" : "world", |
| 603 | pntPos[0], pntPos[1], pntPos[2], pntPos[3]); |
| 604 | } else { |
| 605 | printf("%s: %s(%s:%g,%g,%g) = ", me, airEnumStr(kind->enm, what), |
| 606 | probeSpaceIndex ? "index" : "world", |
| 607 | pntPos[0], pntPos[1], pntPos[2]); |
| 608 | } |
| 609 | printans(stdout__stdoutp, answer, ansLen); |
| 610 | printf("\n"); |
| 611 | /* we're done, get out of here */ |
| 612 | /* except if we're supposed to debug derivatives */ |
| 613 | if (eps && 1 == ansLen) { |
| 614 | double v[3][3][3], fes, ee; |
| 615 | int xo, yo, zo; |
| 616 | if (probeSpaceIndex) { |
| 617 | fprintf(stderr__stderrp, "\n%s: WARNING!!: not probing in world-space (via " |
| 618 | "\"-wsp\") likely leads to errors in estimated " |
| 619 | "derivatives\n\n", me); |
| 620 | } |
| 621 | gageParmSet(ctx, gageParmVerbose, 0); |
| 622 | #define PROBE(x, y, z)((sbp ? gageStackProbeSpace(ctx, x, y, z, posSS, probeSpaceIndex , clamp) : gageProbeSpace(ctx, x, y, z, probeSpaceIndex, clamp )),answer[0]) \ |
| 623 | ((sbp \ |
| 624 | ? gageStackProbeSpace(ctx, x, y, z, posSS, \ |
| 625 | probeSpaceIndex, clamp) \ |
| 626 | : gageProbeSpace(ctx, x, y, z, probeSpaceIndex, \ |
| 627 | clamp)),answer[0]) |
| 628 | for (xo=0; xo<=2; xo++) { |
| 629 | for (yo=0; yo<=2; yo++) { |
| 630 | for (zo=0; zo<=2; zo++) { |
| 631 | v[xo][yo][zo] = PROBE(pntPos[0] + (xo-1)*eps,((sbp ? gageStackProbeSpace(ctx, pntPos[0] + (xo-1)*eps, pntPos [1] + (yo-1)*eps, pntPos[2] + (zo-1)*eps, posSS, probeSpaceIndex , clamp) : gageProbeSpace(ctx, pntPos[0] + (xo-1)*eps, pntPos [1] + (yo-1)*eps, pntPos[2] + (zo-1)*eps, probeSpaceIndex, clamp )),answer[0]) |
| 632 | pntPos[1] + (yo-1)*eps,((sbp ? gageStackProbeSpace(ctx, pntPos[0] + (xo-1)*eps, pntPos [1] + (yo-1)*eps, pntPos[2] + (zo-1)*eps, posSS, probeSpaceIndex , clamp) : gageProbeSpace(ctx, pntPos[0] + (xo-1)*eps, pntPos [1] + (yo-1)*eps, pntPos[2] + (zo-1)*eps, probeSpaceIndex, clamp )),answer[0]) |
| 633 | pntPos[2] + (zo-1)*eps)((sbp ? gageStackProbeSpace(ctx, pntPos[0] + (xo-1)*eps, pntPos [1] + (yo-1)*eps, pntPos[2] + (zo-1)*eps, posSS, probeSpaceIndex , clamp) : gageProbeSpace(ctx, pntPos[0] + (xo-1)*eps, pntPos [1] + (yo-1)*eps, pntPos[2] + (zo-1)*eps, probeSpaceIndex, clamp )),answer[0]); |
| 634 | } |
| 635 | } |
| 636 | } |
| 637 | printf("%s: approx gradient(%s) at (%g,%g,%g) = %f %f %f\n", me, |
| 638 | airEnumStr(kind->enm, what), pntPos[0], pntPos[1], pntPos[2], |
| 639 | (v[2][1][1] - v[0][1][1])/(2*eps), |
| 640 | (v[1][2][1] - v[1][0][1])/(2*eps), |
| 641 | (v[1][1][2] - v[1][1][0])/(2*eps)); |
| 642 | fes = 4*eps*eps; |
| 643 | ee = eps*eps; |
| 644 | printf("%s: approx hessian(%s) at (%g,%g,%g) = \n" |
| 645 | "%f %f %f\n" |
| 646 | " %f %f\n" |
| 647 | " %f\n", me, |
| 648 | airEnumStr(kind->enm, what), pntPos[0], pntPos[1], pntPos[2], |
| 649 | (v[0][1][1] - 2*v[1][1][1] + v[2][1][1])/ee, |
| 650 | (v[2][2][1] - v[0][2][1] - v[2][0][1] + v[0][0][1])/fes, |
| 651 | (v[2][1][2] - v[0][1][2] - v[2][1][0] + v[0][1][0])/fes, |
| 652 | (v[1][2][1] - 2*v[1][1][1] + v[1][0][1])/ee, |
| 653 | (v[1][2][2] - v[1][0][2] - v[1][2][0] + v[1][0][0])/fes, |
| 654 | (v[1][1][2] - 2*v[1][1][1] + v[1][1][0])/ee); |
| 655 | } |
| 656 | airMopOkay(mop); return 0; |
| 657 | } |
| 658 | |
| 659 | if (_npos) { |
| 660 | /* given a nrrd of probe locations */ |
| 661 | double *pos, (*ins)(void *v, size_t I, double d); |
| 662 | size_t II, NN; |
| 663 | unsigned int aidx; |
| 664 | if (!(2 == _npos->dim |
| 665 | && (3 == _npos->axis[0].size || 4 == _npos->axis[0].size))) { |
| 666 | fprintf(stderr__stderrp, "%s: need npos 2-D 3-by-N or 4-by-N " |
| 667 | "(not %u-D %u-by-N)\n", me, _npos->dim, |
| 668 | AIR_UINT(_npos->axis[0].size)((unsigned int)(_npos->axis[0].size))); |
| 669 | airMopError(mop); return 1; |
| 670 | } |
| 671 | if ((sbp && 3 == _npos->axis[0].size) |
| 672 | || (!sbp && 4 == _npos->axis[0].size)) { |
| 673 | fprintf(stderr__stderrp, "%s: have %u point coords but %s using scale-space\n", |
| 674 | me, AIR_UINT(_npos->axis[0].size)((unsigned int)(_npos->axis[0].size)), |
| 675 | sbp ? "are" : "are not"); |
| 676 | airMopError(mop); return 1; |
| 677 | } |
| 678 | NN = _npos->axis[1].size; |
| 679 | npos = nrrdNew(); |
| 680 | airMopAdd(mop, npos, AIR_CAST(airMopper, nrrdNuke)((airMopper)(nrrdNuke)), airMopAlways); |
| 681 | nout = nrrdNew(); |
| 682 | airMopAdd(mop, nout, AIR_CAST(airMopper, nrrdNuke)((airMopper)(nrrdNuke)), airMopAlways); |
| 683 | if (nrrdConvert(npos, _npos, nrrdTypeDouble) |
| 684 | || nrrdMaybeAlloc_va(nout, otype, 2, |
| 685 | AIR_CAST(size_t, ansLen)((size_t)(ansLen)), |
| 686 | AIR_CAST(size_t, NN)((size_t)(NN)))) { |
| 687 | airMopAdd(mop, err = biffGetDone(NRRDnrrdBiffKey), airFree, airMopAlways); |
| 688 | fprintf(stderr__stderrp, "%s: trouble with npos or nout:\n%s\n", me, err); |
| 689 | airMopError(mop); return 1; |
| 690 | } |
| 691 | |
| 692 | pos = AIR_CAST(double *, npos->data)((double *)(npos->data)); |
| 693 | ins = nrrdDInsert[nout->type]; |
| 694 | for (II=0; II<NN; II++) { |
| 695 | if (sbp) { |
| 696 | gageStackProbeSpace(ctx, pos[0], pos[1], pos[2], pos[3], |
| 697 | probeSpaceIndex, clamp); |
| 698 | } else { |
| 699 | gageProbeSpace(ctx, pos[0], pos[1], pos[2], |
| 700 | probeSpaceIndex, clamp); |
| 701 | } |
| 702 | if (1 == ansLen) { |
| 703 | ins(nout->data, II, (ctx->edgeFrac > edgeFracInfo[0] |
| 704 | ? edgeFracInfo[1] |
| 705 | : *answer)); |
| 706 | } else { |
| 707 | for (aidx=0; aidx<ansLen; aidx++) { |
| 708 | ins(nout->data, aidx + ansLen*II, (ctx->edgeFrac > edgeFracInfo[0] |
| 709 | ? edgeFracInfo[1] |
| 710 | : answer[aidx])); |
| 711 | } |
| 712 | } |
| 713 | /* |
| 714 | if (sbp) { |
| 715 | printf("%s: %s(%s:%g,%g,%g,%g) = ", me, airEnumStr(kind->enm, what), |
| 716 | probeSpaceIndex ? "index" : "world", |
| 717 | pos[0], pos[1], pos[2], pos[3]); |
| 718 | } else { |
| 719 | printf("%s: %s(%s:%g,%g,%g) = ", me, airEnumStr(kind->enm, what), |
| 720 | probeSpaceIndex ? "index" : "world", |
| 721 | pos[0], pos[1], pos[2]); |
| 722 | } |
| 723 | printans(stdout, answer, ansLen); |
| 724 | printf("\n"); |
| 725 | */ |
| 726 | pos += _npos->axis[0].size; |
| 727 | } |
| 728 | if (nrrdSave(outS, nout, NULL((void*)0))) { |
| 729 | airMopAdd(mop, err = biffGetDone(NRRDnrrdBiffKey), airFree, airMopAlways); |
| 730 | fprintf(stderr__stderrp, "%s: trouble saving output:\n%s\n", me, err); |
| 731 | airMopError(mop); |
| 732 | return 1; |
| 733 | } |
| 734 | |
| 735 | /* we're done, get out of here */ |
| 736 | airMopOkay(mop); |
| 737 | exit(0); |
| 738 | } |
| 739 | |
| 740 | /* else, we're sampling on some kind of grid */ |
| 741 | ngrid = nrrdNew(); |
| 742 | airMopAdd(mop, ngrid, (airMopper)nrrdNuke, airMopAlways); |
| 743 | iBaseDim = kind->baseDim; |
| 744 | oBaseDim = 1 == ansLen ? 0 : 1; |
| 745 | if (!_ngrid) { |
| 746 | /* did not get a grid, have to use "-s" args to define it */ |
| 747 | double *grid; |
| 748 | size_t gridSize[NRRD_DIM_MAX16]; |
| 749 | six = nin->axis[0+iBaseDim].size; |
| 750 | siy = nin->axis[1+iBaseDim].size; |
| 751 | siz = nin->axis[2+iBaseDim].size; |
| 752 | dsix = AIR_CAST(double, six)((double)(six)); |
| 753 | dsiy = AIR_CAST(double, siy)((double)(siy)); |
| 754 | dsiz = AIR_CAST(double, siz)((double)(siz)); |
| 755 | sox = AIR_CAST(size_t, scale[0]*dsix)((size_t)(scale[0]*dsix)); |
| 756 | soy = AIR_CAST(size_t, scale[1]*dsiy)((size_t)(scale[1]*dsiy)); |
| 757 | soz = AIR_CAST(size_t, scale[2]*dsiz)((size_t)(scale[2]*dsiz)); |
| 758 | dsox = AIR_CAST(double, sox)((double)(sox)); |
| 759 | dsoy = AIR_CAST(double, soy)((double)(soy)); |
| 760 | dsoz = AIR_CAST(double, soz)((double)(soz)); |
| 761 | rscl[0] = dsix/dsox; |
| 762 | rscl[1] = dsiy/dsoy; |
| 763 | rscl[2] = dsiz/dsoz; |
| 764 | if (verbose) { |
| 765 | fprintf(stderr__stderrp, "%s: creating %u x %s x %s x %s output\n", me, |
| 766 | ansLen, |
| 767 | airSprintSize_t(stmp[1], sox), |
| 768 | airSprintSize_t(stmp[2], soy), |
| 769 | airSprintSize_t(stmp[3], soz)); |
| 770 | fprintf(stderr__stderrp, "%s: effective scaling is %g %g %g\n", me, |
| 771 | rscl[0], rscl[1], rscl[2]); |
| 772 | } |
| 773 | gridSize[0] = sbp ? 5 : 4; |
| 774 | gridSize[1] = 4; |
| 775 | if (nrrdMaybeAlloc_nva(ngrid, nrrdTypeDouble, 2, gridSize)) { |
| 776 | airMopAdd(mop, err = biffGetDone(NRRDnrrdBiffKey), airFree, airMopAlways); |
| 777 | fprintf(stderr__stderrp, "%s: trouble making ngrid:\n%s\n", me, err); |
| 778 | airMopError(mop); return 1; |
| 779 | } |
| 780 | grid = AIR_CAST(double *, ngrid->data)((double *)(ngrid->data)); |
| 781 | if (nrrdCenterCell == ctx->shape->center) { |
| 782 | ELL_3V_SET(min, -0.5, -0.5, -0.5)((min)[0] = (-0.5), (min)[1] = (-0.5), (min)[2] = (-0.5)); |
| 783 | ELL_3V_SET(maxOut, dsox-0.5, dsoy-0.5, dsoz-0.5)((maxOut)[0] = (dsox-0.5), (maxOut)[1] = (dsoy-0.5), (maxOut) [2] = (dsoz-0.5)); |
| 784 | ELL_3V_SET(maxIn, dsix-0.5, dsiy-0.5, dsiz-0.5)((maxIn)[0] = (dsix-0.5), (maxIn)[1] = (dsiy-0.5), (maxIn)[2] = (dsiz-0.5)); |
| 785 | } else { |
| 786 | ELL_3V_SET(min, 0, 0, 0)((min)[0] = (0), (min)[1] = (0), (min)[2] = (0)); |
| 787 | ELL_3V_SET(maxOut, dsox-1, dsoy-1, dsoz-1)((maxOut)[0] = (dsox-1), (maxOut)[1] = (dsoy-1), (maxOut)[2] = (dsoz-1)); |
| 788 | ELL_3V_SET(maxIn, dsix-1, dsiy-1, dsiz-1)((maxIn)[0] = (dsix-1), (maxIn)[1] = (dsiy-1), (maxIn)[2] = ( dsiz-1)); |
| 789 | } |
| 790 | ELL_4V_SET(grid + gridSize[0]*0, 3,((grid + gridSize[0]*0)[0] = (3), (grid + gridSize[0]*0)[1] = ((nrrdCenterCell == (ctx->shape->center) ? ( ((double) (((maxIn[0])))-(((min[0]))))*((double)(((0)) + 0.5)-(0)) / (( double)(((sox)))-(0)) + (((min[0])))) : ( ((double)(((maxIn[0 ])))-(((min[0]))))*((double)(((0)))-(0)) / ((double)(((sox))- 1)-(0)) + (((min[0])))))), (grid + gridSize[0]*0)[2] = ((nrrdCenterCell == (ctx->shape->center) ? ( ((double)(((maxIn[1])))-(( (min[1]))))*((double)(((0)) + 0.5)-(0)) / ((double)(((soy)))- (0)) + (((min[1])))) : ( ((double)(((maxIn[1])))-(((min[1]))) )*((double)(((0)))-(0)) / ((double)(((soy))-1)-(0)) + (((min[ 1])))))), (grid + gridSize[0]*0)[3] = ((nrrdCenterCell == (ctx ->shape->center) ? ( ((double)(((maxIn[2])))-(((min[2]) )))*((double)(((0)) + 0.5)-(0)) / ((double)(((soz)))-(0)) + ( ((min[2])))) : ( ((double)(((maxIn[2])))-(((min[2]))))*((double )(((0)))-(0)) / ((double)(((soz))-1)-(0)) + (((min[2]))))))) |
| 791 | NRRD_POS(ctx->shape->center, min[0], maxIn[0], sox, 0),((grid + gridSize[0]*0)[0] = (3), (grid + gridSize[0]*0)[1] = ((nrrdCenterCell == (ctx->shape->center) ? ( ((double) (((maxIn[0])))-(((min[0]))))*((double)(((0)) + 0.5)-(0)) / (( double)(((sox)))-(0)) + (((min[0])))) : ( ((double)(((maxIn[0 ])))-(((min[0]))))*((double)(((0)))-(0)) / ((double)(((sox))- 1)-(0)) + (((min[0])))))), (grid + gridSize[0]*0)[2] = ((nrrdCenterCell == (ctx->shape->center) ? ( ((double)(((maxIn[1])))-(( (min[1]))))*((double)(((0)) + 0.5)-(0)) / ((double)(((soy)))- (0)) + (((min[1])))) : ( ((double)(((maxIn[1])))-(((min[1]))) )*((double)(((0)))-(0)) / ((double)(((soy))-1)-(0)) + (((min[ 1])))))), (grid + gridSize[0]*0)[3] = ((nrrdCenterCell == (ctx ->shape->center) ? ( ((double)(((maxIn[2])))-(((min[2]) )))*((double)(((0)) + 0.5)-(0)) / ((double)(((soz)))-(0)) + ( ((min[2])))) : ( ((double)(((maxIn[2])))-(((min[2]))))*((double )(((0)))-(0)) / ((double)(((soz))-1)-(0)) + (((min[2]))))))) |
| 792 | NRRD_POS(ctx->shape->center, min[1], maxIn[1], soy, 0),((grid + gridSize[0]*0)[0] = (3), (grid + gridSize[0]*0)[1] = ((nrrdCenterCell == (ctx->shape->center) ? ( ((double) (((maxIn[0])))-(((min[0]))))*((double)(((0)) + 0.5)-(0)) / (( double)(((sox)))-(0)) + (((min[0])))) : ( ((double)(((maxIn[0 ])))-(((min[0]))))*((double)(((0)))-(0)) / ((double)(((sox))- 1)-(0)) + (((min[0])))))), (grid + gridSize[0]*0)[2] = ((nrrdCenterCell == (ctx->shape->center) ? ( ((double)(((maxIn[1])))-(( (min[1]))))*((double)(((0)) + 0.5)-(0)) / ((double)(((soy)))- (0)) + (((min[1])))) : ( ((double)(((maxIn[1])))-(((min[1]))) )*((double)(((0)))-(0)) / ((double)(((soy))-1)-(0)) + (((min[ 1])))))), (grid + gridSize[0]*0)[3] = ((nrrdCenterCell == (ctx ->shape->center) ? ( ((double)(((maxIn[2])))-(((min[2]) )))*((double)(((0)) + 0.5)-(0)) / ((double)(((soz)))-(0)) + ( ((min[2])))) : ( ((double)(((maxIn[2])))-(((min[2]))))*((double )(((0)))-(0)) / ((double)(((soz))-1)-(0)) + (((min[2]))))))) |
| 793 | NRRD_POS(ctx->shape->center, min[2], maxIn[2], soz, 0))((grid + gridSize[0]*0)[0] = (3), (grid + gridSize[0]*0)[1] = ((nrrdCenterCell == (ctx->shape->center) ? ( ((double) (((maxIn[0])))-(((min[0]))))*((double)(((0)) + 0.5)-(0)) / (( double)(((sox)))-(0)) + (((min[0])))) : ( ((double)(((maxIn[0 ])))-(((min[0]))))*((double)(((0)))-(0)) / ((double)(((sox))- 1)-(0)) + (((min[0])))))), (grid + gridSize[0]*0)[2] = ((nrrdCenterCell == (ctx->shape->center) ? ( ((double)(((maxIn[1])))-(( (min[1]))))*((double)(((0)) + 0.5)-(0)) / ((double)(((soy)))- (0)) + (((min[1])))) : ( ((double)(((maxIn[1])))-(((min[1]))) )*((double)(((0)))-(0)) / ((double)(((soy))-1)-(0)) + (((min[ 1])))))), (grid + gridSize[0]*0)[3] = ((nrrdCenterCell == (ctx ->shape->center) ? ( ((double)(((maxIn[2])))-(((min[2]) )))*((double)(((0)) + 0.5)-(0)) / ((double)(((soz)))-(0)) + ( ((min[2])))) : ( ((double)(((maxIn[2])))-(((min[2]))))*((double )(((0)))-(0)) / ((double)(((soz))-1)-(0)) + (((min[2]))))))); |
| 794 | ELL_4V_SET(grid + gridSize[0]*1, dsox,((grid + gridSize[0]*1)[0] = (dsox), (grid + gridSize[0]*1)[1 ] = (( ((double)(maxIn[0])-(min[0]))*((double)(1)) / ((double )(maxOut[0])-(min[0])) )), (grid + gridSize[0]*1)[2] = (0), ( grid + gridSize[0]*1)[3] = (0)) |
| 795 | AIR_DELTA(min[0], 1, maxOut[0], min[0], maxIn[0]),((grid + gridSize[0]*1)[0] = (dsox), (grid + gridSize[0]*1)[1 ] = (( ((double)(maxIn[0])-(min[0]))*((double)(1)) / ((double )(maxOut[0])-(min[0])) )), (grid + gridSize[0]*1)[2] = (0), ( grid + gridSize[0]*1)[3] = (0)) |
| 796 | 0,((grid + gridSize[0]*1)[0] = (dsox), (grid + gridSize[0]*1)[1 ] = (( ((double)(maxIn[0])-(min[0]))*((double)(1)) / ((double )(maxOut[0])-(min[0])) )), (grid + gridSize[0]*1)[2] = (0), ( grid + gridSize[0]*1)[3] = (0)) |
| 797 | 0)((grid + gridSize[0]*1)[0] = (dsox), (grid + gridSize[0]*1)[1 ] = (( ((double)(maxIn[0])-(min[0]))*((double)(1)) / ((double )(maxOut[0])-(min[0])) )), (grid + gridSize[0]*1)[2] = (0), ( grid + gridSize[0]*1)[3] = (0)); |
| 798 | ELL_4V_SET(grid + gridSize[0]*2, dsoy,((grid + gridSize[0]*2)[0] = (dsoy), (grid + gridSize[0]*2)[1 ] = (0), (grid + gridSize[0]*2)[2] = (( ((double)(maxIn[1])-( min[1]))*((double)(1)) / ((double)(maxOut[1])-(min[1])) )), ( grid + gridSize[0]*2)[3] = (0)) |
| 799 | 0,((grid + gridSize[0]*2)[0] = (dsoy), (grid + gridSize[0]*2)[1 ] = (0), (grid + gridSize[0]*2)[2] = (( ((double)(maxIn[1])-( min[1]))*((double)(1)) / ((double)(maxOut[1])-(min[1])) )), ( grid + gridSize[0]*2)[3] = (0)) |
| 800 | AIR_DELTA(min[1], 1, maxOut[1], min[1], maxIn[1]),((grid + gridSize[0]*2)[0] = (dsoy), (grid + gridSize[0]*2)[1 ] = (0), (grid + gridSize[0]*2)[2] = (( ((double)(maxIn[1])-( min[1]))*((double)(1)) / ((double)(maxOut[1])-(min[1])) )), ( grid + gridSize[0]*2)[3] = (0)) |
| 801 | 0)((grid + gridSize[0]*2)[0] = (dsoy), (grid + gridSize[0]*2)[1 ] = (0), (grid + gridSize[0]*2)[2] = (( ((double)(maxIn[1])-( min[1]))*((double)(1)) / ((double)(maxOut[1])-(min[1])) )), ( grid + gridSize[0]*2)[3] = (0)); |
| 802 | ELL_4V_SET(grid + gridSize[0]*3, dsoz,((grid + gridSize[0]*3)[0] = (dsoz), (grid + gridSize[0]*3)[1 ] = (0), (grid + gridSize[0]*3)[2] = (0), (grid + gridSize[0] *3)[3] = (( ((double)(maxIn[2])-(min[2]))*((double)(1)) / ((double )(maxOut[2])-(min[2])) ))) |
| 803 | 0,((grid + gridSize[0]*3)[0] = (dsoz), (grid + gridSize[0]*3)[1 ] = (0), (grid + gridSize[0]*3)[2] = (0), (grid + gridSize[0] *3)[3] = (( ((double)(maxIn[2])-(min[2]))*((double)(1)) / ((double )(maxOut[2])-(min[2])) ))) |
| 804 | 0,((grid + gridSize[0]*3)[0] = (dsoz), (grid + gridSize[0]*3)[1 ] = (0), (grid + gridSize[0]*3)[2] = (0), (grid + gridSize[0] *3)[3] = (( ((double)(maxIn[2])-(min[2]))*((double)(1)) / ((double )(maxOut[2])-(min[2])) ))) |
| 805 | AIR_DELTA(min[2], 1, maxOut[2], min[2], maxIn[2]))((grid + gridSize[0]*3)[0] = (dsoz), (grid + gridSize[0]*3)[1 ] = (0), (grid + gridSize[0]*3)[2] = (0), (grid + gridSize[0] *3)[3] = (( ((double)(maxIn[2])-(min[2]))*((double)(1)) / ((double )(maxOut[2])-(min[2])) ))); |
| 806 | if (sbp) { |
| 807 | if (!probeSpaceIndex) { |
| 808 | double idxSS = AIR_NAN(airFloatQNaN.f); |
| 809 | unsigned int vi; |
| 810 | /* there's actually work to do here, weirdly: gageProbe can |
| 811 | either work in index space, or in world space, but the |
| 812 | vprobe-style sampling is index-space-centric, so we have to |
| 813 | convert the world-space stack position to index space, to be |
| 814 | consistent with the way that the grid sampling will be |
| 815 | defined. So, we have to actually replicate work that is done |
| 816 | by _gageProbeSpace() in converting from world to index space */ |
| 817 | /* HEY: the way that idxSS is set is very strange */ |
| 818 | for (vi=0; vi<sbp->num-1; vi++) { |
| 819 | if (AIR_IN_CL(sbp->sigma[vi], posSS, sbp->sigma[vi+1])((sbp->sigma[vi]) <= (posSS) && (posSS) <= ( sbp->sigma[vi+1]))) { |
| 820 | idxSS = vi + AIR_AFFINE(sbp->sigma[vi], posSS, sbp->sigma[vi+1],( ((double)(1)-(0))*((double)(posSS)-(sbp->sigma[vi])) / ( (double)(sbp->sigma[vi+1])-(sbp->sigma[vi])) + (0)) |
| 821 | 0, 1)( ((double)(1)-(0))*((double)(posSS)-(sbp->sigma[vi])) / ( (double)(sbp->sigma[vi+1])-(sbp->sigma[vi])) + (0)); |
| 822 | if (verbose > 1) { |
| 823 | fprintf(stderr__stderrp, "%s: scale pos %g -> idx %g = %u + %g\n", me, |
| 824 | posSS, idxSS, vi, |
| 825 | AIR_AFFINE(sbp->sigma[vi], posSS, sbp->sigma[vi+1],( ((double)(1)-(0))*((double)(posSS)-(sbp->sigma[vi])) / ( (double)(sbp->sigma[vi+1])-(sbp->sigma[vi])) + (0)) |
| 826 | 0, 1)( ((double)(1)-(0))*((double)(posSS)-(sbp->sigma[vi])) / ( (double)(sbp->sigma[vi+1])-(sbp->sigma[vi])) + (0))); |
| 827 | } |
| 828 | break; |
| 829 | } |
| 830 | } |
| 831 | if (vi == sbp->num-1) { |
| 832 | fprintf(stderr__stderrp, "%s: scale pos %g outside range %g=%g, %g=%g\n", me, |
| 833 | posSS, rangeSS[0], sbp->sigma[0], |
| 834 | rangeSS[1], sbp->sigma[sbp->num-1]); |
| 835 | airMopError(mop); return 1; |
| 836 | } |
| 837 | grid[4 + 5*0] = idxSS; |
| 838 | } else { |
| 839 | grid[4 + 5*0] = posSS; |
| 840 | } |
| 841 | grid[4 + 5*1] = 0; |
| 842 | grid[4 + 5*2] = 0; |
| 843 | grid[4 + 5*3] = 0; |
| 844 | } |
| 845 | } else { |
| 846 | /* we did get a grid, here we only copy from _ngrid to ngrid, |
| 847 | and let further massaging be done in gridProbe below */ |
| 848 | six = siy = siz = 0; |
| 849 | sox = soy = soz = 0; |
Value stored to 'sox' is never read | |
| 850 | if (nrrdCopy(ngrid, _ngrid)) { |
| 851 | airMopAdd(mop, err = biffGetDone(NRRDnrrdBiffKey), airFree, airMopAlways); |
| 852 | fprintf(stderr__stderrp, "%s: trouble copying ngrid:\n%s\n", me, err); |
| 853 | airMopError(mop); return 1; |
| 854 | } |
| 855 | } |
| 856 | |
| 857 | /* probe onto grid */ |
| 858 | nout = nrrdNew(); |
| 859 | airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways); |
| 860 | gageParmSet(ctx, gageParmVerbose, verbose); |
| 861 | t0 = airTime(); |
| 862 | if (gridProbe(ctx, pvl, what, nout, otype, ngrid, |
| 863 | (_ngrid |
| 864 | ? probeSpaceIndex /* user specifies grid space */ |
| 865 | : AIR_TRUE1), /* copying vprobe index-space behavior */ |
| 866 | verbose, clamp, scaleIsTau, |
| 867 | edgeFracInfo[0], edgeFracInfo[1])) { |
| 868 | /* note hijacking of GAGE key */ |
| 869 | airMopAdd(mop, err = biffGetDone(GAGEgageBiffKey), airFree, airMopAlways); |
| 870 | fprintf(stderr__stderrp, "%s: trouble probing on grid:\n%s\n", me, err); |
| 871 | airMopError(mop); return 1; |
| 872 | } |
| 873 | t1 = airTime(); |
| 874 | if (verbose) { |
| 875 | fprintf(stderr__stderrp, "probe rate = %g KHz\n", |
| 876 | AIR_CAST(double, nrrdElementNumber(nout)/ansLen)((double)(nrrdElementNumber(nout)/ansLen)) |
| 877 | / (1000.0*(t1-t0))); |
| 878 | } |
| 879 | |
| 880 | /* massage output some */ |
| 881 | nrrdContentSet_va(nout, "gprobe", nin, "%s", airEnumStr(kind->enm, what)); |
| 882 | if (!_ngrid) { |
| 883 | /* did not get a grid, have to emulate vprobe per-axis behavior */ |
| 884 | for (axi=0; axi<3; axi++) { |
| 885 | nout->axis[axi+oBaseDim].label = |
| 886 | airStrdup(nin->axis[axi+iBaseDim].label); |
| 887 | nout->axis[axi+oBaseDim].center = ctx->shape->center; |
| 888 | } |
| 889 | nrrdBasicInfoCopy(nout, nin, (NRRD_BASIC_INFO_DATA_BIT(1<< 1) |
| 890 | | NRRD_BASIC_INFO_TYPE_BIT(1<< 2) |
| 891 | | NRRD_BASIC_INFO_BLOCKSIZE_BIT(1<< 3) |
| 892 | | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4) |
| 893 | | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5) |
| 894 | | NRRD_BASIC_INFO_SPACEORIGIN_BIT(1<<10) |
| 895 | | NRRD_BASIC_INFO_OLDMIN_BIT(1<<12) |
| 896 | | NRRD_BASIC_INFO_OLDMAX_BIT(1<<13) |
| 897 | | NRRD_BASIC_INFO_COMMENTS_BIT(1<<14) |
| 898 | | NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT(1<<15))); |
| 899 | if (ctx->shape->fromOrientation) { |
| 900 | if (nin->space) { |
| 901 | nrrdSpaceSet(nout, nin->space); |
| 902 | } else { |
| 903 | nrrdSpaceDimensionSet(nout, nin->spaceDim); |
| 904 | } |
| 905 | nrrdSpaceVecCopy(nout->spaceOrigin, nin->spaceOrigin); |
| 906 | for (axi=0; axi<3; axi++) { |
| 907 | double tmp; |
| 908 | nrrdSpaceVecScale(nout->axis[axi+oBaseDim].spaceDirection, |
| 909 | rscl[axi], |
| 910 | nin->axis[axi+iBaseDim].spaceDirection); |
| 911 | tmp = AIR_AFFINE(min[axi], 0, maxOut[axi], min[axi], maxIn[axi])( ((double)(maxIn[axi])-(min[axi]))*((double)(0)-(min[axi])) / ((double)(maxOut[axi])-(min[axi])) + (min[axi])); |
| 912 | nrrdSpaceVecScaleAdd2(nout->spaceOrigin, |
| 913 | 1.0, nout->spaceOrigin, |
| 914 | tmp, nin->axis[axi+iBaseDim].spaceDirection); |
| 915 | } |
| 916 | } else { |
| 917 | for (axi=0; axi<3; axi++) { |
| 918 | nout->axis[axi+oBaseDim].spacing = |
| 919 | rscl[axi]*nin->axis[axi+iBaseDim].spacing; |
| 920 | } |
| 921 | } |
| 922 | } |
| 923 | |
| 924 | if (nrrdSave(outS, nout, NULL((void*)0))) { |
| 925 | airMopAdd(mop, err = biffGetDone(NRRDnrrdBiffKey), airFree, airMopAlways); |
| 926 | fprintf(stderr__stderrp, "%s: trouble saving output:\n%s\n", me, err); |
| 927 | airMopError(mop); |
| 928 | return 1; |
| 929 | } |
| 930 | |
| 931 | airMopOkay(mop); |
| 932 | return 0; |
| 933 | } |