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 |
|
|
|
25 |
|
|
#include "ten.h" |
26 |
|
|
#include "privateTen.h" |
27 |
|
|
|
28 |
|
|
int |
29 |
|
|
tenEvecRGB(Nrrd *nout, const Nrrd *nin, |
30 |
|
|
const tenEvecRGBParm *rgbp) { |
31 |
|
|
static const char me[]="tenEvecRGB"; |
32 |
|
|
size_t size[NRRD_DIM_MAX]; |
33 |
|
|
float (*lup)(const void *, size_t), (*ins)(void *, size_t, float); |
34 |
|
|
float ten[7], eval[3], evec[9], RGB[3]; |
35 |
|
|
size_t II, NN; |
36 |
|
|
unsigned char *odataUC; |
37 |
|
|
unsigned short *odataUS; |
38 |
|
|
|
39 |
|
|
if (!(nout && nin)) { |
40 |
|
|
biffAddf(TEN, "%s: got NULL pointer (%p,%p)", |
41 |
|
|
me, AIR_CAST(void *, nout), AIR_CVOIDP(nin)); |
42 |
|
|
return 1; |
43 |
|
|
} |
44 |
|
|
if (tenEvecRGBParmCheck(rgbp)) { |
45 |
|
|
biffAddf(TEN, "%s: RGB parm trouble", me); |
46 |
|
|
return 1; |
47 |
|
|
} |
48 |
|
|
if (!(2 <= nin->dim && 7 == nin->axis[0].size)) { |
49 |
|
|
char stmp[AIR_STRLEN_SMALL]; |
50 |
|
|
biffAddf(TEN, "%s: need nin->dim >= 2 (not %u), axis[0].size == 7 " |
51 |
|
|
"(not %s)", me, nin->dim, |
52 |
|
|
airSprintSize_t(stmp, nin->axis[0].size)); |
53 |
|
|
return 1; |
54 |
|
|
} |
55 |
|
|
|
56 |
|
|
nrrdAxisInfoGet_nva(nin, nrrdAxisInfoSize, size); |
57 |
|
|
size[0] = rgbp->genAlpha ? 4 : 3; |
58 |
|
|
if (nrrdMaybeAlloc_nva(nout, (nrrdTypeDefault == rgbp->typeOut |
59 |
|
|
? nin->type |
60 |
|
|
: rgbp->typeOut), nin->dim, size)) { |
61 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't alloc output", me); |
62 |
|
|
return 1; |
63 |
|
|
} |
64 |
|
|
odataUC = AIR_CAST(unsigned char *, nout->data); |
65 |
|
|
odataUS = AIR_CAST(unsigned short *, nout->data); |
66 |
|
|
|
67 |
|
|
NN = nrrdElementNumber(nin)/7; |
68 |
|
|
lup = nrrdFLookup[nin->type]; |
69 |
|
|
ins = nrrdFInsert[nout->type]; |
70 |
|
|
for (II=0; II<NN; II++) { |
71 |
|
|
TEN_T_SET(ten, lup(nin->data, 0 + 7*II), |
72 |
|
|
lup(nin->data, 1 + 7*II), lup(nin->data, 2 + 7*II), |
73 |
|
|
lup(nin->data, 3 + 7*II), lup(nin->data, 4 + 7*II), |
74 |
|
|
lup(nin->data, 5 + 7*II), lup(nin->data, 6 + 7*II)); |
75 |
|
|
tenEigensolve_f(eval, evec, ten); |
76 |
|
|
tenEvecRGBSingle_f(RGB, ten[0], eval, evec + 3*(rgbp->which), rgbp); |
77 |
|
|
switch (nout->type) { |
78 |
|
|
case nrrdTypeUChar: |
79 |
|
|
odataUC[0 + size[0]*II] = airIndexClamp(0.0, RGB[0], 1.0, 256); |
80 |
|
|
odataUC[1 + size[0]*II] = airIndexClamp(0.0, RGB[1], 1.0, 256); |
81 |
|
|
odataUC[2 + size[0]*II] = airIndexClamp(0.0, RGB[2], 1.0, 256); |
82 |
|
|
if (rgbp->genAlpha) { |
83 |
|
|
odataUC[3 + size[0]*II] = 255; |
84 |
|
|
} |
85 |
|
|
break; |
86 |
|
|
case nrrdTypeUShort: |
87 |
|
|
odataUS[0 + size[0]*II] = airIndexClamp(0.0, RGB[0], 1.0, 65536); |
88 |
|
|
odataUS[1 + size[0]*II] = airIndexClamp(0.0, RGB[1], 1.0, 65536); |
89 |
|
|
odataUS[2 + size[0]*II] = airIndexClamp(0.0, RGB[2], 1.0, 65536); |
90 |
|
|
if (rgbp->genAlpha) { |
91 |
|
|
odataUS[3 + size[0]*II] = 65535; |
92 |
|
|
} |
93 |
|
|
break; |
94 |
|
|
default: |
95 |
|
|
ins(nout->data, 0 + size[0]*II, RGB[0]); |
96 |
|
|
ins(nout->data, 1 + size[0]*II, RGB[1]); |
97 |
|
|
ins(nout->data, 2 + size[0]*II, RGB[2]); |
98 |
|
|
if (rgbp->genAlpha) { |
99 |
|
|
ins(nout->data, 3 + size[0]*II, 1.0); |
100 |
|
|
} |
101 |
|
|
break; |
102 |
|
|
} |
103 |
|
|
} |
104 |
|
|
if (nrrdAxisInfoCopy(nout, nin, NULL, (NRRD_AXIS_INFO_SIZE_BIT))) { |
105 |
|
|
biffMovef(TEN, NRRD, "%s: couldn't copy axis info", me); |
106 |
|
|
return 1; |
107 |
|
|
} |
108 |
|
|
nout->axis[0].kind = nrrdKind3Color; |
109 |
|
|
if (nrrdBasicInfoCopy(nout, nin, |
110 |
|
|
NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) { |
111 |
|
|
biffAddf(TEN, "%s:", me); |
112 |
|
|
return 1; |
113 |
|
|
} |
114 |
|
|
|
115 |
|
|
return 0; |
116 |
|
|
} |
117 |
|
|
|
118 |
|
|
#define SQR(i) ((i)*(i)) |
119 |
|
|
|
120 |
|
|
short |
121 |
|
|
tenEvqSingle(float vec[3], float scl) { |
122 |
|
|
static const char me[]="tenEvqSingle"; |
123 |
|
|
float tmp, L1; |
124 |
|
|
int mi, bins, base, vi, ui; |
125 |
|
|
short ret; |
126 |
|
|
|
127 |
|
|
ELL_3V_NORM_TT(vec, float, vec, tmp); |
128 |
|
|
L1 = AIR_ABS(vec[0]) + AIR_ABS(vec[1]) + AIR_ABS(vec[2]); |
129 |
|
|
ELL_3V_SCALE(vec, 1/L1, vec); |
130 |
|
|
scl = AIR_CLAMP(0.0f, scl, 1.0f); |
131 |
|
|
scl = AIR_CAST(float, pow(scl, 0.75)); |
132 |
|
|
mi = airIndex(0.0, scl, 1.0, 6); |
133 |
|
|
if (mi) { |
134 |
|
|
switch (mi) { |
135 |
|
|
case 1: bins = 16; base = 1; break; |
136 |
|
|
case 2: bins = 32; base = 1+SQR(16); break; |
137 |
|
|
case 3: bins = 48; base = 1+SQR(16)+SQR(32); break; |
138 |
|
|
case 4: bins = 64; base = 1+SQR(16)+SQR(32)+SQR(48); break; |
139 |
|
|
case 5: bins = 80; base = 1+SQR(16)+SQR(32)+SQR(48)+SQR(64); break; |
140 |
|
|
default: |
141 |
|
|
fprintf(stderr, "%s: PANIC: mi = %d\n", me, mi); |
142 |
|
|
exit(0); |
143 |
|
|
} |
144 |
|
|
vi = airIndex(-1, vec[0]+vec[1], 1, bins); |
145 |
|
|
ui = airIndex(-1, vec[0]-vec[1], 1, bins); |
146 |
|
|
ret = vi*bins + ui + base; |
147 |
|
|
} |
148 |
|
|
else { |
149 |
|
|
ret = 0; |
150 |
|
|
} |
151 |
|
|
return ret; |
152 |
|
|
} |
153 |
|
|
|
154 |
|
|
int |
155 |
|
|
tenEvqVolume(Nrrd *nout, |
156 |
|
|
const Nrrd *nin, int which, int aniso, int scaleByAniso) { |
157 |
|
|
static const char me[]="tenEvqVolume"; |
158 |
|
|
int map[3]; |
159 |
|
|
short *qdata; |
160 |
|
|
const float *tdata; |
161 |
|
|
float eval[3], evec[9], an; |
162 |
|
|
size_t N, I, sx, sy, sz; |
163 |
|
|
|
164 |
|
|
if (!(nout && nin)) { |
165 |
|
|
biffAddf(TEN, "%s: got NULL pointer", me); |
166 |
|
|
return 1; |
167 |
|
|
} |
168 |
|
|
if (!(AIR_IN_CL(0, which, 2))) { |
169 |
|
|
biffAddf(TEN, "%s: eigenvector index %d not in range [0..2]", me, which); |
170 |
|
|
return 1; |
171 |
|
|
} |
172 |
|
|
if (scaleByAniso) { |
173 |
|
|
if (airEnumValCheck(tenAniso, aniso)) { |
174 |
|
|
biffAddf(TEN, "%s: anisotropy metric %d not valid", me, aniso); |
175 |
|
|
return 1; |
176 |
|
|
} |
177 |
|
|
} |
178 |
|
|
if (tenTensorCheck(nin, nrrdTypeFloat, AIR_TRUE, AIR_TRUE)) { |
179 |
|
|
biffAddf(TEN, "%s: didn't get a valid DT volume", me); |
180 |
|
|
return 1; |
181 |
|
|
} |
182 |
|
|
sx = nin->axis[1].size; |
183 |
|
|
sy = nin->axis[2].size; |
184 |
|
|
sz = nin->axis[3].size; |
185 |
|
|
if (nrrdMaybeAlloc_va(nout, nrrdTypeShort, 3, |
186 |
|
|
sx, sy, sz)) { |
187 |
|
|
biffMovef(TEN, NRRD, "%s: can't allocate output", me); |
188 |
|
|
return 1; |
189 |
|
|
} |
190 |
|
|
N = sx*sy*sz; |
191 |
|
|
tdata = (float *)nin->data; |
192 |
|
|
qdata = (short *)nout->data; |
193 |
|
|
for (I=0; I<N; I++) { |
194 |
|
|
tenEigensolve_f(eval, evec, tdata); |
195 |
|
|
if (scaleByAniso) { |
196 |
|
|
an = tenAnisoEval_f(eval, aniso); |
197 |
|
|
} else { |
198 |
|
|
an = 1.0; |
199 |
|
|
} |
200 |
|
|
qdata[I] = tenEvqSingle(evec+ 3*which, an); |
201 |
|
|
tdata += 7; |
202 |
|
|
} |
203 |
|
|
ELL_3V_SET(map, 1, 2, 3); |
204 |
|
|
if (nrrdAxisInfoCopy(nout, nin, map, (NRRD_AXIS_INFO_SIZE_BIT |
205 |
|
|
| NRRD_AXIS_INFO_KIND_BIT) )) { |
206 |
|
|
biffMovef(TEN, NRRD, "%s: trouble", me); |
207 |
|
|
return 1; |
208 |
|
|
} |
209 |
|
|
if (nrrdBasicInfoCopy(nout, nin, |
210 |
|
|
NRRD_BASIC_INFO_ALL ^ NRRD_BASIC_INFO_SPACE)) { |
211 |
|
|
biffAddf(TEN, "%s:", me); |
212 |
|
|
return 1; |
213 |
|
|
} |
214 |
|
|
|
215 |
|
|
return 0; |
216 |
|
|
} |
217 |
|
|
|
218 |
|
|
int |
219 |
|
|
tenBMatrixCheck(const Nrrd *nbmat, int type, unsigned int minnum) { |
220 |
|
|
static const char me[]="tenBMatrixCheck"; |
221 |
|
|
|
222 |
|
|
if (nrrdCheck(nbmat)) { |
223 |
|
|
biffMovef(TEN, NRRD, "%s: basic validity check failed", me); |
224 |
|
|
return 1; |
225 |
|
|
} |
226 |
|
|
if (!( 6 == nbmat->axis[0].size && 2 == nbmat->dim )) { |
227 |
|
|
char stmp[AIR_STRLEN_SMALL]; |
228 |
|
|
biffAddf(TEN, "%s: need a 6xN 2-D array (not a %s x? %d-D array)", me, |
229 |
|
|
airSprintSize_t(stmp, nbmat->axis[0].size), nbmat->dim); |
230 |
|
|
return 1; |
231 |
|
|
} |
232 |
|
|
if (nrrdTypeDefault != type && type != nbmat->type) { |
233 |
|
|
biffAddf(TEN, "%s: requested type %s but got type %s", me, |
234 |
|
|
airEnumStr(nrrdType, type), airEnumStr(nrrdType, nbmat->type)); |
235 |
|
|
return 1; |
236 |
|
|
} |
237 |
|
|
if (nrrdTypeBlock == nbmat->type) { |
238 |
|
|
biffAddf(TEN, "%s: sorry, can't use %s type", me, |
239 |
|
|
airEnumStr(nrrdType, nrrdTypeBlock)); |
240 |
|
|
return 1; |
241 |
|
|
} |
242 |
|
|
if (!( minnum <= nbmat->axis[1].size )) { |
243 |
|
|
char stmp[AIR_STRLEN_SMALL]; |
244 |
|
|
biffAddf(TEN, "%s: have only %s B-matrices, need at least %d", me, |
245 |
|
|
airSprintSize_t(stmp, nbmat->axis[1].size), minnum); |
246 |
|
|
return 1; |
247 |
|
|
} |
248 |
|
|
|
249 |
|
|
return 0; |
250 |
|
|
} |
251 |
|
|
|
252 |
|
|
/* |
253 |
|
|
******** _tenFindValley |
254 |
|
|
** |
255 |
|
|
** This is not a general purpose function, and it will take some |
256 |
|
|
** work to make it that way. |
257 |
|
|
** |
258 |
|
|
** the tweak argument implements a cheesy heuristic: threshold should be |
259 |
|
|
** on low side of histogram valley, since stdev for background is much |
260 |
|
|
** narrower then stdev for brain |
261 |
|
|
*/ |
262 |
|
|
int |
263 |
|
|
_tenFindValley(double *valP, const Nrrd *nhist, double tweak, int save) { |
264 |
|
|
static const char me[]="_tenFindValley"; |
265 |
|
|
double gparm[NRRD_KERNEL_PARMS_NUM], dparm[NRRD_KERNEL_PARMS_NUM]; |
266 |
|
|
Nrrd *ntmpA, *ntmpB, *nhistD, *nhistDD; |
267 |
|
|
float *hist, *histD, *histDD; |
268 |
|
|
airArray *mop; |
269 |
|
|
size_t bins, maxbb, bb; |
270 |
|
|
NrrdRange *range; |
271 |
|
|
|
272 |
|
|
/* |
273 |
|
|
tenEMBimodalParm *biparm; |
274 |
|
|
biparm = tenEMBimodalParmNew(); |
275 |
|
|
tenEMBimodal(biparm, nhist); |
276 |
|
|
biparm = tenEMBimodalParmNix(biparm); |
277 |
|
|
*/ |
278 |
|
|
|
279 |
|
|
mop = airMopNew(); |
280 |
|
|
airMopAdd(mop, ntmpA=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
281 |
|
|
airMopAdd(mop, ntmpB=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
282 |
|
|
airMopAdd(mop, nhistD=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
283 |
|
|
airMopAdd(mop, nhistDD=nrrdNew(), (airMopper)nrrdNuke, airMopAlways); |
284 |
|
|
|
285 |
|
|
bins = nhist->axis[0].size; |
286 |
|
|
gparm[0] = bins/128; /* wacky heuristic for gaussian stdev */ |
287 |
|
|
gparm[1] = 3; /* how many stdevs to cut-off at */ |
288 |
|
|
dparm[0] = 1.0; /* unit spacing */ |
289 |
|
|
dparm[1] = 1.0; /* B-Spline kernel */ |
290 |
|
|
dparm[2] = 0.0; |
291 |
|
|
if (nrrdCheapMedian(ntmpA, nhist, AIR_TRUE, AIR_FALSE, 2, 1.0, 1024) |
292 |
|
|
|| nrrdSimpleResample(ntmpB, ntmpA, |
293 |
|
|
nrrdKernelGaussian, gparm, &bins, NULL) |
294 |
|
|
|| nrrdSimpleResample(nhistD, ntmpB, |
295 |
|
|
nrrdKernelBCCubicD, dparm, &bins, NULL) |
296 |
|
|
|| nrrdSimpleResample(nhistDD, ntmpB, |
297 |
|
|
nrrdKernelBCCubicDD, dparm, &bins, NULL)) { |
298 |
|
|
biffMovef(TEN, NRRD, "%s: trouble processing histogram", me); |
299 |
|
|
airMopError(mop); return 1; |
300 |
|
|
} |
301 |
|
|
if (save) { |
302 |
|
|
nrrdSave("tmp-histA.nrrd", ntmpA, NULL); |
303 |
|
|
nrrdSave("tmp-histB.nrrd", ntmpB, NULL); |
304 |
|
|
} |
305 |
|
|
hist = (float*)(ntmpB->data); |
306 |
|
|
histD = (float*)(nhistD->data); |
307 |
|
|
histDD = (float*)(nhistDD->data); |
308 |
|
|
range = nrrdRangeNewSet(ntmpB, nrrdBlind8BitRangeState); |
309 |
|
|
airMopAdd(mop, range, (airMopper)nrrdRangeNix, airMopAlways); |
310 |
|
|
for (bb=0; bb<bins-1; bb++) { |
311 |
|
|
if (hist[bb] == range->max) { |
312 |
|
|
/* first seek to max in histogram */ |
313 |
|
|
break; |
314 |
|
|
} |
315 |
|
|
} |
316 |
|
|
maxbb = bb; |
317 |
|
|
for (; bb<bins-1; bb++) { |
318 |
|
|
if (histD[bb]*histD[bb+1] < 0 && histDD[bb] > 0) { |
319 |
|
|
/* zero-crossing in 1st deriv, positive 2nd deriv */ |
320 |
|
|
break; |
321 |
|
|
} |
322 |
|
|
} |
323 |
|
|
if (bb == bins-1) { |
324 |
|
|
biffAddf(TEN, "%s: never saw a satisfactory zero crossing", me); |
325 |
|
|
airMopError(mop); return 1; |
326 |
|
|
} |
327 |
|
|
|
328 |
|
|
*valP = nrrdAxisInfoPos(nhist, 0, AIR_AFFINE(0, tweak, 1, maxbb, bb)); |
329 |
|
|
|
330 |
|
|
airMopOkay(mop); |
331 |
|
|
return 0; |
332 |
|
|
} |
333 |
|
|
|