Bug Summary

File:src/ten/experSpec.c
Location:line 236, column 10
Description:Access to field 'data' results in a dereference of a null pointer (loaded from variable 'ngrad')

Annotated Source Code

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
27static int
28_experAlloc(tenExperSpec* espec, unsigned int num) {
29 static char me[]="_experAlloc";
30
31 airFree(espec->bval); espec->bval = NULL((void*)0);
32 airFree(espec->grad); espec->grad = NULL((void*)0);
33 /* espec->wght = airFree(espec->wght); */
34 if (!num) {
35 biffAddf(TENtenBiffKey, "%s: need a non-zero number of images", me);
36 return 1;
37 }
38 espec->imgNum = num;
39 espec->bval = AIR_CALLOC(num, double)(double*)(calloc((num), sizeof(double)));
40 espec->grad = AIR_CALLOC(3*num, double)(double*)(calloc((3*num), sizeof(double)));
41 /* espec->wght = AIR_CALLOC(num, double); */
42 if (!( espec->bval && espec->grad /* && espec->wght */ )) {
43 biffAddf(TENtenBiffKey, "%s: couldn't allocate for %u images", me, num);
44 return 1;
45 }
46 return 0;
47}
48
49tenExperSpec*
50tenExperSpecNew(void) {
51 tenExperSpec* espec;
52
53 espec = AIR_CALLOC(1, tenExperSpec)(tenExperSpec*)(calloc((1), sizeof(tenExperSpec)));
54 espec->set = AIR_FALSE0;
55 espec->imgNum = 0;
56 espec->bval = NULL((void*)0);
57 espec->grad = NULL((void*)0);
58 /* espec->wght = NULL; */
59 return espec;
60}
61
62int
63tenExperSpecGradSingleBValSet(tenExperSpec *espec,
64 int insertB0,
65 double bval,
66 const double *grad,
67 unsigned int gradNum) {
68 static const char me[]="tenExperSpecGradSingleBValSet";
69 unsigned int ii, imgNum, ei;
70
71 if (!espec) {
72 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
73 return 1;
74 }
75 if (insertB0 && !ELL_3V_LEN(grad + 3*0)(sqrt((((grad + 3*0))[0]*((grad + 3*0))[0] + ((grad + 3*0))[1
]*((grad + 3*0))[1] + ((grad + 3*0))[2]*((grad + 3*0))[2])))) {
76 biffAddf(TENtenBiffKey, "%s: wanted insertB0 but gradients "
77 "already start with (0,0,0)", me);
78 return 1;
79 }
80 imgNum = gradNum + !!insertB0;
81 if (_experAlloc(espec, imgNum)) {
82 biffAddf(TENtenBiffKey, "%s: couldn't allocate", me);
83 return 1;
84 }
85 if (insertB0) {
86 espec->bval[0] = 0;
87 ELL_3V_SET(espec->grad + 3*0, 1, 0, 0)((espec->grad + 3*0)[0] = (1), (espec->grad + 3*0)[1] =
(0), (espec->grad + 3*0)[2] = (0));
88 ei = 1;
89 } else {
90 ei = 0;
91 }
92 for (ii=0; ii<gradNum; ei++, ii++) {
93 espec->bval[ei] = bval;
94 ELL_3V_COPY(espec->grad + 3*ei, grad + 3*ii)((espec->grad + 3*ei)[0] = (grad + 3*ii)[0], (espec->grad
+ 3*ei)[1] = (grad + 3*ii)[1], (espec->grad + 3*ei)[2] = (
grad + 3*ii)[2]);
95 /* espec->wght[ii] = 1.0; */
96 }
97
98 return 0;
99}
100
101int
102tenExperSpecGradBValSet(tenExperSpec *espec,
103 int insertB0,
104 const double *bval,
105 const double *grad,
106 unsigned int bgNum) {
107 static const char me[]="tenExperSpecGradBValSet";
108 unsigned int ii, imgNum, ei;
109
110 if (!espec) {
111 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
112 return 1;
113 }
114 if (insertB0 && (!ELL_3V_LEN(grad + 3*0)(sqrt((((grad + 3*0))[0]*((grad + 3*0))[0] + ((grad + 3*0))[1
]*((grad + 3*0))[1] + ((grad + 3*0))[2]*((grad + 3*0))[2]))) || !bval[0])) {
115 biffAddf(TENtenBiffKey, "%s: wanted insertB0 but gradients "
116 "already start with (0,0,0) or bvals start with 0", me);
117 return 1;
118 }
119 imgNum = bgNum + !!insertB0;
120 if (_experAlloc(espec, imgNum)) {
121 biffAddf(TENtenBiffKey, "%s: couldn't allocate", me);
122 return 1;
123 }
124 if (insertB0) {
125 espec->bval[0] = 0;
126 ELL_3V_SET(espec->grad + 3*0, 0, 0, 0)((espec->grad + 3*0)[0] = (0), (espec->grad + 3*0)[1] =
(0), (espec->grad + 3*0)[2] = (0));
127 ei = 1;
128 } else {
129 ei = 0;
130 }
131 for (ii=0; ii<bgNum; ei++, ii++) {
132 espec->bval[ei] = bval[ii];
133 ELL_3V_COPY(espec->grad + 3*ei, grad + 3*ii)((espec->grad + 3*ei)[0] = (grad + 3*ii)[0], (espec->grad
+ 3*ei)[1] = (grad + 3*ii)[1], (espec->grad + 3*ei)[2] = (
grad + 3*ii)[2]);
134 /* espec->wght[ii] = 1.0; */
135 }
136
137 return 0;
138}
139
140/*
141int
142tenExperSpecGradBValWghtSet(tenExperSpec *espec,
143 unsigned int imgNum,
144 const double *bval,
145 const double *grad,
146 const double *wght) {
147 static const char me[]="tenExperSpecGradBValWghtSet";
148 unsigned int ii;
149
150 if (!espec) {
151 biffAddf(TEN, "%s: got NULL pointer", me);
152 return 1;
153 }
154 if (_experAlloc(espec, imgNum)) {
155 biffAddf(TEN, "%s: couldn't allocate", me);
156 return 1;
157 }
158 for (ii=0; ii<imgNum; ii++) {
159 espec->bval[ii] = bval[ii];
160 ELL_3V_COPY(espec->grad + 3*ii, grad + 3*ii);
161 espec->wght[ii] = wght[ii];
162 }
163
164 return 0;
165}
166*/
167
168int
169tenExperSpecFromKeyValueSet(tenExperSpec *espec, const Nrrd *ndwi) {
170 static const char me[]="tenExperSpecFromKeyValueSet";
171 unsigned int *skip, skipNum, ii, imgNum, dwiax;
172 Nrrd *ngrad, *nbmat;
173 airArray *mop;
174 double len, singleBval, *bval, *grad;
175
176 for (dwiax=0; dwiax<ndwi->dim; dwiax++) {
1
Loop condition is true. Entering loop body
→
177 if (nrrdKindList == ndwi->axis[dwiax].kind
2
←
Taking true branch
→
178 || nrrdKindVector == ndwi->axis[dwiax].kind) {
179 break;
3
←
 Execution continues on line 182
→
180 }
181 }
182 if (ndwi->dim == dwiax) {
4
←
Taking false branch
→
183 biffAddf(TENtenBiffKey, "%s: need dwis to have a kind %s or %s axis", me,
184 airEnumStr(nrrdKind, nrrdKindList),
185 airEnumStr(nrrdKind, nrrdKindVector));
186 return 1;
187 } else {
188 if (0 != dwiax) {
5
←
Taking false branch
→
189 biffAddf(TENtenBiffKey, "%s: need dwis (kind %s or %s) along axis 0, not %u", me,
190 airEnumStr(nrrdKind, nrrdKindList),
191 airEnumStr(nrrdKind, nrrdKindVector), dwiax);
192 return 1;
193 }
194 }
195 for (ii=dwiax+1; ii<ndwi->dim; ii++) {
6
←
Loop condition is false. Execution continues on line 201
→
196 if (nrrdKindList == ndwi->axis[ii].kind
197 || nrrdKindVector == ndwi->axis[ii].kind) {
198 break;
199 }
200 }
201 if (ii < ndwi->dim) {
7
←
Taking false branch
→
202 biffAddf(TENtenBiffKey, "%s: saw on %u another %s or %s kind axis, after 0", me,
203 ii, airEnumStr(nrrdKind, nrrdKindList),
204 airEnumStr(nrrdKind, nrrdKindVector));
205 return 1;
206 }
207 if (tenDWMRIKeyValueParse(&ngrad, &nbmat, &singleBval,
8
←
Value assigned to 'ngrad'
→
9
←
Taking false branch
→
208 &skip, &skipNum, ndwi)) {
209 biffAddf(TENtenBiffKey, "%s: trouble parsing DWI info from key/value pairs", me);
210 return 1;
211 }
212 mop = airMopNew();
213 if (ngrad) {
10
←
Assuming 'ngrad' is null
→
11
←
Taking false branch
→
214 airMopAdd(mop, ngrad, (airMopper)nrrdNuke, airMopAlways);
215 }
216 if (nbmat) {
12
←
Assuming 'nbmat' is null
→
13
←
Taking false branch
→
217 airMopAdd(mop, nbmat, (airMopper)nrrdNuke, airMopAlways);
218 }
219 if (skip) {
14
←
Assuming 'skip' is null
→
15
←
Taking false branch
→
220 airMopAdd(mop, skip, airFree, airMopAlways);
221 }
222
223 if (nbmat) {
16
←
Taking false branch
→
224 biffAddf(TENtenBiffKey, "%s: sorry, currently can't handle B-matrices here", me);
225 airMopError(mop); return 1;
226 }
227 if (skipNum) {
17
←
Assuming 'skipNum' is 0
→
18
←
Taking false branch
→
228 biffAddf(TENtenBiffKey, "%s: sorry, currently can't handle skipping (%u) here", me,
229 skipNum);
230 airMopError(mop); return 1;
231 }
232
233 imgNum = ngrad->axis[1].size;
234 bval = AIR_CALLOC(imgNum, double)(double*)(calloc((imgNum), sizeof(double)));
235 airMopAdd(mop, bval, airFree, airMopAlways);
236 grad = AIR_CAST(double *, ngrad->data)((double *)(ngrad->data));
19
←
Within the expansion of the macro 'AIR_CAST':
a
Access to field 'data' results in a dereference of a null pointer (loaded from variable 'ngrad')
237 for (ii=0; ii<imgNum; ii++) {
238 len = ELL_3V_LEN(grad + 3*ii)(sqrt((((grad + 3*ii))[0]*((grad + 3*ii))[0] + ((grad + 3*ii)
)[1]*((grad + 3*ii))[1] + ((grad + 3*ii))[2]*((grad + 3*ii))[
2])));
239 bval[ii] = singleBval*len*len;
240 if (len) {
241 ELL_3V_SCALE(grad + 3*ii, 1/len, grad + 3*ii)((grad + 3*ii)[0] = (1/len)*(grad + 3*ii)[0], (grad + 3*ii)[1
] = (1/len)*(grad + 3*ii)[1], (grad + 3*ii)[2] = (1/len)*(grad
+ 3*ii)[2]);
242 } else {
243 ELL_3V_SET(grad + 3*ii, 0, 0, -1)((grad + 3*ii)[0] = (0), (grad + 3*ii)[1] = (0), (grad + 3*ii
)[2] = (-1));
244 }
245 }
246 if (tenExperSpecGradBValSet(espec, AIR_FALSE0, bval, grad, imgNum)) {
247 biffAddf(TENtenBiffKey, "%s: trouble", me);
248 airMopError(mop); return 1;
249 }
250 airMopOkay(mop);
251 return 0;
252}
253
254tenExperSpec*
255tenExperSpecNix(tenExperSpec *espec) {
256
257 if (espec) {
258 airFree(espec->bval);
259 airFree(espec->grad);
260 /* espec->wght = airFree(espec->wght); */
261 airFree(espec);
262 }
263 return NULL((void*)0);
264}
265
266double
267_tenExperSpec_sqe(const double *dwiMeas, const double *dwiSim,
268 const tenExperSpec *espec, int knownB0) {
269 unsigned int ii;
270 double sqe;
271
272 sqe = 0;
273 if (knownB0) {
274 for (ii=0; ii<espec->imgNum; ii++) {
275 double dd;
276 if (!espec->bval[ii]) {
277 continue;
278 }
279 dd = dwiMeas[ii] - dwiSim[ii];
280 sqe += dd*dd;
281 /*
282 fprintf(stderr, "!%s: dwi[%u]: %g - %g -> %g\n", "_tenExperSpec_sqe",
283 ii, dwiMeas[ii], dwiSim[ii], sqe);
284 */
285 }
286 } else {
287 for (ii=0; ii<espec->imgNum; ii++) {
288 double dd;
289 dd = dwiMeas[ii] - dwiSim[ii];
290 sqe += dd*dd;
291 }
292 }
293 return sqe;
294}
295
296double
297_tenExperSpec_nll(const double *dwiMeas, const double *dwiSim,
298 const tenExperSpec *espec,
299 int rician, double sigma, int knownB0) {
300 double nll;
301 unsigned int ii;
302
303 nll = 0;
304 if (rician) {
305 for (ii=0; ii<espec->imgNum; ii++) {
306 if (knownB0 && !espec->bval[ii]) {
307 continue;
308 }
309 nll += -airLogRician(dwiMeas[ii], dwiSim[ii], sigma);
310 }
311 } else {
312 double dd, ladd, denom;
313 ladd = log(sigma*sqrt(2*AIR_PI3.14159265358979323846));
314 denom = 1.0/(2*sigma*sigma);
315 for (ii=0; ii<espec->imgNum; ii++) {
316 if (knownB0 && !espec->bval[ii]) {
317 continue;
318 }
319 dd = dwiMeas[ii] - dwiSim[ii];
320 nll += dd*dd*denom + ladd;
321 }
322 }
323 return nll;
324}
325
326int
327tenDWMRIKeyValueFromExperSpecSet(Nrrd *ndwi, const tenExperSpec *espec) {
328 static char me[]="tenDWMRIKeyValueFromExperSpecSet";
329 char keystr[AIR_STRLEN_MED(256+1)], valstr[AIR_STRLEN_MED(256+1)];
330 double maxb, bb;
331 unsigned int ii;
332
333 if (!(ndwi && espec)) {
334 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
335 return 1;
336 }
337
338 nrrdKeyValueAdd(ndwi, tenDWMRIModalityKey, tenDWMRIModalityVal);
339 maxb = tenExperSpecMaxBGet(espec);
340 sprintf(valstr, "%.17g", maxb)__builtin___sprintf_chk (valstr, 0, __builtin_object_size (valstr
, 2 > 1 ? 1 : 0), "%.17g", maxb);
341 nrrdKeyValueAdd(ndwi, tenDWMRIBValueKey, valstr);
342 for (ii=0; ii<espec->imgNum; ii++) {
343 double vec[3];
344 sprintf(keystr, tenDWMRIGradKeyFmt, ii)__builtin___sprintf_chk (keystr, 0, __builtin_object_size (keystr
, 2 > 1 ? 1 : 0), tenDWMRIGradKeyFmt, ii);
345 ELL_3V_COPY(vec, espec->grad + 3*ii)((vec)[0] = (espec->grad + 3*ii)[0], (vec)[1] = (espec->
grad + 3*ii)[1], (vec)[2] = (espec->grad + 3*ii)[2]);
346 bb = espec->bval[ii];
347 /* Thu Dec 20 03:25:20 CST 2012 this rescaling is not, btw,
348 what is causing the small discrepency between ngrad before
349 and after saving to KVPs */
350 ELL_3V_SCALE(vec, sqrt(bb/maxb), vec)((vec)[0] = (sqrt(bb/maxb))*(vec)[0], (vec)[1] = (sqrt(bb/maxb
))*(vec)[1], (vec)[2] = (sqrt(bb/maxb))*(vec)[2]);
351 sprintf(valstr, "%.17g %.17g %.17g", vec[0], vec[1], vec[2])__builtin___sprintf_chk (valstr, 0, __builtin_object_size (valstr
, 2 > 1 ? 1 : 0), "%.17g %.17g %.17g", vec[0], vec[1], vec
[2]);
352 nrrdKeyValueAdd(ndwi, keystr, valstr);
353 }
354 /* HEY what if its a full B-matrix? */
355
356 return 0;
357}
358
359/*
360** learns B0 from DWIs by simple averaging of all the dwi[ii]
361** without any diffusion weighting, as indicated by espec->bval[ii],
362** or, returns AIR_NAN when there are no such dwi[ii]
363*/
364double
365tenExperSpecKnownB0Get(const tenExperSpec *espec, const double *dwi) {
366 unsigned int ii, nb;
367 double ret, b0;
368
369 if (!( dwi && espec )) {
370 return AIR_NAN(airFloatQNaN.f);
371 }
372
373 nb = 0;
374 b0 = 0.0;
375 for (ii=0; ii<espec->imgNum; ii++) {
376 if (0 == espec->bval[ii]) {
377 b0 += dwi[ii];
378 ++nb;
379 }
380 }
381 if (nb) {
382 ret = b0/nb;
383 } else {
384 ret = AIR_NAN(airFloatQNaN.f);
385 }
386 return ret;
387}
388
389double
390tenExperSpecMaxBGet(const tenExperSpec *espec) {
391 unsigned int ii;
392 double bval;
393
394 if (!( espec )) {
395 return AIR_NAN(airFloatQNaN.f);
396 }
397
398 bval = -1;
399 for (ii=0; ii<espec->imgNum; ii++) {
400 bval = AIR_MAX(bval, espec->bval[ii])((bval) > (espec->bval[ii]) ? (bval) : (espec->bval[
ii]));
401 }
402 return bval;
403}