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 |
|
|
#define CLAMP_HIST_BINS_MIN 512 |
28 |
|
|
#define CLAMP_PERC_MAX 30.0 |
29 |
|
|
|
30 |
|
|
#define DEBUG 0 |
31 |
|
|
|
32 |
|
|
/* TODO: |
33 |
|
|
* |
34 |
|
|
* make sure all outout values exist (test with d/xiao/todering/xing) |
35 |
|
|
* |
36 |
|
|
* permit controlling the extent of correction (don't subtract out |
37 |
|
|
* all the estimated ring signal) |
38 |
|
|
* |
39 |
|
|
* valgrind |
40 |
|
|
* |
41 |
|
|
* make it multi-threaded |
42 |
|
|
* |
43 |
|
|
* try fix for round object boundaries being confused for rings |
44 |
|
|
- (relies on properties of discrete gauss for radial blurring) |
45 |
|
|
- after initial ptxf |
46 |
|
|
- make a few radial blurs (up to scale of radial high-pass) |
47 |
|
|
- measuring scale-normalized |df/dr| on each one |
48 |
|
|
- do non-maximal suppression along radius, zeroing out edges below some |
49 |
|
|
percentile of gradient strength |
50 |
|
|
- where the gradients are high on the most-blurring, and had similar |
51 |
|
|
grad mag at smaller scales, that's a real edge: make a mask for this |
52 |
|
|
- blur this along theta just like the ring map, then multiply w/ ring map |
53 |
|
|
* |
54 |
|
|
* try fix for low-theta-frequency rings |
55 |
|
|
* with high thetaNum, find mode along theta |
56 |
|
|
*/ |
57 |
|
|
|
58 |
|
|
NrrdDeringContext * |
59 |
|
|
nrrdDeringContextNew(void) { |
60 |
|
|
NrrdDeringContext *drc; |
61 |
|
|
unsigned int pi; |
62 |
|
|
|
63 |
|
|
drc = AIR_CALLOC(1, NrrdDeringContext); |
64 |
|
|
if (!drc) { |
65 |
|
|
return NULL; |
66 |
|
|
} |
67 |
|
|
drc->verbose = 0; |
68 |
|
|
drc->linearInterp = AIR_FALSE; |
69 |
|
|
drc->verticalSeam = AIR_FALSE; |
70 |
|
|
drc->nin = NULL; |
71 |
|
|
drc->center[0] = AIR_NAN; |
72 |
|
|
drc->center[1] = AIR_NAN; |
73 |
|
|
drc->clampPerc[0] = 0.0; |
74 |
|
|
drc->clampPerc[1] = 0.0; |
75 |
|
|
drc->radiusScale = 1.0; |
76 |
|
|
drc->thetaNum = 0; |
77 |
|
|
drc->clampHistoBins = CLAMP_HIST_BINS_MIN*4; |
78 |
|
|
drc->rkernel = NULL; |
79 |
|
|
drc->tkernel = NULL; |
80 |
|
|
for (pi=0; pi<NRRD_KERNEL_PARMS_NUM; pi++) { |
81 |
|
|
drc->rkparm[pi] = drc->tkparm[pi] = AIR_NAN; |
82 |
|
|
} |
83 |
|
|
drc->cdataIn = NULL; |
84 |
|
|
drc->cdataOut = NULL; |
85 |
|
|
drc->sliceSize = 0; |
86 |
|
|
drc->clampDo = AIR_FALSE; |
87 |
|
|
drc->clamp[0] = AIR_NAN; |
88 |
|
|
drc->clamp[1] = AIR_NAN; |
89 |
|
|
drc->ringMagnitude = AIR_NAN; |
90 |
|
|
return drc; |
91 |
|
|
} |
92 |
|
|
|
93 |
|
|
NrrdDeringContext * |
94 |
|
|
nrrdDeringContextNix(NrrdDeringContext *drc) { |
95 |
|
|
|
96 |
|
|
if (drc) { |
97 |
|
|
free(drc); |
98 |
|
|
} |
99 |
|
|
return NULL; |
100 |
|
|
} |
101 |
|
|
|
102 |
|
|
int |
103 |
|
|
nrrdDeringVerboseSet(NrrdDeringContext *drc, int verbose) { |
104 |
|
|
static const char me[]="nrrdDeringVerboseSet"; |
105 |
|
|
|
106 |
|
|
if (!drc) { |
107 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
108 |
|
|
return 1; |
109 |
|
|
} |
110 |
|
|
|
111 |
|
|
drc->verbose = verbose; |
112 |
|
|
return 0; |
113 |
|
|
} |
114 |
|
|
|
115 |
|
|
int |
116 |
|
|
nrrdDeringLinearInterpSet(NrrdDeringContext *drc, int linterp) { |
117 |
|
|
static const char me[]="nrrdDeringLinearInterpSet"; |
118 |
|
|
|
119 |
|
|
if (!drc) { |
120 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
121 |
|
|
return 1; |
122 |
|
|
} |
123 |
|
|
|
124 |
|
|
drc->linearInterp = linterp; |
125 |
|
|
return 0; |
126 |
|
|
} |
127 |
|
|
|
128 |
|
|
int |
129 |
|
|
nrrdDeringVerticalSeamSet(NrrdDeringContext *drc, int vertSeam) { |
130 |
|
|
static const char me[]="nrrdDeringVerticalSeamSet"; |
131 |
|
|
|
132 |
|
|
if (!drc) { |
133 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
134 |
|
|
return 1; |
135 |
|
|
} |
136 |
|
|
|
137 |
|
|
drc->verticalSeam = vertSeam; |
138 |
|
|
return 0; |
139 |
|
|
} |
140 |
|
|
|
141 |
|
|
int |
142 |
|
|
nrrdDeringInputSet(NrrdDeringContext *drc, const Nrrd *nin) { |
143 |
|
|
static const char me[]="nrrdDeringInputSet"; |
144 |
|
|
|
145 |
|
|
if (!( drc && nin )) { |
146 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
147 |
|
|
return 1; |
148 |
|
|
} |
149 |
|
|
if (nrrdCheck(nin)) { |
150 |
|
|
biffAddf(NRRD, "%s: problems with given nrrd", me); |
151 |
|
|
return 1; |
152 |
|
|
} |
153 |
|
|
if (nrrdTypeBlock == nin->type) { |
154 |
|
|
biffAddf(NRRD, "%s: can't resample from type %s", me, |
155 |
|
|
airEnumStr(nrrdType, nrrdTypeBlock)); |
156 |
|
|
return 1; |
157 |
|
|
} |
158 |
|
|
if (!( 2 == nin->dim || 3 == nin->dim )) { |
159 |
|
|
biffAddf(NRRD, "%s: need 2 or 3 dim nrrd (not %u)", me, nin->dim); |
160 |
|
|
return 1; |
161 |
|
|
} |
162 |
|
|
|
163 |
|
|
if (drc->verbose > 2) { |
164 |
|
|
fprintf(stderr, "%s: hi\n", me); |
165 |
|
|
} |
166 |
|
|
drc->nin = nin; |
167 |
|
|
drc->cdataIn = AIR_CAST(const char *, nin->data); |
168 |
|
|
drc->sliceSize = (nin->axis[0].size |
169 |
|
|
* nin->axis[1].size |
170 |
|
|
* nrrdElementSize(nin)); |
171 |
|
|
if (drc->verbose > 2) { |
172 |
|
|
fprintf(stderr, "%s: sliceSize = %u\n", me, |
173 |
|
|
AIR_CAST(unsigned int, drc->sliceSize)); |
174 |
|
|
} |
175 |
|
|
|
176 |
|
|
return 0; |
177 |
|
|
} |
178 |
|
|
|
179 |
|
|
int |
180 |
|
|
nrrdDeringCenterSet(NrrdDeringContext *drc, double cx, double cy) { |
181 |
|
|
static const char me[]="nrrdDeringCenterSet"; |
182 |
|
|
|
183 |
|
|
if (!drc) { |
184 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
185 |
|
|
return 1; |
186 |
|
|
} |
187 |
|
|
if (!( AIR_EXISTS(cx) && AIR_EXISTS(cy) )) { |
188 |
|
|
biffAddf(NRRD, "%s: center (%g,%g) doesn't exist", me, cx, cy); |
189 |
|
|
return 1; |
190 |
|
|
} |
191 |
|
|
|
192 |
|
|
drc->center[0] = cx; |
193 |
|
|
drc->center[1] = cy; |
194 |
|
|
|
195 |
|
|
return 0; |
196 |
|
|
} |
197 |
|
|
|
198 |
|
|
int |
199 |
|
|
nrrdDeringClampPercSet(NrrdDeringContext *drc, |
200 |
|
|
double lo, double hi) { |
201 |
|
|
static const char me[]="nrrdDeringClampPercSet"; |
202 |
|
|
|
203 |
|
|
if (!drc) { |
204 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
205 |
|
|
return 1; |
206 |
|
|
} |
207 |
|
|
if (!( AIR_EXISTS(lo) && AIR_EXISTS(hi) |
208 |
|
|
&& lo >= 0 && lo < CLAMP_PERC_MAX |
209 |
|
|
&& hi >= 0 && hi < CLAMP_PERC_MAX)) { |
210 |
|
|
biffAddf(NRRD, "%s: need finite lo and hi both in [0.0, %g), not %g, %g", |
211 |
|
|
me, CLAMP_PERC_MAX, lo, hi); |
212 |
|
|
return 1; |
213 |
|
|
} |
214 |
|
|
|
215 |
|
|
drc->clampPerc[0] = lo; |
216 |
|
|
drc->clampPerc[1] = hi; |
217 |
|
|
|
218 |
|
|
return 0; |
219 |
|
|
} |
220 |
|
|
|
221 |
|
|
int |
222 |
|
|
nrrdDeringClampHistoBinsSet(NrrdDeringContext *drc, |
223 |
|
|
unsigned int bins) { |
224 |
|
|
static const char me[]="nrrdDeringClampHistoBinsSet"; |
225 |
|
|
|
226 |
|
|
if (!drc) { |
227 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
228 |
|
|
return 1; |
229 |
|
|
} |
230 |
|
|
if (!( bins >= CLAMP_HIST_BINS_MIN )) { |
231 |
|
|
biffAddf(NRRD, "%s: given bins %u not >= reasonable min %u", |
232 |
|
|
me, bins, CLAMP_HIST_BINS_MIN); |
233 |
|
|
return 1; |
234 |
|
|
} |
235 |
|
|
|
236 |
|
|
drc->clampHistoBins = bins; |
237 |
|
|
|
238 |
|
|
return 0; |
239 |
|
|
} |
240 |
|
|
|
241 |
|
|
int |
242 |
|
|
nrrdDeringRadiusScaleSet(NrrdDeringContext *drc, double rsc) { |
243 |
|
|
static const char me[]="nrrdDeringRadiusScaleSet"; |
244 |
|
|
|
245 |
|
|
if (!drc) { |
246 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
247 |
|
|
return 1; |
248 |
|
|
} |
249 |
|
|
if (!( AIR_EXISTS(rsc) && rsc > 0.0 )) { |
250 |
|
|
biffAddf(NRRD, "%s: need finite positive radius scale, not %g", me, rsc); |
251 |
|
|
return 1; |
252 |
|
|
} |
253 |
|
|
|
254 |
|
|
drc->radiusScale = rsc; |
255 |
|
|
|
256 |
|
|
return 0; |
257 |
|
|
} |
258 |
|
|
|
259 |
|
|
int |
260 |
|
|
nrrdDeringThetaNumSet(NrrdDeringContext *drc, unsigned int thetaNum) { |
261 |
|
|
static const char me[]="nrrdDeringThetaNumSet"; |
262 |
|
|
|
263 |
|
|
if (!drc) { |
264 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
265 |
|
|
return 1; |
266 |
|
|
} |
267 |
|
|
if (!thetaNum) { |
268 |
|
|
biffAddf(NRRD, "%s: need non-zero thetaNum", me); |
269 |
|
|
return 1; |
270 |
|
|
} |
271 |
|
|
|
272 |
|
|
drc->thetaNum = thetaNum; |
273 |
|
|
|
274 |
|
|
return 0; |
275 |
|
|
} |
276 |
|
|
|
277 |
|
|
int |
278 |
|
|
nrrdDeringRadialKernelSet(NrrdDeringContext *drc, |
279 |
|
|
const NrrdKernel *rkernel, |
280 |
|
|
const double rkparm[NRRD_KERNEL_PARMS_NUM]) { |
281 |
|
|
static const char me[]="nrrdDeringRadialKernelSet"; |
282 |
|
|
unsigned int pi; |
283 |
|
|
|
284 |
|
|
if (!( drc && rkernel && rkparm )) { |
285 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
286 |
|
|
return 1; |
287 |
|
|
} |
288 |
|
|
|
289 |
|
|
drc->rkernel = rkernel; |
290 |
|
|
for (pi=0; pi<NRRD_KERNEL_PARMS_NUM; pi++) { |
291 |
|
|
drc->rkparm[pi] = rkparm[pi]; |
292 |
|
|
} |
293 |
|
|
|
294 |
|
|
return 0; |
295 |
|
|
} |
296 |
|
|
|
297 |
|
|
int |
298 |
|
|
nrrdDeringThetaKernelSet(NrrdDeringContext *drc, |
299 |
|
|
const NrrdKernel *tkernel, |
300 |
|
|
const double tkparm[NRRD_KERNEL_PARMS_NUM]) { |
301 |
|
|
static const char me[]="nrrdDeringThetaKernelSet"; |
302 |
|
|
unsigned int pi; |
303 |
|
|
|
304 |
|
|
if (!( drc && tkernel && tkparm )) { |
305 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
306 |
|
|
return 1; |
307 |
|
|
} |
308 |
|
|
|
309 |
|
|
drc->tkernel = tkernel; |
310 |
|
|
for (pi=0; pi<NRRD_KERNEL_PARMS_NUM; pi++) { |
311 |
|
|
drc->tkparm[pi] = tkparm[pi]; |
312 |
|
|
} |
313 |
|
|
|
314 |
|
|
return 0; |
315 |
|
|
} |
316 |
|
|
|
317 |
|
|
/* |
318 |
|
|
** per-thread state for deringing |
319 |
|
|
*/ |
320 |
|
|
#define ORIG 0 |
321 |
|
|
#define WGHT 1 |
322 |
|
|
#define BLRR 2 |
323 |
|
|
#define DIFF 3 |
324 |
|
|
#define RSHP 4 |
325 |
|
|
#define CROP 5 |
326 |
|
|
#define CBLR 6 |
327 |
|
|
#define RING 7 |
328 |
|
|
#define PTXF_NUM 8 |
329 |
|
|
typedef struct { |
330 |
|
|
unsigned int zi; |
331 |
|
|
double radMax; |
332 |
|
|
size_t radNum; |
333 |
|
|
airArray *mop; |
334 |
|
|
Nrrd *nsliceOrig, /* wrapped slice of nin, sneakily non-const */ |
335 |
|
|
*nslice, /* slice of nin, converted to double */ |
336 |
|
|
*nsliceOut, /* (may be different type than nslice) */ |
337 |
|
|
*nptxf[PTXF_NUM]; |
338 |
|
|
double *slice, *ptxf, *wght, *ring; |
339 |
|
|
NrrdResampleContext *rsmc[2]; /* rsmc[0]: radial blurring ORIG -> BLRR |
340 |
|
|
rsmc[1]: theta blurring DIFF -> RING */ |
341 |
|
|
double ringMag; |
342 |
|
|
} deringBag; |
343 |
|
|
|
344 |
|
|
static deringBag * |
345 |
|
|
deringBagNew(NrrdDeringContext *drc, double radMax) { |
346 |
|
|
deringBag *dbg; |
347 |
|
|
unsigned int pi; |
348 |
|
|
|
349 |
|
|
dbg = AIR_CALLOC(1, deringBag); |
350 |
|
|
dbg->radMax = radMax; |
351 |
|
|
dbg->radNum = AIR_ROUNDUP(drc->radiusScale*radMax); |
352 |
|
|
|
353 |
|
|
dbg->mop = airMopNew(); |
354 |
|
|
dbg->nsliceOrig = nrrdNew(); |
355 |
|
|
airMopAdd(dbg->mop, dbg->nsliceOrig, (airMopper)nrrdNix /* not Nuke! */, |
356 |
|
|
airMopAlways); |
357 |
|
|
dbg->nslice = nrrdNew(); |
358 |
|
|
airMopAdd(dbg->mop, dbg->nslice, (airMopper)nrrdNuke, airMopAlways); |
359 |
|
|
dbg->nsliceOut = nrrdNew(); |
360 |
|
|
airMopAdd(dbg->mop, dbg->nsliceOut, (airMopper)nrrdNix /* not Nuke! */, |
361 |
|
|
airMopAlways); |
362 |
|
|
for (pi=0; pi<PTXF_NUM; pi++) { |
363 |
|
|
dbg->nptxf[pi] = nrrdNew(); |
364 |
|
|
airMopAdd(dbg->mop, dbg->nptxf[pi], (airMopper)nrrdNuke, airMopAlways); |
365 |
|
|
} |
366 |
|
|
|
367 |
|
|
dbg->rsmc[0] = nrrdResampleContextNew(); |
368 |
|
|
airMopAdd(dbg->mop, dbg->rsmc[0], (airMopper)nrrdResampleContextNix, |
369 |
|
|
airMopAlways); |
370 |
|
|
dbg->rsmc[1] = nrrdResampleContextNew(); |
371 |
|
|
airMopAdd(dbg->mop, dbg->rsmc[1], (airMopper)nrrdResampleContextNix, |
372 |
|
|
airMopAlways); |
373 |
|
|
dbg->ringMag = 0.0; |
374 |
|
|
|
375 |
|
|
return dbg; |
376 |
|
|
} |
377 |
|
|
|
378 |
|
|
static deringBag * |
379 |
|
|
deringBagNix(deringBag *dbg) { |
380 |
|
|
|
381 |
|
|
airMopOkay(dbg->mop); |
382 |
|
|
airFree(dbg); |
383 |
|
|
return NULL; |
384 |
|
|
} |
385 |
|
|
|
386 |
|
|
static int |
387 |
|
|
deringPtxfAlloc(NrrdDeringContext *drc, deringBag *dbg) { |
388 |
|
|
static const char me[]="deringPtxfAlloc"; |
389 |
|
|
unsigned int pi, ri; |
390 |
|
|
int E; |
391 |
|
|
|
392 |
|
|
/* polar transform setup */ |
393 |
|
|
for (pi=0; pi<PTXF_NUM; pi++) { |
394 |
|
|
if (drc->verticalSeam && (CROP == pi || CBLR == pi)) { |
395 |
|
|
E = nrrdMaybeAlloc_va(dbg->nptxf[pi], nrrdTypeDouble, 3, |
396 |
|
|
dbg->radNum, |
397 |
|
|
AIR_CAST(size_t, drc->thetaNum/2 - 1), |
398 |
|
|
AIR_CAST(size_t, 2)); |
399 |
|
|
} else { |
400 |
|
|
E = nrrdMaybeAlloc_va(dbg->nptxf[pi], nrrdTypeDouble, 2, |
401 |
|
|
dbg->radNum, |
402 |
|
|
AIR_CAST(size_t, drc->thetaNum)); |
403 |
|
|
} |
404 |
|
|
if (E) { |
405 |
|
|
biffAddf(NRRD, "%s: polar transform allocation problem", me); |
406 |
|
|
return 1; |
407 |
|
|
} |
408 |
|
|
} |
409 |
|
|
dbg->ptxf = AIR_CAST(double *, dbg->nptxf[ORIG]->data); |
410 |
|
|
dbg->wght = AIR_CAST(double *, dbg->nptxf[WGHT]->data); |
411 |
|
|
dbg->ring = AIR_CAST(double *, dbg->nptxf[RING]->data); |
412 |
|
|
|
413 |
|
|
E = AIR_FALSE; |
414 |
|
|
for (ri=0; ri<2; ri++) { |
415 |
|
|
char kstr[AIR_STRLEN_LARGE]; |
416 |
|
|
if (0 == ri) { |
417 |
|
|
if (!E) E |= nrrdResampleInputSet(dbg->rsmc[0], dbg->nptxf[ORIG]); |
418 |
|
|
nrrdKernelSprint(kstr, drc->rkernel, drc->rkparm); |
419 |
|
|
if (drc->verbose > 1) { |
420 |
|
|
fprintf(stderr, "%s: rsmc[0] axis 0 kernel: %s\n", me, kstr); |
421 |
|
|
} |
422 |
|
|
if (!E) E |= nrrdResampleKernelSet(dbg->rsmc[0], 0, |
423 |
|
|
drc->rkernel, drc->rkparm); |
424 |
|
|
if (!E) E |= nrrdResampleKernelSet(dbg->rsmc[0], 1, NULL, NULL); |
425 |
|
|
if (!E) E |= nrrdResampleSamplesSet(dbg->rsmc[0], 1, drc->thetaNum); |
426 |
|
|
if (!E) E |= nrrdResampleBoundarySet(dbg->rsmc[0], nrrdBoundaryWrap); |
427 |
|
|
} else { |
428 |
|
|
if (drc->verticalSeam) { |
429 |
|
|
if (!E) E |= nrrdResampleInputSet(dbg->rsmc[1], dbg->nptxf[CROP]); |
430 |
|
|
} else { |
431 |
|
|
if (!E) E |= nrrdResampleInputSet(dbg->rsmc[1], dbg->nptxf[DIFF]); |
432 |
|
|
} |
433 |
|
|
nrrdKernelSprint(kstr, drc->tkernel, drc->tkparm); |
434 |
|
|
if (drc->verbose > 1) { |
435 |
|
|
fprintf(stderr, "%s: rsmc[1] axis 1 kernel: %s\n", me, kstr); |
436 |
|
|
} |
437 |
|
|
if (!E) E |= nrrdResampleKernelSet(dbg->rsmc[1], 0, NULL, NULL); |
438 |
|
|
if (!E) E |= nrrdResampleKernelSet(dbg->rsmc[1], 1, |
439 |
|
|
drc->tkernel, drc->tkparm); |
440 |
|
|
if (drc->verticalSeam) { |
441 |
|
|
if (!E) E |= nrrdResampleSamplesSet(dbg->rsmc[1], 1, drc->thetaNum/2 - 1); |
442 |
|
|
if (!E) E |= nrrdResampleKernelSet(dbg->rsmc[1], 2, NULL, NULL); |
443 |
|
|
if (!E) E |= nrrdResampleSamplesSet(dbg->rsmc[1], 2, 2); |
444 |
|
|
if (!E) E |= nrrdResampleBoundarySet(dbg->rsmc[1], nrrdBoundaryPad); |
445 |
|
|
if (!E) E |= nrrdResamplePadValueSet(dbg->rsmc[1], 0.0); |
446 |
|
|
} else { |
447 |
|
|
if (!E) E |= nrrdResampleSamplesSet(dbg->rsmc[1], 1, drc->thetaNum); |
448 |
|
|
if (!E) E |= nrrdResampleBoundarySet(dbg->rsmc[1], nrrdBoundaryWrap); |
449 |
|
|
} |
450 |
|
|
} |
451 |
|
|
if (!E) E |= nrrdResampleDefaultCenterSet(dbg->rsmc[ri], nrrdCenterCell); |
452 |
|
|
if (!E) E |= nrrdResampleSamplesSet(dbg->rsmc[ri], 0, dbg->radNum); |
453 |
|
|
if (!E) E |= nrrdResampleTypeOutSet(dbg->rsmc[ri], nrrdTypeDefault); |
454 |
|
|
if (!E) E |= nrrdResampleRenormalizeSet(dbg->rsmc[ri], AIR_TRUE); |
455 |
|
|
if (!E) E |= nrrdResampleNonExistentSet(dbg->rsmc[ri], |
456 |
|
|
nrrdResampleNonExistentRenormalize); |
457 |
|
|
if (!E) E |= nrrdResampleRangeFullSet(dbg->rsmc[ri], 0); |
458 |
|
|
if (!E) E |= nrrdResampleRangeFullSet(dbg->rsmc[ri], 1); |
459 |
|
|
} |
460 |
|
|
if (E) { |
461 |
|
|
biffAddf(NRRD, "%s: couldn't set up resampler", me); |
462 |
|
|
return 1; |
463 |
|
|
} |
464 |
|
|
|
465 |
|
|
|
466 |
|
|
return 0; |
467 |
|
|
} |
468 |
|
|
|
469 |
|
|
static int |
470 |
|
|
deringSliceGet(NrrdDeringContext *drc, deringBag *dbg, unsigned int zi) { |
471 |
|
|
static const char me[]="deringSliceGet"; |
472 |
|
|
void *nonconstdata; |
473 |
|
|
const char *cdata; |
474 |
|
|
|
475 |
|
|
/* HEY: bypass of const-ness of drc->cdataIn; should think about how |
476 |
|
|
to do this without breaking const-correctness... */ |
477 |
|
|
cdata = drc->cdataIn + zi*(drc->sliceSize); |
478 |
|
|
memcpy(&nonconstdata, &cdata, sizeof(void*)); |
479 |
|
|
/* slice setup */ |
480 |
|
|
if (nrrdWrap_va(dbg->nsliceOrig, |
481 |
|
|
nonconstdata, |
482 |
|
|
drc->nin->type, 2, |
483 |
|
|
drc->nin->axis[0].size, |
484 |
|
|
drc->nin->axis[1].size) |
485 |
|
|
|| (nrrdTypeDouble == drc->nin->type |
486 |
|
|
? nrrdCopy(dbg->nslice, dbg->nsliceOrig) |
487 |
|
|
: nrrdConvert(dbg->nslice, dbg->nsliceOrig, nrrdTypeDouble))) { |
488 |
|
|
biffAddf(NRRD, "%s: slice setup trouble", me); |
489 |
|
|
return 1; |
490 |
|
|
} |
491 |
|
|
dbg->slice = AIR_CAST(double *, dbg->nslice->data); |
492 |
|
|
dbg->zi = zi; |
493 |
|
|
|
494 |
|
|
return 0; |
495 |
|
|
} |
496 |
|
|
|
497 |
|
|
static int |
498 |
|
|
deringSliceSet(NrrdDeringContext *drc, deringBag *dbg, |
499 |
|
|
Nrrd *nout, unsigned int zi) { |
500 |
|
|
static const char me[]="deringSliceSet"; |
501 |
|
|
|
502 |
|
|
if (nrrdWrap_va(dbg->nsliceOut, |
503 |
|
|
AIR_CAST(void *, drc->cdataOut + zi*(drc->sliceSize)), |
504 |
|
|
nout->type, 2, |
505 |
|
|
nout->axis[0].size, |
506 |
|
|
nout->axis[1].size) |
507 |
|
|
|| (nrrdTypeDouble == nout->type |
508 |
|
|
? nrrdCopy(dbg->nsliceOut, dbg->nslice) |
509 |
|
|
: nrrdConvert(dbg->nsliceOut, dbg->nslice, nout->type))) { |
510 |
|
|
biffAddf(NRRD, "%s: slice output trouble", me); |
511 |
|
|
return 1; |
512 |
|
|
} |
513 |
|
|
|
514 |
|
|
return 0; |
515 |
|
|
} |
516 |
|
|
|
517 |
|
|
#define EPS 0.000001 |
518 |
|
|
|
519 |
|
|
static void |
520 |
|
|
deringXYtoRT(NrrdDeringContext *drc, deringBag *dbg, |
521 |
|
|
unsigned int xi, unsigned int yi, |
522 |
|
|
unsigned int *rrIdx, unsigned int *thIdx, |
523 |
|
|
double *rrFrc, double *thFrc) { |
524 |
|
|
double dx, dy, rr, th, rrScl, thScl; |
525 |
|
|
dx = xi - drc->center[0]; |
526 |
|
|
dy = yi - drc->center[1]; |
527 |
|
|
rr = sqrt(dx*dx + dy*dy); |
528 |
|
|
th = atan2(-dx, dy); |
529 |
|
|
rrScl = AIR_AFFINE(-EPS, rr, dbg->radMax+EPS, 0.0, dbg->radNum-1); |
530 |
|
|
*rrIdx = AIR_CAST(unsigned int, 0.5 + rrScl); |
531 |
|
|
thScl = AIR_AFFINE(-AIR_PI-EPS, th, AIR_PI+EPS, 0.0, drc->thetaNum); |
532 |
|
|
*thIdx = AIR_CAST(unsigned int, 0.5 + thScl); |
533 |
|
|
if (rrFrc && thFrc) { |
534 |
|
|
*rrFrc = rrScl - *rrIdx; |
535 |
|
|
*thFrc = thScl - *thIdx; |
536 |
|
|
if (*rrFrc < 0) { |
537 |
|
|
*rrIdx -= 1; |
538 |
|
|
*rrFrc += 1; |
539 |
|
|
} |
540 |
|
|
if (*thFrc < 0) { |
541 |
|
|
*thIdx -= 1; |
542 |
|
|
*thFrc += 1; |
543 |
|
|
} |
544 |
|
|
} |
545 |
|
|
return; |
546 |
|
|
} |
547 |
|
|
|
548 |
|
|
static int |
549 |
|
|
deringPtxfDo(NrrdDeringContext *drc, deringBag *dbg) { |
550 |
|
|
/* static const char me[]="deringPtxfDo"; */ |
551 |
|
|
unsigned int sx, sy, xi, yi, rrIdx, thIdx; |
552 |
|
|
|
553 |
|
|
nrrdSetZero(dbg->nptxf[ORIG]); |
554 |
|
|
nrrdSetZero(dbg->nptxf[WGHT]); |
555 |
|
|
sx = AIR_CAST(unsigned int, drc->nin->axis[0].size); |
556 |
|
|
sy = AIR_CAST(unsigned int, drc->nin->axis[1].size); |
557 |
|
|
for (yi=0; yi<sy; yi++) { |
558 |
|
|
for (xi=0; xi<sx; xi++) { |
559 |
|
|
double rrFrc, thFrc, val; |
560 |
|
|
if (drc->linearInterp) { |
561 |
|
|
unsigned int bidx, bidxPlus; |
562 |
|
|
deringXYtoRT(drc, dbg, xi, yi, &rrIdx, &thIdx, &rrFrc, &thFrc); |
563 |
|
|
bidx = AIR_UINT(rrIdx + dbg->radNum*thIdx); |
564 |
|
|
bidxPlus = AIR_UINT(thIdx < drc->thetaNum-1 |
565 |
|
|
? bidx + dbg->radNum |
566 |
|
|
: rrIdx); |
567 |
|
|
val = dbg->slice[xi + sx*yi]; |
568 |
|
|
if (drc->clampDo) { |
569 |
|
|
val = AIR_CLAMP(drc->clamp[0], val, drc->clamp[1]); |
570 |
|
|
} |
571 |
|
|
dbg->ptxf[bidx ] += (1-rrFrc)*(1-thFrc)*val; |
572 |
|
|
dbg->ptxf[bidx + 1] += rrFrc*(1-thFrc)*val; |
573 |
|
|
dbg->ptxf[bidxPlus ] += (1-rrFrc)*thFrc*val; |
574 |
|
|
dbg->ptxf[bidxPlus + 1] += rrFrc*thFrc*val; |
575 |
|
|
dbg->wght[bidx ] += (1-rrFrc)*(1-thFrc); |
576 |
|
|
dbg->wght[bidx + 1] += rrFrc*(1-thFrc); |
577 |
|
|
dbg->wght[bidxPlus ] += (1-rrFrc)*thFrc; |
578 |
|
|
dbg->wght[bidxPlus + 1] += rrFrc*thFrc; |
579 |
|
|
} else { |
580 |
|
|
deringXYtoRT(drc, dbg, xi, yi, &rrIdx, &thIdx, NULL, NULL); |
581 |
|
|
thIdx = thIdx % drc->thetaNum; |
582 |
|
|
dbg->ptxf[rrIdx + dbg->radNum*thIdx] += dbg->slice[xi + sx*yi]; |
583 |
|
|
dbg->wght[rrIdx + dbg->radNum*thIdx] += 1; |
584 |
|
|
} |
585 |
|
|
} |
586 |
|
|
} |
587 |
|
|
for (thIdx=0; thIdx<drc->thetaNum; thIdx++) { |
588 |
|
|
for (rrIdx=0; rrIdx<dbg->radNum; rrIdx++) { |
589 |
|
|
double tmpW; |
590 |
|
|
tmpW = dbg->wght[rrIdx + dbg->radNum*thIdx]; |
591 |
|
|
if (tmpW) { |
592 |
|
|
dbg->ptxf[rrIdx + dbg->radNum*thIdx] /= tmpW; |
593 |
|
|
} else { |
594 |
|
|
dbg->ptxf[rrIdx + dbg->radNum*thIdx] = AIR_NAN; |
595 |
|
|
} |
596 |
|
|
} |
597 |
|
|
} |
598 |
|
|
if (DEBUG) { |
599 |
|
|
char fname[AIR_STRLEN_SMALL]; |
600 |
|
|
sprintf(fname, "wght-%02u.nrrd", dbg->zi); |
601 |
|
|
nrrdSave(fname, dbg->nptxf[WGHT], NULL); |
602 |
|
|
} |
603 |
|
|
|
604 |
|
|
return 0; |
605 |
|
|
} |
606 |
|
|
|
607 |
|
|
static int |
608 |
|
|
deringPtxfFilter(NrrdDeringContext *drc, deringBag *dbg, unsigned int zi) { |
609 |
|
|
static const char me[]="deringPtxfFilter"; |
610 |
|
|
|
611 |
|
|
if ((!zi ? 0 : nrrdResampleInputSet(dbg->rsmc[0], dbg->nptxf[ORIG])) |
612 |
|
|
|| nrrdResampleExecute(dbg->rsmc[0], dbg->nptxf[BLRR]) |
613 |
|
|
|| nrrdArithBinaryOp(dbg->nptxf[DIFF], nrrdBinaryOpSubtract, |
614 |
|
|
dbg->nptxf[ORIG], dbg->nptxf[BLRR])) { |
615 |
|
|
biffAddf(NRRD, "%s: trouble with radial blur", me); |
616 |
|
|
return 1; |
617 |
|
|
} |
618 |
|
|
if (!drc->verticalSeam) { |
619 |
|
|
if ((!zi ? 0 : nrrdResampleInputSet(dbg->rsmc[1], dbg->nptxf[DIFF])) |
620 |
|
|
|| nrrdResampleExecute(dbg->rsmc[1], dbg->nptxf[RING])) { |
621 |
|
|
biffAddf(NRRD, "%s: trouble", me); |
622 |
|
|
return 1; |
623 |
|
|
} |
624 |
|
|
} else { |
625 |
|
|
size_t cmin[3], cmax[3]; |
626 |
|
|
ptrdiff_t pmin[3], pmax[3]; |
627 |
|
|
cmin[0] = 0; cmin[1] = 1; cmin[2] = 0; |
628 |
|
|
pmin[0] = 0; pmin[1] = -1; pmin[2] = 0; |
629 |
|
|
pmax[0] = cmax[0] = dbg->radNum - 1; |
630 |
|
|
cmax[1] = drc->thetaNum/2 - 1; |
631 |
|
|
pmax[1] = drc->thetaNum/2 - 2; /* max sample # of cropped axis 1 */ |
632 |
|
|
pmax[2] = cmax[2] = 1; |
633 |
|
|
if (nrrdAxesSplit(dbg->nptxf[RSHP], dbg->nptxf[DIFF], 1, |
634 |
|
|
drc->thetaNum/2, 2) |
635 |
|
|
|| nrrdCrop(dbg->nptxf[CROP], dbg->nptxf[RSHP], cmin, cmax) |
636 |
|
|
|| (!zi ? 0 : nrrdResampleInputSet(dbg->rsmc[1], dbg->nptxf[CROP])) |
637 |
|
|
|| nrrdResampleExecute(dbg->rsmc[1], dbg->nptxf[CBLR]) |
638 |
|
|
|| nrrdPad_nva(dbg->nptxf[RSHP], dbg->nptxf[CBLR], pmin, pmax, |
639 |
|
|
nrrdBoundaryPad, 0) |
640 |
|
|
|| nrrdAxesMerge(dbg->nptxf[RING], dbg->nptxf[RSHP], 1)) { |
641 |
|
|
biffAddf(NRRD, "%s: trouble with vertical seam", me); |
642 |
|
|
return 1; |
643 |
|
|
} |
644 |
|
|
if (DEBUG) { |
645 |
|
|
char fn[AIR_STRLEN_SMALL]; |
646 |
|
|
sprintf(fn, "rshp-%02u.nrrd", dbg->zi); nrrdSave(fn, dbg->nptxf[RSHP],NULL); |
647 |
|
|
sprintf(fn, "crop-%02u.nrrd", dbg->zi); nrrdSave(fn, dbg->nptxf[CROP],NULL); |
648 |
|
|
} |
649 |
|
|
} |
650 |
|
|
if (DEBUG) { |
651 |
|
|
char fn[AIR_STRLEN_SMALL]; |
652 |
|
|
sprintf(fn, "orig-%02u.nrrd", dbg->zi); nrrdSave(fn, dbg->nptxf[ORIG],NULL); |
653 |
|
|
sprintf(fn, "blrr-%02u.nrrd", dbg->zi); nrrdSave(fn, dbg->nptxf[BLRR],NULL); |
654 |
|
|
sprintf(fn, "diff-%02u.nrrd", dbg->zi); nrrdSave(fn, dbg->nptxf[DIFF],NULL); |
655 |
|
|
sprintf(fn, "ring-%02u.nrrd", dbg->zi); nrrdSave(fn, dbg->nptxf[RING],NULL); |
656 |
|
|
} |
657 |
|
|
|
658 |
|
|
return 0; |
659 |
|
|
} |
660 |
|
|
|
661 |
|
|
static int |
662 |
|
|
deringRingMagMeasure(NrrdDeringContext *drc, deringBag *dbg) { |
663 |
|
|
static const char me[]="deringRingMagMeasure"; |
664 |
|
|
airArray *mop; |
665 |
|
|
Nrrd *ntmp[2]; |
666 |
|
|
|
667 |
|
|
AIR_UNUSED(drc); |
668 |
|
|
mop = airMopNew(); |
669 |
|
|
ntmp[0] = nrrdNew(); |
670 |
|
|
airMopAdd(mop, ntmp[0], (airMopper)nrrdNuke, airMopAlways); |
671 |
|
|
ntmp[1] = nrrdNew(); |
672 |
|
|
airMopAdd(mop, ntmp[1], (airMopper)nrrdNuke, airMopAlways); |
673 |
|
|
if (nrrdReshape_va(ntmp[0], dbg->nptxf[RING], 2, |
674 |
|
|
(dbg->nptxf[RING]->axis[0].size |
675 |
|
|
* dbg->nptxf[RING]->axis[1].size), |
676 |
|
|
AIR_CAST(size_t, 1)) |
677 |
|
|
|| nrrdProject(ntmp[1], ntmp[0], 0, nrrdMeasureL2, nrrdTypeDouble)) { |
678 |
|
|
biffAddf(NRRD, "%s: trouble", me); |
679 |
|
|
airMopError(mop); return 1; |
680 |
|
|
} |
681 |
|
|
dbg->ringMag = *(AIR_CAST(double *, ntmp[1]->data)); |
682 |
|
|
|
683 |
|
|
airMopOkay(mop); |
684 |
|
|
return 0; |
685 |
|
|
} |
686 |
|
|
|
687 |
|
|
static int |
688 |
|
|
deringSubtract(NrrdDeringContext *drc, deringBag *dbg) { |
689 |
|
|
/* static const char me[]="deringSubtract"; */ |
690 |
|
|
unsigned int sx, sy, xi, yi, rrIdx, thIdx; |
691 |
|
|
|
692 |
|
|
sx = AIR_CAST(unsigned int, drc->nin->axis[0].size); |
693 |
|
|
sy = AIR_CAST(unsigned int, drc->nin->axis[1].size); |
694 |
|
|
for (yi=0; yi<sy; yi++) { |
695 |
|
|
for (xi=0; xi<sx; xi++) { |
696 |
|
|
double rrFrc, thFrc, val; |
697 |
|
|
if (drc->linearInterp) { |
698 |
|
|
unsigned int bidx, bidxPlus; |
699 |
|
|
deringXYtoRT(drc, dbg, xi, yi, &rrIdx, &thIdx, &rrFrc, &thFrc); |
700 |
|
|
bidx = AIR_UINT(rrIdx + dbg->radNum*thIdx); |
701 |
|
|
bidxPlus = AIR_UINT(thIdx < drc->thetaNum-1 |
702 |
|
|
? bidx + dbg->radNum |
703 |
|
|
: rrIdx); |
704 |
|
|
val = (dbg->ring[bidx ]*(1-rrFrc)*(1-thFrc) + |
705 |
|
|
dbg->ring[bidx + 1]*rrFrc*(1-thFrc) + |
706 |
|
|
dbg->ring[bidxPlus ]*(1-rrFrc)*thFrc + |
707 |
|
|
dbg->ring[bidxPlus + 1]*rrFrc*thFrc); |
708 |
|
|
dbg->slice[xi + sx*yi] -= val; |
709 |
|
|
} else { |
710 |
|
|
deringXYtoRT(drc, dbg, xi, yi, &rrIdx, &thIdx, NULL, NULL); |
711 |
|
|
thIdx = thIdx % drc->thetaNum; |
712 |
|
|
dbg->slice[xi + sx*yi] -= dbg->ring[rrIdx + dbg->radNum*thIdx]; |
713 |
|
|
} |
714 |
|
|
} |
715 |
|
|
} |
716 |
|
|
if (DEBUG) { |
717 |
|
|
char fname[AIR_STRLEN_SMALL]; |
718 |
|
|
sprintf(fname, "ring2-%02u.nrrd", dbg->zi); nrrdSave(fname, dbg->nptxf[RING], NULL); |
719 |
|
|
sprintf(fname, "drng-%02u.nrrd", dbg->zi); nrrdSave(fname, dbg->nslice, NULL); |
720 |
|
|
} |
721 |
|
|
return 0; |
722 |
|
|
} |
723 |
|
|
|
724 |
|
|
static int |
725 |
|
|
deringDo(NrrdDeringContext *drc, deringBag *dbg, |
726 |
|
|
Nrrd *nout, unsigned int zi) { |
727 |
|
|
static const char me[]="deringDo"; |
728 |
|
|
|
729 |
|
|
if (deringSliceGet(drc, dbg, zi) |
730 |
|
|
|| deringPtxfDo(drc, dbg) |
731 |
|
|
|| deringPtxfFilter(drc, dbg, zi) |
732 |
|
|
|| deringRingMagMeasure(drc, dbg) |
733 |
|
|
|| deringSubtract(drc, dbg) |
734 |
|
|
|| deringSliceSet(drc, dbg, nout, zi)) { |
735 |
|
|
biffAddf(NRRD, "%s: trouble", me); |
736 |
|
|
return 1; |
737 |
|
|
} |
738 |
|
|
|
739 |
|
|
return 0; |
740 |
|
|
} |
741 |
|
|
|
742 |
|
|
static int |
743 |
|
|
deringCheck(NrrdDeringContext *drc) { |
744 |
|
|
static const char me[]="deringCheck"; |
745 |
|
|
|
746 |
|
|
if (!(drc->nin)) { |
747 |
|
|
biffAddf(NRRD, "%s: no input set", me); |
748 |
|
|
return 1; |
749 |
|
|
} |
750 |
|
|
if (!( AIR_EXISTS(drc->center[0]) && AIR_EXISTS(drc->center[1]) )) { |
751 |
|
|
biffAddf(NRRD, "%s: no center set", me); |
752 |
|
|
return 1; |
753 |
|
|
} |
754 |
|
|
if (!(drc->thetaNum)) { |
755 |
|
|
biffAddf(NRRD, "%s: no thetaNum set", me); |
756 |
|
|
return 1; |
757 |
|
|
} |
758 |
|
|
if (!( drc->rkernel && drc->tkernel )) { |
759 |
|
|
biffAddf(NRRD, "%s: R and T kernels not both set", me); |
760 |
|
|
return 1; |
761 |
|
|
} |
762 |
|
|
if (drc->verticalSeam && !(0 == (drc->thetaNum % 2))) { |
763 |
|
|
biffAddf(NRRD, "%s: need even thetaNum (not %u) if wanting verticalSeam", |
764 |
|
|
me, drc->thetaNum); |
765 |
|
|
return 1; |
766 |
|
|
} |
767 |
|
|
return 0; |
768 |
|
|
} |
769 |
|
|
|
770 |
|
|
int |
771 |
|
|
nrrdDeringExecute(NrrdDeringContext *drc, Nrrd *nout) { |
772 |
|
|
static const char me[]="nrrdDeringExecute"; |
773 |
|
|
unsigned int sx, sy, sz, zi; |
774 |
|
|
double dx, dy, radLen, len; |
775 |
|
|
deringBag *dbg; |
776 |
|
|
airArray *mop; |
777 |
|
|
|
778 |
|
|
if (!( drc && nout )) { |
779 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
780 |
|
|
return 1; |
781 |
|
|
} |
782 |
|
|
if (deringCheck(drc)) { |
783 |
|
|
biffAddf(NRRD, "%s: trouble with setup", me); |
784 |
|
|
return 1; |
785 |
|
|
} |
786 |
|
|
|
787 |
|
|
/* initialize output */ |
788 |
|
|
if (nrrdCopy(nout, drc->nin)) { |
789 |
|
|
biffAddf(NRRD, "%s: trouble initializing output with input", me); |
790 |
|
|
return 1; |
791 |
|
|
} |
792 |
|
|
drc->cdataOut = AIR_CAST(char *, nout->data); |
793 |
|
|
|
794 |
|
|
mop = airMopNew(); |
795 |
|
|
|
796 |
|
|
/* set radLen: radial length of polar transform of data */ |
797 |
|
|
radLen = 0; |
798 |
|
|
sx = AIR_CAST(unsigned int, drc->nin->axis[0].size); |
799 |
|
|
sy = AIR_CAST(unsigned int, drc->nin->axis[1].size); |
800 |
|
|
dx = 0 - drc->center[0]; |
801 |
|
|
dy = 0 - drc->center[1]; |
802 |
|
|
len = sqrt(dx*dx + dy*dy); |
803 |
|
|
radLen = AIR_MAX(radLen, len); |
804 |
|
|
dx = sx-1 - drc->center[0]; |
805 |
|
|
dy = 0 - drc->center[1]; |
806 |
|
|
len = sqrt(dx*dx + dy*dy); |
807 |
|
|
radLen = AIR_MAX(radLen, len); |
808 |
|
|
dx = sx-1 - drc->center[0]; |
809 |
|
|
dy = sy-1 - drc->center[1]; |
810 |
|
|
len = sqrt(dx*dx + dy*dy); |
811 |
|
|
radLen = AIR_MAX(radLen, len); |
812 |
|
|
dx = 0 - drc->center[0]; |
813 |
|
|
dy = sy-1 - drc->center[1]; |
814 |
|
|
len = sqrt(dx*dx + dy*dy); |
815 |
|
|
radLen = AIR_MAX(radLen, len); |
816 |
|
|
if (drc->verbose) { |
817 |
|
|
fprintf(stderr, "%s: radLen = %g\n", me, radLen); |
818 |
|
|
} |
819 |
|
|
|
820 |
|
|
/* determine clamping, if any */ |
821 |
|
|
if (drc->clampPerc[0] > 0.0 || drc->clampPerc[1] > 0.0) { |
822 |
|
|
Nrrd *nhist; |
823 |
|
|
double *hist, total, sum; |
824 |
|
|
unsigned int hi; |
825 |
|
|
nhist = nrrdNew(); |
826 |
|
|
airMopAdd(mop, nhist, (airMopper)nrrdNuke, airMopAlways); |
827 |
|
|
if (nrrdHisto(nhist, drc->nin, NULL, NULL, drc->clampHistoBins, |
828 |
|
|
nrrdTypeDouble)) { |
829 |
|
|
biffAddf(NRRD, "%s: trouble making histogram", me); |
830 |
|
|
return 1; |
831 |
|
|
} |
832 |
|
|
hist = AIR_CAST(double *, nhist->data); |
833 |
|
|
total = AIR_CAST(double, nrrdElementNumber(drc->nin)); |
834 |
|
|
sum = 0; |
835 |
|
|
for (hi=0; hi<drc->clampHistoBins; hi++) { |
836 |
|
|
sum += hist[hi]; |
837 |
|
|
if (sum >= drc->clampPerc[0]*total/100.0) { |
838 |
|
|
drc->clamp[0] = AIR_AFFINE(0, hi, drc->clampHistoBins-1, |
839 |
|
|
nhist->axis[0].min, nhist->axis[0].max); |
840 |
|
|
break; |
841 |
|
|
} |
842 |
|
|
} |
843 |
|
|
if (hi == drc->clampHistoBins) { |
844 |
|
|
biffAddf(NRRD, "%s: failed to find lower %g-percentile value", me, |
845 |
|
|
drc->clampPerc[0]); |
846 |
|
|
return 1; |
847 |
|
|
} |
848 |
|
|
sum = 0; |
849 |
|
|
for (hi=drc->clampHistoBins; hi; hi--) { |
850 |
|
|
sum += hist[hi-1]; |
851 |
|
|
if (sum >= drc->clampPerc[1]*total/100.0) { |
852 |
|
|
drc->clamp[1] = AIR_AFFINE(0, hi-1, drc->clampHistoBins-1, |
853 |
|
|
nhist->axis[0].min, nhist->axis[0].max); |
854 |
|
|
break; |
855 |
|
|
} |
856 |
|
|
} |
857 |
|
|
if (!hi) { |
858 |
|
|
biffAddf(NRRD, "%s: failed to find upper %g-percentile value", me, |
859 |
|
|
drc->clampPerc[1]); |
860 |
|
|
return 1; |
861 |
|
|
} |
862 |
|
|
if (drc->verbose) { |
863 |
|
|
fprintf(stderr, "%s: [%g,%g]-%%ile clamping of [%g,%g] --> [%g,%g]\n", |
864 |
|
|
me, drc->clampPerc[0], drc->clampPerc[1], |
865 |
|
|
nhist->axis[0].min, nhist->axis[0].max, |
866 |
|
|
drc->clamp[0], drc->clamp[1]); |
867 |
|
|
} |
868 |
|
|
drc->clampDo = AIR_TRUE; |
869 |
|
|
} else { |
870 |
|
|
drc->clamp[0] = drc->clamp[1] = AIR_NAN; |
871 |
|
|
drc->clampDo = AIR_FALSE; |
872 |
|
|
} |
873 |
|
|
|
874 |
|
|
/* create deringBag(s) */ |
875 |
|
|
dbg = deringBagNew(drc, radLen); |
876 |
|
|
airMopAdd(mop, dbg, (airMopper)deringBagNix, airMopAlways); |
877 |
|
|
|
878 |
|
|
if (deringPtxfAlloc(drc, dbg)) { |
879 |
|
|
biffAddf(NRRD, "%s: trouble on setup", me); |
880 |
|
|
return 1; |
881 |
|
|
} |
882 |
|
|
|
883 |
|
|
sz = (2 == drc->nin->dim |
884 |
|
|
? 1 |
885 |
|
|
: AIR_CAST(unsigned int, drc->nin->axis[2].size)); |
886 |
|
|
drc->ringMagnitude = 0.0; |
887 |
|
|
for (zi=0; zi<sz; zi++) { |
888 |
|
|
if (drc->verbose) { |
889 |
|
|
fprintf(stderr, "%s: slice %u of %u ...\n", me, zi, sz); |
890 |
|
|
} |
891 |
|
|
if (deringDo(drc, dbg, nout, zi)) { |
892 |
|
|
biffAddf(NRRD, "%s: trouble on slice %u", me, zi); |
893 |
|
|
return 1; |
894 |
|
|
} |
895 |
|
|
drc->ringMagnitude += dbg->ringMag; |
896 |
|
|
if (drc->verbose) { |
897 |
|
|
fprintf(stderr, "%s: ... %u done\n", me, zi); |
898 |
|
|
} |
899 |
|
|
} |
900 |
|
|
if (drc->verbose) { |
901 |
|
|
fprintf(stderr, "%s: ring magnitude = %g\n", me, drc->ringMagnitude); |
902 |
|
|
} |
903 |
|
|
|
904 |
|
|
airMopOkay(mop); |
905 |
|
|
return 0; |
906 |
|
|
} |