Bug Summary

File:src/gage/optimsig.c
Location:line 328, column 5
Description:Potential leak of memory pointed to by 'oscx'

Annotated Source Code

1/*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 2013, 2012, 2011, 2010, 2009 University of Chicago
4 Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5 Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6
7 This library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public License
9 (LGPL) as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
11 The terms of redistributing and/or modifying this software also
12 include exceptions to the LGPL that facilitate static linking.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with this library; if not, write to Free Software Foundation, Inc.,
21 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22*/
23
24#include "gage.h"
25#include "privateGage.h"
26
27/*
28static int debugging = 0;
29static int debugii;
30*/
31
32static airArray *debugReconErrArr = NULL((void*)0);
33static double *debugReconErr = NULL((void*)0);
34static char *debugReconErrName = NULL((void*)0);
35
36/*
37** learned:
38**
39** -- debug high/discontinuous errors at the low sigmas: was because
40** cut-off was insufficient to prevent some discontinuous change in
41** kernel values: increased minimum support in the kernel itself, and
42** now using larger cut-offs.
43**
44** -- also, separately from this problem, you can have minima in the
45** inf error (in imgMeasr) *not* at sample points, apparently simply
46** because of how the hermite interpolation works (but this is
47** troubling)
48**
49** -- do now have a different minimization scheme for allMeasr=Linf,
50** but this may still be a work in progress. Recognizing that this is
51** essentially seeking to find a uniform re-parameterization of
52** something with a hidden non-uniform parameterization, we could
53** probably implement a simple global warping of control points within
54** the implied non-uniform domain.
55*/
56
57/* this limits how big the kernel can get with a single evaluation
58 of nrrdKernelDiscreteGaussian; there are some numerical issues
59 with large kernels that need ironing out */
60#define GOOD_SIGMA_MAX5 5
61
62#define N-1 -1
63
64/*
65** NOTE: The idea for this table originated with Raul San Jose Estepar;
66** GLK recomputed it optimizing for 3D recon, but
67** NOTE: there are probably still be bugs in this; look at the
68** "HEY: bug?" notes below, the same problem occurs elsewhere
69**
70** Basic indexing idea: [sigma max][total # samples][which sample]
71**
72** "sigma max" can't be 0; smallest value is 1
73** ==> index with (sigma max)-1
74** biggest value is GAGE_OPTIMSIG_SIGMA_MAX,
75** ==> biggest index is GAGE_OPTIMSIG_SIGMA_MAX-1
76** ==> allocate for GAGE_OPTIMSIG_SIGMA_MAX
77**
78** "total # samples" can't be 0, or 1, smallest value is 2
79** ==> index with (total # samples)-2
80** biggest value is GAGE_OPTIMSIG_SAMPLES_MAXNUM
81** ==> biggest index is GAGE_OPTIMSIG_SAMPLES_MAXNUM-2
82** ==> allocate for GAGE_OPTIMSIG_SAMPLES_MAXNUM-1
83**
84** "which sample" ranges from 0 to GAGE_OPTIMSIG_SAMPLES_MAXNUM-1
85** ==> allocate for GAGE_OPTIMSIG_SAMPLES_MAXNUM
86*/
87static double
88_optimSigTable[GAGE_OPTIMSIG_SIGMA_MAX11][GAGE_OPTIMSIG_SAMPLES_MAXNUM11-1][GAGE_OPTIMSIG_SAMPLES_MAXNUM11] = {
89 {
90 {0,1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
91 {0,0.5279398,1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
92 {0,0.30728838,0.59967405,1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
93 {0,0.25022203,0.47050092,0.69525677,1,N-1,N-1,N-1,N-1,N-1,N-1},
94 {0,0.17127343,0.39234546,0.56356072,0.75660759,1,N-1,N-1,N-1,N-1,N-1},
95 {0,0.16795139,0.37100673,0.51324213,0.65655005,0.81952846,1,N-1,N-1,N-1,N-1},
96 {0,0.1662873,0.34969759,0.46556041,0.55324608,0.68717259,0.83465695,1,N-1,N-1,N-1},
97 {0,0.12720504,0.22565289,0.28316727,0.44209728,0.58615023,0.75034028,0.87391609,1,N-1,N-1},
98 {0,0.12836272 /* HEY: bug? should be < 0.12720504 */,0.22926401,0.27715567,0.43546647,0.56471503,0.69411868,0.80830419,0.89314467,1,N-1},
99 {0,0.13169055 /* HEY: bug? should be < 0.12720504 */,0.23498112,0.26570156,0.42672107,0.54272485,0.62969965,0.73375762,0.76996493,0.89293921,1}
100 }, {
101 {0,2,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
102 {0,0.75118494,2,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
103 {0,0.55478472,1.1535828,2,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
104 {0,0.49007216,0.8412028,1.308665,2,N-1,N-1,N-1,N-1,N-1,N-1},
105 {0,0.29460263,0.57445061,0.93797231,1.368475,2,N-1,N-1,N-1,N-1,N-1},
106 {0,0.2506085,0.49080029,0.73882496,1.069332,1.4497081,2,N-1,N-1,N-1,N-1},
107 {0,0.18255657,0.42056954,0.62766695,0.87999368,1.1692151,1.5175625,2,N-1,N-1,N-1},
108 {0,0.17582123,0.40522173,0.58696139,0.79624867,1.0485514,1.2950466,1.5977446,2,N-1,N-1},
109 {0,0.17304537,0.39376548,0.56427032,0.75127059,0.96672511,1.187861,1.4141362,1.6921321,2,N-1},
110 {0,0.16970521,0.38116929,0.53575242,0.69498152,0.88430929,1.0844854,1.2899524,1.5211773,1.7645421,2}
111 }, {
112 {0,3,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
113 {0,0.92324787,3,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
114 {0,0.59671402,1.3871731,3,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
115 {0,0.53303385,1.0274624,1.6725048,3,N-1,N-1,N-1,N-1,N-1,N-1},
116 {0,0.47298154,0.79659319,1.2379739,1.8434249,3,N-1,N-1,N-1,N-1,N-1},
117 {0,0.29337707,0.56664073,0.94871783,1.3666322,1.949043,3,N-1,N-1,N-1,N-1},
118 {0,0.25583801,0.52919179,0.78387552,1.1250161,1.516176,2.0632432,3,N-1,N-1,N-1},
119 {0,0.25013804,0.48255014,0.72428173,1.0308567,1.3638159,1.7629964,2.2885511,3,N-1,N-1},
120 {0,0.25038671,0.46448985,0.67336935,0.94502586,1.2324173,1.5780864,1.9883285,2.5002999,3,N-1},
121 {0,0.25034565,0.44725224,0.63590652,0.8669008,1.1130947,1.3942779,1.7180597,2.1408446,2.5466051,3}
122 }, {
123 {0,4,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
124 {0,1.0342592,4,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
125 {0,0.6341188,1.5414433,4,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
126 {0,0.5523203,1.1400089,1.9595566,4,N-1,N-1,N-1,N-1,N-1,N-1},
127 {0,0.51082283,0.91567439,1.4275582,2.2504199,4,N-1,N-1,N-1,N-1,N-1},
128 {0,0.46390373,0.76406777,1.1620381,1.6579833,2.470386,4,N-1,N-1,N-1,N-1},
129 {0,0.29957226,0.58226484,0.90447241,1.318499,1.8011117,2.5972142,4,N-1,N-1,N-1},
130 {0,0.29072434,0.5657317,0.8687849,1.2413157,1.7351674,2.2752147,3.1038468,4,N-1,N-1},
131 {0,0.25000414,0.5027808,0.75375289,1.0744231,1.4267329,1.8665372,2.4665236,3.2203004,4,N-1},
132 {0,0.19010291,0.44269502,0.66081244,0.95829803,1.2627038,1.6005131,2.0043969,2.6440792,3.2979164,4}
133 }, {
134 {0,5,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
135 {0,1.1176668,5,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
136 {0,0.66791451,1.6688319,5,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
137 {0,0.56513244,1.2151262,2.2046661,5,N-1,N-1,N-1,N-1,N-1,N-1},
138 {0,0.51955444,0.96157616,1.5293243,2.5639,5,N-1,N-1,N-1,N-1,N-1},
139 {0,0.50639188,0.83235806,1.2596023,1.8475783,2.8751452,5,N-1,N-1,N-1,N-1},
140 {0,0.30821687,0.60048282,1.0057166,1.4351804,2.0372179,3.0747592,5,N-1,N-1,N-1},
141 {0,0.28437388,0.560866,0.92278755,1.3049414,1.7620444,2.4607313,3.5198457,5,N-1,N-1},
142 {0,0.26883101,0.53947717,0.84076571,1.1986721,1.6077875,2.165575,2.9591467,3.931181,5,N-1},
143 {0,0.25029126,0.50162876,0.75587535,1.0861237,1.4452776,1.8865763,2.5002809,3.2476835,4.0337272,5}
144 }, {
145 {0,6,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
146 {0,1.185726,6,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
147 {0,0.69637311,1.7772807,6,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
148 {0,0.57470578,1.2709187,2.4227901,6,N-1,N-1,N-1,N-1,N-1,N-1},
149 {0,0.52996641,1.0128419,1.632214,2.8718762,6,N-1,N-1,N-1,N-1,N-1},
150 {0,0.50426048,0.87729794,1.3428378,2.0053113,3.2981832,6,N-1,N-1,N-1,N-1},
151 {0,0.46658435,0.76617205,1.1726109,1.6950468,2.5514688,4.1463666,6,N-1,N-1,N-1},
152 {0,0.50030917,0.78596908,1.1486269,1.5887094,2.2150676,3.2805684,4.4828262,6,N-1,N-1},
153 {0,0.27919531,0.56878412,0.88591647,1.2631332,1.7201432,2.3851209,3.392889,4.6255312,6,N-1},
154 {0,0.25088972,0.50369233,0.78494686,1.1030188,1.482311,1.9812444,2.6906328,3.734978,4.7532525,6}
155 }, {
156 {0,7,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
157 {0,1.2437892,7,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
158 {0,0.72099203,1.8771845,7,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
159 {0,0.58251196,1.3139123,2.6157444,7,N-1,N-1,N-1,N-1,N-1,N-1},
160 {0,0.5371021,1.0473768,1.7166929,3.1448426,7,N-1,N-1,N-1,N-1,N-1},
161 {0,0.51312029,0.92989367,1.4221185,2.2125893,3.6739931,7,N-1,N-1,N-1,N-1},
162 {0,0.50083971,0.84841007,1.2561073,1.8532455,2.8668625,4.7535434,7,N-1,N-1,N-1},
163 {0,0.3375614,0.63945627,1.0301709,1.4884938,2.073761,3.1614799,5.0744987,7,N-1,N-1},
164 {0,0.29428458,0.58668923,0.93714356,1.3736334,1.8300356,2.6405344,3.9042048,5.3097196,7,N-1},
165 {0,0.25234449,0.52068585,0.79422623,1.1273863,1.5991755,2.1453006,2.8984315,4.1899557,5.4597921,7}
166 }, {
167 {0,8,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
168 {0,1.2942501,8,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
169 {0,0.74332041,1.9693407,8,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
170 {0,0.58823597,1.3483386,2.7880962,8,N-1,N-1,N-1,N-1,N-1,N-1},
171 {0,0.56661958,1.2263036,1.9593971,3.6037345,8,N-1,N-1,N-1,N-1,N-1},
172 {0,0.52106231,0.97026396,1.486012,2.3670862,4.1632919,8,N-1,N-1,N-1,N-1},
173 {0,0.50727636,0.86810225,1.3293955,2.0115428,3.1358411,5.3943086,8,N-1,N-1,N-1},
174 {0,0.47202346,0.77812189,1.1608884,1.6648751,2.4694417,3.9094045,5.7665443,8,N-1,N-1},
175 {0,0.37446901,0.66116196,1.038642,1.4625595,2.0528309,2.9814169,4.4429126,5.9815373,8,N-1},
176 {0,0.26310974,0.54373014,0.84282249,1.2090484,1.6551158,2.3275802,3.3196113,4.7216973,6.1578932,8}
177 }, {
178 {0,9,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
179 {0,1.3413963,9,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
180 {0,0.76222414,2.0542119,9,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
181 {0,0.59559792,1.3777219,2.946173,9,N-1,N-1,N-1,N-1,N-1,N-1},
182 {0,0.56240517,1.1527119,1.9145473,3.6841569,9,N-1,N-1,N-1,N-1,N-1},
183 {0,0.52387071,0.98832464,1.5376476,2.5417714,4.4669261,9,N-1,N-1,N-1,N-1},
184 {0,0.50359035,0.87327009,1.3558764,2.0646384,3.3180211,5.9420524,9,N-1,N-1,N-1},
185 {0,0.50140077,0.83020425,1.256588,1.7709454,2.7100441,4.4434023,6.3934889,9,N-1,N-1},
186 {0,0.36521655,0.65757704,1.0627806,1.5081434,2.1497617,3.1920822,4.870122,6.6418982,9,N-1},
187 {0,0.31160679,0.59032226,0.94745982,1.3620865,1.8115216,2.6007423,3.8324564,5.2064519,6.8468728,9}
188 }, {
189 {0,10,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
190 {0,1.3838946,10,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
191 {0,0.77946955,2.1342247,10,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
192 {0,0.60070014,1.4040204,3.0944126,10,N-1,N-1,N-1,N-1,N-1,N-1},
193 {0,0.55609542,1.1508646,1.9495349,3.9375696,10,N-1,N-1,N-1,N-1,N-1},
194 {0,0.5350194,1.031119,1.6607633,2.8520992,5.4718146,10,N-1,N-1,N-1,N-1},
195 {0,0.5083549,0.90783268,1.4059756,2.1796026,3.571064,6.5497985,10,N-1,N-1,N-1},
196 {0,0.50199872,0.85233968,1.2647815,1.8777326,2.8592849,4.7821364,7.0110598,10,N-1,N-1},
197 {0,0.46663594,0.75212663,1.1302133,1.6134665,2.3560972,3.6558499,5.3189116,7.2945781,10,N-1},
198 {0,0.3789258,0.64023608,1.0374272,1.4685256,2.0717783,3.0241971,4.2591534,5.6669927,7.5286098,10}
199 }, {
200 {0,11,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
201 {0,1.4234025,11,N-1,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
202 {0,0.79513794,2.2098076,11,N-1,N-1,N-1,N-1,N-1,N-1,N-1},
203 {0,0.60728961,1.4287171,3.2358651,11,N-1,N-1,N-1,N-1,N-1,N-1},
204 {0,0.55890071,1.165283,2.0149148,4.1530919,11,N-1,N-1,N-1,N-1,N-1},
205 {0,0.55071467,1.0660659,1.7177736,3.0094495,6.0395317,11,N-1,N-1,N-1,N-1},
206 {0,0.5066433,0.89661205,1.4050072,2.2117786,3.7080047,7.0954437,11,N-1,N-1,N-1},
207 {0,0.50242329,0.86727452,1.3264461,1.9118301,2.9509099,5.1184769,7.624764,11,N-1,N-1},
208 {0,0.47785854,0.78873962,1.1769236,1.6880652,2.4978926,4.0288033,5.7288432,7.9420485,11,N-1},
209 {0,0.50532979,0.79486167,1.1706896,1.6148115,2.2648265,3.3499777,4.5595574,6.116312,8.2049971,11}
210 }
211};
212
213/*
214** this is only for retreiving part of the table above
215*/
216int
217gageOptimSigSet(double *scale, unsigned int num, unsigned int sigmaMax) {
218 static const char me[]="gageOptimSigSet";
219 unsigned int si;
220
221 if (!scale) {
222 biffAddf(GAGEgageBiffKey, "%s: got NULL pointer", me);
223 return 1;
224 }
225 if (!AIR_IN_CL(2, num, GAGE_OPTIMSIG_SAMPLES_MAXNUM)((2) <= (num) && (num) <= (11))) {
226 biffAddf(GAGEgageBiffKey,
227 "%s: requested # sigma samples %u not in known range [2,%u]",
228 me, num, GAGE_OPTIMSIG_SAMPLES_MAXNUM11);
229 return 1;
230 }
231 if (!AIR_IN_CL(1, sigmaMax, GAGE_OPTIMSIG_SIGMA_MAX)((1) <= (sigmaMax) && (sigmaMax) <= (11))) {
232 biffAddf(GAGEgageBiffKey, "%s: requested sigma max %u not in known range [1,%u]",
233 me, sigmaMax, GAGE_OPTIMSIG_SIGMA_MAX11);
234 return 1;
235 }
236
237 for (si=0; si<num; si++) {
238 scale[si] = _optimSigTable[sigmaMax-1][num-2][si];
239 }
240 return 0;
241}
242
243/* ------- from here down is the stuff for computing the table ------ */
244
245/* rho is a stand-in for tau - and something that will likely change
246 based on the findings from using this code; the idea is that it
247 reflects the needed density of samples for optimal scale-space
248 reconstruction. In order to be used for the internal workings of
249 the sigma optimization, its important that the conversion between
250 sigma and rho be accurately invertible. */
251
252/*
253**
254** This is a decent approximation of tau(sigma), made of a slightly
255** scaled version of the taylor expansion of tau(sigma=0) which meets
256** up with the large-scale approximation of tau from Lindeberg.
257** However, because its so flat at sigma=0, its not really invertible
258** there, so its a poor basis for computations that are parameterized
259** by rho. Keeping it around for reference.
260
261static double
262_RhoOfSig(double sig) {
263 double rho;
264
265 if (sig < 1.05189095) {
266 rho = sig*sig*(0.2775733212544225 + 0.13078298856958057*sig*sig);
267 } else {
268 double tee;
269 tee = sig*sig;
270 rho = 0.53653222368715360118 + log(tee)/2.0 + log(1.0 - 1.0/(8.0*tee));
271 }
272 return rho;
273}
274
275static double
276_SigOfRho(double rho) {
277 double sig;
278
279 if (rho < 0.46724360022171363) {
280 sig = 0.00033978812426865065 *
281 sqrt(-9.191366355042886e6 + 245.3752559286824 *
282 sqrt(1.403132301e9 + 9.526961876920057e9*rho));
283 } else {
284 double ee, tee;
285 ee = exp(2.0*rho);
286 tee = 0.0063325739776461107152*(27.0*ee + 2*AIR_PI*AIR_PI
287 + 3.0*sqrt(81.0*ee*ee
288 + 12*ee*AIR_PI*AIR_PI));
289 sig = sqrt(tee);
290 }
291 return sig;
292}
293*/
294
295static double
296_RhoOfSig(double sig) {
297
298 return log(sig + 1);
299}
300
301static double
302_SigOfRho(double rho) {
303
304 return exp(rho)-1;
305}
306
307
308/*
309** allocates context, with error checking
310** does use biff
311*/
312gageOptimSigContext *gageOptimSigContextNew(unsigned int dim,
313 unsigned int sampleNumMax,
314 unsigned int trueImgNum,
315 double sigmaMin, double sigmaMax,
316 double cutoff) {
317 static const char me[]="gageOptimSigContextNew";
318 gageOptimSigContext *oscx;
319 unsigned int support, ii;
320 double kparm[2];
321
322 oscx = AIR_CALLOC(1, gageOptimSigContext)(gageOptimSigContext*)(calloc((1), sizeof(gageOptimSigContext
)));
1
Within the expansion of the macro 'AIR_CALLOC':
→
a
Memory is allocated
323 if (!oscx) {
2
←
Assuming 'oscx' is non-null
→
3
←
Taking false branch
→
324 biffAddf(GAGEgageBiffKey, "%s: couldn't allocate context", me);
325 return NULL((void*)0);
326 }
327 if (!AIR_IN_CL(1, dim, 3)((1) <= (dim) && (dim) <= (3))) {
4
←
Taking true branch
→
328 biffAddf(GAGEgageBiffKey, "%s: dim %u not 1, 2, or 3", me, dim);
5
←
Potential leak of memory pointed to by 'oscx'
329 return NULL((void*)0);
330 }
331 if (!(sampleNumMax >= 3)) {
332 biffAddf(GAGEgageBiffKey, "%s: sampleNumMax %u not >= 3", me, sampleNumMax);
333 return NULL((void*)0);
334 }
335 if (!(trueImgNum >= 3)) {
336 biffAddf(GAGEgageBiffKey, "%s: trueImgNum %u not >= 3", me, trueImgNum);
337 return NULL((void*)0);
338 }
339 if (!(sigmaMin >= 0 && sigmaMax > sigmaMin && cutoff > 0)) {
340 biffAddf(GAGEgageBiffKey, "%s: sigmaMin %g, sigmaMax %g, cutoff %g not valid", me,
341 sigmaMin, sigmaMax, cutoff);
342 return NULL((void*)0);
343 }
344
345 oscx->dim = dim;
346 oscx->sampleNumMax = sampleNumMax;
347 oscx->trueImgNum = trueImgNum;
348 oscx->sigmaRange[0] = sigmaMin;
349 oscx->sigmaRange[1] = sigmaMax;
350 oscx->cutoff = cutoff;
351
352 /* will be set later */
353 oscx->kssSpec = NULL((void*)0);
354 oscx->sampleNum = 0;
355 oscx->maxIter = 0;
356 oscx->imgMeasr = nrrdMeasureUnknown;
357 oscx->allMeasr = nrrdMeasureUnknown;
358 oscx->convEps = AIR_NAN(airFloatQNaN.f);
359
360 /* allocate internal buffers based on arguments */
361 kparm[0] = oscx->sigmaRange[1];
362 kparm[1] = oscx->cutoff;
363 support = AIR_ROUNDUP(nrrdKernelDiscreteGaussian->support(kparm))((int)(floor((nrrdKernelDiscreteGaussian->support(kparm))+
0.5)));
364 oscx->sx = 2*support - 1;
365 oscx->sy = dim >= 2 ? 2*support - 1 : 1;
366 oscx->sz = dim >= 3 ? 2*support - 1 : 1;
367 /*
368 fprintf(stderr, "%s: max sigma = %g, cutoff %g ==> support=%u, "
369 "3D vol size=%u x %u x %u\n", me,
370 sigmaMax, cutoff, support, oscx->sx, oscx->sy, oscx->sz);
371 */
372 oscx->nerr = nrrdNew();
373 oscx->ninterp = nrrdNew();
374 oscx->ndiff = nrrdNew();
375 if (nrrdMaybeAlloc_va(oscx->nerr, nrrdTypeDouble, 1,
376 AIR_CAST(size_t, oscx->trueImgNum)((size_t)(oscx->trueImgNum)))
377 || nrrdMaybeAlloc_va(oscx->ninterp, nrrdTypeDouble, 3,
378 AIR_CAST(size_t, oscx->sx)((size_t)(oscx->sx)),
379 AIR_CAST(size_t, oscx->sy)((size_t)(oscx->sy)),
380 AIR_CAST(size_t, oscx->sz)((size_t)(oscx->sz)))
381 || nrrdMaybeAlloc_va(oscx->ndiff, nrrdTypeDouble, 3,
382 AIR_CAST(size_t, oscx->sx)((size_t)(oscx->sx)),
383 AIR_CAST(size_t, oscx->sy)((size_t)(oscx->sy)),
384 AIR_CAST(size_t, oscx->sz)((size_t)(oscx->sz)))) {
385 biffMovef(GAGEgageBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate buffers", me);
386 return NULL((void*)0);
387 }
388 nrrdAxisInfoSet_va(oscx->ninterp, nrrdAxisInfoSpacing,
389 1.0, 1.0, 1.0);
390 nrrdAxisInfoSet_va(oscx->ndiff, nrrdAxisInfoSpacing,
391 1.0, 1.0, 1.0);
392 oscx->rhoRange[0] = _RhoOfSig(oscx->sigmaRange[0]);
393 oscx->rhoRange[1] = _RhoOfSig(oscx->sigmaRange[1]);
394 nrrdAxisInfoSet_va(oscx->nerr, nrrdAxisInfoMin,
395 oscx->rhoRange[0]);
396 nrrdAxisInfoSet_va(oscx->nerr, nrrdAxisInfoMax,
397 oscx->rhoRange[1]);
398
399 fprintf(stderr__stderrp, "!%s: sigma [%g,%g] -> rho [%g,%g]\n", me,
400 oscx->sigmaRange[0], oscx->sigmaRange[1],
401 oscx->rhoRange[0], oscx->rhoRange[1]);
402 fprintf(stderr__stderrp, "!%s: rho %g -- %g\n", me,
403 oscx->rhoRange[0], oscx->rhoRange[1]);
404 for (ii=0; ii<oscx->trueImgNum; ii++) {
405 double rr, ss, rc, eps=1e-13;
406 rr = AIR_AFFINE(0, ii, oscx->trueImgNum-1,( ((double)(oscx->rhoRange[1])-(oscx->rhoRange[0]))*((double
)(ii)-(0)) / ((double)(oscx->trueImgNum-1)-(0)) + (oscx->
rhoRange[0]))
407 oscx->rhoRange[0], oscx->rhoRange[1])( ((double)(oscx->rhoRange[1])-(oscx->rhoRange[0]))*((double
)(ii)-(0)) / ((double)(oscx->trueImgNum-1)-(0)) + (oscx->
rhoRange[0]));
408 ss = _SigOfRho(rr);
409 rc = _RhoOfSig(ss);
410 /*
411 fprintf(stderr, "!%s: (%u) rho %.17g -> sig %.17g -> rho %.17g (%.17g)\n",
412 me, ii, rr, ss, rc, AIR_ABS(rr - rc)/(rr + eps/2));
413 */
414 if ( AIR_ABS(rr - rc)((rr - rc) > 0.0f ? (rr - rc) : -(rr - rc))/(rr + eps) > eps ) {
415 biffAddf(GAGEgageBiffKey, "%s: rho %g -> sig %g -> rho %g has error %g > %g; "
416 "_SigOfRho() and _RhoOfSig() not invertible",
417 me, rr, ss, rc, AIR_ABS(rr - rc)((rr - rc) > 0.0f ? (rr - rc) : -(rr - rc))/(rr + eps/2), eps);
418 return NULL((void*)0);
419 }
420 }
421
422 oscx->kloc = AIR_CALLOC(oscx->sx, double)(double*)(calloc((oscx->sx), sizeof(double)));
423 oscx->kern = AIR_CALLOC(oscx->sx, double)(double*)(calloc((oscx->sx), sizeof(double)));
424 oscx->ktmp1 = AIR_CALLOC(oscx->sx, double)(double*)(calloc((oscx->sx), sizeof(double)));
425 oscx->ktmp2 = AIR_CALLOC(oscx->sx, double)(double*)(calloc((oscx->sx), sizeof(double)));
426 if (!(oscx->kloc && oscx->kern && oscx->ktmp1 && oscx->ktmp2)) {
427 biffAddf(GAGEgageBiffKey, "%s: couldn't allocate kernel buffers", me);
428 return NULL((void*)0);
429 }
430 for (ii=0; ii<oscx->sx; ii++) {
431 oscx->kloc[ii] = AIR_CAST(double, ii)((double)(ii)) - ((oscx->sx + 1)/2 - 1);
432 }
433 oscx->kone[0] = 1.0;
434
435 oscx->gctx = NULL((void*)0);
436 oscx->pvlBase = NULL((void*)0);
437 oscx->pvlSS = AIR_CALLOC(oscx->sampleNumMax, gagePerVolume *)(gagePerVolume **)(calloc((oscx->sampleNumMax), sizeof(gagePerVolume
*)));
438 oscx->nsampleImg = AIR_CALLOC(oscx->sampleNumMax, Nrrd *)(Nrrd **)(calloc((oscx->sampleNumMax), sizeof(Nrrd *)));
439 oscx->sampleSigma = AIR_CALLOC(oscx->sampleNumMax, double)(double*)(calloc((oscx->sampleNumMax), sizeof(double)));
440 oscx->sampleRho = AIR_CALLOC(oscx->sampleNumMax, double)(double*)(calloc((oscx->sampleNumMax), sizeof(double)));
441 oscx->sampleTmp = AIR_CALLOC(oscx->sampleNumMax, double)(double*)(calloc((oscx->sampleNumMax), sizeof(double)));
442 oscx->sampleErrMax = AIR_CALLOC(oscx->sampleNumMax, double)(double*)(calloc((oscx->sampleNumMax), sizeof(double)));
443 oscx->step = AIR_CALLOC(oscx->sampleNumMax, double)(double*)(calloc((oscx->sampleNumMax), sizeof(double)));
444 if (!(oscx->pvlSS && oscx->nsampleImg
445 && oscx->sampleSigma && oscx->sampleRho
446 && oscx->sampleTmp && oscx->sampleErrMax && oscx->step)) {
447 biffAddf(GAGEgageBiffKey, "%s: couldn't allocate per-sample arrays", me);
448 return NULL((void*)0);
449 }
450 for (ii=0; ii<oscx->sampleNumMax; ii++) {
451 oscx->nsampleImg[ii] = nrrdNew();
452 if (nrrdMaybeAlloc_va(oscx->nsampleImg[ii], nrrdTypeDouble, 3,
453 AIR_CAST(size_t, oscx->sx)((size_t)(oscx->sx)),
454 AIR_CAST(size_t, oscx->sy)((size_t)(oscx->sy)),
455 AIR_CAST(size_t, oscx->sz)((size_t)(oscx->sz)))) {
456 biffMovef(GAGEgageBiffKey, NRRDnrrdBiffKey, "%s: couldn't allocate vol[%u]", me, ii);
457 return NULL((void*)0);
458 }
459 nrrdAxisInfoSet_va(oscx->nsampleImg[ii], nrrdAxisInfoSpacing,
460 1.0, 1.0, 1.0);
461 }
462
463 /* implementation not started
464 oscx->nsampleHist = nrrdNew();
465 */
466
467 return oscx;
468}
469
470gageOptimSigContext *
471gageOptimSigContextNix(gageOptimSigContext *oscx) {
472
473 if (oscx) {
474 unsigned int si;
475 nrrdKernelSpecNix(oscx->kssSpec);
476 nrrdNuke(oscx->nerr);
477 nrrdNuke(oscx->ninterp);
478 nrrdNuke(oscx->ndiff);
479 airFree(oscx->kloc);
480 airFree(oscx->kern);
481 airFree(oscx->ktmp1);
482 airFree(oscx->ktmp2);
483 gageContextNix(oscx->gctx);
484 /* airFree(oscx->pvlSS); needed? */
485 for (si=0; si<oscx->sampleNumMax; si++) {
486 nrrdNuke(oscx->nsampleImg[si]);
487 }
488 airFree(oscx->nsampleImg);
489 airFree(oscx->sampleSigma);
490 airFree(oscx->sampleRho);
491 airFree(oscx->sampleTmp);
492 airFree(oscx->sampleErrMax);
493 airFree(oscx->step);
494 /* nrrdNuke(oscx->nsampleHist); */
495 airFree(oscx);
496 }
497 return NULL((void*)0);
498}
499
500static int
501_volInterp(Nrrd *ninterp, double rho, gageOptimSigContext *oscx) {
502 static const char me[]="_volInterp";
503 double *interp, scaleIdx, sigma;
504 const double *answer;
505 unsigned int xi, yi, zi;
506 int outside;
507
508 /*
509 debugging = rho > 1.197;
510 gageParmSet(oscx->gctx, gageParmVerbose, 2*debugging);
511 */
512 sigma = _SigOfRho(rho);
513 scaleIdx = gageStackWtoI(oscx->gctx, sigma, &outside);
514 /* Because of limited numerical precision, _SigOfRho(rhoRange[1])
515 can end up "outside" stack, which should really be a bug.
516 However, since the use of gage is pretty straight-forward here,
517 we're okay with ignoring the "outside" here, and also clamping
518 the probe below */
519 answer = gageAnswerPointer(oscx->gctx, oscx->pvlBase, gageSclValue);
520 interp = AIR_CAST(double *, ninterp->data)((double *)(ninterp->data));
521 for (zi=0; zi<oscx->sz; zi++) {
522 for (yi=0; yi<oscx->sy; yi++) {
523 for (xi=0; xi<oscx->sx; xi++) {
524 if (gageStackProbeSpace(oscx->gctx, xi, yi, zi, scaleIdx,
525 AIR_TRUE1 /* index space */,
526 AIR_TRUE1 /* clamping */)) {
527 biffAddf(GAGEgageBiffKey, "%s: probe error at (%u,%u,%u,%.17g): %s (%d)", me,
528 xi, yi, zi, scaleIdx,
529 oscx->gctx->errStr, oscx->gctx->errNum);
530 return 1;
531 }
532 interp[xi + oscx->sx*(yi + oscx->sy*zi)] = answer[0];
533 }
534 }
535 }
536 /*
537 gageParmSet(oscx->gctx, gageParmVerbose, 0);
538 */
539 /*
540 if (debugging) {
541 char fname[AIR_STRLEN_SMALL];
542 sprintf(fname, "interp-%04u.nrrd", debugii);
543 nrrdSave(fname, ninterp, NULL);
544 }
545 */
546 return 0;
547}
548
549static void
550_kernset(double **kzP, double **kyP, double **kxP,
551 gageOptimSigContext *oscx,
552 double rho) {
553 NrrdKernel *dg;
554 double sig, kparm[NRRD_KERNEL_PARMS_NUM8],
555 *kloc, *kern, *ktmp1, *ktmp2;
556 unsigned int ki, kj, kk, sx;
557
558 kern = oscx->kern;
559 kloc = oscx->kloc;
560 ktmp1 = oscx->ktmp1;
561 ktmp2 = oscx->ktmp2;
562 sx = oscx->sx;
563 sig = _SigOfRho(rho);
564 dg = nrrdKernelDiscreteGaussian;
565 kparm[1] = oscx->cutoff;
566 if (sig < GOOD_SIGMA_MAX5) {
567 /* for small sigma, can evaluate directly into kern */
568 kparm[0] = sig;
569 dg->evalN_d(kern, kloc, sx, kparm);
570 } else {
571 double timeleft, tdelta;
572 unsigned int rx;
573 rx = (sx + 1)/2 - 1;
574 /* we have to iteratively blur */
575 kparm[0] = GOOD_SIGMA_MAX5;
576 dg->evalN_d(kern, kloc, sx, kparm);
577 timeleft = sig*sig - GOOD_SIGMA_MAX5*GOOD_SIGMA_MAX5;
578 do {
579 tdelta = AIR_MIN( GOOD_SIGMA_MAX*GOOD_SIGMA_MAX, timeleft )((5*5) < (timeleft) ? (5*5) : (timeleft));
580 kparm[0] = sqrt(tdelta);
581 dg->evalN_d(ktmp1, kloc, sx, kparm);
582 for (ki=0; ki<sx; ki++) {
583 double csum = 0.0;
584 for (kj=0; kj<sx; kj++) {
585 kk = ki - kj + rx;
586 if (kk < sx) {
587 csum += kern[kk]*ktmp1[kj];
588 }
589 }
590 ktmp2[ki] = csum;
591 }
592 for (ki=0; ki<sx; ki++) {
593 kern[ki] = ktmp2[ki];
594 }
595 timeleft -= tdelta;
596 } while (timeleft);
597 }
598 *kzP = oscx->dim >= 3 ? kern : oscx->kone;
599 *kyP = oscx->dim >= 2 ? kern : oscx->kone;
600 *kxP = kern;
601 return;
602}
603
604/*
605** sets one of the sampleImg, to be used as a sample in scale-space interp
606*/
607static void
608_sampleSet(gageOptimSigContext *oscx, unsigned int si, double rho) {
609 double *vol, *kz, *ky, *kx;
610 unsigned int ii, xi, yi, zi;
611
612 oscx->sampleSigma[si] = _SigOfRho(rho);
613 oscx->sampleRho[si] = rho;
614 vol = AIR_CAST(double *, oscx->nsampleImg[si]->data)((double *)(oscx->nsampleImg[si]->data));
615 _kernset(&kz, &ky, &kx, oscx, rho);
616 ii = 0;
617 for (zi=0; zi<oscx->sz; zi++) {
618 for (yi=0; yi<oscx->sy; yi++) {
619 for (xi=0; xi<oscx->sx; xi++) {
620 vol[ii] = kz[zi]*ky[yi]*kx[xi];
621 ii++;
622 }
623 }
624 }
625 if (oscx->gctx) {
626 /* the gage stack needs to know new scale pos */
627 oscx->gctx->stackPos[si] = oscx->sampleSigma[si];
628 /* HEY: GLK forgets why this is needed,
629 but remembers it was a tricky bug to find */
630 gagePointReset(&(oscx->gctx->point));
631 }
632 return;
633}
634
635static int
636_errSingle(double *retP, gageOptimSigContext *oscx, double rho) {
637 static const char me[]="_errSingle";
638 double *diff, *kx, *ky, *kz;
639 const double *interp;
640 unsigned int ii, xi, yi, zi;
641
642 if (_volInterp(oscx->ninterp, rho, oscx)) {
643 biffAddf(GAGEgageBiffKey, "%s: trouble at rho %.17g\n", me, rho);
644 return 1;
645 }
646 /*
647 if (debugging) {
648 char fname[AIR_STRLEN_SMALL];
649 sprintf(fname, "interp-%04u.nrrd", debugii);
650 nrrdSave(fname, oscx->ninterp, NULL);
651 }
652 */
653 interp = AIR_CAST(const double *, oscx->ninterp->data)((const double *)(oscx->ninterp->data));
654 diff = AIR_CAST(double *, oscx->ndiff->data)((double *)(oscx->ndiff->data));
655 _kernset(&kz, &ky, &kx, oscx, rho);
656 ii = 0;
657 for (zi=0; zi<oscx->sz; zi++) {
658 for (yi=0; yi<oscx->sy; yi++) {
659 for (xi=0; xi<oscx->sx; xi++) {
660 double tru;
661 tru = kz[zi]*ky[yi]*kx[xi];
662 diff[ii] = interp[ii] - tru;
663 if (debugReconErrArr) {
664 unsigned int idx = airArrayLenIncr(debugReconErrArr, 2);
665 debugReconErr[idx + 0] = tru;
666 debugReconErr[idx + 1] = interp[ii];
667 }
668 ii++;
669 }
670 }
671 }
672 nrrdMeasureLine[oscx->imgMeasr](retP, nrrdTypeDouble,
673 diff, nrrdTypeDouble,
674 ii /* HEY: cleverness */,
675 AIR_NAN(airFloatQNaN.f), AIR_NAN(airFloatQNaN.f));
676 return 0;
677}
678
679static int
680_errTotal(double *retP, gageOptimSigContext *oscx) {
681 static const char me[]="_errTotal";
682 unsigned int ii;
683 double *err;
684
685 err = AIR_CAST(double *, oscx->nerr->data)((double *)(oscx->nerr->data));
686 for (ii=0; ii<oscx->trueImgNum; ii++) {
687 /* debugii = ii; */
688 if (_errSingle(err + ii, oscx, AIR_AFFINE(0, ii, oscx->trueImgNum-1,( ((double)(oscx->rhoRange[1])-(oscx->rhoRange[0]))*((double
)(ii)-(0)) / ((double)(oscx->trueImgNum-1)-(0)) + (oscx->
rhoRange[0]))
689 oscx->rhoRange[0],( ((double)(oscx->rhoRange[1])-(oscx->rhoRange[0]))*((double
)(ii)-(0)) / ((double)(oscx->trueImgNum-1)-(0)) + (oscx->
rhoRange[0]))
690 oscx->rhoRange[1])( ((double)(oscx->rhoRange[1])-(oscx->rhoRange[0]))*((double
)(ii)-(0)) / ((double)(oscx->trueImgNum-1)-(0)) + (oscx->
rhoRange[0])))) {
691 biffAddf(GAGEgageBiffKey, "%s: trouble at ii %u", me, ii);
692 return 1;
693 }
694 }
695 nrrdMeasureLine[oscx->allMeasr](retP, nrrdTypeDouble,
696 err, nrrdTypeDouble,
697 oscx->trueImgNum,
698 AIR_NAN(airFloatQNaN.f), AIR_NAN(airFloatQNaN.f));
699 /*
700 if (debugging) {
701 char fname[AIR_STRLEN_SMALL];
702 unsigned int ni;
703 for (ni=0; ni<oscx->sampleNum; ni++) {
704 sprintf(fname, "sample-%04u.nrrd", ni);
705 nrrdSave(fname, oscx->nsampleImg[ni], NULL);
706 }
707 }
708 */
709 if (0) {
710 static unsigned int call=0;
711 char fname[AIR_STRLEN_SMALL(128+1)];
712 sprintf(fname, "err-%04u.nrrd", call)__builtin___sprintf_chk (fname, 0, __builtin_object_size (fname
, 2 > 1 ? 1 : 0), "err-%04u.nrrd", call);
713 nrrdSave(fname, oscx->nerr, NULL((void*)0));
714 call++;
715 }
716 return 0;
717}
718
719static int
720_errTotalLinf(double *retP, gageOptimSigContext *oscx,
721 unsigned int mmIdx[2], double mmErr[2]) {
722 static const char me[]="_errTotalLinf";
723 unsigned int ii, pi;
724 double *err, *sem, *rr, rho, sig, pid;
725 int outside;
726
727 err = AIR_CAST(double *, oscx->nerr->data)((double *)(oscx->nerr->data));
728 sem = oscx->sampleErrMax;
729 rr = oscx->rhoRange;
730 for (pi=0; pi<=oscx->sampleNum-2; pi++) {
731 sem[pi] = 0;
732 }
733 /* NOTE: we don't bother with last "true image": it will always be a
734 low error, and not meaningfully associated with a gap */
735 for (ii=0; ii<oscx->trueImgNum-1; ii++) {
736 /* debugii = ii; */
737 rho = AIR_AFFINE(0, ii, oscx->trueImgNum-1, rr[0], rr[1])( ((double)(rr[1])-(rr[0]))*((double)(ii)-(0)) / ((double)(oscx
->trueImgNum-1)-(0)) + (rr[0]));
738 if (_errSingle(err + ii, oscx, rho)) {
739 biffAddf(GAGEgageBiffKey, "%s: trouble at ii %u", me, ii);
740 return 1;
741 }
742 sig = _SigOfRho(rho);
743 pid = gageStackWtoI(oscx->gctx, sig, &outside);
744 pi = AIR_CAST(unsigned int, pid)((unsigned int)(pid));
745 if (outside || !(pi <= oscx->sampleNum-2)) {
746 biffAddf(GAGEgageBiffKey, "%s: ii %u -> rho %g -> sig %g -> pi %u-> OUTSIDE",
747 me, ii, rho, sig, pi);
748 return 1;
749 }
750 sem[pi] = AIR_MAX(sem[pi], err[ii])((sem[pi]) > (err[ii]) ? (sem[pi]) : (err[ii]));
751 }
752 mmIdx[0] = mmIdx[1] = 0;
753 mmErr[0] = mmErr[1] = sem[0];
754 for (pi=1; pi<=oscx->sampleNum-2; pi++) {
755 if (sem[pi] < mmErr[0]) {
756 mmIdx[0] = pi;
757 mmErr[0] = sem[pi];
758 }
759 if (sem[pi] > mmErr[1]) {
760 mmIdx[1] = pi;
761 mmErr[1] = sem[pi];
762 }
763 }
764 /* returned error invented, not really Linf, but minimizing this
765 does the right thing */
766 *retP = 1000*oscx->sampleNum*(mmErr[1] - mmErr[0])/(rr[1] - rr[0]);
767 if (0) {
768 static unsigned int call=0;
769 char fname[AIR_STRLEN_SMALL(128+1)];
770 sprintf(fname, "err-%04u.nrrd", call)__builtin___sprintf_chk (fname, 0, __builtin_object_size (fname
, 2 > 1 ? 1 : 0), "err-%04u.nrrd", call);
771 nrrdSave(fname, oscx->nerr, NULL((void*)0));
772 call++;
773 }
774 return 0;
775}
776
777/*
778** this can be called repeatedly
779*/
780static int
781_gageSetup(gageOptimSigContext *oscx) {
782 static const char me[]="_gageSetup";
783 double kparm[NRRD_KERNEL_PARMS_NUM8];
784 int E;
785
786 if (oscx->gctx) {
787 gageContextNix(oscx->gctx);
788 }
789 oscx->gctx = gageContextNew();
790 gageParmSet(oscx->gctx, gageParmVerbose, 0);
791 gageParmSet(oscx->gctx, gageParmRenormalize, AIR_FALSE0);
792 gageParmSet(oscx->gctx, gageParmCheckIntegrals, AIR_FALSE0);
793 gageParmSet(oscx->gctx, gageParmOrientationFromSpacing, AIR_TRUE1);
794 gageParmSet(oscx->gctx, gageParmStackUse, AIR_TRUE1);
795 E = 0;
796 if (!E) E |= !(oscx->pvlBase = gagePerVolumeNew(oscx->gctx,
797 oscx->nsampleImg[0],
798 gageKindScl));
799 if (!E) E |= gageStackPerVolumeNew(oscx->gctx, oscx->pvlSS,
800 AIR_CAST(const Nrrd*const*,((const Nrrd*const*)(oscx->nsampleImg))
801 oscx->nsampleImg)((const Nrrd*const*)(oscx->nsampleImg)),
802 oscx->sampleNum, gageKindScl);
803 if (!E) E |= gageStackPerVolumeAttach(oscx->gctx,
804 oscx->pvlBase, oscx->pvlSS,
805 oscx->sampleSigma, oscx->sampleNum);
806 kparm[0] = 1;
807 if (!E) E |= gageKernelSet(oscx->gctx, gageKernel00,
808 nrrdKernelBox, kparm);
809 if (!E) E |= gageKernelSet(oscx->gctx, gageKernelStack,
810 oscx->kssSpec->kernel, oscx->kssSpec->parm);
811 if (!E) E |= gageQueryItemOn(oscx->gctx, oscx->pvlBase, gageSclValue);
812 if (!E) E |= gageUpdate(oscx->gctx);
813 if (E) {
814 biffAddf(GAGEgageBiffKey, "%s: problem setting up gage", me);
815 return 1;
816 }
817 return 0;
818}
819
820static char *
821_timefmt(char tstr[AIR_STRLEN_MED(256+1)], double deltim) {
822
823 if (deltim < 60) {
824 sprintf(tstr, "%g secs", deltim)__builtin___sprintf_chk (tstr, 0, __builtin_object_size (tstr
, 2 > 1 ? 1 : 0), "%g secs", deltim);
825 return tstr;
826 }
827 deltim /= 60;
828 if (deltim < 60) {
829 sprintf(tstr, "%g mins", deltim)__builtin___sprintf_chk (tstr, 0, __builtin_object_size (tstr
, 2 > 1 ? 1 : 0), "%g mins", deltim);
830 return tstr;
831 }
832 deltim /= 60;
833 if (deltim < 24) {
834 sprintf(tstr, "%g hours", deltim)__builtin___sprintf_chk (tstr, 0, __builtin_object_size (tstr
, 2 > 1 ? 1 : 0), "%g hours", deltim);
835 return tstr;
836 }
837 deltim /= 24;
838 if (deltim < 7) {
839 sprintf(tstr, "%g days", deltim)__builtin___sprintf_chk (tstr, 0, __builtin_object_size (tstr
, 2 > 1 ? 1 : 0), "%g days", deltim);
840 return tstr;
841 }
842 deltim /= 7;
843 sprintf(tstr, "%g weeks", deltim)__builtin___sprintf_chk (tstr, 0, __builtin_object_size (tstr
, 2 > 1 ? 1 : 0), "%g weeks", deltim);
844 return tstr;
845}
846
847static int
848_optsigrun(gageOptimSigContext *oscx) {
849 static const char me[]="_optsigrun";
850 char tstr[AIR_STRLEN_MED(256+1)];
851 unsigned int iter, pnt;
852 double lastErr, newErr, rhoeps, oppor, lastPos, backoff, decavg, time0;
853 int badStep;
854
855 /* used to debug failure to measure error meaningfully
856 {
857 unsigned int hi, hn=200;
858 double rr;
859 for (hi=0; hi<hn; hi++) {
860 rr = AIR_AFFINE(-20, hi, hn+20, oscx->rhoRange[0], oscx->rhoRange[1]);
861 _sampleSet(oscx, 1, rr);
862 _errTotal(oscx);
863 }
864 }
865 */
866
867 time0 = airTime();
868 if (_errTotal(&lastErr, oscx)) {
869 biffAddf(GAGEgageBiffKey, "%s: first error measurement", me);
870 return 1;
871 }
872 fprintf(stderr__stderrp, "%s: (%s for initial error measr)\n", me,
873 _timefmt(tstr, airTime() - time0));
874 newErr = AIR_NAN(airFloatQNaN.f);
875 decavg = oscx->sampleNum; /* hack */
876 /* meaningful discrete difference for looking at error gradient is
877 bounded by the resolution of the sampling we're doing along scale */
878 rhoeps = 0.1*(oscx->rhoRange[1]-oscx->rhoRange[0])/oscx->trueImgNum;
879 oppor = 1.3333;
880 backoff = 0.25;
881 /* initialize step for the moving samples: 1 through oscx->sampleNum-2 */
882 for (pnt=1; pnt<oscx->sampleNum-1; pnt++) {
883 oscx->step[pnt] = 10;
884 }
885 for (iter=0; iter<oscx->maxIter; iter++) {
886 double limit, err1, grad, delta;
887 unsigned int tryi;
888 int zerodelta, esgn;
889 esgn = 2*AIR_CAST(int, airRandInt(2))((int)(airRandInt(2))) - 1;
890 pnt = 1 + (iter % (oscx->sampleNum-2));
891 lastPos = oscx->sampleRho[pnt];
892 fprintf(stderr__stderrp, "%s: ***** iter %u; [[ err %g ]] %s\n",
893 me, iter, lastErr, _timefmt(tstr, airTime() - time0));
894 limit = AIR_MIN((oscx->sampleRho[pnt] - oscx->sampleRho[pnt-1])/3,(((oscx->sampleRho[pnt] - oscx->sampleRho[pnt-1])/3) <
((oscx->sampleRho[pnt+1] - oscx->sampleRho[pnt])/3) ? (
(oscx->sampleRho[pnt] - oscx->sampleRho[pnt-1])/3) : ((
oscx->sampleRho[pnt+1] - oscx->sampleRho[pnt])/3))
895 (oscx->sampleRho[pnt+1] - oscx->sampleRho[pnt])/3)(((oscx->sampleRho[pnt] - oscx->sampleRho[pnt-1])/3) <
((oscx->sampleRho[pnt+1] - oscx->sampleRho[pnt])/3) ? (
(oscx->sampleRho[pnt] - oscx->sampleRho[pnt-1])/3) : ((
oscx->sampleRho[pnt+1] - oscx->sampleRho[pnt])/3));
896 fprintf(stderr__stderrp, ". pnt %u: pos %g, step %g\n",
897 pnt, lastPos, oscx->step[pnt]);
898 fprintf(stderr__stderrp, ". limit = min((%g-%g)/3,(%g-%g)/3) = %g\n",
899 oscx->sampleRho[pnt], oscx->sampleRho[pnt-1],
900 oscx->sampleRho[pnt+1], oscx->sampleRho[pnt], limit);
901 _sampleSet(oscx, pnt, lastPos + esgn*rhoeps);
902 if (_errTotal(&err1, oscx)) {
903 biffAddf(GAGEgageBiffKey, "%s: for err1 (%u -> %.17g)",
904 me, pnt, lastPos + esgn*rhoeps);
905 return 1;
906 }
907 _sampleSet(oscx, pnt, lastPos);
908 grad = (err1 - lastErr)/(esgn*rhoeps);
909 fprintf(stderr__stderrp, ". grad = %g\n", grad);
910 delta = -grad*oscx->step[pnt];
911 if (!AIR_EXISTS(delta)(((int)(!((delta) - (delta)))))) {
912 biffAddf(GAGEgageBiffKey, "%s: got non-exist delta %g on iter %u (pnt %u) err %g",
913 me, delta, iter, pnt, lastErr);
914 return 1;
915 }
916 if (AIR_ABS(delta)((delta) > 0.0f ? (delta) : -(delta)) > limit) {
917 oscx->step[pnt] *= limit/AIR_ABS(delta)((delta) > 0.0f ? (delta) : -(delta));
918 fprintf(stderr__stderrp, ". step *= %g/%g -> %g\n",
919 limit, AIR_ABS(delta)((delta) > 0.0f ? (delta) : -(delta)), oscx->step[pnt]);
920 delta = -grad*oscx->step[pnt];
921 }
922 fprintf(stderr__stderrp, ". delta = %g\n", delta);
923 tryi = 0;
924 badStep = AIR_FALSE0;
925 do {
926 if (tryi == oscx->maxIter) {
927 biffAddf(GAGEgageBiffKey, "%s: confusion (tryi %u) on iter %u (pnt %u) err %g",
928 me, tryi, iter, pnt, lastErr);
929 return 1;
930 }
931 if (!delta) {
932 fprintf(stderr__stderrp, "... try %u: delta = 0; nothing to do\n", tryi);
933 newErr = lastErr;
934 zerodelta = AIR_TRUE1;
935 } else {
936 zerodelta = AIR_FALSE0;
937 _sampleSet(oscx, pnt, lastPos + delta);
938 if (_errTotal(&newErr, oscx)) {
939 biffAddf(GAGEgageBiffKey, "%s: for newErr (%u -> %.17g)", me,
940 pnt, lastPos + delta);
941 return 1;
942 }
943 badStep = newErr > lastErr;
944 fprintf(stderr__stderrp, "... try %u: pos[%u] %g + %g = %g;\n"
945 "%s: err %g %s %g\n",
946 tryi, pnt, lastPos, delta,
947 oscx->sampleRho[pnt],
948 badStep ? "*BAD*" : "good",
949 newErr, newErr > lastErr ? ">" : "<=", lastErr);
950 if (badStep) {
951 oscx->step[pnt] *= backoff;
952 if (oscx->step[pnt] < rhoeps/1000) {
953 /* step got so small its stupid to be moving this point */
954 fprintf(stderr__stderrp, "... !! step %g < %g pointlessly small, "
955 "moving on\n", oscx->step[pnt], rhoeps/1000);
956 _sampleSet(oscx, pnt, lastPos);
957 newErr = lastErr;
958 badStep = AIR_FALSE0;
959 } else {
960 delta = -grad*oscx->step[pnt];
961 }
962 }
963 }
964 tryi++;
965 } while (badStep);
966 if (!zerodelta) {
967 /* don't update decavg if we moved on because slope was EXACTLY zero */
968 decavg = AIR_AFFINE(0, 1, oscx->sampleNum,( ((double)((lastErr - newErr)/lastErr)-(decavg))*((double)(1
)-(0)) / ((double)(oscx->sampleNum)-(0)) + (decavg))
969 decavg, (lastErr - newErr)/lastErr)( ((double)((lastErr - newErr)/lastErr)-(decavg))*((double)(1
)-(0)) / ((double)(oscx->sampleNum)-(0)) + (decavg));
970 oscx->step[pnt] *= oppor;
971 }
972 if (decavg <= oscx->convEps) {
973 fprintf(stderr__stderrp, "%s: converged (%g <= %g) after %u iters\n", me,
974 decavg, oscx->convEps, iter);
975 break;
976 } else {
977 fprintf(stderr__stderrp, "%s: _____ iter %u done; decavg = %g > eps %g\n", me,
978 iter, decavg, oscx->convEps);
979 }
980 lastErr = newErr;
981 }
982 if (iter == oscx->maxIter && decavg > oscx->convEps) {
983 biffAddf(GAGEgageBiffKey, "%s: failed to converge (%g > %g) after %u iters\n", me,
984 decavg, oscx->convEps, iter);
985 return 1;
986 }
987 oscx->finalErr = lastErr;
988
989 return 0;
990}
991
992static int
993_optsigrunLinf(gageOptimSigContext *oscx) {
994 static const char me[]="_optsigrunLinf";
995 char tstr[AIR_STRLEN_MED(256+1)];
996 double *srho, *stmp, time0, lastErr, newErr, decavg,
997 step, oppor, backoff, ceps, mmErr[2];
998 unsigned int iter, si, sn, mmIdx[2];
999 int shrink;
1000
1001 time0 = airTime();
1002 if (_errTotalLinf(&lastErr, oscx, mmIdx, mmErr)) {
1003 biffAddf(GAGEgageBiffKey, "%s: first error measurement", me);
1004 return 1;
1005 }
1006 fprintf(stderr__stderrp, "%s: (init) min %u %g max %u %g\n", me,
1007 mmIdx[0], mmErr[0], mmIdx[1], mmErr[1]);
1008 fprintf(stderr__stderrp, "%s: (%s for initial error measr)\n", me,
1009 _timefmt(tstr, airTime() - time0));
1010 newErr = AIR_NAN(airFloatQNaN.f);
1011
1012 /* shorcuts */
1013 sn = oscx->sampleNum;
1014 srho = oscx->sampleRho;
1015 stmp = oscx->sampleTmp;
1016
1017 /* Linf uses a single scalar step istead of oscx->step array */
1018 step = 0.1;
1019 oppor = 1.1;
1020 backoff = 0.5;
1021
1022 /* more demanding for more points */
1023 ceps = oscx->convEps/sn;
1024
1025 decavg = 2*ceps;
1026 for (iter=0; iter<oscx->maxIter; iter++) {
1027 double midp, wlo, whi, glo, ghi, gerr;
1028 unsigned int gap, tryi;
1029 int badStep;
1030
1031 if (iter % 2) {
1032 /* we will grow gap around smallest peak */
1033 gap = mmIdx[0];
1034 gerr = mmErr[0];
1035 shrink = AIR_FALSE0;
1036 } else {
1037 /* we will shrink gap around tallest peak */
1038 gap = mmIdx[1];
1039 gerr = mmErr[1];
1040 shrink = AIR_TRUE1;
1041 }
1042 midp = (srho[gap] + srho[gap+1])/2;
1043 fprintf(stderr__stderrp, "%s: ---- iter %u (step %g): gap [%u]/%g (%s)\n",
1044 me, iter, step, gap, gerr,
1045 shrink ? "shrinking tallest" : "growing lowest");
1046 /* save last set of positions to restore after bad step */
1047 for (si=1; si<sn-1; si++) {
1048 stmp[si] = srho[si];
1049 }
1050 tryi = 0;
1051 badStep = AIR_FALSE0;
1052 do {
1053 if (tryi == oscx->maxIter) {
1054 biffAddf(GAGEgageBiffKey, "%s: confusion (tryi %u) on iter %u err %g",
1055 me, tryi, iter, lastErr);
1056 return 1;
1057 }
1058 if (shrink) {
1059 wlo = AIR_AFFINE(0, step, 1, srho[gap], midp)( ((double)(midp)-(srho[gap]))*((double)(step)-(0)) / ((double
)(1)-(0)) + (srho[gap]));
1060 whi = AIR_AFFINE(0, step, 1, srho[gap+1], midp)( ((double)(midp)-(srho[gap+1]))*((double)(step)-(0)) / ((double
)(1)-(0)) + (srho[gap+1]));
1061 } else {
1062 wlo = AIR_AFFINE(0, step, -2, srho[gap], midp)( ((double)(midp)-(srho[gap]))*((double)(step)-(0)) / ((double
)(-2)-(0)) + (srho[gap]));
1063 whi = AIR_AFFINE(0, step, -2, srho[gap+1], midp)( ((double)(midp)-(srho[gap+1]))*((double)(step)-(0)) / ((double
)(-2)-(0)) + (srho[gap+1]));
1064 }
1065 glo = srho[gap];
1066 ghi = srho[gap+1];
1067 fprintf(stderr__stderrp, "%s: rho[%u] %g | %g -- rho[%u] %g | %g\n", me,
1068 gap, srho[gap], wlo, gap+1, srho[gap+1], whi);
1069 for (si=1; si<sn-1; si++) {
1070 _sampleSet(oscx, si, (si <= gap
1071 ? AIR_AFFINE(srho[0], srho[si], glo,( ((double)(wlo)-(srho[0]))*((double)(srho[si])-(srho[0])) / (
(double)(glo)-(srho[0])) + (srho[0]))
1072 srho[0], wlo)( ((double)(wlo)-(srho[0]))*((double)(srho[si])-(srho[0])) / (
(double)(glo)-(srho[0])) + (srho[0]))
1073 : AIR_AFFINE(ghi, srho[si], srho[sn-1],( ((double)(srho[sn-1])-(whi))*((double)(srho[si])-(ghi)) / (
(double)(srho[sn-1])-(ghi)) + (whi))
1074 whi, srho[sn-1])( ((double)(srho[sn-1])-(whi))*((double)(srho[si])-(ghi)) / (
(double)(srho[sn-1])-(ghi)) + (whi))));
1075 }
1076 if (_errTotalLinf(&newErr,
1077 oscx, mmIdx, mmErr)) {
1078 biffAddf(GAGEgageBiffKey, "%s: iter %u", me, iter);
1079 return 1;
1080 }
1081 fprintf(stderr__stderrp, "%s: min %u %g max %u %g\n", me,
1082 mmIdx[0], mmErr[0], mmIdx[1], mmErr[1]);
1083 if (iter % 3) {
1084 badStep = newErr > lastErr;
1085 fprintf(stderr__stderrp, "... try %u [%u] step %g -> newErr %g %s "
1086 "lastErr %g %s\n",
1087 tryi, gap, step, newErr,
1088 badStep ? ">" : "<=",
1089 lastErr,
1090 badStep ? "*BAD*" : "good");
1091 if (badStep) {
1092 step *= backoff;
1093 for (si=1; si<sn-1; si++) {
1094 srho[si] = stmp[si];
1095 }
1096 }
1097 tryi++;
1098 }
1099 } while (badStep);
1100 step *= oppor;
1101 decavg = (decavg + (lastErr - newErr))/2;
1102 if (AIR_IN_OP(0, decavg, ceps)((0) < (decavg) && (decavg) < (ceps))) {
1103 fprintf(stderr__stderrp, "%s: converged (%g <= %g) after %u iters\n", me,
1104 decavg, ceps, iter);
1105 break;
1106 } else {
1107 fprintf(stderr__stderrp, "%s: iter %u done; decavg = %g > eps %g\n", me,
1108 iter, decavg, ceps);
1109 }
1110 lastErr = newErr;
1111 } /* for iter */
1112 if (oscx->maxIter
1113 && iter == oscx->maxIter
1114 && decavg > ceps) {
1115 biffAddf(GAGEgageBiffKey, "%s: failed to converge (%g > %g) after %u iters\n", me,
1116 decavg, ceps, iter);
1117 return 1;
1118 }
1119 oscx->finalErr = lastErr;
1120
1121 return 0;
1122}
1123
1124int
1125gageOptimSigCalculate(gageOptimSigContext *oscx,
1126 /* output */ double *sigma,
1127 unsigned int sigmaNum,
1128 const NrrdKernelSpec *kssSpec,
1129 int imgMeasr, int allMeasr,
1130 unsigned int maxIter, double convEps) {
1131 static const char me[]="gageOptimSigCalculate";
1132 unsigned int ii;
1133
1134 if (!( oscx && sigma && kssSpec )) {
1135 biffAddf(GAGEgageBiffKey, "%s: got NULL pointer", me);
1136 return 1;
1137 }
1138 if (sigmaNum > oscx->sampleNumMax) {
1139 biffAddf(GAGEgageBiffKey, "%s: initialized for max %u samples, not %u", me,
1140 oscx->sampleNumMax, sigmaNum);
1141 return 1;
1142 }
1143
1144 /* initialize to uniform samples in rho */
1145 oscx->sampleNum = sigmaNum;
1146 fprintf(stderr__stderrp, "%s: initializing %u samples ... ", me, oscx->sampleNum);
1147 fflush(stderr__stderrp);
1148 for (ii=0; ii<oscx->sampleNum; ii++) {
1149 _sampleSet(oscx, ii, AIR_AFFINE(0, ii, oscx->sampleNum-1,( ((double)(oscx->rhoRange[1])-(oscx->rhoRange[0]))*((double
)(ii)-(0)) / ((double)(oscx->sampleNum-1)-(0)) + (oscx->
rhoRange[0]))
1150 oscx->rhoRange[0], oscx->rhoRange[1])( ((double)(oscx->rhoRange[1])-(oscx->rhoRange[0]))*((double
)(ii)-(0)) / ((double)(oscx->sampleNum-1)-(0)) + (oscx->
rhoRange[0])));
1151 }
1152 fprintf(stderr__stderrp, "done.\n");
1153
1154 /* copy remaining input parameters */
1155 nrrdKernelSpecNix(oscx->kssSpec);
1156 oscx->kssSpec = nrrdKernelSpecCopy(kssSpec);
1157 oscx->imgMeasr = imgMeasr;
1158 oscx->allMeasr = allMeasr;
1159 oscx->convEps = convEps;
1160 oscx->maxIter = maxIter;
1161 oscx->convEps = convEps;
1162
1163 /* set up gage */
1164 fprintf(stderr__stderrp, "%s: setting up gage ... \n", me);
1165 if (_gageSetup(oscx)) {
1166 biffAddf(GAGEgageBiffKey, "%s: problem setting up gage", me);
1167 return 1;
1168 }
1169 fprintf(stderr__stderrp, "%s: ... gage setup done.\n", me);
1170
1171 /* run the optimization */
1172 if (oscx->sampleNum > 2) {
1173 int EE;
1174 EE = (nrrdMeasureLinf == oscx->allMeasr
1175 ? _optsigrunLinf(oscx)
1176 : _optsigrun(oscx));
1177 if (EE) {
1178 biffAddf(GAGEgageBiffKey, "%s: trouble", me);
1179 return 1;
1180 }
1181 } else {
1182 fprintf(stderr__stderrp, "%s: num == 2, no optimization, finding error ... ", me);
1183 fflush(stderr__stderrp);
1184 if (_errTotal(&(oscx->finalErr), oscx)) {
1185 biffAddf(GAGEgageBiffKey, "%s: for finalErr", me);
1186 return 1;
1187 }
1188 fprintf(stderr__stderrp, "done.\n");
1189 }
1190
1191 /* save output */
1192 for (ii=0; ii<oscx->sampleNum; ii++) {
1193 sigma[ii] = oscx->sampleSigma[ii];
1194 }
1195
1196 return 0;
1197}
1198
1199int
1200gageOptimSigErrorPlot(gageOptimSigContext *oscx, Nrrd *nout,
1201 const double *sigma,
1202 unsigned int sigmaNum,
1203 const NrrdKernelSpec *kssSpec,
1204 int imgMeasr) {
1205 static const char me[]="gageOptimSigErrorPlot";
1206 char doneStr[AIR_STRLEN_SMALL(128+1)];
1207 double *out;
1208 unsigned int ii;
1209
1210 if (!(oscx && nout && sigma)) {
1211 biffAddf(GAGEgageBiffKey, "%s: got NULL pointer", me);
1212 return 1;
1213 }
1214 if (!( sigmaNum >= 2 )) {
1215 biffAddf(GAGEgageBiffKey, "%s: need sigmaNum >= 2 (not %u)", me, sigmaNum);
1216 return 1;
1217 }
1218 if (sigmaNum > oscx->sampleNumMax) {
1219 biffAddf(GAGEgageBiffKey, "%s: initialized for max %u samples, not %u", me,
1220 oscx->sampleNumMax, sigmaNum);
1221 return 1;
1222 }
1223
1224 /* copy remaining input parms */
1225 nrrdKernelSpecNix(oscx->kssSpec);
1226 oscx->kssSpec = nrrdKernelSpecCopy(kssSpec);
1227 oscx->sampleNum = sigmaNum;
1228 oscx->maxIter = 0;
1229 oscx->imgMeasr = imgMeasr;
1230 oscx->allMeasr = nrrdMeasureUnknown;
1231 oscx->convEps = AIR_NAN(airFloatQNaN.f);
1232 if (nrrdMaybeAlloc_va(nout, nrrdTypeDouble, 2,
1233 AIR_CAST(size_t, 2)((size_t)(2)),
1234 AIR_CAST(size_t, oscx->trueImgNum)((size_t)(oscx->trueImgNum)))) {
1235 biffMovef(GAGEgageBiffKey, NRRDnrrdBiffKey, "%s: trouble allocating output", me);
1236 return 1;
1237 }
1238 out = AIR_CAST(double *, nout->data)((double *)(nout->data));
1239
1240 /* set up requested samples */
1241 for (ii=0; ii<oscx->sampleNum; ii++) {
1242 _sampleSet(oscx, ii, _RhoOfSig(sigma[ii]));
1243 }
1244 if (_gageSetup(oscx)) {
1245 biffAddf(GAGEgageBiffKey, "%s: problem setting up gage", me);
1246 return 1;
1247 }
1248 fprintf(stderr__stderrp, "%s: plotting ... ", me); fflush(stderr__stderrp);
1249 for (ii=0; ii<oscx->trueImgNum; ii++) {
1250 double rho, err;
1251 fprintf(stderr__stderrp, "%s", airDoneStr(0, ii, oscx->trueImgNum, doneStr));
1252 fflush(stderr__stderrp);
1253 rho = AIR_AFFINE(0, ii, oscx->trueImgNum-1,( ((double)(oscx->rhoRange[1])-(oscx->rhoRange[0]))*((double
)(ii)-(0)) / ((double)(oscx->trueImgNum-1)-(0)) + (oscx->
rhoRange[0]))
1254 oscx->rhoRange[0], oscx->rhoRange[1])( ((double)(oscx->rhoRange[1])-(oscx->rhoRange[0]))*((double
)(ii)-(0)) / ((double)(oscx->trueImgNum-1)-(0)) + (oscx->
rhoRange[0]));
1255 out[0 + 2*ii] = rho;
1256 /* debugii = ii; */
1257 if (_errSingle(&err, oscx, rho)) {
1258 biffAddf(GAGEgageBiffKey, "%s: plotting %u", me, ii);
1259 return 1;
1260 }
1261 out[1 + 2*ii] = err;
1262 /*
1263 if (AIR_IN_CL(69,ii,71)) {
1264 fprintf(stderr, "!%s: ----- %u : %g\n", me, ii, err);
1265 }
1266 */
1267 }
1268 fprintf(stderr__stderrp, "%s\n", airDoneStr(0, ii, oscx->trueImgNum, doneStr));
1269
1270 /*
1271 if (0) {
1272 static unsigned int call=0;
1273 char fname[AIR_STRLEN_SMALL];
1274 unsigned int ni;
1275 if (0) {
1276 sprintf(fname, "err-%04u.nrrd", call);
1277 nrrdSave(fname, oscx->nerr, NULL);
1278 }
1279 if (debugging) {
1280 for (ni=0; ni<oscx->sampleNum; ni++) {
1281 sprintf(fname, "sample-%04u.nrrd", ni);
1282 nrrdSave(fname, oscx->nsampleImg[ni], NULL);
1283 }
1284 }
1285 call++;
1286 debugging = (2 == call);
1287 }
1288 */
1289
1290 return 0;
1291}
1292
1293int
1294gageOptimSigErrorPlotSliding(gageOptimSigContext *oscx, Nrrd *nout,
1295 double windowRho,
1296 unsigned int sampleNum,
1297 const NrrdKernelSpec *kssSpec,
1298 int imgMeasr) {
1299 static const char me[]="gageOptimSigRecondErrorPlotSliding";
1300 char doneStr[AIR_STRLEN_SMALL(128+1)];
1301 unsigned int ii;
1302 double *out;
1303 char hackKeyStr[]="TEEM_OPTSIG_RECONERR";
1304
1305 if (!( oscx && nout && kssSpec )) {
1306 biffAddf(GAGEgageBiffKey, "%s: got NULL pointer", me);
1307 return 1;
1308 }
1309 if (windowRho <= 0) {
1310 biffAddf(GAGEgageBiffKey, "%s: need positive windowRho (not %g)", me, windowRho);
1311 return 1;
1312 }
1313 if (windowRho > oscx->rhoRange[1] - oscx->rhoRange[0]) {
1314 biffAddf(GAGEgageBiffKey, "%s: window %g > rhorange %g-%g=%g", me,
1315 windowRho, oscx->rhoRange[1], oscx->rhoRange[0],
1316 oscx->rhoRange[1] - oscx->rhoRange[0]);
1317 return 1;
1318 }
1319
1320 if (nrrdGetenvString(&debugReconErrName, hackKeyStr)) {
1321 fprintf(stderr__stderrp, "%s: %s hack on: will save recon results to %s\n",
1322 me, hackKeyStr, debugReconErrName);
1323 debugReconErrArr = airArrayNew(AIR_CAST(void **, &(debugReconErr))((void **)(&(debugReconErr))), NULL((void*)0),
1324 sizeof(double), 1000);
1325 }
1326
1327 /* copy remaining input parms */
1328 nrrdKernelSpecNix(oscx->kssSpec);
1329 oscx->kssSpec = nrrdKernelSpecCopy(kssSpec);
1330 oscx->sampleNum = 3; /* hacky */
1331 oscx->maxIter = 0;
1332 oscx->imgMeasr = imgMeasr;
1333 oscx->allMeasr = nrrdMeasureUnknown;
1334 oscx->convEps = AIR_NAN(airFloatQNaN.f);
1335 oscx->sampleSigma[0] = oscx->sigmaRange[0]; /* just for gage setup */
1336 oscx->sampleSigma[1] = oscx->sigmaRange[1]; /* just for gage setup */
1337 oscx->sampleSigma[2] = oscx->sigmaRange[1]+1;
1338 if (_gageSetup(oscx)) {
1339 biffAddf(GAGEgageBiffKey, "%s: problem setting up gage", me);
1340 return 1;
1341 }
1342 if (nrrdMaybeAlloc_va(nout, nrrdTypeDouble, 2,
1343 AIR_CAST(size_t, 2)((size_t)(2)),
1344 AIR_CAST(size_t, sampleNum)((size_t)(sampleNum)))) {
1345 biffMovef(GAGEgageBiffKey, NRRDnrrdBiffKey, "%s: trouble allocating output", me);
1346 return 1;
1347 }
1348 out = AIR_CAST(double *, nout->data)((double *)(nout->data));
1349 fprintf(stderr__stderrp, "%s: plotting ... ", me); fflush(stderr__stderrp);
1350 for (ii=0; ii<sampleNum; ii++) {
1351 double rho, rlo, rhi, err;
1352 fprintf(stderr__stderrp, "%s", airDoneStr(0, ii, oscx->trueImgNum, doneStr));
1353 fflush(stderr__stderrp);
1354 rho = AIR_AFFINE(0, ii, sampleNum,( ((double)(oscx->rhoRange[1])-(oscx->rhoRange[0]))*((double
)(ii)-(0)) / ((double)(sampleNum)-(0)) + (oscx->rhoRange[0
]))
1355 oscx->rhoRange[0], oscx->rhoRange[1])( ((double)(oscx->rhoRange[1])-(oscx->rhoRange[0]))*((double
)(ii)-(0)) / ((double)(sampleNum)-(0)) + (oscx->rhoRange[0
]));
1356 rlo = rho - windowRho/2;
1357 rhi = rho + windowRho/2;
1358 if (rlo < oscx->rhoRange[0]
1359 || rhi > oscx->rhoRange[1]) {
1360 if (debugReconErrArr) {
1361 /* we have to simulate the updates to debugReconErrArr
1362 that would happen with a call to _errSingle */
1363 airArrayLenIncr(debugReconErrArr, 2*oscx->sz*oscx->sy*oscx->sx);
1364 }
1365 out[0 + 2*ii] = _SigOfRho(rho);
1366 out[1 + 2*ii] = AIR_NAN(airFloatQNaN.f);
1367 continue;
1368 }
1369 /* required samples will slide along with plotting */
1370 _sampleSet(oscx, 0, rlo);
1371 _sampleSet(oscx, 1, rhi);
1372 if (_errSingle(&err, oscx, rho)) {
1373 biffAddf(GAGEgageBiffKey, "%s: plotting/sliding %u", me, ii);
1374 return 1;
1375 }
1376 out[0 + 2*ii] = _SigOfRho(rho);
1377 out[1 + 2*ii] = err;
1378 }
1379 fprintf(stderr__stderrp, "%s\n", airDoneStr(0, ii, oscx->trueImgNum, doneStr));
1380
1381 if (debugReconErrArr) {
1382 Nrrd *nre = nrrdNew();
1383 nrrdWrap_va(nre, debugReconErr, nrrdTypeDouble, 3,
1384 AIR_CAST(size_t, 2)((size_t)(2)),
1385 AIR_CAST(size_t, oscx->sz*oscx->sy*oscx->sx)((size_t)(oscx->sz*oscx->sy*oscx->sx)),
1386 AIR_CAST(size_t, sampleNum)((size_t)(sampleNum)));
1387 nrrdSave(debugReconErrName, nre, NULL((void*)0));
1388 nrrdNix(nre);
1389 }
1390
1391
1392 return 0;
1393}
1394