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 "Create joint histogram of two or more nrrds" |
28 |
|
|
static const char *_unrrdu_jhistoInfoL = |
29 |
|
|
(INFO |
30 |
|
|
". Each axis of the output corresponds to one of the " |
31 |
|
|
"input nrrds, and each bin in the output records the " |
32 |
|
|
"number of corresponding positions in the inputs with " |
33 |
|
|
"a combination of values represented by the coordinates " |
34 |
|
|
"of the bin.\n " |
35 |
|
|
"* Uses nrrdHistoJoint"); |
36 |
|
|
|
37 |
|
|
int |
38 |
|
|
unrrdu_jhistoMain(int argc, const char **argv, const char *me, |
39 |
|
|
hestParm *hparm) { |
40 |
|
2 |
hestOpt *opt = NULL; |
41 |
|
1 |
char *out, *err; |
42 |
|
1 |
Nrrd **nin, **npass; |
43 |
|
1 |
Nrrd *nout, *nwght; |
44 |
|
1 |
size_t *bin; |
45 |
|
1 |
int type, clamp[NRRD_DIM_MAX], pret; |
46 |
|
1 |
unsigned int minLen, maxLen, ninLen, binLen, ai, diceax; |
47 |
|
|
airArray *mop; |
48 |
|
1 |
double *min, *max; |
49 |
|
|
NrrdRange **range; |
50 |
|
|
|
51 |
|
1 |
hestOptAdd(&opt, "b,bin", "bins0 bins1", airTypeSize_t, 2, -1, &bin, NULL, |
52 |
|
|
"bins<i> is the number of bins to use along axis i (of joint " |
53 |
|
|
"histogram), which represents the values of nin<i> ", |
54 |
|
|
&binLen); |
55 |
|
1 |
hestOptAdd(&opt, "w,weight", "nweight", airTypeOther, 1, 1, &nwght, "", |
56 |
|
|
"how to weigh contributions to joint histogram. By default " |
57 |
|
|
"(not using this option), the increment is one bin count per " |
58 |
|
|
"sample, but by giving a nrrd, the value in the nrrd at the " |
59 |
|
|
"corresponding location will be the bin count increment ", |
60 |
|
1 |
NULL, NULL, nrrdHestNrrd); |
61 |
|
1 |
hestOptAdd(&opt, "min,minimum", "min0 min1", airTypeDouble, 2, -1, |
62 |
|
|
&min, "nan nan", |
63 |
|
|
"min<i> is the low range of values to be quantized along " |
64 |
|
|
"axis i; use \"nan\" to represent lowest value present ", |
65 |
|
|
&minLen); |
66 |
|
1 |
hestOptAdd(&opt, "max,maximum", "max0 max1", airTypeDouble, 2, -1, |
67 |
|
|
&max, "nan nan", |
68 |
|
|
"max<i> is the high range of values to be quantized along " |
69 |
|
|
"axis i; use \"nan\" to represent highest value present ", |
70 |
|
|
&maxLen); |
71 |
|
1 |
OPT_ADD_TYPE(type, "type to use for output (the type used to store hit " |
72 |
|
|
"counts in the joint histogram). Clamping is done on hit " |
73 |
|
|
"counts so that they never overflow a fixed-point type", |
74 |
|
|
"uint"); |
75 |
|
1 |
hestOptAdd(&opt, "i,input", "nin0 [nin1]", airTypeOther, 1, -1, &nin, "-", |
76 |
|
|
"list of nrrds (one for each axis of joint histogram), " |
77 |
|
|
"or, single nrrd that will be sliced along specified axis.", |
78 |
|
1 |
&ninLen, NULL, nrrdHestNrrd); |
79 |
|
1 |
hestOptAdd(&opt, "a,axis", "axis", airTypeUInt, 1, 1, &diceax, "0", |
80 |
|
|
"axis to slice along when working with single nrrd. "); |
81 |
|
1 |
OPT_ADD_NOUT(out, "output nrrd"); |
82 |
|
|
|
83 |
|
1 |
mop = airMopNew(); |
84 |
|
1 |
airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways); |
85 |
|
|
|
86 |
✓✗ |
2 |
USAGE(_unrrdu_jhistoInfoL); |
87 |
|
|
PARSE(); |
88 |
|
|
airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways); |
89 |
|
|
|
90 |
|
|
if (ninLen == 1) { |
91 |
|
|
/* Slice a nrrd on the fly */ |
92 |
|
|
size_t asize; |
93 |
|
|
if (!( diceax <= nin[0]->dim-1 )) { |
94 |
|
|
fprintf(stderr, "%s: slice axis %u not valid for single %u-D nrrd", |
95 |
|
|
me, diceax, nin[0]->dim); |
96 |
|
|
airMopError(mop); |
97 |
|
|
return 1; |
98 |
|
|
} |
99 |
|
|
asize = nin[0]->axis[diceax].size; |
100 |
|
|
if (asize != binLen) { |
101 |
|
|
fprintf(stderr, |
102 |
|
|
"%s: size (%u) of slice axis %u != # bins given (%u)\n", me, |
103 |
|
|
AIR_CAST(unsigned int, asize), diceax, |
104 |
|
|
AIR_CAST(unsigned int, binLen)); |
105 |
|
|
airMopError(mop); |
106 |
|
|
return 1; |
107 |
|
|
} |
108 |
|
|
/* create buffer for slices */ |
109 |
|
|
if (!( npass = AIR_CALLOC(binLen, Nrrd*) )) { |
110 |
|
|
fprintf(stderr, "%s: couldn't allocate nrrd array (size %u)\n", |
111 |
|
|
me, binLen); |
112 |
|
|
airMopError(mop); return 1; |
113 |
|
|
} |
114 |
|
|
airMopMem(mop, &npass, airMopAlways); |
115 |
|
|
/* slice this nrrd, allocate new nrrds, and store the slices in nin */ |
116 |
|
|
for (ai=0; ai<binLen; ai++) { |
117 |
|
|
/* Allocate each slice nrrd */ |
118 |
|
|
if (!( npass[ai] = nrrdNew() )) { |
119 |
|
|
fprintf(stderr, "%s: couldn't allocate npass[%u]\n", me, ai); |
120 |
|
|
airMopError(mop); return 1; |
121 |
|
|
} |
122 |
|
|
airMopAdd(mop, npass[ai], (airMopper)nrrdNuke, airMopAlways); |
123 |
|
|
if (nrrdSlice(npass[ai], nin[0], diceax, ai)) { |
124 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
125 |
|
|
fprintf(stderr, "%s: error slicing:\n%s", me, err); |
126 |
|
|
airMopError(mop); return 1; |
127 |
|
|
} |
128 |
|
|
} |
129 |
|
|
} else { |
130 |
|
|
/* we were given multiple nrrds */ |
131 |
|
|
if (ninLen != binLen) { |
132 |
|
|
fprintf(stderr, |
133 |
|
|
"%s: # input nrrds (%u) != # bin specifications (%u)\n", me, |
134 |
|
|
AIR_CAST(unsigned int, ninLen), AIR_CAST(unsigned int, binLen)); |
135 |
|
|
airMopError(mop); |
136 |
|
|
return 1; |
137 |
|
|
} |
138 |
|
|
/* create buffer for slices (HEY copy and paste) */ |
139 |
|
|
if (!( npass = AIR_CALLOC(binLen, Nrrd*) )) { |
140 |
|
|
fprintf(stderr, "%s: couldn't allocate nrrd array (size %u)\n", |
141 |
|
|
me, binLen); |
142 |
|
|
airMopError(mop); return 1; |
143 |
|
|
} |
144 |
|
|
for (ai=0; ai<binLen; ai++) { |
145 |
|
|
npass[ai] = nin[ai]; |
146 |
|
|
} |
147 |
|
|
} |
148 |
|
|
range = AIR_CALLOC(binLen, NrrdRange*); |
149 |
|
|
airMopAdd(mop, range, airFree, airMopAlways); |
150 |
|
|
for (ai=0; ai<binLen; ai++) { |
151 |
|
|
range[ai] = nrrdRangeNew(AIR_NAN, AIR_NAN); |
152 |
|
|
airMopAdd(mop, range[ai], (airMopper)nrrdRangeNix, airMopAlways); |
153 |
|
|
} |
154 |
|
|
if (2 != minLen || (AIR_EXISTS(min[0]) || AIR_EXISTS(min[1]))) { |
155 |
|
|
if (minLen != binLen) { |
156 |
|
|
fprintf(stderr, "%s: # mins (%d) != # input nrrds (%d)\n", me, |
157 |
|
|
minLen, binLen); |
158 |
|
|
airMopError(mop); return 1; |
159 |
|
|
} |
160 |
|
|
for (ai=0; ai<binLen; ai++) { |
161 |
|
|
range[ai]->min = min[ai]; |
162 |
|
|
} |
163 |
|
|
} |
164 |
|
|
if (2 != maxLen || (AIR_EXISTS(max[0]) || AIR_EXISTS(max[1]))) { |
165 |
|
|
if (maxLen != binLen) { |
166 |
|
|
fprintf(stderr, "%s: # maxs (%d) != # input nrrds (%d)\n", me, |
167 |
|
|
maxLen, binLen); |
168 |
|
|
airMopError(mop); return 1; |
169 |
|
|
} |
170 |
|
|
for (ai=0; ai<binLen; ai++) { |
171 |
|
|
range[ai]->max = max[ai]; |
172 |
|
|
} |
173 |
|
|
} |
174 |
|
|
for (ai=0; ai<binLen; ai++) { |
175 |
|
|
clamp[ai] = 0; |
176 |
|
|
} |
177 |
|
|
|
178 |
|
|
nout = nrrdNew(); |
179 |
|
|
airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways); |
180 |
|
|
|
181 |
|
|
if (nrrdHistoJoint(nout, (const Nrrd*const*)npass, |
182 |
|
|
(const NrrdRange*const*)range, |
183 |
|
|
binLen, nwght, bin, type, clamp)) { |
184 |
|
|
airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways); |
185 |
|
|
fprintf(stderr, "%s: error doing joint histogram:\n%s", me, err); |
186 |
|
|
airMopError(mop); |
187 |
|
|
return 1; |
188 |
|
|
} |
189 |
|
|
|
190 |
|
|
SAVE(out, nout, NULL); |
191 |
|
|
|
192 |
|
|
airMopOkay(mop); |
193 |
|
|
return 0; |
194 |
|
1 |
} |
195 |
|
|
|
196 |
|
|
UNRRDU_CMD(jhisto, INFO); |