GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/unrrdu/vidicon.c Lines: 23 127 18.1 %
Date: 2017-05-26 Branches: 1 90 1.1 %

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 "Try to create the look of early 80s analog B+W video"
28
static const char *_unrrdu_vidiconInfoL =
29
  (INFO
30
   ". Does various things, some more justified than others.\n "
31
   "* (as yet there's no single nrrd function which does all this)");
32
33
int
34
unrrdu_vidiconMain(int argc, const char **argv, const char *me,
35
                   hestParm *hparm) {
36
2
  hestOpt *opt = NULL;
37
  airArray *mop, *submop=NULL;
38
39
1
  unsigned int vsize[2], vpadding[2], rpadding[2];
40
1
  double rescale, rperc;
41
1
  Nrrd *nin, *nrescale, *npad, *nvbase, *ntmp, *nout;
42
1
  char *out, *err, *stpfx, stname[AIR_STRLEN_SMALL];
43
  int pret;
44
  NrrdResampleContext *rsmc;
45
1
  NrrdKernelSpec *rescaleKsp, *vdsmp[2];
46
  NrrdRange *b8range;
47
48
1
  hparm->elideSingleOtherDefault = AIR_FALSE;
49
50
1
  hestOptAdd(&opt, "i", "input", airTypeOther, 1, 1, &nin, NULL,
51
             "input image. Should be grayscale PNG.",
52
1
             NULL, NULL, nrrdHestNrrd);
53
1
  hestOptAdd(&opt, "rs", "rescale", airTypeDouble, 1, 1, &rescale, "0.75",
54
             "how to rescale (downsample) the image prior to processing, "
55
             "just to get a better representation of the floating-point "
56
             "range of image values (overcoming 8-bit quantization effects)");
57
1
  hestOptAdd(&opt, "rsk", "kern", airTypeOther, 1, 1, &rescaleKsp, "hann:5",
58
             "kernel for rescaling.",
59
1
             NULL, NULL, nrrdHestKernelSpec);
60
1
  hestOptAdd(&opt, "rsp", "percentile", airTypeDouble, 1, 1, &rperc, "1.5",
61
             "after rescaling, the highest and lowest percentiles are mapped "
62
             "to 0.0 and 255.0, just to have a uniform range of intensities "
63
             "in subsequent processing. This option determines how big those "
64
             "percentiles are.");
65
1
  hestOptAdd(&opt, "vs", "sx sy", airTypeUInt, 2, 2, vsize, "550 525",
66
             "the lowest (\"video\") resolution to which the image is "
67
             "down-sampled, reflecting the limited resolution of the "
68
             "vidicon tubes");
69
1
  hestOptAdd(&opt, "pad", "padX padY", airTypeUInt, 2, 2, vpadding, "10 10",
70
             "at the lowest resolution, there should be this much padding "
71
             "by black, to reflect the fact the signal outside the tube "
72
             "(e.g. between scanlines is black)");
73
2
  hestOptAdd(&opt, "vk", "kernX kernY", airTypeOther, 2, 2, vdsmp,
74
             "hann:1,4 cubic:0,0.5",
75
             "kernels for downsampling to video resolution; the horizontal "
76
             "and vertical kernels are different",
77
1
             NULL, NULL, nrrdHestKernelSpec);
78
1
  hestOptAdd(&opt, "stp", "prefix", airTypeString, 1, 1, &stpfx, "",
79
             "if a string is given here, a series of images are saved, "
80
             "representing the various stages of processing");
81
1
  hestOptAdd(&opt, "o", "output", airTypeString, 1, 1, &out, NULL,
82
             "output nrrd");
83
84
1
  mop = airMopNew();
85
1
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);
86
2
  USAGE(_unrrdu_vidiconInfoL);
87
  PARSE();
88
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);
89
  ntmp = nrrdNew();
90
  airMopAdd(mop, ntmp, (airMopper)nrrdNuke, airMopAlways);
91
  nout = nrrdNew();
92
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);
93
  b8range = nrrdRangeNew(0.0, 255.0);
94
  airMopAdd(mop, b8range, (airMopper)nrrdRangeNix, airMopAlways);
95
96
  if (!( 2 == nin->dim && nrrdTypeBlock != nin->type )) {
97
    fprintf(stderr, "%s: need input as 2D grayscale image (not %u-d %s)\n",
98
            me, nin->dim, airEnumStr(nrrdType, nin->type));
99
    airMopError(mop); return 1;
100
  }
101
  nrescale = nrrdNew();
102
  airMopAdd(mop, nrescale, (airMopper)nrrdNuke, airMopAlways);
103
104
  fprintf(stderr, "%s: rescaling by %g ... \n", me, rescale);
105
  rsmc = nrrdResampleContextNew();
106
  airMopAdd(mop, rsmc, (airMopper)nrrdResampleContextNix, airMopAlways);
107
  if (nrrdResampleDefaultCenterSet(rsmc, nrrdCenterCell)
108
      || nrrdResampleInputSet(rsmc, nin)
109
      || nrrdResampleKernelSet(rsmc, 0, rescaleKsp->kernel, rescaleKsp->parm)
110
      || nrrdResampleKernelSet(rsmc, 1, rescaleKsp->kernel, rescaleKsp->parm)
111
      || nrrdResampleSamplesSet(rsmc, 0, AIR_CAST(size_t,
112
                                                  rescale*nin->axis[0].size))
113
      || nrrdResampleSamplesSet(rsmc, 1, AIR_CAST(size_t,
114
                                                  rescale*nin->axis[1].size))
115
      || nrrdResampleRangeFullSet(rsmc, 0)
116
      || nrrdResampleRangeFullSet(rsmc, 1)
117
      || nrrdResampleTypeOutSet(rsmc, nrrdTypeFloat)
118
      || nrrdResampleRenormalizeSet(rsmc, AIR_TRUE)
119
      || nrrdResampleExecute(rsmc, nrescale)) {
120
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
121
    fprintf(stderr, "%s: problem rescaling:\n%s", me, err);
122
    airMopError(mop); return 1;
123
  }
124
125
#define SAVE_TMP(name, nrrd)                                            \
126
  if (airStrlen(stpfx)) {                                               \
127
    sprintf(stname, "%s-" #name ".png", stpfx);                         \
128
    if (nrrdQuantize(ntmp, nrrd, b8range, 8)                            \
129
        || nrrdSave(stname, ntmp, 0)) {                                 \
130
      airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);   \
131
      fprintf(stderr, "%s: problem saving %s:\n%s", me, stname, err);   \
132
      airMopError(mop); return 1;                                       \
133
    }                                                                   \
134
  }
135
  SAVE_TMP(rescale, nrescale);
136
137
  /* rescaling values to 0.0 -- 255.0 based on percentile rperc */
138
  {
139
    Nrrd *nhist;
140
    double *hist, sum, total, minval, maxval;
141
    unsigned int hi, hbins;
142
    float *rescaled;
143
    size_t ii, nn;
144
145
    submop = airMopNew();
146
    nhist = nrrdNew();
147
    hbins = 3000;
148
    airMopAdd(submop, nhist, (airMopper)nrrdNuke, airMopAlways);
149
    if (nrrdHisto(nhist, nrescale, NULL, NULL, hbins, nrrdTypeDouble)) {
150
      airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
151
      fprintf(stderr, "%s: trouble making histogram:\n%s", me, err);
152
      airMopError(submop); airMopError(mop); return 1;
153
    }
154
    hist = AIR_CAST(double *, nhist->data);
155
    total = AIR_CAST(double, nrrdElementNumber(nrescale));
156
    minval = AIR_NAN;
157
    sum = 0;
158
    for (hi=0; hi<hbins; hi++) {
159
      sum += hist[hi];
160
      if (sum >= rperc*total/100.0) {
161
        minval = AIR_AFFINE(0, hi, hbins-1,
162
                            nhist->axis[0].min, nhist->axis[0].max);
163
        break;
164
      }
165
    }
166
    if (hi == hbins || !AIR_EXISTS(minval)) {
167
      fprintf(stderr, "%s: failed to find lower %g-percentile value",
168
              me, rperc);
169
      airMopError(submop); airMopError(mop); return 1;
170
    }
171
    maxval = AIR_NAN;
172
    sum = 0;
173
    for (hi=hbins; hi; hi--) {
174
      sum += hist[hi-1];
175
      if (sum >= rperc*total/100.0) {
176
        maxval = AIR_AFFINE(0, hi-1, hbins-1,
177
                            nhist->axis[0].min, nhist->axis[0].max);
178
        break;
179
      }
180
    }
181
    if (!hi || !AIR_EXISTS(maxval)) {
182
      fprintf(stderr, "%s: failed to find upper %g-percentile value",
183
              me, rperc);
184
      airMopError(submop); airMopError(mop); return 1;
185
    }
186
    fprintf(stderr, "%s: min %g --> 0, max %g --> 255\n", me, minval, maxval);
187
    nn = nrrdElementNumber(nrescale);
188
    rescaled = AIR_CAST(float *, nrescale->data);
189
    for (ii=0; ii<nn; ii++) {
190
      rescaled[ii] = AIR_CAST(float, AIR_AFFINE(minval, rescaled[ii],
191
                                                maxval, 0.0, 255.0));
192
    }
193
    airMopOkay(submop);
194
    submop = NULL;
195
  }
196
197
  /* padding rescaled image with black */
198
  rpadding[0] = AIR_ROUNDUP(AIR_CAST(double, vpadding[0])
199
                            * nrescale->axis[0].size / vsize[0]);
200
  rpadding[1] = AIR_ROUNDUP(AIR_CAST(double, vpadding[1])
201
                            * nrescale->axis[1].size / vsize[1]);
202
  fprintf(stderr, "%s: padding in rescaled image: %u x %u\n",
203
          me, rpadding[0], rpadding[1]);
204
  npad = nrrdNew();
205
  airMopAdd(mop, npad, (airMopper)nrrdNuke, airMopAlways);
206
  {
207
    ptrdiff_t pmin[2], pmax[2];
208
    pmin[0] = -AIR_CAST(ptrdiff_t, rpadding[0]);
209
    pmin[1] = -AIR_CAST(ptrdiff_t, rpadding[1]);
210
    pmax[0] = nrescale->axis[0].size + rpadding[0];
211
    pmax[1] = nrescale->axis[1].size + rpadding[1];
212
    if (nrrdPad_nva(npad, nrescale, pmin, pmax, nrrdBoundaryPad, 0.0)) {
213
      airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
214
      fprintf(stderr, "%s: problem padding:\n%s", me, err);
215
      airMopError(mop); return 1;
216
    }
217
  }
218
219
  /* rescaling down to "video" resolution */
220
  fprintf(stderr, "%s: downsampling to %u x %u\n", me,
221
          AIR_CAST(unsigned int, vsize[0] + 2*vpadding[0]),
222
          AIR_CAST(unsigned int, vsize[1] + 2*vpadding[1]));
223
  nvbase = nrrdNew();
224
  airMopAdd(mop, nvbase, (airMopper)nrrdNuke, airMopAlways);
225
  if (nrrdResampleDefaultCenterSet(rsmc, nrrdCenterCell)
226
      || nrrdResampleInputSet(rsmc, npad)
227
      || nrrdResampleKernelSet(rsmc, 0, vdsmp[0]->kernel, vdsmp[0]->parm)
228
      || nrrdResampleKernelSet(rsmc, 1, vdsmp[1]->kernel, vdsmp[1]->parm)
229
      || nrrdResampleSamplesSet(rsmc, 0, vsize[0] + 2*vpadding[0])
230
      || nrrdResampleSamplesSet(rsmc, 1, vsize[1] + 2*vpadding[1])
231
      || nrrdResampleRangeFullSet(rsmc, 0)
232
      || nrrdResampleRangeFullSet(rsmc, 1)
233
      || nrrdResampleTypeOutSet(rsmc, nrrdTypeFloat)
234
      || nrrdResampleRenormalizeSet(rsmc, AIR_TRUE)
235
      || nrrdResampleExecute(rsmc, nvbase)) {
236
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
237
    fprintf(stderr, "%s: problem downsampling to video resolution:\n%s", me, err);
238
    airMopError(mop); return 1;
239
  }
240
241
  /* halo, lowfilt, windowing, noise, filt, interlace, noise, fuzz, upsample */
242
243
  nrrdCopy(nout, nvbase);
244
245
  SAVE(out, nout, NULL);
246
  airMopOkay(mop);
247
  return 0;
248
1
}
249
250
UNRRDU_CMD_HIDE(vidicon, INFO);