1 /*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 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
tenEvecRGB(Nrrd * nout,const Nrrd * nin,const tenEvecRGBParm * rgbp)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
tenEvqSingle(float vec[3],float scl)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
tenEvqVolume(Nrrd * nout,const Nrrd * nin,int which,int aniso,int scaleByAniso)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
tenBMatrixCheck(const Nrrd * nbmat,int type,unsigned int minnum)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
_tenFindValley(double * valP,const Nrrd * nhist,double tweak,int save)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
334