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 |
|
|
/* |
28 |
|
|
** learned: if you have globals, such as _nrrdCC_verb, which are |
29 |
|
|
** defined and declared here, but which are NOT initialized, then |
30 |
|
|
** C++ apps which are linking against Teem will have problems!!! |
31 |
|
|
** This was first seen on the mac. |
32 |
|
|
*/ |
33 |
|
|
int _nrrdCC_verb = 0; |
34 |
|
|
int _nrrdCC_EqvIncr = 10000; /* HEY: this has to be big so that ccfind is not |
35 |
|
|
stuck constantly re-allocating the eqv array. |
36 |
|
|
This is will be one of the best places to |
37 |
|
|
test the new multiplicative reallocation |
38 |
|
|
strategy, planned for Teem 2.0 */ |
39 |
|
|
|
40 |
|
|
int |
41 |
|
|
_nrrdCCFind_1(Nrrd *nout, unsigned int *numid, const Nrrd *nin) { |
42 |
|
|
/* static const char me[]="_nrrdCCFind_1"; */ |
43 |
|
|
unsigned int sx, I, id, lval, val, *out, (*lup)(const void *, size_t); |
44 |
|
|
|
45 |
|
|
lup = nrrdUILookup[nin->type]; |
46 |
|
|
out = AIR_CAST(unsigned int*, nout->data); |
47 |
|
|
out[0] = id = 0; |
48 |
|
|
*numid = 1; |
49 |
|
|
|
50 |
|
|
sx = AIR_CAST(unsigned int, nin->axis[0].size); |
51 |
|
|
lval = lup(nin->data, 0); |
52 |
|
|
for (I=1; I<sx; I++) { |
53 |
|
|
val = lup(nin->data, I); |
54 |
|
|
if (lval != val) { |
55 |
|
|
id++; |
56 |
|
|
(*numid)++; |
57 |
|
|
} |
58 |
|
|
out[I] = id; |
59 |
|
|
lval = val; |
60 |
|
|
} |
61 |
|
|
|
62 |
|
|
return 0; |
63 |
|
|
} |
64 |
|
|
|
65 |
|
|
/* |
66 |
|
|
** layout of value (pvl) and index (pid) cache: |
67 |
|
|
** |
68 |
|
|
** 2 3 4 --> X |
69 |
|
|
** 1 . . oddly, index 0 is never used |
70 |
|
|
** . . . |
71 |
|
|
** | |
72 |
|
|
** v Y |
73 |
|
|
*/ |
74 |
|
|
int |
75 |
|
|
_nrrdCCFind_2(Nrrd *nout, unsigned int *numid, airArray *eqvArr, |
76 |
|
|
const Nrrd *nin, unsigned int conny) { |
77 |
|
|
static const char me[]="_nrrdCCFind_2"; |
78 |
|
|
double vl=0, pvl[5]={0,0,0,0,0}; |
79 |
|
|
unsigned int id, pid[5]={0,0,0,0,0}, (*lup)(const void *, size_t), *out; |
80 |
|
|
unsigned int p, x, y, sx, sy; |
81 |
|
|
|
82 |
|
|
id = 0; /* sssh! compiler warnings */ |
83 |
|
|
lup = nrrdUILookup[nin->type]; |
84 |
|
|
out = AIR_CAST(unsigned int*, nout->data); |
85 |
|
|
sx = AIR_CAST(unsigned int, nin->axis[0].size); |
86 |
|
|
sy = AIR_CAST(unsigned int, nin->axis[1].size); |
87 |
|
|
#define GETV_2(x,y) ((AIR_IN_CL(0, AIR_CAST(int, x), AIR_CAST(int, sx-1)) \ |
88 |
|
|
&& AIR_IN_CL(0, AIR_CAST(int, y), AIR_CAST(int, sy-1))) \ |
89 |
|
|
? lup(nin->data, (x) + sx*(y)) \ |
90 |
|
|
: 0.5) /* value that can't come from an array of uints */ |
91 |
|
|
#define GETI_2(x,y) ((AIR_IN_CL(0, AIR_CAST(int, x), AIR_CAST(int, sx-1)) \ |
92 |
|
|
&& AIR_IN_CL(0, AIR_CAST(int, y), AIR_CAST(int, sy-1))) \ |
93 |
|
|
? out[(x) + sx*(y)] \ |
94 |
|
|
: AIR_CAST(unsigned int, -1)) /* CC index (probably!) |
95 |
|
|
never assigned */ |
96 |
|
|
|
97 |
|
|
*numid = 0; |
98 |
|
|
for (y=0; y<sy; y++) { |
99 |
|
|
for (x=0; x<sx; x++) { |
100 |
|
|
if (_nrrdCC_verb) { |
101 |
|
|
fprintf(stderr, "%s(%d,%d) -----------\n", me, x, y); |
102 |
|
|
} |
103 |
|
|
if (!x) { |
104 |
|
|
pvl[1] = GETV_2(-1, y); pid[1] = GETI_2(-1, y); |
105 |
|
|
pvl[2] = GETV_2(-1, y-1); pid[2] = GETI_2(-1, y-1); |
106 |
|
|
pvl[3] = GETV_2(0, y-1); pid[3] = GETI_2(0, y-1); |
107 |
|
|
|
108 |
|
|
} else { |
109 |
|
|
pvl[1] = vl; pid[1] = id; |
110 |
|
|
pvl[2] = pvl[3]; pid[2] = pid[3]; |
111 |
|
|
pvl[3] = pvl[4]; pid[3] = pid[4]; |
112 |
|
|
} |
113 |
|
|
pvl[4] = GETV_2(x+1, y-1); pid[4] = GETI_2(x+1, y-1); |
114 |
|
|
vl = GETV_2(x, y); |
115 |
|
|
p = 0; |
116 |
|
|
if (vl == pvl[1]) { |
117 |
|
|
id = pid[p=1]; |
118 |
|
|
} |
119 |
|
|
#define TEST(P) \ |
120 |
|
|
if (vl == pvl[(P)]) { \ |
121 |
|
|
if (p) { /* we already had a value match */ \ |
122 |
|
|
if (id != pid[(P)]) { \ |
123 |
|
|
airEqvAdd(eqvArr, pid[(P)], id); \ |
124 |
|
|
} \ |
125 |
|
|
} else { \ |
126 |
|
|
id = pid[p=(P)]; \ |
127 |
|
|
} \ |
128 |
|
|
} |
129 |
|
|
TEST(3); |
130 |
|
|
if (2 == conny) { |
131 |
|
|
TEST(2); |
132 |
|
|
TEST(4); |
133 |
|
|
} |
134 |
|
|
if (!p) { |
135 |
|
|
/* didn't match anything previous */ |
136 |
|
|
id = *numid; |
137 |
|
|
(*numid)++; |
138 |
|
|
} |
139 |
|
|
if (_nrrdCC_verb) { |
140 |
|
|
fprintf(stderr, "%s: pvl: %g %g %g %g (vl = %g)\n", me, |
141 |
|
|
pvl[1], pvl[2], pvl[3], pvl[4], vl); |
142 |
|
|
fprintf(stderr, " pid: %d %d %d %d\n", |
143 |
|
|
pid[1], pid[2], pid[3], pid[4]); |
144 |
|
|
fprintf(stderr, " --> p = %d, id = %d, *numid = %d\n", |
145 |
|
|
p, id, *numid); |
146 |
|
|
} |
147 |
|
|
out[x + sx*y] = id; |
148 |
|
|
} |
149 |
|
|
} |
150 |
|
|
|
151 |
|
|
return 0; |
152 |
|
|
} |
153 |
|
|
|
154 |
|
|
/* |
155 |
|
|
** |
156 |
|
|
** 5 6 7 --> X |
157 |
|
|
** 8 9 10 |
158 |
|
|
** 11 12 13 |
159 |
|
|
** | |
160 |
|
|
** v Y |
161 |
|
|
** 2 3 4 |
162 |
|
|
** / 1 . . again, 0 index never used, for reasons forgotten |
163 |
|
|
** Z . . . |
164 |
|
|
*/ |
165 |
|
|
int |
166 |
|
|
_nrrdCCFind_3(Nrrd *nout, unsigned int *numid, airArray *eqvArr, |
167 |
|
|
const Nrrd *nin, unsigned int conny) { |
168 |
|
|
/* static const char me[]="_nrrdCCFind_3" ; */ |
169 |
|
|
double pvl[14]={0,0,0,0,0,0,0,0,0,0,0,0,0,0}, vl=0; |
170 |
|
|
unsigned int id, *out, (*lup)(const void *, size_t), |
171 |
|
|
pid[14]={0,0,0,0,0,0,0,0,0,0,0,0,0,0}; |
172 |
|
|
unsigned int p, x, y, z, sx, sy, sz; |
173 |
|
|
|
174 |
|
|
id = 0; /* sssh! compiler warnings */ |
175 |
|
|
lup = nrrdUILookup[nin->type]; |
176 |
|
|
out = AIR_CAST(unsigned int*, nout->data); |
177 |
|
|
sx = AIR_CAST(unsigned int, nin->axis[0].size); |
178 |
|
|
sy = AIR_CAST(unsigned int, nin->axis[1].size); |
179 |
|
|
sz = AIR_CAST(unsigned int, nin->axis[2].size); |
180 |
|
|
#define GETV_3(x,y,z) ((AIR_IN_CL(0, AIR_CAST(int, x), AIR_CAST(int, sx-1)) \ |
181 |
|
|
&& AIR_IN_CL(0, AIR_CAST(int, y), AIR_CAST(int, sy-1)) \ |
182 |
|
|
&& AIR_IN_CL(0, AIR_CAST(int, z), AIR_CAST(int, sz-1)))\ |
183 |
|
|
? lup(nin->data, (x) + sx*((y) + sy*(z))) \ |
184 |
|
|
: 0.5) |
185 |
|
|
#define GETI_3(x,y,z) ((AIR_IN_CL(0, AIR_CAST(int, x), AIR_CAST(int, sx-1)) \ |
186 |
|
|
&& AIR_IN_CL(0, AIR_CAST(int, y), AIR_CAST(int, sy-1)) \ |
187 |
|
|
&& AIR_IN_CL(0, AIR_CAST(int, z), AIR_CAST(int, sz-1)))\ |
188 |
|
|
? out[(x) + sx*((y) + sy*(z))] \ |
189 |
|
|
: AIR_CAST(unsigned int, -1)) |
190 |
|
|
|
191 |
|
|
*numid = 0; |
192 |
|
|
for (z=0; z<sz; z++) { |
193 |
|
|
for (y=0; y<sy; y++) { |
194 |
|
|
for (x=0; x<sx; x++) { |
195 |
|
|
if (!x) { |
196 |
|
|
pvl[ 1] = GETV_3( -1, y, z); pid[ 1] = GETI_3( -1, y, z); |
197 |
|
|
pvl[ 2] = GETV_3( -1, y-1, z); pid[ 2] = GETI_3( -1, y-1, z); |
198 |
|
|
pvl[ 3] = GETV_3( 0, y-1, z); pid[ 3] = GETI_3( 0, y-1, z); |
199 |
|
|
pvl[ 5] = GETV_3( -1, y-1, z-1); pid[ 5] = GETI_3( -1, y-1, z-1); |
200 |
|
|
pvl[ 8] = GETV_3( -1, y, z-1); pid[ 8] = GETI_3( -1, y, z-1); |
201 |
|
|
pvl[11] = GETV_3( -1, y+1, z-1); pid[11] = GETI_3( -1, y+1, z-1); |
202 |
|
|
pvl[ 6] = GETV_3( 0, y-1, z-1); pid[ 6] = GETI_3( 0, y-1, z-1); |
203 |
|
|
pvl[ 9] = GETV_3( 0, y, z-1); pid[ 9] = GETI_3( 0, y, z-1); |
204 |
|
|
pvl[12] = GETV_3( 0, y+1, z-1); pid[12] = GETI_3( 0, y+1, z-1); |
205 |
|
|
} else { |
206 |
|
|
pvl[ 1] = vl; pid[ 1] = id; |
207 |
|
|
pvl[ 2] = pvl[ 3]; pid[ 2] = pid[ 3]; |
208 |
|
|
pvl[ 3] = pvl[ 4]; pid[ 3] = pid[ 4]; |
209 |
|
|
pvl[ 5] = pvl[ 6]; pid[ 5] = pid[ 6]; |
210 |
|
|
pvl[ 8] = pvl[ 9]; pid[ 8] = pid[ 9]; |
211 |
|
|
pvl[11] = pvl[12]; pid[11] = pid[12]; |
212 |
|
|
pvl[ 6] = pvl[ 7]; pid[ 6] = pid[ 7]; |
213 |
|
|
pvl[ 9] = pvl[10]; pid[ 9] = pid[10]; |
214 |
|
|
pvl[12] = pvl[13]; pid[12] = pid[13]; |
215 |
|
|
} |
216 |
|
|
pvl[ 4] = GETV_3(x+1, y-1, z); pid[ 4] = GETI_3(x+1, y-1, z); |
217 |
|
|
pvl[ 7] = GETV_3(x+1, y-1, z-1); pid[ 7] = GETI_3(x+1, y-1, z-1); |
218 |
|
|
pvl[10] = GETV_3(x+1, y, z-1); pid[10] = GETI_3(x+1, y, z-1); |
219 |
|
|
pvl[13] = GETV_3(x+1, y+1, z-1); pid[13] = GETI_3(x+1, y+1, z-1); |
220 |
|
|
vl = GETV_3(x, y, z); |
221 |
|
|
p = 0; |
222 |
|
|
if (vl == pvl[1]) { |
223 |
|
|
id = pid[p=1]; |
224 |
|
|
} |
225 |
|
|
TEST(3); |
226 |
|
|
TEST(9); |
227 |
|
|
if (2 <= conny) { |
228 |
|
|
TEST(2); TEST(4); |
229 |
|
|
TEST(6); TEST(8); TEST(10); TEST(12); |
230 |
|
|
if (3 == conny) { |
231 |
|
|
TEST(5); TEST(7); TEST(11); TEST(13); |
232 |
|
|
} |
233 |
|
|
} |
234 |
|
|
if (!p) { |
235 |
|
|
/* didn't match anything previous */ |
236 |
|
|
id = *numid; |
237 |
|
|
(*numid)++; |
238 |
|
|
} |
239 |
|
|
out[x + sx*(y + sy*z)] = id; |
240 |
|
|
} |
241 |
|
|
} |
242 |
|
|
} |
243 |
|
|
|
244 |
|
|
return 0; |
245 |
|
|
} |
246 |
|
|
|
247 |
|
|
int |
248 |
|
|
_nrrdCCFind_N(Nrrd *nout, unsigned int *numid, airArray *eqvArr, |
249 |
|
|
const Nrrd *nin, unsigned int conny) { |
250 |
|
|
static const char me[]="_nrrdCCFind_N"; |
251 |
|
|
|
252 |
|
|
AIR_UNUSED(nout); |
253 |
|
|
AIR_UNUSED(numid); |
254 |
|
|
AIR_UNUSED(eqvArr); |
255 |
|
|
AIR_UNUSED(nin); |
256 |
|
|
AIR_UNUSED(conny); |
257 |
|
|
biffAddf(NRRD, "%s: sorry, not implemented yet", me); |
258 |
|
|
return 1; |
259 |
|
|
} |
260 |
|
|
|
261 |
|
|
/* |
262 |
|
|
******** nrrdCCFind |
263 |
|
|
** |
264 |
|
|
** finds connected components (CCs) in given integral type nrrd "nin", |
265 |
|
|
** according to connectivity "conny", putting the results in "nout". |
266 |
|
|
** The "type" argument controls what type the output will be. If |
267 |
|
|
** type == nrrdTypeDefault, the type used will be the smallest that |
268 |
|
|
** can contain the CC id values. Otherwise, the specified type "type" |
269 |
|
|
** will be used, assuming that it is large enough to hold the CC ids. |
270 |
|
|
** |
271 |
|
|
** "conny": the number of coordinates that need to varied together in |
272 |
|
|
** order to reach all the samples that are to consitute the neighborhood |
273 |
|
|
** around a sample. For 2-D, conny==1 specifies the 4 edge-connected |
274 |
|
|
** pixels, and 2 specifies the 8 edge- and corner-connected. |
275 |
|
|
** |
276 |
|
|
** The caller can get a record of the values in each CC by passing a |
277 |
|
|
** non-NULL nval, which will be allocated to an array of the same type |
278 |
|
|
** as nin, so that nval->data[I] is the value in nin inside CC #I. |
279 |
|
|
*/ |
280 |
|
|
int |
281 |
|
|
nrrdCCFind(Nrrd *nout, Nrrd **nvalP, const Nrrd *nin, int type, |
282 |
|
|
unsigned int conny) { |
283 |
|
|
static const char me[]="nrrdCCFind", func[]="ccfind"; |
284 |
|
|
Nrrd *nfpid; /* first-pass IDs */ |
285 |
|
|
airArray *mop, *eqvArr; |
286 |
|
|
unsigned int *fpid, numid, numsettleid, *map, |
287 |
|
|
(*lup)(const void *, size_t), (*ins)(void *, size_t, unsigned int); |
288 |
|
|
int ret; |
289 |
|
|
size_t I, NN; |
290 |
|
|
void *val; |
291 |
|
|
|
292 |
|
|
if (!(nout && nin)) { |
293 |
|
|
/* NULL nvalP okay */ |
294 |
|
|
biffAddf(NRRD, "%s: got NULL pointer", me); |
295 |
|
|
return 1; |
296 |
|
|
} |
297 |
|
|
if (nout == nin) { |
298 |
|
|
biffAddf(NRRD, "%s: nout == nin disallowed", me); |
299 |
|
|
return 1; |
300 |
|
|
} |
301 |
|
|
if (!( nrrdTypeIsIntegral[nin->type] |
302 |
|
|
&& nrrdTypeIsUnsigned[nin->type] |
303 |
|
|
&& nrrdTypeSize[nin->type] <= 4 )) { |
304 |
|
|
biffAddf(NRRD, "%s: can only find connected components in " |
305 |
|
|
"1, 2, or 4 byte unsigned integral values (not %s)", |
306 |
|
|
me, airEnumStr(nrrdType, nin->type)); |
307 |
|
|
return 1; |
308 |
|
|
} |
309 |
|
|
if (nrrdTypeDefault != type) { |
310 |
|
|
if (!( AIR_IN_OP(nrrdTypeUnknown, type, nrrdTypeLast) )) { |
311 |
|
|
biffAddf(NRRD, "%s: got invalid target type %d", me, type); |
312 |
|
|
return 1; |
313 |
|
|
} |
314 |
|
|
if (!( nrrdTypeIsIntegral[type] |
315 |
|
|
&& nrrdTypeIsUnsigned[nin->type] |
316 |
|
|
&& nrrdTypeSize[type] <= 4 )) { |
317 |
|
|
biffAddf(NRRD, |
318 |
|
|
"%s: can only save connected components to 1, 2, or 4 byte " |
319 |
|
|
"unsigned integral values (not %s)", |
320 |
|
|
me, airEnumStr(nrrdType, type)); |
321 |
|
|
return 1; |
322 |
|
|
} |
323 |
|
|
} |
324 |
|
|
if (!( conny <= nin->dim )) { |
325 |
|
|
biffAddf(NRRD, "%s: connectivity value must be in [1..%d] for %d-D " |
326 |
|
|
"data (not %d)", me, nin->dim, nin->dim, conny); |
327 |
|
|
return 1; |
328 |
|
|
} |
329 |
|
|
if (nrrdConvert(nfpid=nrrdNew(), nin, nrrdTypeUInt)) { |
330 |
|
|
biffAddf(NRRD, "%s: couldn't allocate fpid %s array to match input size", |
331 |
|
|
me, airEnumStr(nrrdType, nrrdTypeUInt)); |
332 |
|
|
return 1; |
333 |
|
|
} |
334 |
|
|
|
335 |
|
|
mop = airMopNew(); |
336 |
|
|
airMopAdd(mop, nfpid, (airMopper)nrrdNuke, airMopAlways); |
337 |
|
|
eqvArr = airArrayNew(NULL, NULL, 2*sizeof(unsigned int), _nrrdCC_EqvIncr); |
338 |
|
|
airMopAdd(mop, eqvArr, (airMopper)airArrayNuke, airMopAlways); |
339 |
|
|
ret = 0; |
340 |
|
|
switch(nin->dim) { |
341 |
|
|
case 1: |
342 |
|
|
ret = _nrrdCCFind_1(nfpid, &numid, nin); |
343 |
|
|
break; |
344 |
|
|
case 2: |
345 |
|
|
ret = _nrrdCCFind_2(nfpid, &numid, eqvArr, nin, conny); |
346 |
|
|
break; |
347 |
|
|
case 3: |
348 |
|
|
ret = _nrrdCCFind_3(nfpid, &numid, eqvArr, nin, conny); |
349 |
|
|
break; |
350 |
|
|
default: |
351 |
|
|
ret = _nrrdCCFind_N(nfpid, &numid, eqvArr, nin, conny); |
352 |
|
|
break; |
353 |
|
|
} |
354 |
|
|
if (ret) { |
355 |
|
|
biffAddf(NRRD, "%s: initial pass failed", me); |
356 |
|
|
airMopError(mop); return 1; |
357 |
|
|
} |
358 |
|
|
|
359 |
|
|
map = AIR_MALLOC(numid, unsigned int); |
360 |
|
|
airMopAdd(mop, map, airFree, airMopAlways); |
361 |
|
|
numsettleid = airEqvMap(eqvArr, map, numid); |
362 |
|
|
/* convert fpid values to final id values */ |
363 |
|
|
fpid = (unsigned int*)(nfpid->data); |
364 |
|
|
NN = nrrdElementNumber(nfpid); |
365 |
|
|
for (I=0; I<NN; I++) { |
366 |
|
|
fpid[I] = map[fpid[I]]; |
367 |
|
|
} |
368 |
|
|
if (nvalP) { |
369 |
|
|
if (!(*nvalP)) { |
370 |
|
|
*nvalP = nrrdNew(); |
371 |
|
|
} |
372 |
|
|
if (nrrdMaybeAlloc_va(*nvalP, nin->type, 1, |
373 |
|
|
AIR_CAST(size_t, numsettleid))) { |
374 |
|
|
biffAddf(NRRD, "%s: couldn't allocate output value list", me); |
375 |
|
|
airMopError(mop); return 1; |
376 |
|
|
} |
377 |
|
|
airMopAdd(mop, nvalP, (airMopper)airSetNull, airMopOnError); |
378 |
|
|
airMopAdd(mop, *nvalP, (airMopper)nrrdNuke, airMopOnError); |
379 |
|
|
val = (*nvalP)->data; |
380 |
|
|
lup = nrrdUILookup[nin->type]; |
381 |
|
|
ins = nrrdUIInsert[nin->type]; |
382 |
|
|
/* I'm not sure if its more work to do all the redundant assignments |
383 |
|
|
or to check whether or not to do them */ |
384 |
|
|
for (I=0; I<NN; I++) { |
385 |
|
|
ins(val, fpid[I], lup(nin->data, I)); |
386 |
|
|
} |
387 |
|
|
} |
388 |
|
|
|
389 |
|
|
if (nrrdTypeDefault != type) { |
390 |
|
|
if (numsettleid-1 > nrrdTypeMax[type]) { |
391 |
|
|
biffAddf(NRRD, |
392 |
|
|
"%s: max cc id %u is too large to fit in output type %s", |
393 |
|
|
me, numsettleid-1, airEnumStr(nrrdType, type)); |
394 |
|
|
airMopError(mop); return 1; |
395 |
|
|
} |
396 |
|
|
} else { |
397 |
|
|
type = (numsettleid-1 <= nrrdTypeMax[nrrdTypeUChar] |
398 |
|
|
? nrrdTypeUChar |
399 |
|
|
: (numsettleid-1 <= nrrdTypeMax[nrrdTypeUShort] |
400 |
|
|
? nrrdTypeUShort |
401 |
|
|
: nrrdTypeUInt)); |
402 |
|
|
} |
403 |
|
|
if (nrrdConvert(nout, nfpid, type)) { |
404 |
|
|
biffAddf(NRRD, "%s: trouble converting to final output", me); |
405 |
|
|
airMopError(mop); return 1; |
406 |
|
|
} |
407 |
|
|
if (nrrdContentSet_va(nout, func, nin, "%s,%d", |
408 |
|
|
airEnumStr(nrrdType, type), conny)) { |
409 |
|
|
biffAddf(NRRD, "%s:", me); |
410 |
|
|
return 1; |
411 |
|
|
} |
412 |
|
|
if (nout != nin) { |
413 |
|
|
nrrdAxisInfoCopy(nout, nin, NULL, NRRD_AXIS_INFO_NONE); |
414 |
|
|
} |
415 |
|
|
/* basic info handled by nrrdConvert */ |
416 |
|
|
|
417 |
|
|
airMopOkay(mop); |
418 |
|
|
return 0; |
419 |
|
|
} |
420 |
|
|
|
421 |
|
|
int |
422 |
|
|
_nrrdCCAdj_1(unsigned char *out, int numid, const Nrrd *nin) { |
423 |
|
|
|
424 |
|
|
AIR_UNUSED(out); |
425 |
|
|
AIR_UNUSED(numid); |
426 |
|
|
AIR_UNUSED(nin); |
427 |
|
|
return 0; |
428 |
|
|
} |
429 |
|
|
|
430 |
|
|
int |
431 |
|
|
_nrrdCCAdj_2(unsigned char *out, unsigned int numid, const Nrrd *nin, |
432 |
|
|
unsigned int conny) { |
433 |
|
|
unsigned int (*lup)(const void *, size_t), x, y, sx, sy, id=0; |
434 |
|
|
double pid[5]={0,0,0,0,0}; |
435 |
|
|
|
436 |
|
|
lup = nrrdUILookup[nin->type]; |
437 |
|
|
sx = AIR_CAST(unsigned int, nin->axis[0].size); |
438 |
|
|
sy = AIR_CAST(unsigned int, nin->axis[1].size); |
439 |
|
|
for (y=0; y<sy; y++) { |
440 |
|
|
for (x=0; x<sx; x++) { |
441 |
|
|
if (!x) { |
442 |
|
|
pid[1] = GETV_2(-1, y); |
443 |
|
|
pid[2] = GETV_2(-1, y-1); |
444 |
|
|
pid[3] = GETV_2(0, y-1); |
445 |
|
|
} else { |
446 |
|
|
pid[1] = id; |
447 |
|
|
pid[2] = pid[3]; |
448 |
|
|
pid[3] = pid[4]; |
449 |
|
|
} |
450 |
|
|
pid[4] = GETV_2(x+1, y-1); |
451 |
|
|
id = AIR_CAST(unsigned int, GETV_2(x, y)); |
452 |
|
|
#define TADJ(P) \ |
453 |
|
|
if (pid[(P)] != 0.5 && id != pid[(P)]) { \ |
454 |
|
|
out[id + numid*AIR_CAST(unsigned int, pid[(P)])] = \ |
455 |
|
|
out[AIR_CAST(unsigned int, pid[(P)]) + numid*id] = 1; \ |
456 |
|
|
} |
457 |
|
|
TADJ(1); |
458 |
|
|
TADJ(3); |
459 |
|
|
if (2 == conny) { |
460 |
|
|
TADJ(2); |
461 |
|
|
TADJ(4); |
462 |
|
|
} |
463 |
|
|
} |
464 |
|
|
} |
465 |
|
|
|
466 |
|
|
return 0; |
467 |
|
|
} |
468 |
|
|
|
469 |
|
|
int |
470 |
|
|
_nrrdCCAdj_3(unsigned char *out, int numid, const Nrrd *nin, |
471 |
|
|
unsigned int conny) { |
472 |
|
|
unsigned int (*lup)(const void *, size_t), x, y, z, sx, sy, sz, id=0; |
473 |
|
|
double pid[14]={0,0,0,0,0,0,0,0,0,0,0,0,0,0}; |
474 |
|
|
|
475 |
|
|
lup = nrrdUILookup[nin->type]; |
476 |
|
|
sx = AIR_CAST(unsigned int, nin->axis[0].size); |
477 |
|
|
sy = AIR_CAST(unsigned int, nin->axis[1].size); |
478 |
|
|
sz = AIR_CAST(unsigned int, nin->axis[2].size); |
479 |
|
|
for (z=0; z<sz; z++) { |
480 |
|
|
for (y=0; y<sy; y++) { |
481 |
|
|
for (x=0; x<sx; x++) { |
482 |
|
|
if (!x) { |
483 |
|
|
pid[ 1] = GETV_3(-1, y, z); |
484 |
|
|
pid[ 2] = GETV_3(-1, y-1, z); |
485 |
|
|
pid[ 3] = GETV_3( 0, y-1, z); |
486 |
|
|
pid[ 5] = GETV_3(-1, y-1, z-1); |
487 |
|
|
pid[ 8] = GETV_3(-1, y, z-1); |
488 |
|
|
pid[11] = GETV_3(-1, y+1, z-1); |
489 |
|
|
pid[ 6] = GETV_3( 0, y-1, z-1); |
490 |
|
|
pid[ 9] = GETV_3( 0, y, z-1); |
491 |
|
|
pid[12] = GETV_3( 0, y+1, z-1); |
492 |
|
|
} else { |
493 |
|
|
pid[ 1] = id; |
494 |
|
|
pid[ 2] = pid[ 3]; |
495 |
|
|
pid[ 3] = pid[ 4]; |
496 |
|
|
pid[ 5] = pid[ 6]; |
497 |
|
|
pid[ 8] = pid[ 9]; |
498 |
|
|
pid[11] = pid[12]; |
499 |
|
|
pid[ 6] = pid[ 7]; |
500 |
|
|
pid[ 9] = pid[10]; |
501 |
|
|
pid[12] = pid[13]; |
502 |
|
|
} |
503 |
|
|
pid[ 4] = GETV_3(x+1, y-1, z); |
504 |
|
|
pid[ 7] = GETV_3(x+1, y-1, z-1); |
505 |
|
|
pid[10] = GETV_3(x+1, y, z-1); |
506 |
|
|
pid[13] = GETV_3(x+1, y+1, z-1); |
507 |
|
|
id = AIR_CAST(unsigned int, GETV_3(x, y, z)); |
508 |
|
|
TADJ(1); |
509 |
|
|
TADJ(3); |
510 |
|
|
TADJ(9); |
511 |
|
|
if (2 <= conny) { |
512 |
|
|
TADJ(2); TADJ(4); |
513 |
|
|
TADJ(6); TADJ(8); TADJ(10); TADJ(12); |
514 |
|
|
if (3 == conny) { |
515 |
|
|
TADJ(5); TADJ(7); TADJ(11); TADJ(13); |
516 |
|
|
} |
517 |
|
|
} |
518 |
|
|
} |
519 |
|
|
} |
520 |
|
|
} |
521 |
|
|
|
522 |
|
|
return 0; |
523 |
|
|
} |
524 |
|
|
|
525 |
|
|
int |
526 |
|
|
_nrrdCCAdj_N(unsigned char *out, int numid, const Nrrd *nin, |
527 |
|
|
unsigned int conny) { |
528 |
|
|
static const char me[]="_nrrdCCAdj_N"; |
529 |
|
|
|
530 |
|
|
AIR_UNUSED(out); |
531 |
|
|
AIR_UNUSED(numid); |
532 |
|
|
AIR_UNUSED(nin); |
533 |
|
|
AIR_UNUSED(conny); |
534 |
|
|
biffAddf(NRRD, "%s: sorry, not implemented", me); |
535 |
|
|
return 1; |
536 |
|
|
} |
537 |
|
|
|
538 |
|
|
int |
539 |
|
|
nrrdCCAdjacency(Nrrd *nout, const Nrrd *nin, unsigned int conny) { |
540 |
|
|
static const char me[]="nrrdCCAdjacency", func[]="ccadj"; |
541 |
|
|
int ret; |
542 |
|
|
unsigned int maxid; |
543 |
|
|
unsigned char *out; |
544 |
|
|
|
545 |
|
|
if (!( nout && nrrdCCValid(nin) )) { |
546 |
|
|
biffAddf(NRRD, "%s: invalid args", me); |
547 |
|
|
return 1; |
548 |
|
|
} |
549 |
|
|
if (nout == nin) { |
550 |
|
|
biffAddf(NRRD, "%s: nout == nin disallowed", me); |
551 |
|
|
return 1; |
552 |
|
|
} |
553 |
|
|
if (!( AIR_IN_CL(1, conny, nin->dim) )) { |
554 |
|
|
biffAddf(NRRD, "%s: connectivity value must be in [1..%d] for %d-D " |
555 |
|
|
"data (not %d)", me, nin->dim, nin->dim, conny); |
556 |
|
|
return 1; |
557 |
|
|
} |
558 |
|
|
maxid = nrrdCCMax(nin); |
559 |
|
|
if (nrrdMaybeAlloc_va(nout, nrrdTypeUChar, 2, |
560 |
|
|
AIR_CAST(size_t, maxid+1), |
561 |
|
|
AIR_CAST(size_t, maxid+1))) { |
562 |
|
|
biffAddf(NRRD, "%s: trouble allocating output", me); |
563 |
|
|
return 1; |
564 |
|
|
} |
565 |
|
|
out = (unsigned char *)(nout->data); |
566 |
|
|
|
567 |
|
|
switch(nin->dim) { |
568 |
|
|
case 1: |
569 |
|
|
ret = _nrrdCCAdj_1(out, maxid+1, nin); |
570 |
|
|
break; |
571 |
|
|
case 2: |
572 |
|
|
ret = _nrrdCCAdj_2(out, maxid+1, nin, conny); |
573 |
|
|
break; |
574 |
|
|
case 3: |
575 |
|
|
ret = _nrrdCCAdj_3(out, maxid+1, nin, conny); |
576 |
|
|
break; |
577 |
|
|
default: |
578 |
|
|
ret = _nrrdCCAdj_N(out, maxid+1, nin, conny); |
579 |
|
|
break; |
580 |
|
|
} |
581 |
|
|
if (ret) { |
582 |
|
|
biffAddf(NRRD, "%s: trouble", me); |
583 |
|
|
return 1; |
584 |
|
|
} |
585 |
|
|
/* this goofiness is just so that histo-based projections |
586 |
|
|
return the sorts of values that we expect */ |
587 |
|
|
nout->axis[0].center = nout->axis[1].center = nrrdCenterCell; |
588 |
|
|
nout->axis[0].min = nout->axis[1].min = -0.5; |
589 |
|
|
nout->axis[0].max = nout->axis[1].max = maxid + 0.5; |
590 |
|
|
if (nrrdContentSet_va(nout, func, nin, "%d", conny)) { |
591 |
|
|
biffAddf(NRRD, "%s:", me); |
592 |
|
|
return 1; |
593 |
|
|
} |
594 |
|
|
|
595 |
|
|
return 0; |
596 |
|
|
} |
597 |
|
|
|
598 |
|
|
/* |
599 |
|
|
******** nrrdCCMerge |
600 |
|
|
** |
601 |
|
|
** Slightly-too-multi-purpose tool for merging small connected components |
602 |
|
|
** (CCs) into larger ones, according to a number of possible different |
603 |
|
|
** constraints, as explained below. |
604 |
|
|
** |
605 |
|
|
** valDir: (value direction) uses information about the original values |
606 |
|
|
** in the CC to constrain whether darker gets merged into brighter, or vice |
607 |
|
|
** versa, or neither. For non-zero valDir values, a non-NULL _nval (from |
608 |
|
|
** nrrdCCFind) must be passed. |
609 |
|
|
** valDir > 0 : merge dark CCs into bright, but not vice versa |
610 |
|
|
** valDir = 0 : merge either way, values are irrelevant |
611 |
|
|
** valDir < 0 : merge bright CCs into dark, but not vice versa |
612 |
|
|
** When merging with multiple neighbors (maxNeighbor > 1), the value |
613 |
|
|
** of the largest neighbor is considered. |
614 |
|
|
** |
615 |
|
|
** maxSize: a cap on how large "small" is- CCs any larger than maxSize are |
616 |
|
|
** not merged, as they are deemed too significant. Or, a maxSize of 0 says |
617 |
|
|
** size is no object for merging CCs. |
618 |
|
|
** |
619 |
|
|
** maxNeighbor: a maximum number of neighbors that a CC can have (either |
620 |
|
|
** bigger than the CC or not) if it is to be merged. Use 1 to merge |
621 |
|
|
** isolated islands into their surrounds, 2 to merge CC with the larger |
622 |
|
|
** of their two neighbors, etc., or 0 to allow any number of neighbors. |
623 |
|
|
** |
624 |
|
|
** conny: passed to nrrdCCAdjacency() when determining neighbors |
625 |
|
|
** |
626 |
|
|
** In order to prevent weirdness, the merging done in one call to this |
627 |
|
|
** function is not transitive: if A is merged to B, then B will not be |
628 |
|
|
** merged to anything else, even if meets all the requirements defined |
629 |
|
|
** by the given parameters. This is accomplished by working from the |
630 |
|
|
** smallest CCs to the largest. Iterated calls may be needed to acheive |
631 |
|
|
** the desired effect. |
632 |
|
|
** |
633 |
|
|
** Note: the output of this is not "settled"- the CC id values are not |
634 |
|
|
** shiftward downwards to their lowest possible values, since this would |
635 |
|
|
** needlessly invalidate the nval value store. |
636 |
|
|
*/ |
637 |
|
|
int |
638 |
|
|
nrrdCCMerge(Nrrd *nout, const Nrrd *nin, Nrrd *_nval, |
639 |
|
|
int valDir, unsigned int maxSize, unsigned int maxNeighbor, |
640 |
|
|
unsigned int conny) { |
641 |
|
|
static const char me[]="nrrdCCMerge", func[]="ccmerge"; |
642 |
|
|
const char *valcnt; |
643 |
|
|
unsigned int _i, i, j, bigi=0, numid, *size, *sizeId, |
644 |
|
|
*nn, /* number of neighbors */ |
645 |
|
|
*val=NULL, *hit, |
646 |
|
|
(*lup)(const void *, size_t), (*ins)(void *, size_t, unsigned int); |
647 |
|
|
Nrrd *nadj, *nsize, *nval=NULL, *nnn; |
648 |
|
|
unsigned char *adj; |
649 |
|
|
unsigned int *map, *id; |
650 |
|
|
airArray *mop; |
651 |
|
|
size_t I, NN; |
652 |
|
|
|
653 |
|
|
mop = airMopNew(); |
654 |
|
|
if (!( nout && nrrdCCValid(nin) )) { |
655 |
|
|
/* _nval can be NULL */ |
656 |
|
|
biffAddf(NRRD, "%s: invalid args", me); |
657 |
|
|
airMopError(mop); return 1; |
658 |
|
|
} |
659 |
|
|
if (valDir) { |
660 |
|
|
airMopAdd(mop, nval = nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
661 |
|
|
if (nrrdConvert(nval, _nval, nrrdTypeUInt)) { |
662 |
|
|
biffAddf(NRRD, "%s: value-directed merging needs usable nval", me); |
663 |
|
|
airMopError(mop); return 1; |
664 |
|
|
} |
665 |
|
|
val = (unsigned int*)(nval->data); |
666 |
|
|
} |
667 |
|
|
if (nout != nin) { |
668 |
|
|
if (nrrdCopy(nout, nin)) { |
669 |
|
|
biffAddf(NRRD, "%s:", me); |
670 |
|
|
airMopError(mop); return 1; |
671 |
|
|
} |
672 |
|
|
} |
673 |
|
|
airMopAdd(mop, nadj = nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
674 |
|
|
airMopAdd(mop, nsize = nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
675 |
|
|
airMopAdd(mop, nnn = nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
676 |
|
|
|
677 |
|
|
if (nrrdCCSize(nsize, nin) |
678 |
|
|
|| nrrdCopy(nnn, nsize) /* just to allocate to right size and type */ |
679 |
|
|
|| nrrdCCAdjacency(nadj, nin, conny)) { |
680 |
|
|
biffAddf(NRRD, "%s:", me); |
681 |
|
|
airMopError(mop); return 1; |
682 |
|
|
} |
683 |
|
|
size = (unsigned int*)(nsize->data); |
684 |
|
|
adj = (unsigned char*)(nadj->data); |
685 |
|
|
nn = (unsigned int*)(nnn->data); |
686 |
|
|
numid = AIR_CAST(unsigned int, nsize->axis[0].size); |
687 |
|
|
for (i=0; i<numid; i++) { |
688 |
|
|
nn[i] = 0; |
689 |
|
|
for (j=0; j<numid; j++) { |
690 |
|
|
nn[i] += adj[j + numid*i]; |
691 |
|
|
} |
692 |
|
|
} |
693 |
|
|
map = AIR_MALLOC(numid, unsigned int); |
694 |
|
|
id = AIR_MALLOC(numid, unsigned int); |
695 |
|
|
hit = AIR_MALLOC(numid, unsigned int); |
696 |
|
|
sizeId = AIR_MALLOC(2*numid, unsigned int); |
697 |
|
|
/* we add to the mops BEFORE error checking so that anything non-NULL |
698 |
|
|
will get airFree'd, and happily airFree is a no-op on NULL */ |
699 |
|
|
airMopAdd(mop, map, airFree, airMopAlways); |
700 |
|
|
airMopAdd(mop, id, airFree, airMopAlways); |
701 |
|
|
airMopAdd(mop, hit, airFree, airMopAlways); |
702 |
|
|
airMopAdd(mop, sizeId, airFree, airMopAlways); |
703 |
|
|
if (!(map && id && hit && sizeId)) { |
704 |
|
|
biffAddf(NRRD, "%s: couldn't allocate buffers", me); |
705 |
|
|
airMopError(mop); return 1; |
706 |
|
|
} |
707 |
|
|
|
708 |
|
|
/* store and sort size/id pairs */ |
709 |
|
|
for (i=0; i<numid; i++) { |
710 |
|
|
sizeId[0 + 2*i] = size[i]; |
711 |
|
|
sizeId[1 + 2*i] = i; |
712 |
|
|
} |
713 |
|
|
qsort(sizeId, numid, 2*sizeof(unsigned int), nrrdValCompare[nrrdTypeUInt]); |
714 |
|
|
for (i=0; i<numid; i++) { |
715 |
|
|
id[i] = sizeId[1 + 2*i]; |
716 |
|
|
} |
717 |
|
|
|
718 |
|
|
/* initialize arrays */ |
719 |
|
|
for (i=0; i<numid; i++) { |
720 |
|
|
map[i] = i; |
721 |
|
|
hit[i] = AIR_FALSE; |
722 |
|
|
} |
723 |
|
|
/* _i goes through 0 to numid-1, |
724 |
|
|
i goes through the CC ids in ascending order of size */ |
725 |
|
|
for (_i=0; _i<numid; _i++) { |
726 |
|
|
i = id[_i]; |
727 |
|
|
if (hit[i]) { |
728 |
|
|
continue; |
729 |
|
|
} |
730 |
|
|
if (maxSize && (size[i] > maxSize)) { |
731 |
|
|
continue; |
732 |
|
|
} |
733 |
|
|
if (maxNeighbor && (nn[i] > maxNeighbor)) { |
734 |
|
|
continue; |
735 |
|
|
} |
736 |
|
|
/* find biggest neighbor, exploiting the fact that we already |
737 |
|
|
sorted CC ids on size. j descends through indices of id[], |
738 |
|
|
bigi goes through CC ids which are larger than CC i */ |
739 |
|
|
for (j=numid-1; j>_i; j--) { |
740 |
|
|
bigi = id[j]; |
741 |
|
|
if (adj[bigi + numid*i]) |
742 |
|
|
break; |
743 |
|
|
} |
744 |
|
|
if (j == _i) { |
745 |
|
|
continue; /* we had no neighbors ?!?! */ |
746 |
|
|
} |
747 |
|
|
if (valDir && (AIR_CAST(int, val[bigi]) |
748 |
|
|
- AIR_CAST(int, val[i]))*valDir < 0 ) { |
749 |
|
|
continue; |
750 |
|
|
} |
751 |
|
|
/* else all criteria for merging have been met */ |
752 |
|
|
map[i] = bigi; |
753 |
|
|
hit[bigi] = AIR_TRUE; |
754 |
|
|
} |
755 |
|
|
lup = nrrdUILookup[nin->type]; |
756 |
|
|
ins = nrrdUIInsert[nout->type]; |
757 |
|
|
NN = nrrdElementNumber(nin); |
758 |
|
|
for (I=0; I<NN; I++) { |
759 |
|
|
ins(nout->data, I, map[lup(nin->data, I)]); |
760 |
|
|
} |
761 |
|
|
|
762 |
|
|
valcnt = ((_nval && _nval->content) |
763 |
|
|
? _nval->content |
764 |
|
|
: nrrdStateUnknownContent); |
765 |
|
|
if ( (valDir && nrrdContentSet_va(nout, func, nin, "%c(%s),%d,%d,%d", |
766 |
|
|
(valDir > 0 ? '+' : '-'), valcnt, |
767 |
|
|
maxSize, maxNeighbor, conny)) |
768 |
|
|
|| |
769 |
|
|
(!valDir && nrrdContentSet_va(nout, func, nin, ".,%d,%d,%d", |
770 |
|
|
maxSize, maxNeighbor, conny)) ) { |
771 |
|
|
biffAddf(NRRD, "%s:", me); |
772 |
|
|
airMopError(mop); return 1; |
773 |
|
|
} |
774 |
|
|
/* basic info handled by nrrdCopy */ |
775 |
|
|
airMopOkay(mop); |
776 |
|
|
return 0; |
777 |
|
|
} |
778 |
|
|
|
779 |
|
|
/* |
780 |
|
|
******** nrrdCCRevalue() |
781 |
|
|
** |
782 |
|
|
** assigns the original values back to the connected components |
783 |
|
|
** obviously, this could be subsumed by nrrdApply1DLut(), but this |
784 |
|
|
** is so special purpose that it seemed simpler to code from scratch |
785 |
|
|
*/ |
786 |
|
|
int |
787 |
|
|
nrrdCCRevalue (Nrrd *nout, const Nrrd *nin, const Nrrd *nval) { |
788 |
|
|
static const char me[]="nrrdCCRevalue"; |
789 |
|
|
size_t I, NN; |
790 |
|
|
unsigned int (*vlup)(const void *, size_t), (*ilup)(const void *, size_t), |
791 |
|
|
(*ins)(void *, size_t, unsigned int); |
792 |
|
|
|
793 |
|
|
if (!( nout && nrrdCCValid(nin) && nval )) { |
794 |
|
|
biffAddf(NRRD, "%s: invalid args", me); |
795 |
|
|
return 1; |
796 |
|
|
} |
797 |
|
|
if (nrrdConvert(nout, nin, nval->type)) { |
798 |
|
|
biffAddf(NRRD, "%s: couldn't initialize output", me); |
799 |
|
|
return 1; |
800 |
|
|
} |
801 |
|
|
NN = nrrdElementNumber(nin); |
802 |
|
|
vlup = nrrdUILookup[nval->type]; |
803 |
|
|
ilup = nrrdUILookup[nin->type]; |
804 |
|
|
ins = nrrdUIInsert[nout->type]; |
805 |
|
|
for (I=0; I<NN; I++) { |
806 |
|
|
ins(nout->data, I, vlup(nval->data, ilup(nin->data, I))); |
807 |
|
|
} |
808 |
|
|
/* basic info handled by nrrdConvert */ |
809 |
|
|
|
810 |
|
|
return 0; |
811 |
|
|
} |
812 |
|
|
|
813 |
|
|
int |
814 |
|
|
nrrdCCSettle(Nrrd *nout, Nrrd **nvalP, const Nrrd *nin) { |
815 |
|
|
static const char me[]="nrrdCCSettle", func[]="ccsettle"; |
816 |
|
|
unsigned int numid, maxid, jd, id, *map, |
817 |
|
|
(*lup)(const void *, size_t), (*ins)(void *, size_t, unsigned int); |
818 |
|
|
size_t I, NN; |
819 |
|
|
airArray *mop; |
820 |
|
|
|
821 |
|
|
mop = airMopNew(); |
822 |
|
|
if (!( nout && nrrdCCValid(nin) )) { |
823 |
|
|
/* nvalP can be NULL */ |
824 |
|
|
biffAddf(NRRD, "%s: invalid args", me); |
825 |
|
|
airMopError(mop); return 1; |
826 |
|
|
} |
827 |
|
|
if (nrrdCopy(nout, nin)) { |
828 |
|
|
biffAddf(NRRD, "%s: initial copy failed", me); |
829 |
|
|
airMopError(mop); return 1; |
830 |
|
|
} |
831 |
|
|
maxid = nrrdCCMax(nin); |
832 |
|
|
lup = nrrdUILookup[nin->type]; |
833 |
|
|
ins = nrrdUIInsert[nin->type]; |
834 |
|
|
NN = nrrdElementNumber(nin); |
835 |
|
|
map = AIR_CALLOC(maxid+1, unsigned int); /* we do need it zeroed out */ |
836 |
|
|
if (!map) { |
837 |
|
|
biffAddf(NRRD, "%s: couldn't allocate internal LUT", me); |
838 |
|
|
airMopError(mop); return 1; |
839 |
|
|
} |
840 |
|
|
airMopAdd(mop, map, airFree, airMopAlways); |
841 |
|
|
for (I=0; I<NN; I++) { |
842 |
|
|
map[lup(nin->data, I)] = 1; |
843 |
|
|
} |
844 |
|
|
numid = 0; |
845 |
|
|
for (jd=0; jd<=maxid; jd++) { |
846 |
|
|
numid += map[jd]; |
847 |
|
|
} |
848 |
|
|
|
849 |
|
|
if (nvalP) { |
850 |
|
|
if (!(*nvalP)) { |
851 |
|
|
*nvalP = nrrdNew(); |
852 |
|
|
} |
853 |
|
|
if (nrrdMaybeAlloc_va(*nvalP, nin->type, 1, |
854 |
|
|
AIR_CAST(size_t, numid))) { |
855 |
|
|
biffAddf(NRRD, "%s: couldn't allocate output value list", me); |
856 |
|
|
airMopError(mop); return 1; |
857 |
|
|
} |
858 |
|
|
airMopAdd(mop, nvalP, (airMopper)airSetNull, airMopOnError); |
859 |
|
|
airMopAdd(mop, *nvalP, (airMopper)nrrdNuke, airMopOnError); |
860 |
|
|
} |
861 |
|
|
|
862 |
|
|
id = 0; |
863 |
|
|
for (jd=0; jd<=maxid; jd++) { |
864 |
|
|
if (map[jd]) { |
865 |
|
|
map[jd] = id; |
866 |
|
|
if (nvalP) { |
867 |
|
|
ins((*nvalP)->data, id, jd); |
868 |
|
|
} |
869 |
|
|
id++; |
870 |
|
|
} |
871 |
|
|
} |
872 |
|
|
for (I=0; I<NN; I++) { |
873 |
|
|
ins(nout->data, I, map[lup(nin->data, I)]); |
874 |
|
|
} |
875 |
|
|
|
876 |
|
|
if (nrrdContentSet_va(nout, func, nin, "")) { |
877 |
|
|
biffAddf(NRRD, "%s:", me); |
878 |
|
|
airMopError(mop); return 1; |
879 |
|
|
} |
880 |
|
|
/* basic info handled by nrrdCopy */ |
881 |
|
|
airMopOkay(mop); |
882 |
|
|
return 0; |
883 |
|
|
} |