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); |
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_MAX], coordOut[NRRD_DIM_MAX], II, NN; |
59 |
|
|
double (*ins)(void *v, size_t I, double d); |
60 |
|
|
char stmp[2][AIR_STRLEN_SMALL]; |
61 |
|
|
|
62 |
|
|
if (!(ctx && pvl && nout && _ngrid)) { |
63 |
|
|
biffAddf(GAGE, "%s: got NULL pointer", me); |
64 |
|
|
return 1; |
65 |
|
|
} |
66 |
|
|
if (airEnumValCheck(nrrdType, typeOut)) { |
67 |
|
|
biffAddf(GAGE, "%s: type %d not valid", me, typeOut); |
68 |
|
|
return 1; |
69 |
|
|
} |
70 |
|
|
if (!gagePerVolumeIsAttached(ctx, pvl)) { |
71 |
|
|
biffAddf(GAGE, "%s: given pvl not attached to context", me); |
72 |
|
|
return 1; |
73 |
|
|
} |
74 |
|
|
if (!(2 == _ngrid->dim)) { |
75 |
|
|
biffAddf(GAGE, "%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(GAGE, "%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)); |
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(GAGE, NRRD, "%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); /* 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(GAGE, NRRD, "%s: trouble converting/padding ngrid", me); |
108 |
|
|
airMopError(mop); return 1; |
109 |
|
|
} |
110 |
|
|
} |
111 |
|
|
grid = AIR_CAST(double *, ngrid->data); |
112 |
|
|
gridDim = AIR_ROUNDUP_UI(grid[0]); |
113 |
|
|
if (gridDim + 1 != ngrid->axis[1].size) { |
114 |
|
|
biffAddf(GAGE, "%s: ngrid->axis[1].size = %u but expected %u = 1 + %u", |
115 |
|
|
me, AIR_UINT(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_MAX) { |
124 |
|
|
biffAddf(GAGE, "%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)]); |
134 |
|
|
NN *= sizeOut[aidx + baseDim]; |
135 |
|
|
coordOut[aidx + baseDim] = 0; |
136 |
|
|
} |
137 |
|
|
if (nrrdMaybeAlloc_nva(nout, typeOut, dim, sizeOut)) { |
138 |
|
|
biffMovef(GAGE, NRRD, "%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, "z = "); |
148 |
|
|
} |
149 |
|
|
fprintf(stderr, " %s/%s", |
150 |
|
|
airSprintSize_t(stmp[0], coordOut[2+baseDim]), |
151 |
|
|
airSprintSize_t(stmp[1], sizeOut[2+baseDim])); |
152 |
|
|
fflush(stderr); |
153 |
|
|
if (verbose > 1) { |
154 |
|
|
fprintf(stderr, "\n"); |
155 |
|
|
} |
156 |
|
|
} |
157 |
|
|
ELL_4V_COPY(pos, grid + 1 + 5*0); |
158 |
|
|
for (aidx=0; aidx<gridDim; aidx++) { |
159 |
|
|
ELL_4V_SCALE_ADD2(pos, 1, pos, |
160 |
|
|
AIR_CAST(double, coordOut[aidx + baseDim]), |
161 |
|
|
grid + 1 + 5*(1+aidx)); |
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(GAGE, "%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); |
199 |
|
|
} |
200 |
|
|
if (verbose && verbose <= 1) { |
201 |
|
|
fprintf(stderr, "\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; |
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; |
240 |
|
|
Nrrd *ngrad=NULL, *nbmat=NULL; |
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; |
246 |
|
|
double t0, t1, rscl[3], min[3], maxOut[3], maxIn[3]; |
247 |
|
|
airArray *mop; |
248 |
|
|
#define NON_SBP_OPT_NUM 5 |
249 |
|
|
unsigned int ansLen, *skip, skipNum, pntPosNum, |
250 |
|
|
nonSbpOpi[NON_SBP_OPT_NUM], nsi; |
251 |
|
|
gageStackBlurParm *sbpIN, *sbpCL, *sbp; |
252 |
|
|
int otype, clamp, scaleIsTau; |
253 |
|
|
char stmp[4][AIR_STRLEN_SMALL]; |
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, "******************************************\n"); |
263 |
|
|
fprintf(stderr, "******************************************\n"); |
264 |
|
|
fprintf(stderr, "\n"); |
265 |
|
|
fprintf(stderr, " %s: nrrd sanity check FAILED.\n", me); |
266 |
|
|
fprintf(stderr, "\n"); |
267 |
|
|
fprintf(stderr, " This means that either nrrd can't work on this " |
268 |
|
|
"platform, or (more likely)\n"); |
269 |
|
|
fprintf(stderr, " there was an error in the compilation options " |
270 |
|
|
"and variable definitions\n"); |
271 |
|
|
fprintf(stderr, " for how Teem was built here.\n"); |
272 |
|
|
fprintf(stderr, "\n"); |
273 |
|
|
fprintf(stderr, " %s\n", err = biffGetDone(NRRD)); |
274 |
|
|
fprintf(stderr, "\n"); |
275 |
|
|
fprintf(stderr, "******************************************\n"); |
276 |
|
|
fprintf(stderr, "******************************************\n"); |
277 |
|
|
free(err); |
278 |
|
|
return 1; |
279 |
|
|
} |
280 |
|
|
|
281 |
|
|
mop = airMopNew(); |
282 |
|
|
hparm = hestParmNew(); |
283 |
|
|
airMopAdd(mop, hparm, AIR_CAST(airMopper, hestParmFree), airMopAlways); |
284 |
|
|
hparm->elideSingleOtherType = AIR_TRUE; |
285 |
|
|
hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &nin, NULL, |
286 |
|
|
"input volume", NULL, NULL, nrrdHestNrrd); |
287 |
|
|
hestOptAdd(&hopt, "k", "kind", airTypeOther, 1, 1, &kind, NULL, |
288 |
|
|
"\"kind\" of volume (\"scalar\", \"vector\", " |
289 |
|
|
"\"tensor\", or \"dwi\")", |
290 |
|
|
NULL, NULL, meetHestGageKind); |
291 |
|
|
hestOptAdd(&hopt, "v", "verbosity", airTypeInt, 1, 1, &verbose, "1", |
292 |
|
|
"verbosity level"); |
293 |
|
|
hestOptAdd(&hopt, "q", "query", airTypeString, 1, 1, &whatS, NULL, |
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, "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, NULL, nrrdHestKernelSpec); |
321 |
|
|
hestOptAdd(&hopt, "k11", "kern11", airTypeOther, 1, 1, &k11, |
322 |
|
|
"cubicd:1,0", "kernel for gageKernel11", |
323 |
|
|
NULL, NULL, nrrdHestKernelSpec); |
324 |
|
|
hestOptAdd(&hopt, "k22", "kern22", airTypeOther, 1, 1, &k22, |
325 |
|
|
"cubicdd:1,0", "kernel for gageKernel22", |
326 |
|
|
NULL, NULL, nrrdHestKernelSpec); |
327 |
|
|
hestOptAdd(&hopt, "rn", NULL, airTypeInt, 0, 0, &renorm, NULL, |
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, airTypeInt, 0, 0, &uniformSS, NULL, |
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, airTypeInt, 0, 0, &optimSS, NULL, |
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, NULL, nrrdHestKernelSpec); |
352 |
|
|
if (nsi != NON_SBP_OPT_NUM) { |
353 |
|
|
fprintf(stderr, "%s: PANIC nsi %u != %u", me, nsi, NON_SBP_OPT_NUM); |
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, NULL, 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, airTypeInt, 0, 0, &normdSS, NULL, |
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, NULL, 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, NULL, 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, NULL, 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, 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_TRUE, AIR_TRUE, AIR_TRUE); |
417 |
|
|
airMopAdd(mop, hopt, AIR_CAST(airMopper, hestOptFree), airMopAlways); |
418 |
|
|
airMopAdd(mop, hopt, AIR_CAST(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, "%s: couldn't parse \"%s\" as measure of \"%s\" volume\n", |
424 |
|
|
me, whatS, kind->name); |
425 |
|
|
hestUsage(stderr, hopt, me, hparm); |
426 |
|
|
hestGlossary(stderr, 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, kind->name)) { |
432 |
|
|
if (tenDWMRIKeyValueParse(&ngrad, &nbmat, &bval, &skip, &skipNum, nin)) { |
433 |
|
|
airMopAdd(mop, err = biffGetDone(TEN), airFree, airMopAlways); |
434 |
|
|
fprintf(stderr, "%s: trouble parsing DWI info:\n%s\n", me, err); |
435 |
|
|
airMopError(mop); return 1; |
436 |
|
|
} |
437 |
|
|
if (skipNum) { |
438 |
|
|
fprintf(stderr, "%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(TEN), airFree, airMopAlways); |
446 |
|
|
fprintf(stderr, "%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_FALSE; |
459 |
|
|
for (nsi=0; nsi<NON_SBP_OPT_NUM; nsi++) { |
460 |
|
|
gotOld |= (hestSourceUser == hopt[nonSbpOpi[nsi]].source); |
461 |
|
|
} |
462 |
|
|
if (gotOld) { |
463 |
|
|
fprintf(stderr, "%s: with new -sbp option; can't also use older " |
464 |
|
|
"scale-space options (used", me); |
465 |
|
|
for (nsi=0; nsi<NON_SBP_OPT_NUM; nsi++) { |
466 |
|
|
if (hestSourceUser == hopt[nonSbpOpi[nsi]].source) { |
467 |
|
|
fprintf(stderr, " -%s", hopt[nonSbpOpi[nsi]].flag); |
468 |
|
|
} |
469 |
|
|
} |
470 |
|
|
fprintf(stderr, ")\n"); |
471 |
|
|
airMopError(mop); return 1; |
472 |
|
|
} |
473 |
|
|
if (gageStackBlurManage(&ninSS, &recompute, sbpCL, |
474 |
|
|
stackFnameFormat, AIR_TRUE, NULL, |
475 |
|
|
nin, kind)) { |
476 |
|
|
airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways); |
477 |
|
|
fprintf(stderr, "%s: trouble getting volume stack:\n%s\n", me, err); |
478 |
|
|
airMopError(mop); return 1; |
479 |
|
|
} |
480 |
|
|
if (sbpCL->verbose > 2) { |
481 |
|
|
fprintf(stderr, "%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, " 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_TRUE) |
499 |
|
|
|| gageStackBlurParmBoundarySet(sbpIN, nrrdBoundaryBleed, AIR_NAN) |
500 |
|
|
|| gageStackBlurManage(&ninSS, &recompute, sbpIN, |
501 |
|
|
stackFnameFormat, AIR_TRUE, NULL, |
502 |
|
|
nin, kind)) { |
503 |
|
|
airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways); |
504 |
|
|
fprintf(stderr, "%s: trouble getting volume stack:\n%s\n", me, err); |
505 |
|
|
airMopError(mop); return 1; |
506 |
|
|
} |
507 |
|
|
if (verbose > 2) { |
508 |
|
|
fprintf(stderr, "%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, " 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; |
519 |
|
|
sbpIN = NULL; |
520 |
|
|
sbp = NULL; |
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), airMopAlways); |
530 |
|
|
gageParmSet(ctx, gageParmGradMagCurvMin, gmc); |
531 |
|
|
gageParmSet(ctx, gageParmVerbose, verbose); |
532 |
|
|
gageParmSet(ctx, gageParmTwoDimZeroZ, zeroZ); |
533 |
|
|
gageParmSet(ctx, gageParmRenormalize, renorm ? AIR_TRUE : AIR_FALSE); |
534 |
|
|
gageParmSet(ctx, gageParmCheckIntegrals, AIR_TRUE); |
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_TRUE); |
544 |
|
|
gageParmSet(ctx, gageParmStackNormalizeDeriv, normdSS); |
545 |
|
|
gageParmSet(ctx, gageParmStackNormalizeDerivBias, biasSS); |
546 |
|
|
if (!E) E |= !(pvlSS = AIR_CAST(gagePerVolume **, |
547 |
|
|
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), |
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(GAGE), airFree, airMopAlways); |
562 |
|
|
fprintf(stderr, "%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, "%s: kernel support = %d^3 samples\n", me, |
572 |
|
|
2*ctx->radius); |
573 |
|
|
} |
574 |
|
|
|
575 |
|
|
if (ELL_3V_EXISTS(pntPos)) { |
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))) { |
580 |
|
|
fprintf(stderr, "%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, "%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, "%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, 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, "\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) \ |
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, |
632 |
|
|
pntPos[1] + (yo-1)*eps, |
633 |
|
|
pntPos[2] + (zo-1)*eps); |
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, "%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)); |
669 |
|
|
airMopError(mop); return 1; |
670 |
|
|
} |
671 |
|
|
if ((sbp && 3 == _npos->axis[0].size) |
672 |
|
|
|| (!sbp && 4 == _npos->axis[0].size)) { |
673 |
|
|
fprintf(stderr, "%s: have %u point coords but %s using scale-space\n", |
674 |
|
|
me, AIR_UINT(_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), airMopAlways); |
681 |
|
|
nout = nrrdNew(); |
682 |
|
|
airMopAdd(mop, nout, AIR_CAST(airMopper, nrrdNuke), airMopAlways); |
683 |
|
|
if (nrrdConvert(npos, _npos, nrrdTypeDouble) |
684 |
|
|
|| nrrdMaybeAlloc_va(nout, otype, 2, |
685 |
|
|
AIR_CAST(size_t, ansLen), |
686 |
|
|
AIR_CAST(size_t, NN))) { |
687 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
688 |
|
|
fprintf(stderr, "%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); |
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)) { |
729 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
730 |
|
|
fprintf(stderr, "%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_MAX]; |
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); |
753 |
|
|
dsiy = AIR_CAST(double, siy); |
754 |
|
|
dsiz = AIR_CAST(double, siz); |
755 |
|
|
sox = AIR_CAST(size_t, scale[0]*dsix); |
756 |
|
|
soy = AIR_CAST(size_t, scale[1]*dsiy); |
757 |
|
|
soz = AIR_CAST(size_t, scale[2]*dsiz); |
758 |
|
|
dsox = AIR_CAST(double, sox); |
759 |
|
|
dsoy = AIR_CAST(double, soy); |
760 |
|
|
dsoz = AIR_CAST(double, soz); |
761 |
|
|
rscl[0] = dsix/dsox; |
762 |
|
|
rscl[1] = dsiy/dsoy; |
763 |
|
|
rscl[2] = dsiz/dsoz; |
764 |
|
|
if (verbose) { |
765 |
|
|
fprintf(stderr, "%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, "%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(NRRD), airFree, airMopAlways); |
777 |
|
|
fprintf(stderr, "%s: trouble making ngrid:\n%s\n", me, err); |
778 |
|
|
airMopError(mop); return 1; |
779 |
|
|
} |
780 |
|
|
grid = AIR_CAST(double *, ngrid->data); |
781 |
|
|
if (nrrdCenterCell == ctx->shape->center) { |
782 |
|
|
ELL_3V_SET(min, -0.5, -0.5, -0.5); |
783 |
|
|
ELL_3V_SET(maxOut, dsox-0.5, dsoy-0.5, dsoz-0.5); |
784 |
|
|
ELL_3V_SET(maxIn, dsix-0.5, dsiy-0.5, dsiz-0.5); |
785 |
|
|
} else { |
786 |
|
|
ELL_3V_SET(min, 0, 0, 0); |
787 |
|
|
ELL_3V_SET(maxOut, dsox-1, dsoy-1, dsoz-1); |
788 |
|
|
ELL_3V_SET(maxIn, dsix-1, dsiy-1, dsiz-1); |
789 |
|
|
} |
790 |
|
|
ELL_4V_SET(grid + gridSize[0]*0, 3, |
791 |
|
|
NRRD_POS(ctx->shape->center, min[0], maxIn[0], sox, 0), |
792 |
|
|
NRRD_POS(ctx->shape->center, min[1], maxIn[1], soy, 0), |
793 |
|
|
NRRD_POS(ctx->shape->center, min[2], maxIn[2], soz, 0)); |
794 |
|
|
ELL_4V_SET(grid + gridSize[0]*1, dsox, |
795 |
|
|
AIR_DELTA(min[0], 1, maxOut[0], min[0], maxIn[0]), |
796 |
|
|
0, |
797 |
|
|
0); |
798 |
|
|
ELL_4V_SET(grid + gridSize[0]*2, dsoy, |
799 |
|
|
0, |
800 |
|
|
AIR_DELTA(min[1], 1, maxOut[1], min[1], maxIn[1]), |
801 |
|
|
0); |
802 |
|
|
ELL_4V_SET(grid + gridSize[0]*3, dsoz, |
803 |
|
|
0, |
804 |
|
|
0, |
805 |
|
|
AIR_DELTA(min[2], 1, maxOut[2], min[2], maxIn[2])); |
806 |
|
|
if (sbp) { |
807 |
|
|
if (!probeSpaceIndex) { |
808 |
|
|
double idxSS = AIR_NAN; |
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])) { |
820 |
|
|
idxSS = vi + AIR_AFFINE(sbp->sigma[vi], posSS, sbp->sigma[vi+1], |
821 |
|
|
0, 1); |
822 |
|
|
if (verbose > 1) { |
823 |
|
|
fprintf(stderr, "%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], |
826 |
|
|
0, 1)); |
827 |
|
|
} |
828 |
|
|
break; |
829 |
|
|
} |
830 |
|
|
} |
831 |
|
|
if (vi == sbp->num-1) { |
832 |
|
|
fprintf(stderr, "%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; |
850 |
|
|
if (nrrdCopy(ngrid, _ngrid)) { |
851 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
852 |
|
|
fprintf(stderr, "%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_TRUE), /* 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(GAGE), airFree, airMopAlways); |
870 |
|
|
fprintf(stderr, "%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, "probe rate = %g KHz\n", |
876 |
|
|
AIR_CAST(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 |
890 |
|
|
| NRRD_BASIC_INFO_TYPE_BIT |
891 |
|
|
| NRRD_BASIC_INFO_BLOCKSIZE_BIT |
892 |
|
|
| NRRD_BASIC_INFO_DIMENSION_BIT |
893 |
|
|
| NRRD_BASIC_INFO_CONTENT_BIT |
894 |
|
|
| NRRD_BASIC_INFO_SPACEORIGIN_BIT |
895 |
|
|
| NRRD_BASIC_INFO_OLDMIN_BIT |
896 |
|
|
| NRRD_BASIC_INFO_OLDMAX_BIT |
897 |
|
|
| NRRD_BASIC_INFO_COMMENTS_BIT |
898 |
|
|
| NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT)); |
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]); |
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)) { |
925 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
926 |
|
|
fprintf(stderr, "%s: trouble saving output:\n%s\n", me, err); |
927 |
|
|
airMopError(mop); |
928 |
|
|
return 1; |
929 |
|
|
} |
930 |
|
|
|
931 |
|
|
airMopOkay(mop); |
932 |
|
|
return 0; |
933 |
|
|
} |