Bug Summary

File:src/ten/bvec.c
Location:line 157, column 5
Description:Value stored to 'error' is never read

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
27double
28tenBVecNonLinearFit_error(double *bb, double *ss, double *ww, int len,
29 double amp, double dec) {
30 int ii;
31 double err, tmp;
32
33 err = 0;
34 for (ii=0; ii<len; ii++) {
35 tmp = ww[ii]*(amp*exp(-dec*bb[ii]) - ss[ii]);
36 err += tmp*tmp;
37 }
38 return err;
39}
40
41void
42tenBVecNonLinearFit_linear(double *amp, double *dec,
43 double *bb, double *ss, double *ww, int len) {
44 double x, y, wi=0, xi=0, yi=0, xiyi=0, xisq=0, det;
45 int ii;
46
47 for (ii=0; ii<len; ii++) {
48 x = bb[ii];
49 y = log(AIR_MAX(ss[ii], 0.01)((ss[ii]) > (0.01) ? (ss[ii]) : (0.01)));
50 xi += ww[ii]*x;
51 yi += ww[ii]*y;
52 xiyi += ww[ii]*x*y;
53 xisq += ww[ii]*x*x;
54 wi += ww[ii];
55 }
56 det = xisq*wi - xi*xi;
57 *dec = -(wi*xiyi - xi*yi)/det; /* negative sign assumed in model */
58 *amp = exp((-xi*xiyi + xisq*yi)/det);
59 return;
60}
61
62void
63tenBVecNonLinearFit_GNstep(double *d_amp, double *d_dec,
64 double *bb, double *ss, double *ww, int len,
65 double amp, double dec) {
66 double tmp, ff, dfdx1, dfdx2, AA=0, BB=0, CC=0, JTf[2], det;
67 int ii;
68
69 JTf[0] = JTf[1] = 0;
70 for (ii=0; ii<len; ii++) {
71 tmp = exp(-dec*bb[ii]);
72 ff = ww[ii]*(amp*tmp - ss[ii]);
73 dfdx1 = ww[ii]*tmp;
74 dfdx2 = -ww[ii]*amp*bb[ii]*tmp;
75 AA += dfdx1*dfdx1;
76 BB += dfdx1*dfdx2;
77 CC += dfdx2*dfdx2;
78 JTf[0] += dfdx1*ff;
79 JTf[1] += dfdx2*ff;
80 }
81 det = AA*CC - BB*BB;
82 *d_amp = -(CC*JTf[0] - BB*JTf[1])/det;
83 *d_dec = -(-BB*JTf[0] + AA*JTf[1])/det;
84 return;
85}
86
87
88/*
89******** tenBVecNonLinearFit
90**
91** Assuming that axis 0 represents a sequence of DWI measurements at a
92** range of b values (as described by bb[i]), do non-linear least-squares
93** fitting of those measurements, governed by weights ww[i] (with at
94** most iterMax interations, or terminated when L2 norm change < eps).
95**
96** Based on model fit amp*exp(-b*dec), output nrrd's axis 0 has three values:
97** 0: amp
98** 1: dec
99** 2: error of fit
100** and all other axes are unchanged from input. Output type is always double.
101*/
102int
103tenBVecNonLinearFit(Nrrd *nout, const Nrrd *nin,
104 double *bb, double *ww, int iterMax, double eps) {
105 static const char me[]="tenBVecNonLinearFit";
106 int map[NRRD_DIM_MAX16], vecSize, iter;
107 size_t ii, size[NRRD_DIM_MAX16], vecI, vecNum;
108 char *vec;
109 double *out, ss[AIR_STRLEN_SMALL(128+1)], amp, dec, d_amp, d_dec, error, diff,
110 (*vecLup)(const void *v, size_t I);
111
112 if (!( nout && nin && bb && ww )) {
113 biffAddf(TENtenBiffKey, "%s: got NULL pointer", me);
114 return 1;
115 }
116
117 if (!( nin->dim >= 2 )) {
118 biffAddf(TENtenBiffKey, "%s: nin->dim (%d) not >= 2", me, nin->dim);
119 return 1;
120 }
121 if (!( nin->axis[0].size < AIR_STRLEN_SMALL(128+1) )) {
122 char stmp[AIR_STRLEN_SMALL(128+1)];
123 biffAddf(TENtenBiffKey, "%s: sorry need nin->axis[0].size (%s) < %d", me,
124 airSprintSize_t(stmp, nin->axis[0].size), AIR_STRLEN_SMALL(128+1));
125 return 1;
126 }
127
128 /* allocate/set-up output */
129 nrrdAxisInfoGet_nva(nin, nrrdAxisInfoSize, size);
130 size[0] = 3;
131 if (nrrdMaybeAlloc_nva(nout, nrrdTypeDouble, nin->dim, size)) {
132 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate output", me);
133 return 1;
134 }
135 for (ii=1; ii<nin->dim; ii++) {
136 map[ii] = ii;
137 }
138 map[0] = -1;
139 if (nrrdAxisInfoCopy(nout, nin, map, NRRD_AXIS_INFO_NONE0)) {
140 biffMovef(TENtenBiffKey, NRRDnrrdBiffKey, "%s: couldn't copy axis info", me);
141 return 1;
142 }
143
144 /* process all b vectors */
145 vecSize = nin->axis[0].size*nrrdTypeSize[nin->type];
146 vecNum = nrrdElementNumber(nin)/nin->axis[0].size;
147 vecLup = nrrdDLookup[nin->type];
148 vec = (char*)nin->data;
149 out = (double*)nout->data;
150 for (vecI=0; vecI<vecNum; vecI++) {
151 /* copy DWI signal values */
152 for (ii=0; ii<nin->axis[0].size; ii++) {
153 ss[ii] = vecLup(vec, ii);
154 }
155 /* start with linear fit */
156 tenBVecNonLinearFit_linear(&amp, &dec, bb, ss, ww, nin->axis[0].size);
157 error = tenBVecNonLinearFit_error(bb, ss, ww, nin->axis[0].size, amp, dec);
Value stored to 'error' is never read
158 /* possibly refine with gauss-newton */
159 if (iterMax > 0) {
160 iter = 0;
161 do {
162 iter++;
163 tenBVecNonLinearFit_GNstep(&d_amp, &d_dec,
164 bb, ss, ww, nin->axis[0].size, amp, dec);
165 amp += 0.3*d_amp;
166 dec += 0.3*d_dec;
167 diff = d_amp*d_amp + d_dec*d_dec;
168 } while (iter < iterMax && diff > eps);
169 }
170 error = tenBVecNonLinearFit_error(bb, ss, ww, nin->axis[0].size, amp, dec);
171 out[0] = amp;
172 out[1] = dec;
173 out[2] = error;
174 vec += vecSize;
175 out += 3;
176 }
177
178 return 0;
179}
180