GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/ten/tendEllipse.c Lines: 24 189 12.7 %
Date: 2017-05-26 Branches: 1 84 1.2 %

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 "ten.h"
25
#include "privateTen.h"
26
27
#define INFO "Generate postscript renderings of 2D glyphs"
28
static const char *_tend_ellipseInfoL =
29
  (INFO
30
   ".  Not much to look at here.");
31
32
int
33
tend_ellipseDoit(FILE *file, Nrrd *nten, Nrrd *npos, Nrrd *nstn,
34
                 float min[2], float max[2],
35
                 float gscale, float dotRad, float lineWidth, float cthresh,
36
                 int invert) {
37
  size_t sx=0, sy=0, ti, nt;
38
  int x, y, vi, *sdata;
39
  double aspect, minX, minY, maxX, maxY, conf, Dxx, Dxy, Dyy, px, py, spx, spy;
40
  float *tdata, *pdata;
41
42
  if (npos) {
43
    nt = npos->axis[1].size;
44
    aspect = (max[0] - min[0])/(max[1] - min[1]);
45
  } else {
46
    spx = (AIR_EXISTS(nten->axis[1].spacing)
47
           ? nten->axis[1].spacing
48
           : 1);
49
    spy = (AIR_EXISTS(nten->axis[2].spacing)
50
           ? nten->axis[2].spacing
51
           : 1);
52
    sx = nten->axis[1].size;
53
    sy = nten->axis[2].size;
54
    nt = sx*sy;
55
    aspect = sx*spx/(sy*spy);
56
  }
57
58
  if (aspect > 7.5/10) {
59
    /* image has a wider aspect ratio than safely printable page area */
60
    minX = 0.5;
61
    maxX = 8.0;
62
    minY = 5.50 - 7.5/2/aspect;
63
    maxY = 5.50 + 7.5/2/aspect;
64
  } else {
65
    /* image is taller ... */
66
    minX = 4.25 - 10.0/2*aspect;
67
    maxX = 4.25 + 10.0/2*aspect;
68
    minY = 0.5;
69
    maxY = 10.5;
70
  }
71
  minX *= 72; minY *= 72;
72
  maxX *= 72; maxY *= 72;
73
  if (npos) {
74
    gscale *= AIR_CAST(float, (maxX - minX)/(max[0] - min[0]));
75
    dotRad *= AIR_CAST(float, (maxX - minX)/(max[0] - min[0]));
76
    lineWidth *= AIR_CAST(float, (maxX - minX)/(max[0] - min[0]));
77
  }
78
79
  fprintf(file, "%%!PS-Adobe-3.0 EPSF-3.0\n");
80
  fprintf(file, "%%%%Creator: tend ellipse\n");
81
  fprintf(file, "%%%%Title: blah blah blah\n");
82
  fprintf(file, "%%%%Pages: 1\n");
83
  fprintf(file, "%%%%BoundingBox: %d %d %d %d\n",
84
          AIR_CAST(int, floor(minX)), AIR_CAST(int, floor(minY)),
85
          AIR_CAST(int, ceil(maxX)), AIR_CAST(int, ceil(maxY)));
86
  fprintf(file, "%%%%HiResBoundingBox: %g %g %g %g\n",
87
          minX, minY, maxX, maxY);
88
  fprintf(file, "%%%%EndComments\n");
89
  fprintf(file, "%%%%BeginProlog\n");
90
  fprintf(file, "%%%%EndProlog\n");
91
  fprintf(file, "%%%%Page: 1 1\n");
92
93
  fprintf(file, "gsave\n");
94
95
  if (invert) {
96
    fprintf(file, "0 setgray\n");
97
    fprintf(file, "%g %g moveto\n", minX, minY);
98
    fprintf(file, "%g %g lineto\n", maxX, minY);
99
    fprintf(file, "%g %g lineto\n", maxX, maxY);
100
    fprintf(file, "%g %g lineto\n", minX, maxY);
101
    fprintf(file, "closepath fill\n");
102
  }
103
104
  fprintf(file, "gsave\n");
105
  fprintf(file, "0.5 setgray\n");
106
  tdata = (float*)nten->data;
107
  pdata = npos ? (float*)npos->data : NULL;
108
  for (ti=0; ti<nt; ti++) {
109
    if (npos) {
110
      px = AIR_AFFINE(min[0], pdata[0], max[0], minX, maxX);
111
      py = AIR_AFFINE(min[1], pdata[1], max[1], maxY, minY);
112
      pdata += 2;
113
    } else {
114
      x = ti % sx;
115
      y = ti / sx;
116
      px = NRRD_CELL_POS(minX, maxX, sx, x);
117
      py = NRRD_CELL_POS(minY, maxY, sy, sy-1-y);
118
    }
119
    conf = tdata[0];
120
    if (conf > cthresh) {
121
      double eval0, eval1, dd;
122
      Dxx = tdata[1];
123
      Dxy = tdata[2];
124
      Dyy = tdata[3];
125
      dd = Dxx - Dyy;
126
      eval0 = 0.5*(-Dxx + sqrt(4*Dxy*Dxy + dd*dd) - Dyy);
127
      eval1 = 0.5*(-Dxx - sqrt(4*Dxy*Dxy + dd*dd) - Dyy);
128
      fprintf(file, "gsave\n");
129
      fprintf(file, "matrix currentmatrix\n");
130
      fprintf(file, "[%g %g %g %g %g %g] concat\n",
131
              Dxx, -Dxy, -Dxy, Dyy, px, py);
132
      fprintf(file, "0 0 %g 0 360 arc closepath\n", gscale);
133
      fprintf(file, "setmatrix\n");
134
      if (eval0 * eval1 < 0) {
135
        fprintf(file, "gsave\n");
136
        fprintf(file, "0.15 setgray\n");
137
        fprintf(file, "fill\n");
138
        fprintf(file, "grestore\n");
139
      } else {
140
        fprintf(file, "fill\n");
141
      }
142
      fprintf(file, "grestore\n");
143
    }
144
    tdata += 4;
145
  }
146
  fprintf(file, "grestore\n");
147
148
  if (dotRad && !nstn) {
149
    fprintf(file, "gsave\n");
150
    tdata = (float*)nten->data;
151
    pdata = npos ? (float*)npos->data : NULL;
152
    fprintf(file, "%g setgray\n", invert ? 1.0 : 0.0);
153
    for (ti=0; ti<nt; ti++) {
154
      if (npos) {
155
        px = AIR_AFFINE(min[0], pdata[0], max[0], minX, maxX);
156
        py = AIR_AFFINE(min[1], pdata[1], max[1], maxY, minY);
157
        pdata += 2;
158
      } else {
159
        x = ti % sx;
160
        y = ti / sx;
161
        px = NRRD_CELL_POS(minX, maxX, sx, x);
162
        py = NRRD_CELL_POS(minY, maxY, sy, sy-1-y);
163
      }
164
      conf = tdata[0];
165
      if (conf > cthresh) {
166
        fprintf(file, "%g %g %g 0 360 arc closepath fill\n", px, py, dotRad);
167
      }
168
      tdata += 4;
169
    }
170
    fprintf(file, "grestore\n");
171
  }
172
173
  if ((dotRad || lineWidth) && npos && nstn) {
174
    fprintf(file, "gsave\n");
175
    tdata = (float*)nten->data;
176
    pdata = npos ? (float*)npos->data : NULL;
177
    sdata = nstn ? (int*)nstn->data : NULL;
178
    fprintf(file, "%g setlinewidth\n", lineWidth);
179
    fprintf(file, "%g setgray\n", invert ? 1.0 : 0.0);
180
    fprintf(file, "1 setlinecap\n");
181
    fprintf(file, "1 setlinejoin\n");
182
    for (ti=0; ti<nstn->axis[1].size; ti++) {
183
      if (1 == sdata[1 + 3*ti]) {
184
        vi = sdata[0 + 3*ti];
185
        px = AIR_AFFINE(min[0], pdata[0 + 2*vi], max[0], minX, maxX);
186
        py = AIR_AFFINE(min[1], pdata[1 + 2*vi], max[1], maxY, minY);
187
        if (tdata[0 + 4*vi] > cthresh) {
188
          fprintf(file, "%g %g %g 0 360 arc closepath fill\n", px, py, dotRad);
189
        }
190
      } else {
191
        fprintf(file, "newpath\n");
192
        for (vi = sdata[0 + 3*ti];
193
             vi < sdata[0 + 3*ti] + sdata[1 + 3*ti];
194
             vi++) {
195
          px = AIR_AFFINE(min[0], pdata[0 + 2*vi], max[0], minX, maxX);
196
          py = AIR_AFFINE(min[1], pdata[1 + 2*vi], max[1], maxY, minY);
197
          fprintf(file, "%g %g %s\n", px, py,
198
                  vi == sdata[0 + 3*ti] ? "moveto" : "lineto");
199
        }
200
        fprintf(file, "stroke\n");
201
        vi = sdata[0 + 3*ti] + sdata[2 + 3*ti];
202
        px = AIR_AFFINE(min[0], pdata[0 + 2*vi], max[0], minX, maxX);
203
        py = AIR_AFFINE(min[1], pdata[1 + 2*vi], max[1], maxY, minY);
204
        fprintf(file, "%g %g %g 0 360 arc closepath fill\n",
205
                px, py, dotRad + lineWidth);
206
      }
207
    }
208
    fprintf(file, "grestore\n");
209
  }
210
211
  fprintf(file, "grestore\n");
212
213
  return 0;
214
}
215
216
int
217
tend_ellipseMain(int argc, const char **argv, const char *me,
218
                 hestParm *hparm) {
219
  int pret;
220
2
  hestOpt *hopt = NULL;
221
1
  char *perr;
222
  airArray *mop;
223
224
1
  Nrrd *nten, *npos, *nstn;
225
1
  char *outS;
226
1
  float gscale, dotRad, lineWidth, cthresh, min[2], max[2];
227
  FILE *fout;
228
1
  int invert;
229
230
1
  mop = airMopNew();
231
232
1
  hestOptAdd(&hopt, "ctr", "conf thresh", airTypeFloat, 1, 1, &cthresh, "0.5",
233
             "Glyphs will be drawn only for tensors with confidence "
234
             "values greater than this threshold");
235
1
  hestOptAdd(&hopt, "gsc", "scale", airTypeFloat, 1, 1, &gscale, "1",
236
             "over-all glyph size");
237
1
  hestOptAdd(&hopt, "dot", "radius", airTypeFloat, 1, 1, &dotRad, "0.0",
238
             "radius of little dot to put in middle of ellipse, or \"0\" "
239
             "for no such dot");
240
1
  hestOptAdd(&hopt, "wid", "width", airTypeFloat, 1, 1, &lineWidth, "0.0",
241
             "with of lines for tractlets");
242
1
  hestOptAdd(&hopt, "inv", NULL, airTypeInt, 0, 0, &invert, NULL,
243
             "use white ellipses on black background, instead of reverse");
244
1
  hestOptAdd(&hopt, "min", "minX minY", airTypeFloat, 2, 2, min, "-1 -1",
245
             "when using \"-p\", minimum corner");
246
1
  hestOptAdd(&hopt, "max", "maxX maxY", airTypeFloat, 2, 2, max, "1 1",
247
             "when using \"-p\", maximum corner");
248
249
  /* input/output */
250
1
  hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &nten, "-",
251
1
             "image of 2D tensors", NULL, NULL, nrrdHestNrrd);
252
1
  hestOptAdd(&hopt, "p", "pos array", airTypeOther, 1, 1, &npos, "",
253
             "Instead of being on a grid, tensors are at arbitrary locations, "
254
             "as defined by this 2-by-N array of floats", NULL, NULL,
255
1
             nrrdHestNrrd);
256
1
  hestOptAdd(&hopt, "s", "stn array", airTypeOther, 1, 1, &nstn, "",
257
             "Locations given by \"-p\" have this connectivity", NULL, NULL,
258
1
             nrrdHestNrrd);
259
1
  hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-",
260
             "output PostScript file");
261
262
1
  airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
263
2
  USAGE(_tend_ellipseInfoL);
264
  JUSTPARSE();
265
  airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
266
267
  if (npos) {
268
    if (!( 2 == nten->dim && 4 == nten->axis[0].size
269
           && 2 == npos->dim && 2 == npos->axis[0].size
270
           && nten->axis[1].size == npos->axis[1].size )) {
271
      fprintf(stderr, "%s: didn't get matching lists of tensors and pos's\n",
272
              me);
273
      airMopError(mop); return 1;
274
    }
275
    if (!( nrrdTypeFloat == npos->type )) {
276
      fprintf(stderr, "%s: didn't get float type positions\n", me);
277
      airMopError(mop); return 1;
278
    }
279
  } else {
280
    if (!(3 == nten->dim && 4 == nten->axis[0].size)) {
281
      fprintf(stderr, "%s: didn't get a 3-D 4-by-X-by-Y 2D tensor array\n",
282
              me);
283
      airMopError(mop); return 1;
284
    }
285
  }
286
  if (!( nrrdTypeFloat == nten->type )) {
287
    fprintf(stderr, "%s: didn't get float type tensors\n", me);
288
    airMopError(mop); return 1;
289
  }
290
  if (nstn) {
291
    if (!( nrrdTypeUInt == nstn->type
292
           && 2 == nstn->dim
293
           && 3 == nstn->axis[0].size )) {
294
      fprintf(stderr, "%s: connectivity isn't 2-D 3-by-N array of %ss\n",
295
              me, airEnumStr(nrrdType, nrrdTypeInt));
296
      airMopError(mop); return 1;
297
    }
298
  }
299
  if (!(fout = airFopen(outS, stdout, "wb"))) {
300
    fprintf(stderr, "%s: couldn't open \"%s\" for writing\n", me, outS);
301
    airMopError(mop); return 1;
302
  }
303
  airMopAdd(mop, fout, (airMopper)airFclose, airMopAlways);
304
305
  tend_ellipseDoit(fout, nten, npos, nstn, min, max,
306
                   gscale, dotRad, lineWidth, cthresh, invert);
307
308
  airMopOkay(mop);
309
  return 0;
310
1
}
311
TEND_CMD(ellipse, INFO);