Bug Summary

File:src/nrrd/filt.c
Location:line 718, column 58
Description:Array access (from variable 'dd') results in a null pointer dereference

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
27int
28_nrrdCM_median(const float *hist, float half) {
29 float sum = 0;
30 const float *hpt;
31
32 hpt = hist;
33
34 while(sum < half)
35 sum += *hpt++;
36
37 return AIR_CAST(int, hpt - 1 - hist)((int)(hpt - 1 - hist));
38}
39
40int
41_nrrdCM_mode(const float *hist, int bins) {
42 float max;
43 int i, mi;
44
45 mi = -1;
46 max = 0;
47 for (i=0; i<bins; i++) {
48 if (hist[i] && (!max || hist[i] > max) ) {
49 max = hist[i];
50 mi = i;
51 }
52 }
53
54 return mi;
55}
56
57#define INDEX(nin, range, lup, idxIn, bins, val)( val = (lup)((nin)->data, (idxIn)), airIndex((range)->
min, (val), (range)->max, bins) ) ( \
58 val = (lup)((nin)->data, (idxIn)), \
59 airIndex((range)->min, (val), (range)->max, bins) \
60)
61
62void
63_nrrdCM_printhist(const float *hist, int bins, const char *desc) {
64 int i;
65
66 printf("%s:\n", desc);
67 for (i=0; i<bins; i++) {
68 if (hist[i]) {
69 printf(" %d: %g\n", i, hist[i]);
70 }
71 }
72}
73
74float *
75_nrrdCM_wtAlloc(int radius, float wght) {
76 float *wt, sum;
77 int diam, r;
78
79 diam = 2*radius + 1;
80 wt = (float *)calloc(diam, sizeof(float));
81 wt[radius] = 1.0;
82 for (r=1; r<=radius; r++) {
83 wt[radius+r] = AIR_CAST(float, pow(1.0/wght, r))((float)(pow(1.0/wght, r)));
84 wt[radius-r] = AIR_CAST(float, pow(1.0/wght, r))((float)(pow(1.0/wght, r)));
85 }
86 sum = 0.0;
87 for (r=0; r<diam; r++) {
88 sum += wt[r];
89 }
90 for (r=0; r<diam; r++) {
91 wt[r] /= sum;
92 }
93 /*
94 for (r=0; r<diam; r++) {
95 fprintf(stderr, "!%s: wt[%d] = %g\n", "_nrrdCM_wtAlloc", r, wt[r]);
96 }
97 */
98 return wt;
99}
100
101void
102_nrrdCheapMedian1D(Nrrd *nout, const Nrrd *nin, const NrrdRange *range,
103 int radius, float wght,
104 int bins, int mode, float *hist) {
105 /* static const char me[]="_nrrdCheapMedian1D"; */
106 size_t num;
107 int X, I, idx, diam;
108 float half, *wt;
109 double val, (*lup)(const void *, size_t);
110
111 diam = 2*radius + 1;
112 lup = nrrdDLookup[nin->type];
113 num = nrrdElementNumber(nin);
114 if (1 == wght) {
115 /* uniform weighting-> can do sliding histogram optimization */
116 /* initialize histogram */
117 half = AIR_CAST(float, diam/2 + 1)((float)(diam/2 + 1));
118 memset(hist, 0, bins*sizeof(float))__builtin___memset_chk (hist, 0, bins*sizeof(float), __builtin_object_size
(hist, 0));
119 for (X=0; X<diam; X++) {
120 hist[INDEX(nin, range, lup, X, bins, val)( val = (lup)((nin)->data, (X)), airIndex((range)->min,
(val), (range)->max, bins) )]++;
121 }
122 /* _nrrdCM_printhist(hist, bins, "after init"); */
123 /* find median at each point using existing histogram */
124 for (X=radius; X<(int)num-radius; X++) {
125 /* _nrrdCM_printhist(hist, bins, "----------"); */
126 idx = mode ? _nrrdCM_mode(hist, bins) : _nrrdCM_median(hist, half);
127 val = NRRD_NODE_POS(range->min, range->max, bins, idx)( ((double)((range->max))-((range->min)))*((double)((idx
))-(0)) / ((double)((bins)-1)-(0)) + ((range->min)));
128 /* printf(" median idx = %d -> val = %g\n", idx, val); */
129 nrrdDInsert[nout->type](nout->data, X, val);
130 /* probably update histogram for next iteration */
131 if (X < (int)num-radius-1) {
132 hist[INDEX(nin, range, lup, X+radius+1, bins, val)( val = (lup)((nin)->data, (X+radius+1)), airIndex((range)
->min, (val), (range)->max, bins) )]++;
133 hist[INDEX(nin, range, lup, X-radius, bins, val)( val = (lup)((nin)->data, (X-radius)), airIndex((range)->
min, (val), (range)->max, bins) )]--;
134 }
135 }
136 } else {
137 /* non-uniform weighting --> slow and stupid */
138 wt = _nrrdCM_wtAlloc(radius, wght);
139 half = 0.5;
140 for (X=radius; X<(int)num-radius; X++) {
141 memset(hist, 0, bins*sizeof(float))__builtin___memset_chk (hist, 0, bins*sizeof(float), __builtin_object_size
(hist, 0));
142 for (I=-radius; I<=radius; I++) {
143 hist[INDEX(nin, range, lup, I+X, bins, val)( val = (lup)((nin)->data, (I+X)), airIndex((range)->min
, (val), (range)->max, bins) )] += wt[I+radius];
144 }
145 idx = mode ? _nrrdCM_mode(hist, bins) : _nrrdCM_median(hist, half);
146 val = NRRD_NODE_POS(range->min, range->max, bins, idx)( ((double)((range->max))-((range->min)))*((double)((idx
))-(0)) / ((double)((bins)-1)-(0)) + ((range->min)));
147 nrrdDInsert[nout->type](nout->data, X, val);
148 }
149 free(wt);
150 }
151}
152
153void
154_nrrdCheapMedian2D(Nrrd *nout, const Nrrd *nin, const NrrdRange *range,
155 int radius, float wght,
156 int bins, int mode, float *hist) {
157 /* static const char me[]="_nrrdCheapMedian2D"; */
158 int X, Y, I, J;
159 int sx, sy, idx, diam;
160 float half, *wt;
161 double val, (*lup)(const void *, size_t);
162
163 diam = 2*radius + 1;
164 sx = AIR_CAST(unsigned int, nin->axis[0].size)((unsigned int)(nin->axis[0].size));
165 sy = AIR_CAST(unsigned int, nin->axis[1].size)((unsigned int)(nin->axis[1].size));
166 lup = nrrdDLookup[nin->type];
167 if (1 == wght) {
168 /* uniform weighting-> can do sliding histogram optimization */
169 half = AIR_CAST(float, diam*diam/2 + 1)((float)(diam*diam/2 + 1));
170 for (Y=radius; Y<sy-radius; Y++) {
171 /* initialize histogram */
172 memset(hist, 0, bins*sizeof(float))__builtin___memset_chk (hist, 0, bins*sizeof(float), __builtin_object_size
(hist, 0));
173 X = radius;
174 for (J=-radius; J<=radius; J++) {
175 for (I=-radius; I<=radius; I++) {
176 hist[INDEX(nin, range, lup, X+I + sx*(J+Y), bins, val)( val = (lup)((nin)->data, (X+I + sx*(J+Y))), airIndex((range
)->min, (val), (range)->max, bins) )]++;
177 }
178 }
179 /* _nrrdCM_printhist(hist, bins, "after init"); */
180 /* find median at each point using existing histogram */
181 for (X=radius; X<sx-radius; X++) {
182 idx = mode ? _nrrdCM_mode(hist, bins) : _nrrdCM_median(hist, half);
183 val = NRRD_NODE_POS(range->min, range->max, bins, idx)( ((double)((range->max))-((range->min)))*((double)((idx
))-(0)) / ((double)((bins)-1)-(0)) + ((range->min)));
184 nrrdDInsert[nout->type](nout->data, X + sx*Y, val);
185 /* probably update histogram for next iteration */
186 if (X < sx-radius-1) {
187 for (J=-radius; J<=radius; J++) {
188 hist[INDEX(nin, range, lup, X+radius+1 + sx*(J+Y), bins, val)( val = (lup)((nin)->data, (X+radius+1 + sx*(J+Y))), airIndex
((range)->min, (val), (range)->max, bins) )]++;
189 hist[INDEX(nin, range, lup, X-radius + sx*(J+Y), bins, val)( val = (lup)((nin)->data, (X-radius + sx*(J+Y))), airIndex
((range)->min, (val), (range)->max, bins) )]--;
190 }
191 }
192 }
193 }
194 } else {
195 /* non-uniform weighting --> slow and stupid */
196 wt = _nrrdCM_wtAlloc(radius, wght);
197 half = 0.5;
198 for (Y=radius; Y<sy-radius; Y++) {
199 for (X=radius; X<sx-radius; X++) {
200 memset(hist, 0, bins*sizeof(float))__builtin___memset_chk (hist, 0, bins*sizeof(float), __builtin_object_size
(hist, 0));
201 for (J=-radius; J<=radius; J++) {
202 for (I=-radius; I<=radius; I++) {
203 hist[INDEX(nin, range, lup, I+X + sx*(J+Y),( val = (lup)((nin)->data, (I+X + sx*(J+Y))), airIndex((range
)->min, (val), (range)->max, bins) )
204 bins, val)( val = (lup)((nin)->data, (I+X + sx*(J+Y))), airIndex((range
)->min, (val), (range)->max, bins) )] += wt[I+radius]*wt[J+radius];
205 }
206 }
207 idx = mode ? _nrrdCM_mode(hist, bins) : _nrrdCM_median(hist, half);
208 val = NRRD_NODE_POS(range->min, range->max, bins, idx)( ((double)((range->max))-((range->min)))*((double)((idx
))-(0)) / ((double)((bins)-1)-(0)) + ((range->min)));
209 nrrdDInsert[nout->type](nout->data, X + sx*Y, val);
210 }
211 }
212 free(wt);
213 }
214}
215
216void
217_nrrdCheapMedian3D(Nrrd *nout, const Nrrd *nin, const NrrdRange *range,
218 int radius, float wght,
219 int bins, int mode, float *hist) {
220 static const char me[]="_nrrdCheapMedian3D";
221 char done[13];
222 int X, Y, Z, I, J, K;
223 int sx, sy, sz, idx, diam;
224 float half, *wt;
225 double val, (*lup)(const void *, size_t);
226
227 diam = 2*radius + 1;
228 sx = AIR_CAST(int, nin->axis[0].size)((int)(nin->axis[0].size));
229 sy = AIR_CAST(int, nin->axis[1].size)((int)(nin->axis[1].size));
230 sz = AIR_CAST(int, nin->axis[2].size)((int)(nin->axis[2].size));
231 lup = nrrdDLookup[nin->type];
232 fprintf(stderr__stderrp, "%s: ... ", me);
233 if (1 == wght) {
234 /* uniform weighting-> can do sliding histogram optimization */
235 half = AIR_CAST(float, diam*diam*diam/2 + 1)((float)(diam*diam*diam/2 + 1));
236 fflush(stderr__stderrp);
237 for (Z=radius; Z<sz-radius; Z++) {
238 fprintf(stderr__stderrp, "%s", airDoneStr(radius, Z, sz-radius-1, done));
239 fflush(stderr__stderrp);
240 for (Y=radius; Y<sy-radius; Y++) {
241 /* initialize histogram */
242 memset(hist, 0, bins*sizeof(float))__builtin___memset_chk (hist, 0, bins*sizeof(float), __builtin_object_size
(hist, 0));
243 X = radius;
244 for (K=-radius; K<=radius; K++) {
245 for (J=-radius; J<=radius; J++) {
246 for (I=-radius; I<=radius; I++) {
247 hist[INDEX(nin, range, lup, I+X + sx*(J+Y + sy*(K+Z)),( val = (lup)((nin)->data, (I+X + sx*(J+Y + sy*(K+Z)))), airIndex
((range)->min, (val), (range)->max, bins) )
248 bins, val)( val = (lup)((nin)->data, (I+X + sx*(J+Y + sy*(K+Z)))), airIndex
((range)->min, (val), (range)->max, bins) )]++;
249 }
250 }
251 }
252 /* find median at each point using existing histogram */
253 for (X=radius; X<sx-radius; X++) {
254 idx = mode ? _nrrdCM_mode(hist, bins) : _nrrdCM_median(hist, half);
255 val = NRRD_NODE_POS(range->min, range->max, bins, idx)( ((double)((range->max))-((range->min)))*((double)((idx
))-(0)) / ((double)((bins)-1)-(0)) + ((range->min)));
256 nrrdDInsert[nout->type](nout->data, X + sx*(Y + sy*Z), val);
257 /* probably update histogram for next iteration */
258 if (X < sx-radius-1) {
259 for (K=-radius; K<=radius; K++) {
260 for (J=-radius; J<=radius; J++) {
261 hist[INDEX(nin, range, lup, X+radius+1 + sx*(J+Y + sy*(K+Z)),( val = (lup)((nin)->data, (X+radius+1 + sx*(J+Y + sy*(K+Z
)))), airIndex((range)->min, (val), (range)->max, bins)
)
262 bins, val)( val = (lup)((nin)->data, (X+radius+1 + sx*(J+Y + sy*(K+Z
)))), airIndex((range)->min, (val), (range)->max, bins)
)]++;
263 hist[INDEX(nin, range, lup, X-radius + sx*(J+Y + sy*(K+Z)),( val = (lup)((nin)->data, (X-radius + sx*(J+Y + sy*(K+Z))
)), airIndex((range)->min, (val), (range)->max, bins) )
264 bins, val)( val = (lup)((nin)->data, (X-radius + sx*(J+Y + sy*(K+Z))
)), airIndex((range)->min, (val), (range)->max, bins) )]--;
265 }
266 }
267 }
268 }
269 }
270 }
271 } else {
272 /* non-uniform weighting --> slow and stupid */
273 wt = _nrrdCM_wtAlloc(radius, wght);
274 half = 0.5;
275 for (Z=radius; Z<sz-radius; Z++) {
276 fprintf(stderr__stderrp, "%s", airDoneStr(radius, Z, sz-radius-1, done));
277 fflush(stderr__stderrp);
278 for (Y=radius; Y<sy-radius; Y++) {
279 for (X=radius; X<sx-radius; X++) {
280 memset(hist, 0, bins*sizeof(float))__builtin___memset_chk (hist, 0, bins*sizeof(float), __builtin_object_size
(hist, 0));
281 for (K=-radius; K<=radius; K++) {
282 for (J=-radius; J<=radius; J++) {
283 for (I=-radius; I<=radius; I++) {
284 hist[INDEX(nin, range, lup, I+X + sx*(J+Y + sy*(K+Z)),( val = (lup)((nin)->data, (I+X + sx*(J+Y + sy*(K+Z)))), airIndex
((range)->min, (val), (range)->max, bins) )
285 bins, val)( val = (lup)((nin)->data, (I+X + sx*(J+Y + sy*(K+Z)))), airIndex
((range)->min, (val), (range)->max, bins) )]
286 += wt[I+radius]*wt[J+radius]*wt[K+radius];
287 }
288 }
289 }
290 idx = mode ? _nrrdCM_mode(hist, bins) : _nrrdCM_median(hist, half);
291 val = NRRD_NODE_POS(range->min, range->max, bins, idx)( ((double)((range->max))-((range->min)))*((double)((idx
))-(0)) / ((double)((bins)-1)-(0)) + ((range->min)));
292 nrrdDInsert[nout->type](nout->data, X + sx*(Y + sy*Z), val);
293 }
294 }
295 }
296 free(wt);
297 }
298 fprintf(stderr__stderrp, "\b\b\b\b\b\b done\n");
299}
300
301/* HEY: total copy-and-paste from _nrrdCheapMedian3D */
302void
303_nrrdCheapMedian4D(Nrrd *nout, const Nrrd *nin, const NrrdRange *range,
304 int radius, float wght,
305 int bins, int mode, float *hist) {
306 static const char me[]="_nrrdCheapMedian4D";
307 char done[13];
308 int W, X, Y, Z, H, I, J, K;
309 int sw, sx, sy, sz, idx, diam;
310 float half, *wt;
311 double val, (*lup)(const void *, size_t);
312
313 diam = 2*radius + 1;
314 sw = AIR_CAST(int, nin->axis[0].size)((int)(nin->axis[0].size));
315 sx = AIR_CAST(int, nin->axis[1].size)((int)(nin->axis[1].size));
316 sy = AIR_CAST(int, nin->axis[2].size)((int)(nin->axis[2].size));
317 sz = AIR_CAST(int, nin->axis[3].size)((int)(nin->axis[3].size));
318 lup = nrrdDLookup[nin->type];
319 fprintf(stderr__stderrp, "%s: ... ", me);
320 if (1 == wght) {
321 /* uniform weighting-> can do sliding histogram optimization */
322 half = AIR_CAST(float, diam*diam*diam*diam/2 + 1)((float)(diam*diam*diam*diam/2 + 1));
323 fflush(stderr__stderrp);
324 for (Z=radius; Z<sz-radius; Z++) {
325 fprintf(stderr__stderrp, "%s", airDoneStr(radius, Z, sz-radius-1, done));
326 fflush(stderr__stderrp);
327 for (Y=radius; Y<sy-radius; Y++) {
328 for (X=radius; X<sx-radius; X++) {
329 /* initialize histogram */
330 memset(hist, 0, bins*sizeof(float))__builtin___memset_chk (hist, 0, bins*sizeof(float), __builtin_object_size
(hist, 0));
331 W = radius;
332 for (K=-radius; K<=radius; K++) {
333 for (J=-radius; J<=radius; J++) {
334 for (I=-radius; I<=radius; I++) {
335 for (H=-radius; H<=radius; H++) {
336 hist[INDEX(nin, range, lup,( val = (lup)((nin)->data, (H+W + sw*(I+X + sx*(J+Y + sy*(
K+Z))))), airIndex((range)->min, (val), (range)->max, bins
) )
337 H+W + sw*(I+X + sx*(J+Y + sy*(K+Z))),( val = (lup)((nin)->data, (H+W + sw*(I+X + sx*(J+Y + sy*(
K+Z))))), airIndex((range)->min, (val), (range)->max, bins
) )
338 bins, val)( val = (lup)((nin)->data, (H+W + sw*(I+X + sx*(J+Y + sy*(
K+Z))))), airIndex((range)->min, (val), (range)->max, bins
) )]++;
339 }
340 }
341 }
342 }
343 /* find median at each point using existing histogram */
344 for (W=radius; W<sw-radius; W++) {
345 idx = mode ? _nrrdCM_mode(hist, bins) : _nrrdCM_median(hist, half);
346 val = NRRD_NODE_POS(range->min, range->max, bins, idx)( ((double)((range->max))-((range->min)))*((double)((idx
))-(0)) / ((double)((bins)-1)-(0)) + ((range->min)));
347 nrrdDInsert[nout->type](nout->data, W + sw*(X + sx*(Y + sy*Z)), val);
348 /* probably update histogram for next iteration */
349 if (W < sw-radius-1) {
350 for (K=-radius; K<=radius; K++) {
351 for (J=-radius; J<=radius; J++) {
352 for (I=-radius; I<=radius; I++) {
353 hist[INDEX(nin, range, lup, W+radius+1 + sw*(I+X + sx*(J+Y + sy*(K+Z))),( val = (lup)((nin)->data, (W+radius+1 + sw*(I+X + sx*(J+Y
+ sy*(K+Z))))), airIndex((range)->min, (val), (range)->
max, bins) )
354 bins, val)( val = (lup)((nin)->data, (W+radius+1 + sw*(I+X + sx*(J+Y
+ sy*(K+Z))))), airIndex((range)->min, (val), (range)->
max, bins) )]++;
355 hist[INDEX(nin, range, lup, W-radius + sw*(I+X + sx*(J+Y + sy*(K+Z))),( val = (lup)((nin)->data, (W-radius + sw*(I+X + sx*(J+Y +
sy*(K+Z))))), airIndex((range)->min, (val), (range)->max
, bins) )
356 bins, val)( val = (lup)((nin)->data, (W-radius + sw*(I+X + sx*(J+Y +
sy*(K+Z))))), airIndex((range)->min, (val), (range)->max
, bins) )]--;
357 }
358 }
359 }
360 }
361 }
362 }
363 }
364 }
365 } else {
366 /* non-uniform weighting --> slow and stupid */
367 wt = _nrrdCM_wtAlloc(radius, wght);
368 half = 0.5;
369 for (Z=radius; Z<sz-radius; Z++) {
370 fprintf(stderr__stderrp, "%s", airDoneStr(radius, Z, sz-radius-1, done));
371 fflush(stderr__stderrp);
372 for (Y=radius; Y<sy-radius; Y++) {
373 for (X=radius; X<sx-radius; X++) {
374 for (W=radius; W<sw-radius; W++) {
375 memset(hist, 0, bins*sizeof(float))__builtin___memset_chk (hist, 0, bins*sizeof(float), __builtin_object_size
(hist, 0));
376 for (K=-radius; K<=radius; K++) {
377 for (J=-radius; J<=radius; J++) {
378 for (I=-radius; I<=radius; I++) {
379 for (H=-radius; H<=radius; H++) {
380 hist[INDEX(nin, range, lup,( val = (lup)((nin)->data, (H+W + sw*(I+X + sx*(J+Y + sy*(
K+Z))))), airIndex((range)->min, (val), (range)->max, bins
) )
381 H+W + sw*(I+X + sx*(J+Y + sy*(K+Z))),( val = (lup)((nin)->data, (H+W + sw*(I+X + sx*(J+Y + sy*(
K+Z))))), airIndex((range)->min, (val), (range)->max, bins
) )
382 bins, val)( val = (lup)((nin)->data, (H+W + sw*(I+X + sx*(J+Y + sy*(
K+Z))))), airIndex((range)->min, (val), (range)->max, bins
) )]
383 += wt[H+radius]*wt[I+radius]*wt[J+radius]*wt[K+radius];
384 }
385 }
386 }
387 }
388 idx = mode ? _nrrdCM_mode(hist, bins) : _nrrdCM_median(hist, half);
389 val = NRRD_NODE_POS(range->min, range->max, bins, idx)( ((double)((range->max))-((range->min)))*((double)((idx
))-(0)) / ((double)((bins)-1)-(0)) + ((range->min)));
390 nrrdDInsert[nout->type](nout->data, W + sw*(X + sx*(Y + sy*Z)), val);
391 }
392 }
393 }
394 }
395 free(wt);
396 }
397 fprintf(stderr__stderrp, "\b\b\b\b\b\b done\n");
398}
399
400/*
401******** nrrdCheapMedian
402**
403** histogram-based median or mode filtering
404** !mode: median filtering
405** mode: mode filtering
406*/
407int
408nrrdCheapMedian(Nrrd *_nout, const Nrrd *_nin,
409 int pad, int mode,
410 unsigned int radius, float wght, unsigned int bins) {
411 static const char me[]="nrrdCheapMedian", func[]="cmedian";
412 NrrdRange *range;
413 float *hist;
414 Nrrd *nout, *nin;
415 airArray *mop;
416 unsigned int minsize;
417
418 if (!(_nin && _nout)) {
419 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
420 return 1;
421 }
422 if (!(radius >= 1)) {
423 biffAddf(NRRDnrrdBiffKey, "%s: need radius >= 1 (got %d)", me, radius);
424 return 1;
425 }
426 if (!(bins >= 1)) {
427 biffAddf(NRRDnrrdBiffKey, "%s: need bins >= 1 (got %d)", me, bins);
428 return 1;
429 }
430 if (!(AIR_IN_CL(1, _nin->dim, 4)((1) <= (_nin->dim) && (_nin->dim) <= (4)
))) {
431 biffAddf(NRRDnrrdBiffKey, "%s: sorry, can only handle dim 1, 2, 3 (not %d)",
432 me, _nin->dim);
433 return 1;
434 }
435 minsize = AIR_CAST(unsigned int, _nin->axis[0].size)((unsigned int)(_nin->axis[0].size));
436 if (_nin->dim > 1) {
437 minsize = AIR_MIN(minsize, AIR_CAST(unsigned int, _nin->axis[1].size))((minsize) < (((unsigned int)(_nin->axis[1].size))) ? (
minsize) : (((unsigned int)(_nin->axis[1].size))));
438 }
439 if (_nin->dim > 2) {
440 minsize = AIR_MIN(minsize, AIR_CAST(unsigned int, _nin->axis[2].size))((minsize) < (((unsigned int)(_nin->axis[2].size))) ? (
minsize) : (((unsigned int)(_nin->axis[2].size))));
441 }
442 if (!pad && minsize < 2*radius+1) {
443 biffAddf(NRRDnrrdBiffKey, "%s: minimum nrrd size (%d) smaller than filtering "
444 "window size (%d) with radius %d; must enable padding", me,
445 minsize, 2*radius+1, radius);
446 return 1;
447 }
448 if (_nout == _nin) {
449 biffAddf(NRRDnrrdBiffKey, "%s: nout==nin disallowed", me);
450 return 1;
451 }
452 if (nrrdTypeBlock == _nin->type) {
453 biffAddf(NRRDnrrdBiffKey, "%s: can't filter nrrd type %s", me,
454 airEnumStr(nrrdType, nrrdTypeBlock));
455 return 1;
456 }
457
458 mop = airMopNew();
459 /* set nin based on _nin */
460 airMopAdd(mop, nin=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
461 if (pad) {
462 airMopAdd(mop, nout=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
463 if (nrrdSimplePad_va(nin, _nin, radius, nrrdBoundaryBleed)) {
464 biffAddf(NRRDnrrdBiffKey, "%s: trouble padding input", me);
465 airMopError(mop); return 1;
466 }
467 } else {
468 if (nrrdCopy(nin, _nin)) {
469 biffAddf(NRRDnrrdBiffKey, "%s: trouble copying input", me);
470 airMopError(mop); return 1;
471 }
472 nout = _nout;
473 }
474 if (nrrdCopy(nout, nin)) {
475 biffAddf(NRRDnrrdBiffKey, "%s: failed to create initial copy of input", me);
476 airMopError(mop); return 1;
477 }
478 range = nrrdRangeNewSet(nin, nrrdBlind8BitRangeFalse);
479 airMopAdd(mop, range, (airMopper)nrrdRangeNix, airMopAlways);
480 if (!(hist = (float*)calloc(bins, sizeof(float)))) {
481 biffAddf(NRRDnrrdBiffKey, "%s: couldn't allocate histogram (%d bins)", me, bins);
482 airMopError(mop); return 1;
483 }
484 airMopAdd(mop, hist, airFree, airMopAlways);
485 if (!AIR_EXISTS(wght)(((int)(!((wght) - (wght)))))) {
486 wght = 1.0;
487 }
488 switch (nin->dim) {
489 case 1:
490 _nrrdCheapMedian1D(nout, nin, range, radius, wght, bins, mode, hist);
491 break;
492 case 2:
493 _nrrdCheapMedian2D(nout, nin, range, radius, wght, bins, mode, hist);
494 break;
495 case 3:
496 _nrrdCheapMedian3D(nout, nin, range, radius, wght, bins, mode, hist);
497 break;
498 case 4:
499 _nrrdCheapMedian4D(nout, nin, range, radius, wght, bins, mode, hist);
500 break;
501 default:
502 biffAddf(NRRDnrrdBiffKey, "%s: sorry, %d-dimensional median unimplemented",
503 me, nin->dim);
504 airMopError(mop); return 1;
505 }
506
507 nrrdAxisInfoCopy(nout, nin, NULL((void*)0), NRRD_AXIS_INFO_NONE0);
508 if (nrrdContentSet_va(nout, func, nin, "%d,%d,%g,%d",
509 mode, radius, wght, bins)) {
510 biffAddf(NRRDnrrdBiffKey, "%s:", me);
511 airMopError(mop); return 1;
512 }
513 /* basic info handled by nrrdCopy above */
514
515 /* set _nout based on nout */
516 if (pad) {
517 if (nrrdSimpleCrop(_nout, nout, radius)) {
518 biffAddf(NRRDnrrdBiffKey, "%s: trouble cropping output", me);
519 airMopError(mop); return 1;
520 }
521 } else {
522 /* we've already set output in _nout == nout */
523 }
524
525 airMopOkay(mop);
526 return 0;
527}
528
529/*
530** returns intersection of parabolas c(x) = spc^2 (x - xi)^2 + yi
531*/
532static double
533intx(double xx0, double yy0, double xx1, double yy1, double spc) {
534 double ss;
535
536 ss = spc*spc;
537 return (yy1/ss + xx1*xx1 - (yy0/ss + xx0*xx0))/(2*(xx1 - xx0));
538}
539
540/*
541** squared-distance transform of single scanline
542**
543** based on code published with:
544** Distance Transforms of Sampled Functions
545** Pedro F. Felzenszwalb and Daniel P. Huttenlocher
546** Cornell Computing and Information Science TR2004-1963
547**
548** dd: output (pre-allocated for "len")
549** ff: input function (pre-allocated for "len")
550** zz: buffer (pre-allocated for "len"+1) for locations of
551** boundaries between parabolas
552** vv: buffer (pre-allocated for "len") for locations of
553** parabolas in lower envelope
554**
555** The major modification from the published method is the inclusion
556** of the "spc" parameter that gives the inter-sample spacing, so that
557** the multi-dimensional version can work on non-isotropic samples.
558**
559** FLT_MIN and FLT_MAX are used to be consistent with the
560** initialization in nrrdDistanceL2(), which uses FLT_MAX
561** to be compatible with the case of using floats
562*/
563static void
564distanceL2Sqrd1D(double *dd, const double *ff,
565 double *zz, unsigned int *vv,
566 size_t len, double spc) {
567 size_t kk, qq;
568
569 if (!( dd && ff && zz && vv && len > 0 )) {
570 /* error */
571 return;
572 }
573
574 kk = 0;
575 vv[0] = 0;
576 zz[0] = -FLT_MAX3.40282347e+38F;
577 zz[1] = +FLT_MAX3.40282347e+38F;
578 for (qq=1; qq<len; qq++) {
579 double ss;
580 ss = intx(AIR_CAST(double, qq)((double)(qq)), ff[qq], vv[kk], ff[vv[kk]], spc);
581 /* while (ss <= zz[kk]) {
582 ** HEY this can have kk going to -1 and into memory errors!
583 */
584 while (ss <= zz[kk] && kk) {
585 kk--;
586 ss = intx(AIR_CAST(double, qq)((double)(qq)), ff[qq], vv[kk], ff[vv[kk]], spc);
587 }
588 kk++;
589 vv[kk] = AIR_CAST(unsigned int, qq)((unsigned int)(qq));
590 zz[kk] = ss;
591 zz[kk+1] = +FLT_MAX3.40282347e+38F;
592 }
593
594 kk = 0;
595 for (qq=00; qq<len; qq++) {
596 double dx;
597 while (zz[kk+1] < qq) {
598 kk++;
599 }
600 /* cast to avoid overflow weirdness on the unsigned ints */
601 dx = AIR_CAST(double, qq)((double)(qq)) - vv[kk];
602 dd[qq] = spc*spc*dx*dx + ff[vv[kk]];
603 }
604
605 return;
606}
607
608static int
609distanceL2Sqrd(Nrrd *ndist, double *spcMean) {
610 static const char me[]="distanceL2Sqrd";
611 size_t sizeMax; /* max size of all axes */
612 Nrrd *ntmpA, *ntmpB, *npass[NRRD_DIM_MAX16+1];
613 int spcSomeExist, spcSomeNonExist;
614 unsigned int di, *vv;
615 double *dd, *ff, *zz;
616 double spc[NRRD_DIM_MAX16], vector[NRRD_SPACE_DIM_MAX8];
617 double (*lup)(const void *, size_t), (*ins)(void *, size_t, double);
618 airArray *mop;
619
620 if (!( nrrdTypeFloat == ndist->type || nrrdTypeDouble == ndist->type )) {
1
Taking true branch
→
621 /* MWC: This error condition was/is being ignored. */
622 /* biffAddf(NRRD, "%s: sorry, can only process type %s or %s (not %s)",
623 me,
624 airEnumStr(nrrdType, nrrdTypeFloat),
625 airEnumStr(nrrdType, nrrdTypeDouble),
626 airEnumStr(nrrdType, ndist->type));
627 */
628 }
629
630 spcSomeExist = AIR_FALSE0;
631 spcSomeNonExist = AIR_FALSE0;
632 for (di=0; di<ndist->dim; di++) {
2
←
Loop condition is true. Entering loop body
→
3
←
Loop condition is false. Execution continues on line 637
→
633 nrrdSpacingCalculate(ndist, di, spc + di, vector);
634 spcSomeExist |= AIR_EXISTS(spc[di])(((int)(!((spc[di]) - (spc[di])))));
635 spcSomeNonExist |= !AIR_EXISTS(spc[di])(((int)(!((spc[di]) - (spc[di])))));
636 }
637 if (spcSomeExist && spcSomeNonExist) {
4
←
Taking false branch
→
638 biffAddf(NRRDnrrdBiffKey, "%s: axis spacings must all exist or all non-exist", me);
639 return 1;
640 }
641 if (!spcSomeExist) {
5
←
Taking false branch
→
642 for (di=0; di<ndist->dim; di++) {
643 spc[di] = 1.0;
644 }
645 }
646 *spcMean = 0;
647 for (di=0; di<ndist->dim; di++) {
6
←
Loop condition is true. Entering loop body
→
7
←
Loop condition is false. Execution continues on line 650
→
648 *spcMean += spc[di];
649 }
650 *spcMean /= ndist->dim;
651
652 sizeMax = 0;
653 for (di=0; di<ndist->dim; di++) {
8
←
Loop condition is true. Entering loop body
→
9
←
Loop condition is false. Execution continues on line 658
→
654 sizeMax = AIR_MAX(sizeMax, ndist->axis[di].size)((sizeMax) > (ndist->axis[di].size) ? (sizeMax) : (ndist
->axis[di].size));
655 }
656
657 /* create mop and allocate tmp buffers */
658 mop = airMopNew();
659 ntmpA = nrrdNew();
660 airMopAdd(mop, ntmpA, (airMopper)nrrdNuke, airMopAlways);
661 if (ndist->dim > 2) {
10
←
Taking false branch
→
662 ntmpB = nrrdNew();
663 airMopAdd(mop, ntmpB, (airMopper)nrrdNuke, airMopAlways);
664 } else {
665 ntmpB = NULL((void*)0);
666 }
667 if (nrrdCopy(ntmpA, ndist)
668 || (ndist->dim > 2 && nrrdCopy(ntmpB, ndist))) {
669 biffAddf(NRRDnrrdBiffKey, "%s: couldn't allocate image buffers", me);
670 }
671 dd = AIR_CAST(double *, calloc(sizeMax, sizeof(double)))((double *)(calloc(sizeMax, sizeof(double))));
11
←
Value assigned to 'dd'
→
672 ff = AIR_CAST(double *, calloc(sizeMax, sizeof(double)))((double *)(calloc(sizeMax, sizeof(double))));
673 zz = AIR_CAST(double *, calloc(sizeMax+1, sizeof(double)))((double *)(calloc(sizeMax+1, sizeof(double))));
674 vv = AIR_CAST(unsigned int *, calloc(sizeMax, sizeof(unsigned int)))((unsigned int *)(calloc(sizeMax, sizeof(unsigned int))));
675 airMopAdd(mop, dd, airFree, airMopAlways);
676 airMopAdd(mop, ff, airFree, airMopAlways);
677 airMopAdd(mop, zz, airFree, airMopAlways);
678 airMopAdd(mop, vv, airFree, airMopAlways);
679 if (!( dd && ff && zz && vv )) {
12
←
Assuming pointer value is null
→
13
←
Taking true branch
→
680 /* MWC: This error condition was/is being ignored. */
681 /* biffAddf(NRRD, "%s: couldn't allocate scanline buffers", me); */
682 }
683
684 /* set up array of buffers */
685 npass[0] = ndist;
686 for (di=1; di<ndist->dim; di++) {
14
←
Loop condition is false. Execution continues on line 689
→
687 npass[di] = (di % 2) ? ntmpA : ntmpB;
688 }
689 npass[ndist->dim] = ndist;
690
691 /* run the multiple passes */
692 /* what makes the indexing here so simple is that by assuming that
693 we're processing every axis, the input to any given pass can be
694 logically considered a 2-D image (a list of scanlines), where the
695 second axis is the merge of all input axes but the first. With
696 the rotational shuffle of axes through passes, the initial axis
697 and the set of other axes swap places, so its like the 2-D image
698 is being transposed. NOTE: the Nrrds that were allocated as
699 buffers are really being mis-used, in that the axis sizes and
700 raster ordering of what we're storing there is *not* the same as
701 told by axis[].size */
702 lup = nrrdDLookup[ndist->type];
703 ins = nrrdDInsert[ndist->type];
704 for (di=0; di<ndist->dim; di++) {
15
←
Loop condition is true. Entering loop body
→
705 size_t lineIdx, lineNum, valIdx, valNum;
706
707 valNum = ndist->axis[di].size;
708 lineNum = nrrdElementNumber(ndist)/valNum;
709 for (lineIdx=0; lineIdx<lineNum; lineIdx++) {
16
←
Assuming 'lineIdx' is < 'lineNum'
→
17
←
Loop condition is true. Entering loop body
→
710 /* read input scanline into ff */
711 for (valIdx=0; valIdx<valNum; valIdx++) {
18
←
Loop condition is true. Entering loop body
→
19
←
Assuming 'valIdx' is >= 'valNum'
→
20
←
Loop condition is false. Execution continues on line 715
→
712 ff[valIdx] = lup(npass[di]->data, valIdx + valNum*lineIdx);
713 }
714 /* do the transform */
715 distanceL2Sqrd1D(dd, ff, zz, vv, valNum, spc[di]);
716 /* write dd to output scanline */
717 for (valIdx=0; valIdx<valNum; valIdx++) {
21
←
Loop condition is true. Entering loop body
→
718 ins(npass[di+1]->data, lineIdx + lineNum*valIdx, dd[valIdx]);
22
←
Array access (from variable 'dd') results in a null pointer dereference
719 }
720 }
721 }
722
723 airMopOkay(mop);
724 return 0;
725}
726
727/*
728** helper function for distance transforms, is called by things that want to do
729** specific kinds of transforms.
730*/
731static int
732_distanceBase(Nrrd *nout, const Nrrd *nin,
733 int typeOut, const int *axisDo,
734 double thresh, double bias, int insideHigher) {
735 static const char me[]="_distanceBase";
736 size_t ii, nn;
737 double (*lup)(const void *, size_t), (*ins)(void *, size_t, double);
738 double spcMean;
739
740 if (!( nout && nin )) {
741 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
742 return 1;
743 }
744 if (nrrdTypeBlock == nin->type) {
745 biffAddf(NRRDnrrdBiffKey, "%s: need scalar type for distance transform (not %s)",
746 me, airEnumStr(nrrdType, nrrdTypeBlock));
747 return 1;
748 }
749 if (!( nrrdTypeDouble == typeOut || nrrdTypeFloat == typeOut )) {
750 biffAddf(NRRDnrrdBiffKey, "%s: sorry, can only transform to type %s or %s "
751 "(not %s)", me,
752 airEnumStr(nrrdType, nrrdTypeFloat),
753 airEnumStr(nrrdType, nrrdTypeDouble),
754 airEnumStr(nrrdType, typeOut));
755 return 1;
756 }
757 if (axisDo) {
758 biffAddf(NRRDnrrdBiffKey, "%s: sorry, selective axis transform not implemented",
759 me);
760 return 1;
761 }
762 if (!AIR_EXISTS(thresh)(((int)(!((thresh) - (thresh)))))) {
763 biffAddf(NRRDnrrdBiffKey, "%s: threshold (%g) doesn't exist", me, thresh);
764 return 1;
765 }
766
767 if (nrrdConvert(nout, nin, typeOut)) {
768 biffAddf(NRRDnrrdBiffKey, "%s: couldn't allocate output", me);
769 return 1;
770 }
771 lup = nrrdDLookup[nout->type];
772 ins = nrrdDInsert[nout->type];
773
774 nn = nrrdElementNumber(nout);
775 for (ii=0; ii<nn; ii++) {
776 double val, bb;
777 val = lup(nout->data, ii);
778 if (insideHigher) {
779 bb = bias*(val - thresh);
780 ins(nout->data, ii, val > thresh ? bb*bb : FLT_MAX3.40282347e+38F);
781 } else {
782 bb = bias*(thresh - val);
783 ins(nout->data, ii, val <= thresh ? bb*bb : FLT_MAX3.40282347e+38F);
784 }
785 }
786
787 if (distanceL2Sqrd(nout, &spcMean)) {
788 biffAddf(NRRDnrrdBiffKey, "%s: trouble doing transform", me);
789 return 1;
790 }
791
792 for (ii=0; ii<nn; ii++) {
793 double val;
794 val = sqrt(lup(nout->data, ii));
795 /* here's where the distance is tweaked downwards by half a sample width */
796 ins(nout->data, ii, AIR_MAX(0, val - spcMean/2)((0) > (val - spcMean/2) ? (0) : (val - spcMean/2)));
797 }
798
799 return 0;
800}
801
802/*
803******** nrrdDistanceL2
804**
805** computes euclidean (L2) distance transform of input image, after
806** thresholding at "thresh".
807**
808** NOTE: the output of this is slightly offset from what one might
809** expect; decreased by half of the average (over all axes) sample
810** spacing. The reason for this is so that when the transform is
811** applied to the inverted image and negated, to create a full
812** signed distance map, the transition from interior to exterior
813** distance values is smooth. Without this trick, there is a
814** small little plateau at the transition.
815*/
816int
817nrrdDistanceL2(Nrrd *nout, const Nrrd *nin,
818 int typeOut, const int *axisDo,
819 double thresh, int insideHigher) {
820 static const char me[]="nrrdDistanceL2";
821
822 if (_distanceBase(nout, nin, typeOut, axisDo, thresh, 0, insideHigher)) {
823 biffAddf(NRRDnrrdBiffKey, "%s: trouble doing distance transform", me);
824 return 1;
825 }
826 return 0;
827}
828
829int
830nrrdDistanceL2Biased(Nrrd *nout, const Nrrd *nin,
831 int typeOut, const int *axisDo,
832 double thresh, double bias, int insideHigher) {
833 static const char me[]="nrrdDistanceL2Biased";
834
835 if (_distanceBase(nout, nin, typeOut, axisDo, thresh, bias, insideHigher)) {
836 biffAddf(NRRDnrrdBiffKey, "%s: trouble doing distance transform", me);
837 return 1;
838 }
839 return 0;
840}
841
842int
843nrrdDistanceL2Signed(Nrrd *nout, const Nrrd *nin,
844 int typeOut, const int *axisDo,
845 double thresh, int insideHigher) {
846 static const char me[]="nrrdDistanceL2Signed";
847 airArray *mop;
848 Nrrd *ninv;
849
850 if (!(nout && nin)) {
851 biffAddf(NRRDnrrdBiffKey, "%s: got NULL pointer", me);
852 return 1;
853 }
854
855 mop = airMopNew();
856 ninv = nrrdNew();
857 airMopAdd(mop, ninv, (airMopper)nrrdNuke, airMopAlways);
858
859 if (nrrdDistanceL2(nout, nin, typeOut, axisDo, thresh, insideHigher)
860 || nrrdDistanceL2(ninv, nin, typeOut, axisDo, thresh, !insideHigher)
861 || nrrdArithUnaryOp(ninv, nrrdUnaryOpNegative, ninv)
862 || nrrdArithBinaryOp(nout, nrrdBinaryOpAdd, nout, ninv)) {
863 biffAddf(NRRDnrrdBiffKey, "%s: trouble doing or combining transforms", me);
864 airMopError(mop); return 1;
865 }
866
867 airMopOkay(mop);
868 return 0;
869}