Bug Summary

File:src/nrrd/axis.c
Location:line 640, column 28
Description:Assigned value is garbage or undefined

Annotated Source Code

1/*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago
4 Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5 Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6
7 This library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public License
9 (LGPL) as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11 The terms of redistributing and/or modifying this software also
12 include exceptions to the LGPL that facilitate static linking.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with this library; if not, write to Free Software Foundation, Inc.,
21 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22*/
23
24#include "nrrd.h"
25#include "privateNrrd.h"
26
27/* ------------------------------------------------------------ */
28
29void
30_nrrdAxisInfoInit(NrrdAxisInfo *axis) {
31 int dd;
32
33 if (axis) {
34 axis->size = 0;
35 axis->spacing = axis->thickness = AIR_NAN(airFloatQNaN.f);
36 axis->min = axis->max = AIR_NAN(airFloatQNaN.f);
37 for (dd=0; dd<NRRD_SPACE_DIM_MAX8; dd++) {
38 axis->spaceDirection[dd] = AIR_NAN(airFloatQNaN.f);
39 }
40 axis->center = nrrdCenterUnknown;
41 axis->kind = nrrdKindUnknown;
42 axis->label = (char *)airFree(axis->label);
43 axis->units = (char *)airFree(axis->units);
44 }
45}
46
47void
48_nrrdAxisInfoNewInit(NrrdAxisInfo *axis) {
49
50 if (axis) {
51 axis->label = NULL((void*)0);
52 axis->units = NULL((void*)0);
53 _nrrdAxisInfoInit(axis);
54 }
55}
56
57/* ------------------------------------------------------------ */
58
59/*
60******** nrrdKindIsDomain
61**
62** returns non-zero for kinds (from nrrdKind* enum) that are domain
63** axes, or independent variable axes, or resample-able axes, all
64** different ways of describing the same thing
65*/
66int
67nrrdKindIsDomain(int kind) {
68
69 return (nrrdKindDomain == kind
70 || nrrdKindSpace == kind
71 || nrrdKindTime == kind);
72}
73
74/*
75******** nrrdKindSize
76**
77** returns suggested size (length) of an axis with the given kind, or,
78** 0 if either (1) there is no suggested size because the axis is the
79** kind of an independent or domain variable or (2) the kind is invalid
80*/
81unsigned int
82nrrdKindSize(int kind) {
83 static const char me[]="nrrdKindSize";
84 unsigned int ret;
85
86 if (!( AIR_IN_OP(nrrdKindUnknown, kind, nrrdKindLast)((nrrdKindUnknown) < (kind) && (kind) < (nrrdKindLast
)) )) {
87 /* they gave us invalid or unknown kind */
88 return 0;
89 }
90
91 switch (kind) {
92 case nrrdKindDomain:
93 case nrrdKindSpace:
94 case nrrdKindTime:
95 case nrrdKindList:
96 case nrrdKindPoint:
97 case nrrdKindVector:
98 case nrrdKindCovariantVector:
99 case nrrdKindNormal:
100 ret = 0;
101 break;
102 case nrrdKindStub:
103 case nrrdKindScalar:
104 ret = 1;
105 break;
106 case nrrdKindComplex:
107 case nrrdKind2Vector:
108 ret = 2;
109 break;
110 case nrrdKind3Color:
111 case nrrdKindRGBColor:
112 case nrrdKindHSVColor:
113 case nrrdKindXYZColor:
114 ret = 3;
115 break;
116 case nrrdKind4Color:
117 case nrrdKindRGBAColor:
118 ret = 4;
119 break;
120 case nrrdKind3Vector:
121 case nrrdKind3Normal:
122 ret = 3;
123 break;
124 case nrrdKind4Vector:
125 case nrrdKindQuaternion:
126 ret = 4;
127 break;
128 case nrrdKind2DSymMatrix:
129 ret = 3;
130 break;
131 case nrrdKind2DMaskedSymMatrix:
132 ret = 4;
133 break;
134 case nrrdKind2DMatrix:
135 ret = 4;
136 break;
137 case nrrdKind2DMaskedMatrix:
138 ret = 5;
139 break;
140 case nrrdKind3DSymMatrix:
141 ret = 6;
142 break;
143 case nrrdKind3DMaskedSymMatrix:
144 ret = 7;
145 break;
146 case nrrdKind3DMatrix:
147 ret = 9;
148 break;
149 case nrrdKind3DMaskedMatrix:
150 ret = 10;
151 break;
152 default:
153 fprintf(stderr__stderrp, "%s: PANIC: nrrdKind %d not implemented!\n", me, kind);
154 ret = UINT_MAX(2147483647 *2U +1U);
155 }
156
157 return ret;
158}
159
160/*
161** _nrrdKindAltered:
162**
163** implements logic for how kind should be updated when samples
164** along the axis are altered
165*/
166int
167_nrrdKindAltered(int kindIn, int resampling) {
168 int kindOut;
169
170 if (nrrdStateKindNoop) {
171 kindOut = nrrdKindUnknown;
172 /* HEY: setting the kindOut to unknown is arguably not a no-op.
173 It is more like pointedly and stubbornly simplistic. So maybe
174 nrrdStateKindNoop could be renamed .. */
175 } else {
176 if (nrrdKindIsDomain(kindIn)
177 || (0 == nrrdKindSize(kindIn) && !resampling)) {
178 kindOut = kindIn;
179 } else {
180 kindOut = nrrdKindUnknown;
181 }
182 }
183 return kindOut;
184}
185
186/*
187** _nrrdAxisInfoCopy
188**
189** HEY: we have a void return even though this function potentially
190** involves calling airStrdup!!
191*/
192void
193_nrrdAxisInfoCopy(NrrdAxisInfo *dest, const NrrdAxisInfo *src, int bitflag) {
194 int ii;
195
196 if (!(NRRD_AXIS_INFO_SIZE_BIT(1<< 1) & bitflag)) {
197 dest->size = src->size;
198 }
199 if (!(NRRD_AXIS_INFO_SPACING_BIT(1<< 2) & bitflag)) {
200 dest->spacing = src->spacing;
201 }
202 if (!(NRRD_AXIS_INFO_THICKNESS_BIT(1<< 3) & bitflag)) {
203 dest->thickness = src->thickness;
204 }
205 if (!(NRRD_AXIS_INFO_MIN_BIT(1<< 4) & bitflag)) {
206 dest->min = src->min;
207 }
208 if (!(NRRD_AXIS_INFO_MAX_BIT(1<< 5) & bitflag)) {
209 dest->max = src->max;
210 }
211 if (!(NRRD_AXIS_INFO_SPACEDIRECTION_BIT(1<< 6) & bitflag)) {
212 for (ii=0; ii<NRRD_SPACE_DIM_MAX8; ii++) {
213 dest->spaceDirection[ii] = src->spaceDirection[ii];
214 }
215 }
216 if (!(NRRD_AXIS_INFO_CENTER_BIT(1<< 7) & bitflag)) {
217 dest->center = src->center;
218 }
219 if (!(NRRD_AXIS_INFO_KIND_BIT(1<< 8) & bitflag)) {
220 dest->kind = src->kind;
221 }
222 if (!(NRRD_AXIS_INFO_LABEL_BIT(1<< 9) & bitflag)) {
223 if (dest->label != src->label) {
224 dest->label = (char *)airFree(dest->label);
225 dest->label = (char *)airStrdup(src->label);
226 }
227 }
228 if (!(NRRD_AXIS_INFO_UNITS_BIT(1<<10) & bitflag)) {
229 if (dest->units != src->units) {
230 dest->units = (char *)airFree(dest->units);
231 dest->units = (char *)airStrdup(src->units);
232 }
233 }
234
235 return;
236}
237
238/*
239******** nrrdAxisInfoCopy()
240**
241** For copying all the per-axis peripheral information. Takes a
242** permutation "map"; map[d] tells from which axis in input should the
243** output axis d copy its information. The length of this permutation
244** array is nout->dim. If map is NULL, the identity permutation is
245** assumed. If map[i]==-1 for any i in [0,dim-1], then nothing is
246** copied into axis i of output. The "bitflag" field controls which
247** per-axis fields will NOT be copied; if bitflag==0, then all fields
248** are copied. The value of bitflag should be |'s of NRRD_AXIS_INFO_*
249** defines.
250**
251** Decided to Not use Biff, since many times map will be NULL, in
252** which case the only error is getting a NULL nrrd, or an invalid map
253** permutation, which will probably be unlikely given the contexts in
254** which this is called. For the paranoid, the integer return value
255** indicates error.
256**
257** Sun Feb 27 21:12:57 EST 2005: decided to allow nout==nin, so now
258** use a local array of NrrdAxisInfo as buffer.
259*/
260int
261nrrdAxisInfoCopy(Nrrd *nout, const Nrrd *nin, const int *axmap, int bitflag) {
262 NrrdAxisInfo axisBuffer[NRRD_DIM_MAX16];
263 const NrrdAxisInfo *axis;
264 unsigned int from, axi;
265
266 if (!(nout && nin)) {
267 return 1;
268 }
269 if (axmap) {
270 for (axi=0; axi<nout->dim; axi++) {
271 if (-1 == axmap[axi]) {
272 continue;
273 }
274 if (!AIR_IN_CL(0, axmap[axi], (int)nin->dim-1)((0) <= (axmap[axi]) && (axmap[axi]) <= ((int)nin
->dim-1))) {
275 return 3;
276 }
277 }
278 }
279 if (nout == nin) {
280 /* copy axis info to local buffer */
281 for (axi=0; axi<nin->dim; axi++) {
282 _nrrdAxisInfoNewInit(axisBuffer + axi);
283 _nrrdAxisInfoCopy(axisBuffer + axi, nin->axis + axi, bitflag);
284 }
285 axis = axisBuffer;
286 } else {
287 axis = nin->axis;
288 }
289 for (axi=0; axi<nout->dim; axi++) {
290 if (axmap && -1 == axmap[axi]) {
291 /* for this axis, we don't touch a thing */
292 continue;
293 }
294 from = axmap ? (unsigned int)axmap[axi] : axi;
295 _nrrdAxisInfoCopy(nout->axis + axi, axis + from, bitflag);
296 }
297 if (nout == nin) {
298 /* free dynamically allocated stuff */
299 for (axi=0; axi<nin->dim; axi++) {
300 _nrrdAxisInfoInit(axisBuffer + axi);
301 }
302 }
303 return 0;
304}
305
306/*
307******** nrrdAxisInfoSet_nva()
308**
309** Simple means of setting fields of the axis array in the nrrd.
310**
311** type to pass for third argument:
312** nrrdAxisInfoSize: size_t*
313** nrrdAxisInfoSpacing: double*
314** nrrdAxisInfoThickness: double*
315** nrrdAxisInfoMin: double*
316** nrrdAxisInfoMax: double*
317** nrrdAxisInfoSpaceDirection: double (*var)[NRRD_SPACE_DIM_MAX]
318** nrrdAxisInfoCenter: int*
319** nrrdAxisInfoKind: int*
320** nrrdAxisInfoLabel: char**
321** nrrdAxisInfoUnits: char**
322**
323** Note that in the case of nrrdAxisInfoSpaceDirection, we only access
324** spaceDim elements of info.V[ai] (so caller can allocate it for less
325** than NRRD_SPACE_DIM_MAX if they know what they're doing)
326*/
327void
328nrrdAxisInfoSet_nva(Nrrd *nrrd, int axInfo, const void *_info) {
329 _nrrdAxisInfoSetPtrs info;
330 int exists;
331 unsigned int ai, si, minsi;
332
333 if (!( nrrd
334 && AIR_IN_CL(1, nrrd->dim, NRRD_DIM_MAX)((1) <= (nrrd->dim) && (nrrd->dim) <= (16
))
335 && AIR_IN_OP(nrrdAxisInfoUnknown, axInfo, nrrdAxisInfoLast)((nrrdAxisInfoUnknown) < (axInfo) && (axInfo) <
(nrrdAxisInfoLast))
336 && _info )) {
337 return;
338 }
339 info.P = _info;
340
341 for (ai=0; ai<nrrd->dim; ai++) {
342 switch (axInfo) {
343 case nrrdAxisInfoSize:
344 nrrd->axis[ai].size = info.ST[ai];
345 break;
346 case nrrdAxisInfoSpacing:
347 nrrd->axis[ai].spacing = info.D[ai];
348 break;
349 case nrrdAxisInfoThickness:
350 nrrd->axis[ai].thickness = info.D[ai];
351 break;
352 case nrrdAxisInfoMin:
353 nrrd->axis[ai].min = info.D[ai];
354 break;
355 case nrrdAxisInfoMax:
356 nrrd->axis[ai].max = info.D[ai];
357 break;
358 case nrrdAxisInfoSpaceDirection:
359 /* we won't allow setting an invalid direction */
360 exists = AIR_EXISTS(info.V[ai][0])(((int)(!((info.V[ai][0]) - (info.V[ai][0])))));
361 minsi = nrrd->spaceDim;
362 for (si=0; si<nrrd->spaceDim; si++) {
363 nrrd->axis[ai].spaceDirection[si] = info.V[ai][si];
364 if (exists ^ AIR_EXISTS(info.V[ai][si])(((int)(!((info.V[ai][si]) - (info.V[ai][si])))))) {
365 minsi = 0;
366 break;
367 }
368 }
369 for (si=minsi; si<NRRD_SPACE_DIM_MAX8; si++) {
370 nrrd->axis[ai].spaceDirection[si] = AIR_NAN(airFloatQNaN.f);
371 }
372 break;
373 case nrrdAxisInfoCenter:
374 nrrd->axis[ai].center = info.I[ai];
375 break;
376 case nrrdAxisInfoKind:
377 nrrd->axis[ai].kind = info.I[ai];
378 break;
379 case nrrdAxisInfoLabel:
380 nrrd->axis[ai].label = (char *)airFree(nrrd->axis[ai].label);
381 nrrd->axis[ai].label = (char *)airStrdup(info.CP[ai]);
382 break;
383 case nrrdAxisInfoUnits:
384 nrrd->axis[ai].units = (char *)airFree(nrrd->axis[ai].units);
385 nrrd->axis[ai].units = (char *)airStrdup(info.CP[ai]);
386 break;
387 }
388 }
389 if (nrrdAxisInfoSpaceDirection == axInfo) {
390 for (ai=nrrd->dim; ai<NRRD_DIM_MAX16; ai++) {
391 for (si=0; si<NRRD_SPACE_DIM_MAX8; si++) {
392 nrrd->axis[ai].spaceDirection[si] = AIR_NAN(airFloatQNaN.f);
393 }
394 }
395 }
396 return;
397}
398
399/*
400******** nrrdAxisInfoSet_va()
401**
402** var args front-end for nrrdAxisInfoSet_nva
403**
404** types to pass, one for each dimension:
405** nrrdAxisInfoSize: size_t
406** nrrdAxisInfoSpacing: double
407** nrrdAxisInfoThickness: double
408** nrrdAxisInfoMin: double
409** nrrdAxisInfoMax: double
410** nrrdAxisInfoSpaceDirection: double*
411** nrrdAxisInfoCenter: int
412** nrrdAxisInfoKind: int
413** nrrdAxisInfoLabel: char*
414** nrrdAxisInfoUnits: char*
415*/
416void
417nrrdAxisInfoSet_va(Nrrd *nrrd, int axInfo, ...) {
418 NRRD_TYPE_BIGGESTdouble *buffer[NRRD_DIM_MAX16];
419 _nrrdAxisInfoSetPtrs info;
420 unsigned int ai, si;
421 va_list ap;
422 double *dp, svec[NRRD_DIM_MAX16][NRRD_SPACE_DIM_MAX8];
423
424 if (!( nrrd
425 && AIR_IN_CL(1, nrrd->dim, NRRD_DIM_MAX)((1) <= (nrrd->dim) && (nrrd->dim) <= (16
))
426 && AIR_IN_OP(nrrdAxisInfoUnknown, axInfo, nrrdAxisInfoLast)((nrrdAxisInfoUnknown) < (axInfo) && (axInfo) <
(nrrdAxisInfoLast)) )) {
427 return;
428 }
429
430 info.P = buffer;
431 va_start(ap, axInfo)__builtin_va_start(ap, axInfo);
432 for (ai=0; ai<nrrd->dim; ai++) {
433 switch (axInfo) {
434 case nrrdAxisInfoSize:
435 info.ST[ai] = va_arg(ap, size_t)__builtin_va_arg(ap, size_t);
436 /*
437 printf("!%s: got int[%d] = %d\n", "nrrdAxisInfoSet", d, info.I[ai]);
438 */
439 break;
440 case nrrdAxisInfoSpaceDirection:
441 dp = va_arg(ap, double*)__builtin_va_arg(ap, double*); /* punting on using info enum */
442 /*
443 printf("!%s: got dp = %lu\n", "nrrdAxisInfoSet",
444 (unsigned long)(dp));
445 */
446 for (si=0; si<nrrd->spaceDim; si++) {
447 /* nrrd->axis[ai].spaceDirection[si] = dp[si]; */
448 svec[ai][si] = dp[si];
449 }
450 for (si=nrrd->spaceDim; si<NRRD_SPACE_DIM_MAX8; si++) {
451 /* nrrd->axis[ai].spaceDirection[si] = AIR_NAN; */
452 svec[ai][si] = dp[si];
453 }
454 break;
455 case nrrdAxisInfoCenter:
456 case nrrdAxisInfoKind:
457 info.I[ai] = va_arg(ap, int)__builtin_va_arg(ap, int);
458 /*
459 printf("!%s: got int[%d] = %d\n",
460 "nrrdAxisInfoSet", d, info.I[ai]);
461 */
462 break;
463 case nrrdAxisInfoSpacing:
464 case nrrdAxisInfoThickness:
465 case nrrdAxisInfoMin:
466 case nrrdAxisInfoMax:
467 info.D[ai] = va_arg(ap, double)__builtin_va_arg(ap, double);
468 /*
469 printf("!%s: got double[%d] = %g\n",
470 "nrrdAxisInfoSet", d, info.D[ai]);
471 */
472 break;
473 case nrrdAxisInfoLabel:
474 /* we DO NOT do the airStrdup() here because this pointer value is
475 just going to be handed to nrrdAxisInfoSet_nva(), which WILL do the
476 airStrdup(); we're not violating the rules for axis labels */
477 info.CP[ai] = va_arg(ap, char *)__builtin_va_arg(ap, char *);
478 /*
479 printf("!%s: got char*[%d] = |%s|\n",
480 "nrrdAxisInfoSet", d, info.CP[ai]);
481 */
482 break;
483 case nrrdAxisInfoUnits:
484 /* see not above */
485 info.CP[ai] = va_arg(ap, char *)__builtin_va_arg(ap, char *);
486 break;
487 }
488 }
489 va_end(ap)__builtin_va_end(ap);
490
491 if (nrrdAxisInfoSpaceDirection != axInfo) {
492 /* now set the quantities which we've gotten from the var args */
493 nrrdAxisInfoSet_nva(nrrd, axInfo, info.P);
494 } else {
495 nrrdAxisInfoSet_nva(nrrd, axInfo, svec);
496 }
497
498 return;
499}
500
501/*
502******** nrrdAxisInfoGet_nva()
503**
504** get any of the axis fields into an array
505**
506** Note that getting axes labels involves implicitly allocating space
507** for them, due to the action of airStrdup(). The user is
508** responsible for free()ing these strings when done with them.
509**
510** type to pass for third argument:
511** nrrdAxisInfoSize: size_t*
512** nrrdAxisInfoSpacing: double*
513** nrrdAxisInfoThickness: double*
514** nrrdAxisInfoMin: double*
515** nrrdAxisInfoMax: double*
516** nrrdAxisInfoSpaceDirection: double (*var)[NRRD_SPACE_DIM_MAX]
517** nrrdAxisInfoCenter: int*
518** nrrdAxisInfoKind: int*
519** nrrdAxisInfoLabel: char**
520** nrrdAxisInfoUnits: char**
521*/
522void
523nrrdAxisInfoGet_nva(const Nrrd *nrrd, int axInfo, void *_info) {
524 _nrrdAxisInfoGetPtrs info;
525 unsigned int ai, si;
526
527 if (!( nrrd
528 && AIR_IN_CL(1, nrrd->dim, NRRD_DIM_MAX)((1) <= (nrrd->dim) && (nrrd->dim) <= (16
))
529 && AIR_IN_OP(nrrdAxisInfoUnknown, axInfo, nrrdAxisInfoLast)((nrrdAxisInfoUnknown) < (axInfo) && (axInfo) <
(nrrdAxisInfoLast)) )) {
530 return;
531 }
532
533 info.P = _info;
534 for (ai=0; ai<nrrd->dim; ai++) {
535 switch (axInfo) {
536 case nrrdAxisInfoSize:
537 info.ST[ai] = nrrd->axis[ai].size;
538 break;
539 case nrrdAxisInfoSpacing:
540 info.D[ai] = nrrd->axis[ai].spacing;
541 break;
542 case nrrdAxisInfoThickness:
543 info.D[ai] = nrrd->axis[ai].thickness;
544 break;
545 case nrrdAxisInfoMin:
546 info.D[ai] = nrrd->axis[ai].min;
547 break;
548 case nrrdAxisInfoMax:
549 info.D[ai] = nrrd->axis[ai].max;
550 break;
551 case nrrdAxisInfoSpaceDirection:
552 for (si=0; si<nrrd->spaceDim; si++) {
553 info.V[ai][si] = nrrd->axis[ai].spaceDirection[si];
554 }
555 for (si=nrrd->spaceDim; si<NRRD_SPACE_DIM_MAX8; si++) {
556 info.V[ai][si] = AIR_NAN(airFloatQNaN.f);
557 }
558 break;
559 case nrrdAxisInfoCenter:
560 info.I[ai] = nrrd->axis[ai].center;
561 break;
562 case nrrdAxisInfoKind:
563 info.I[ai] = nrrd->axis[ai].kind;
564 break;
565 case nrrdAxisInfoLabel:
566 /* note airStrdup()! */
567 info.CP[ai] = airStrdup(nrrd->axis[ai].label);
568 break;
569 case nrrdAxisInfoUnits:
570 /* note airStrdup()! */
571 info.CP[ai] = airStrdup(nrrd->axis[ai].units);
572 break;
573 }
574 }
575 if (nrrdAxisInfoSpaceDirection == axInfo) {
576 for (ai=nrrd->dim; ai<NRRD_DIM_MAX16; ai++) {
577 for (si=0; si<NRRD_SPACE_DIM_MAX8; si++) {
578 info.V[ai][si] = AIR_NAN(airFloatQNaN.f);
579 }
580 }
581 }
582 return;
583}
584
585/*
586** types to pass, one for each dimension:
587** nrrdAxisInfoSize: size_t*
588** nrrdAxisInfoSpacing: double*
589** nrrdAxisInfoThickness: double*
590** nrrdAxisInfoMin: double*
591** nrrdAxisInfoMax: double*
592** nrrdAxisInfoSpaceDirection: double*
593** nrrdAxisInfoCenter: int*
594** nrrdAxisInfoKind: int*
595** nrrdAxisInfoLabel: char**
596** nrrdAxisInfoUnits: char**
597*/
598void
599nrrdAxisInfoGet_va(const Nrrd *nrrd, int axInfo, ...) {
600 void *buffer[NRRD_DIM_MAX16], *ptr;
601 _nrrdAxisInfoGetPtrs info;
602 unsigned int ai, si;
603 va_list ap;
604 double svec[NRRD_DIM_MAX16][NRRD_SPACE_DIM_MAX8];
605
606 if (!( nrrd
1
Taking false branch
→
607 && AIR_IN_CL(1, nrrd->dim, NRRD_DIM_MAX)((1) <= (nrrd->dim) && (nrrd->dim) <= (16
))
608 && AIR_IN_OP(nrrdAxisInfoUnknown, axInfo, nrrdAxisInfoLast)((nrrdAxisInfoUnknown) < (axInfo) && (axInfo) <
(nrrdAxisInfoLast)) )) {
609 return;
610 }
611
612 if (nrrdAxisInfoSpaceDirection != axInfo) {
2
←
Assuming 'axInfo' is not equal to nrrdAxisInfoSpaceDirection
→
3
←
Taking true branch
→
613 info.P = buffer;
614 nrrdAxisInfoGet_nva(nrrd, axInfo, info.P);
615 } else {
616 nrrdAxisInfoGet_nva(nrrd, axInfo, svec);
617 }
618
619 va_start(ap, axInfo)__builtin_va_start(ap, axInfo);
620 for (ai=0; ai<nrrd->dim; ai++) {
4
←
Loop condition is true. Entering loop body
→
621 ptr = va_arg(ap, void*)__builtin_va_arg(ap, void*);
622 /*
623 printf("!%s(%d): ptr = %lu\n",
624 "nrrdAxisInfoGet", d, (unsigned long)ptr);
625 */
626 switch (axInfo) {
5
←
Control jumps to 'case nrrdAxisInfoSpaceDirection:' at line 638
→
627 case nrrdAxisInfoSize:
628 *((size_t*)ptr) = info.ST[ai];
629 break;
630 case nrrdAxisInfoSpacing:
631 case nrrdAxisInfoThickness:
632 case nrrdAxisInfoMin:
633 case nrrdAxisInfoMax:
634 *((double*)ptr) = info.D[ai];
635 /* printf("!%s: got double[%d] = %lg\n", "nrrdAxisInfoGet", d,
636 *((double*)ptr)); */
637 break;
638 case nrrdAxisInfoSpaceDirection:
639 for (si=0; si<nrrd->spaceDim; si++) {
6
←
Loop condition is true. Entering loop body
→
640 ((double*)ptr)[si] = svec[ai][si];
7
←
Assigned value is garbage or undefined
641 }
642 for (si=nrrd->spaceDim; si<NRRD_SPACE_DIM_MAX8; si++) {
643 ((double*)ptr)[si] = AIR_NAN(airFloatQNaN.f);
644 }
645 break;
646 case nrrdAxisInfoCenter:
647 case nrrdAxisInfoKind:
648 *((int*)ptr) = info.I[ai];
649 /* printf("!%s: got int[%d] = %d\n",
650 "nrrdAxisInfoGet", d, *((int*)ptr)); */
651 break;
652 case nrrdAxisInfoLabel:
653 case nrrdAxisInfoUnits:
654 /* we DO NOT do the airStrdup() here because this pointer value just
655 came from nrrdAxisInfoGet_nva(), which already did the airStrdup() */
656 *((char**)ptr) = info.CP[ai];
657 /* printf("!%s: got char*[%d] = |%s|\n", "nrrdAxisInfoSet", d,
658 *((char**)ptr)); */
659 break;
660 }
661 }
662 va_end(ap)__builtin_va_end(ap);
663
664 return;
665}
666
667/*
668** _nrrdCenter()
669**
670** for nrrdCenterCell and nrrdCenterNode, return will be the same
671** as input. Converts nrrdCenterUnknown into nrrdDefaultCenter,
672** and then clamps to (nrrdCenterUnknown+1, nrrdCenterLast-1).
673**
674** Thus, this ALWAYS returns nrrdCenterNode or nrrdCenterCell
675** (as long as those are the only two centering schemes).
676*/
677int
678_nrrdCenter(int center) {
679
680 center = (nrrdCenterUnknown == center
681 ? nrrdDefaultCenter
682 : center);
683 center = AIR_CLAMP(nrrdCenterUnknown+1, center, nrrdCenterLast-1)((center) < (nrrdCenterUnknown+1) ? (nrrdCenterUnknown+1) :
((center) > (nrrdCenterLast-1) ? (nrrdCenterLast-1) : (center
)));
684 return center;
685}
686
687int
688_nrrdCenter2(int center, int defCenter) {
689
690 center = (nrrdCenterUnknown == center
691 ? defCenter
692 : center);
693 center = AIR_CLAMP(nrrdCenterUnknown+1, center, nrrdCenterLast-1)((center) < (nrrdCenterUnknown+1) ? (nrrdCenterUnknown+1) :
((center) > (nrrdCenterLast-1) ? (nrrdCenterLast-1) : (center
)));
694 return center;
695}
696
697
698/*
699******** nrrdAxisInfoPos()
700**
701** given a nrrd, an axis, and a (floating point) index space position,
702** return the position implied the axis's min, max, and center
703** Does the opposite of nrrdAxisIdx().
704**
705** does not use biff
706*/
707double
708nrrdAxisInfoPos(const Nrrd *nrrd, unsigned int ax, double idx) {
709 int center;
710 size_t size;
711 double min, max;
712
713 if (!( nrrd && ax <= nrrd->dim-1 )) {
714 return AIR_NAN(airFloatQNaN.f);
715 }
716 center = _nrrdCenter(nrrd->axis[ax].center);
717 min = nrrd->axis[ax].min;
718 max = nrrd->axis[ax].max;
719 size = nrrd->axis[ax].size;
720
721 return NRRD_POS(center, min, max, size, idx)(nrrdCenterCell == (center) ? ( ((double)(((max)))-(((min))))
*((double)(((idx)) + 0.5)-(0)) / ((double)(((size)))-(0)) + (
((min)))) : ( ((double)(((max)))-(((min))))*((double)(((idx))
)-(0)) / ((double)(((size))-1)-(0)) + (((min)))));
722}
723
724/*
725******** nrrdAxisInfoIdx()
726**
727** given a nrrd, an axis, and a (floating point) world space position,
728** return the index implied the axis's min, max, and center.
729** Does the opposite of nrrdAxisPos().
730**
731** does not use biff
732*/
733double
734nrrdAxisInfoIdx(const Nrrd *nrrd, unsigned int ax, double pos) {
735 int center;
736 size_t size;
737 double min, max;
738
739 if (!( nrrd && ax <= nrrd->dim-1 )) {
740 return AIR_NAN(airFloatQNaN.f);
741 }
742 center = _nrrdCenter(nrrd->axis[ax].center);
743 min = nrrd->axis[ax].min;
744 max = nrrd->axis[ax].max;
745 size = nrrd->axis[ax].size;
746
747 return NRRD_IDX(center, min, max, size, pos)(nrrdCenterCell == (center) ? (( ((double)(((size)))-(0))*((double
)(((pos)))-(((min)))) / ((double)(((max)))-(((min)))) + (0)) -
0.5) : ( ((double)(((size))-1)-(0))*((double)(((pos)))-(((min
)))) / ((double)(((max)))-(((min)))) + (0)));
748}
749
750/*
751******** nrrdAxisInfoPosRange()
752**
753** given a nrrd, an axis, and two (floating point) index space positions,
754** return the range of positions implied the axis's min, max, and center
755** The opposite of nrrdAxisIdxRange()
756*/
757void
758nrrdAxisInfoPosRange(double *loP, double *hiP,
759 const Nrrd *nrrd, unsigned int ax,
760 double loIdx, double hiIdx) {
761 int center, flip = 0;
762 size_t size;
763 double min, max, tmp;
764
765 if (!( loP && hiP && nrrd && ax <= nrrd->dim-1 )) {
766 if (loP) *loP = AIR_NAN(airFloatQNaN.f);
767 if (hiP) *hiP = AIR_NAN(airFloatQNaN.f);
768 return;
769 }
770 center = _nrrdCenter(nrrd->axis[ax].center);
771 min = nrrd->axis[ax].min;
772 max = nrrd->axis[ax].max;
773 size = nrrd->axis[ax].size;
774
775 if (loIdx > hiIdx) {
776 flip = 1;
777 tmp = loIdx; loIdx = hiIdx; hiIdx = tmp;
778 }
779 if (nrrdCenterCell == center) {
780 *loP = AIR_AFFINE(0, loIdx, size, min, max)( ((double)(max)-(min))*((double)(loIdx)-(0)) / ((double)(size
)-(0)) + (min));
781 *hiP = AIR_AFFINE(0, hiIdx+1, size, min, max)( ((double)(max)-(min))*((double)(hiIdx+1)-(0)) / ((double)(size
)-(0)) + (min));
782 } else {
783 *loP = AIR_AFFINE(0, loIdx, size-1, min, max)( ((double)(max)-(min))*((double)(loIdx)-(0)) / ((double)(size
-1)-(0)) + (min));
784 *hiP = AIR_AFFINE(0, hiIdx, size-1, min, max)( ((double)(max)-(min))*((double)(hiIdx)-(0)) / ((double)(size
-1)-(0)) + (min));
785 }
786 if (flip) {
787 tmp = *loP; *loP = *hiP; *hiP = tmp;
788 }
789
790 return;
791}
792
793/*
794******** nrrdAxisInfoIdxRange()
795**
796** given a nrrd, an axis, and two (floating point) world space positions,
797** return the range of index space implied the axis's min, max, and center
798** The opposite of nrrdAxisPosRange().
799**
800** Actually- there are situations where sending an interval through
801** nrrdAxisIdxRange -> nrrdAxisPosRange -> nrrdAxisIdxRange
802** such as in cell centering, when the range of positions given does
803** not even span one sample. Such as:
804** axis->size = 4, axis->min = -4, axis->max = 4, loPos = 0, hiPos = 1
805** --> nrrdAxisIdxRange == (2, 1.5) --> nrrdAxisPosRange == (2, -1)
806** The basic problem is that because of the 0.5 offset inherent in
807** cell centering, there are situations where (in terms of the arguments
808** to nrrdAxisIdxRange()) loPos < hiPos, but *loP > *hiP.
809*/
810void
811nrrdAxisInfoIdxRange(double *loP, double *hiP,
812 const Nrrd *nrrd, unsigned int ax,
813 double loPos, double hiPos) {
814 int center, flip = 0;
815 size_t size;
816 double min, max, tmp;
817
818 if (!( loP && hiP && nrrd && ax <= nrrd->dim-1 )) {
819 *loP = *hiP = AIR_NAN(airFloatQNaN.f);
820 return;
821 }
822 center = _nrrdCenter(nrrd->axis[ax].center);
823 min = nrrd->axis[ax].min;
824 max = nrrd->axis[ax].max;
825 size = nrrd->axis[ax].size;
826
827 if (loPos > hiPos) {
828 flip = 1;
829 tmp = loPos; loPos = hiPos; hiPos = tmp;
830 }
831 if (nrrdCenterCell == center) {
832 if (min < max) {
833 *loP = AIR_AFFINE(min, loPos, max, 0, size)( ((double)(size)-(0))*((double)(loPos)-(min)) / ((double)(max
)-(min)) + (0));
834 *hiP = AIR_AFFINE(min, hiPos, max, -1, size-1)( ((double)(size-1)-(-1))*((double)(hiPos)-(min)) / ((double)
(max)-(min)) + (-1));
835 } else {
836 *loP = AIR_AFFINE(min, loPos, max, -1, size-1)( ((double)(size-1)-(-1))*((double)(loPos)-(min)) / ((double)
(max)-(min)) + (-1));
837 *hiP = AIR_AFFINE(min, hiPos, max, 0, size)( ((double)(size)-(0))*((double)(hiPos)-(min)) / ((double)(max
)-(min)) + (0));
838 }
839 } else {
840 *loP = AIR_AFFINE(min, loPos, max, 0, size-1)( ((double)(size-1)-(0))*((double)(loPos)-(min)) / ((double)(
max)-(min)) + (0));
841 *hiP = AIR_AFFINE(min, hiPos, max, 0, size-1)( ((double)(size-1)-(0))*((double)(hiPos)-(min)) / ((double)(
max)-(min)) + (0));
842 }
843 if (flip) {
844 tmp = *loP; *loP = *hiP; *hiP = tmp;
845 }
846
847 return;
848}
849
850void
851nrrdAxisInfoSpacingSet(Nrrd *nrrd, unsigned int ax) {
852 int sign;
853 double min, max, tmp;
854
855 if (!( nrrd && ax <= nrrd->dim-1 )) {
856 return;
857 }
858
859 min = nrrd->axis[ax].min;
860 max = nrrd->axis[ax].max;
861 if (!( AIR_EXISTS(min)(((int)(!((min) - (min))))) && AIR_EXISTS(max)(((int)(!((max) - (max))))) )) {
862 /* there's no actual basis on which to set the spacing information,
863 but we have to set it something, so here goes .. */
864 nrrd->axis[ax].spacing = nrrdDefaultSpacing;
865 return;
866 }
867
868 if (min > max) {
869 tmp = min; min = max; max = tmp;
870 sign = -1;
871 } else {
872 sign = 1;
873 }
874
875 /* the skinny */
876 nrrd->axis[ax].spacing = NRRD_SPACING(_nrrdCenter(nrrd->axis[ax].center),(nrrdCenterCell == _nrrdCenter(nrrd->axis[ax].center) ? ((
max) - (min))/((double)(nrrd->axis[ax].size)) : ((max) - (
min))/(((double)((nrrd->axis[ax].size)- 1))))
877 min, max, nrrd->axis[ax].size)(nrrdCenterCell == _nrrdCenter(nrrd->axis[ax].center) ? ((
max) - (min))/((double)(nrrd->axis[ax].size)) : ((max) - (
min))/(((double)((nrrd->axis[ax].size)- 1))));
878 nrrd->axis[ax].spacing *= sign;
879
880 return;
881}
882
883void
884nrrdAxisInfoMinMaxSet(Nrrd *nrrd, unsigned int ax, int defCenter) {
885 int center;
886 double spacing;
887
888 if (!( nrrd && ax <= nrrd->dim-1 )) {
889 return;
890 }
891
892 center = _nrrdCenter2(nrrd->axis[ax].center, defCenter);
893 spacing = nrrd->axis[ax].spacing;
894 if (!AIR_EXISTS(spacing)(((int)(!((spacing) - (spacing))))))
895 spacing = nrrdDefaultSpacing;
896 if (nrrdCenterCell == center) {
897 nrrd->axis[ax].min = 0;
898 nrrd->axis[ax].max = spacing*AIR_CAST(double, nrrd->axis[ax].size)((double)(nrrd->axis[ax].size));
899 } else {
900 nrrd->axis[ax].min = 0;
901 nrrd->axis[ax].max = spacing*AIR_CAST(double, nrrd->axis[ax].size - 1)((double)(nrrd->axis[ax].size - 1));
902 }
903
904 return;
905}
906/* ---- BEGIN non-NrrdIO */
907
908/*
909** not using the value comparators in accessors.c because of their
910** slightly strange behavior WRT infinity (+inf < -42). This code
911** may eventually warrant wider availability, for now its here but
912** accessible to nrrd files via privateNrrd.h
913*/
914int
915_nrrdDblcmp(double aa, double bb) {
916 int nna, nnb, ret;
917
918 nna = AIR_EXISTS(aa)(((int)(!((aa) - (aa))))) || !airIsNaN(aa);
919 nnb = AIR_EXISTS(bb)(((int)(!((bb) - (bb))))) || !airIsNaN(bb);
920 if (nna && nnb) {
921 /* both either exist or are an infinity */
922 ret = (aa < bb
923 ? -1
924 : (aa > bb
925 ? 1
926 : 0));
927 } else {
928 /* one or the other is NaN */
929 ret = (nna < nnb
930 ? -1
931 : (nna > nnb
932 ? 1
933 : 0));
934 }
935 return ret;
936}
937
938/*
939******** nrrdAxisInfoCompare
940**
941** compares all fields in the NrrdAxisInfoCompare
942**
943** See comment about logic of return value above nrrdCompare()
944**
945** NOTE: the structure of this code is very similar to that of
946** nrrdCompare, and any improvements here should be reflected there
947*/
948int
949nrrdAxisInfoCompare(const NrrdAxisInfo *axisA, const NrrdAxisInfo *axisB,
950 int *differ, char explain[AIR_STRLEN_LARGE(512+1)]) {
951 static const char me[]="nrrdAxisInfoCompare";
952 unsigned int saxi;
953
954 if (!(axisA && axisB && differ)) {
955 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer (%p, %p, or %p)", me,
956 AIR_CVOIDP(axisA)((const void *)(axisA)), AIR_CVOIDP(axisB)((const void *)(axisB)), AIR_VOIDP(differ)((void *)(differ)));
957 return 1;
958 }
959
960 if (explain) {
961 strcpy(explain, "")__builtin___strcpy_chk (explain, "", __builtin_object_size (explain
, 2 > 1 ? 1 : 0));
962 }
963 if (axisA->size != axisB->size) {
964 char stmp1[AIR_STRLEN_SMALL(128+1)], stmp2[AIR_STRLEN_SMALL(128+1)];
965 *differ = axisA->size < axisB->size ? -1 : 1;
966 if (explain) {
967 sprintf(explain, "axisA->size=%s %s axisB->size=%s",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->size=%s %s axisB->size=%s"
, airSprintSize_t(stmp1, axisA->size), *differ < 0 ? "<"
: ">", airSprintSize_t(stmp2, axisB->size))
968 airSprintSize_t(stmp1, axisA->size),__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->size=%s %s axisB->size=%s"
, airSprintSize_t(stmp1, axisA->size), *differ < 0 ? "<"
: ">", airSprintSize_t(stmp2, axisB->size))
969 *differ < 0 ? "<" : ">",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->size=%s %s axisB->size=%s"
, airSprintSize_t(stmp1, axisA->size), *differ < 0 ? "<"
: ">", airSprintSize_t(stmp2, axisB->size))
970 airSprintSize_t(stmp2, axisB->size))__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->size=%s %s axisB->size=%s"
, airSprintSize_t(stmp1, axisA->size), *differ < 0 ? "<"
: ">", airSprintSize_t(stmp2, axisB->size));
971 }
972 return 0;
973 }
974
975#define DOUBLE_COMPARE(VAL, STR) \
976 *differ = _nrrdDblcmp(axisA->VAL, axisB->VAL); \
977 if (*differ) { \
978 if (explain) { \
979 sprintf(explain, "axisA->%s %.17g %s axisB->%s %.17g", \__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->%s %.17g %s axisB->%s %.17g"
, STR, axisA->VAL, *differ < 0 ? "<" : ">", STR, axisB
->VAL)
980 STR, axisA->VAL, *differ < 0 ? "<" : ">", \__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->%s %.17g %s axisB->%s %.17g"
, STR, axisA->VAL, *differ < 0 ? "<" : ">", STR, axisB
->VAL)
981 STR, axisB->VAL)__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->%s %.17g %s axisB->%s %.17g"
, STR, axisA->VAL, *differ < 0 ? "<" : ">", STR, axisB
->VAL); \
982 } \
983 return 0; \
984 }
985
986 DOUBLE_COMPARE(spacing, "spacing");
987 DOUBLE_COMPARE(thickness, "thickness");
988 DOUBLE_COMPARE(min, "min");
989 DOUBLE_COMPARE(max, "max");
990 for (saxi=0; saxi<NRRD_SPACE_DIM_MAX8; saxi++) {
991 char stmp[AIR_STRLEN_SMALL(128+1)];
992 sprintf(stmp, "spaceDirection[%u]", saxi)__builtin___sprintf_chk (stmp, 0, __builtin_object_size (stmp
, 2 > 1 ? 1 : 0), "spaceDirection[%u]", saxi);
993 DOUBLE_COMPARE(spaceDirection[saxi], stmp);
994 }
995#undef DOUBLE_COMPARE
996
997 if (axisA->center != axisB->center) {
998 *differ = axisA->center < axisB->center ? -1 : 1;
999 if (explain) {
1000 sprintf(explain, "axisA->center %s %s axisB->center %s",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->center %s %s axisB->center %s"
, airEnumStr(nrrdCenter, axisA->center), *differ < 0 ? "<"
: ">", airEnumStr(nrrdCenter, axisB->center))
1001 airEnumStr(nrrdCenter, axisA->center),__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->center %s %s axisB->center %s"
, airEnumStr(nrrdCenter, axisA->center), *differ < 0 ? "<"
: ">", airEnumStr(nrrdCenter, axisB->center))
1002 *differ < 0 ? "<" : ">",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->center %s %s axisB->center %s"
, airEnumStr(nrrdCenter, axisA->center), *differ < 0 ? "<"
: ">", airEnumStr(nrrdCenter, axisB->center))
1003 airEnumStr(nrrdCenter, axisB->center))__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->center %s %s axisB->center %s"
, airEnumStr(nrrdCenter, axisA->center), *differ < 0 ? "<"
: ">", airEnumStr(nrrdCenter, axisB->center));
1004 }
1005 return 0;
1006 }
1007 if (axisA->kind != axisB->kind) {
1008 *differ = axisA->kind < axisB->kind ? -1 : 1;
1009 if (explain) {
1010 sprintf(explain, "axisA->kind %s %s axisB->kind %s",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->kind %s %s axisB->kind %s"
, airEnumStr(nrrdKind, axisA->kind), *differ < 0 ? "<"
: ">", airEnumStr(nrrdKind, axisB->kind))
1011 airEnumStr(nrrdKind, axisA->kind),__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->kind %s %s axisB->kind %s"
, airEnumStr(nrrdKind, axisA->kind), *differ < 0 ? "<"
: ">", airEnumStr(nrrdKind, axisB->kind))
1012 *differ < 0 ? "<" : ">",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->kind %s %s axisB->kind %s"
, airEnumStr(nrrdKind, axisA->kind), *differ < 0 ? "<"
: ">", airEnumStr(nrrdKind, axisB->kind))
1013 airEnumStr(nrrdKind, axisB->kind))__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->kind %s %s axisB->kind %s"
, airEnumStr(nrrdKind, axisA->kind), *differ < 0 ? "<"
: ">", airEnumStr(nrrdKind, axisB->kind));
1014 }
1015 return 0;
1016 }
1017 *differ = airStrcmp(axisA->label, axisB->label);
1018 if (*differ) {
1019 if (explain) {
1020 /* can't safely print whole labels because of fixed-size of explain */
1021 sprintf(explain, "axisA->label %s axisB->label",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->label %s axisB->label", *differ
< 0 ? "<" : ">")
1022 *differ < 0 ? "<" : ">")__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->label %s axisB->label", *differ
< 0 ? "<" : ">");
1023 if (strlen(explain) + airStrlen(axisA->label)
1024 + airStrlen(axisB->label)
1025 + 2*strlen(" \"\" ") + 1 < AIR_STRLEN_LARGE(512+1)) {
1026 /* ok, we can print them */
1027 sprintf(explain, "axisA->label \"%s\" %s axisB->label \"%s\"",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->label \"%s\" %s axisB->label \"%s\""
, axisA->label ? axisA->label : "", *differ < 0 ? "<"
: ">", axisB->label ? axisB->label : "")
1028 axisA->label ? axisA->label : "",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->label \"%s\" %s axisB->label \"%s\""
, axisA->label ? axisA->label : "", *differ < 0 ? "<"
: ">", axisB->label ? axisB->label : "")
1029 *differ < 0 ? "<" : ">",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->label \"%s\" %s axisB->label \"%s\""
, axisA->label ? axisA->label : "", *differ < 0 ? "<"
: ">", axisB->label ? axisB->label : "")
1030 axisB->label ? axisB->label : "")__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->label \"%s\" %s axisB->label \"%s\""
, axisA->label ? axisA->label : "", *differ < 0 ? "<"
: ">", axisB->label ? axisB->label : "");
1031 }
1032 }
1033 return 0;
1034 }
1035 *differ = airStrcmp(axisA->units, axisB->units);
1036 if (*differ) {
1037 if (explain) {
1038 /* can't print whole string because of fixed-size of explain */
1039 sprintf(explain, "axisA->units %s axisB->units",__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->units %s axisB->units", *differ
< 0 ? "<" : ">")
1040 *differ < 0 ? "<" : ">")__builtin___sprintf_chk (explain, 0, __builtin_object_size (explain
, 2 > 1 ? 1 : 0), "axisA->units %s axisB->units", *differ
< 0 ? "<" : ">");
1041 }
1042 return 0;
1043 }
1044
1045 return 0;
1046}
1047/* ---- END non-NrrdIO */
1048
1049/*
1050******** nrrdDomainAxesGet
1051**
1052** Based on the per-axis "kind" field, learns which are the domain
1053** (resample-able) axes of an image, in other words, the axes which
1054** correspond to independent variables. The return value is the
1055** number of domain axes, and that many values are set in the given
1056** axisIdx[] array
1057**
1058** NOTE: this takes a wild guess that an unset (nrrdKindUnknown) kind
1059** is a domain axis.
1060*/
1061unsigned int
1062nrrdDomainAxesGet(const Nrrd *nrrd, unsigned int axisIdx[NRRD_DIM_MAX16]) {
1063 unsigned int domAxi, axi;
1064
1065 if (!( nrrd && axisIdx )) {
1066 return 0;
1067 }
1068 domAxi = 0;
1069 for (axi=0; axi<nrrd->dim; axi++) {
1070 if (nrrdKindUnknown == nrrd->axis[axi].kind
1071 || nrrdKindIsDomain(nrrd->axis[axi].kind)) {
1072 axisIdx[domAxi++] = axi;
1073 }
1074 }
1075 return domAxi;
1076}
1077
1078int
1079_nrrdSpaceVecExists(const Nrrd *nrrd, unsigned int axi) {
1080 unsigned int sai;
1081 int ret;
1082
1083 if (!( nrrd && axi < nrrd->dim && nrrd->spaceDim )) {
1084 ret = AIR_FALSE0;
1085 } else {
1086 ret = AIR_TRUE1;
1087 for (sai=0; sai<nrrd->spaceDim; sai++) {
1088 ret &= AIR_EXISTS(nrrd->axis[axi].spaceDirection[sai])(((int)(!((nrrd->axis[axi].spaceDirection[sai]) - (nrrd->
axis[axi].spaceDirection[sai])))));
1089 }
1090 }
1091 return ret;
1092}
1093
1094unsigned int
1095nrrdSpatialAxesGet(const Nrrd *nrrd, unsigned int axisIdx[NRRD_DIM_MAX16]) {
1096 unsigned int spcAxi, axi;
1097
1098 if (!( nrrd && axisIdx && nrrd->spaceDim)) {
1099 return 0;
1100 }
1101 spcAxi = 0;
1102 for (axi=0; axi<nrrd->dim; axi++) {
1103 if (_nrrdSpaceVecExists(nrrd, axi)) {
1104 axisIdx[spcAxi++] = axi;
1105 }
1106 }
1107 return spcAxi;
1108}
1109
1110/*
1111******** nrrdRangeAxesGet
1112**
1113** Based on the per-axis "kind" field, learns which are the range
1114** (non-resample-able) axes of an image, in other words, the axes
1115** which correspond to dependent variables. The return value is the
1116** number of range axes; that number of values are set in the given
1117** axisIdx[] array
1118**
1119** Note: this really is as simple as returning the complement of the
1120** axis selected by nrrdDomainAxesGet()
1121*/
1122unsigned int
1123nrrdRangeAxesGet(const Nrrd *nrrd, unsigned int axisIdx[NRRD_DIM_MAX16]) {
1124 unsigned int domNum, domIdx[NRRD_DIM_MAX16], rngAxi, axi, ii, isDom;
1125
1126 if (!( nrrd && axisIdx )) {
1127 return 0;
1128 }
1129 domNum = nrrdDomainAxesGet(nrrd, domIdx);
1130 rngAxi = 0;
1131 for (axi=0; axi<nrrd->dim; axi++) {
1132 isDom = AIR_FALSE0;
1133 for (ii=0; ii<domNum; ii++) { /* yes, inefficient */
1134 isDom |= axi == domIdx[ii];
1135 }
1136 if (!isDom) {
1137 axisIdx[rngAxi++] = axi;
1138 }
1139 }
1140 return rngAxi;
1141}
1142
1143unsigned int
1144nrrdNonSpatialAxesGet(const Nrrd *nrrd, unsigned int axisIdx[NRRD_DIM_MAX16]) {
1145 unsigned int spcNum, spcIdx[NRRD_DIM_MAX16], nspAxi, axi, ii, isSpc;
1146
1147 if (!( nrrd && axisIdx )) {
1148 return 0;
1149 }
1150 /* HEY: copy and paste, should refactor with above */
1151 spcNum = nrrdSpatialAxesGet(nrrd, spcIdx);
1152 nspAxi = 0;
1153 for (axi=0; axi<nrrd->dim; axi++) {
1154 isSpc = AIR_FALSE0;
1155 for (ii=0; ii<spcNum; ii++) { /* yes, inefficient */
1156 isSpc |= axi == spcIdx[ii];
1157 }
1158 if (!isSpc) {
1159 axisIdx[nspAxi++] = axi;
1160 }
1161 }
1162 return nspAxi;
1163}
1164
1165
1166/*
1167******** nrrdSpacingCalculate
1168**
1169** Determine nrrdSpacingStatus, and whatever can be calculated about
1170** spacing for a given axis. Takes a nrrd, an axis, a double pointer
1171** (for returning a scalar), a space vector, and an int pointer for
1172** returning the known length of the space vector.
1173**
1174** The behavior of what has been set by the function is determined by
1175** the return value, which takes values from the nrrdSpacingStatus*
1176** enum, as follows:
1177**
1178** returned status value: what it means, and what it set
1179** ---------------------------------------------------------------------------
1180** nrrdSpacingStatusUnknown Something about the given arguments is
1181** invalid.
1182** *spacing = NaN,
1183** vector = all NaNs
1184**
1185** nrrdSpacingStatusNone There is no spacing info at all:
1186** *spacing = NaN,
1187** vector = all NaNs
1188**
1189** nrrdSpacingStatusScalarNoSpace There is no surrounding space, but the
1190** axis's spacing was known.
1191** *spacing = axis->spacing,
1192** vector = all NaNs
1193**
1194** nrrdSpacingStatusScalarWithSpace There *is* a surrounding space, but the
1195** given axis does not live in that space,
1196** because it has no space direction. Caller
1197** may want to think about what's going on.
1198** *spacing = axis->spacing,
1199** vector = all NaNs
1200**
1201** nrrdSpacingStatusDirection There is a surrounding space, in which
1202** this axis has a direction V:
1203** *spacing = |V| (length of direction),
1204** vector = V/|V| (normalized direction)
1205** NOTE: it is still possible for both
1206** *spacing and vector to be all NaNs!!
1207*/
1208int
1209nrrdSpacingCalculate(const Nrrd *nrrd, unsigned int ax,
1210 double *spacing, double vector[NRRD_SPACE_DIM_MAX8]) {
1211 int ret;
1212
1213 if (!( nrrd && spacing && vector
1214 && ax <= nrrd->dim-1
1215 && !_nrrdCheck(nrrd, AIR_FALSE0, AIR_FALSE0) )) {
1216 /* there's a problem with the arguments. Note: the _nrrdCheck()
1217 call does not check on non-NULL-ity of nrrd->data */
1218 ret = nrrdSpacingStatusUnknown;
1219 if (spacing) {
1220 *spacing = AIR_NAN(airFloatQNaN.f);
1221 }
1222 if (vector) {
1223 nrrdSpaceVecSetNaN(vector);
1224 }
1225 } else {
1226 if (AIR_EXISTS(nrrd->axis[ax].spacing)(((int)(!((nrrd->axis[ax].spacing) - (nrrd->axis[ax].spacing
)))))) {
1227 if (nrrd->spaceDim > 0) {
1228 ret = nrrdSpacingStatusScalarWithSpace;
1229 } else {
1230 ret = nrrdSpacingStatusScalarNoSpace;
1231 }
1232 *spacing = nrrd->axis[ax].spacing;
1233 nrrdSpaceVecSetNaN(vector);
1234 } else {
1235 if (nrrd->spaceDim > 0 && _nrrdSpaceVecExists(nrrd, ax)) {
1236 ret = nrrdSpacingStatusDirection;
1237 *spacing = nrrdSpaceVecNorm(nrrd->spaceDim,
1238 nrrd->axis[ax].spaceDirection);
1239 nrrdSpaceVecScale(vector, 1.0/(*spacing),
1240 nrrd->axis[ax].spaceDirection);
1241 } else {
1242 ret = nrrdSpacingStatusNone;
1243 *spacing = AIR_NAN(airFloatQNaN.f);
1244 nrrdSpaceVecSetNaN(vector);
1245 }
1246 }
1247 }
1248 return ret;
1249}
1250
1251int
1252nrrdOrientationReduce(Nrrd *nout, const Nrrd *nin,
1253 int setMinsFromOrigin) {
1254 static const char me[]="nrrdOrientationReduce";
1255 unsigned int spatialAxisNum, spatialAxisIdx[NRRD_DIM_MAX16], saxii;
1256 NrrdAxisInfo *axis;
1257
1258 if (!(nout && nin)) {
1259 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
1260 return 1;
1261 }
1262
1263 if (nout != nin) {
1264 if (nrrdCopy(nout, nin)) {
1265 biffAddf(NRRDnrrdBiffKey, "%s: trouble doing initial copying", me);
1266 return 1;
1267 }
1268 }
1269 if (!nout->spaceDim) {
1270 /* we're done! */
1271 return 0;
1272 }
1273 spatialAxisNum = nrrdSpatialAxesGet(nout, spatialAxisIdx);
1274 for (saxii=0; saxii<spatialAxisNum; saxii++) {
1275 axis = nout->axis + spatialAxisIdx[saxii];
1276 axis->spacing = nrrdSpaceVecNorm(nout->spaceDim,
1277 axis->spaceDirection);
1278 if (setMinsFromOrigin) {
1279 axis->min = (saxii < nout->spaceDim
1280 ? nout->spaceOrigin[saxii]
1281 : AIR_NAN(airFloatQNaN.f));
1282 }
1283 }
1284 nrrdSpaceSet(nout, nrrdSpaceUnknown);
1285
1286 return 0;
1287}
1288
1289/*
1290******** nrrdMetaData
1291**
1292** The brains of "unu dnorm" (for Diderot normalization): put all meta-data
1293** of a nrrd into some simpler canonical form.
1294**
1295** This function probably doesn't belong in this file, but it is kind
1296** the opposite of nrrdOrientationReduce (above), so here it is
1297*/
1298int
1299nrrdMetaDataNormalize(Nrrd *nout, const Nrrd *nin,
1300 int version,
1301 int trivialOrient,
1302 int permuteComponentAxisFastest,
1303 int recenterGrid,
1304 double sampleSpacing,
1305 int *lostMeasurementFrame) {
1306 static const char me[]="nrrdMetaDataNormalize";
1307 size_t size[NRRD_DIM_MAX16];
1308 int kindIn, kindOut, haveMM, gotmf;
1309 unsigned int kindAxis, axi, si, sj;
1310 Nrrd *ntmp;
1311 airArray *mop;
1312
1313 if (!(nout && nin)) {
1314 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
1315 return 1;
1316 }
1317 if (airEnumValCheck(nrrdMetaDataCanonicalVersion, version)) {
1318 biffAddf(NRRDnrrdBiffKey, "%s: version %d not valid %s", me,
1319 version, nrrdMetaDataCanonicalVersion->name);
1320 return 1;
1321 }
1322 if (nrrdMetaDataCanonicalVersionAlpha != version) {
1323 biffAddf(NRRDnrrdBiffKey, "%s: sorry, %s %s not implemented (only %s)", me,
1324 nrrdMetaDataCanonicalVersion->name,
1325 airEnumStr(nrrdMetaDataCanonicalVersion, version),
1326 airEnumStr(nrrdMetaDataCanonicalVersion,
1327 nrrdMetaDataCanonicalVersionAlpha));
1328 return 1;
1329 }
1330
1331 if (_nrrdCheck(nin, AIR_FALSE0 /* checkData */, AIR_TRUE1 /* useBiff */)) {
1332 biffAddf(NRRDnrrdBiffKey, "%s: basic check failed", me);
1333 return 1;
1334 }
1335 /* but can't deal with block type */
1336 if (nrrdTypeBlock == nin->type) {
1337 biffAddf(NRRDnrrdBiffKey, "%s: can only have scalar types (not %s)", me,
1338 airEnumStr(nrrdType, nrrdTypeBlock));
1339 return 1;
1340 }
1341
1342 /* look at all per-axis kinds */
1343 /* see if there's a range kind, verify that there's only one */
1344 /* set haveMM */
1345 haveMM = AIR_TRUE1;
1346 kindIn = nrrdKindUnknown;
1347 kindAxis = 0; /* only means something if kindIn != nrrdKindUnknown */
1348 for (axi=0; axi<nin->dim; axi++) {
1349 if (nrrdKindUnknown == nin->axis[axi].kind
1350 || nrrdKindIsDomain(nin->axis[axi].kind)) {
1351 haveMM &= AIR_EXISTS(nin->axis[axi].min)(((int)(!((nin->axis[axi].min) - (nin->axis[axi].min)))
));
1352 haveMM &= AIR_EXISTS(nin->axis[axi].max)(((int)(!((nin->axis[axi].max) - (nin->axis[axi].max)))
));
1353 } else {
1354 if (nrrdKindUnknown != kindIn) {
1355 biffAddf(NRRDnrrdBiffKey, "%s: got non-domain kind %s on axis %u, but already "
1356 "have kind %s on previous axis %u", me,
1357 airEnumStr(nrrdKind, nin->axis[axi].kind), axi,
1358 airEnumStr(nrrdKind, kindIn), kindAxis);
1359 return 1;
1360 }
1361 kindIn = nin->axis[axi].kind;
1362 kindAxis = axi;
1363 }
1364 }
1365
1366 if (nrrdKindUnknown != kindIn && kindAxis) {
1367 /* have a non-domain axis, and it isn't the fastest */
1368 if (permuteComponentAxisFastest) {
1369 if (nout == nin) {
1370 biffAddf(NRRDnrrdBiffKey, "%s: can't permute non-domain axis %u (kind %s) "
1371 "to axis 0 with nout == nin", me,
1372 kindAxis, airEnumStr(nrrdKind, kindIn));
1373 return 1;
1374 }
1375 biffAddf(NRRDnrrdBiffKey, "%s: sorry, permuting non-domain axis %u (kind %s) "
1376 "to axis 0 not yet implemented", me,
1377 kindAxis, airEnumStr(nrrdKind, kindIn));
1378 return 1;
1379 } else {
1380 /* caller thinks its okay for non-domain axis to be on
1381 something other than fastest axis */
1382 if (nrrdMetaDataCanonicalVersionAlpha == version) {
1383 biffAddf(NRRDnrrdBiffKey, "%s: (%s) non-domain axis %u (kind %s) "
1384 "must be fastest axis", me,
1385 airEnumStr(nrrdMetaDataCanonicalVersion, version),
1386 kindAxis, airEnumStr(nrrdKind, kindIn));
1387 return 1;
1388 }
1389 /* maybe with nrrdMetaDataCanonicalVersionAlpha != version
1390 it is okay to have non-domain axis on non-fastest axis? */
1391 }
1392 }
1393
1394 /* HEY: would be nice to handle a stub "scalar" axis by deleting it */
1395
1396 /* see if the non-domain kind is something we can interpret as a tensor */
1397 if (nrrdKindUnknown != kindIn) {
1398 switch (kindIn) {
1399 /* ======= THESE are the kinds that we can possibly output ======= */
1400 case nrrdKind2Vector:
1401 case nrrdKind3Vector:
1402 case nrrdKind4Vector:
1403 case nrrdKind2DSymMatrix:
1404 case nrrdKind2DMatrix:
1405 case nrrdKind3DSymMatrix:
1406 case nrrdKind3DMatrix:
1407 /* =============================================================== */
1408 kindOut = kindIn;
1409 break;
1410 /* Some other kinds are mapped to those above */
1411 case nrrdKind3Color:
1412 case nrrdKindRGBColor:
1413 kindOut = nrrdKind3Vector;
1414 break;
1415 case nrrdKind4Color:
1416 case nrrdKindRGBAColor:
1417 kindOut = nrrdKind4Vector;
1418 break;
1419 default:
1420 biffAddf(NRRDnrrdBiffKey, "%s: got non-conforming kind %s on axis %u", me,
1421 airEnumStr(nrrdKind, kindIn), kindAxis);
1422 return 1;
1423 }
1424 } else {
1425 /* kindIn is nrrdKindUnknown, so its a simple scalar image,
1426 and that's what the output will be too; kindOut == nrrdKindUnknown
1427 is used in the code below to say "its a scalar image" */
1428 kindOut = nrrdKindUnknown;
1429 }
1430
1431 /* initialize output by copying meta-data from nin to ntmp */
1432 mop = airMopNew();
1433 ntmp = nrrdNew();
1434 airMopAdd(mop, ntmp, (airMopper)nrrdNix, airMopAlways);
1435 /* HEY this is doing the work of a shallow copy, which isn't
1436 available in the API. You can pass nrrdCopy() a nin with NULL
1437 nin->data, which implements a shallow copy, but we can't set
1438 nin->data=NULL here because of const correctness */
1439 nrrdAxisInfoGet_nva(nin, nrrdAxisInfoSize, size);
1440 if (nrrdWrap_nva(ntmp, NULL((void*)0), nin->type, nin->dim, size)) {
1441 biffAddf(NRRDnrrdBiffKey, "%s: couldn't wrap buffer nrrd around NULL", me);
1442 airMopError(mop); return 1;
1443 }
1444 /* so ntmp->data == NULL */
1445 nrrdAxisInfoCopy(ntmp, nin, NULL((void*)0), NRRD_AXIS_INFO_SIZE_BIT(1<< 1));
1446 if (nrrdBasicInfoCopy(ntmp, nin, NRRD_BASIC_INFO_DATA_BIT(1<< 1))) {
1447 biffAddf(NRRDnrrdBiffKey, "%s: trouble copying basic info", me);
1448 airMopError(mop); return 1;
1449 }
1450
1451 /* no comments */
1452 nrrdCommentClear(ntmp);
1453
1454 /* no measurement frame */
1455 gotmf = AIR_FALSE0;
1456 for (si=0; si<NRRD_SPACE_DIM_MAX8; si++) {
1457 for (sj=0; sj<NRRD_SPACE_DIM_MAX8; sj++) {
1458 gotmf |= AIR_EXISTS(ntmp->measurementFrame[si][sj])(((int)(!((ntmp->measurementFrame[si][sj]) - (ntmp->measurementFrame
[si][sj])))));
1459 }
1460 }
1461 if (lostMeasurementFrame) {
1462 *lostMeasurementFrame = gotmf;
1463 }
1464 for (si=0; si<NRRD_SPACE_DIM_MAX8; si++) {
1465 for (sj=0; sj<NRRD_SPACE_DIM_MAX8; sj++) {
1466 ntmp->measurementFrame[si][sj] = AIR_NAN(airFloatQNaN.f);
1467 }
1468 }
1469
1470 /* no key/value pairs */
1471 nrrdKeyValueClear(ntmp);
1472
1473 /* no content field */
1474 ntmp->content = airFree(ntmp->content);
1475
1476 /* normalize domain kinds to "space" */
1477 /* HEY: if Diderot supports time-varying fields, this will have to change */
1478 /* turn off centers (current Diderot semantics don't expose centering) */
1479 /* turn off thickness */
1480 /* turn off labels and units */
1481 for (axi=0; axi<ntmp->dim; axi++) {
1482 if (nrrdKindUnknown == kindOut) {
1483 ntmp->axis[axi].kind = nrrdKindSpace;
1484 } else {
1485 ntmp->axis[axi].kind = (kindAxis == axi
1486 ? kindOut
1487 : nrrdKindSpace);
1488 }
1489 ntmp->axis[axi].center = nrrdCenterUnknown;
1490 ntmp->axis[axi].thickness = AIR_NAN(airFloatQNaN.f);
1491 ntmp->axis[axi].label = airFree(ntmp->axis[axi].label);
1492 ntmp->axis[axi].units = airFree(ntmp->axis[axi].units);
1493 ntmp->axis[axi].min = AIR_NAN(airFloatQNaN.f);
1494 ntmp->axis[axi].max = AIR_NAN(airFloatQNaN.f);
1495 ntmp->axis[axi].spacing = AIR_NAN(airFloatQNaN.f);
1496 }
1497
1498 /* logic of orientation definition:
1499 If space dimension is known:
1500 set origin to zero if not already set
1501 set space direction to unit vector if not already set
1502 Else if have per-axis min and max:
1503 set spae origin and directions to communicate same intent
1504 as original per-axis min and max and original centering
1505 Else
1506 set origin to zero and all space directions to units.
1507 (It might be nice to use gage's logic for mapping from world to index,
1508 but we have to accept a greater variety of kinds and dimensions
1509 than gage ever has to process.)
1510 The result is that space origin and space directions are set.
1511 the "space" field is not used, only "spaceDim"
1512 */
1513 /* no named space */
1514 ntmp->space = nrrdSpaceUnknown;
1515 if (ntmp->spaceDim && !trivialOrient) {
1516 int saxi = 0;
1517 if (!nrrdSpaceVecExists(ntmp->spaceDim, ntmp->spaceOrigin)) {
1518 nrrdSpaceVecSetZero(ntmp->spaceOrigin);
1519 }
1520 for (axi=0; axi<ntmp->dim; axi++) {
1521 if (nrrdKindUnknown == kindOut || kindAxis != axi) {
1522 /* its a domain axis of output */
1523 if (!nrrdSpaceVecExists(ntmp->spaceDim,
1524 ntmp->axis[axi].spaceDirection)) {
1525 nrrdSpaceVecSetZero(ntmp->axis[axi].spaceDirection);
1526 ntmp->axis[axi].spaceDirection[saxi] = sampleSpacing;
1527 }
1528 /* else we leave existing space vector as is */
1529 saxi++;
1530 } else {
1531 /* else its a range (non-domain, component) axis */
1532 nrrdSpaceVecSetNaN(ntmp->axis[axi].spaceDirection);
1533 }
1534 }
1535 } else if (haveMM && !trivialOrient) {
1536 int saxi = 0;
1537 size_t N;
1538 double rng;
1539 for (axi=0; axi<ntmp->dim; axi++) {
1540 if (nrrdKindUnknown == kindOut || kindAxis != axi) {
1541 /* its a domain axis of output */
1542 nrrdSpaceVecSetZero(ntmp->axis[axi].spaceDirection);
1543 rng = nin->axis[axi].max - nin->axis[axi].min;
1544 if (nrrdCenterNode == nin->axis[axi].center) {
1545 ntmp->spaceOrigin[saxi] = nin->axis[axi].min;
1546 N = nin->axis[axi].size;
1547 ntmp->axis[axi].spaceDirection[saxi] = rng/(N-1);
1548 } else {
1549 /* unknown centering treated as cell */
1550 N = nin->axis[axi].size;
1551 ntmp->spaceOrigin[saxi] = nin->axis[axi].min + (rng/N)/2;
1552 ntmp->axis[axi].spaceDirection[saxi] = rng/N;
1553 }
1554 saxi++;
1555 } else {
1556 /* else its a range axis */
1557 nrrdSpaceVecSetNaN(ntmp->axis[axi].spaceDirection);
1558 }
1559 }
1560 ntmp->spaceDim = saxi;
1561 } else {
1562 /* either trivialOrient, or, not spaceDim and not haveMM */
1563 int saxi = 0;
1564 nrrdSpaceVecSetZero(ntmp->spaceOrigin);
1565 for (axi=0; axi<ntmp->dim; axi++) {
1566 if (nrrdKindUnknown == kindOut || kindAxis != axi) {
1567 /* its a domain axis of output */
1568 nrrdSpaceVecSetZero(ntmp->axis[axi].spaceDirection);
1569 ntmp->axis[axi].spaceDirection[saxi]
1570 = (AIR_EXISTS(nin->axis[axi].spacing)(((int)(!((nin->axis[axi].spacing) - (nin->axis[axi].spacing
)))))
1571 ? nin->axis[axi].spacing
1572 : sampleSpacing);
1573 saxi++;
1574 } else {
1575 /* else its a range axis */
1576 nrrdSpaceVecSetNaN(ntmp->axis[axi].spaceDirection);
1577 }
1578 }
1579 ntmp->spaceDim = saxi;
1580 }
1581
1582 /* space dimension has to match the number of domain axes */
1583 if (ntmp->dim != ntmp->spaceDim + !!kindOut) {
1584 biffAddf(NRRDnrrdBiffKey, "%s: output dim %d != spaceDim %d + %d %s%s%s%s",
1585 me, ntmp->dim, ntmp->spaceDim, !!kindOut,
1586 kindOut ? "for non-scalar (" : "(scalar data)",
1587 kindOut ? airEnumStr(nrrdKind, kindOut) : "",
1588 kindOut ? ") data" : "",
1589 kindOut ? "" : "; a non-domain axis in the input "
1590 "may be missing an informative \"kind\", leading to the "
1591 "false assumption of a scalar array");
1592 airMopError(mop); return 1;
1593 }
1594
1595 if (recenterGrid) {
1596 /* sets field's origin so field is centered on the origin. capiche? */
1597 /* this code was tacked on later than the stuff above, so its
1598 logic could probably be moved up there, but it seems cleaner to
1599 have it as a separate post-process */
1600 double mean[NRRD_SPACE_DIM_MAX8];
1601 nrrdSpaceVecSetZero(mean);
1602 for (axi=0; axi<ntmp->dim; axi++) {
1603 if (nrrdKindUnknown == kindOut || kindAxis != axi) {
1604 nrrdSpaceVecScaleAdd2(mean, 1.0, mean,
1605 0.5*(ntmp->axis[axi].size - 1),
1606 ntmp->axis[axi].spaceDirection);
1607 }
1608 }
1609 nrrdSpaceVecScaleAdd2(mean, 1.0, mean,
1610 1.0, ntmp->spaceOrigin);
1611 /* now mean is the center of the field */
1612 nrrdSpaceVecScaleAdd2(ntmp->spaceOrigin,
1613 1.0, ntmp->spaceOrigin,
1614 -1.0, mean);
1615 }
1616
1617 /* with that all done, now copy from ntmp to nout */
1618 if (nout != nin) {
1619 /* have to copy data */
1620 ntmp->data = nin->data;
1621 if (nrrdCopy(nout, ntmp)) {
1622 biffAddf(NRRDnrrdBiffKey, "%s: problem copying (with data) to output", me);
1623 airMopError(mop); return 1;
1624 }
1625 } else {
1626 /* nout == nin; have to copy only meta-data, leave data as is */
1627 void *data = nin->data;
1628 /* ntmp->data == NULL, so this is a shallow copy */
1629 if (nrrdCopy(nout, ntmp)) {
1630 biffAddf(NRRDnrrdBiffKey, "%s: problem copying meta-data to output", me);
1631 airMopError(mop); return 1;
1632 }
1633 nout->data = data;
1634 }
1635
1636 airMopOkay(mop);
1637 return 0;
1638}
1639