Bug Summary

File:src/gage/sclanswer.c
Location:line 226, column 7
Description:Value stored to 'pness' 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 "gage.h"
25#include "privateGage.h"
26
27/* HEY copied from ten.h */
28#define TEN_M2T(t, m) ( \
29 (t)[1] = (m)[0], \
30 (t)[2] = ((m)[1]+(m)[3])/2.0, \
31 (t)[3] = ((m)[2]+(m)[6])/2.0, \
32 (t)[4] = (m)[4], \
33 (t)[5] = ((m)[5]+(m)[7])/2.0, \
34 (t)[6] = (m)[8])
35#define TEN_T_SCALE(a, s, b) ( \
36 (a)[0] = (b)[0], \
37 (a)[1] = (s)*(b)[1], \
38 (a)[2] = (s)*(b)[2], \
39 (a)[3] = (s)*(b)[3], \
40 (a)[4] = (s)*(b)[4], \
41 (a)[5] = (s)*(b)[5], \
42 (a)[6] = (s)*(b)[6])
43
44void
45_gageSclAnswer(gageContext *ctx, gagePerVolume *pvl) {
46 char me[]="_gageSclAnswer";
47 double gmag=0, *hess, *norm, *gvec, *gten, *k1, *k2, curv=0,
48 sHess[9]={0,0,0,0,0,0,0,0,0};
49 double tmpMat[9], tmpVec[3], hevec[9], heval[3];
50 double len, gp1[3], gp2[3], *nPerp, ncTen[9], nProj[9]={0,0,0,0,0,0,0,0,0};
51 double alpha = 0.5;
52 double beta = 0.5;
53 double _gamma = 5;
54 double cc = 1e-6;
55#define FD_MEDIAN_MAX16 16
56 int fd, nidx, xi, yi, zi;
57 double *fw, iv3wght[2*FD_MEDIAN_MAX16*FD_MEDIAN_MAX16*FD_MEDIAN_MAX16],
58 wghtSum, wght;
59
60 /* convenience pointers for work below */
61 hess = pvl->directAnswer[gageSclHessian];
62 gvec = pvl->directAnswer[gageSclGradVec];
63 norm = pvl->directAnswer[gageSclNormal];
64 nPerp = pvl->directAnswer[gageSclNPerp];
65 gten = pvl->directAnswer[gageSclGeomTens];
66 k1 = pvl->directAnswer[gageSclK1];
67 k2 = pvl->directAnswer[gageSclK2];
68
69 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclValue)(pvl->query[gageSclValue/8] & (1 << (gageSclValue
% 8)))) {
70 /* done if doV */
71 if (ctx->verbose > 2) {
72 fprintf(stderr__stderrp, "%s: val = % 15.7f\n", me,
73 (double)(pvl->directAnswer[gageSclValue][0]));
74 }
75 }
76 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclGradVec)(pvl->query[gageSclGradVec/8] & (1 << (gageSclGradVec
% 8)))) {
77 /* done if doD1 */
78 if (ctx->verbose > 2) {
79 fprintf(stderr__stderrp, "%s: gvec = ", me);
80 ell_3v_print_d(stderr__stderrp, gvec);
81 }
82 }
83 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclGradMag)(pvl->query[gageSclGradMag/8] & (1 << (gageSclGradMag
% 8)))) {
84 /* this is the true value of gradient magnitude */
85 gmag = pvl->directAnswer[gageSclGradMag][0] = sqrt(ELL_3V_DOT(gvec, gvec)((gvec)[0]*(gvec)[0] + (gvec)[1]*(gvec)[1] + (gvec)[2]*(gvec)
[2]));
86 }
87
88 /* NB: it would seem that gageParmGradMagMin is completely ignored . . . */
89
90 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclNormal)(pvl->query[gageSclNormal/8] & (1 << (gageSclNormal
% 8)))) {
91 if (gmag) {
92 ELL_3V_SCALE(norm, 1/gmag, gvec)((norm)[0] = (1/gmag)*(gvec)[0], (norm)[1] = (1/gmag)*(gvec)[
1], (norm)[2] = (1/gmag)*(gvec)[2]);
93 /* polishing
94 len = sqrt(ELL_3V_DOT(norm, norm));
95 ELL_3V_SCALE(norm, 1/len, norm);
96 */
97 } else {
98 ELL_3V_COPY(norm, gageZeroNormal)((norm)[0] = (gageZeroNormal)[0], (norm)[1] = (gageZeroNormal
)[1], (norm)[2] = (gageZeroNormal)[2]);
99 }
100 }
101 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclNProj)(pvl->query[gageSclNProj/8] & (1 << (gageSclNProj
% 8)))) {
102 ELL_3MV_OUTER(pvl->directAnswer[gageSclNProj], norm, norm)((((pvl->directAnswer[gageSclNProj])+0)[0] = ((norm)[0])*(
(norm))[0], ((pvl->directAnswer[gageSclNProj])+0)[1] = ((norm
)[0])*((norm))[1], ((pvl->directAnswer[gageSclNProj])+0)[2
] = ((norm)[0])*((norm))[2]), (((pvl->directAnswer[gageSclNProj
])+3)[0] = ((norm)[1])*((norm))[0], ((pvl->directAnswer[gageSclNProj
])+3)[1] = ((norm)[1])*((norm))[1], ((pvl->directAnswer[gageSclNProj
])+3)[2] = ((norm)[1])*((norm))[2]), (((pvl->directAnswer[
gageSclNProj])+6)[0] = ((norm)[2])*((norm))[0], ((pvl->directAnswer
[gageSclNProj])+6)[1] = ((norm)[2])*((norm))[1], ((pvl->directAnswer
[gageSclNProj])+6)[2] = ((norm)[2])*((norm))[2]));
103 }
104 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclNPerp)(pvl->query[gageSclNPerp/8] & (1 << (gageSclNPerp
% 8)))) {
105 /* nPerp = I - outer(norm, norm) */
106 /* NB: this sets both nPerp and nProj */
107 ELL_3MV_OUTER(nProj, norm, norm)((((nProj)+0)[0] = ((norm)[0])*((norm))[0], ((nProj)+0)[1] = (
(norm)[0])*((norm))[1], ((nProj)+0)[2] = ((norm)[0])*((norm))
[2]), (((nProj)+3)[0] = ((norm)[1])*((norm))[0], ((nProj)+3)[
1] = ((norm)[1])*((norm))[1], ((nProj)+3)[2] = ((norm)[1])*((
norm))[2]), (((nProj)+6)[0] = ((norm)[2])*((norm))[0], ((nProj
)+6)[1] = ((norm)[2])*((norm))[1], ((nProj)+6)[2] = ((norm)[2
])*((norm))[2]));
108 ELL_3M_SCALE(nPerp, -1, nProj)((((nPerp)+0)[0] = ((-1))*((nProj)+0)[0], ((nPerp)+0)[1] = ((
-1))*((nProj)+0)[1], ((nPerp)+0)[2] = ((-1))*((nProj)+0)[2]),
(((nPerp)+3)[0] = ((-1))*((nProj)+3)[0], ((nPerp)+3)[1] = ((
-1))*((nProj)+3)[1], ((nPerp)+3)[2] = ((-1))*((nProj)+3)[2]),
(((nPerp)+6)[0] = ((-1))*((nProj)+6)[0], ((nPerp)+6)[1] = ((
-1))*((nProj)+6)[1], ((nPerp)+6)[2] = ((-1))*((nProj)+6)[2]));
109 nPerp[0] += 1;
110 nPerp[4] += 1;
111 nPerp[8] += 1;
112 }
113 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessian)(pvl->query[gageSclHessian/8] & (1 << (gageSclHessian
% 8)))) {
114 /* done if doD2 */
115 if (ctx->verbose > 2) {
116 fprintf(stderr__stderrp, "%s: hess = \n", me);
117 ell_3m_print_d(stderr__stderrp, hess);
118 }
119 }
120 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessianTen)(pvl->query[gageSclHessianTen/8] & (1 << (gageSclHessianTen
% 8)))) {
121 pvl->directAnswer[gageSclHessianTen][0] = 1.0;
122 TEN_M2T(pvl->directAnswer[gageSclHessianTen],
123 pvl->directAnswer[gageSclHessian]);
124 }
125 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclLaplacian)(pvl->query[gageSclLaplacian/8] & (1 << (gageSclLaplacian
% 8)))) {
126 pvl->directAnswer[gageSclLaplacian][0] = hess[0] + hess[4] + hess[8];
127 if (ctx->verbose > 2) {
128 fprintf(stderr__stderrp, "%s: lapl = %g + %g + %g = %g\n", me,
129 hess[0], hess[4], hess[8],
130 pvl->directAnswer[gageSclLaplacian][0]);
131 }
132 }
133 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessFrob)(pvl->query[gageSclHessFrob/8] & (1 << (gageSclHessFrob
% 8)))) {
134 pvl->directAnswer[gageSclHessFrob][0] = ELL_3M_FROB(hess)(sqrt((((hess)+0)[0]*((hess)+0)[0] + ((hess)+0)[1]*((hess)+0)
[1] + ((hess)+0)[2]*((hess)+0)[2]) + (((hess)+3)[0]*((hess)+3
)[0] + ((hess)+3)[1]*((hess)+3)[1] + ((hess)+3)[2]*((hess)+3)
[2]) + (((hess)+6)[0]*((hess)+6)[0] + ((hess)+6)[1]*((hess)+6
)[1] + ((hess)+6)[2]*((hess)+6)[2])));
135 }
136 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessEval)(pvl->query[gageSclHessEval/8] & (1 << (gageSclHessEval
% 8)))) {
137 if (ctx->parm.twoDimZeroZ) {
138 ell_3m2sub_eigensolve_d(heval, hevec, hess);
139 } else {
140 /* HEY: look at the return value for root multiplicity? */
141 ell_3m_eigensolve_d(heval, hevec, hess, AIR_TRUE1);
142 }
143 ELL_3V_COPY(pvl->directAnswer[gageSclHessEval], heval)((pvl->directAnswer[gageSclHessEval])[0] = (heval)[0], (pvl
->directAnswer[gageSclHessEval])[1] = (heval)[1], (pvl->
directAnswer[gageSclHessEval])[2] = (heval)[2]);
144 }
145 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessEvec)(pvl->query[gageSclHessEvec/8] & (1 << (gageSclHessEvec
% 8)))) {
146 ELL_3M_COPY(pvl->directAnswer[gageSclHessEvec], hevec)((((pvl->directAnswer[gageSclHessEvec])+0)[0] = ((hevec)+0
)[0], ((pvl->directAnswer[gageSclHessEvec])+0)[1] = ((hevec
)+0)[1], ((pvl->directAnswer[gageSclHessEvec])+0)[2] = ((hevec
)+0)[2]), (((pvl->directAnswer[gageSclHessEvec])+3)[0] = (
(hevec)+3)[0], ((pvl->directAnswer[gageSclHessEvec])+3)[1]
= ((hevec)+3)[1], ((pvl->directAnswer[gageSclHessEvec])+3
)[2] = ((hevec)+3)[2]), (((pvl->directAnswer[gageSclHessEvec
])+6)[0] = ((hevec)+6)[0], ((pvl->directAnswer[gageSclHessEvec
])+6)[1] = ((hevec)+6)[1], ((pvl->directAnswer[gageSclHessEvec
])+6)[2] = ((hevec)+6)[2]));
147 }
148#if 1
149 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessRidgeness)(pvl->query[gageSclHessRidgeness/8] & (1 << (gageSclHessRidgeness
% 8)))) {
150 double A, B, S;
151 if (heval[1] >0 || heval[2]>0) {
152 pvl->directAnswer[gageSclHessRidgeness][0] = 0;
153 }
154 else if (AIR_ABS(heval[1])((heval[1]) > 0.0f ? (heval[1]) : -(heval[1]))<1e-10 || AIR_ABS(heval[2])((heval[2]) > 0.0f ? (heval[2]) : -(heval[2]))<1e-10) {
155 pvl->directAnswer[gageSclHessRidgeness][0] = 0;
156 }
157 else {
158 double *ans;
159 A = AIR_ABS(heval[1])((heval[1]) > 0.0f ? (heval[1]) : -(heval[1]))/AIR_ABS(heval[2])((heval[2]) > 0.0f ? (heval[2]) : -(heval[2]));
160 B = AIR_ABS(heval[0])((heval[0]) > 0.0f ? (heval[0]) : -(heval[0]))/sqrt(AIR_ABS(heval[1]*heval[2])((heval[1]*heval[2]) > 0.0f ? (heval[1]*heval[2]) : -(heval
[1]*heval[2])));
161 S = sqrt(heval[0]*heval[0] + heval[1]*heval[1] + heval[2]*heval[2]);
162 ans = pvl->directAnswer[gageSclHessRidgeness];
163 ans[0] = (1-exp(-A*A/(2*alpha*alpha))) *
164 exp(-B*B/(2*beta*beta)) *
165 (1-exp(-S*S/(2*_gamma*_gamma))) *
166 exp(-2*cc*cc/(AIR_ABS(heval[1])((heval[1]) > 0.0f ? (heval[1]) : -(heval[1]))*heval[2]*heval[2]));
167 }
168 }
169#else
170 /* alternative implementation by GLK, based on directly following
171 Frangi text. Only significant difference from above is a
172 discontinuity at heval[0] = -heval[1] */
173 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessRidgeness)(pvl->query[gageSclHessRidgeness/8] & (1 << (gageSclHessRidgeness
% 8)))) {
174 double ev[4], tmp;
175 ELL_3V_COPY(ev+1, heval)((ev+1)[0] = (heval)[0], (ev+1)[1] = (heval)[1], (ev+1)[2] = (
heval)[2]);
176 if (AIR_ABS(ev[2])((ev[2]) > 0.0f ? (ev[2]) : -(ev[2])) > AIR_ABS(ev[3])((ev[3]) > 0.0f ? (ev[3]) : -(ev[3]))) { ELL_SWAP2(ev[2], ev[3], tmp)((tmp)=(ev[2]),(ev[2])=(ev[3]),(ev[3])=(tmp)); }
177 if (AIR_ABS(ev[1])((ev[1]) > 0.0f ? (ev[1]) : -(ev[1])) > AIR_ABS(ev[2])((ev[2]) > 0.0f ? (ev[2]) : -(ev[2]))) { ELL_SWAP2(ev[1], ev[2], tmp)((tmp)=(ev[1]),(ev[1])=(ev[2]),(ev[2])=(tmp)); }
178 if (AIR_ABS(ev[2])((ev[2]) > 0.0f ? (ev[2]) : -(ev[2])) > AIR_ABS(ev[3])((ev[3]) > 0.0f ? (ev[3]) : -(ev[3]))) { ELL_SWAP2(ev[2], ev[3], tmp)((tmp)=(ev[2]),(ev[2])=(ev[3]),(ev[3])=(tmp)); }
179 if (ev[2] > 0 || ev[3] > 0) {
180 pvl->directAnswer[gageSclHessRidgeness][0] = 0;
181 } else {
182 double a1, a2, a3, RB, RA, SS, fa, fb, fg, aa, bb, gg, frangi;
183 a1 = AIR_ABS(ev[1])((ev[1]) > 0.0f ? (ev[1]) : -(ev[1]));
184 a2 = AIR_ABS(ev[2])((ev[2]) > 0.0f ? (ev[2]) : -(ev[2]));
185 a3 = AIR_ABS(ev[3])((ev[3]) > 0.0f ? (ev[3]) : -(ev[3]));
186 RB = a1/sqrt(a2*a3);
187 RA = a2/a3;
188 SS = sqrt(a1*a1 + a2*a2 + a3*a3);
189 aa = bb = 0.5;
190 gg = 1;
191 fa = 1 - exp(-RA*RA/(2*aa*aa));
192 fb = exp(-RB*RB/(2*bb*bb));
193 fg = 1 - exp(-SS*SS/(2*gg*gg));
194 frangi = fa*fb*fg;
195 if (!AIR_EXISTS(frangi)(((int)(!((frangi) - (frangi)))))) { frangi = 0.0; }
196 pvl->directAnswer[gageSclHessRidgeness][0] = frangi;
197 }
198 }
199#endif
200 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessValleyness)(pvl->query[gageSclHessValleyness/8] & (1 << (gageSclHessValleyness
% 8)))) {
201 double A, B, S;
202 if (heval[0] <0 || heval[1]<0) {
203 pvl->directAnswer[gageSclHessValleyness][0] = 0;
204 }
205 else if (AIR_ABS(heval[0])((heval[0]) > 0.0f ? (heval[0]) : -(heval[0]))<1e-10 || AIR_ABS(heval[1])((heval[1]) > 0.0f ? (heval[1]) : -(heval[1]))<1e-10) {
206 pvl->directAnswer[gageSclHessValleyness][0] = 0;
207 }
208 else {
209 double *ans;
210 A = AIR_ABS(heval[1])((heval[1]) > 0.0f ? (heval[1]) : -(heval[1]))/AIR_ABS(heval[0])((heval[0]) > 0.0f ? (heval[0]) : -(heval[0]));
211 B = AIR_ABS(heval[2])((heval[2]) > 0.0f ? (heval[2]) : -(heval[2]))/sqrt(AIR_ABS(heval[1]*heval[0])((heval[1]*heval[0]) > 0.0f ? (heval[1]*heval[0]) : -(heval
[1]*heval[0])));
212 S = sqrt(heval[0]*heval[0] + heval[1]*heval[1] + heval[2]*heval[2]);
213 ans = pvl->directAnswer[gageSclHessValleyness];
214 ans[0] = (1-exp(-A*A/(2*alpha*alpha))) *
215 exp(-B*B/(2*beta*beta)) *
216 (1-exp(-S*S/(2*_gamma*_gamma))) *
217 exp(-2*cc*cc/(AIR_ABS(heval[1])((heval[1]) > 0.0f ? (heval[1]) : -(heval[1]))*heval[0]*heval[0]));
218 }
219 }
220 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessDotPeakness)(pvl->query[gageSclHessDotPeakness/8] & (1 << (gageSclHessDotPeakness
% 8)))) {
221#define OSQT 0.57735026918962576451
222 double neval[3], hlen, pness,
223 peak[3] = {-OSQT, -OSQT, -OSQT};
224 ELL_3V_NORM(neval, heval, hlen)(hlen = (sqrt((((heval))[0]*((heval))[0] + ((heval))[1]*((heval
))[1] + ((heval))[2]*((heval))[2]))), ((neval)[0] = (1.0/hlen
)*(heval)[0], (neval)[1] = (1.0/hlen)*(heval)[1], (neval)[2] =
(1.0/hlen)*(heval)[2]));
225 if (!hlen) {
226 pness = 0;
Value stored to 'pness' is never read
227 } else {
228 pness = ELL_3V_DOT(peak, neval)((peak)[0]*(neval)[0] + (peak)[1]*(neval)[1] + (peak)[2]*(neval
)[2]);
229 pness = AIR_AFFINE(-1, pness, 1, 0, 1)( ((double)(1)-(0))*((double)(pness)-(-1)) / ((double)(1)-(-1
)) + (0));
230 pvl->directAnswer[gageSclHessDotPeakness][0] = pness;
231 }
232#undef OSQT
233 }
234 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclHessMode)(pvl->query[gageSclHessMode/8] & (1 << (gageSclHessMode
% 8)))) {
235 pvl->directAnswer[gageSclHessMode][0] = airMode3_d(heval);
236 }
237
238 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageScl2ndDD)(pvl->query[gageScl2ndDD/8] & (1 << (gageScl2ndDD
% 8)))) {
239 ELL_3MV_MUL(tmpVec, hess, norm)((tmpVec)[0] = (hess)[0]*(norm)[0] + (hess)[1]*(norm)[1] + (hess
)[2]*(norm)[2], (tmpVec)[1] = (hess)[3]*(norm)[0] + (hess)[4]
*(norm)[1] + (hess)[5]*(norm)[2], (tmpVec)[2] = (hess)[6]*(norm
)[0] + (hess)[7]*(norm)[1] + (hess)[8]*(norm)[2]);
240 pvl->directAnswer[gageScl2ndDD][0] = ELL_3V_DOT(norm, tmpVec)((norm)[0]*(tmpVec)[0] + (norm)[1]*(tmpVec)[1] + (norm)[2]*(tmpVec
)[2]);
241 }
242 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclGeomTens)(pvl->query[gageSclGeomTens/8] & (1 << (gageSclGeomTens
% 8)))) {
243 if (gmag > ctx->parm.gradMagCurvMin) {
244 /* parm.curvNormalSide applied here to determine the sense of the
245 normal when doing all curvature calculations */
246 ELL_3M_SCALE(sHess, -(ctx->parm.curvNormalSide)/gmag, hess)((((sHess)+0)[0] = ((-(ctx->parm.curvNormalSide)/gmag))*((
hess)+0)[0], ((sHess)+0)[1] = ((-(ctx->parm.curvNormalSide
)/gmag))*((hess)+0)[1], ((sHess)+0)[2] = ((-(ctx->parm.curvNormalSide
)/gmag))*((hess)+0)[2]), (((sHess)+3)[0] = ((-(ctx->parm.curvNormalSide
)/gmag))*((hess)+3)[0], ((sHess)+3)[1] = ((-(ctx->parm.curvNormalSide
)/gmag))*((hess)+3)[1], ((sHess)+3)[2] = ((-(ctx->parm.curvNormalSide
)/gmag))*((hess)+3)[2]), (((sHess)+6)[0] = ((-(ctx->parm.curvNormalSide
)/gmag))*((hess)+6)[0], ((sHess)+6)[1] = ((-(ctx->parm.curvNormalSide
)/gmag))*((hess)+6)[1], ((sHess)+6)[2] = ((-(ctx->parm.curvNormalSide
)/gmag))*((hess)+6)[2]));
247
248 /* gten = nPerp * sHess * nPerp */
249 ELL_3M_MUL(tmpMat, sHess, nPerp)((tmpMat)[0] = (sHess)[0]*(nPerp)[0] + (sHess)[1]*(nPerp)[3] +
(sHess)[2]*(nPerp)[6], (tmpMat)[1] = (sHess)[0]*(nPerp)[1] +
(sHess)[1]*(nPerp)[4] + (sHess)[2]*(nPerp)[7], (tmpMat)[2] =
(sHess)[0]*(nPerp)[2] + (sHess)[1]*(nPerp)[5] + (sHess)[2]*(
nPerp)[8], (tmpMat)[3] = (sHess)[3]*(nPerp)[0] + (sHess)[4]*(
nPerp)[3] + (sHess)[5]*(nPerp)[6], (tmpMat)[4] = (sHess)[3]*(
nPerp)[1] + (sHess)[4]*(nPerp)[4] + (sHess)[5]*(nPerp)[7], (tmpMat
)[5] = (sHess)[3]*(nPerp)[2] + (sHess)[4]*(nPerp)[5] + (sHess
)[5]*(nPerp)[8], (tmpMat)[6] = (sHess)[6]*(nPerp)[0] + (sHess
)[7]*(nPerp)[3] + (sHess)[8]*(nPerp)[6], (tmpMat)[7] = (sHess
)[6]*(nPerp)[1] + (sHess)[7]*(nPerp)[4] + (sHess)[8]*(nPerp)[
7], (tmpMat)[8] = (sHess)[6]*(nPerp)[2] + (sHess)[7]*(nPerp)[
5] + (sHess)[8]*(nPerp)[8]);
250 ELL_3M_MUL(gten, nPerp, tmpMat)((gten)[0] = (nPerp)[0]*(tmpMat)[0] + (nPerp)[1]*(tmpMat)[3] +
(nPerp)[2]*(tmpMat)[6], (gten)[1] = (nPerp)[0]*(tmpMat)[1] +
(nPerp)[1]*(tmpMat)[4] + (nPerp)[2]*(tmpMat)[7], (gten)[2] =
(nPerp)[0]*(tmpMat)[2] + (nPerp)[1]*(tmpMat)[5] + (nPerp)[2]
*(tmpMat)[8], (gten)[3] = (nPerp)[3]*(tmpMat)[0] + (nPerp)[4]
*(tmpMat)[3] + (nPerp)[5]*(tmpMat)[6], (gten)[4] = (nPerp)[3]
*(tmpMat)[1] + (nPerp)[4]*(tmpMat)[4] + (nPerp)[5]*(tmpMat)[7
], (gten)[5] = (nPerp)[3]*(tmpMat)[2] + (nPerp)[4]*(tmpMat)[5
] + (nPerp)[5]*(tmpMat)[8], (gten)[6] = (nPerp)[6]*(tmpMat)[0
] + (nPerp)[7]*(tmpMat)[3] + (nPerp)[8]*(tmpMat)[6], (gten)[7
] = (nPerp)[6]*(tmpMat)[1] + (nPerp)[7]*(tmpMat)[4] + (nPerp)
[8]*(tmpMat)[7], (gten)[8] = (nPerp)[6]*(tmpMat)[2] + (nPerp)
[7]*(tmpMat)[5] + (nPerp)[8]*(tmpMat)[8]);
251
252 if (ctx->verbose > 2) {
253 fprintf(stderr__stderrp, "%s: gten: \n", me);
254 ell_3m_print_d(stderr__stderrp, gten);
255 ELL_3MV_MUL(tmpVec, gten, norm)((tmpVec)[0] = (gten)[0]*(norm)[0] + (gten)[1]*(norm)[1] + (gten
)[2]*(norm)[2], (tmpVec)[1] = (gten)[3]*(norm)[0] + (gten)[4]
*(norm)[1] + (gten)[5]*(norm)[2], (tmpVec)[2] = (gten)[6]*(norm
)[0] + (gten)[7]*(norm)[1] + (gten)[8]*(norm)[2]);
256 len = ELL_3V_LEN(tmpVec)(sqrt((((tmpVec))[0]*((tmpVec))[0] + ((tmpVec))[1]*((tmpVec))
[1] + ((tmpVec))[2]*((tmpVec))[2])));
257 fprintf(stderr__stderrp, "%s: should be small: %30.15f\n", me, (double)len);
258 ell_3v_perp_d(gp1, norm);
259 ELL_3MV_MUL(tmpVec, gten, gp1)((tmpVec)[0] = (gten)[0]*(gp1)[0] + (gten)[1]*(gp1)[1] + (gten
)[2]*(gp1)[2], (tmpVec)[1] = (gten)[3]*(gp1)[0] + (gten)[4]*(
gp1)[1] + (gten)[5]*(gp1)[2], (tmpVec)[2] = (gten)[6]*(gp1)[0
] + (gten)[7]*(gp1)[1] + (gten)[8]*(gp1)[2]);
260 len = ELL_3V_LEN(tmpVec)(sqrt((((tmpVec))[0]*((tmpVec))[0] + ((tmpVec))[1]*((tmpVec))
[1] + ((tmpVec))[2]*((tmpVec))[2])));
261 fprintf(stderr__stderrp, "%s: should be bigger: %30.15f\n", me, (double)len);
262 ELL_3V_CROSS(gp2, gp1, norm)((gp2)[0] = (gp1)[1]*(norm)[2] - (gp1)[2]*(norm)[1], (gp2)[1]
= (gp1)[2]*(norm)[0] - (gp1)[0]*(norm)[2], (gp2)[2] = (gp1)[
0]*(norm)[1] - (gp1)[1]*(norm)[0]);
263 ELL_3MV_MUL(tmpVec, gten, gp2)((tmpVec)[0] = (gten)[0]*(gp2)[0] + (gten)[1]*(gp2)[1] + (gten
)[2]*(gp2)[2], (tmpVec)[1] = (gten)[3]*(gp2)[0] + (gten)[4]*(
gp2)[1] + (gten)[5]*(gp2)[2], (tmpVec)[2] = (gten)[6]*(gp2)[0
] + (gten)[7]*(gp2)[1] + (gten)[8]*(gp2)[2]);
264 len = ELL_3V_LEN(tmpVec)(sqrt((((tmpVec))[0]*((tmpVec))[0] + ((tmpVec))[1]*((tmpVec))
[1] + ((tmpVec))[2]*((tmpVec))[2])));
265 fprintf(stderr__stderrp, "%s: should (also) be bigger: %30.15f\n",
266 me, (double)len);
267 }
268 } else {
269 ELL_3M_ZERO_SET(gten)((((gten)+0)[0] = (0), ((gten)+0)[1] = (0), ((gten)+0)[2] = (
0)), (((gten)+3)[0] = (0), ((gten)+3)[1] = (0), ((gten)+3)[2]
= (0)), (((gten)+6)[0] = (0), ((gten)+6)[1] = (0), ((gten)+6
)[2] = (0)));
270 }
271 }
272 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclGeomTensTen)(pvl->query[gageSclGeomTensTen/8] & (1 << (gageSclGeomTensTen
% 8)))) {
273 pvl->directAnswer[gageSclGeomTensTen][0] = 1.0;
274 TEN_M2T(pvl->directAnswer[gageSclGeomTensTen],
275 pvl->directAnswer[gageSclGeomTens]);
276 TEN_T_SCALE(pvl->directAnswer[gageSclGeomTensTen], -1,
277 pvl->directAnswer[gageSclGeomTensTen]);
278 }
279 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclTotalCurv)(pvl->query[gageSclTotalCurv/8] & (1 << (gageSclTotalCurv
% 8)))) {
280 curv = pvl->directAnswer[gageSclTotalCurv][0] = ELL_3M_FROB(gten)(sqrt((((gten)+0)[0]*((gten)+0)[0] + ((gten)+0)[1]*((gten)+0)
[1] + ((gten)+0)[2]*((gten)+0)[2]) + (((gten)+3)[0]*((gten)+3
)[0] + ((gten)+3)[1]*((gten)+3)[1] + ((gten)+3)[2]*((gten)+3)
[2]) + (((gten)+6)[0]*((gten)+6)[0] + ((gten)+6)[1]*((gten)+6
)[1] + ((gten)+6)[2]*((gten)+6)[2])));
281 }
282 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclShapeTrace)(pvl->query[gageSclShapeTrace/8] & (1 << (gageSclShapeTrace
% 8)))) {
283 pvl->directAnswer[gageSclShapeTrace][0] = (curv
284 ? ELL_3M_TRACE(gten)((gten)[0] + (gten)[4] + (gten)[8])/curv
285 : 0);
286 }
287 if ( (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclK1)(pvl->query[gageSclK1/8] & (1 << (gageSclK1 % 8)
))) ||
288 (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclK2)(pvl->query[gageSclK2/8] & (1 << (gageSclK2 % 8)
))) ){
289 double T, N, D;
290 T = ELL_3M_TRACE(gten)((gten)[0] + (gten)[4] + (gten)[8]);
291 N = curv;
292 D = 2*N*N - T*T;
293 /*
294 if (D < -0.0000001) {
295 fprintf(stderr, "%s: %g %g\n", me, T, N);
296 fprintf(stderr, "%s: !!! D curv determinant % 22.10f < 0.0\n", me, D);
297 fprintf(stderr, "%s: gten: \n", me);
298 ell_3m_print_d(stderr, gten);
299 }
300 */
301 D = AIR_MAX(D, 0)((D) > (0) ? (D) : (0));
302 D = sqrt(D);
303 k1[0] = 0.5*(T + D);
304 k2[0] = 0.5*(T - D);
305 }
306 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclMeanCurv)(pvl->query[gageSclMeanCurv/8] & (1 << (gageSclMeanCurv
% 8)))) {
307 pvl->directAnswer[gageSclMeanCurv][0] = (*k1 + *k2)/2;
308 }
309 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclGaussCurv)(pvl->query[gageSclGaussCurv/8] & (1 << (gageSclGaussCurv
% 8)))) {
310 pvl->directAnswer[gageSclGaussCurv][0] = (*k1)*(*k2);
311 }
312 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclShapeIndex)(pvl->query[gageSclShapeIndex/8] & (1 << (gageSclShapeIndex
% 8)))) {
313 pvl->directAnswer[gageSclShapeIndex][0] =
314 -(2/AIR_PI3.14159265358979323846)*atan2(*k1 + *k2, *k1 - *k2);
315 }
316 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclCurvDir1)(pvl->query[gageSclCurvDir1/8] & (1 << (gageSclCurvDir1
% 8)))) {
317 /* HEY: this only works when K1, K2, 0 are all well mutually distinct,
318 since these are the eigenvalues of the geometry tensor, and this
319 code assumes that the eigenspaces are all one-dimensional */
320 ELL_3M_COPY(tmpMat, gten)((((tmpMat)+0)[0] = ((gten)+0)[0], ((tmpMat)+0)[1] = ((gten)+
0)[1], ((tmpMat)+0)[2] = ((gten)+0)[2]), (((tmpMat)+3)[0] = (
(gten)+3)[0], ((tmpMat)+3)[1] = ((gten)+3)[1], ((tmpMat)+3)[2
] = ((gten)+3)[2]), (((tmpMat)+6)[0] = ((gten)+6)[0], ((tmpMat
)+6)[1] = ((gten)+6)[1], ((tmpMat)+6)[2] = ((gten)+6)[2]));
321 ELL_3M_DIAG_SET(tmpMat, gten[0] - *k1, gten[4]- *k1, gten[8] - *k1)((tmpMat)[0] = (gten[0] - *k1), (tmpMat)[4] = (gten[4]- *k1),
(tmpMat)[8] = (gten[8] - *k1));
322 ell_3m_1d_nullspace_d(tmpVec, tmpMat);
323 ELL_3V_COPY(pvl->directAnswer[gageSclCurvDir1], tmpVec)((pvl->directAnswer[gageSclCurvDir1])[0] = (tmpVec)[0], (pvl
->directAnswer[gageSclCurvDir1])[1] = (tmpVec)[1], (pvl->
directAnswer[gageSclCurvDir1])[2] = (tmpVec)[2]);
324 }
325 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclCurvDir2)(pvl->query[gageSclCurvDir2/8] & (1 << (gageSclCurvDir2
% 8)))) {
326 /* HEY: this only works when K1, K2, 0 are all well mutually distinct,
327 since these are the eigenvalues of the geometry tensor, and this
328 code assumes that the eigenspaces are all one-dimensional */
329 ELL_3M_COPY(tmpMat, gten)((((tmpMat)+0)[0] = ((gten)+0)[0], ((tmpMat)+0)[1] = ((gten)+
0)[1], ((tmpMat)+0)[2] = ((gten)+0)[2]), (((tmpMat)+3)[0] = (
(gten)+3)[0], ((tmpMat)+3)[1] = ((gten)+3)[1], ((tmpMat)+3)[2
] = ((gten)+3)[2]), (((tmpMat)+6)[0] = ((gten)+6)[0], ((tmpMat
)+6)[1] = ((gten)+6)[1], ((tmpMat)+6)[2] = ((gten)+6)[2]));
330 ELL_3M_DIAG_SET(tmpMat, gten[0] - *k2, gten[4] - *k2, gten[8] - *k2)((tmpMat)[0] = (gten[0] - *k2), (tmpMat)[4] = (gten[4] - *k2)
, (tmpMat)[8] = (gten[8] - *k2));
331 ell_3m_1d_nullspace_d(tmpVec, tmpMat);
332 ELL_3V_COPY(pvl->directAnswer[gageSclCurvDir2], tmpVec)((pvl->directAnswer[gageSclCurvDir2])[0] = (tmpVec)[0], (pvl
->directAnswer[gageSclCurvDir2])[1] = (tmpVec)[1], (pvl->
directAnswer[gageSclCurvDir2])[2] = (tmpVec)[2]);
333 }
334 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclFlowlineCurv)(pvl->query[gageSclFlowlineCurv/8] & (1 << (gageSclFlowlineCurv
% 8)))) {
335 if (gmag >= ctx->parm.gradMagCurvMin) {
336 /* because of the gageSclGeomTens prerequisite, sHess, nPerp, and
337 nProj are all already set */
338 /* ncTen = nPerp * sHess * nProj */
339 ELL_3M_MUL(tmpMat, sHess, nProj)((tmpMat)[0] = (sHess)[0]*(nProj)[0] + (sHess)[1]*(nProj)[3] +
(sHess)[2]*(nProj)[6], (tmpMat)[1] = (sHess)[0]*(nProj)[1] +
(sHess)[1]*(nProj)[4] + (sHess)[2]*(nProj)[7], (tmpMat)[2] =
(sHess)[0]*(nProj)[2] + (sHess)[1]*(nProj)[5] + (sHess)[2]*(
nProj)[8], (tmpMat)[3] = (sHess)[3]*(nProj)[0] + (sHess)[4]*(
nProj)[3] + (sHess)[5]*(nProj)[6], (tmpMat)[4] = (sHess)[3]*(
nProj)[1] + (sHess)[4]*(nProj)[4] + (sHess)[5]*(nProj)[7], (tmpMat
)[5] = (sHess)[3]*(nProj)[2] + (sHess)[4]*(nProj)[5] + (sHess
)[5]*(nProj)[8], (tmpMat)[6] = (sHess)[6]*(nProj)[0] + (sHess
)[7]*(nProj)[3] + (sHess)[8]*(nProj)[6], (tmpMat)[7] = (sHess
)[6]*(nProj)[1] + (sHess)[7]*(nProj)[4] + (sHess)[8]*(nProj)[
7], (tmpMat)[8] = (sHess)[6]*(nProj)[2] + (sHess)[7]*(nProj)[
5] + (sHess)[8]*(nProj)[8]);
340 ELL_3M_MUL(ncTen, nPerp, tmpMat)((ncTen)[0] = (nPerp)[0]*(tmpMat)[0] + (nPerp)[1]*(tmpMat)[3]
+ (nPerp)[2]*(tmpMat)[6], (ncTen)[1] = (nPerp)[0]*(tmpMat)[1
] + (nPerp)[1]*(tmpMat)[4] + (nPerp)[2]*(tmpMat)[7], (ncTen)[
2] = (nPerp)[0]*(tmpMat)[2] + (nPerp)[1]*(tmpMat)[5] + (nPerp
)[2]*(tmpMat)[8], (ncTen)[3] = (nPerp)[3]*(tmpMat)[0] + (nPerp
)[4]*(tmpMat)[3] + (nPerp)[5]*(tmpMat)[6], (ncTen)[4] = (nPerp
)[3]*(tmpMat)[1] + (nPerp)[4]*(tmpMat)[4] + (nPerp)[5]*(tmpMat
)[7], (ncTen)[5] = (nPerp)[3]*(tmpMat)[2] + (nPerp)[4]*(tmpMat
)[5] + (nPerp)[5]*(tmpMat)[8], (ncTen)[6] = (nPerp)[6]*(tmpMat
)[0] + (nPerp)[7]*(tmpMat)[3] + (nPerp)[8]*(tmpMat)[6], (ncTen
)[7] = (nPerp)[6]*(tmpMat)[1] + (nPerp)[7]*(tmpMat)[4] + (nPerp
)[8]*(tmpMat)[7], (ncTen)[8] = (nPerp)[6]*(tmpMat)[2] + (nPerp
)[7]*(tmpMat)[5] + (nPerp)[8]*(tmpMat)[8]);
341 } else {
342 ELL_3M_ZERO_SET(ncTen)((((ncTen)+0)[0] = (0), ((ncTen)+0)[1] = (0), ((ncTen)+0)[2] =
(0)), (((ncTen)+3)[0] = (0), ((ncTen)+3)[1] = (0), ((ncTen)+
3)[2] = (0)), (((ncTen)+6)[0] = (0), ((ncTen)+6)[1] = (0), ((
ncTen)+6)[2] = (0)));
343 }
344 /* there used to be a wrong extra sqrt() here */
345 pvl->directAnswer[gageSclFlowlineCurv][0] = ELL_3M_FROB(ncTen)(sqrt((((ncTen)+0)[0]*((ncTen)+0)[0] + ((ncTen)+0)[1]*((ncTen
)+0)[1] + ((ncTen)+0)[2]*((ncTen)+0)[2]) + (((ncTen)+3)[0]*((
ncTen)+3)[0] + ((ncTen)+3)[1]*((ncTen)+3)[1] + ((ncTen)+3)[2]
*((ncTen)+3)[2]) + (((ncTen)+6)[0]*((ncTen)+6)[0] + ((ncTen)+
6)[1]*((ncTen)+6)[1] + ((ncTen)+6)[2]*((ncTen)+6)[2])));
346 }
347 if (GAGE_QUERY_ITEM_TEST(pvl->query, gageSclMedian)(pvl->query[gageSclMedian/8] & (1 << (gageSclMedian
% 8)))) {
348 /* this item is currently a complete oddball in that it does not
349 benefit from anything done in the "filter" stage, which is in
350 fact a waste of time if the query consists only of this item */
351 fd = 2*ctx->radius;
352 if (fd > FD_MEDIAN_MAX16) {
353 fprintf(stderr__stderrp, "%s: PANIC: current filter diameter = %d "
354 "> FD_MEDIAN_MAX = %d\n", me, fd, FD_MEDIAN_MAX16);
355 exit(1);
356 }
357 fw = ctx->fw + fd*3*gageKernel00;
358 /* HEY: this needs some optimization help */
359 wghtSum = 0;
360 nidx = 0;
361 for (xi=0; xi<fd; xi++) {
362 for (yi=0; yi<fd; yi++) {
363 for (zi=0; zi<fd; zi++) {
364 iv3wght[0 + 2*nidx] = pvl->iv3[nidx];
365 iv3wght[1 + 2*nidx] = fw[xi + 0*fd]*fw[yi + 1*fd]*fw[zi + 2*fd];
366 wghtSum += iv3wght[1 + 2*nidx];
367 nidx++;
368 }
369 }
370 }
371 qsort(iv3wght, fd*fd*fd, 2*sizeof(double), nrrdValCompare[nrrdTypeDouble]);
372 wght = 0;
373 for (nidx=0; nidx<fd*fd*fd; nidx++) {
374 wght += iv3wght[1 + 2*nidx];
375 if (wght > wghtSum/2) {
376 break;
377 }
378 }
379 pvl->directAnswer[gageSclMedian][0] = iv3wght[0 + 2*nidx];
380 }
381 return;
382}
383
384#undef TEN_M2T
385#undef TEN_T_SCALE