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 "teem/gage.h" |
25 |
|
|
|
26 |
|
|
#define PROBE "probePolynomial" |
27 |
|
|
|
28 |
|
|
/* |
29 |
|
|
** Tests: |
30 |
|
|
** |
31 |
|
|
*/ |
32 |
|
|
|
33 |
|
|
/* the hope is that without too much work it would be possible to |
34 |
|
|
increase this to 4th-order polynomials and 4th-order derivatives */ |
35 |
|
|
#define POWER_MAX 3 |
36 |
|
|
|
37 |
|
|
/* |
38 |
|
|
** polyeval takes a polynomial defined by coef -- |
39 |
|
|
** coef[i][j][k] is the coefficient for (x^i)*(y^j)*(z^k) -- |
40 |
|
|
** and evaluates at position pos[] the dx-th derivative along X, |
41 |
|
|
** the dy-th derivative along Y, and the dz-th derivative along Z, |
42 |
|
|
** dx = dy = dz = 0 to evaluate the polynomial itself |
43 |
|
|
*/ |
44 |
|
|
static double |
45 |
|
|
polyeval(double coef[POWER_MAX+1][POWER_MAX+1][POWER_MAX+1], |
46 |
|
|
unsigned int dx, unsigned int dy, unsigned int dz, |
47 |
|
|
const double pos[3]) { |
48 |
|
|
unsigned int xi, yi, zi; |
49 |
|
28597632 |
double tmp, ret, |
50 |
|
|
xp[POWER_MAX+1], yp[POWER_MAX+1], zp[POWER_MAX+1], |
51 |
|
14298816 |
dc[POWER_MAX+1][POWER_MAX+1] = { |
52 |
|
|
{1.0, 1.0, 1.0, 1.0}, /* how coeffs are scaled for 0th derivative */ |
53 |
|
|
{0.0, 1.0, 2.0, 3.0}, /* how coeffs are scaled for 1st derivative */ |
54 |
|
|
{0.0, 0.0, 2.0, 6.0}, /* how coeffs are scaled for 2nd derivative */ |
55 |
|
|
{0.0, 0.0, 0.0, 6.0}}; /* how coeffs are scaled for 3rd derivative */ |
56 |
|
|
|
57 |
|
|
tmp = 1.0; |
58 |
✓✓ |
142988160 |
for (xi=0; xi<=POWER_MAX; xi++) { |
59 |
|
57195264 |
xp[xi] = tmp; |
60 |
|
57195264 |
tmp *= pos[0]; |
61 |
|
|
} |
62 |
|
|
tmp = 1.0; |
63 |
✓✓ |
142988160 |
for (yi=0; yi<=POWER_MAX; yi++) { |
64 |
|
57195264 |
yp[yi] = tmp; |
65 |
|
57195264 |
tmp *= pos[1]; |
66 |
|
|
} |
67 |
|
|
tmp = 1.0; |
68 |
✓✓ |
142988160 |
for (zi=0; zi<=POWER_MAX; zi++) { |
69 |
|
57195264 |
zp[zi] = tmp; |
70 |
|
57195264 |
tmp *= pos[2]; |
71 |
|
|
} |
72 |
|
|
|
73 |
|
|
ret = 0.0; |
74 |
✓✓ |
142976810 |
for (xi=dx; xi<=POWER_MAX; xi++) { |
75 |
✓✓ |
571852760 |
for (yi=dy; yi<=POWER_MAX; yi++) { |
76 |
✓✓ |
2287204470 |
for (zi=dz; zi<=POWER_MAX; zi++) { |
77 |
|
1829730888 |
ret += (coef[xi][yi][zi]*xp[xi-dx]*yp[yi-dy]*zp[zi-dz] |
78 |
|
914865444 |
*dc[dx][xi]*dc[dy][yi]*dc[dz][zi]); |
79 |
|
|
} |
80 |
|
|
} |
81 |
|
|
} |
82 |
|
|
|
83 |
|
14298816 |
return ret; |
84 |
|
14298816 |
} |
85 |
|
|
|
86 |
|
|
static void |
87 |
|
|
randvec(double vec[NRRD_SPACE_DIM_MAX], double lenexp, double lensig, |
88 |
|
|
airRandMTState *rng) { |
89 |
|
224 |
double tmp, clen[2], len; |
90 |
|
|
|
91 |
|
112 |
airNormalRand_r(vec + 0, vec + 1, rng); |
92 |
|
112 |
airNormalRand_r(vec + 2, NULL, rng); |
93 |
|
112 |
ELL_3V_NORM(vec, vec, tmp); |
94 |
|
112 |
airNormalRand_r(clen + 0, clen + 1, rng); |
95 |
|
|
/* rician noise, actually */ |
96 |
|
112 |
clen[0] = lenexp + lensig*clen[0]; |
97 |
|
112 |
clen[1] = lensig*clen[1]; |
98 |
|
112 |
len = ELL_2V_LEN(clen); |
99 |
|
112 |
ELL_3V_SCALE(vec, len, vec); |
100 |
|
112 |
} |
101 |
|
|
|
102 |
|
|
/* |
103 |
|
|
** makeVolume allocates inside nin a new randomly oriented volume of |
104 |
|
|
** size sx-by-sy-sz; the spaceDirection vectors are random but enforced |
105 |
|
|
** to not be too parallel |
106 |
|
|
*/ |
107 |
|
|
static gageContext * |
108 |
|
|
makeVolume(Nrrd *nin, unsigned int sx, unsigned int sy, unsigned int sz, |
109 |
|
|
double coef[POWER_MAX+1][POWER_MAX+1][POWER_MAX+1], |
110 |
|
|
airRandMTState *rng) { |
111 |
|
|
static const char me[]="makeVolume"; |
112 |
|
60 |
double spcOrig[NRRD_SPACE_DIM_MAX], spcVec[3][NRRD_SPACE_DIM_MAX], |
113 |
|
|
angle10, angle20, angle21, ooff[3], aperm=0.6, lexp=0.1, lstdv=0.04, |
114 |
|
|
kparm[NRRD_KERNEL_PARMS_NUM]; |
115 |
|
|
gageContext *gctx; |
116 |
|
|
gagePerVolume *gpvl; |
117 |
|
|
unsigned int xi, yi, zi; |
118 |
|
|
airArray *submop; |
119 |
|
|
int EE; |
120 |
|
|
|
121 |
|
30 |
submop = airMopNew(); |
122 |
|
30 |
randvec(spcVec[0], lexp, lstdv, rng); |
123 |
|
30 |
do { |
124 |
|
32 |
randvec(spcVec[1], lexp, lstdv, rng); |
125 |
|
32 |
angle10 = ell_3v_angle_d(spcVec[1], spcVec[0]); |
126 |
✓✓✓✓
|
95 |
} while (!( AIR_IN_OP(AIR_PI*(1-aperm)/2, angle10, AIR_PI*(1+aperm)/2) )); |
127 |
|
|
do { |
128 |
|
50 |
randvec(spcVec[2], lexp, lstdv, rng); |
129 |
|
50 |
angle20 = ell_3v_angle_d(spcVec[2], spcVec[0]); |
130 |
|
50 |
angle21 = ell_3v_angle_d(spcVec[2], spcVec[1]); |
131 |
✓✓✓✓
✓✓ |
146 |
} while (!( AIR_IN_OP(AIR_PI*(1-aperm)/2, angle20, AIR_PI*(1+aperm)/2) && |
132 |
✓✓ |
76 |
AIR_IN_OP(AIR_PI*(1-aperm)/2, angle21, AIR_PI*(1+aperm)/2) )); |
133 |
|
|
|
134 |
|
|
/* center volume near origin */ |
135 |
|
30 |
airNormalRand_r(ooff + 0, ooff + 1, rng); |
136 |
|
30 |
airNormalRand_r(ooff + 2, NULL, rng); |
137 |
|
30 |
ELL_3V_SET(spcOrig, 0.0, 0.0, 0.0); |
138 |
|
30 |
ELL_3V_SCALE_INCR(spcOrig, -AIR_CAST(double, sx)/2.0 + 2*ooff[0], spcVec[0]); |
139 |
|
30 |
ELL_3V_SCALE_INCR(spcOrig, -AIR_CAST(double, sy)/2.0 + 2*ooff[1], spcVec[1]); |
140 |
|
30 |
ELL_3V_SCALE_INCR(spcOrig, -AIR_CAST(double, sz)/2.0 + 2*ooff[1], spcVec[2]); |
141 |
|
|
|
142 |
|
|
/* allocate data */ |
143 |
✗✓ |
60 |
if (nrrdMaybeAlloc_va(nin, nrrdTypeDouble, 3, |
144 |
|
30 |
AIR_CAST(size_t, sx), |
145 |
|
30 |
AIR_CAST(size_t, sy), |
146 |
|
30 |
AIR_CAST(size_t, sz)) |
147 |
✓✗ |
60 |
|| nrrdSpaceSet(nin, nrrdSpaceRightAnteriorSuperior) |
148 |
✓✗ |
60 |
|| nrrdSpaceOriginSet(nin, spcOrig)) { |
149 |
|
|
biffMovef(PROBE, NRRD, "%s: trouble setting volume", me); |
150 |
|
|
airMopError(submop); return NULL; |
151 |
|
|
} |
152 |
|
30 |
nrrdAxisInfoSet_va(nin, nrrdAxisInfoSpaceDirection, |
153 |
|
|
spcVec[0], |
154 |
|
|
spcVec[1], |
155 |
|
|
spcVec[2]); |
156 |
|
30 |
nrrdAxisInfoSet_va(nin, nrrdAxisInfoCenter, |
157 |
|
|
nrrdCenterCell, |
158 |
|
|
nrrdCenterCell, |
159 |
|
|
nrrdCenterCell); |
160 |
|
|
|
161 |
|
|
/* set data */ |
162 |
|
|
{ |
163 |
|
30 |
double *ddata = AIR_CAST(double *, nin->data); |
164 |
✓✓ |
5236 |
for (zi=0; zi<sz; zi++) { |
165 |
|
2588 |
double pos[3]; |
166 |
✓✓ |
405998 |
for (yi=0; yi<sy; yi++) { |
167 |
✓✓ |
28975754 |
for (xi=0; xi<sx; xi++) { |
168 |
|
14287466 |
ELL_3V_SCALE_ADD2(pos, 1.0, spcOrig, xi, spcVec[0]); |
169 |
|
14287466 |
ELL_3V_SCALE_ADD2(pos, 1.0, pos, yi, spcVec[1]); |
170 |
|
14287466 |
ELL_3V_SCALE_ADD2(pos, 1.0, pos, zi, spcVec[2]); |
171 |
|
14287466 |
ddata[xi + sx*(yi + sy*zi)] = polyeval(coef, 0, 0, 0, pos); |
172 |
|
|
} |
173 |
|
|
} |
174 |
|
2588 |
} |
175 |
|
|
} |
176 |
|
|
|
177 |
|
|
/* create context, using the c4hexic kernel and its derivatives. c4hexic |
178 |
|
|
is can reconstruct cubics without pre-filtering, and its analytical |
179 |
|
|
derivatives are readily available. tmf:n,3,4 can also reconstruct |
180 |
|
|
cubics but its not clear which other kernels are its derivatives */ |
181 |
|
30 |
gctx = gageContextNew(); |
182 |
|
30 |
airMopAdd(submop, gctx, (airMopper)gageContextNix, airMopOnError); |
183 |
|
30 |
gageParmSet(gctx, gageParmRenormalize, AIR_FALSE); |
184 |
|
30 |
gageParmSet(gctx, gageParmCheckIntegrals, AIR_TRUE); |
185 |
|
|
EE = 0; |
186 |
✓✗ |
60 |
if (!EE) EE |= gageKernelSet(gctx, gageKernel00, nrrdKernelC4Hexic, kparm); |
187 |
✓✗ |
60 |
if (!EE) EE |= gageKernelSet(gctx, gageKernel11, nrrdKernelC4HexicD, kparm); |
188 |
✓✗ |
60 |
if (!EE) EE |= gageKernelSet(gctx, gageKernel22, nrrdKernelC4HexicDD, kparm); |
189 |
✓✗ |
60 |
if (!EE) EE |= !(gpvl = gagePerVolumeNew(gctx, nin, gageKindScl)); |
190 |
✓✗ |
60 |
if (!EE) EE |= gagePerVolumeAttach(gctx, gpvl); |
191 |
✓✗ |
60 |
if (!EE) EE |= gageQueryItemOn(gctx, gpvl, gageSclValue); |
192 |
✓✗ |
60 |
if (!EE) EE |= gageQueryItemOn(gctx, gpvl, gageSclGradVec); |
193 |
✓✗ |
60 |
if (!EE) EE |= gageQueryItemOn(gctx, gpvl, gageSclHessian); |
194 |
✓✗ |
60 |
if (!EE) EE |= gageUpdate(gctx); |
195 |
✗✓ |
30 |
if (EE) { |
196 |
|
|
biffMovef(PROBE, GAGE, "%s: trouble setting up gage", me); |
197 |
|
|
airMopError(submop); return NULL; |
198 |
|
|
} |
199 |
|
|
|
200 |
|
30 |
airMopOkay(submop); |
201 |
|
30 |
return gctx; |
202 |
|
30 |
} |
203 |
|
|
|
204 |
|
|
int |
205 |
|
|
main() { |
206 |
|
|
airArray *mop; |
207 |
|
|
char *err; |
208 |
|
|
|
209 |
|
|
airRandMTState *rng; |
210 |
|
|
Nrrd *nin; |
211 |
|
|
unsigned int xi, yi, zi, runNum, runIdx; |
212 |
|
2 |
double coef[POWER_MAX+1][POWER_MAX+1][POWER_MAX+1], pos[3], epsilon; |
213 |
|
|
const double *vmeas, *gmeas, *hmeas; |
214 |
|
|
gageContext *gctx; |
215 |
|
|
|
216 |
|
1 |
mop = airMopNew(); |
217 |
|
|
|
218 |
|
1 |
rng = airRandMTStateNew(429); |
219 |
|
1 |
airMopAdd(mop, rng, (airMopper)airRandMTStateNix, airMopAlways); |
220 |
|
|
|
221 |
|
1 |
nin = nrrdNew(); |
222 |
|
1 |
airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopAlways); |
223 |
|
|
epsilon = 1e-8; |
224 |
|
|
runNum = 30; |
225 |
|
|
|
226 |
✓✓ |
62 |
for (runIdx=0; runIdx<runNum; runIdx++) { |
227 |
|
|
unsigned int sx, sy, sz, probeNum; |
228 |
|
|
double avgDiff; |
229 |
|
|
airArray *submop; |
230 |
|
|
|
231 |
|
30 |
submop = airMopNew(); |
232 |
|
|
/* set random coefficients in polynomial */ |
233 |
✓✓ |
300 |
for (xi=0; xi<=POWER_MAX; xi++) { |
234 |
✓✓ |
1200 |
for (yi=0; yi<=POWER_MAX; yi++) { |
235 |
✓✓ |
4800 |
for (zi=0; zi<=POWER_MAX; zi++) { |
236 |
✓✓✓✓
|
3840 |
if (xi + yi + zi > POWER_MAX) { |
237 |
|
3240 |
coef[xi][yi][zi] = 0.0; |
238 |
|
1320 |
} else { |
239 |
|
600 |
airNormalRand_r(&(coef[xi][yi][zi]), NULL, rng); |
240 |
|
|
} |
241 |
|
|
} |
242 |
|
|
} |
243 |
|
|
} |
244 |
|
|
|
245 |
|
30 |
sx = 20 + airRandInt_r(rng, 120); |
246 |
|
30 |
sy = 20 + airRandInt_r(rng, 120); |
247 |
|
30 |
sz = 20 + airRandInt_r(rng, 120); |
248 |
|
30 |
fprintf(stderr, "%u: %u %u %u: ", runIdx, sx, sy, sz); fflush(stderr); |
249 |
✗✓ |
30 |
if (!(gctx = makeVolume(nin, sx, sy, sz, coef, rng))) { |
250 |
|
|
airMopAdd(mop, err = biffGetDone(PROBE), airFree, airMopAlways); |
251 |
|
|
fprintf(stderr, "trouble creating volume:\n%s", err); |
252 |
|
|
airMopError(submop); airMopError(mop); return 1; |
253 |
|
|
} |
254 |
|
30 |
airMopAdd(submop, gctx, (airMopper)gageContextNix, airMopAlways); |
255 |
|
30 |
vmeas = gageAnswerPointer(gctx, gctx->pvl[0], gageSclValue); |
256 |
|
30 |
gmeas = gageAnswerPointer(gctx, gctx->pvl[0], gageSclGradVec); |
257 |
|
30 |
hmeas = gageAnswerPointer(gctx, gctx->pvl[0], gageSclHessian); |
258 |
|
30 |
ELL_3V_SET(pos, 0, 0, 0); |
259 |
|
30 |
gageProbeSpace(gctx, pos[0], pos[1], pos[2], |
260 |
|
|
AIR_FALSE /* indexSpace */, |
261 |
|
|
AIR_TRUE /* clamp */); |
262 |
|
|
probeNum = 0; |
263 |
|
|
avgDiff = 0.0; |
264 |
|
|
/* take a random walk starting at (0,0,0), making sure that at the |
265 |
|
|
visited positions the gage-measured values, gradients, and |
266 |
|
|
hessians are the same as the analytical ones */ |
267 |
|
30 |
do { |
268 |
|
1135 |
double vanal, ganal[3], hanal[9], vdiff, gdiff[3], hdiff[9], |
269 |
|
|
step[3], stepSize = 0.2; |
270 |
|
1135 |
probeNum++; |
271 |
|
1135 |
vanal = polyeval(coef, 0, 0, 0, pos); |
272 |
|
1135 |
ganal[0] = polyeval(coef, 1, 0, 0, pos); |
273 |
|
1135 |
ganal[1] = polyeval(coef, 0, 1, 0, pos); |
274 |
|
1135 |
ganal[2] = polyeval(coef, 0, 0, 1, pos); |
275 |
|
1135 |
hanal[0] = polyeval(coef, 2, 0, 0, pos); |
276 |
|
1135 |
hanal[1] = polyeval(coef, 1, 1, 0, pos); |
277 |
|
1135 |
hanal[2] = polyeval(coef, 1, 0, 1, pos); |
278 |
|
|
hanal[3] = hanal[1]; |
279 |
|
1135 |
hanal[4] = polyeval(coef, 0, 2, 0, pos); |
280 |
|
1135 |
hanal[5] = polyeval(coef, 0, 1, 1, pos); |
281 |
|
|
hanal[6] = hanal[2]; |
282 |
|
|
hanal[7] = hanal[5]; |
283 |
|
1135 |
hanal[8] = polyeval(coef, 0, 0, 2, pos); |
284 |
|
1135 |
vdiff = fabs(vmeas[0] - vanal); |
285 |
|
1135 |
ELL_3V_SUB(gdiff, gmeas, ganal); |
286 |
|
1135 |
ELL_3M_SUB(hdiff, hmeas, hanal); |
287 |
✓✗✗✓
|
2270 |
if (vdiff > epsilon || |
288 |
✓✗ |
1135 |
ELL_3V_LEN(gdiff) > epsilon || |
289 |
|
1135 |
ELL_3M_FROB(hdiff) > epsilon) { |
290 |
|
|
fprintf(stderr, "at (%g,%g,%g) one diff %.17g %.17g %.17g " |
291 |
|
|
"> eps %.17g\n", |
292 |
|
|
pos[0], pos[1], pos[2], vdiff, |
293 |
|
|
ELL_3V_LEN(gdiff), ELL_3M_FROB(hdiff), epsilon); |
294 |
|
|
airMopError(submop); airMopError(mop); return 1; |
295 |
|
|
} |
296 |
|
1135 |
avgDiff += vdiff + ELL_3V_LEN(gdiff) + ELL_3M_FROB(hdiff); |
297 |
|
1135 |
airNormalRand_r(step + 0, step + 1, rng); |
298 |
|
1135 |
airNormalRand_r(step + 2, NULL, rng); |
299 |
|
1135 |
ELL_3V_SCALE_INCR(pos, stepSize, step); |
300 |
|
1135 |
gageProbeSpace(gctx, pos[0], pos[1], pos[2], |
301 |
|
|
AIR_FALSE /* indexSpace */, |
302 |
|
|
AIR_TRUE /* clamp */); |
303 |
✓✗✓✓
|
3405 |
} while (!gctx->edgeFrac); |
304 |
|
30 |
avgDiff /= probeNum; |
305 |
|
30 |
fprintf(stderr, "%u (%.17g)\n", probeNum, avgDiff); |
306 |
|
30 |
airMopOkay(submop); |
307 |
|
30 |
} |
308 |
|
|
|
309 |
|
1 |
airMopOkay(mop); |
310 |
|
1 |
return 0; |
311 |
|
1 |
} |