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 |
|
|
|
27 |
|
|
typedef struct { |
28 |
|
|
double *buffTen, *buffWght; |
29 |
|
|
tenInterpParm *tip; /* sneakiness: using tip->allocLen to record |
30 |
|
|
allocation sizes of buffTen and buffWght, too */ |
31 |
|
|
} _tenGagePvlData; |
32 |
|
|
|
33 |
|
|
gageItemEntry |
34 |
|
|
_tenGageTable[TEN_GAGE_ITEM_MAX+1] = { |
35 |
|
|
/* enum value len,deriv, prereqs, parent item, parent index, needData */ |
36 |
|
|
{tenGageUnknown, 0, 0, {0}, 0, 0, AIR_FALSE}, |
37 |
|
|
{tenGageTensor, 7, 0, {0}, 0, 0, AIR_FALSE}, |
38 |
|
|
{tenGageConfidence, 1, 0, {tenGageTensor}, tenGageTensor, 0, AIR_FALSE}, |
39 |
|
|
|
40 |
|
|
{tenGageTrace, 1, 0, {tenGageTensor}, 0, 0, AIR_FALSE}, |
41 |
|
|
{tenGageNorm, 1, 0, {tenGageTensor}, 0, 0, AIR_FALSE}, |
42 |
|
|
{tenGageB, 1, 0, {tenGageTensor}, 0, 0, AIR_FALSE}, |
43 |
|
|
{tenGageDet, 1, 0, {tenGageTensor}, 0, 0, AIR_FALSE}, |
44 |
|
|
{tenGageS, 1, 0, {tenGageTensor}, 0, 0, AIR_FALSE}, |
45 |
|
|
{tenGageQ, 1, 0, {tenGageS, tenGageB}, 0, 0, AIR_FALSE}, |
46 |
|
|
{tenGageFA, 1, 0, {tenGageQ, tenGageS}, 0, 0, AIR_FALSE}, |
47 |
|
|
{tenGageR, 1, 0, {tenGageTrace, tenGageB, tenGageDet, tenGageS}, 0, 0, AIR_FALSE}, |
48 |
|
|
{tenGageMode, 1, 0, {tenGageR, tenGageQ}, 0, 0, AIR_FALSE}, |
49 |
|
|
{tenGageTheta, 1, 0, {tenGageMode}, 0, 0, AIR_FALSE}, |
50 |
|
|
{tenGageModeWarp, 1, 0, {tenGageMode}, 0, 0, AIR_FALSE}, |
51 |
|
|
{tenGageOmega, 1, 0, {tenGageFA, tenGageMode}, 0, 0, AIR_FALSE}, |
52 |
|
|
|
53 |
|
|
{tenGageEval, 3, 0, {tenGageTensor}, 0, 0, AIR_FALSE}, |
54 |
|
|
{tenGageEval0, 1, 0, {tenGageEval}, tenGageEval, 0, AIR_FALSE}, |
55 |
|
|
{tenGageEval1, 1, 0, {tenGageEval}, tenGageEval, 1, AIR_FALSE}, |
56 |
|
|
{tenGageEval2, 1, 0, {tenGageEval}, tenGageEval, 2, AIR_FALSE}, |
57 |
|
|
{tenGageEvec, 9, 0, {tenGageTensor}, 0, 0, AIR_FALSE}, |
58 |
|
|
{tenGageEvec0, 3, 0, {tenGageEvec}, tenGageEvec, 0, AIR_FALSE}, |
59 |
|
|
{tenGageEvec1, 3, 0, {tenGageEvec}, tenGageEvec, 3, AIR_FALSE}, |
60 |
|
|
{tenGageEvec2, 3, 0, {tenGageEvec}, tenGageEvec, 6, AIR_FALSE}, |
61 |
|
|
|
62 |
|
|
{tenGageDelNormK2, 7, 0, {tenGageTensor}, 0, 0, AIR_FALSE}, |
63 |
|
|
{tenGageDelNormK3, 7, 0, {tenGageTensor}, 0, 0, AIR_FALSE}, |
64 |
|
|
{tenGageDelNormR1, 7, 0, {tenGageTensor}, 0, 0, AIR_FALSE}, |
65 |
|
|
{tenGageDelNormR2, 7, 0, {tenGageTensor}, 0, 0, AIR_FALSE}, |
66 |
|
|
|
67 |
|
|
{tenGageDelNormPhi1, 7, 0, {tenGageEvec}, 0, 0, AIR_FALSE}, |
68 |
|
|
{tenGageDelNormPhi2, 7, 0, {tenGageEvec}, 0, 0, AIR_FALSE}, |
69 |
|
|
{tenGageDelNormPhi3, 7, 0, {tenGageEvec}, 0, 0, AIR_FALSE}, |
70 |
|
|
|
71 |
|
|
{tenGageTensorGrad, 21, 1, {0}, 0, 0, AIR_FALSE}, |
72 |
|
|
{tenGageTensorGradMag, 3, 1, {tenGageTensorGrad}, 0, 0, AIR_FALSE}, |
73 |
|
|
{tenGageTensorGradMagMag, 1, 1, {tenGageTensorGradMag}, 0, 0, AIR_FALSE}, |
74 |
|
|
|
75 |
|
|
{tenGageTraceGradVec, 3, 1, {tenGageTensor, tenGageTensorGrad}, 0, 0, AIR_FALSE}, |
76 |
|
|
{tenGageTraceGradMag, 1, 1, {tenGageTraceGradVec}, 0, 0, AIR_FALSE}, |
77 |
|
|
{tenGageTraceNormal, 3, 1, {tenGageTraceGradVec, tenGageTraceGradMag}, 0, 0, AIR_FALSE}, |
78 |
|
|
|
79 |
|
|
{tenGageNormGradVec, 3, 1, {tenGageNorm, tenGageSGradVec}, 0, 0, AIR_FALSE}, |
80 |
|
|
{tenGageNormGradMag, 1, 1, {tenGageNormGradVec}, 0, 0, AIR_FALSE}, |
81 |
|
|
{tenGageNormNormal, 3, 1, {tenGageNormGradVec, tenGageNormGradMag}, 0, 0, AIR_FALSE}, |
82 |
|
|
|
83 |
|
|
{tenGageBGradVec, 3, 1, {tenGageTensor, tenGageTensorGrad}, 0, 0, AIR_FALSE}, |
84 |
|
|
{tenGageBGradMag, 1, 1, {tenGageBGradVec}, 0, 0, AIR_FALSE}, |
85 |
|
|
{tenGageBNormal, 3, 1, {tenGageBGradVec, tenGageBGradMag}, 0, 0, AIR_FALSE}, |
86 |
|
|
|
87 |
|
|
{tenGageDetGradVec, 3, 1, {tenGageTensor, tenGageTensorGrad}, 0, 0, AIR_FALSE}, |
88 |
|
|
{tenGageDetGradMag, 1, 1, {tenGageDetGradVec}, 0, 0, AIR_FALSE}, |
89 |
|
|
{tenGageDetNormal, 3, 1, {tenGageDetGradVec, tenGageDetGradMag}, 0, 0, AIR_FALSE}, |
90 |
|
|
|
91 |
|
|
{tenGageSGradVec, 3, 1, {tenGageTensor, tenGageTensorGrad}, 0, 0, AIR_FALSE}, |
92 |
|
|
{tenGageSGradMag, 1, 1, {tenGageSGradVec}, 0, 0, AIR_FALSE}, |
93 |
|
|
{tenGageSNormal, 3, 1, {tenGageSGradVec, tenGageSGradMag}, 0, 0, AIR_FALSE}, |
94 |
|
|
|
95 |
|
|
{tenGageQGradVec, 3, 1, {tenGageSGradVec, tenGageBGradVec}, 0, 0, AIR_FALSE}, |
96 |
|
|
{tenGageQGradMag, 1, 1, {tenGageQGradVec}, 0, 0, AIR_FALSE}, |
97 |
|
|
{tenGageQNormal, 3, 1, {tenGageQGradVec, tenGageQGradMag}, 0, 0, AIR_FALSE}, |
98 |
|
|
|
99 |
|
|
{tenGageFAGradVec, 3, 1, {tenGageQGradVec, tenGageSGradVec, tenGageFA}, 0, 0, AIR_FALSE}, |
100 |
|
|
{tenGageFAGradMag, 1, 1, {tenGageFAGradVec}, 0, 0, AIR_FALSE}, |
101 |
|
|
{tenGageFANormal, 3, 1, {tenGageFAGradVec, tenGageFAGradMag}, 0, 0, AIR_FALSE}, |
102 |
|
|
|
103 |
|
|
{tenGageRGradVec, 3, 1, {tenGageR, tenGageTraceGradVec, tenGageBGradVec, |
104 |
|
|
tenGageDetGradVec, tenGageSGradVec}, 0, 0, AIR_FALSE}, |
105 |
|
|
{tenGageRGradMag, 1, 1, {tenGageRGradVec}, 0, 0, AIR_FALSE}, |
106 |
|
|
{tenGageRNormal, 3, 1, {tenGageRGradVec, tenGageRGradMag}, 0, 0, AIR_FALSE}, |
107 |
|
|
|
108 |
|
|
{tenGageModeGradVec, 3, 1, {tenGageRGradVec, tenGageQGradVec, tenGageMode}, 0, 0, AIR_FALSE}, |
109 |
|
|
{tenGageModeGradMag, 1, 1, {tenGageModeGradVec}, 0, 0, AIR_FALSE}, |
110 |
|
|
{tenGageModeNormal, 3, 1, {tenGageModeGradVec, tenGageModeGradMag}, 0, 0, AIR_FALSE}, |
111 |
|
|
|
112 |
|
|
{tenGageThetaGradVec, 3, 1, {tenGageRGradVec, tenGageQGradVec, tenGageTheta}, 0, 0, AIR_FALSE}, |
113 |
|
|
{tenGageThetaGradMag, 1, 1, {tenGageThetaGradVec}, 0, 0, AIR_FALSE}, |
114 |
|
|
{tenGageThetaNormal, 3, 1, {tenGageThetaGradVec, tenGageThetaGradMag}, 0, 0, AIR_FALSE}, |
115 |
|
|
|
116 |
|
|
{tenGageOmegaGradVec, 3, 1, {tenGageFA, tenGageMode, tenGageFAGradVec, tenGageModeGradVec}, 0, 0, AIR_FALSE}, |
117 |
|
|
{tenGageOmegaGradMag, 1, 1, {tenGageOmegaGradVec}, 0, 0, AIR_FALSE}, |
118 |
|
|
{tenGageOmegaNormal, 3, 1, {tenGageOmegaGradVec, tenGageOmegaGradMag}, 0, 0, AIR_FALSE}, |
119 |
|
|
|
120 |
|
|
{tenGageInvarKGrads, 9, 1, {tenGageDelNormK2, tenGageDelNormK3, tenGageTensorGrad}, 0, 0, AIR_FALSE}, |
121 |
|
|
{tenGageInvarKGradMags, 3, 1, {tenGageInvarKGrads}, 0, 0, AIR_FALSE}, |
122 |
|
|
{tenGageInvarRGrads, 9, 1, {tenGageDelNormR1, tenGageDelNormR2, tenGageDelNormK3, |
123 |
|
|
tenGageTensorGrad}, 0, 0, AIR_FALSE}, |
124 |
|
|
{tenGageInvarRGradMags, 3, 1, {tenGageInvarRGrads}, 0, 0, AIR_FALSE}, |
125 |
|
|
|
126 |
|
|
{tenGageRotTans, 9, 1, {tenGageDelNormPhi1, tenGageDelNormPhi2, tenGageDelNormPhi3, |
127 |
|
|
tenGageTensorGrad}, 0, 0, AIR_FALSE}, |
128 |
|
|
{tenGageRotTanMags, 3, 1, {tenGageRotTans}, 0, 0, AIR_FALSE}, |
129 |
|
|
|
130 |
|
|
{tenGageEvalGrads, 9, 1, {tenGageTensorGrad, tenGageEval, tenGageEvec}, 0, 0, AIR_FALSE}, |
131 |
|
|
|
132 |
|
|
{tenGageCl1, 1, 0, {tenGageTensor, tenGageEval0, tenGageEval1, tenGageEval2}, 0, 0, AIR_FALSE}, |
133 |
|
|
{tenGageCp1, 1, 0, {tenGageTensor, tenGageEval0, tenGageEval1, tenGageEval2}, 0, 0, AIR_FALSE}, |
134 |
|
|
{tenGageCa1, 1, 0, {tenGageTensor, tenGageEval0, tenGageEval1, tenGageEval2}, 0, 0, AIR_FALSE}, |
135 |
|
|
{tenGageClpmin1, 1, 0, {tenGageTensor, tenGageEval0, tenGageEval1, tenGageEval2}, 0, 0, AIR_FALSE}, |
136 |
|
|
{tenGageCl2, 1, 0, {tenGageTensor, tenGageEval0, tenGageEval1, tenGageEval2}, 0, 0, AIR_FALSE}, |
137 |
|
|
{tenGageCp2, 1, 0, {tenGageTensor, tenGageEval0, tenGageEval1, tenGageEval2}, 0, 0, AIR_FALSE}, |
138 |
|
|
{tenGageCa2, 1, 0, {tenGageTensor, tenGageEval0, tenGageEval1, tenGageEval2}, 0, 0, AIR_FALSE}, |
139 |
|
|
{tenGageClpmin2, 1, 0, {tenGageTensor, tenGageEval0, tenGageEval1, tenGageEval2}, 0, 0, AIR_FALSE}, |
140 |
|
|
|
141 |
|
|
{tenGageHessian, 63, 2, {0}, 0, 0, AIR_FALSE}, |
142 |
|
|
{tenGageTraceHessian, 9, 2, {tenGageHessian}, 0, 0, AIR_FALSE}, |
143 |
|
|
{tenGageTraceHessianEval, 3, 2, {tenGageTraceHessian}, 0, 0, AIR_FALSE}, |
144 |
|
|
{tenGageTraceHessianEval0, 1, 2, {tenGageTraceHessianEval}, tenGageTraceHessianEval, 0, AIR_FALSE}, |
145 |
|
|
{tenGageTraceHessianEval1, 1, 2, {tenGageTraceHessianEval}, tenGageTraceHessianEval, 1, AIR_FALSE}, |
146 |
|
|
{tenGageTraceHessianEval2, 1, 2, {tenGageTraceHessianEval}, tenGageTraceHessianEval, 2, AIR_FALSE}, |
147 |
|
|
{tenGageTraceHessianEvec, 9, 2, {tenGageTraceHessian}, 0, 0, AIR_FALSE}, |
148 |
|
|
{tenGageTraceHessianEvec0, 3, 2, {tenGageTraceHessianEvec}, tenGageTraceHessianEvec, 0, AIR_FALSE}, |
149 |
|
|
{tenGageTraceHessianEvec1, 3, 2, {tenGageTraceHessianEvec}, tenGageTraceHessianEvec, 3, AIR_FALSE}, |
150 |
|
|
{tenGageTraceHessianEvec2, 3, 2, {tenGageTraceHessianEvec}, tenGageTraceHessianEvec, 6, AIR_FALSE}, |
151 |
|
|
{tenGageTraceHessianFrob, 1, 2, {tenGageTraceHessian}, 0, 0, AIR_FALSE}, |
152 |
|
|
{tenGageBHessian, 9, 2, {tenGageTensor, tenGageTensorGrad, tenGageHessian}, 0, 0, AIR_FALSE}, |
153 |
|
|
{tenGageDetHessian, 9, 2, {tenGageTensor, tenGageTensorGrad, tenGageHessian}, 0, 0, AIR_FALSE}, |
154 |
|
|
{tenGageSHessian, 9, 2, {tenGageTensor, tenGageTensorGrad, tenGageHessian}, 0, 0, AIR_FALSE}, |
155 |
|
|
{tenGageQHessian, 9, 2, {tenGageBHessian, tenGageSHessian}, 0, 0, AIR_FALSE}, |
156 |
|
|
|
157 |
|
|
{tenGageFAHessian, 9, 2, {tenGageSHessian, tenGageQHessian, |
158 |
|
|
tenGageSGradVec, tenGageQGradVec, tenGageFA}, 0, 0, AIR_FALSE}, |
159 |
|
|
{tenGageFAHessianEval, 3, 2, {tenGageFAHessian}, 0, 0, AIR_FALSE}, |
160 |
|
|
{tenGageFAHessianEval0, 1, 2, {tenGageFAHessianEval}, tenGageFAHessianEval, 0, AIR_FALSE}, |
161 |
|
|
{tenGageFAHessianEval1, 1, 2, {tenGageFAHessianEval}, tenGageFAHessianEval, 1, AIR_FALSE}, |
162 |
|
|
{tenGageFAHessianEval2, 1, 2, {tenGageFAHessianEval}, tenGageFAHessianEval, 2, AIR_FALSE}, |
163 |
|
|
{tenGageFAHessianEvec, 9, 2, {tenGageFAHessian}, 0, 0, AIR_FALSE}, |
164 |
|
|
{tenGageFAHessianEvec0, 3, 2, {tenGageFAHessianEvec}, tenGageFAHessianEvec, 0, AIR_FALSE}, |
165 |
|
|
{tenGageFAHessianEvec1, 3, 2, {tenGageFAHessianEvec}, tenGageFAHessianEvec, 3, AIR_FALSE}, |
166 |
|
|
{tenGageFAHessianEvec2, 3, 2, {tenGageFAHessianEvec}, tenGageFAHessianEvec, 6, AIR_FALSE}, |
167 |
|
|
{tenGageFAHessianFrob, 1, 2, {tenGageFAHessian}, 0, 0, AIR_FALSE}, |
168 |
|
|
{tenGageFARidgeSurfaceStrength, 1, 2, {tenGageConfidence, tenGageFAHessianEval}, 0, 0, AIR_FALSE}, |
169 |
|
|
{tenGageFAValleySurfaceStrength, 1, 2, {tenGageConfidence, tenGageFAHessianEval}, 0, 0, AIR_FALSE}, |
170 |
|
|
{tenGageFALaplacian, 1, 2, {tenGageFAHessian}, 0, 0, AIR_FALSE}, |
171 |
|
|
{tenGageFAHessianEvalMode, 1, 2, {tenGageFAHessianEval}, 0, 0, AIR_FALSE}, |
172 |
|
|
{tenGageFARidgeLineAlignment, 1, 2, {tenGageEvec0, tenGageFAHessianEvec0, tenGageFAHessianEvalMode}, 0, 0, AIR_FALSE}, |
173 |
|
|
{tenGageFARidgeSurfaceAlignment, 1, 2, {tenGageEvec0, tenGageFAHessianEvec2, tenGageFAHessianEvalMode}, 0, 0, AIR_FALSE}, |
174 |
|
|
{tenGageFA2ndDD, 1, 2, {tenGageFAHessian, tenGageFANormal}, 0, 0, AIR_FALSE}, |
175 |
|
|
|
176 |
|
|
{tenGageFAGeomTens, 9, 2, {tenGageFAHessian, tenGageFAGradMag, tenGageFANormal}, 0, 0, AIR_FALSE}, |
177 |
|
|
{tenGageFAKappa1, 1, 2, {tenGageFAGeomTens}, 0, 0, AIR_FALSE}, |
178 |
|
|
{tenGageFAKappa2, 1, 2, {tenGageFAGeomTens}, 0, 0, AIR_FALSE}, |
179 |
|
|
{tenGageFATotalCurv, 1, 2, {tenGageFAGeomTens}, 0, 0, AIR_FALSE}, |
180 |
|
|
{tenGageFAShapeIndex, 1, 2, {tenGageFAKappa1, tenGageFAKappa2}, 0, 0, AIR_FALSE}, |
181 |
|
|
{tenGageFAMeanCurv, 1, 2, {tenGageFAKappa1, tenGageFAKappa2}, 0, 0, AIR_FALSE}, |
182 |
|
|
{tenGageFAGaussCurv, 1, 2, {tenGageFAKappa1, tenGageFAKappa2}, 0, 0, AIR_FALSE}, |
183 |
|
|
{tenGageFACurvDir1, 3, 2, {tenGageFAGeomTens, tenGageFAKappa2}, 0, 0, AIR_FALSE}, |
184 |
|
|
{tenGageFACurvDir2, 3, 2, {tenGageFAGeomTens, tenGageFAKappa1}, 0, 0, AIR_FALSE}, |
185 |
|
|
{tenGageFAFlowlineCurv, 1, 2, {tenGageFANormal, tenGageFAHessian, tenGageFAGradMag}, 0, 0, AIR_FALSE}, |
186 |
|
|
|
187 |
|
|
{tenGageRHessian, 9, 2, {tenGageR, tenGageRGradVec, tenGageTraceHessian, |
188 |
|
|
tenGageBHessian, tenGageDetHessian, tenGageSHessian}, 0, 0, AIR_FALSE}, |
189 |
|
|
|
190 |
|
|
{tenGageModeHessian, 9, 2, {tenGageR, tenGageQ, tenGageRGradVec, tenGageQGradVec, |
191 |
|
|
tenGageRHessian, tenGageQHessian}, 0, 0, AIR_FALSE}, |
192 |
|
|
{tenGageModeHessianEval, 3, 2, {tenGageModeHessian}, 0, 0, AIR_FALSE}, |
193 |
|
|
{tenGageModeHessianEval0, 1, 2, {tenGageModeHessianEval}, tenGageModeHessianEval, 0, AIR_FALSE}, |
194 |
|
|
{tenGageModeHessianEval1, 1, 2, {tenGageModeHessianEval}, tenGageModeHessianEval, 1, AIR_FALSE}, |
195 |
|
|
{tenGageModeHessianEval2, 1, 2, {tenGageModeHessianEval}, tenGageModeHessianEval, 2, AIR_FALSE}, |
196 |
|
|
{tenGageModeHessianEvec, 9, 2, {tenGageModeHessian}, 0, 0, AIR_FALSE}, |
197 |
|
|
{tenGageModeHessianEvec0, 3, 2, {tenGageModeHessianEvec}, tenGageModeHessianEvec, 0, AIR_FALSE}, |
198 |
|
|
{tenGageModeHessianEvec1, 3, 2, {tenGageModeHessianEvec}, tenGageModeHessianEvec, 3, AIR_FALSE}, |
199 |
|
|
{tenGageModeHessianEvec2, 3, 2, {tenGageModeHessianEvec}, tenGageModeHessianEvec, 6, AIR_FALSE}, |
200 |
|
|
{tenGageModeHessianFrob, 1, 2, {tenGageModeHessian}, 0, 0, AIR_FALSE}, |
201 |
|
|
|
202 |
|
|
{tenGageOmegaHessian, 9, 2, {tenGageFA, tenGageMode, tenGageFAGradVec, tenGageModeGradVec, |
203 |
|
|
tenGageFAHessian, tenGageModeHessian}, 0, 0, AIR_FALSE}, |
204 |
|
|
{tenGageOmegaHessianEval, 3, 2, {tenGageOmegaHessian}, 0, 0, AIR_FALSE}, |
205 |
|
|
{tenGageOmegaHessianEval0, 1, 2, {tenGageOmegaHessianEval}, tenGageOmegaHessianEval, 0, AIR_FALSE}, |
206 |
|
|
{tenGageOmegaHessianEval1, 1, 2, {tenGageOmegaHessianEval}, tenGageOmegaHessianEval, 1, AIR_FALSE}, |
207 |
|
|
{tenGageOmegaHessianEval2, 1, 2, {tenGageOmegaHessianEval}, tenGageOmegaHessianEval, 2, AIR_FALSE}, |
208 |
|
|
{tenGageOmegaHessianEvec, 9, 2, {tenGageOmegaHessian}, 0, 0, AIR_FALSE}, |
209 |
|
|
{tenGageOmegaHessianEvec0, 3, 2, {tenGageOmegaHessianEvec}, tenGageOmegaHessianEvec, 0, AIR_FALSE}, |
210 |
|
|
{tenGageOmegaHessianEvec1, 3, 2, {tenGageOmegaHessianEvec}, tenGageOmegaHessianEvec, 3, AIR_FALSE}, |
211 |
|
|
{tenGageOmegaHessianEvec2, 3, 2, {tenGageOmegaHessianEvec}, tenGageOmegaHessianEvec, 6, AIR_FALSE}, |
212 |
|
|
{tenGageOmegaLaplacian, 1, 2, {tenGageOmegaHessian}, 0, 0, AIR_FALSE}, |
213 |
|
|
{tenGageOmega2ndDD, 1, 2, {tenGageOmegaHessian, tenGageOmegaNormal}, 0, 0, AIR_FALSE}, |
214 |
|
|
|
215 |
|
|
{tenGageOmegaHessianContrTenEvec0, 1, 2, {tenGageOmegaHessian, tenGageEvec0}, 0, 0, AIR_FALSE}, |
216 |
|
|
{tenGageOmegaHessianContrTenEvec1, 1, 2, {tenGageOmegaHessian, tenGageEvec1}, 0, 0, AIR_FALSE}, |
217 |
|
|
{tenGageOmegaHessianContrTenEvec2, 1, 2, {tenGageOmegaHessian, tenGageEvec2}, 0, 0, AIR_FALSE}, |
218 |
|
|
|
219 |
|
|
{tenGageTraceGradVecDotEvec0, 1, 1, {tenGageTraceGradVec, tenGageEvec0}, 0, 0, AIR_FALSE}, |
220 |
|
|
{tenGageTraceDiffusionAlign, 1, 1, {tenGageTraceNormal, tenGageEvec0}, 0, 0, AIR_FALSE}, |
221 |
|
|
{tenGageTraceDiffusionFraction, 1, 1, {tenGageTraceNormal, tenGageTensor}, 0, 0, AIR_FALSE}, |
222 |
|
|
|
223 |
|
|
{tenGageFAGradVecDotEvec0, 1, 1, {tenGageFAGradVec, tenGageEvec0}, 0, 0, AIR_FALSE}, |
224 |
|
|
{tenGageFADiffusionAlign, 1, 1, {tenGageFANormal, tenGageEvec0}, 0, 0, AIR_FALSE}, |
225 |
|
|
{tenGageFADiffusionFraction, 1, 1, {tenGageFANormal, tenGageTensor}, 0, 0, AIR_FALSE}, |
226 |
|
|
|
227 |
|
|
{tenGageOmegaGradVecDotEvec0, 1, 1, {tenGageOmegaGradVec, tenGageEvec0}, 0, 0, AIR_FALSE}, |
228 |
|
|
{tenGageOmegaDiffusionAlign, 1, 1, {tenGageOmegaNormal, tenGageEvec0}, 0, 0, AIR_FALSE}, |
229 |
|
|
{tenGageOmegaDiffusionFraction, 1, 1, {tenGageOmegaNormal, tenGageTensor}, 0, 0, AIR_FALSE}, |
230 |
|
|
|
231 |
|
|
/* currently don't have tenGageConfGradVec */ |
232 |
|
|
{tenGageConfGradVecDotEvec0, 1, 1, {tenGageTensorGrad, tenGageEvec0}, 0, 0, AIR_FALSE}, |
233 |
|
|
{tenGageConfDiffusionAlign, 1, 1, {tenGageTensorGrad, tenGageEvec0}, 0, 0, AIR_FALSE}, |
234 |
|
|
{tenGageConfDiffusionFraction, 1, 1, {tenGageTensorGrad, tenGageTensor}, 0, 0, AIR_FALSE}, |
235 |
|
|
|
236 |
|
|
{tenGageCovariance, 21, 0, {tenGageTensor}, /* and all the values in iv3 */ 0, 0, AIR_FALSE}, |
237 |
|
|
{tenGageCovarianceRGRT, 21, 0, {tenGageCovariance, |
238 |
|
|
tenGageDelNormR1, tenGageDelNormR2, tenGageDelNormK3, |
239 |
|
|
tenGageDelNormPhi1, tenGageDelNormPhi2, tenGageDelNormPhi3}, 0, 0, AIR_FALSE}, |
240 |
|
|
{tenGageCovarianceKGRT, 21, 0, {tenGageCovariance, |
241 |
|
|
tenGageDelNormK2, tenGageDelNormK3, |
242 |
|
|
tenGageDelNormPhi1, tenGageDelNormPhi2, tenGageDelNormPhi3}, 0, 0, AIR_FALSE}, |
243 |
|
|
|
244 |
|
|
{tenGageTensorLogEuclidean, 7, 0, {0}, 0, 0, AIR_FALSE}, |
245 |
|
|
{tenGageTensorQuatGeoLoxK, 7, 0, {0}, 0, 0, AIR_FALSE}, |
246 |
|
|
{tenGageTensorQuatGeoLoxR, 7, 0, {0}, 0, 0, AIR_FALSE}, |
247 |
|
|
{tenGageTensorRThetaPhiLinear, 7, 0, {0}, 0, 0, AIR_FALSE}, |
248 |
|
|
{tenGageCl1GradVec, 3, 1, {tenGageTrace, tenGageEval, tenGageEvalGrads}, 0, 0, AIR_FALSE}, |
249 |
|
|
{tenGageCl1GradMag, 1, 1, {tenGageCl1GradVec}, 0, 0, AIR_FALSE}, |
250 |
|
|
{tenGageCl1Normal, 3, 1, {tenGageCl1GradVec, tenGageCl1GradMag}, 0, 0, AIR_FALSE}, |
251 |
|
|
{tenGageCp1GradVec, 3, 1, {tenGageTrace, tenGageEval, tenGageEvalGrads}, 0, 0, AIR_FALSE}, |
252 |
|
|
{tenGageCp1GradMag, 1, 1, {tenGageCp1GradVec}, 0, 0, AIR_FALSE}, |
253 |
|
|
{tenGageCp1Normal, 3, 1, {tenGageCp1GradVec, tenGageCp1GradMag}, 0, 0, AIR_FALSE}, |
254 |
|
|
{tenGageCa1GradVec, 3, 1, {tenGageTrace, tenGageEval, tenGageEvalGrads}, 0, 0, AIR_FALSE}, |
255 |
|
|
{tenGageCa1GradMag, 1, 1, {tenGageCa1GradVec}, 0, 0, AIR_FALSE}, |
256 |
|
|
{tenGageCa1Normal, 3, 1, {tenGageCa1GradVec, tenGageCa1GradMag}, 0, 0, AIR_FALSE}, |
257 |
|
|
{tenGageTensorGradRotE, 21, 1, {tenGageTensorGrad, tenGageEval, tenGageEvec}, 0, 0, AIR_FALSE}, |
258 |
|
|
{tenGageEvalHessian, 27, 2, {tenGageTensorGradRotE, tenGageHessian, tenGageEval}, 0, 0, AIR_FALSE }, |
259 |
|
|
{tenGageCl1Hessian, 9, 2, {tenGageTensorGradRotE, tenGageEvalHessian}, 0, 0, AIR_FALSE }, |
260 |
|
|
{tenGageCl1HessianEval, 3, 2, {tenGageCl1Hessian}, 0, 0, AIR_FALSE }, |
261 |
|
|
{tenGageCl1HessianEval0, 1, 2, {tenGageCl1HessianEval}, tenGageCl1HessianEval, 0, AIR_FALSE}, |
262 |
|
|
{tenGageCl1HessianEval1, 1, 2, {tenGageCl1HessianEval}, tenGageCl1HessianEval, 1, AIR_FALSE}, |
263 |
|
|
{tenGageCl1HessianEval2, 1, 2, {tenGageCl1HessianEval}, tenGageCl1HessianEval, 2, AIR_FALSE}, |
264 |
|
|
{tenGageCl1HessianEvec, 9, 2, {tenGageCl1Hessian}, 0, 0, AIR_FALSE }, |
265 |
|
|
{tenGageCl1HessianEvec0, 3, 2, {tenGageCl1HessianEvec}, tenGageCl1HessianEvec, 0, AIR_FALSE}, |
266 |
|
|
{tenGageCl1HessianEvec1, 3, 2, {tenGageCl1HessianEvec}, tenGageCl1HessianEvec, 3, AIR_FALSE}, |
267 |
|
|
{tenGageCl1HessianEvec2, 3, 2, {tenGageCl1HessianEvec}, tenGageCl1HessianEvec, 6, AIR_FALSE}, |
268 |
|
|
{tenGageCp1Hessian, 9, 2, {tenGageTensorGradRotE, tenGageEvalHessian}, 0, 0, AIR_FALSE }, |
269 |
|
|
{tenGageCp1HessianEval, 3, 2, {tenGageCp1Hessian}, 0, 0, AIR_FALSE }, |
270 |
|
|
{tenGageCp1HessianEval0, 1, 2, {tenGageCp1HessianEval}, tenGageCp1HessianEval, 0, AIR_FALSE}, |
271 |
|
|
{tenGageCp1HessianEval1, 1, 2, {tenGageCp1HessianEval}, tenGageCp1HessianEval, 1, AIR_FALSE}, |
272 |
|
|
{tenGageCp1HessianEval2, 1, 2, {tenGageCp1HessianEval}, tenGageCp1HessianEval, 2, AIR_FALSE}, |
273 |
|
|
{tenGageCp1HessianEvec, 9, 2, {tenGageCp1Hessian}, 0, 0, AIR_FALSE }, |
274 |
|
|
{tenGageCp1HessianEvec0, 3, 2, {tenGageCp1HessianEvec}, tenGageCp1HessianEvec, 0, AIR_FALSE}, |
275 |
|
|
{tenGageCp1HessianEvec1, 3, 2, {tenGageCp1HessianEvec}, tenGageCp1HessianEvec, 3, AIR_FALSE}, |
276 |
|
|
{tenGageCp1HessianEvec2, 3, 2, {tenGageCp1HessianEvec}, tenGageCp1HessianEvec, 6, AIR_FALSE}, |
277 |
|
|
{tenGageCa1Hessian, 9, 2, {tenGageTensorGradRotE, tenGageEvalHessian}, 0, 0, AIR_FALSE }, |
278 |
|
|
{tenGageCa1HessianEval, 3, 2, {tenGageCa1Hessian}, 0, 0, AIR_FALSE }, |
279 |
|
|
{tenGageCa1HessianEval0, 1, 2, {tenGageCa1HessianEval}, tenGageCa1HessianEval, 0, AIR_FALSE}, |
280 |
|
|
{tenGageCa1HessianEval1, 1, 2, {tenGageCa1HessianEval}, tenGageCa1HessianEval, 1, AIR_FALSE}, |
281 |
|
|
{tenGageCa1HessianEval2, 1, 2, {tenGageCa1HessianEval}, tenGageCa1HessianEval, 2, AIR_FALSE}, |
282 |
|
|
{tenGageCa1HessianEvec, 9, 2, {tenGageCa1Hessian}, 0, 0, AIR_FALSE }, |
283 |
|
|
{tenGageCa1HessianEvec0, 3, 2, {tenGageCa1HessianEvec}, tenGageCa1HessianEvec, 0, AIR_FALSE}, |
284 |
|
|
{tenGageCa1HessianEvec1, 3, 2, {tenGageCa1HessianEvec}, tenGageCa1HessianEvec, 3, AIR_FALSE}, |
285 |
|
|
{tenGageCa1HessianEvec2, 3, 2, {tenGageCa1HessianEvec}, tenGageCa1HessianEvec, 6, AIR_FALSE}, |
286 |
|
|
{tenGageFiberCurving, 1, 1, {tenGageRotTans, tenGageEvec}, 0, 0, AIR_FALSE }, |
287 |
|
|
{tenGageFiberDispersion, 1, 1, {tenGageRotTans, tenGageEvec}, 0, 0, AIR_FALSE }, |
288 |
|
|
{tenGageAniso, TEN_ANISO_MAX+1, 0, {tenGageEval0, tenGageEval1, tenGageEval2}, 0, 0, AIR_FALSE} |
289 |
|
|
}; |
290 |
|
|
|
291 |
|
|
void |
292 |
|
|
_tenGageIv3Print(FILE *file, gageContext *ctx, gagePerVolume *pvl) { |
293 |
|
|
double *iv3; |
294 |
|
|
int i, fd; |
295 |
|
|
|
296 |
|
|
fd = 2*ctx->radius; |
297 |
|
|
iv3 = pvl->iv3 + fd*fd*fd; |
298 |
|
|
fprintf(file, "iv3[]'s *Dxx* component:\n"); |
299 |
|
|
switch(fd) { |
300 |
|
|
case 2: |
301 |
|
|
fprintf(file, "% 10.4f % 10.4f\n", (float)iv3[6], (float)iv3[7]); |
302 |
|
|
fprintf(file, " % 10.4f % 10.4f\n\n", (float)iv3[4], (float)iv3[5]); |
303 |
|
|
fprintf(file, "% 10.4f % 10.4f\n", (float)iv3[2], (float)iv3[3]); |
304 |
|
|
fprintf(file, " % 10.4f % 10.4f\n", (float)iv3[0], (float)iv3[1]); |
305 |
|
|
break; |
306 |
|
|
case 4: |
307 |
|
|
for (i=3; i>=0; i--) { |
308 |
|
|
fprintf(file, "% 10.4f % 10.4f % 10.4f % 10.4f\n", |
309 |
|
|
(float)iv3[12+16*i], (float)iv3[13+16*i], |
310 |
|
|
(float)iv3[14+16*i], (float)iv3[15+16*i]); |
311 |
|
|
fprintf(file, " % 10.4f %c% 10.4f % 10.4f%c % 10.4f\n", |
312 |
|
|
(float)iv3[ 8+16*i], (i==1||i==2)?'\\':' ', |
313 |
|
|
(float)iv3[ 9+16*i], (float)iv3[10+16*i], (i==1||i==2)?'\\':' ', |
314 |
|
|
(float)iv3[11+16*i]); |
315 |
|
|
fprintf(file, " % 10.4f %c% 10.4f % 10.4f%c % 10.4f\n", |
316 |
|
|
(float)iv3[ 4+16*i], (i==1||i==2)?'\\':' ', |
317 |
|
|
(float)iv3[ 5+16*i], (float)iv3[ 6+16*i], (i==1||i==2)?'\\':' ', |
318 |
|
|
(float)iv3[ 7+16*i]); |
319 |
|
|
fprintf(file, " % 10.4f % 10.4f % 10.4f % 10.4f\n", |
320 |
|
|
(float)iv3[ 0+16*i], (float)iv3[ 1+16*i], |
321 |
|
|
(float)iv3[ 2+16*i], (float)iv3[ 3+16*i]); |
322 |
|
|
if (i) fprintf(file, "\n"); |
323 |
|
|
} |
324 |
|
|
break; |
325 |
|
|
default: |
326 |
|
|
for (i=0; i<fd*fd*fd; i++) { |
327 |
|
|
fprintf(file, " iv3[% 3d,% 3d,% 3d] = % 10.4f\n", |
328 |
|
|
i%fd, (i/fd)%fd, i/(fd*fd), (float)iv3[i]); |
329 |
|
|
} |
330 |
|
|
break; |
331 |
|
|
} |
332 |
|
|
return; |
333 |
|
|
} |
334 |
|
|
|
335 |
|
|
void |
336 |
|
|
_tenGageFilter(gageContext *ctx, gagePerVolume *pvl) { |
337 |
|
59400 |
char me[]="_tenGageFilter"; |
338 |
|
|
double *fw00, *fw11, *fw22, *ten, *tgrad, *thess; |
339 |
|
|
int fd; |
340 |
|
29700 |
gageScl3PFilter_t *filter[5] = {NULL, gageScl3PFilter2, gageScl3PFilter4, |
341 |
|
|
gageScl3PFilter6, gageScl3PFilter8}; |
342 |
|
|
unsigned int valIdx; |
343 |
|
|
|
344 |
|
29700 |
fd = 2*ctx->radius; |
345 |
|
29700 |
ten = pvl->directAnswer[tenGageTensor]; |
346 |
|
29700 |
tgrad = pvl->directAnswer[tenGageTensorGrad]; |
347 |
|
29700 |
thess = pvl->directAnswer[tenGageHessian]; |
348 |
✗✓ |
29700 |
if (!ctx->parm.k3pack) { |
349 |
|
|
fprintf(stderr, "!%s: sorry, 6pack filtering not implemented\n", me); |
350 |
|
|
return; |
351 |
|
|
} |
352 |
|
29700 |
fw00 = ctx->fw + fd*3*gageKernel00; |
353 |
|
29700 |
fw11 = ctx->fw + fd*3*gageKernel11; |
354 |
|
29700 |
fw22 = ctx->fw + fd*3*gageKernel22; |
355 |
|
|
/* perform the filtering */ |
356 |
|
|
/* HEY: we still want trilinear interpolation of confidence, no? */ |
357 |
✓✓ |
29700 |
if (fd <= 8) { |
358 |
✓✓ |
351000 |
for (valIdx=0; valIdx<7; valIdx++) { |
359 |
|
327600 |
filter[ctx->radius](ctx->shape, |
360 |
|
163800 |
pvl->iv3 + valIdx*fd*fd*fd, |
361 |
|
163800 |
pvl->iv2 + valIdx*fd*fd, |
362 |
|
163800 |
pvl->iv1 + valIdx*fd, |
363 |
|
|
fw00, fw11, fw22, |
364 |
|
163800 |
ten + valIdx, tgrad + valIdx*3, thess + valIdx*9, |
365 |
|
163800 |
pvl->needD); |
366 |
|
|
} |
367 |
|
|
} else { |
368 |
✓✓ |
94500 |
for (valIdx=0; valIdx<7; valIdx++) { |
369 |
|
88200 |
gageScl3PFilterN(ctx->shape, fd, |
370 |
|
44100 |
pvl->iv3 + valIdx*fd*fd*fd, |
371 |
|
44100 |
pvl->iv2 + valIdx*fd*fd, |
372 |
|
44100 |
pvl->iv1 + valIdx*fd, |
373 |
|
|
fw00, fw11, fw22, |
374 |
|
44100 |
ten + valIdx, tgrad + valIdx*3, thess + valIdx*9, |
375 |
|
44100 |
pvl->needD); |
376 |
|
|
} |
377 |
|
|
} |
378 |
|
|
|
379 |
|
29700 |
return; |
380 |
|
29700 |
} |
381 |
|
|
|
382 |
|
|
void |
383 |
|
|
_tenGageAnswer(gageContext *ctx, gagePerVolume *pvl) { |
384 |
|
59400 |
char me[]="_tenGageAnswer"; |
385 |
|
29700 |
double *tenAns, *evalAns, *evecAns, *vecTmp=NULL, *matTmp=NULL, |
386 |
|
|
*gradDtA=NULL, *gradDtB=NULL, *gradDtC=NULL, |
387 |
|
|
*gradDtD=NULL, *gradDtE=NULL, *gradDtF=NULL, |
388 |
|
|
*hessDtA=NULL, *hessDtB=NULL, *hessDtC=NULL, |
389 |
|
|
*hessDtD=NULL, *hessDtE=NULL, *hessDtF=NULL, |
390 |
|
|
*gradCbS=NULL, *gradCbB=NULL, *gradCbQ=NULL, *gradCbR=NULL, |
391 |
|
|
*hessCbS=NULL, *hessCbB=NULL, *hessCbQ=NULL, *hessCbR=NULL, |
392 |
|
29700 |
gradDdXYZ[21]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; |
393 |
|
|
double tmp0, tmp1, tmp3, magTmp=0, |
394 |
|
|
dtA=0, dtB=0, dtC=0, dtD=0, dtE=0, dtF=0, |
395 |
|
|
cbQ=0, cbR=0, cbA=0, cbB=0, cbC=0, cbS=0, |
396 |
|
|
gradCbA[3]={0,0,0}, gradCbC[3]={0,0,0}; |
397 |
|
|
double hessCbA[9]={0,0,0,0,0,0,0,0,0}, |
398 |
|
|
hessCbC[9]={0,0,0,0,0,0,0,0,0}; |
399 |
|
|
int ci; |
400 |
|
|
|
401 |
|
29700 |
tenAns = pvl->directAnswer[tenGageTensor]; |
402 |
|
29700 |
evalAns = pvl->directAnswer[tenGageEval]; |
403 |
|
29700 |
evecAns = pvl->directAnswer[tenGageEvec]; |
404 |
✓✗ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensor)) { |
405 |
|
|
/* done if doV */ |
406 |
|
|
/* HEY: this was prohibiting a Deft-related hack |
407 |
|
|
tenAns[0] = AIR_CLAMP(0, tenAns[0], 1); |
408 |
|
|
*/ |
409 |
|
|
/* HEY: and this was botching using 1-conf as potential energy for push |
410 |
|
|
tenAns[0] = AIR_MAX(0, tenAns[0]); |
411 |
|
|
*/ |
412 |
|
29700 |
dtA = tenAns[1]; |
413 |
|
29700 |
dtB = tenAns[2]; |
414 |
|
29700 |
dtC = tenAns[3]; |
415 |
|
29700 |
dtD = tenAns[4]; |
416 |
|
29700 |
dtE = tenAns[5]; |
417 |
|
29700 |
dtF = tenAns[6]; |
418 |
✗✓ |
29700 |
if (ctx->verbose) { |
419 |
|
|
fprintf(stderr, "!%s: tensor = (%g) %g %g %g %g %g %g\n", me, |
420 |
|
|
tenAns[0], dtA, dtB, dtC, dtD, dtE, dtF); |
421 |
|
|
} |
422 |
|
|
} |
423 |
|
|
/* done if doV |
424 |
|
|
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageConfidence)) { |
425 |
|
|
} |
426 |
|
|
*/ |
427 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTrace)) { |
428 |
|
|
cbA = -(pvl->directAnswer[tenGageTrace][0] = dtA + dtD + dtF); |
429 |
|
|
} |
430 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageNorm)) { |
431 |
|
|
pvl->directAnswer[tenGageNorm][0] = |
432 |
|
|
sqrt(dtA*dtA + dtD*dtD + dtF*dtF + 2*dtB*dtB + 2*dtC*dtC + 2*dtE*dtE); |
433 |
|
|
} |
434 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageB)) { |
435 |
|
|
cbB = pvl->directAnswer[tenGageB][0] = |
436 |
|
|
dtA*dtD + dtA*dtF + dtD*dtF - dtB*dtB - dtC*dtC - dtE*dtE; |
437 |
|
|
} |
438 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageDet)) { |
439 |
|
|
cbC = -(pvl->directAnswer[tenGageDet][0] = |
440 |
|
|
2*dtB*dtC*dtE + dtA*dtD*dtF |
441 |
|
|
- dtC*dtC*dtD - dtA*dtE*dtE - dtB*dtB*dtF); |
442 |
|
|
} |
443 |
✓✗ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageS)) { |
444 |
|
29700 |
cbS = (pvl->directAnswer[tenGageS][0] = |
445 |
|
29700 |
dtA*dtA + dtD*dtD + dtF*dtF |
446 |
|
29700 |
+ 2*dtB*dtB + 2*dtC*dtC + 2*dtE*dtE); |
447 |
|
29700 |
} |
448 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageQ)) { |
449 |
|
|
cbQ = pvl->directAnswer[tenGageQ][0] = (cbS - cbB)/9; |
450 |
|
|
} |
451 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFA)) { |
452 |
|
|
tmp0 = (cbS |
453 |
|
|
? cbQ/cbS |
454 |
|
|
: 0); |
455 |
|
|
tmp0 = AIR_MAX(0, tmp0); |
456 |
|
|
pvl->directAnswer[tenGageFA][0] = 3*sqrt(tmp0); |
457 |
|
|
/* |
458 |
|
|
if (!AIR_EXISTS(pvl->directAnswer[tenGageFA][0])) { |
459 |
|
|
fprintf(stderr, "!%s: cbS = %g, cbQ = %g, cbQ/(epsilon + cbS) = %g\n" |
460 |
|
|
"tmp0 = max(0, cbQ/(epsilon + cbS)) = %g\n" |
461 |
|
|
"sqrt(tmp0) = %g --> %g\n", me, |
462 |
|
|
cbS, cbQ, cbQ/(epsilon + cbS), |
463 |
|
|
tmp0, sqrt(tmp0), pvl->directAnswer[tenGageFA][0]); |
464 |
|
|
} |
465 |
|
|
*/ |
466 |
|
|
} |
467 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageR)) { |
468 |
|
|
cbR = pvl->directAnswer[tenGageR][0] = |
469 |
|
|
(5*cbA*cbB - 27*cbC - 2*cbA*cbS)/54; |
470 |
|
|
} |
471 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageMode)) { |
472 |
|
|
double cbQQQ; |
473 |
|
|
cbQQQ = 2*AIR_MAX(0, cbQ*cbQ*cbQ); |
474 |
|
|
tmp0 = 1.41421356237309504880*(cbQQQ ? cbR/(sqrt(cbQQQ)) : 0); |
475 |
|
|
pvl->directAnswer[tenGageMode][0] = AIR_CLAMP(-1, tmp0, 1); |
476 |
|
|
} |
477 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmega)) { |
478 |
|
|
pvl->directAnswer[tenGageOmega][0] = |
479 |
|
|
pvl->directAnswer[tenGageFA][0]*(1+pvl->directAnswer[tenGageMode][0])/2; |
480 |
|
|
} |
481 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTheta)) { |
482 |
|
|
pvl->directAnswer[tenGageTheta][0] = |
483 |
|
|
acos(-pvl->directAnswer[tenGageMode][0])/AIR_PI; |
484 |
|
|
} |
485 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageModeWarp)) { |
486 |
|
|
pvl->directAnswer[tenGageModeWarp][0] = |
487 |
|
|
cos((1-pvl->directAnswer[tenGageMode][0])*AIR_PI/2); |
488 |
|
|
} |
489 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageEvec)) { |
490 |
|
|
/* we do the longer process to get eigenvectors, and in the process |
491 |
|
|
we always find the eigenvalues, whether or not they were asked for */ |
492 |
|
|
tenEigensolve_d(evalAns, evecAns, tenAns); |
493 |
✗✓ |
29700 |
} else if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageEval)) { |
494 |
|
|
/* else eigenvectors are NOT needed, but eigenvalues ARE needed */ |
495 |
|
|
tenEigensolve_d(evalAns, NULL, tenAns); |
496 |
|
|
} |
497 |
|
|
|
498 |
✗✓ |
59400 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageDelNormK2) |
499 |
✓✗ |
59400 |
|| GAGE_QUERY_ITEM_TEST(pvl->query, tenGageDelNormK3)) { |
500 |
|
|
double tmp[7]; |
501 |
|
|
tenInvariantGradientsK_d(tmp, |
502 |
|
|
pvl->directAnswer[tenGageDelNormK2], |
503 |
|
|
pvl->directAnswer[tenGageDelNormK3], |
504 |
|
|
tenAns, 0.0000001); |
505 |
|
|
} |
506 |
✗✓ |
59400 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageDelNormR1) |
507 |
✓✗ |
59400 |
|| GAGE_QUERY_ITEM_TEST(pvl->query, tenGageDelNormR2)) { |
508 |
|
|
tenInvariantGradientsR_d(pvl->directAnswer[tenGageDelNormR1], |
509 |
|
|
pvl->directAnswer[tenGageDelNormR2], |
510 |
|
|
pvl->directAnswer[tenGageDelNormK3], |
511 |
|
|
tenAns, 0.0000001); |
512 |
|
|
} |
513 |
|
|
|
514 |
✗✓ |
59400 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageDelNormPhi1) |
515 |
✓✗ |
59400 |
|| GAGE_QUERY_ITEM_TEST(pvl->query, tenGageDelNormPhi2) |
516 |
✓✗ |
59400 |
|| GAGE_QUERY_ITEM_TEST(pvl->query, tenGageDelNormPhi3)) { |
517 |
|
|
tenRotationTangents_d(pvl->directAnswer[tenGageDelNormPhi1], |
518 |
|
|
pvl->directAnswer[tenGageDelNormPhi2], |
519 |
|
|
pvl->directAnswer[tenGageDelNormPhi3], |
520 |
|
|
evecAns); |
521 |
|
|
} |
522 |
|
|
|
523 |
✓✗ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensorGrad)) { |
524 |
|
|
/* done if doD1 */ |
525 |
|
|
/* still have to set up pointer variables that item answers |
526 |
|
|
below will rely on as short-cuts */ |
527 |
|
29700 |
vecTmp = pvl->directAnswer[tenGageTensorGrad]; |
528 |
|
29700 |
gradDtA = vecTmp + 1*3; |
529 |
|
29700 |
gradDtB = vecTmp + 2*3; |
530 |
|
29700 |
gradDtC = vecTmp + 3*3; |
531 |
|
29700 |
gradDtD = vecTmp + 4*3; |
532 |
|
29700 |
gradDtE = vecTmp + 5*3; |
533 |
|
29700 |
gradDtF = vecTmp + 6*3; |
534 |
|
29700 |
TEN_T_SET(gradDdXYZ + 0*7, tenAns[0], |
535 |
|
|
gradDtA[0], gradDtB[0], gradDtC[0], |
536 |
|
|
gradDtD[0], gradDtE[0], |
537 |
|
|
gradDtF[0]); |
538 |
|
29700 |
TEN_T_SET(gradDdXYZ + 1*7, tenAns[0], |
539 |
|
|
gradDtA[1], gradDtB[1], gradDtC[1], |
540 |
|
|
gradDtD[1], gradDtE[1], |
541 |
|
|
gradDtF[1]); |
542 |
|
29700 |
TEN_T_SET(gradDdXYZ + 2*7, tenAns[0], |
543 |
|
|
gradDtA[2], gradDtB[2], gradDtC[2], |
544 |
|
|
gradDtD[2], gradDtE[2], |
545 |
|
|
gradDtF[2]); |
546 |
|
29700 |
} |
547 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensorGradMag)) { |
548 |
|
|
vecTmp = pvl->directAnswer[tenGageTensorGradMag]; |
549 |
|
|
vecTmp[0] = sqrt(TEN_T_DOT(gradDdXYZ + 0*7, gradDdXYZ + 0*7)); |
550 |
|
|
vecTmp[1] = sqrt(TEN_T_DOT(gradDdXYZ + 1*7, gradDdXYZ + 1*7)); |
551 |
|
|
vecTmp[2] = sqrt(TEN_T_DOT(gradDdXYZ + 2*7, gradDdXYZ + 2*7)); |
552 |
|
|
} |
553 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensorGradMag)) { |
554 |
|
|
pvl->directAnswer[tenGageTensorGradMagMag][0] = ELL_3V_LEN(vecTmp); |
555 |
|
|
} |
556 |
|
|
|
557 |
|
|
/* --- Trace --- */ |
558 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTraceGradVec)) { |
559 |
|
|
vecTmp = pvl->directAnswer[tenGageTraceGradVec]; |
560 |
|
|
ELL_3V_ADD3(vecTmp, gradDtA, gradDtD, gradDtF); |
561 |
|
|
ELL_3V_SCALE(gradCbA, -1, vecTmp); |
562 |
|
|
} |
563 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTraceGradMag)) { |
564 |
|
|
magTmp = pvl->directAnswer[tenGageTraceGradMag][0] = ELL_3V_LEN(vecTmp); |
565 |
|
|
} |
566 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTraceNormal)) { |
567 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageTraceNormal], |
568 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
569 |
|
|
} |
570 |
|
|
|
571 |
|
|
/* ---- Norm stuff handled after S */ |
572 |
|
|
|
573 |
|
|
/* --- B --- */ |
574 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageBGradVec)) { |
575 |
|
|
gradCbB = vecTmp = pvl->directAnswer[tenGageBGradVec]; |
576 |
|
|
ELL_3V_SCALE_ADD6(vecTmp, |
577 |
|
|
dtD + dtF, gradDtA, |
578 |
|
|
dtA + dtF, gradDtD, |
579 |
|
|
dtA + dtD, gradDtF, |
580 |
|
|
-2*dtB, gradDtB, |
581 |
|
|
-2*dtC, gradDtC, |
582 |
|
|
-2*dtE, gradDtE); |
583 |
|
|
} |
584 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageBGradMag)) { |
585 |
|
|
magTmp = pvl->directAnswer[tenGageBGradMag][0] = ELL_3V_LEN(vecTmp); |
586 |
|
|
} |
587 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageBNormal)) { |
588 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageBNormal], |
589 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
590 |
|
|
} |
591 |
|
|
/* --- Det --- */ |
592 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageDetGradVec)) { |
593 |
|
|
vecTmp = pvl->directAnswer[tenGageDetGradVec]; |
594 |
|
|
ELL_3V_SCALE_ADD6(vecTmp, |
595 |
|
|
dtD*dtF - dtE*dtE, gradDtA, |
596 |
|
|
2*(dtC*dtE - dtB*dtF), gradDtB, |
597 |
|
|
2*(dtB*dtE - dtC*dtD), gradDtC, |
598 |
|
|
dtA*dtF - dtC*dtC, gradDtD, |
599 |
|
|
2*(dtB*dtC - dtA*dtE), gradDtE, |
600 |
|
|
dtA*dtD - dtB*dtB, gradDtF); |
601 |
|
|
ELL_3V_SCALE(gradCbC, -1, vecTmp); |
602 |
|
|
} |
603 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageDetGradMag)) { |
604 |
|
|
magTmp = pvl->directAnswer[tenGageDetGradMag][0] = |
605 |
|
|
AIR_CAST(float, ELL_3V_LEN(vecTmp)); |
606 |
|
|
} |
607 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageDetNormal)) { |
608 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageDetNormal], |
609 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
610 |
|
|
} |
611 |
|
|
/* --- S --- */ |
612 |
✓✗ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageSGradVec)) { |
613 |
|
29700 |
gradCbS = vecTmp = pvl->directAnswer[tenGageSGradVec]; |
614 |
|
29700 |
ELL_3V_SCALE_ADD6(vecTmp, |
615 |
|
|
2*dtA, gradDtA, |
616 |
|
|
2*dtD, gradDtD, |
617 |
|
|
2*dtF, gradDtF, |
618 |
|
|
4*dtB, gradDtB, |
619 |
|
|
4*dtC, gradDtC, |
620 |
|
|
4*dtE, gradDtE); |
621 |
|
29700 |
} |
622 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageSGradMag)) { |
623 |
|
|
magTmp = pvl->directAnswer[tenGageSGradMag][0] = ELL_3V_LEN(vecTmp); |
624 |
|
|
} |
625 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageSNormal)) { |
626 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageSNormal], |
627 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
628 |
|
|
} |
629 |
|
|
|
630 |
|
|
/* --- Norm --- */ |
631 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageNormGradVec)) { |
632 |
|
|
double nslc; |
633 |
|
|
nslc = pvl->directAnswer[tenGageNorm][0]; |
634 |
|
|
nslc = nslc ? 1/(2*nslc) : 0.0; |
635 |
|
|
vecTmp = pvl->directAnswer[tenGageNormGradVec]; |
636 |
|
|
ELL_3V_SCALE(vecTmp, nslc, pvl->directAnswer[tenGageSGradVec]); |
637 |
|
|
} |
638 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageNormGradMag)) { |
639 |
|
|
magTmp = pvl->directAnswer[tenGageNormGradMag][0] = ELL_3V_LEN(vecTmp); |
640 |
|
|
} |
641 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageNormNormal)) { |
642 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageNormNormal], |
643 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
644 |
|
|
} |
645 |
|
|
|
646 |
|
|
/* --- Q --- */ |
647 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageQGradVec)) { |
648 |
|
|
gradCbQ = vecTmp = pvl->directAnswer[tenGageQGradVec]; |
649 |
|
|
ELL_3V_SCALE_ADD2(vecTmp, |
650 |
|
|
1.0/9, gradCbS, |
651 |
|
|
-1.0/9, gradCbB); |
652 |
|
|
|
653 |
|
|
} |
654 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageQGradMag)) { |
655 |
|
|
magTmp = pvl->directAnswer[tenGageQGradMag][0] = ELL_3V_LEN(vecTmp); |
656 |
|
|
} |
657 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageQNormal)) { |
658 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageQNormal], |
659 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
660 |
|
|
} |
661 |
|
|
/* --- FA --- */ |
662 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAGradVec)) { |
663 |
|
|
vecTmp = pvl->directAnswer[tenGageFAGradVec]; |
664 |
|
|
tmp3 = AIR_MAX(0, pvl->directAnswer[tenGageFA][0]); |
665 |
|
|
tmp0 = cbQ ? tmp3/(2*cbQ) : 0; |
666 |
|
|
tmp1 = cbS ? -tmp3/(2*cbS) : 0; |
667 |
|
|
ELL_3V_SCALE_ADD2(vecTmp, |
668 |
|
|
tmp0, gradCbQ, |
669 |
|
|
tmp1, gradCbS); |
670 |
|
|
} |
671 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAGradMag)) { |
672 |
|
|
magTmp = pvl->directAnswer[tenGageFAGradMag][0] = ELL_3V_LEN(vecTmp); |
673 |
|
|
} |
674 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFANormal)) { |
675 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageFANormal], |
676 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
677 |
|
|
} |
678 |
|
|
/* --- R --- */ |
679 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageRGradVec)) { |
680 |
|
|
gradCbR = vecTmp = pvl->directAnswer[tenGageRGradVec]; |
681 |
|
|
ELL_3V_SCALE_ADD4(vecTmp, |
682 |
|
|
(5*cbB - 2*cbS)/54, gradCbA, |
683 |
|
|
5*cbA/54, gradCbB, |
684 |
|
|
-0.5, gradCbC, |
685 |
|
|
-cbA/27, gradCbS); |
686 |
|
|
} |
687 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageRGradMag)) { |
688 |
|
|
magTmp = pvl->directAnswer[tenGageRGradMag][0] = ELL_3V_LEN(vecTmp); |
689 |
|
|
} |
690 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageRNormal)) { |
691 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageRNormal], |
692 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
693 |
|
|
} |
694 |
|
|
/* --- Mode --- */ |
695 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageModeGradVec)) { |
696 |
|
|
vecTmp = pvl->directAnswer[tenGageModeGradVec]; |
697 |
|
|
tmp1 = AIR_MAX(0, cbQ*cbQ*cbQ); |
698 |
|
|
tmp1 = tmp1 ? sqrt(1/tmp1) : 0; |
699 |
|
|
tmp0 = cbQ ? -tmp1*3*cbR/(2*cbQ) : 0; |
700 |
|
|
ELL_3V_SCALE_ADD2(vecTmp, |
701 |
|
|
tmp0, gradCbQ, |
702 |
|
|
tmp1, gradCbR); |
703 |
|
|
} |
704 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageModeGradMag)) { |
705 |
|
|
magTmp = pvl->directAnswer[tenGageModeGradMag][0] = ELL_3V_LEN(vecTmp); |
706 |
|
|
} |
707 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageModeNormal)) { |
708 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageModeNormal], |
709 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
710 |
|
|
} |
711 |
|
|
/* --- Theta --- */ |
712 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageThetaGradVec)) { |
713 |
|
|
vecTmp = pvl->directAnswer[tenGageThetaGradVec]; |
714 |
|
|
tmp1 = AIR_MAX(0, cbQ*cbQ*cbQ); |
715 |
|
|
tmp0 = tmp1 ? cbR*cbR/tmp1 : 0; |
716 |
|
|
tmp1 = sqrt(tmp1)*sqrt(1.0 - tmp0); |
717 |
|
|
tmp1 = tmp1 ? 1/(AIR_PI*tmp1) : 0.0; |
718 |
|
|
tmp0 = cbQ ? -tmp1*3*cbR/(2*cbQ) : 0.0; |
719 |
|
|
ELL_3V_SCALE_ADD2(vecTmp, |
720 |
|
|
tmp0, gradCbQ, |
721 |
|
|
tmp1, gradCbR); |
722 |
|
|
} |
723 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageThetaGradMag)) { |
724 |
|
|
magTmp = pvl->directAnswer[tenGageThetaGradMag][0] = ELL_3V_LEN(vecTmp); |
725 |
|
|
} |
726 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageThetaNormal)) { |
727 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageThetaNormal], |
728 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
729 |
|
|
} |
730 |
|
|
/* --- Omega --- */ |
731 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaGradVec)) { |
732 |
|
|
double fa, mode, *faGrad, *modeGrad; |
733 |
|
|
vecTmp = pvl->directAnswer[tenGageOmegaGradVec]; |
734 |
|
|
fa = pvl->directAnswer[tenGageFA][0]; |
735 |
|
|
mode = pvl->directAnswer[tenGageMode][0]; |
736 |
|
|
faGrad = pvl->directAnswer[tenGageFAGradVec]; |
737 |
|
|
modeGrad = pvl->directAnswer[tenGageModeGradVec]; |
738 |
|
|
ELL_3V_SCALE_ADD2(vecTmp, |
739 |
|
|
(1+mode)/2, faGrad, |
740 |
|
|
fa/2, modeGrad); |
741 |
|
|
} |
742 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaGradMag)) { |
743 |
|
|
magTmp = pvl->directAnswer[tenGageOmegaGradMag][0] = ELL_3V_LEN(vecTmp); |
744 |
|
|
} |
745 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaNormal)) { |
746 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageOmegaNormal], |
747 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
748 |
|
|
} |
749 |
|
|
|
750 |
|
|
#define SQRT_1_OVER_3 0.57735026918962576450 |
751 |
|
|
|
752 |
|
|
/* --- Invariant gradients + rotation tangents --- */ |
753 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageInvarKGrads)) { |
754 |
|
|
double mu1Grad[7], *mu2Grad, *skwGrad; |
755 |
|
|
|
756 |
|
|
TEN_T_SET(mu1Grad, 1, |
757 |
|
|
SQRT_1_OVER_3, 0, 0, |
758 |
|
|
SQRT_1_OVER_3, 0, |
759 |
|
|
SQRT_1_OVER_3); |
760 |
|
|
mu2Grad = pvl->directAnswer[tenGageDelNormK2]; |
761 |
|
|
skwGrad = pvl->directAnswer[tenGageDelNormK3]; |
762 |
|
|
|
763 |
|
|
ELL_3V_SET(pvl->directAnswer[tenGageInvarKGrads] + 0*3, |
764 |
|
|
TEN_T_DOT(mu1Grad, gradDdXYZ + 0*7), |
765 |
|
|
TEN_T_DOT(mu1Grad, gradDdXYZ + 1*7), |
766 |
|
|
TEN_T_DOT(mu1Grad, gradDdXYZ + 2*7)); |
767 |
|
|
ELL_3V_SET(pvl->directAnswer[tenGageInvarKGrads] + 1*3, |
768 |
|
|
TEN_T_DOT(mu2Grad, gradDdXYZ + 0*7), |
769 |
|
|
TEN_T_DOT(mu2Grad, gradDdXYZ + 1*7), |
770 |
|
|
TEN_T_DOT(mu2Grad, gradDdXYZ + 2*7)); |
771 |
|
|
ELL_3V_SET(pvl->directAnswer[tenGageInvarKGrads] + 2*3, |
772 |
|
|
TEN_T_DOT(skwGrad, gradDdXYZ + 0*7), |
773 |
|
|
TEN_T_DOT(skwGrad, gradDdXYZ + 1*7), |
774 |
|
|
TEN_T_DOT(skwGrad, gradDdXYZ + 2*7)); |
775 |
|
|
} |
776 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageInvarKGradMags)) { |
777 |
|
|
ELL_3V_SET(pvl->directAnswer[tenGageInvarKGradMags], |
778 |
|
|
ELL_3V_LEN(pvl->directAnswer[tenGageInvarKGrads] + 0*3), |
779 |
|
|
ELL_3V_LEN(pvl->directAnswer[tenGageInvarKGrads] + 1*3), |
780 |
|
|
ELL_3V_LEN(pvl->directAnswer[tenGageInvarKGrads] + 2*3)); |
781 |
|
|
} |
782 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageInvarRGrads)) { |
783 |
|
|
double *R1Grad, *R2Grad, *R3Grad; |
784 |
|
|
|
785 |
|
|
R1Grad = pvl->directAnswer[tenGageDelNormR1]; |
786 |
|
|
R2Grad = pvl->directAnswer[tenGageDelNormR2]; |
787 |
|
|
R3Grad = pvl->directAnswer[tenGageDelNormK3]; |
788 |
|
|
|
789 |
|
|
ELL_3V_SET(pvl->directAnswer[tenGageInvarRGrads] + 0*3, |
790 |
|
|
TEN_T_DOT(R1Grad, gradDdXYZ + 0*7), |
791 |
|
|
TEN_T_DOT(R1Grad, gradDdXYZ + 1*7), |
792 |
|
|
TEN_T_DOT(R1Grad, gradDdXYZ + 2*7)); |
793 |
|
|
ELL_3V_SET(pvl->directAnswer[tenGageInvarRGrads] + 1*3, |
794 |
|
|
TEN_T_DOT(R2Grad, gradDdXYZ + 0*7), |
795 |
|
|
TEN_T_DOT(R2Grad, gradDdXYZ + 1*7), |
796 |
|
|
TEN_T_DOT(R2Grad, gradDdXYZ + 2*7)); |
797 |
|
|
ELL_3V_SET(pvl->directAnswer[tenGageInvarRGrads] + 2*3, |
798 |
|
|
TEN_T_DOT(R3Grad, gradDdXYZ + 0*7), |
799 |
|
|
TEN_T_DOT(R3Grad, gradDdXYZ + 1*7), |
800 |
|
|
TEN_T_DOT(R3Grad, gradDdXYZ + 2*7)); |
801 |
|
|
} |
802 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageInvarRGradMags)) { |
803 |
|
|
ELL_3V_SET(pvl->directAnswer[tenGageInvarRGradMags], |
804 |
|
|
ELL_3V_LEN(pvl->directAnswer[tenGageInvarRGrads] + 0*3), |
805 |
|
|
ELL_3V_LEN(pvl->directAnswer[tenGageInvarRGrads] + 1*3), |
806 |
|
|
ELL_3V_LEN(pvl->directAnswer[tenGageInvarRGrads] + 2*3)); |
807 |
|
|
} |
808 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageEvalGrads)) { |
809 |
|
|
double matOut[9], tenOut[9], tmpRes[9], |
810 |
|
|
rel1, rel2, w1, w2, eps; |
811 |
|
|
unsigned int evi; |
812 |
|
|
for (evi=0; evi<=2; evi++) { |
813 |
|
|
ELL_3MV_OUTER(matOut, evecAns + evi*3, evecAns + evi*3); |
814 |
|
|
TEN_M2T(tenOut, matOut); |
815 |
|
|
ELL_3V_SET(tmpRes + evi*3, |
816 |
|
|
TEN_T_DOT(tenOut, gradDdXYZ + 0*7), |
817 |
|
|
TEN_T_DOT(tenOut, gradDdXYZ + 1*7), |
818 |
|
|
TEN_T_DOT(tenOut, gradDdXYZ + 2*7)); |
819 |
|
|
} |
820 |
|
|
|
821 |
|
|
/* Added 2008-06-27: In case there are duplicate eigenvalues, |
822 |
|
|
* average their derivatives to avoid visible artifacs in edge |
823 |
|
|
* maps. Provide a smooth transition to the ill-defined case */ |
824 |
|
|
eps=0.05; /* threshold at which we start the transition */ |
825 |
|
|
|
826 |
|
|
/* interpolation weights from relative eigenvalue distance */ |
827 |
|
|
rel1=(evalAns[0]-evalAns[1])/(fabs(evalAns[0])+fabs(evalAns[1])); |
828 |
|
|
rel2=(evalAns[1]-evalAns[2])/(fabs(evalAns[1])+fabs(evalAns[2])); |
829 |
|
|
w1=rel1/eps-1; |
830 |
|
|
w1*=w1; |
831 |
|
|
w2=rel2/eps-1; |
832 |
|
|
w2*=w2; |
833 |
|
|
|
834 |
|
|
if (rel1>eps) { |
835 |
|
|
ELL_3V_COPY(pvl->directAnswer[tenGageEvalGrads], tmpRes); |
836 |
|
|
} else { |
837 |
|
|
if (rel2>eps) { |
838 |
|
|
ELL_3V_SCALE_ADD2(pvl->directAnswer[tenGageEvalGrads], |
839 |
|
|
1-0.5*w1, tmpRes, 0.5*w1, tmpRes+3); |
840 |
|
|
} else { |
841 |
|
|
ELL_3V_SCALE_ADD3(pvl->directAnswer[tenGageEvalGrads], |
842 |
|
|
1-0.5*w1-w1*w2/6.0, tmpRes, |
843 |
|
|
0.5*w1-w1*w2/6.0, tmpRes+3, |
844 |
|
|
w1*w2/3.0, tmpRes+6); |
845 |
|
|
} |
846 |
|
|
} |
847 |
|
|
|
848 |
|
|
if (rel2>eps) { |
849 |
|
|
ELL_3V_COPY(pvl->directAnswer[tenGageEvalGrads]+6, tmpRes+6); |
850 |
|
|
} else { |
851 |
|
|
if (rel1>eps) { |
852 |
|
|
ELL_3V_SCALE_ADD2(pvl->directAnswer[tenGageEvalGrads]+6, |
853 |
|
|
1-0.5*w2, tmpRes+6, 0.5*w2, tmpRes+3); |
854 |
|
|
} else { |
855 |
|
|
ELL_3V_SCALE_ADD3(pvl->directAnswer[tenGageEvalGrads]+6, |
856 |
|
|
1-0.5*w2-w1*w2/6.0, tmpRes+6, |
857 |
|
|
0.5*w2-w1*w2/6.0, tmpRes+3, |
858 |
|
|
w1*w2/3.0, tmpRes); |
859 |
|
|
} |
860 |
|
|
} |
861 |
|
|
|
862 |
|
|
ELL_3V_ADD3(pvl->directAnswer[tenGageEvalGrads]+3, |
863 |
|
|
tmpRes, tmpRes+3, tmpRes+6); |
864 |
|
|
ELL_3V_ADD2(tmpRes, pvl->directAnswer[tenGageEvalGrads], |
865 |
|
|
pvl->directAnswer[tenGageEvalGrads]+6); |
866 |
|
|
ELL_3V_SUB(pvl->directAnswer[tenGageEvalGrads]+3, |
867 |
|
|
pvl->directAnswer[tenGageEvalGrads]+3, |
868 |
|
|
tmpRes); /* l2 = trace - l1 - l3 */ |
869 |
|
|
} |
870 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageRotTans)) { |
871 |
|
|
double phi1[7], phi2[7], phi3[7]; |
872 |
|
|
|
873 |
|
|
tenRotationTangents_d(phi1, phi2, phi3, evecAns); |
874 |
|
|
ELL_3V_SET(pvl->directAnswer[tenGageRotTans] + 0*3, |
875 |
|
|
TEN_T_DOT(phi1, gradDdXYZ + 0*7), |
876 |
|
|
TEN_T_DOT(phi1, gradDdXYZ + 1*7), |
877 |
|
|
TEN_T_DOT(phi1, gradDdXYZ + 2*7)); |
878 |
|
|
ELL_3V_SET(pvl->directAnswer[tenGageRotTans] + 1*3, |
879 |
|
|
TEN_T_DOT(phi2, gradDdXYZ + 0*7), |
880 |
|
|
TEN_T_DOT(phi2, gradDdXYZ + 1*7), |
881 |
|
|
TEN_T_DOT(phi2, gradDdXYZ + 2*7)); |
882 |
|
|
ELL_3V_SET(pvl->directAnswer[tenGageRotTans] + 2*3, |
883 |
|
|
TEN_T_DOT(phi3, gradDdXYZ + 0*7), |
884 |
|
|
TEN_T_DOT(phi3, gradDdXYZ + 1*7), |
885 |
|
|
TEN_T_DOT(phi3, gradDdXYZ + 2*7)); |
886 |
|
|
} |
887 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageRotTanMags)) { |
888 |
|
|
ELL_3V_SET(pvl->directAnswer[tenGageRotTanMags], |
889 |
|
|
ELL_3V_LEN(pvl->directAnswer[tenGageRotTans] + 0*3), |
890 |
|
|
ELL_3V_LEN(pvl->directAnswer[tenGageRotTans] + 1*3), |
891 |
|
|
ELL_3V_LEN(pvl->directAnswer[tenGageRotTans] + 2*3)); |
892 |
|
|
} |
893 |
|
|
/* --- C{l,p,a,lpmin}1 --- */ |
894 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCl1)) { |
895 |
|
|
tmp0 = tenAnisoEval_d(evalAns, tenAniso_Cl1); |
896 |
|
|
pvl->directAnswer[tenGageCl1][0] = AIR_CLAMP(0, tmp0, 1); |
897 |
|
|
} |
898 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCp1)) { |
899 |
|
|
tmp0 = tenAnisoEval_d(evalAns, tenAniso_Cp1); |
900 |
|
|
pvl->directAnswer[tenGageCp1][0] = AIR_CLAMP(0, tmp0, 1); |
901 |
|
|
} |
902 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCa1)) { |
903 |
|
|
tmp0 = tenAnisoEval_d(evalAns, tenAniso_Ca1); |
904 |
|
|
pvl->directAnswer[tenGageCa1][0] = AIR_CLAMP(0, tmp0, 1); |
905 |
|
|
} |
906 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageClpmin1)) { |
907 |
|
|
tmp0 = tenAnisoEval_d(evalAns, tenAniso_Clpmin1); |
908 |
|
|
pvl->directAnswer[tenGageClpmin1][0] = AIR_CLAMP(0, tmp0, 1); |
909 |
|
|
} |
910 |
|
|
/* --- C{l,p,a,lpmin}2 --- */ |
911 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCl2)) { |
912 |
|
|
tmp0 = tenAnisoEval_d(evalAns, tenAniso_Cl2); |
913 |
|
|
pvl->directAnswer[tenGageCl2][0] = AIR_CLAMP(0, tmp0, 1); |
914 |
|
|
} |
915 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCp2)) { |
916 |
|
|
tmp0 = tenAnisoEval_d(evalAns, tenAniso_Cp2); |
917 |
|
|
pvl->directAnswer[tenGageCp2][0] = AIR_CLAMP(0, tmp0, 1); |
918 |
|
|
} |
919 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCa2)) { |
920 |
|
|
tmp0 = tenAnisoEval_d(evalAns, tenAniso_Ca2); |
921 |
|
|
pvl->directAnswer[tenGageCa2][0] = AIR_CLAMP(0, tmp0, 1); |
922 |
|
|
} |
923 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageClpmin2)) { |
924 |
|
|
tmp0 = tenAnisoEval_d(evalAns, tenAniso_Clpmin2); |
925 |
|
|
pvl->directAnswer[tenGageClpmin2][0] = AIR_CLAMP(0, tmp0, 1); |
926 |
|
|
} |
927 |
|
|
/* --- Hessian madness (the derivative, not the soldier) --- */ |
928 |
✓✗ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageHessian)) { |
929 |
|
|
/* done if doD2; still have to set up pointers */ |
930 |
|
29700 |
matTmp = pvl->directAnswer[tenGageHessian]; |
931 |
|
29700 |
hessDtA = matTmp + 1*9; |
932 |
|
29700 |
hessDtB = matTmp + 2*9; |
933 |
|
29700 |
hessDtC = matTmp + 3*9; |
934 |
|
29700 |
hessDtD = matTmp + 4*9; |
935 |
|
29700 |
hessDtE = matTmp + 5*9; |
936 |
|
29700 |
hessDtF = matTmp + 6*9; |
937 |
|
29700 |
} |
938 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTraceHessian)) { |
939 |
|
|
ELL_3M_SCALE_ADD3(pvl->directAnswer[tenGageTraceHessian], |
940 |
|
|
1.0, hessDtA, |
941 |
|
|
1.0, hessDtD, |
942 |
|
|
1.0, hessDtF); |
943 |
|
|
ELL_3M_SCALE(hessCbA, -1, pvl->directAnswer[tenGageTraceHessian]); |
944 |
|
|
} |
945 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTraceHessianEvec)) { |
946 |
|
|
ell_3m_eigensolve_d(pvl->directAnswer[tenGageTraceHessianEval], |
947 |
|
|
pvl->directAnswer[tenGageTraceHessianEvec], |
948 |
|
|
pvl->directAnswer[tenGageTraceHessian], AIR_TRUE); |
949 |
✗✓ |
29700 |
} else if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTraceHessianEval)) { |
950 |
|
|
ell_3m_eigenvalues_d(pvl->directAnswer[tenGageTraceHessianEval], |
951 |
|
|
pvl->directAnswer[tenGageTraceHessian], AIR_TRUE); |
952 |
|
|
} |
953 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTraceHessianFrob)) { |
954 |
|
|
pvl->directAnswer[tenGageTraceHessianFrob][0] |
955 |
|
|
= ELL_3M_FROB(pvl->directAnswer[tenGageTraceHessian]); |
956 |
|
|
} |
957 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageBHessian)) { |
958 |
|
|
hessCbB = matTmp = pvl->directAnswer[tenGageBHessian]; |
959 |
|
|
ELL_3M_ZERO_SET(matTmp); |
960 |
|
|
ELL_3M_SCALE_INCR(matTmp, dtB, hessDtB); |
961 |
|
|
ELL_3M_SCALE_INCR(matTmp, dtC, hessDtC); |
962 |
|
|
ELL_3M_SCALE_INCR(matTmp, dtE, hessDtE); |
963 |
|
|
ELL_3MV_OUTER_INCR(matTmp, gradDtB, gradDtB); |
964 |
|
|
ELL_3MV_OUTER_INCR(matTmp, gradDtC, gradDtC); |
965 |
|
|
ELL_3MV_OUTER_INCR(matTmp, gradDtE, gradDtE); |
966 |
|
|
ELL_3M_SCALE(matTmp, -2, matTmp); |
967 |
|
|
ELL_3MV_OUTER_INCR(matTmp, gradDtD, gradDtA); |
968 |
|
|
ELL_3MV_OUTER_INCR(matTmp, gradDtF, gradDtA); |
969 |
|
|
ELL_3MV_OUTER_INCR(matTmp, gradDtA, gradDtD); |
970 |
|
|
ELL_3MV_OUTER_INCR(matTmp, gradDtF, gradDtD); |
971 |
|
|
ELL_3MV_OUTER_INCR(matTmp, gradDtA, gradDtF); |
972 |
|
|
ELL_3MV_OUTER_INCR(matTmp, gradDtD, gradDtF); |
973 |
|
|
ELL_3M_SCALE_INCR(matTmp, dtD + dtF, hessDtA); |
974 |
|
|
ELL_3M_SCALE_INCR(matTmp, dtA + dtF, hessDtD); |
975 |
|
|
ELL_3M_SCALE_INCR(matTmp, dtA + dtD, hessDtF); |
976 |
|
|
} |
977 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageDetHessian)) { |
978 |
|
|
double tmp[3]; |
979 |
|
|
matTmp = pvl->directAnswer[tenGageDetHessian]; |
980 |
|
|
ELL_3M_ZERO_SET(matTmp); |
981 |
|
|
ELL_3V_SCALE_ADD3(tmp, dtD, gradDtF, |
982 |
|
|
dtF, gradDtD, |
983 |
|
|
-2*dtE, gradDtE); |
984 |
|
|
ELL_3MV_OUTER_INCR(matTmp, tmp, gradDtA); |
985 |
|
|
ELL_3M_SCALE_INCR(matTmp, dtD*dtF - dtE*dtE, hessDtA); |
986 |
|
|
ELL_3V_SCALE_ADD4(tmp, 2*dtC, gradDtE, |
987 |
|
|
2*dtE, gradDtC, |
988 |
|
|
-2*dtB, gradDtF, |
989 |
|
|
-2*dtF, gradDtB); |
990 |
|
|
ELL_3MV_OUTER_INCR(matTmp, tmp, gradDtB); |
991 |
|
|
ELL_3M_SCALE_INCR(matTmp, 2*(dtC*dtE - dtB*dtF), hessDtB); |
992 |
|
|
ELL_3V_SCALE_ADD4(tmp, 2*dtB, gradDtE, |
993 |
|
|
2*dtE, gradDtB, |
994 |
|
|
-2*dtC, gradDtD, |
995 |
|
|
-2*dtD, gradDtC); |
996 |
|
|
ELL_3MV_OUTER_INCR(matTmp, tmp, gradDtC); |
997 |
|
|
ELL_3M_SCALE_INCR(matTmp, 2*(dtB*dtE - dtC*dtD), hessDtC); |
998 |
|
|
ELL_3V_SCALE_ADD3(tmp, dtA, gradDtF, |
999 |
|
|
dtF, gradDtA, |
1000 |
|
|
-2*dtC, gradDtC); |
1001 |
|
|
ELL_3MV_OUTER_INCR(matTmp, tmp, gradDtD); |
1002 |
|
|
ELL_3M_SCALE_INCR(matTmp, dtA*dtF - dtC*dtC, hessDtD); |
1003 |
|
|
ELL_3V_SCALE_ADD4(tmp, 2*dtB, gradDtC, |
1004 |
|
|
2*dtC, gradDtB, |
1005 |
|
|
-2*dtA, gradDtE, |
1006 |
|
|
-2*dtE, gradDtA); |
1007 |
|
|
ELL_3MV_OUTER_INCR(matTmp, tmp, gradDtE); |
1008 |
|
|
ELL_3M_SCALE_INCR(matTmp, 2*(dtB*dtC - dtA*dtE), hessDtE); |
1009 |
|
|
ELL_3V_SCALE_ADD3(tmp, dtA, gradDtD, |
1010 |
|
|
dtD, gradDtA, |
1011 |
|
|
-2*dtB, gradDtB); |
1012 |
|
|
ELL_3MV_OUTER_INCR(matTmp, tmp, gradDtF); |
1013 |
|
|
ELL_3M_SCALE_INCR(matTmp, dtA*dtD - dtB*dtB, hessDtF); |
1014 |
|
|
ELL_3M_SCALE(hessCbC, -1, pvl->directAnswer[tenGageDetHessian]); |
1015 |
|
|
} |
1016 |
✓✗ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageSHessian)) { |
1017 |
|
29700 |
hessCbS = matTmp = pvl->directAnswer[tenGageSHessian]; |
1018 |
|
29700 |
ELL_3M_ZERO_SET(matTmp); |
1019 |
|
29700 |
ELL_3M_SCALE_INCR(matTmp, dtB, hessDtB); |
1020 |
|
29700 |
ELL_3MV_OUTER_INCR(matTmp, gradDtB, gradDtB); |
1021 |
|
29700 |
ELL_3M_SCALE_INCR(matTmp, dtC, hessDtC); |
1022 |
|
29700 |
ELL_3MV_OUTER_INCR(matTmp, gradDtC, gradDtC); |
1023 |
|
29700 |
ELL_3M_SCALE_INCR(matTmp, dtE, hessDtE); |
1024 |
|
29700 |
ELL_3MV_OUTER_INCR(matTmp, gradDtE, gradDtE); |
1025 |
|
29700 |
ELL_3M_SCALE(matTmp, 2, matTmp); |
1026 |
|
29700 |
ELL_3M_SCALE_INCR(matTmp, dtA, hessDtA); |
1027 |
|
29700 |
ELL_3MV_OUTER_INCR(matTmp, gradDtA, gradDtA); |
1028 |
|
29700 |
ELL_3M_SCALE_INCR(matTmp, dtD, hessDtD); |
1029 |
|
29700 |
ELL_3MV_OUTER_INCR(matTmp, gradDtD, gradDtD); |
1030 |
|
29700 |
ELL_3M_SCALE_INCR(matTmp, dtF, hessDtF); |
1031 |
|
29700 |
ELL_3MV_OUTER_INCR(matTmp, gradDtF, gradDtF); |
1032 |
|
29700 |
ELL_3M_SCALE(matTmp, 2, matTmp); |
1033 |
|
29700 |
} |
1034 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageQHessian)) { |
1035 |
|
|
hessCbQ = pvl->directAnswer[tenGageQHessian]; |
1036 |
|
|
ELL_3M_SCALE_ADD2(hessCbQ, |
1037 |
|
|
1.0/9, hessCbS, |
1038 |
|
|
-1.0/9, hessCbB); |
1039 |
|
|
} |
1040 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAHessian)) { |
1041 |
|
|
double tmpQ, rQ, orQ, oQ, tmpS, rS, orS, oS; |
1042 |
|
|
tmpQ = AIR_MAX(0, cbQ); |
1043 |
|
|
tmpS = AIR_MAX(0, cbS); |
1044 |
|
|
oQ = tmpQ ? 1/tmpQ : 0; |
1045 |
|
|
oS = tmpS ? 1/tmpS : 0; |
1046 |
|
|
rQ = sqrt(tmpQ); |
1047 |
|
|
rS = sqrt(tmpS); |
1048 |
|
|
orQ = rQ ? 1/rQ : 0; |
1049 |
|
|
orS = rS ? 1/rS : 0; |
1050 |
|
|
matTmp = pvl->directAnswer[tenGageFAHessian]; |
1051 |
|
|
ELL_3M_ZERO_SET(matTmp); |
1052 |
|
|
ELL_3M_SCALE_INCR(matTmp, orS*orQ, hessCbQ); |
1053 |
|
|
ELL_3M_SCALE_INCR(matTmp, -rQ*orS*oS, hessCbS); |
1054 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, -orS*orQ*oQ/2, gradCbQ, gradCbQ); |
1055 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, 3*rQ*orS*oS*oS/2, gradCbS, gradCbS); |
1056 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, -orS*oS*orQ/2, gradCbS, gradCbQ); |
1057 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, -orQ*orS*oS/2, gradCbQ, gradCbS); |
1058 |
|
|
ELL_3M_SCALE(matTmp, 3.0/2, matTmp); |
1059 |
|
|
} |
1060 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAHessianEvec)) { |
1061 |
|
|
ell_3m_eigensolve_d(pvl->directAnswer[tenGageFAHessianEval], |
1062 |
|
|
pvl->directAnswer[tenGageFAHessianEvec], |
1063 |
|
|
pvl->directAnswer[tenGageFAHessian], AIR_TRUE); |
1064 |
✗✓ |
29700 |
} else if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAHessianEval)) { |
1065 |
|
|
ell_3m_eigenvalues_d(pvl->directAnswer[tenGageFAHessianEval], |
1066 |
|
|
pvl->directAnswer[tenGageFAHessian], AIR_TRUE); |
1067 |
|
|
} |
1068 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAHessianFrob)) { |
1069 |
|
|
pvl->directAnswer[tenGageFAHessianFrob][0] |
1070 |
|
|
= ELL_3M_FROB(pvl->directAnswer[tenGageFAHessian]); |
1071 |
|
|
} |
1072 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFARidgeSurfaceStrength)) { |
1073 |
|
|
double ev; |
1074 |
|
|
ev = -pvl->directAnswer[tenGageFAHessianEval][2]; |
1075 |
|
|
ev = AIR_MAX(0, ev); |
1076 |
|
|
pvl->directAnswer[tenGageFARidgeSurfaceStrength][0] = tenAns[0]*ev; |
1077 |
|
|
} |
1078 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAValleySurfaceStrength)) { |
1079 |
|
|
double ev; |
1080 |
|
|
ev = pvl->directAnswer[tenGageFAHessianEval][0]; |
1081 |
|
|
ev = AIR_MAX(0, ev); |
1082 |
|
|
pvl->directAnswer[tenGageFAValleySurfaceStrength][0] = tenAns[0]*ev; |
1083 |
|
|
} |
1084 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFALaplacian)) { |
1085 |
|
|
double *hess; |
1086 |
|
|
hess = pvl->directAnswer[tenGageFAHessian]; |
1087 |
|
|
pvl->directAnswer[tenGageFALaplacian][0] = hess[0] + hess[4] + hess[8]; |
1088 |
|
|
} |
1089 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAHessianEvalMode)) { |
1090 |
|
|
double *heval; |
1091 |
|
|
heval = pvl->directAnswer[tenGageFAHessianEval]; |
1092 |
|
|
pvl->directAnswer[tenGageFAHessianEvalMode][0] = airMode3_d(heval); |
1093 |
|
|
} |
1094 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFARidgeLineAlignment)) { |
1095 |
|
|
double *hev0, *dev0, dot, mde; |
1096 |
|
|
hev0 = pvl->directAnswer[tenGageFAHessianEvec0]; |
1097 |
|
|
dev0 = pvl->directAnswer[tenGageEvec0]; |
1098 |
|
|
dot = ELL_3V_DOT(hev0, dev0); |
1099 |
|
|
mde = pvl->directAnswer[tenGageFAHessianEvalMode][0]; |
1100 |
|
|
mde = AIR_AFFINE(-1, mde, 1, 0, 1); |
1101 |
|
|
pvl->directAnswer[tenGageFARidgeLineAlignment][0] = mde*dot*dot; |
1102 |
|
|
} |
1103 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFARidgeSurfaceAlignment)) { |
1104 |
|
|
double *hev2, *dev0, dot, mde; |
1105 |
|
|
hev2 = pvl->directAnswer[tenGageFAHessianEvec2]; |
1106 |
|
|
dev0 = pvl->directAnswer[tenGageEvec0]; |
1107 |
|
|
dot = ELL_3V_DOT(hev2, dev0); |
1108 |
|
|
mde = pvl->directAnswer[tenGageFAHessianEvalMode][0]; |
1109 |
|
|
mde = AIR_AFFINE(-1, mde, 1, 1, 0); |
1110 |
|
|
pvl->directAnswer[tenGageFARidgeSurfaceAlignment][0]= mde*(1-dot*dot); |
1111 |
|
|
} |
1112 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFA2ndDD)) { |
1113 |
|
|
double *hess, *norm, tmpv[3]; |
1114 |
|
|
hess = pvl->directAnswer[tenGageFAHessian]; |
1115 |
|
|
norm = pvl->directAnswer[tenGageFANormal]; |
1116 |
|
|
ELL_3MV_MUL(tmpv, hess, norm); |
1117 |
|
|
pvl->directAnswer[tenGageFA2ndDD][0] = ELL_3V_DOT(norm, tmpv); |
1118 |
|
|
} |
1119 |
|
|
|
1120 |
|
|
/* HEY: lots of this is copy/paste from gage/sclanswer.c */ |
1121 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAGeomTens)) { |
1122 |
|
|
double denom, *fahess, *fagmag, tmpMat[9], *fnorm, nPerp[9], sHess[9]; |
1123 |
|
|
|
1124 |
|
|
fahess = pvl->directAnswer[tenGageFAHessian]; |
1125 |
|
|
fagmag = pvl->directAnswer[tenGageFAGradMag]; |
1126 |
|
|
fnorm = pvl->directAnswer[tenGageFANormal]; |
1127 |
|
|
denom = (*fagmag) ? 1.0/(*fagmag) : 0.0; |
1128 |
|
|
ELL_3M_SCALE(sHess, denom, fahess); |
1129 |
|
|
ELL_3M_IDENTITY_SET(nPerp); |
1130 |
|
|
ELL_3MV_SCALE_OUTER_INCR(nPerp, -1, fnorm, fnorm); |
1131 |
|
|
/* gten = nPerp * sHess * nPerp */ |
1132 |
|
|
ELL_3M_MUL(tmpMat, sHess, nPerp); |
1133 |
|
|
ELL_3M_MUL(pvl->directAnswer[tenGageFAGeomTens], nPerp, tmpMat); |
1134 |
|
|
} |
1135 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFATotalCurv)) { |
1136 |
|
|
pvl->directAnswer[tenGageFATotalCurv][0] = |
1137 |
|
|
ELL_3M_FROB(pvl->directAnswer[tenGageFAGeomTens]); |
1138 |
|
|
} |
1139 |
✓✗✗✓
|
59400 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAKappa1) || |
1140 |
|
29700 |
GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAKappa2)) { |
1141 |
|
|
double *k1, *k2, T, N, D; |
1142 |
|
|
k1 = pvl->directAnswer[tenGageFAKappa1]; |
1143 |
|
|
k2 = pvl->directAnswer[tenGageFAKappa2]; |
1144 |
|
|
T = ELL_3M_TRACE(pvl->directAnswer[tenGageFAGeomTens]); |
1145 |
|
|
N = pvl->directAnswer[tenGageFATotalCurv][0]; |
1146 |
|
|
D = 2*N*N - T*T; |
1147 |
|
|
D = AIR_MAX(D, 0); |
1148 |
|
|
D = sqrt(D); |
1149 |
|
|
k1[0] = 0.5*(T + D); |
1150 |
|
|
k2[0] = 0.5*(T - D); |
1151 |
|
|
} |
1152 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAMeanCurv)) { |
1153 |
|
|
double k1, k2; |
1154 |
|
|
k1 = pvl->directAnswer[tenGageFAKappa1][0]; |
1155 |
|
|
k2 = pvl->directAnswer[tenGageFAKappa2][0]; |
1156 |
|
|
pvl->directAnswer[tenGageFAMeanCurv][0] = (k1 + k2)/2; |
1157 |
|
|
} |
1158 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAShapeIndex)) { |
1159 |
|
|
double k1, k2; |
1160 |
|
|
k1 = pvl->directAnswer[tenGageFAKappa1][0]; |
1161 |
|
|
k2 = pvl->directAnswer[tenGageFAKappa2][0]; |
1162 |
|
|
pvl->directAnswer[tenGageFAShapeIndex][0] = |
1163 |
|
|
-(2/AIR_PI)*atan2(k1 + k2, k1 - k2); |
1164 |
|
|
} |
1165 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAGaussCurv)) { |
1166 |
|
|
double k1, k2; |
1167 |
|
|
k1 = pvl->directAnswer[tenGageFAKappa1][0]; |
1168 |
|
|
k2 = pvl->directAnswer[tenGageFAKappa2][0]; |
1169 |
|
|
pvl->directAnswer[tenGageFAGaussCurv][0] = k1*k2; |
1170 |
|
|
} |
1171 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFACurvDir1)) { |
1172 |
|
|
double kk, tmpMat[9], tmpVec[3]; |
1173 |
|
|
kk = pvl->directAnswer[tenGageFAKappa2][0]; |
1174 |
|
|
ELL_3M_COPY(tmpMat, pvl->directAnswer[tenGageFAGeomTens]); |
1175 |
|
|
tmpMat[0] -= kk; tmpMat[4] -= kk; tmpMat[8] -= kk; |
1176 |
|
|
ell_3m_1d_nullspace_d(tmpVec, tmpMat); |
1177 |
|
|
ELL_3V_COPY(pvl->directAnswer[tenGageFACurvDir1], tmpVec); |
1178 |
|
|
} |
1179 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFACurvDir2)) { |
1180 |
|
|
double kk, tmpMat[9], tmpVec[3]; |
1181 |
|
|
kk = pvl->directAnswer[tenGageFAKappa1][0]; |
1182 |
|
|
ELL_3M_COPY(tmpMat, pvl->directAnswer[tenGageFAGeomTens]); |
1183 |
|
|
tmpMat[0] -= kk; tmpMat[4] -= kk; tmpMat[8] -= kk; |
1184 |
|
|
ell_3m_1d_nullspace_d(tmpVec, tmpMat); |
1185 |
|
|
ELL_3V_COPY(pvl->directAnswer[tenGageFACurvDir1], tmpVec); |
1186 |
|
|
} |
1187 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAFlowlineCurv)) { |
1188 |
|
|
double *fahess, *fnorm, *fagmag, denom, nPerp[9], nProj[9], |
1189 |
|
|
ncTen[9], sHess[9], tmpMat[9]; |
1190 |
|
|
fnorm = pvl->directAnswer[tenGageFANormal]; |
1191 |
|
|
fahess = pvl->directAnswer[tenGageFAHessian]; |
1192 |
|
|
fagmag = pvl->directAnswer[tenGageFAGradMag]; |
1193 |
|
|
ELL_3MV_OUTER(nProj, fnorm, fnorm); |
1194 |
|
|
ELL_3M_IDENTITY_SET(nPerp); |
1195 |
|
|
ELL_3M_SCALE_INCR(nPerp, -1, nProj); |
1196 |
|
|
denom = (*fagmag) ? 1.0/(*fagmag) : 0.0; |
1197 |
|
|
ELL_3M_SCALE(sHess, denom, fahess); |
1198 |
|
|
/* ncTen = nPerp * sHess * nProj */ |
1199 |
|
|
ELL_3M_MUL(tmpMat, sHess, nProj); |
1200 |
|
|
ELL_3M_MUL(ncTen, nPerp, tmpMat); |
1201 |
|
|
pvl->directAnswer[gageSclFlowlineCurv][0] = ELL_3M_FROB(ncTen); |
1202 |
|
|
} |
1203 |
|
|
|
1204 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageRHessian)) { |
1205 |
|
|
hessCbR = matTmp = pvl->directAnswer[tenGageRHessian]; |
1206 |
|
|
ELL_3M_ZERO_SET(matTmp); |
1207 |
|
|
ELL_3M_SCALE_INCR(matTmp, 5*cbB - 2*cbS, hessCbA); |
1208 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, 5, gradCbB, gradCbA); |
1209 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, -2, gradCbS, gradCbA); |
1210 |
|
|
ELL_3M_SCALE_INCR(matTmp, 5*cbA, hessCbB); |
1211 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, 5, gradCbA, gradCbB); |
1212 |
|
|
ELL_3M_SCALE_INCR(matTmp, -27, hessCbC); |
1213 |
|
|
ELL_3M_SCALE_INCR(matTmp, -2*cbA, hessCbS); |
1214 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, -2, gradCbA, gradCbS); |
1215 |
|
|
ELL_3M_SCALE(matTmp, 1.0/54, matTmp); |
1216 |
|
|
} |
1217 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageModeHessian)) { |
1218 |
|
|
double tmpQ, oQ, rQ; |
1219 |
|
|
tmpQ = AIR_MAX(0, cbQ); |
1220 |
|
|
rQ = sqrt(tmpQ); |
1221 |
|
|
oQ = tmpQ ? 1/tmpQ : 0; |
1222 |
|
|
matTmp = pvl->directAnswer[tenGageModeHessian]; |
1223 |
|
|
ELL_3M_ZERO_SET(matTmp); |
1224 |
|
|
ELL_3M_SCALE_INCR(matTmp, -(3.0/2)*cbR, hessCbQ); |
1225 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, (15.0/4)*cbR*oQ, gradCbQ, gradCbQ); |
1226 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, -(3.0/2), gradCbR, gradCbQ); |
1227 |
|
|
ELL_3M_SCALE_INCR(matTmp, cbQ, hessCbR); |
1228 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, -(3.0/2), gradCbQ, gradCbR); |
1229 |
|
|
tmp0 = (tmpQ && rQ) ? 1/(tmpQ*tmpQ*rQ) : 0.0; |
1230 |
|
|
ELL_3M_SCALE(matTmp, tmp0, matTmp); |
1231 |
|
|
} |
1232 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageModeHessianEvec)) { |
1233 |
|
|
ell_3m_eigensolve_d(pvl->directAnswer[tenGageModeHessianEval], |
1234 |
|
|
pvl->directAnswer[tenGageModeHessianEvec], |
1235 |
|
|
pvl->directAnswer[tenGageModeHessian], AIR_TRUE); |
1236 |
✗✓ |
29700 |
} else if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageModeHessianEval)) { |
1237 |
|
|
ell_3m_eigenvalues_d(pvl->directAnswer[tenGageModeHessianEval], |
1238 |
|
|
pvl->directAnswer[tenGageModeHessian], AIR_TRUE); |
1239 |
|
|
} |
1240 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageModeHessianFrob)) { |
1241 |
|
|
pvl->directAnswer[tenGageModeHessianFrob][0] |
1242 |
|
|
= ELL_3M_FROB(pvl->directAnswer[tenGageModeHessian]); |
1243 |
|
|
} |
1244 |
|
|
|
1245 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaHessian)) { |
1246 |
|
|
double fa, mode, *modeGrad, *faGrad, *modeHess, *faHess; |
1247 |
|
|
fa = pvl->directAnswer[tenGageFA][0]; |
1248 |
|
|
mode = pvl->directAnswer[tenGageMode][0]; |
1249 |
|
|
faGrad = pvl->directAnswer[tenGageFAGradVec]; |
1250 |
|
|
modeGrad = pvl->directAnswer[tenGageModeGradVec]; |
1251 |
|
|
faHess = pvl->directAnswer[tenGageFAHessian]; |
1252 |
|
|
modeHess = pvl->directAnswer[tenGageModeHessian]; |
1253 |
|
|
matTmp = pvl->directAnswer[tenGageOmegaHessian]; |
1254 |
|
|
ELL_3M_ZERO_SET(matTmp); |
1255 |
|
|
ELL_3M_SCALE_INCR(matTmp, (1+mode)/2, faHess); |
1256 |
|
|
ELL_3M_SCALE_INCR(matTmp, fa/2, modeHess); |
1257 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, 0.5, modeGrad, faGrad); |
1258 |
|
|
ELL_3MV_SCALE_OUTER_INCR(matTmp, 0.5, faGrad, modeGrad); |
1259 |
|
|
} |
1260 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaHessianEvec)) { |
1261 |
|
|
ell_3m_eigensolve_d(pvl->directAnswer[tenGageOmegaHessianEval], |
1262 |
|
|
pvl->directAnswer[tenGageOmegaHessianEvec], |
1263 |
|
|
pvl->directAnswer[tenGageOmegaHessian], AIR_TRUE); |
1264 |
✗✓ |
29700 |
} else if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaHessianEval)) { |
1265 |
|
|
ell_3m_eigenvalues_d(pvl->directAnswer[tenGageOmegaHessianEval], |
1266 |
|
|
pvl->directAnswer[tenGageOmegaHessian], AIR_TRUE); |
1267 |
|
|
} |
1268 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaLaplacian)) { |
1269 |
|
|
double *hess; |
1270 |
|
|
hess = pvl->directAnswer[tenGageOmegaHessian]; |
1271 |
|
|
pvl->directAnswer[tenGageOmegaLaplacian][0] = hess[0] + hess[4] + hess[8]; |
1272 |
|
|
} |
1273 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmega2ndDD)) { |
1274 |
|
|
double *hess, *norm, tmpv[3]; |
1275 |
|
|
hess = pvl->directAnswer[tenGageOmegaHessian]; |
1276 |
|
|
norm = pvl->directAnswer[tenGageOmegaNormal]; |
1277 |
|
|
ELL_3MV_MUL(tmpv, hess, norm); |
1278 |
|
|
pvl->directAnswer[tenGageOmega2ndDD][0] = ELL_3V_DOT(norm, tmpv); |
1279 |
|
|
} |
1280 |
|
|
|
1281 |
|
|
/* the copy-and-paste nature of this is really getting out of control ... */ |
1282 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaHessianContrTenEvec0)) { |
1283 |
|
|
double *hess, *evec, tmpv[3]; |
1284 |
|
|
hess = pvl->directAnswer[tenGageOmegaHessian]; |
1285 |
|
|
evec = pvl->directAnswer[tenGageEvec0]; |
1286 |
|
|
ELL_3MV_MUL(tmpv, hess, evec); |
1287 |
|
|
pvl->directAnswer[tenGageOmegaHessianContrTenEvec0][0] = ELL_3V_DOT(evec, tmpv); |
1288 |
|
|
} |
1289 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaHessianContrTenEvec1)) { |
1290 |
|
|
double *hess, *evec, tmpv[3]; |
1291 |
|
|
hess = pvl->directAnswer[tenGageOmegaHessian]; |
1292 |
|
|
evec = pvl->directAnswer[tenGageEvec1]; |
1293 |
|
|
ELL_3MV_MUL(tmpv, hess, evec); |
1294 |
|
|
pvl->directAnswer[tenGageOmegaHessianContrTenEvec1][0] = ELL_3V_DOT(evec, tmpv); |
1295 |
|
|
} |
1296 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaHessianContrTenEvec2)) { |
1297 |
|
|
double *hess, *evec, tmpv[3]; |
1298 |
|
|
hess = pvl->directAnswer[tenGageOmegaHessian]; |
1299 |
|
|
evec = pvl->directAnswer[tenGageEvec2]; |
1300 |
|
|
ELL_3MV_MUL(tmpv, hess, evec); |
1301 |
|
|
pvl->directAnswer[tenGageOmegaHessianContrTenEvec2][0] = ELL_3V_DOT(evec, tmpv); |
1302 |
|
|
} |
1303 |
|
|
|
1304 |
|
|
/* --- evec0 dot products */ |
1305 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTraceGradVecDotEvec0)) { |
1306 |
|
|
tmp0 = ELL_3V_DOT(evecAns + 0*3, pvl->directAnswer[tenGageTraceGradVec]); |
1307 |
|
|
pvl->directAnswer[tenGageTraceGradVecDotEvec0][0] = AIR_ABS(tmp0); |
1308 |
|
|
} |
1309 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTraceDiffusionAlign)) { |
1310 |
|
|
tmp0 = ELL_3V_DOT(evecAns + 0*3, pvl->directAnswer[tenGageTraceNormal]); |
1311 |
|
|
tmp0 = 1 - 2*acos(AIR_ABS(tmp0))/AIR_PI; |
1312 |
|
|
pvl->directAnswer[tenGageTraceDiffusionAlign][0] = tmp0; |
1313 |
|
|
} |
1314 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTraceDiffusionFraction)) { |
1315 |
|
|
double tmpv[3]; |
1316 |
|
|
TEN_TV_MUL(tmpv, tenAns, pvl->directAnswer[tenGageTraceNormal]); |
1317 |
|
|
tmp0 = ELL_3V_DOT(tmpv, pvl->directAnswer[tenGageTraceNormal]); |
1318 |
|
|
tmp0 /= TEN_T_TRACE(tenAns) ? TEN_T_TRACE(tenAns) : 1; |
1319 |
|
|
pvl->directAnswer[tenGageTraceDiffusionFraction][0] = tmp0; |
1320 |
|
|
} |
1321 |
|
|
|
1322 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFAGradVecDotEvec0)) { |
1323 |
|
|
tmp0 = ELL_3V_DOT(evecAns + 0*3, pvl->directAnswer[tenGageFAGradVec]); |
1324 |
|
|
pvl->directAnswer[tenGageFAGradVecDotEvec0][0] = AIR_ABS(tmp0); |
1325 |
|
|
} |
1326 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFADiffusionAlign)) { |
1327 |
|
|
tmp0 = ELL_3V_DOT(evecAns + 0*3, pvl->directAnswer[tenGageFANormal]); |
1328 |
|
|
tmp0 = 1 - 2*acos(AIR_ABS(tmp0))/AIR_PI; |
1329 |
|
|
pvl->directAnswer[tenGageFADiffusionAlign][0] = tmp0; |
1330 |
|
|
} |
1331 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFADiffusionFraction)) { |
1332 |
|
|
double tmpv[3]; |
1333 |
|
|
TEN_TV_MUL(tmpv, tenAns, pvl->directAnswer[tenGageFANormal]); |
1334 |
|
|
tmp0 = ELL_3V_DOT(tmpv, pvl->directAnswer[tenGageFANormal]); |
1335 |
|
|
tmp0 /= TEN_T_TRACE(tenAns) ? TEN_T_TRACE(tenAns) : 1; |
1336 |
|
|
pvl->directAnswer[tenGageFADiffusionFraction][0] = tmp0; |
1337 |
|
|
} |
1338 |
|
|
|
1339 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaGradVecDotEvec0)) { |
1340 |
|
|
tmp0 = ELL_3V_DOT(evecAns + 0*3, pvl->directAnswer[tenGageOmegaGradVec]); |
1341 |
|
|
pvl->directAnswer[tenGageOmegaGradVecDotEvec0][0] = AIR_ABS(tmp0); |
1342 |
|
|
} |
1343 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaDiffusionAlign)) { |
1344 |
|
|
tmp0 = ELL_3V_DOT(evecAns + 0*3, pvl->directAnswer[tenGageOmegaNormal]); |
1345 |
|
|
tmp0 = 1 - 2*acos(AIR_ABS(tmp0))/AIR_PI; |
1346 |
|
|
pvl->directAnswer[tenGageOmegaDiffusionAlign][0] = tmp0; |
1347 |
|
|
} |
1348 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageOmegaDiffusionFraction)) { |
1349 |
|
|
double tmpv[3]; |
1350 |
|
|
TEN_TV_MUL(tmpv, tenAns, pvl->directAnswer[tenGageOmegaNormal]); |
1351 |
|
|
tmp0 = ELL_3V_DOT(tmpv, pvl->directAnswer[tenGageOmegaNormal]); |
1352 |
|
|
tmp0 /= TEN_T_TRACE(tenAns) ? TEN_T_TRACE(tenAns) : 1; |
1353 |
|
|
pvl->directAnswer[tenGageOmegaDiffusionFraction][0] = tmp0; |
1354 |
|
|
} |
1355 |
|
|
|
1356 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageConfGradVecDotEvec0)) { |
1357 |
|
|
double *confGrad; |
1358 |
|
|
confGrad = pvl->directAnswer[tenGageTensorGrad]; |
1359 |
|
|
tmp0 = ELL_3V_DOT(evecAns + 0*3, confGrad); |
1360 |
|
|
pvl->directAnswer[tenGageConfGradVecDotEvec0][0] = AIR_ABS(tmp0); |
1361 |
|
|
} |
1362 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageConfDiffusionAlign)) { |
1363 |
|
|
double *confGrad, confNorm[3], tmp; |
1364 |
|
|
confGrad = pvl->directAnswer[tenGageTensorGrad]; |
1365 |
|
|
ELL_3V_NORM(confNorm, confGrad, tmp); |
1366 |
|
|
tmp0 = ELL_3V_DOT(evecAns + 0*3, confNorm); |
1367 |
|
|
tmp0 = 1 - 2*acos(AIR_ABS(tmp0))/AIR_PI; |
1368 |
|
|
pvl->directAnswer[tenGageConfDiffusionAlign][0] = tmp0; |
1369 |
|
|
} |
1370 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageConfDiffusionFraction)) { |
1371 |
|
|
double *confGrad, confNorm[3], tmp, tmpv[3]; |
1372 |
|
|
confGrad = pvl->directAnswer[tenGageTensorGrad]; |
1373 |
|
|
ELL_3V_NORM(confNorm, confGrad, tmp); |
1374 |
|
|
TEN_TV_MUL(tmpv, tenAns, confNorm); |
1375 |
|
|
tmp0 = ELL_3V_DOT(tmpv, confNorm); |
1376 |
|
|
tmp0 /= TEN_T_TRACE(tenAns) ? TEN_T_TRACE(tenAns) : 1; |
1377 |
|
|
pvl->directAnswer[tenGageConfDiffusionFraction][0] = tmp0; |
1378 |
|
|
} |
1379 |
|
|
|
1380 |
|
|
|
1381 |
|
|
/* --- Covariance --- */ |
1382 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCovariance)) { |
1383 |
|
|
unsigned int cc, tt, taa, tbb, |
1384 |
|
|
vijk, vii, vjj, vkk, fd, fddd; |
1385 |
|
|
double *cov, ww, wxx, wyy, wzz, ten[7]; |
1386 |
|
|
|
1387 |
|
|
cov = pvl->directAnswer[tenGageCovariance]; |
1388 |
|
|
/* HEY: casting because radius signed (shouldn't be) */ |
1389 |
|
|
fd = AIR_CAST(unsigned int, 2*ctx->radius); |
1390 |
|
|
fddd = fd*fd*fd; |
1391 |
|
|
|
1392 |
|
|
/* reset answer */ |
1393 |
|
|
for (cc=0; cc<21; cc++) { |
1394 |
|
|
cov[cc] = 0; |
1395 |
|
|
} |
1396 |
|
|
|
1397 |
|
|
ten[0] = 1; /* never used anyway */ |
1398 |
|
|
for (vijk=0; vijk<fddd; vijk++) { |
1399 |
|
|
vii = vijk % fd; |
1400 |
|
|
vjj = (vijk/fd) % fd; |
1401 |
|
|
vkk = vijk/fd/fd; |
1402 |
|
|
wxx = ctx->fw[vii + fd*(0 + 3*gageKernel00)]; |
1403 |
|
|
wyy = ctx->fw[vjj + fd*(1 + 3*gageKernel00)]; |
1404 |
|
|
wzz = ctx->fw[vkk + fd*(2 + 3*gageKernel00)]; |
1405 |
|
|
ww = wxx*wyy*wzz; |
1406 |
|
|
for (tt=1; tt<7; tt++) { |
1407 |
|
|
ten[tt] = ww*(pvl->iv3[vijk + fddd*tt] - tenAns[tt]); |
1408 |
|
|
} |
1409 |
|
|
|
1410 |
|
|
cc = 0; |
1411 |
|
|
for (taa=0; taa<6; taa++) { |
1412 |
|
|
for (tbb=taa; tbb<6; tbb++) { |
1413 |
|
|
/* HEY: do I really mean to have this factor in here? */ |
1414 |
|
|
/* it probably meant that the units in the IGRT TMI paper were |
1415 |
|
|
wrong by this factor ... */ |
1416 |
|
|
cov[cc] += 100000*ten[taa+1]*ten[tbb+1]; |
1417 |
|
|
cc++; |
1418 |
|
|
} |
1419 |
|
|
} |
1420 |
|
|
} |
1421 |
|
|
} |
1422 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCovarianceRGRT)) { |
1423 |
|
|
double *cov, *covr, *igrt[6]; |
1424 |
|
|
unsigned int taa, tbb, cc; |
1425 |
|
|
|
1426 |
|
|
cov = pvl->directAnswer[tenGageCovariance]; |
1427 |
|
|
covr = pvl->directAnswer[tenGageCovarianceRGRT]; |
1428 |
|
|
igrt[0] = pvl->directAnswer[tenGageDelNormR1]; |
1429 |
|
|
igrt[1] = pvl->directAnswer[tenGageDelNormR2]; |
1430 |
|
|
igrt[2] = pvl->directAnswer[tenGageDelNormK3]; |
1431 |
|
|
igrt[3] = pvl->directAnswer[tenGageDelNormPhi1]; |
1432 |
|
|
igrt[4] = pvl->directAnswer[tenGageDelNormPhi2]; |
1433 |
|
|
igrt[5] = pvl->directAnswer[tenGageDelNormPhi3]; |
1434 |
|
|
|
1435 |
|
|
cc = 0; |
1436 |
|
|
for (taa=0; taa<6; taa++) { |
1437 |
|
|
for (tbb=taa; tbb<6; tbb++) { |
1438 |
|
|
covr[cc] = tenDoubleContract_d(igrt[tbb], cov, igrt[taa]); |
1439 |
|
|
cc++; |
1440 |
|
|
} |
1441 |
|
|
} |
1442 |
|
|
} |
1443 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCovarianceKGRT)) { |
1444 |
|
|
double *cov, *covk, *igrt[6], delnormk1[7]; |
1445 |
|
|
unsigned int taa, tbb, cc; |
1446 |
|
|
|
1447 |
|
|
cov = pvl->directAnswer[tenGageCovariance]; |
1448 |
|
|
covk = pvl->directAnswer[tenGageCovarianceKGRT]; |
1449 |
|
|
TEN_T_SET(delnormk1, 1, 0.57735026, 0, 0, 0.57735026, 0, 0.57735026); |
1450 |
|
|
igrt[0] = delnormk1; |
1451 |
|
|
igrt[1] = pvl->directAnswer[tenGageDelNormK2]; |
1452 |
|
|
igrt[2] = pvl->directAnswer[tenGageDelNormK3]; |
1453 |
|
|
igrt[3] = pvl->directAnswer[tenGageDelNormPhi1]; |
1454 |
|
|
igrt[4] = pvl->directAnswer[tenGageDelNormPhi2]; |
1455 |
|
|
igrt[5] = pvl->directAnswer[tenGageDelNormPhi3]; |
1456 |
|
|
|
1457 |
|
|
cc = 0; |
1458 |
|
|
for (taa=0; taa<6; taa++) { |
1459 |
|
|
for (tbb=taa; tbb<6; tbb++) { |
1460 |
|
|
covk[cc] = tenDoubleContract_d(igrt[tbb], cov, igrt[taa]); |
1461 |
|
|
cc++; |
1462 |
|
|
} |
1463 |
|
|
} |
1464 |
|
|
} |
1465 |
|
|
|
1466 |
|
|
/* these are items that somewhat bypass the convolution result |
1467 |
|
|
(in tenGageTensor) because it has to do something else fancy |
1468 |
|
|
with the constituent tensors. This is young and hacky code; |
1469 |
|
|
and it may be that facilitating this kind of processing should |
1470 |
|
|
be better supported by the gage API ... */ |
1471 |
✓✗✗✓
|
59400 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensorLogEuclidean) || |
1472 |
✓✗ |
29700 |
GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensorQuatGeoLoxK) || |
1473 |
✓✗ |
29700 |
GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensorQuatGeoLoxR) || |
1474 |
|
29700 |
GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensorRThetaPhiLinear)) { |
1475 |
|
|
unsigned int vijk, vii, vjj, vkk, fd, fddd; |
1476 |
|
|
_tenGagePvlData *pvlData; |
1477 |
|
|
double *ans; |
1478 |
|
|
int qret; |
1479 |
|
|
|
1480 |
|
|
pvlData = AIR_CAST(_tenGagePvlData *, pvl->data); |
1481 |
|
|
/* HEY: casting because radius is signed (shouldn't be) */ |
1482 |
|
|
fd = AIR_CAST(unsigned int, 2*ctx->radius); |
1483 |
|
|
fddd = fd*fd*fd; |
1484 |
|
|
for (vijk=0; vijk<fddd; vijk++) { |
1485 |
|
|
double wxx, wyy, wzz; |
1486 |
|
|
unsigned int tt; |
1487 |
|
|
vii = vijk % fd; |
1488 |
|
|
vjj = (vijk/fd) % fd; |
1489 |
|
|
vkk = vijk/fd/fd; |
1490 |
|
|
wxx = ctx->fw[vii + fd*(0 + 3*gageKernel00)]; |
1491 |
|
|
wyy = ctx->fw[vjj + fd*(1 + 3*gageKernel00)]; |
1492 |
|
|
wzz = ctx->fw[vkk + fd*(2 + 3*gageKernel00)]; |
1493 |
|
|
pvlData->buffWght[vijk] = wxx*wyy*wzz; |
1494 |
|
|
for (tt=0; tt<7; tt++) { |
1495 |
|
|
pvlData->buffTen[tt + 7*vijk] = pvl->iv3[vijk + fddd*tt]; |
1496 |
|
|
} |
1497 |
|
|
} |
1498 |
|
|
|
1499 |
|
|
qret = 0; |
1500 |
|
|
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensorLogEuclidean)) { |
1501 |
|
|
ans = pvl->directAnswer[tenGageTensorLogEuclidean]; |
1502 |
|
|
qret = tenInterpN_d(ans, pvlData->buffTen, pvlData->buffWght, fddd, |
1503 |
|
|
tenInterpTypeLogLinear, pvlData->tip); |
1504 |
|
|
} |
1505 |
|
|
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensorQuatGeoLoxK)) { |
1506 |
|
|
ans = pvl->directAnswer[tenGageTensorQuatGeoLoxK]; |
1507 |
|
|
qret = tenInterpN_d(ans, pvlData->buffTen, pvlData->buffWght, fddd, |
1508 |
|
|
tenInterpTypeQuatGeoLoxK, pvlData->tip); |
1509 |
|
|
} |
1510 |
|
|
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensorQuatGeoLoxR)) { |
1511 |
|
|
ans = pvl->directAnswer[tenGageTensorQuatGeoLoxR]; |
1512 |
|
|
qret= tenInterpN_d(ans, pvlData->buffTen, pvlData->buffWght, fddd, |
1513 |
|
|
tenInterpTypeQuatGeoLoxR, pvlData->tip); |
1514 |
|
|
} |
1515 |
|
|
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensorRThetaPhiLinear)) { |
1516 |
|
|
ans = pvl->directAnswer[tenGageTensorRThetaPhiLinear]; |
1517 |
|
|
qret= tenInterpN_d(ans, pvlData->buffTen, pvlData->buffWght, fddd, |
1518 |
|
|
tenInterpTypeRThetaPhiLinear, pvlData->tip); |
1519 |
|
|
} |
1520 |
|
|
if (qret) { |
1521 |
|
|
char *lerr; |
1522 |
|
|
fprintf(stderr, "!%s: problem!!!\n %s", me, lerr = biffGet(TEN)); |
1523 |
|
|
free(lerr); |
1524 |
|
|
} |
1525 |
|
|
} |
1526 |
|
|
|
1527 |
|
|
/* --- cl/cp/ca gradients --- */ |
1528 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCl1GradVec)) { |
1529 |
|
|
vecTmp = pvl->directAnswer[tenGageCl1GradVec]; |
1530 |
|
|
|
1531 |
|
|
ELL_3V_SET(vecTmp, |
1532 |
|
|
(evalAns[0]*(-2*pvl->directAnswer[tenGageEvalGrads][3]-pvl->directAnswer[tenGageEvalGrads][6]) |
1533 |
|
|
+evalAns[1]*(2*pvl->directAnswer[tenGageEvalGrads][0]+pvl->directAnswer[tenGageEvalGrads][6]) |
1534 |
|
|
+evalAns[2]*(pvl->directAnswer[tenGageEvalGrads][0]-pvl->directAnswer[tenGageEvalGrads][3])) |
1535 |
|
|
/(pvl->directAnswer[tenGageTrace][0]*pvl->directAnswer[tenGageTrace][0]), |
1536 |
|
|
(evalAns[0]*(-2*pvl->directAnswer[tenGageEvalGrads][4]-pvl->directAnswer[tenGageEvalGrads][7]) |
1537 |
|
|
+evalAns[1]*(2*pvl->directAnswer[tenGageEvalGrads][1]+pvl->directAnswer[tenGageEvalGrads][7]) |
1538 |
|
|
+evalAns[2]*(pvl->directAnswer[tenGageEvalGrads][1]-pvl->directAnswer[tenGageEvalGrads][4])) |
1539 |
|
|
/(pvl->directAnswer[tenGageTrace][0]*pvl->directAnswer[tenGageTrace][0]), |
1540 |
|
|
(evalAns[0]*(-2*pvl->directAnswer[tenGageEvalGrads][5]-pvl->directAnswer[tenGageEvalGrads][8]) |
1541 |
|
|
+evalAns[1]*(2*pvl->directAnswer[tenGageEvalGrads][2]+pvl->directAnswer[tenGageEvalGrads][8]) |
1542 |
|
|
+evalAns[2]*(pvl->directAnswer[tenGageEvalGrads][2]-pvl->directAnswer[tenGageEvalGrads][5])) |
1543 |
|
|
/(pvl->directAnswer[tenGageTrace][0]*pvl->directAnswer[tenGageTrace][0])); |
1544 |
|
|
} |
1545 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCl1GradMag)) { |
1546 |
|
|
magTmp = pvl->directAnswer[tenGageCl1GradMag][0] = ELL_3V_LEN(vecTmp); |
1547 |
|
|
} |
1548 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCl1Normal)) { |
1549 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageCl1Normal], |
1550 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
1551 |
|
|
} |
1552 |
|
|
|
1553 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCp1GradVec)) { |
1554 |
|
|
vecTmp = pvl->directAnswer[tenGageCp1GradVec]; |
1555 |
|
|
|
1556 |
|
|
ELL_3V_SET(vecTmp, |
1557 |
|
|
2*(evalAns[0]*(pvl->directAnswer[tenGageEvalGrads][3]-pvl->directAnswer[tenGageEvalGrads][6]) |
1558 |
|
|
+evalAns[1]*(-pvl->directAnswer[tenGageEvalGrads][0]-2*pvl->directAnswer[tenGageEvalGrads][6]) |
1559 |
|
|
+evalAns[2]*(pvl->directAnswer[tenGageEvalGrads][0]+2*pvl->directAnswer[tenGageEvalGrads][3])) |
1560 |
|
|
/(pvl->directAnswer[tenGageTrace][0]*pvl->directAnswer[tenGageTrace][0]), |
1561 |
|
|
2*(evalAns[0]*(pvl->directAnswer[tenGageEvalGrads][4]-pvl->directAnswer[tenGageEvalGrads][7]) |
1562 |
|
|
+evalAns[1]*(-pvl->directAnswer[tenGageEvalGrads][1]-2*pvl->directAnswer[tenGageEvalGrads][7]) |
1563 |
|
|
+evalAns[2]*(pvl->directAnswer[tenGageEvalGrads][1]+2*pvl->directAnswer[tenGageEvalGrads][4])) |
1564 |
|
|
/(pvl->directAnswer[tenGageTrace][0]*pvl->directAnswer[tenGageTrace][0]), |
1565 |
|
|
2*(evalAns[0]*(pvl->directAnswer[tenGageEvalGrads][5]-pvl->directAnswer[tenGageEvalGrads][8]) |
1566 |
|
|
+evalAns[1]*(-pvl->directAnswer[tenGageEvalGrads][2]-2*pvl->directAnswer[tenGageEvalGrads][8]) |
1567 |
|
|
+evalAns[2]*(pvl->directAnswer[tenGageEvalGrads][2]+2*pvl->directAnswer[tenGageEvalGrads][5])) |
1568 |
|
|
/(pvl->directAnswer[tenGageTrace][0]*pvl->directAnswer[tenGageTrace][0])); |
1569 |
|
|
} |
1570 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCp1GradMag)) { |
1571 |
|
|
magTmp = pvl->directAnswer[tenGageCp1GradMag][0] = ELL_3V_LEN(vecTmp); |
1572 |
|
|
} |
1573 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCp1Normal)) { |
1574 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageCp1Normal], |
1575 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
1576 |
|
|
} |
1577 |
|
|
|
1578 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCa1GradVec)) { |
1579 |
|
|
vecTmp = pvl->directAnswer[tenGageCa1GradVec]; |
1580 |
|
|
|
1581 |
|
|
ELL_3V_SET(vecTmp, |
1582 |
|
|
-3*((evalAns[0]+evalAns[1])*pvl->directAnswer[tenGageEvalGrads][6] |
1583 |
|
|
-evalAns[2]*(pvl->directAnswer[tenGageEvalGrads][0]+pvl->directAnswer[tenGageEvalGrads][3])) |
1584 |
|
|
/(pvl->directAnswer[tenGageTrace][0]*pvl->directAnswer[tenGageTrace][0]), |
1585 |
|
|
-3*((evalAns[0]+evalAns[1])*pvl->directAnswer[tenGageEvalGrads][7] |
1586 |
|
|
-evalAns[2]*(pvl->directAnswer[tenGageEvalGrads][1]+pvl->directAnswer[tenGageEvalGrads][4])) |
1587 |
|
|
/(pvl->directAnswer[tenGageTrace][0]*pvl->directAnswer[tenGageTrace][0]), |
1588 |
|
|
-3*((evalAns[0]+evalAns[1])*pvl->directAnswer[tenGageEvalGrads][8] |
1589 |
|
|
-evalAns[2]*(pvl->directAnswer[tenGageEvalGrads][2]+pvl->directAnswer[tenGageEvalGrads][5])) |
1590 |
|
|
/(pvl->directAnswer[tenGageTrace][0]*pvl->directAnswer[tenGageTrace][0])); |
1591 |
|
|
} |
1592 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCa1GradMag)) { |
1593 |
|
|
magTmp = pvl->directAnswer[tenGageCa1GradMag][0] = ELL_3V_LEN(vecTmp); |
1594 |
|
|
} |
1595 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCa1Normal)) { |
1596 |
|
|
ELL_3V_SCALE(pvl->directAnswer[tenGageCa1Normal], |
1597 |
|
|
magTmp ? 1/magTmp : 0, vecTmp); |
1598 |
|
|
} |
1599 |
|
|
|
1600 |
|
|
/* --- tensor gradient, rotated into eigenframe of the tensor itself --- */ |
1601 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageTensorGradRotE)) { |
1602 |
|
|
/* confidence not affected by rotation */ |
1603 |
|
|
double evecsT[9], evecs[9], tmp[9], tmp2[9]; |
1604 |
|
|
double diff0, diff1, diff2, diffthresh; |
1605 |
|
|
double rdiff0, rdiff1, rdiff2; |
1606 |
|
|
unsigned int evi, tci; |
1607 |
|
|
|
1608 |
|
|
ELL_3V_COPY(pvl->directAnswer[tenGageTensorGradRotE], |
1609 |
|
|
pvl->directAnswer[tenGageTensorGrad]); |
1610 |
|
|
|
1611 |
|
|
/* pre-compute relative eval diffs to detect ill-conditioned case */ |
1612 |
|
|
diff0=(evalAns[0]-evalAns[1])/(fabs(evalAns[0])+fabs(evalAns[1])); |
1613 |
|
|
diff1=(evalAns[1]-evalAns[2])/(fabs(evalAns[1])+fabs(evalAns[2])); |
1614 |
|
|
diff2=(evalAns[0]-evalAns[2])/(fabs(evalAns[0])+fabs(evalAns[2])); |
1615 |
|
|
diffthresh=0.05; |
1616 |
|
|
|
1617 |
|
|
if (diff2>diffthresh) rdiff2=1.0; else rdiff2=diff2/diffthresh; |
1618 |
|
|
if (diff1>diffthresh) rdiff1=1.0; else rdiff1=diff1/diffthresh; |
1619 |
|
|
if (diff0>diffthresh) rdiff0=1.0; else rdiff0=diff0/diffthresh; |
1620 |
|
|
|
1621 |
|
|
ELL_3M_COPY(evecs,evecAns); |
1622 |
|
|
ELL_3M_TRANSPOSE(evecsT,evecs); |
1623 |
|
|
for (evi=0; evi<3; evi++) { |
1624 |
|
|
char sign; |
1625 |
|
|
/* first, simply rotate derivatives into eigenframe of value */ |
1626 |
|
|
TEN_T2M(tmp,gradDdXYZ + evi*7); |
1627 |
|
|
ell_3m_mul_d(tmp2,tmp,evecsT); |
1628 |
|
|
ell_3m_mul_d(tmp,evecs,tmp2); |
1629 |
|
|
|
1630 |
|
|
/* If necessary, perform a number of additional rotations to |
1631 |
|
|
* distribute eigenvalue derivatives equally in ill-defined |
1632 |
|
|
* cases. Explanation in Schultz and Seidel, "Using Eigenvalue |
1633 |
|
|
* Derivatives for Edge Detection in DT-MRI Data", DAGM 2008*/ |
1634 |
|
|
if (rdiff0<1.0) { |
1635 |
|
|
/* the goal is to find the smallest angle phi such that |
1636 |
|
|
* rotation by phi around z will result in tmp[0]=tmp[4], i.e.: |
1637 |
|
|
* |
1638 |
|
|
* cos(phi)^2*tmp[0]-2*sin(phi)*cos(phi)*tmp[1]+sin(phi)^2*tmp[4] = |
1639 |
|
|
* sin(phi)^2*tmp[0]+2*sin(phi)*cos(phi)*tmp[1]+cos(phi)^2*tmp[4] |
1640 |
|
|
* => |
1641 |
|
|
* tan(2*phi)=(tmp[0]-tmp[4])/(2*tmp[1]) |
1642 |
|
|
* |
1643 |
|
|
* we use atan2 to avoid potential problems with tmp[1]==0, |
1644 |
|
|
* but manipulate the signs of the arguments s.t. the result |
1645 |
|
|
* is always in [-pi/2,pi/2] (i.e., the smallest solution of |
1646 |
|
|
* the above equality) |
1647 |
|
|
*/ |
1648 |
|
|
|
1649 |
|
|
/* rotate around z axis */ |
1650 |
|
|
double phi, R[9], RT[9]; |
1651 |
|
|
sign = (tmp[0]-tmp[4])*tmp[1]>0?1:-1; |
1652 |
|
|
phi=0.5*atan2(sign*fabs(tmp[0]-tmp[4]),fabs(2*tmp[1])); |
1653 |
|
|
ELL_3M_ROTATE_Z_SET(R, (1.0-rdiff0)*phi); |
1654 |
|
|
ELL_3M_TRANSPOSE(RT,R); |
1655 |
|
|
ell_3m_mul_d(tmp2,tmp,RT); |
1656 |
|
|
ell_3m_mul_d(tmp,R,tmp2); |
1657 |
|
|
} |
1658 |
|
|
if (rdiff1<1.0) { |
1659 |
|
|
/* rotate around x axis */ |
1660 |
|
|
double phi, R[9], RT[9]; |
1661 |
|
|
sign = (tmp[4]-tmp[8])*tmp[5]>0?1:-1; |
1662 |
|
|
phi=0.5*atan2(sign*fabs(tmp[4]-tmp[8]),fabs(2*tmp[5])); |
1663 |
|
|
ELL_3M_ROTATE_X_SET(R, (1.0-rdiff1)*phi); |
1664 |
|
|
ELL_3M_TRANSPOSE(RT,R); |
1665 |
|
|
ell_3m_mul_d(tmp2,tmp,RT); |
1666 |
|
|
ell_3m_mul_d(tmp,R,tmp2); |
1667 |
|
|
} |
1668 |
|
|
if (rdiff2<1.0) { |
1669 |
|
|
double mean, submatrix[3], isoPhi, _gamma, beta, A, C, R[9],RT[9]; |
1670 |
|
|
int axis, midaxis, smallest; |
1671 |
|
|
|
1672 |
|
|
mean=(tmp[0]+tmp[4]+tmp[8])/3.0; |
1673 |
|
|
/* what's the median? */ |
1674 |
|
|
midaxis=0; |
1675 |
|
|
if ((tmp[0]>tmp[4] && tmp[4]>tmp[8])|| |
1676 |
|
|
(tmp[0]<tmp[4] && tmp[4]<tmp[8])) |
1677 |
|
|
midaxis=1; |
1678 |
|
|
else if ((tmp[4]>tmp[8] && tmp[8]>tmp[0])|| |
1679 |
|
|
(tmp[4]<tmp[8] && tmp[8]<tmp[0])) |
1680 |
|
|
midaxis=2; |
1681 |
|
|
/* do we first rotate around smallest or largest? */ |
1682 |
|
|
smallest = 0; |
1683 |
|
|
if (mean>tmp[4*midaxis]) { |
1684 |
|
|
smallest=1; |
1685 |
|
|
} |
1686 |
|
|
|
1687 |
|
|
sign=1; |
1688 |
|
|
if ((smallest && (tmp[0]<tmp[4] && tmp[0]<tmp[8])) || |
1689 |
|
|
(!smallest && (tmp[0]>tmp[4] && tmp[0]>tmp[8]))) { |
1690 |
|
|
axis=0; |
1691 |
|
|
submatrix[0]=tmp[4]; |
1692 |
|
|
submatrix[1]=tmp[5]; |
1693 |
|
|
submatrix[2]=tmp[8]; |
1694 |
|
|
if (midaxis!=1) sign=-1; |
1695 |
|
|
} else if ((smallest && (tmp[4]<tmp[0] && tmp[4]<tmp[8])) || |
1696 |
|
|
(!smallest && (tmp[4]>tmp[0] && tmp[4]>tmp[8]))) { |
1697 |
|
|
axis=1; |
1698 |
|
|
submatrix[0]=tmp[8]; |
1699 |
|
|
submatrix[1]=tmp[2]; |
1700 |
|
|
submatrix[2]=tmp[0]; |
1701 |
|
|
if (midaxis!=2) sign=-1; |
1702 |
|
|
} else { |
1703 |
|
|
axis=2; |
1704 |
|
|
submatrix[0]=tmp[0]; |
1705 |
|
|
submatrix[1]=tmp[1]; |
1706 |
|
|
submatrix[2]=tmp[4]; |
1707 |
|
|
if (midaxis!=0) sign=-1; |
1708 |
|
|
} |
1709 |
|
|
isoPhi=0.0f; |
1710 |
|
|
_gamma=sign*(submatrix[0]-submatrix[2]); |
1711 |
|
|
beta=-sign*2*submatrix[1]; |
1712 |
|
|
A=sqrt(_gamma*_gamma+beta*beta); |
1713 |
|
|
C=atan2(_gamma,beta); |
1714 |
|
|
isoPhi=0.5*(asin(2.0/A*(mean-0.5*(submatrix[0]+submatrix[2])))-C); |
1715 |
|
|
/* make sure we use the minimal rotation */ |
1716 |
|
|
isoPhi=asin(2.0/A*(mean-0.5*(submatrix[0]+submatrix[2]))); |
1717 |
|
|
if (isoPhi>0) { |
1718 |
|
|
if (fabs(AIR_PI-isoPhi-C)<fabs(isoPhi-C)) |
1719 |
|
|
isoPhi=0.5*(AIR_PI-isoPhi-C); |
1720 |
|
|
else |
1721 |
|
|
isoPhi=0.5*(isoPhi-C); |
1722 |
|
|
} else if (isoPhi<0) { |
1723 |
|
|
if (fabs(-AIR_PI-isoPhi-C)<fabs(isoPhi-C)) |
1724 |
|
|
isoPhi=0.5*(-AIR_PI-isoPhi-C); |
1725 |
|
|
else |
1726 |
|
|
isoPhi=0.5*(isoPhi-C); |
1727 |
|
|
} |
1728 |
|
|
|
1729 |
|
|
/* perform the rotation */ |
1730 |
|
|
switch (axis) { |
1731 |
|
|
case 0: |
1732 |
|
|
ELL_3M_ROTATE_X_SET(R, (1.0-rdiff2)*isoPhi); |
1733 |
|
|
break; |
1734 |
|
|
case 1: |
1735 |
|
|
ELL_3M_ROTATE_Y_SET(R, (1.0-rdiff2)*isoPhi); |
1736 |
|
|
break; |
1737 |
|
|
default: |
1738 |
|
|
ELL_3M_ROTATE_Z_SET(R, (1.0-rdiff2)*isoPhi); |
1739 |
|
|
break; |
1740 |
|
|
} |
1741 |
|
|
ELL_3M_TRANSPOSE(RT,R); |
1742 |
|
|
ell_3m_mul_d(tmp2,tmp,RT); |
1743 |
|
|
ell_3m_mul_d(tmp,R,tmp2); |
1744 |
|
|
|
1745 |
|
|
/* rotate around the now corrected evec */ |
1746 |
|
|
axis=midaxis; |
1747 |
|
|
switch (midaxis) { |
1748 |
|
|
case 0: |
1749 |
|
|
sign = (tmp[0]-tmp[4])*tmp[1]>0?1:-1; |
1750 |
|
|
ELL_3M_ROTATE_X_SET(R, (1.0-rdiff2)*0.5*atan2(sign*fabs(tmp[0]-tmp[4]),fabs(2*tmp[1]))); |
1751 |
|
|
break; |
1752 |
|
|
case 1: |
1753 |
|
|
sign = (tmp[8]-tmp[0])*tmp[2]>0?1:-1; |
1754 |
|
|
ELL_3M_ROTATE_Y_SET(R, (1.0-rdiff2)*0.5*atan2(sign*fabs(tmp[8]-tmp[0]),fabs(2*tmp[2]))); |
1755 |
|
|
break; |
1756 |
|
|
case 2: |
1757 |
|
|
sign = (tmp[4]-tmp[8])*tmp[5]>0?1:-1; |
1758 |
|
|
ELL_3M_ROTATE_Z_SET(R, (1.0-rdiff2)*0.5*atan2(sign*fabs(tmp[4]-tmp[8]),fabs(2*tmp[5]))); |
1759 |
|
|
break; |
1760 |
|
|
} |
1761 |
|
|
ELL_3M_TRANSPOSE(RT,R); |
1762 |
|
|
ell_3m_mul_d(tmp2,tmp,RT); |
1763 |
|
|
ell_3m_mul_d(tmp,R,tmp2); |
1764 |
|
|
} |
1765 |
|
|
/* Now, we can set the answer */ |
1766 |
|
|
TEN_M2T(tmp2,tmp); |
1767 |
|
|
for (tci=1; tci<7; tci++) { |
1768 |
|
|
pvl->directAnswer[tenGageTensorGradRotE][3*tci+evi] = tmp2[tci]; |
1769 |
|
|
} |
1770 |
|
|
} |
1771 |
|
|
} |
1772 |
|
|
|
1773 |
|
|
/* --- Eigenvalue Hessians: rotate partial second derivatives into |
1774 |
|
|
* eigenframe and take into account correction factor based on |
1775 |
|
|
* eigenvector derivatives --- */ |
1776 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageEvalHessian)) { |
1777 |
|
|
/* some aliases for convenience: */ |
1778 |
|
|
double *thess = pvl->directAnswer[tenGageHessian]; |
1779 |
|
|
double *evhess = pvl->directAnswer[tenGageEvalHessian]; |
1780 |
|
|
double *tgradE = pvl->directAnswer[tenGageTensorGradRotE]; |
1781 |
|
|
int dira, dirb; /* directions of first/second partial derivative */ |
1782 |
|
|
int k; /* number of the eigenvalue */ |
1783 |
|
|
double evecsT[9], tmp[9]; |
1784 |
|
|
/* HEY: some copy-paste from tenGageEvalGrads */ |
1785 |
|
|
double eps=0.05; /* treshold for handling degeneracies */ |
1786 |
|
|
double rel1, rel2, rel3, w1, w2; |
1787 |
|
|
/* interpolation weights from relative eigenvalue distance */ |
1788 |
|
|
rel1=evalAns[0]*evalAns[0]+evalAns[1]*evalAns[1]; |
1789 |
|
|
if (rel1>1e-10) |
1790 |
|
|
rel1=(evalAns[0]-evalAns[1])*(evalAns[0]-evalAns[1])/rel1; |
1791 |
|
|
rel2=evalAns[0]*evalAns[0]+evalAns[2]*evalAns[2]; |
1792 |
|
|
if (rel2>1e-10) |
1793 |
|
|
rel2=(evalAns[0]-evalAns[2])*(evalAns[0]-evalAns[2])/rel2; |
1794 |
|
|
rel3=evalAns[1]*evalAns[1]+evalAns[2]*evalAns[2]; |
1795 |
|
|
if (rel3>1e-10) |
1796 |
|
|
rel3=(evalAns[1]-evalAns[2])*(evalAns[1]-evalAns[2])/rel3; |
1797 |
|
|
w1=rel1/eps-1; |
1798 |
|
|
w1*=w1; |
1799 |
|
|
w2=rel2/eps-1; |
1800 |
|
|
w2*=w2; |
1801 |
|
|
|
1802 |
|
|
ELL_3M_TRANSPOSE(evecsT,evecAns); |
1803 |
|
|
for (dira=0; dira<3; dira++) { |
1804 |
|
|
for (dirb=0; dirb<=dira; dirb++) { /* exploit symmetry of Hessian */ |
1805 |
|
|
double rdiff1,rdiff2; |
1806 |
|
|
double l1res, l2res, l3res; |
1807 |
|
|
/* collect second partial derivatives in dira,dirb */ |
1808 |
|
|
double H[9]; |
1809 |
|
|
ELL_3V_SET(H, thess[9+3*dirb+dira], thess[18+3*dirb+dira], |
1810 |
|
|
thess[27+3*dirb+dira]); |
1811 |
|
|
ELL_3V_SET(H+3, thess[18+3*dirb+dira], thess[36+3*dirb+dira], |
1812 |
|
|
thess[45+3*dirb+dira]); |
1813 |
|
|
ELL_3V_SET(H+6, thess[27+3*dirb+dira], thess[45+3*dirb+dira], |
1814 |
|
|
thess[54+3*dirb+dira]); |
1815 |
|
|
/* rotate into eigenframe of value */ |
1816 |
|
|
ell_3m_mul_d(tmp,H,evecsT); |
1817 |
|
|
ell_3m_mul_d(H,evecAns,tmp); |
1818 |
|
|
|
1819 |
|
|
/* we have to divide by rdiff=lambda_1-lambda_2; the following |
1820 |
|
|
* is a heuristic to avoid numerical problems in case rdiff is |
1821 |
|
|
* near zero */ |
1822 |
|
|
if (rel1>eps) rdiff1=1.0/(evalAns[0]-evalAns[1]); |
1823 |
|
|
else if (rel1>1e-10) rdiff1=(rel1/eps)/(evalAns[0]-evalAns[1]); |
1824 |
|
|
else rdiff1=0; |
1825 |
|
|
|
1826 |
|
|
if (rel2>eps) rdiff2=1.0/(evalAns[0]-evalAns[2]); |
1827 |
|
|
else if (rel2>1e-10) rdiff2=(rel2/eps)/(evalAns[0]-evalAns[2]); |
1828 |
|
|
else rdiff2=0; |
1829 |
|
|
|
1830 |
|
|
l1res=H[0]+2*(tgradE[6+dira]*tgradE[6+dirb]*rdiff1+ |
1831 |
|
|
tgradE[9+dira]*tgradE[9+dirb]*rdiff2); |
1832 |
|
|
|
1833 |
|
|
if (rel1>eps) rdiff1=1.0/(evalAns[1]-evalAns[0]); |
1834 |
|
|
else if (rel1>1e-10) rdiff1=(rel1/eps)/(evalAns[1]-evalAns[0]); |
1835 |
|
|
else rdiff1=0; |
1836 |
|
|
|
1837 |
|
|
if (rel3>eps) rdiff2=1.0/(evalAns[1]-evalAns[2]); |
1838 |
|
|
else if (rel3>1e-10) rdiff2=(rel3/eps)/(evalAns[1]-evalAns[2]); |
1839 |
|
|
else rdiff2=0; |
1840 |
|
|
l2res=H[4]+2*(tgradE[6+dira]*tgradE[6+dirb]*rdiff1+ |
1841 |
|
|
tgradE[15+dira]*tgradE[15+dirb]*rdiff2); |
1842 |
|
|
|
1843 |
|
|
if (rel2>eps) rdiff1=1.0/(evalAns[2]-evalAns[0]); |
1844 |
|
|
else if (rel2>1e-10) rdiff1=(rel2/eps)/(evalAns[2]-evalAns[0]); |
1845 |
|
|
else rdiff1=0; |
1846 |
|
|
|
1847 |
|
|
if (rel3>eps) rdiff2=1.0/(evalAns[2]-evalAns[1]); |
1848 |
|
|
else if (rel3>1e-10) rdiff2=(rel3/eps)/(evalAns[2]-evalAns[1]); |
1849 |
|
|
else rdiff2=0; |
1850 |
|
|
|
1851 |
|
|
l3res=H[8]+2*(tgradE[9+dira]*tgradE[9+dirb]*rdiff1+ |
1852 |
|
|
tgradE[15+dira]*tgradE[15+dirb]*rdiff2); |
1853 |
|
|
|
1854 |
|
|
if (rel1>eps) |
1855 |
|
|
evhess[3*dirb+dira]=l1res; |
1856 |
|
|
else { |
1857 |
|
|
if (rel2>eps) |
1858 |
|
|
evhess[3*dirb+dira]=(1-0.5*w1)*l1res+0.5*w1*l2res; |
1859 |
|
|
else |
1860 |
|
|
evhess[3*dirb+dira]=(1-0.5*w1-w1*w2/6.0)*l1res+ |
1861 |
|
|
(0.5*w1-w1*w2/6.0)*l2res+ |
1862 |
|
|
w1*w2/3.0*l3res; |
1863 |
|
|
} |
1864 |
|
|
|
1865 |
|
|
if (rel2>eps) |
1866 |
|
|
evhess[18+3*dirb+dira]=l3res; |
1867 |
|
|
else { |
1868 |
|
|
if (rel1>eps) |
1869 |
|
|
evhess[18+3*dirb+dira]=(1-0.5*w2)*l3res+0.5*w2*l2res; |
1870 |
|
|
else |
1871 |
|
|
evhess[18+3*dirb+dira]=(1-0.5*w2-w1*w2/6.0)*l3res+ |
1872 |
|
|
(0.5*w2-w1*w2/6.0)*l2res+ |
1873 |
|
|
w1*w2/3.0*l1res; |
1874 |
|
|
} |
1875 |
|
|
evhess[9+3*dirb+dira]=l1res+l2res+l3res - evhess[3*dirb+dira] - |
1876 |
|
|
evhess[18+3*dirb+dira]; |
1877 |
|
|
} |
1878 |
|
|
} |
1879 |
|
|
for (dira=0; dira<2; dira++) { |
1880 |
|
|
for (dirb=dira+1; dirb<3; dirb++) { /* copy over symmetric values */ |
1881 |
|
|
for (k=0; k<3; k++) { |
1882 |
|
|
evhess[9*k+3*dirb+dira]=evhess[9*k+3*dira+dirb]; |
1883 |
|
|
} |
1884 |
|
|
} |
1885 |
|
|
} |
1886 |
|
|
} |
1887 |
|
|
|
1888 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCl1Hessian)) { |
1889 |
|
|
int dira, dirb; |
1890 |
|
|
double *cl1hess = pvl->directAnswer[tenGageCl1Hessian]; |
1891 |
|
|
double *tgradE = pvl->directAnswer[tenGageTensorGradRotE]; |
1892 |
|
|
double *evhess = pvl->directAnswer[tenGageEvalHessian]; |
1893 |
|
|
/* A and B come out of the quotient rule; cf. appendix of Schultz |
1894 |
|
|
* et al., TVCG 2010 for more details */ |
1895 |
|
|
double B = (evalAns[0]+evalAns[1]+evalAns[2])* |
1896 |
|
|
(evalAns[0]+evalAns[1]+evalAns[2]); |
1897 |
|
|
for (dira=0; dira<3; dira++) { |
1898 |
|
|
for (dirb=0; dirb<=dira; dirb++) { /* again, exploit Hessian symmetry */ |
1899 |
|
|
double A = evalAns[0]*(-2*tgradE[12+dira]-tgradE[18+dira])+ |
1900 |
|
|
evalAns[1]*(2*tgradE[3+dira]+tgradE[18+dira])+ |
1901 |
|
|
evalAns[2]*(tgradE[3+dira]-tgradE[12+dira]); |
1902 |
|
|
double Ad = tgradE[3+dirb]*(-2*tgradE[12+dira]-tgradE[18+dira])+ |
1903 |
|
|
evalAns[0]*(-2*evhess[9+3*dirb+dira]-evhess[18+3*dirb+dira])+ |
1904 |
|
|
tgradE[12+dirb]*(2*tgradE[3+dira]+tgradE[18+dira])+ |
1905 |
|
|
evalAns[1]*(2*evhess[3*dirb+dira]+evhess[18+3*dirb+dira])+ |
1906 |
|
|
tgradE[18+dirb]*(tgradE[3+dira]-tgradE[12+dira])+ |
1907 |
|
|
evalAns[2]*(evhess[3*dirb+dira]-evhess[9+3*dirb+dira]); |
1908 |
|
|
double Bd = 2*(evalAns[0]+evalAns[1]+evalAns[2])* |
1909 |
|
|
(tgradE[3+dirb]+tgradE[12+dirb]+tgradE[18+dirb]); |
1910 |
|
|
cl1hess[3*dirb+dira]=Ad/B-A/B*Bd/B; |
1911 |
|
|
} |
1912 |
|
|
} |
1913 |
|
|
for (dira=0; dira<2; dira++) { |
1914 |
|
|
for (dirb=dira+1; dirb<3; dirb++) { /* copy over symmetric values */ |
1915 |
|
|
cl1hess[3*dirb+dira]=cl1hess[3*dira+dirb]; |
1916 |
|
|
} |
1917 |
|
|
} |
1918 |
|
|
} |
1919 |
|
|
|
1920 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCl1HessianEvec)) { |
1921 |
|
|
ell_3m_eigensolve_d(pvl->directAnswer[tenGageCl1HessianEval], |
1922 |
|
|
pvl->directAnswer[tenGageCl1HessianEvec], |
1923 |
|
|
pvl->directAnswer[tenGageCl1Hessian], AIR_TRUE); |
1924 |
✗✓ |
29700 |
} else if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCl1HessianEval)) { |
1925 |
|
|
ell_3m_eigenvalues_d(pvl->directAnswer[tenGageCl1HessianEval], |
1926 |
|
|
pvl->directAnswer[tenGageCl1Hessian], AIR_TRUE); |
1927 |
|
|
} |
1928 |
|
|
|
1929 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCp1Hessian)) { |
1930 |
|
|
int dira, dirb; |
1931 |
|
|
double *cp1hess = pvl->directAnswer[tenGageCp1Hessian]; |
1932 |
|
|
double *tgradE = pvl->directAnswer[tenGageTensorGradRotE]; |
1933 |
|
|
double *evhess = pvl->directAnswer[tenGageEvalHessian]; |
1934 |
|
|
double B = (evalAns[0]+evalAns[1]+evalAns[2])* |
1935 |
|
|
(evalAns[0]+evalAns[1]+evalAns[2]); |
1936 |
|
|
for (dira=0; dira<3; dira++) { |
1937 |
|
|
for (dirb=0; dirb<=dira; dirb++) { /* again, exploit Hessian symmetry */ |
1938 |
|
|
double A = 2*(evalAns[0]*(tgradE[12+dira]-tgradE[18+dira])+ |
1939 |
|
|
evalAns[1]*(-tgradE[3+dira]-2*tgradE[18+dira])+ |
1940 |
|
|
evalAns[2]*(tgradE[3+dira]+2*tgradE[12+dira])); |
1941 |
|
|
double Ad=2*(evalAns[0]*(evhess[9+3*dirb+dira]-evhess[18+3*dirb+dira])+ |
1942 |
|
|
evalAns[1]*(-evhess[3*dirb+dira]-2*evhess[18+3*dirb+dira])+ |
1943 |
|
|
evalAns[2]*(evhess[3*dirb+dira]+2*evhess[9+3*dirb+dira])+ |
1944 |
|
|
tgradE[3+dirb]*(tgradE[12+dira]-tgradE[18+dira])+ |
1945 |
|
|
tgradE[12+dirb]*(-tgradE[3+dira]-2*tgradE[18+dira])+ |
1946 |
|
|
tgradE[18+dirb]*(tgradE[3+dira]+2*tgradE[12+dira])); |
1947 |
|
|
double Bd = 2*(evalAns[0]+evalAns[1]+evalAns[2])* |
1948 |
|
|
(tgradE[3+dirb]+tgradE[12+dirb]+tgradE[18+dirb]); |
1949 |
|
|
cp1hess[3*dirb+dira]=Ad/B-A/B*Bd/B; |
1950 |
|
|
} |
1951 |
|
|
} |
1952 |
|
|
for (dira=0; dira<2; dira++) { |
1953 |
|
|
for (dirb=dira+1; dirb<3; dirb++) { /* copy over symmetric values */ |
1954 |
|
|
cp1hess[3*dirb+dira]=cp1hess[3*dira+dirb]; |
1955 |
|
|
} |
1956 |
|
|
} |
1957 |
|
|
} |
1958 |
|
|
|
1959 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCp1HessianEvec)) { |
1960 |
|
|
ell_3m_eigensolve_d(pvl->directAnswer[tenGageCp1HessianEval], |
1961 |
|
|
pvl->directAnswer[tenGageCp1HessianEvec], |
1962 |
|
|
pvl->directAnswer[tenGageCp1Hessian], AIR_TRUE); |
1963 |
✗✓ |
29700 |
} else if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCp1HessianEval)) { |
1964 |
|
|
ell_3m_eigenvalues_d(pvl->directAnswer[tenGageCp1HessianEval], |
1965 |
|
|
pvl->directAnswer[tenGageCp1Hessian], AIR_TRUE); |
1966 |
|
|
} |
1967 |
|
|
|
1968 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCa1Hessian)) { |
1969 |
|
|
int dira, dirb; |
1970 |
|
|
double *ca1hess = pvl->directAnswer[tenGageCa1Hessian]; |
1971 |
|
|
double *tgradE = pvl->directAnswer[tenGageTensorGradRotE]; |
1972 |
|
|
double *evhess = pvl->directAnswer[tenGageEvalHessian]; |
1973 |
|
|
double B = (evalAns[0]+evalAns[1]+evalAns[2])* |
1974 |
|
|
(evalAns[0]+evalAns[1]+evalAns[2]); |
1975 |
|
|
for (dira=0; dira<3; dira++) { |
1976 |
|
|
for (dirb=0; dirb<=dira; dirb++) { /* again, exploit Hessian symmetry */ |
1977 |
|
|
double A = 3*(evalAns[0]*tgradE[18+dira]+evalAns[1]*tgradE[18+dira]+ |
1978 |
|
|
evalAns[2]*(-tgradE[3+dira]-tgradE[12+dira])); |
1979 |
|
|
double Ad = 3*(tgradE[3+dirb]*tgradE[18+dira]+ |
1980 |
|
|
tgradE[12+dirb]*tgradE[18+dira]+ |
1981 |
|
|
tgradE[18+dirb]*(-tgradE[3+dira]-tgradE[12+dira])+ |
1982 |
|
|
evalAns[0]*evhess[18+3*dirb+dira]+ |
1983 |
|
|
evalAns[1]*evhess[18+3*dirb+dira]+ |
1984 |
|
|
evalAns[2]*(-evhess[3*dirb+dira]-evhess[9+3*dirb+dira])); |
1985 |
|
|
double Bd = 2*(evalAns[0]+evalAns[1]+evalAns[2])* |
1986 |
|
|
(tgradE[3+dirb]+tgradE[12+dirb]+tgradE[18+dirb]); |
1987 |
|
|
/* above formulas are true for cs, so flip sign here */ |
1988 |
|
|
ca1hess[3*dirb+dira]=-Ad/B+A/B*Bd/B; |
1989 |
|
|
} |
1990 |
|
|
} |
1991 |
|
|
for (dira=0; dira<2; dira++) { |
1992 |
|
|
for (dirb=dira+1; dirb<3; dirb++) { /* copy over symmetric values */ |
1993 |
|
|
ca1hess[3*dirb+dira]=ca1hess[3*dira+dirb]; |
1994 |
|
|
} |
1995 |
|
|
} |
1996 |
|
|
} |
1997 |
|
|
|
1998 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCa1HessianEvec)) { |
1999 |
|
|
ell_3m_eigensolve_d(pvl->directAnswer[tenGageCa1HessianEval], |
2000 |
|
|
pvl->directAnswer[tenGageCa1HessianEvec], |
2001 |
|
|
pvl->directAnswer[tenGageCa1Hessian], AIR_TRUE); |
2002 |
✗✓ |
29700 |
} else if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageCa1HessianEval)) { |
2003 |
|
|
ell_3m_eigenvalues_d(pvl->directAnswer[tenGageCa1HessianEval], |
2004 |
|
|
pvl->directAnswer[tenGageCa1Hessian], AIR_TRUE); |
2005 |
|
|
} |
2006 |
|
|
|
2007 |
✗✓ |
59400 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFiberCurving) |
2008 |
✓✗ |
59400 |
|| GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFiberDispersion)) { |
2009 |
|
|
double rtout[7], evout[7]; |
2010 |
|
|
TEN_T3V_OUTER(rtout, pvl->directAnswer[tenGageRotTans] + 1*3); |
2011 |
|
|
TEN_T3V_OUTER_INCR(rtout, pvl->directAnswer[tenGageRotTans] + 2*3); |
2012 |
|
|
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFiberCurving)) { |
2013 |
|
|
TEN_T3V_OUTER(evout, pvl->directAnswer[tenGageEvec0]); |
2014 |
|
|
pvl->directAnswer[tenGageFiberCurving][0] = TEN_T_DOT(rtout, evout); |
2015 |
|
|
/* pvl->directAnswer[tenGageFiberCurving][0] *= 100000; */ |
2016 |
|
|
} |
2017 |
|
|
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageFiberDispersion)) { |
2018 |
|
|
TEN_T3V_OUTER(evout, pvl->directAnswer[tenGageEvec1]); |
2019 |
|
|
TEN_T3V_OUTER_INCR(evout, pvl->directAnswer[tenGageEvec2]); |
2020 |
|
|
pvl->directAnswer[tenGageFiberDispersion][0] = TEN_T_DOT(rtout, evout); |
2021 |
|
|
/* pvl->directAnswer[tenGageFiberDispersion][0] *= 100000; */ |
2022 |
|
|
} |
2023 |
|
|
} |
2024 |
|
|
|
2025 |
|
|
/* --- Aniso --- */ |
2026 |
✗✓ |
29700 |
if (GAGE_QUERY_ITEM_TEST(pvl->query, tenGageAniso)) { |
2027 |
|
|
for (ci=tenAnisoUnknown+1; ci<=TEN_ANISO_MAX; ci++) { |
2028 |
|
|
pvl->directAnswer[tenGageAniso][ci] = tenAnisoEval_d(evalAns, ci); |
2029 |
|
|
} |
2030 |
|
|
} |
2031 |
|
|
|
2032 |
|
|
return; |
2033 |
|
29700 |
} |
2034 |
|
|
|
2035 |
|
|
|
2036 |
|
|
void * |
2037 |
|
|
_tenGagePvlDataNew(const struct gageKind_t *kind) { |
2038 |
|
|
_tenGagePvlData *pvlData; |
2039 |
|
|
|
2040 |
|
|
AIR_UNUSED(kind); |
2041 |
|
16 |
pvlData = AIR_CALLOC(1, _tenGagePvlData); |
2042 |
✓✗ |
8 |
if (pvlData) { |
2043 |
|
8 |
pvlData->buffTen = NULL; |
2044 |
|
8 |
pvlData->buffWght = NULL; |
2045 |
|
8 |
pvlData->tip = tenInterpParmNew(); |
2046 |
|
8 |
} |
2047 |
|
8 |
return pvlData; |
2048 |
|
|
} |
2049 |
|
|
|
2050 |
|
|
void * |
2051 |
|
|
_tenGagePvlDataCopy(const struct gageKind_t *kind, |
2052 |
|
|
const void *_pvlDataOld) { |
2053 |
|
|
_tenGagePvlData *pvlDataNew, *pvlDataOld; |
2054 |
|
|
unsigned int num; |
2055 |
|
|
|
2056 |
|
|
AIR_UNUSED(kind); |
2057 |
|
16 |
pvlDataOld = AIR_CAST(_tenGagePvlData *, _pvlDataOld); |
2058 |
|
8 |
num = pvlDataOld->tip->allocLen; |
2059 |
|
8 |
pvlDataNew = AIR_CALLOC(1, _tenGagePvlData); |
2060 |
✓✗ |
8 |
if (pvlDataNew) { |
2061 |
|
8 |
pvlDataNew->buffTen = AIR_CALLOC(7*num, double); |
2062 |
|
8 |
pvlDataNew->buffWght = AIR_CALLOC(num, double); |
2063 |
|
8 |
pvlDataNew->tip = tenInterpParmCopy(pvlDataOld->tip); |
2064 |
|
8 |
} |
2065 |
|
8 |
return pvlDataNew; |
2066 |
|
|
} |
2067 |
|
|
|
2068 |
|
|
void * |
2069 |
|
|
_tenGagePvlDataNix(const struct gageKind_t *kind, |
2070 |
|
|
void *_pvlData) { |
2071 |
|
|
_tenGagePvlData *pvlData; |
2072 |
|
|
|
2073 |
|
|
AIR_UNUSED(kind); |
2074 |
|
32 |
pvlData = AIR_CAST(_tenGagePvlData *, _pvlData); |
2075 |
|
16 |
airFree(pvlData->buffTen); |
2076 |
|
16 |
airFree(pvlData->buffWght); |
2077 |
|
16 |
tenInterpParmNix(pvlData->tip); |
2078 |
|
16 |
airFree(pvlData); |
2079 |
|
16 |
return NULL; |
2080 |
|
|
} |
2081 |
|
|
|
2082 |
|
|
int |
2083 |
|
|
_tenGagePvlDataUpdate(const struct gageKind_t *kind, |
2084 |
|
|
const gageContext *ctx, const gagePerVolume *pvl, |
2085 |
|
|
const void *_pvlData) { |
2086 |
|
|
_tenGagePvlData *pvlData; |
2087 |
|
|
unsigned int fd, num; |
2088 |
|
|
|
2089 |
|
|
AIR_UNUSED(kind); |
2090 |
|
|
AIR_UNUSED(pvl); |
2091 |
|
16 |
pvlData = AIR_CAST(_tenGagePvlData *, _pvlData); |
2092 |
|
8 |
fd = AIR_CAST(unsigned int, 2*ctx->radius); |
2093 |
|
8 |
num = fd*fd*fd; |
2094 |
✓✗ |
8 |
if (num != pvlData->tip->allocLen) { |
2095 |
|
|
/* HEY: no error checking */ |
2096 |
|
8 |
airFree(pvlData->buffTen); |
2097 |
|
8 |
pvlData->buffTen = NULL; |
2098 |
|
8 |
airFree(pvlData->buffWght); |
2099 |
|
8 |
pvlData->buffWght = NULL; |
2100 |
|
8 |
pvlData->buffTen = AIR_CALLOC(7*num, double); |
2101 |
|
8 |
pvlData->buffWght = AIR_CALLOC(num, double); |
2102 |
|
8 |
tenInterpParmBufferAlloc(pvlData->tip, num); |
2103 |
|
8 |
} |
2104 |
|
8 |
return 0; |
2105 |
|
|
} |
2106 |
|
|
|
2107 |
|
|
|
2108 |
|
|
gageKind |
2109 |
|
|
_tenGageKind = { |
2110 |
|
|
AIR_FALSE, /* statically allocated */ |
2111 |
|
|
"tensor", |
2112 |
|
|
&_tenGage, |
2113 |
|
|
1, |
2114 |
|
|
7, |
2115 |
|
|
TEN_GAGE_ITEM_MAX, |
2116 |
|
|
_tenGageTable, |
2117 |
|
|
_tenGageIv3Print, |
2118 |
|
|
_tenGageFilter, |
2119 |
|
|
_tenGageAnswer, |
2120 |
|
|
_tenGagePvlDataNew, |
2121 |
|
|
_tenGagePvlDataCopy, |
2122 |
|
|
_tenGagePvlDataNix, |
2123 |
|
|
_tenGagePvlDataUpdate, |
2124 |
|
|
NULL |
2125 |
|
|
}; |
2126 |
|
|
gageKind * |
2127 |
|
|
tenGageKind = &_tenGageKind; |