Bug Summary

File:src/nrrd/subset.c
Location:line 246, column 22
Description:The left operand of '*' is a garbage value

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******** nrrdSlice()
29**
30** slices a nrrd along a given axis, at a given position.
31**
32** This is a newer version of the procedure, which is simpler, faster,
33** and requires less memory overhead than the first one. It is based
34** on the observation that any slice is a periodic square-wave pattern
35** in the original data (viewed as a one- dimensional array). The
36** characteristics of that periodic pattern are how far from the
37** beginning it starts (offset), the length of the "on" part (length),
38** the period (period), and the number of periods (numper).
39*/
40int
41nrrdSlice(Nrrd *nout, const Nrrd *cnin, unsigned int saxi, size_t pos) {
42 static const char me[]="nrrdSlice", func[]="slice";
43 size_t
44 I,
45 rowLen, /* length of segment */
46 colStep, /* distance between start of each segment */
47 colLen, /* number of periods */
48 szOut[NRRD_DIM_MAX16];
49 unsigned int ai, outdim;
50 int map[NRRD_DIM_MAX16];
51 const char *src;
52 char *dest, stmp[2][AIR_STRLEN_SMALL(128+1)];
53 airArray *mop;
54 Nrrd *nin;
55
56 if (!(cnin && nout)) {
57 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
58 return 1;
59 }
60 if (nout == cnin) {
61 biffAddf(NRRDnrrdBiffKey, "%s: nout==nin disallowed", me);
62 return 1;
63 }
64 if (1 == cnin->dim) {
65 if (0 != saxi) {
66 biffAddf(NRRDnrrdBiffKey, "%s: slice axis must be 0, not %u, for 1-D array",
67 me, saxi);
68 return 1;
69 }
70 } else {
71 if (!( saxi < cnin->dim )) {
72 biffAddf(NRRDnrrdBiffKey, "%s: slice axis %d out of bounds (0 to %d)",
73 me, saxi, cnin->dim-1);
74 return 1;
75 }
76 }
77 if (!( pos < cnin->axis[saxi].size )) {
78 biffAddf(NRRDnrrdBiffKey, "%s: position %s out of bounds (0 to %s)", me,
79 airSprintSize_t(stmp[0], pos),
80 airSprintSize_t(stmp[1], cnin->axis[saxi].size-1));
81 return 1;
82 }
83 /* this shouldn't actually be necessary .. */
84 if (!nrrdElementSize(cnin)) {
85 biffAddf(NRRDnrrdBiffKey, "%s: nrrd reports zero element size!", me);
86 return 1;
87 }
88
89 /* HEY: copy and paste from measure.c/nrrdProject */
90 mop = airMopNew();
91 if (1 == cnin->dim) {
92 /* There are more efficient ways of dealing with this case; this way is
93 easy to implement because it leaves most of the established code below
94 only superficially changed; uniformly replacing nin with (nin ? nin :
95 cnin), even if pointlessly so; this expression that can't be assigned
96 to a new variable because of the difference in const. */
97 nin = nrrdNew();
98 airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopAlways);
99 if (nrrdAxesInsert(nin, cnin, 1)) {
100 biffAddf(NRRDnrrdBiffKey, "%s: trouble inserting axis on 1-D array", me);
101 airMopError(mop); return 1;
102 }
103 } else {
104 nin = NULL((void*)0);
105 }
106
107 /* set up control variables */
108 rowLen = colLen = 1;
109 for (ai=0; ai<(nin ? nin : cnin)->dim; ai++) {
110 if (ai < saxi) {
111 rowLen *= (nin ? nin : cnin)->axis[ai].size;
112 } else if (ai > saxi) {
113 colLen *= (nin ? nin : cnin)->axis[ai].size;
114 }
115 }
116 rowLen *= nrrdElementSize(nin ? nin : cnin);
117 colStep = rowLen*(nin ? nin : cnin)->axis[saxi].size;
118
119 outdim = (nin ? nin : cnin)->dim-1;
120 for (ai=0; ai<outdim; ai++) {
121 map[ai] = AIR_INT(ai)((int)(ai)) + (ai >= saxi);
122 szOut[ai] = (nin ? nin : cnin)->axis[map[ai]].size;
123 }
124 nout->blockSize = (nin ? nin : cnin)->blockSize;
125 if (nrrdMaybeAlloc_nva(nout, (nin ? nin : cnin)->type, outdim, szOut)) {
126 biffAddf(NRRDnrrdBiffKey, "%s: failed to create slice", me);
127 airMopError(mop); return 1;
128 }
129
130 /* the skinny */
131 src = AIR_CAST(const char *, (nin ? nin : cnin)->data)((const char *)((nin ? nin : cnin)->data));
132 dest = AIR_CAST(char *, nout->data)((char *)(nout->data));
133 src += rowLen*pos;
134 for (I=0; I<colLen; I++) {
135 /* HEY: replace with AIR_MEMCPY() or similar, when applicable */
136 memcpy(dest, src, rowLen)__builtin___memcpy_chk (dest, src, rowLen, __builtin_object_size
(dest, 0));
137 src += colStep;
138 dest += rowLen;
139 }
140
141 /* copy the peripheral information */
142 if (nrrdAxisInfoCopy(nout, (nin ? nin : cnin), map, NRRD_AXIS_INFO_NONE0)) {
143 biffAddf(NRRDnrrdBiffKey, "%s:", me);
144 airMopError(mop); return 1;
145 }
146 if (nrrdContentSet_va(nout, func, cnin /* hide possible axinsert*/,
147 "%d,%d", saxi, pos)) {
148 biffAddf(NRRDnrrdBiffKey, "%s:", me);
149 airMopError(mop); return 1;
150 }
151 if (nrrdBasicInfoCopy(nout, (nin ? nin : cnin),
152 NRRD_BASIC_INFO_DATA_BIT(1<< 1)
153 | NRRD_BASIC_INFO_TYPE_BIT(1<< 2)
154 | NRRD_BASIC_INFO_BLOCKSIZE_BIT(1<< 3)
155 | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4)
156 | NRRD_BASIC_INFO_SPACEORIGIN_BIT(1<<10)
157 | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5)
158 | NRRD_BASIC_INFO_COMMENTS_BIT(1<<14)
159 | (nrrdStateKeyValuePairsPropagate
160 ? 0
161 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT(1<<15)))) {
162 biffAddf(NRRDnrrdBiffKey, "%s:", me);
163 airMopError(mop); return 1;
164 }
165 /* translate origin if this was a spatial axis, otherwise copy */
166 /* note that if there is no spatial info at all, this is all harmless */
167 if (AIR_EXISTS((nin ? nin : cnin)->axis[saxi].spaceDirection[0])(((int)(!(((nin ? nin : cnin)->axis[saxi].spaceDirection[0
]) - ((nin ? nin : cnin)->axis[saxi].spaceDirection[0]))))
)) {
168 nrrdSpaceVecScaleAdd2(nout->spaceOrigin,
169 1.0, (nin ? nin : cnin)->spaceOrigin,
170 AIR_CAST(double, pos)((double)(pos)),
171 (nin ? nin : cnin)->axis[saxi].spaceDirection);
172 } else {
173 nrrdSpaceVecCopy(nout->spaceOrigin, (nin ? nin : cnin)->spaceOrigin);
174 }
175 airMopOkay(mop);
176 return 0;
177}
178
179/*
180******** nrrdCrop()
181**
182** select some sub-volume inside a given nrrd, producing an output
183** nrrd with the same dimensions, but with equal or smaller sizes
184** along each axis.
185*/
186int
187nrrdCrop(Nrrd *nout, const Nrrd *nin, size_t *min, size_t *max) {
188 static const char me[]="nrrdCrop", func[] = "crop";
189 char buff1[NRRD_DIM_MAX16*30], buff2[AIR_STRLEN_SMALL(128+1)];
190 unsigned int ai;
191 size_t I,
192 lineSize, /* #bytes in one scanline to be copied */
193 typeSize, /* size of data type */
194 cIn[NRRD_DIM_MAX16], /* coords for line start, in input */
195 cOut[NRRD_DIM_MAX16], /* coords for line start, in output */
196 szIn[NRRD_DIM_MAX16],
197 szOut[NRRD_DIM_MAX16],
198 idxIn, idxOut, /* linear indices for input and output */
199 numLines; /* number of scanlines in output nrrd */
200 char *dataIn, *dataOut, stmp[3][AIR_STRLEN_SMALL(128+1)];
201
202 /* errors */
203 if (!(nout && nin && min && max)) {
1
Taking false branch
→
204 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
205 return 1;
206 }
207 if (nout == nin) {
2
←
Assuming 'nout' is not equal to 'nin'
→
3
←
Taking false branch
→
208 biffAddf(NRRDnrrdBiffKey, "%s: nout==nin disallowed", me);
209 return 1;
210 }
211 for (ai=0; ai<nin->dim; ai++) {
4
←
Loop condition is false. Execution continues on line 227
→
212 if (!(min[ai] <= max[ai])) {
213 biffAddf(NRRDnrrdBiffKey, "%s: axis %d min (%s) not <= max (%s)", me, ai,
214 airSprintSize_t(stmp[0], min[ai]),
215 airSprintSize_t(stmp[1], max[ai]));
216 return 1;
217 }
218 if (!( min[ai] < nin->axis[ai].size && max[ai] < nin->axis[ai].size )) {
219 biffAddf(NRRDnrrdBiffKey, "%s: axis %d min (%s) or max (%s) out of bounds [0,%s]",
220 me, ai, airSprintSize_t(stmp[0], min[ai]),
221 airSprintSize_t(stmp[1], max[ai]),
222 airSprintSize_t(stmp[2], nin->axis[ai].size-1));
223 return 1;
224 }
225 }
226 /* this shouldn't actually be necessary .. */
227 if (!nrrdElementSize(nin)) {
5
←
Taking false branch
→
228 biffAddf(NRRDnrrdBiffKey, "%s: nrrd reports zero element size!", me);
229 return 1;
230 }
231
232 /* allocate */
233 nrrdAxisInfoGet_nva(nin, nrrdAxisInfoSize, szIn);
234 numLines = 1;
235 for (ai=0; ai<nin->dim; ai++) {
6
←
Loop condition is false. Execution continues on line 241
→
236 szOut[ai] = max[ai] - min[ai] + 1;
237 if (ai) {
238 numLines *= szOut[ai];
239 }
240 }
241 nout->blockSize = nin->blockSize;
242 if (nrrdMaybeAlloc_nva(nout, nin->type, nin->dim, szOut)) {
7
←
Taking false branch
→
243 biffAddf(NRRDnrrdBiffKey, "%s:", me);
244 return 1;
245 }
246 lineSize = szOut[0]*nrrdElementSize(nin);
8
←
The left operand of '*' is a garbage value
247
248 /* the skinny */
249 typeSize = nrrdElementSize(nin);
250 dataIn = (char *)nin->data;
251 dataOut = (char *)nout->data;
252 memset(cOut, 0, NRRD_DIM_MAX*sizeof(*cOut))__builtin___memset_chk (cOut, 0, 16*sizeof(*cOut), __builtin_object_size
(cOut, 0));
253 /*
254 printf("!%s: nin->dim = %d\n", me, nin->dim);
255 printf("!%s: min = %d %d %d\n", me, min[0], min[1], min[2]);
256 printf("!%s: szIn = %d %d %d\n", me, szIn[0], szIn[1], szIn[2]);
257 printf("!%s: szOut = %d %d %d\n", me, szOut[0], szOut[1], szOut[2]);
258 printf("!%s: lineSize = %d\n", me, lineSize);
259 printf("!%s: typeSize = %d\n", me, typeSize);
260 printf("!%s: numLines = %d\n", me, (int)numLines);
261 */
262 for (I=0; I<numLines; I++) {
263 for (ai=0; ai<nin->dim; ai++) {
264 cIn[ai] = cOut[ai] + min[ai];
265 }
266 NRRD_INDEX_GEN(idxOut, cOut, szOut, nin->dim){ unsigned int ddd = (nin->dim); (idxOut) = 0; while (ddd)
{ ddd--; (idxOut) = (cOut)[ddd] + (szOut)[ddd]*(idxOut); } };
267 NRRD_INDEX_GEN(idxIn, cIn, szIn, nin->dim){ unsigned int ddd = (nin->dim); (idxIn) = 0; while (ddd) {
ddd--; (idxIn) = (cIn)[ddd] + (szIn)[ddd]*(idxIn); } };
268 /*
269 printf("!%s: %5d: cOut=(%3d,%3d,%3d) --> idxOut = %5d\n",
270 me, (int)I, cOut[0], cOut[1], cOut[2], (int)idxOut);
271 printf("!%s: %5d: cIn=(%3d,%3d,%3d) --> idxIn = %5d\n",
272 me, (int)I, cIn[0], cIn[1], cIn[2], (int)idxIn);
273 */
274 memcpy(dataOut + idxOut*typeSize, dataIn + idxIn*typeSize, lineSize)__builtin___memcpy_chk (dataOut + idxOut*typeSize, dataIn + idxIn
*typeSize, lineSize, __builtin_object_size (dataOut + idxOut*
typeSize, 0));
275 /* the lowest coordinate in cOut[] will stay zero, since we are
276 copying one (1-D) scanline at a time */
277 NRRD_COORD_INCR(cOut, szOut, nin->dim, 1)if ((1) < (nin->dim)) { (cOut)[(1)]++; { unsigned int ddd
; for (ddd=0; ddd+1 < ((nin->dim)-(1)) && ((cOut
)+(1))[ddd] >= ((szOut)+(1))[ddd]; ddd++) { ((cOut)+(1))[ddd
] = 0; ((cOut)+(1))[ddd+1]++; } if ((nin->dim)-(1)) { ((cOut
)+(1))[((nin->dim)-(1))-1] = ((((cOut)+(1))[((nin->dim)
-(1))-1]) < (((szOut)+(1))[((nin->dim)-(1))-1]-1) ? (((
cOut)+(1))[((nin->dim)-(1))-1]) : (((szOut)+(1))[((nin->
dim)-(1))-1]-1)); } }; };
278 }
279 if (nrrdAxisInfoCopy(nout, nin, NULL((void*)0), (NRRD_AXIS_INFO_SIZE_BIT(1<< 1) |
280 NRRD_AXIS_INFO_MIN_BIT(1<< 4) |
281 NRRD_AXIS_INFO_MAX_BIT(1<< 5) ))) {
282 biffAddf(NRRDnrrdBiffKey, "%s:", me);
283 return 1;
284 }
285 for (ai=0; ai<nin->dim; ai++) {
286 nrrdAxisInfoPosRange(&(nout->axis[ai].min), &(nout->axis[ai].max),
287 nin, ai, AIR_CAST(double, min[ai])((double)(min[ai])),
288 AIR_CAST(double, max[ai])((double)(max[ai])));
289 /* do the safe thing first */
290 nout->axis[ai].kind = _nrrdKindAltered(nin->axis[ai].kind, AIR_FALSE0);
291 /* try cleverness */
292 if (!nrrdStateKindNoop) {
293 if (nout->axis[ai].size == nin->axis[ai].size) {
294 /* we can safely copy kind; the samples didn't change */
295 nout->axis[ai].kind = nin->axis[ai].kind;
296 } else if (nrrdKind4Color == nin->axis[ai].kind
297 && 3 == szOut[ai]) {
298 nout->axis[ai].kind = nrrdKind3Color;
299 } else if (nrrdKind4Vector == nin->axis[ai].kind
300 && 3 == szOut[ai]) {
301 nout->axis[ai].kind = nrrdKind3Vector;
302 } else if ((nrrdKind4Vector == nin->axis[ai].kind
303 || nrrdKind3Vector == nin->axis[ai].kind)
304 && 2 == szOut[ai]) {
305 nout->axis[ai].kind = nrrdKind2Vector;
306 } else if (nrrdKindRGBAColor == nin->axis[ai].kind
307 && 0 == min[ai]
308 && 2 == max[ai]) {
309 nout->axis[ai].kind = nrrdKindRGBColor;
310 } else if (nrrdKind2DMaskedSymMatrix == nin->axis[ai].kind
311 && 1 == min[ai]
312 && max[ai] == szIn[ai]-1) {
313 nout->axis[ai].kind = nrrdKind2DSymMatrix;
314 } else if (nrrdKind2DMaskedMatrix == nin->axis[ai].kind
315 && 1 == min[ai]
316 && max[ai] == szIn[ai]-1) {
317 nout->axis[ai].kind = nrrdKind2DMatrix;
318 } else if (nrrdKind3DMaskedSymMatrix == nin->axis[ai].kind
319 && 1 == min[ai]
320 && max[ai] == szIn[ai]-1) {
321 nout->axis[ai].kind = nrrdKind3DSymMatrix;
322 } else if (nrrdKind3DMaskedMatrix == nin->axis[ai].kind
323 && 1 == min[ai]
324 && max[ai] == szIn[ai]-1) {
325 nout->axis[ai].kind = nrrdKind3DMatrix;
326 }
327 }
328 }
329 strcpy(buff1, "")__builtin___strcpy_chk (buff1, "", __builtin_object_size (buff1
, 2 > 1 ? 1 : 0));
330 for (ai=0; ai<nin->dim; ai++) {
331 sprintf(buff2, "%s[%s,%s]", (ai ? "x" : ""),__builtin___sprintf_chk (buff2, 0, __builtin_object_size (buff2
, 2 > 1 ? 1 : 0), "%s[%s,%s]", (ai ? "x" : ""), airSprintSize_t
(stmp[0], min[ai]), airSprintSize_t(stmp[1], max[ai]))
332 airSprintSize_t(stmp[0], min[ai]),__builtin___sprintf_chk (buff2, 0, __builtin_object_size (buff2
, 2 > 1 ? 1 : 0), "%s[%s,%s]", (ai ? "x" : ""), airSprintSize_t
(stmp[0], min[ai]), airSprintSize_t(stmp[1], max[ai]))
333 airSprintSize_t(stmp[1], max[ai]))__builtin___sprintf_chk (buff2, 0, __builtin_object_size (buff2
, 2 > 1 ? 1 : 0), "%s[%s,%s]", (ai ? "x" : ""), airSprintSize_t
(stmp[0], min[ai]), airSprintSize_t(stmp[1], max[ai]));
334 strcat(buff1, buff2)__builtin___strcat_chk (buff1, buff2, __builtin_object_size (
buff1, 2 > 1 ? 1 : 0));
335 }
336 if (nrrdContentSet_va(nout, func, nin, "%s", buff1)) {
337 biffAddf(NRRDnrrdBiffKey, "%s:", me);
338 return 1;
339 }
340 if (nrrdBasicInfoCopy(nout, nin,
341 NRRD_BASIC_INFO_DATA_BIT(1<< 1)
342 | NRRD_BASIC_INFO_TYPE_BIT(1<< 2)
343 | NRRD_BASIC_INFO_BLOCKSIZE_BIT(1<< 3)
344 | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4)
345 | NRRD_BASIC_INFO_SPACEORIGIN_BIT(1<<10)
346 | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5)
347 | NRRD_BASIC_INFO_COMMENTS_BIT(1<<14)
348 | (nrrdStateKeyValuePairsPropagate
349 ? 0
350 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT(1<<15)))) {
351 biffAddf(NRRDnrrdBiffKey, "%s:", me);
352 return 1;
353 }
354 /* copy origin, then shift it along the spatial axes */
355 nrrdSpaceVecCopy(nout->spaceOrigin, nin->spaceOrigin);
356 for (ai=0; ai<nin->dim; ai++) {
357 if (AIR_EXISTS(nin->axis[ai].spaceDirection[0])(((int)(!((nin->axis[ai].spaceDirection[0]) - (nin->axis
[ai].spaceDirection[0])))))) {
358 nrrdSpaceVecScaleAdd2(nout->spaceOrigin,
359 1.0, nout->spaceOrigin,
360 AIR_CAST(double, min[ai])((double)(min[ai])),
361 nin->axis[ai].spaceDirection);
362 }
363 }
364
365
366 return 0;
367}
368
369/* ---- BEGIN non-NrrdIO */
370
371/*
372******** nrrdSliceSelect
373**
374** selects slices along axis "axi" from "nin", according to whether
375** line[i] is above or below thresh:
376**
377** line[i] >= thresh: slice i goes into noutAbove
378** line[i] < thresh: slice i goes into noutBelow
379**
380** Either noutAbove or noutBelow (but not both) can be passed
381** as NULL if you aren't interested in the output. It is a
382** biff-able error if the threshold is outside the range of
383** errors and a non-Null nrrd was passed for the correspondingly
384** empty output.
385*/
386int
387nrrdSliceSelect(Nrrd *noutAbove, Nrrd *noutBelow, const Nrrd *nin,
388 unsigned int saxi, Nrrd *_nline, double thresh) {
389 static const char me[]="nrrdSliceSelect";
390 airArray *mop;
391 Nrrd *nline, *nslice;
392 NrrdRange *rng;
393 double *line;
394 size_t II, LL, numAbove, numBelow, stride,
395 sizeAbove[NRRD_DIM_MAX16], sizeBelow[NRRD_DIM_MAX16];
396 unsigned int aa, bb;
397 int axmap[NRRD_DIM_MAX16];
398 char *above, *below, stmp[2][AIR_STRLEN_SMALL(128+1)];
399
400 if (!( (noutAbove || noutBelow) && nin && _nline )) {
401 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
402 return 1;
403 }
404 if (!AIR_EXISTS(thresh)(((int)(!((thresh) - (thresh)))))) {
405 biffAddf(NRRDnrrdBiffKey, "%s: thresh %g doesn't exist", me, thresh);
406 return 1;
407 }
408 if (!(saxi < nin->dim)) {
409 biffAddf(NRRDnrrdBiffKey, "%s: can't select axis %u of a %u-D input nrrd", me,
410 saxi, nin->dim);
411 return 1;
412 }
413 if (nrrdCheck(nin) || nrrdCheck(_nline)) {
414 biffAddf(NRRDnrrdBiffKey, "%s: basic problems with nin or nline", me);
415 return 1;
416 }
417 if (nrrdTypeBlock == nin->type) {
418 /* no good reason for this except that GLK forgets out to
419 set the blocksize in output */
420 biffAddf(NRRDnrrdBiffKey, "%s: sorry, can't handle type %s input", me,
421 airEnumStr(nrrdType, nrrdTypeBlock));
422 return 1;
423 }
424 if (!( nrrdTypeBlock != _nline->type
425 && 1 == _nline->dim)) {
426 biffAddf(NRRDnrrdBiffKey, "%s: nline must be 1-D array of scalars (not %u-D of %s)",
427 me, _nline->dim, airEnumStr(nrrdType, _nline->type));
428 return 1;
429 }
430 if (!( _nline->axis[0].size == nin->axis[saxi].size )) {
431 biffAddf(NRRDnrrdBiffKey, "%s: line length (%s) != axis[%u].size (%s)", me,
432 airSprintSize_t(stmp[0], _nline->axis[0].size), saxi,
433 airSprintSize_t(stmp[1], nin->axis[saxi].size));
434 return 1;
435 }
436 if (1 == nin->dim) {
437 biffAddf(NRRDnrrdBiffKey, "%s: sorry, slice-based implementation requires input "
438 "dimension > 1", me);
439 return 1;
440 }
441
442 mop = airMopNew();
443 rng = nrrdRangeNewSet(_nline, AIR_FALSE0);
444 airMopAdd(mop, rng, (airMopper)nrrdRangeNix, airMopAlways);
445 if (rng->hasNonExist) {
446 biffAddf(NRRDnrrdBiffKey, "%s: had non-existent values in line", me);
447 airMopError(mop); return 1;
448 }
449
450 nslice = nrrdNew();
451 airMopAdd(mop, nslice, (airMopper)nrrdNix, airMopAlways);
452 nline = nrrdNew();
453 airMopAdd(mop, nline, (airMopper)nrrdNuke, airMopAlways);
454 if (nrrdConvert(nline, _nline, nrrdTypeDouble)) {
455 biffAddf(NRRDnrrdBiffKey, "%s: problem copying line", me);
456 airMopError(mop); return 1;
457 }
458
459 line = AIR_CAST(double *, nline->data)((double *)(nline->data));
460 LL = nline->axis[0].size;
461 numAbove = numBelow = 0;
462 for (II=0; II<LL; II++) {
463 if (line[II] >= thresh) {
464 numAbove++;
465 } else {
466 numBelow++;
467 }
468 }
469 if (noutAbove && !numAbove) {
470 biffAddf(NRRDnrrdBiffKey, "%s: want slices for val >= thresh %g, "
471 "but highest value is %g < %g\n", me, thresh,
472 rng->max, thresh);
473 airMopError(mop); return 1;
474 }
475 if (noutBelow && !numBelow) {
476 biffAddf(NRRDnrrdBiffKey, "%s: want slices for val < thresh %g, "
477 "but lowest value is %g >= %g\n", me, thresh,
478 rng->min, thresh);
479 airMopError(mop); return 1;
480 }
481
482 nslice->dim = nin->dim-1;
483 nslice->type = nin->type;
484 bb = 0;
485 stride = nrrdElementSize(nin);
486 for (aa=0; aa<nin->dim; aa++) {
487 sizeAbove[aa] = sizeBelow[aa] = nin->axis[aa].size;
488 if (aa != saxi) {
489 axmap[aa] = aa;
490 nslice->axis[bb].size = nin->axis[aa].size;
491 if (aa < saxi) {
492 stride *= nin->axis[aa].size;
493 }
494 bb++;
495 } else {
496 axmap[aa] = -1;
497 }
498 }
499 sizeAbove[saxi] = numAbove;
500 sizeBelow[saxi] = numBelow;
501 if ((noutAbove
502 && nrrdMaybeAlloc_nva(noutAbove, nin->type, nin->dim, sizeAbove))
503 ||
504 (noutBelow
505 && nrrdMaybeAlloc_nva(noutBelow, nin->type, nin->dim, sizeBelow))) {
506 biffAddf(NRRDnrrdBiffKey, "%s: trouble allocating output", me);
507 airMopError(mop); return 1;
508 }
509 if (noutAbove) {
510 above = AIR_CAST(char *, noutAbove->data)((char *)(noutAbove->data));
511 } else {
512 above = NULL((void*)0);
513 }
514 if (noutBelow) {
515 below = AIR_CAST(char *, noutBelow->data)((char *)(noutBelow->data));
516 } else {
517 below = NULL((void*)0);
518 }
519
520 /* the skinny */
521 for (II=0; II<LL; II++) {
522 if (line[II] >= thresh) {
523 if (noutAbove) {
524 nslice->data = above;
525 if (nrrdSlice(nslice, nin, saxi, II)) {
526 biffAddf(NRRDnrrdBiffKey, "%s: trouble slicing (above) at %s", me,
527 airSprintSize_t(stmp[0], II));
528 airMopError(mop); return 1;
529 }
530 above += stride;
531 }
532 } else {
533 if (noutBelow) {
534 nslice->data = below;
535 if (nrrdSlice(nslice, nin, saxi, II)) {
536 biffAddf(NRRDnrrdBiffKey, "%s: trouble slicing (below) at %s", me,
537 airSprintSize_t(stmp[0], II));
538 airMopError(mop); return 1;
539 }
540 below += stride;
541 }
542 }
543 }
544
545 if (noutAbove) {
546 nrrdAxisInfoCopy(noutAbove, nin, axmap, NRRD_AXIS_INFO_NONE0);
547 if (nrrdBasicInfoCopy(noutAbove, nin,
548 NRRD_BASIC_INFO_DATA_BIT(1<< 1)
549 | NRRD_BASIC_INFO_TYPE_BIT(1<< 2)
550 | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4)
551 | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5)
552 | NRRD_BASIC_INFO_COMMENTS_BIT(1<<14)
553 | (nrrdStateKeyValuePairsPropagate
554 ? 0
555 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT(1<<15)))) {
556 biffAddf(NRRDnrrdBiffKey, "%s:", me);
557 return 1;
558 }
559 }
560 if (noutBelow) {
561 nrrdAxisInfoCopy(noutBelow, nin, axmap, NRRD_AXIS_INFO_NONE0);
562 if (nrrdBasicInfoCopy(noutBelow, nin,
563 NRRD_BASIC_INFO_DATA_BIT(1<< 1)
564 | NRRD_BASIC_INFO_TYPE_BIT(1<< 2)
565 | NRRD_BASIC_INFO_DIMENSION_BIT(1<< 4)
566 | NRRD_BASIC_INFO_CONTENT_BIT(1<< 5)
567 | NRRD_BASIC_INFO_COMMENTS_BIT(1<<14)
568 | (nrrdStateKeyValuePairsPropagate
569 ? 0
570 : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT(1<<15)))) {
571 biffAddf(NRRDnrrdBiffKey, "%s:", me);
572 return 1;
573 }
574 }
575
576 airMopOkay(mop);
577 return 0;
578}
579
580/*
581******** nrrdSample_nva()
582**
583** given coordinates within a nrrd, copies the
584** single element into given *val
585*/
586int
587nrrdSample_nva(void *val, const Nrrd *nrrd, const size_t *coord) {
588 static const char me[]="nrrdSample_nva";
589 size_t I, size[NRRD_DIM_MAX16], typeSize;
590 unsigned int ai;
591 char stmp[2][AIR_STRLEN_SMALL(128+1)];
592
593 if (!(nrrd && coord && val)) {
594 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
595 return 1;
596 }
597 /* this shouldn't actually be necessary .. */
598 if (!nrrdElementSize(nrrd)) {
599 biffAddf(NRRDnrrdBiffKey, "%s: nrrd reports zero element size!", me);
600 return 1;
601 }
602
603 typeSize = nrrdElementSize(nrrd);
604 nrrdAxisInfoGet_nva(nrrd, nrrdAxisInfoSize, size);
605 for (ai=0; ai<nrrd->dim; ai++) {
606 if (!( coord[ai] < size[ai] )) {
607 biffAddf(NRRDnrrdBiffKey, "%s: coordinate %s on axis %d out of bounds (0 to %s)",
608 me, airSprintSize_t(stmp[0], coord[ai]),
609 ai, airSprintSize_t(stmp[1], size[ai]-1));
610 return 1;
611 }
612 }
613
614 NRRD_INDEX_GEN(I, coord, size, nrrd->dim){ unsigned int ddd = (nrrd->dim); (I) = 0; while (ddd) { ddd
--; (I) = (coord)[ddd] + (size)[ddd]*(I); } };
615
616 memcpy(val, (char*)(nrrd->data) + I*typeSize, typeSize)__builtin___memcpy_chk (val, (char*)(nrrd->data) + I*typeSize
, typeSize, __builtin_object_size (val, 0));
617 return 0;
618}
619
620/*
621******** nrrdSample_va()
622**
623** var-args version of nrrdSample_nva()
624*/
625int
626nrrdSample_va(void *val, const Nrrd *nrrd, ...) {
627 static const char me[]="nrrdSample_va";
628 unsigned int ai;
629 size_t coord[NRRD_DIM_MAX16];
630 va_list ap;
631
632 if (!(nrrd && val)) {
633 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
634 return 1;
635 }
636
637 va_start(ap, nrrd)__builtin_va_start(ap, nrrd);
638 for (ai=0; ai<nrrd->dim; ai++) {
639 coord[ai] = va_arg(ap, size_t)__builtin_va_arg(ap, size_t);
640 }
641 va_end(ap)__builtin_va_end(ap);
642
643 if (nrrdSample_nva(val, nrrd, coord)) {
644 biffAddf(NRRDnrrdBiffKey, "%s:", me);
645 return 1;
646 }
647 return 0;
648}
649
650/*
651******** nrrdSimpleCrop()
652**
653*/
654int
655nrrdSimpleCrop(Nrrd *nout, const Nrrd *nin, unsigned int crop) {
656 static const char me[]="nrrdSimpleCrop";
657 unsigned int ai;
658 size_t min[NRRD_DIM_MAX16], max[NRRD_DIM_MAX16];
659
660 if (!(nout && nin)) {
661 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
662 return 1;
663 }
664 for (ai=0; ai<nin->dim; ai++) {
665 min[ai] = crop;
666 max[ai] = nin->axis[ai].size-1 - crop;
667 }
668 if (nrrdCrop(nout, nin, min, max)) {
669 biffAddf(NRRDnrrdBiffKey, "%s:", me);
670 return 1;
671 }
672 return 0;
673}
674
675int
676nrrdCropAuto(Nrrd *nout, const Nrrd *nin,
677 size_t _min[NRRD_DIM_MAX16], size_t _max[NRRD_DIM_MAX16],
678 const unsigned int *keep, unsigned int keepNum,
679 int measr, double frac, int offset) {
680 static const char me[]="nrrdCropAuto";
681 size_t min[NRRD_DIM_MAX16], max[NRRD_DIM_MAX16], NN, II;
682 int cropdo[NRRD_DIM_MAX16];
683 airArray *mop;
684 Nrrd *nperm, *nline;
685 unsigned int axi;
686 double *line;
687
688 if (!( nout && nin )) {
689 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
690 return 1;
691 }
692 if (keepNum && !keep) {
693 biffAddf(NRRDnrrdBiffKey, "%s: non-zero keepNum %u but NULL keep", me, keepNum);
694 return 1;
695 }
696 if (airEnumValCheck(nrrdMeasure, measr)) {
697 biffAddf(NRRDnrrdBiffKey, "%s: invalid %s measr %d", me,
698 nrrdMeasure->name, measr);
699 return 1;
700 }
701 if (!( AIR_EXISTS(frac)(((int)(!((frac) - (frac))))) && frac >= 0.0 && frac < 0.5 )) {
702 biffAddf(NRRDnrrdBiffKey, "%s: frac %g not in interval [0.0,0.5)", me, frac);
703 return 1;
704 }
705 for (axi=0; axi<nin->dim; axi++) {
706 cropdo[axi] = AIR_TRUE1;
707 }
708 for (axi=0; axi<keepNum; axi++) {
709 if (!( keep[axi] < nin->dim )) {
710 biffAddf(NRRDnrrdBiffKey, "%s: keep[%u] %u out of range [0,%u]", me,
711 axi, keep[axi], nin->dim-1);
712 return 1;
713 }
714 if (!cropdo[keep[axi]]) {
715 biffAddf(NRRDnrrdBiffKey, "%s: keep[%u] %u redundant", me, axi, keep[axi]);
716 return 1;
717 }
718 cropdo[keep[axi]] = AIR_FALSE0;
719 }
720 if (keepNum == nin->dim) {
721 /* weird- wanted to keep all axes and crop none; that's easy */
722 if (nrrdCopy(nout, nin)) {
723 biffAddf(NRRDnrrdBiffKey, "%s: trouble copying for trivial case", me);
724 return 1;
725 }
726 return 0;
727 }
728
729 /* else there's work to do */
730 mop = airMopNew();
731 nperm = nrrdNew();
732 airMopAdd(mop, nperm, (airMopper)nrrdNuke, airMopAlways);
733 nline = nrrdNew();
734 airMopAdd(mop, nline, (airMopper)nrrdNuke, airMopAlways);
735 for (axi=0; axi<nin->dim; axi++) {
736 double wsum, part;
737 min[axi] = 0;
738 max[axi] = nin->axis[axi].size-1;
739 if (!cropdo[axi]) {
740 continue;
741 }
742 /* else some analysis is required for this axis */
743 /* NN is product of axes NOT being cropped */
744 NN = nrrdElementNumber(nin)/nin->axis[axi].size;
745 if (nrrdAxesSwap(nperm, nin, axi, nin->dim-1)
746 || nrrdReshape_va(nperm, nperm, 2, NN, nin->axis[axi].size)
747 || nrrdProject(nline, nperm, 0, measr, nrrdTypeDouble)) {
748 biffAddf(NRRDnrrdBiffKey, "%s: trouble forming projection line", me);
749 airMopError(mop); return 1;
750 }
751 /* find sum of array */
752 line = AIR_CAST(double *, nline->data)((double *)(nline->data));
753 wsum = part = 0.0;
754 for (II=0; II<nin->axis[axi].size; II++) {
755 wsum += line[II];
756 }
757 /* sum bottom of array until hit fraction */
758 for (II=0; II<nin->axis[axi].size; II++) {
759 part += line[II];
760 if (part/wsum > frac) {
761 min[axi] = II;
762 break;
763 }
764 }
765 if (II == nin->axis[axi].size) {
766 biffAddf(NRRDnrrdBiffKey, "%s: confusion on bottom of axis %u", me, axi);
767 airMopError(mop); return 1;
768 }
769 /* sum top of array until hit fraction */
770 part = 0.0;
771 for (II=nin->axis[axi].size; II>0; II--) {
772 part += line[II-1];
773 if (part/wsum > frac) {
774 max[axi] = II-1;
775 break;
776 }
777 }
778 if (II == 0) {
779 biffAddf(NRRDnrrdBiffKey, "%s: confusion on top of axis %u", me, axi);
780 airMopError(mop); return 1;
781 }
782 /*
783 fprintf(stderr, "!%s: axis %u [%u,%u] --> ", me, axi,
784 AIR_CAST(unsigned int, min[axi]),
785 AIR_CAST(unsigned int, max[axi]));
786 */
787 /* adjust based on offset */
788 if (offset > 0) {
789 if (min[axi] < AIR_CAST(size_t, offset)((size_t)(offset))) {
790 /* desired outwards offset is more than cropping set */
791 min[axi] = 0;
792 } else {
793 min[axi] -= offset;
794 }
795 max[axi] += offset;
796 max[axi] = AIR_MIN(max[axi], nin->axis[axi].size-1)((max[axi]) < (nin->axis[axi].size-1) ? (max[axi]) : (nin
->axis[axi].size-1));
797 }
798 /*
799 fprintf(stderr, "[%u,%u]\n",
800 AIR_CAST(unsigned int, min[axi]),
801 AIR_CAST(unsigned int, max[axi]));
802 */
803 }
804
805 /* can now do the crop */
806 if (nrrdCrop(nout, nin, min, max)) {
807 biffAddf(NRRDnrrdBiffKey, "%s: trouble doing the crop", me);
808 return 1;
809 }
810 /* save the extents */
811 for (axi=0; axi<nin->dim; axi++) {
812 if (_min) {
813 _min[axi] = min[axi];
814 }
815 if (_max) {
816 _max[axi] = max[axi];
817 }
818 }
819 airMopOkay(mop);
820 return 0;
821}
822
823/* ---- END non-NrrdIO */