GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/nrrd/deringNrrd.c Lines: 0 451 0.0 %
Date: 2017-05-26 Branches: 0 294 0.0 %

Line Branch Exec Source
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
}