GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/unrrdu/cmedian.c Lines: 18 61 29.5 %
Date: 2017-05-26 Branches: 1 28 3.6 %

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 "unrrdu.h"
25
#include "privateUnrrdu.h"
26
27
#define INFO "Cheap histogram-based median/mode filtering"
28
static const char *_unrrdu_cmedianInfoL =
29
(INFO
30
 ". Only works on 1, 2, or 3 dimensions.  The window "
31
 "over which filtering is done is always square, and "
32
 "only a simplistic weighting scheme is available. "
33
 "The method is cheap because it does the median "
34
 "or mode based on a histogram, which enforces a quantization to the number "
35
 "of bins in the histogram, which probably means a loss of precision for "
36
 "anything except 8-bit data.  Also, integral values can be recovered "
37
 "exactly only when the number of bins is exactly min-max+1 (as reported "
38
 "by \"unu minmax\").\n "
39
 "* Uses nrrdCheapMedian, plus nrrdSlice and nrrdJoin in "
40
 "case of \"-c\"");
41
42
int
43
unrrdu_cmedianMain(int argc, const char **argv, const char *me,
44
                   hestParm *hparm) {
45
2
  hestOpt *opt = NULL;
46
1
  char *out, *err;
47
1
  Nrrd *nin, *nout, *ntmp, **mnout;
48
1
  int pad, pret, mode, chan, ni, nsize;
49
1
  unsigned int bins, radius;
50
  airArray *mop;
51
1
  float wght;
52
53
1
  hestOptAdd(&opt, "r,radius", "radius", airTypeUInt, 1, 1, &radius, NULL,
54
             "how big a window to filter over. \"-r 1\" leads to a "
55
             "3x3 window in an image, and a 3x3x3 window in a volume");
56
1
  hestOptAdd(&opt, "mode", NULL, airTypeInt, 0, 0, &mode, NULL,
57
             "By default, median filtering is done.  Using this option "
58
             "enables mode filtering, in which the most common value is "
59
             "used as output");
60
1
  hestOptAdd(&opt, "b,bins", "num", airTypeUInt, 1, 1, &bins, "256",
61
             "# of bins in histogram.  It is in your interest to minimize "
62
             "this number, since big histograms mean slower execution "
63
             "times.  8-bit data needs at most 256 bins.");
64
1
  hestOptAdd(&opt, "w,weight", "weight", airTypeFloat, 1, 1, &wght, "1.0",
65
             "How much higher to preferentially weight samples that are "
66
             "closer to the center of the window.  \"1.0\" weight means that "
67
             "all samples are uniformly weighted over the window, which "
68
             "facilitates a simple speed-up. ");
69
  /* FYI: these are weights which are just high enough to preserve
70
     an island of N contiguous high pixels in a row:
71
     1: 7.695
72
     2: 6.160
73
     3: 4.829 (actually only the middle pixel remains */
74
1
  hestOptAdd(&opt, "p,pad", NULL, airTypeInt, 0, 0, &pad, NULL,
75
             "Pad the input (with boundary method \"bleed\"), "
76
             "and crop the output, so as to "
77
             "overcome our cheapness and correctly "
78
             "handle the border.  Obviously, this takes more memory.");
79
1
  hestOptAdd(&opt, "c,channel", NULL, airTypeInt, 0, 0, &chan, NULL,
80
             "Slice the input along axis 0, run filtering on all slices, "
81
             "and join the results back together.  This is the way you'd "
82
             "want to process color (multi-channel) images or volumes.");
83
1
  OPT_ADD_NIN(nin, "input nrrd");
84
1
  OPT_ADD_NOUT(out, "output nrrd");
85
86
1
  mop = airMopNew();
87
1
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);
88
89
2
  USAGE(_unrrdu_cmedianInfoL);
90
  PARSE();
91
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);
92
93
  nout = nrrdNew();
94
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);
95
96
  if (chan) {
97
    nsize = AIR_UINT(nin->axis[0].size);
98
    mnout = AIR_CALLOC(nsize, Nrrd*);
99
    airMopAdd(mop, mnout, airFree, airMopAlways);
100
    ntmp = nrrdNew();
101
    airMopAdd(mop, ntmp, (airMopper)nrrdNuke, airMopAlways);
102
    for (ni=0; ni<nsize; ni++) {
103
      if (nrrdSlice(ntmp, nin, 0, ni)) {
104
        airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
105
        fprintf(stderr, "%s: error slicing input at pos = %d:\n%s",
106
                me, ni, err);
107
        airMopError(mop);
108
        return 1;
109
      }
110
      airMopAdd(mop, mnout[ni] = nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
111
      if (nrrdCheapMedian(mnout[ni], ntmp, pad, mode, radius, wght, bins)) {
112
        airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
113
        fprintf(stderr, "%s: error doing cheap median:\n%s", me, err);
114
        airMopError(mop);
115
        return 1;
116
      }
117
    }
118
    if (nrrdJoin(nout, (const Nrrd*const*)mnout, nsize, 0, AIR_TRUE)) {
119
      airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
120
      fprintf(stderr, "%s: error doing final join:\n%s", me, err);
121
      airMopError(mop);
122
      return 1;
123
    }
124
    nrrdAxisInfoCopy(nout, nin, NULL, NRRD_AXIS_INFO_NONE);
125
    if (nrrdBasicInfoCopy(nout, nin,
126
                          NRRD_BASIC_INFO_DATA_BIT
127
                          | NRRD_BASIC_INFO_TYPE_BIT
128
                          | NRRD_BASIC_INFO_BLOCKSIZE_BIT
129
                          | NRRD_BASIC_INFO_DIMENSION_BIT
130
                          | NRRD_BASIC_INFO_CONTENT_BIT
131
                          | NRRD_BASIC_INFO_COMMENTS_BIT
132
                          | (nrrdStateKeyValuePairsPropagate
133
                             ? 0
134
                             : NRRD_BASIC_INFO_KEYVALUEPAIRS_BIT))) {
135
      airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
136
      fprintf(stderr, "%s: error copying basic info:\n%s", me, err);
137
      airMopError(mop);
138
      return 1;
139
    }
140
  } else {
141
    if (nrrdCheapMedian(nout, nin, pad, mode, radius, wght, bins)) {
142
      airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
143
      fprintf(stderr, "%s: error doing cheap median:\n%s", me, err);
144
      airMopError(mop);
145
      return 1;
146
    }
147
  }
148
149
  SAVE(out, nout, NULL);
150
151
  airMopOkay(mop);
152
  return 0;
153
1
}
154
155
UNRRDU_CMD(cmedian, INFO);