CR2RE Pipeline Reference Manual 1.7.0
hdrl_fit-test.c
1/*
2 * This file is part of the HDRL
3 * Copyright (C) 2014 European Southern Observatory
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20#ifdef HAVE_CONFIG_H
21#include <config.h>
22#endif
23
24/*-----------------------------------------------------------------------------
25 Includes
26 -----------------------------------------------------------------------------*/
27
28#include "hdrl_fit.c"
29#include "hdrl_imagelist.h"
30#include "hdrl_test.h"
31
32#include <cpl.h>
33
34#include <math.h>
35
36#ifndef ARRAY_LEN
37#define ARRAY_LEN(a) sizeof((a))/sizeof((a)[0])
38#endif
39
40/*----------------------------------------------------------------------------*/
44/*----------------------------------------------------------------------------*/
45
46#define matrix_eq(m, exp, eps) \
47 do { \
48 for (size_t i = 0; i < (size_t)cpl_matrix_get_nrow(m); i++) { \
49 for (size_t j = 0; j < (size_t)cpl_matrix_get_ncol(m); j++) { \
50 cpl_test_abs(cpl_matrix_get(m, i, j), exp[i * \
51 cpl_matrix_get_ncol(m) + j], eps); \
52 } \
53 } \
54 } while (0)
55
56static void test_vander1d(void)
57{
58 {
59 double p[] = {1, 2, 3, 4};
60 double exp[] = {1., 1., 1.,
61 1., 2., 4.,
62 1., 3., 9.,
63 1., 4., 16.};
64 cpl_vector * pv = cpl_vector_wrap(ARRAY_LEN(p), p);
65 cpl_matrix * v = polyvander1d(pv, 2);
66 cpl_test_error(CPL_ERROR_NONE);
67 matrix_eq(v, exp, DBL_EPSILON * 5);
68 cpl_matrix_delete(v);
69
70 v = polyvander1d(pv, 2);
71 cpl_test_error(CPL_ERROR_NONE);
72 matrix_eq(v, exp, DBL_EPSILON * 5);
73 cpl_matrix_delete(v);
74
75 cpl_vector_unwrap(pv);
76 }
77 /* could be used for weight_design function */
78#if 0
79 {
80 double p[] = {1, 2, 3, 4};
81 double e[] = {1, 1, 4, 1};
82 double exp[] = {1., 1., 1.,
83 1., 2., 4.,
84 0.25, 0.75, 2.25,
85 1., 4., 16.};
86 cpl_vector * pv = cpl_vector_wrap(ARRAY_LEN(p), p);
87 cpl_vector * ev = cpl_vector_wrap(ARRAY_LEN(e), e);
88 cpl_matrix * v = polyvander1d(pv, 2);
89 cpl_test_error(CPL_ERROR_NONE);
90
91 matrix_eq(v, exp, DBL_EPSILON * 5);
92 cpl_vector_unwrap(pv);
93 cpl_vector_unwrap(ev);
94 cpl_matrix_delete(v);
95 }
96#endif
97}
98
99static void test_fit(void)
100{
101 {
102 double s[] = {1, 2, 3, 4};
103 double p[] = {2, 2.5, 3, 3.5};
104 double exp[] = {1.5, 0.5, 0.};
105 double eres[] = {0., 0., 0., 0.};
106 cpl_vector * sv = cpl_vector_wrap(ARRAY_LEN(s), s);
107 cpl_vector * pv = cpl_vector_wrap(ARRAY_LEN(p), p);
108 cpl_vector * vre = cpl_vector_wrap(ARRAY_LEN(eres), eres);
109 cpl_matrix * v = polyvander1d(sv, 2);
110 hdrl_ls_fit_result * r = fit(v, pv, NULL);
111 cpl_test_error(CPL_ERROR_NONE);
112 matrix_eq(r->coef, exp, DBL_EPSILON * 10);
113 cpl_vector * res = hdrl_ls_fit_result_get_residuals(r, pv);
114 cpl_test_vector_abs(res, vre, DBL_EPSILON * 5);
115 cpl_vector_unwrap(vre);
116 cpl_vector_delete(res);
117 cpl_matrix_delete(v);
118 cpl_vector_unwrap(pv);
119 cpl_vector_unwrap(sv);
120 hdrl_ls_fit_result_delete(r);
121 }
122 {
123 double s[] = {1, 2, 3, 4, 5};
124 double p[] = {1.1, 2.5, 3.4, 3.8, 7};
125 double e[] = {0.3, 0.2, 0.2, 0.1, 0.5};
126 double exp[] = {0.54529, 0.858981};
127 double cexp[] = {0.0756216, -0.0206541, -0.0206541, 0.00613226};
128 cpl_vector * vfit = cpl_vector_new(ARRAY_LEN(s));
129 cpl_vector * vres = cpl_vector_new(ARRAY_LEN(s));
130 for (size_t i = 0; i < ARRAY_LEN(s); i++) {
131 cpl_vector_set(vfit, i, exp[0] + exp[1] * s[i]);
132 cpl_vector_set(vres, i, p[i] - (exp[0] + exp[1] * s[i]));
133 }
134
135 cpl_vector * sv = cpl_vector_wrap(ARRAY_LEN(s), s);
136 cpl_vector * pv = cpl_vector_wrap(ARRAY_LEN(p), p);
137 cpl_vector * ev = cpl_vector_wrap(ARRAY_LEN(e), e);
138 cpl_matrix * v = polyvander1d(sv, 1);
139
140 hdrl_ls_fit_result * r = fit(v, pv, ev);
141 cpl_test_error(CPL_ERROR_NONE);
142 matrix_eq(r->coef, exp, DBL_EPSILON * 1e10);
143 matrix_eq(r->cov, cexp, DBL_EPSILON * 1e10);
144
145 cpl_vector * values = hdrl_ls_fit_result_get_fitted_values(r);
146 cpl_test_vector_abs(values, vfit, DBL_EPSILON * 1e10);
147 cpl_vector_delete(vfit);
148 cpl_vector_delete(values);
149 cpl_vector * resi = hdrl_ls_fit_result_get_residuals(r, pv);
150 cpl_test_vector_abs(resi, vres, DBL_EPSILON * 1e10);
151 cpl_vector_delete(vres);
152 cpl_vector_delete(resi);
153 cpl_matrix_delete(v);
154
155 hdrl_ls_fit_result_delete(r);
156
157 r = polyfit1d(sv, pv, ev, 1);
158 cpl_test_error(CPL_ERROR_NONE);
159 matrix_eq(r->coef, exp, DBL_EPSILON * 1e10);
160 hdrl_ls_fit_result_delete(r);
161
162 cpl_vector_unwrap(pv);
163 cpl_vector_unwrap(ev);
164 cpl_vector_unwrap(sv);
165 }
166}
167
168/* ---------------------------------------------------------------------------*/
182/* ---------------------------------------------------------------------------*/
183static inline cpl_vector * hdrl_ls_fit_result_get_fit_interval(
184 const hdrl_ls_fit_result * r,
185 const cpl_vector * data,
186 cpl_vector * errors)
187{
188 /* mse = sum(sqrt(weights) * residuals ** 2) / df */
189 double mse = hdrl_ls_fit_result_get_chi2(r, data, errors) /
190 hdrl_ls_fit_result_get_residual_dof(r);
191 /* var = mse / weights */
192 cpl_vector * serror = cpl_vector_duplicate(errors);
193 cpl_vector_multiply(serror, serror);
194 cpl_vector_multiply_scalar(serror, mse);
195 cpl_vector_power(serror, 0.5);
196
197 /* TODO: accounting for covariance missing
198 * + (exog * np.dot(covb, exog.T).T).sum(axis=1) */
199
200 return serror;
201}
202
203
204void test_poisson(void)
205{
206 double x[] = { 10. , 62.1 , 114.2, 166.3, 218.4, 270.5, 322.6,
207 374.7, 426.8, 478.9, 531.1, 583.2, 635.3, 687.4, 739.5,
208 791.6, 843.7, 895.8, 947.9, 1000. };
209 /* poisson data with lambda: x / 10 (== variance) and 100 offset */
210 double y[] = { 103, 107, 111, 112, 117, 127, 126, 125, 139, 150, 157, 162,
211 153, 158, 162, 184, 191, 195, 182, 196 };
212 cpl_vector * vx = cpl_vector_wrap(ARRAY_LEN(x), x);
213 cpl_vector * vy = cpl_vector_wrap(ARRAY_LEN(y), y);
214 /* relative model errors (poisson model ~ sqrt(x)) */
215 cpl_vector * ve_model = cpl_vector_duplicate(vx);
216 cpl_vector_power(ve_model, 0.5);
217 /* real absolute error of the population != relative model errors */
218 cpl_vector * ve_real = cpl_vector_duplicate(vx);
219 cpl_vector_divide_scalar(ve_real, 10);
220 cpl_vector_power(ve_real, 0.5);
221
222
223 double exp_c[] = {101.4164, 0.0919476};
224 cpl_matrix * design = polyvander1d(vx, 1);
225 hdrl_ls_fit_result * res = polyfit1d(vx, vy, ve_model, 1);
226 matrix_eq(res->coef, exp_c, DBL_EPSILON * 2e12);
227 /* get sample error */
228 cpl_vector * pred_e =
229 hdrl_ls_fit_result_get_fit_interval(res, vy, ve_model);
230 /* should be less than < 20% deviation to population error */
231 cpl_test_vector_abs(pred_e, ve_real, 0.7);
232 cpl_vector_delete(pred_e);
233
234 cpl_vector_unwrap(vx);
235 cpl_vector_unwrap(vy);
236 cpl_vector_delete(ve_model);
237 cpl_vector_delete(ve_real);
238 cpl_matrix_delete(design);
239 hdrl_ls_fit_result_delete(res);
240}
241
242void test_imglistfit(void)
243{
244 {
245 size_t n = 5;
246 hdrl_imagelist * input = hdrl_imagelist_new();
247 cpl_vector * samp = cpl_vector_new(n);
248 hdrl_imagelist * out_coef = NULL;
249 cpl_image * out_chi2 = NULL;
250 cpl_image * out_dof = NULL;
251 hdrl_fit_polynomial_imagelist(NULL, NULL, 0, NULL, NULL, NULL);
252 cpl_test_error(CPL_ERROR_NULL_INPUT);
253
254 hdrl_fit_polynomial_imagelist(input, NULL, 0, NULL, NULL, NULL);
255 cpl_test_error(CPL_ERROR_NULL_INPUT);
256
257 hdrl_fit_polynomial_imagelist(NULL, samp, 0, NULL, NULL, NULL);
258 cpl_test_error(CPL_ERROR_NULL_INPUT);
259
260 hdrl_fit_polynomial_imagelist(input, samp, 0, NULL, NULL, NULL);
261 cpl_test_error(CPL_ERROR_NULL_INPUT);
262
263 hdrl_fit_polynomial_imagelist(input, samp, -1, &out_coef, &out_chi2, &out_dof);
264 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
265
266 hdrl_fit_polynomial_imagelist(input, samp, n + 2, &out_coef, &out_chi2, &out_dof);
267 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
268
269 for (size_t i = 0; i < n; i++) {
270 double t = (i + 1) * 100;
271 hdrl_image * img = hdrl_image_new(10, 10);
272 /* slight deviation from linear to get nonzero scale */
274 (hdrl_value){0.5 * t + 100 - i, sqrt(0.5 * t)});
275 if (i == 3) {
276 hdrl_image_reject(img, 3, 4);
277 }
278 hdrl_imagelist_set(input, img, i);
279 cpl_vector_set(samp, i, t);
280 }
281
282 hdrl_fit_polynomial_imagelist(input, samp, 1, &out_coef, &out_chi2, &out_dof);
283 cpl_test_error(CPL_ERROR_NONE);
284 cpl_test_eq(cpl_image_get_type(out_dof), HDRL_TYPE_DATA);
285 cpl_test_eq(cpl_image_get_type(out_chi2), HDRL_TYPE_DATA);
286 cpl_test_eq(hdrl_imagelist_get_size(out_coef), 2);
287
288 hdrl_image * coef0 = hdrl_imagelist_get(out_coef, 0);
289 hdrl_image * expect = hdrl_image_new(10, 10);
290 hdrl_image_add_scalar(expect, (hdrl_value){101, 9.0045});
291 hdrl_image_set_pixel(expect, 3, 4, (hdrl_value){101, 9.29448});
292 hdrl_test_image_abs(coef0, expect, HDRL_EPS_DATA * 1e11);
293 hdrl_image_delete(expect);
294
295 hdrl_image * coef1 = hdrl_imagelist_get(out_coef, 1);
296 expect = hdrl_image_new(10, 10);
297 hdrl_image_add_scalar(expect, (hdrl_value){0.49, 0.0351317});
298 hdrl_image_set_pixel(expect, 3, 4, (hdrl_value){0.49, 0.0399607});
299 hdrl_test_image_abs(coef1, expect, HDRL_EPS_DATA * 1e11);
300 hdrl_image_delete(expect);
301
302 cpl_image * cexpect = cpl_image_new(10, 10, HDRL_TYPE_DATA);
303 /* dof*scaling between weights and measured error */
304 cpl_image_add_scalar(cexpect, 1.831e-29);
305 cpl_test_image_abs(out_chi2, cexpect, DBL_EPSILON * 1e9);
306 cpl_image_delete(cexpect);
307
308 cexpect = cpl_image_new(10, 10, HDRL_TYPE_DATA);
309 cpl_image_add_scalar(cexpect, 3);
310 cpl_image_set(cexpect, 3, 4, 2);
311 cpl_test_image_abs(out_dof, cexpect, 0);
312 cpl_image_delete(cexpect);
313
314 hdrl_imagelist_delete(out_coef);
315 cpl_image_delete(out_chi2);
316 cpl_image_delete(out_dof);
317
318
319 /* smoke test less good pixels than fit degree */
320 hdrl_image_reject(hdrl_imagelist_get(input, 0), 2, 2);
321 hdrl_image_reject(hdrl_imagelist_get(input, 1), 2, 2);
322 hdrl_image_reject(hdrl_imagelist_get(input, 2), 2, 2);
323 hdrl_image_reject(hdrl_imagelist_get(input, 3), 2, 2);
324 hdrl_fit_polynomial_imagelist(input, samp, 3, &out_coef, &out_chi2, &out_dof);
325
327 cpl_vector_delete(samp);
328 hdrl_imagelist_delete(out_coef);
329 cpl_image_delete(out_chi2);
330 cpl_image_delete(out_dof);
331 }
332}
333
334void test_imglistfit2(void)
335{
336 {
337 size_t n = 5;
338 hdrl_imagelist * input = hdrl_imagelist_new();
339 cpl_imagelist * samp = cpl_imagelist_new();
340 hdrl_imagelist * out_coef = NULL;
341 cpl_image * out_chi2 = NULL;
342 cpl_image * out_dof = NULL;
343 hdrl_fit_polynomial_imagelist2(NULL, NULL, 0, NULL, NULL, NULL);
344 cpl_test_error(CPL_ERROR_NULL_INPUT);
345
346 hdrl_fit_polynomial_imagelist2(input, NULL, 0, NULL, NULL, NULL);
347 cpl_test_error(CPL_ERROR_NULL_INPUT);
348
349 hdrl_fit_polynomial_imagelist2(NULL, samp, 0, NULL, NULL, NULL);
350 cpl_test_error(CPL_ERROR_NULL_INPUT);
351
352 hdrl_fit_polynomial_imagelist2(input, samp, 0, NULL, NULL, NULL);
353 cpl_test_error(CPL_ERROR_NULL_INPUT);
354
355 hdrl_fit_polynomial_imagelist2(input, samp, -1, &out_coef, &out_chi2, &out_dof);
356 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
357
358 hdrl_fit_polynomial_imagelist2(input, samp, n + 2, &out_coef, &out_chi2, &out_dof);
359 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
360
361 for (size_t i = 0; i < n; i++) {
362 double t = (i + 1) * 100;
363 hdrl_image * img = hdrl_image_new(10, 10);
364 /* slight deviation from linear to get nonzero scale */
366 (hdrl_value){0.5 * t + 100 - i, sqrt(0.5 * t)});
367 hdrl_imagelist_set(input, img, i);
368 cpl_image * sampi = cpl_image_new(10, 10, HDRL_TYPE_DATA);
369 cpl_image_add_scalar(sampi, t);
370 if (i == 3) {
371 cpl_image_reject(sampi, 3, 4);
372 }
373 cpl_imagelist_set(samp, sampi, i);
374 }
375
376 hdrl_fit_polynomial_imagelist2(input, samp, 1, &out_coef, &out_chi2, &out_dof);
377 cpl_test_error(CPL_ERROR_NONE);
378 cpl_test_eq(cpl_image_get_type(out_dof), HDRL_TYPE_DATA);
379 cpl_test_eq(cpl_image_get_type(out_chi2), HDRL_TYPE_DATA);
380 cpl_test_eq(hdrl_imagelist_get_size(out_coef), 2);
381
382 hdrl_image * coef0 = hdrl_imagelist_get(out_coef, 0);
383 hdrl_image * expect = hdrl_image_new(10, 10);
384 hdrl_image_add_scalar(expect, (hdrl_value){101, 9.0045});
385 hdrl_image_set_pixel(expect, 3, 4, (hdrl_value){101, 9.29448});
386 hdrl_test_image_abs(coef0, expect, HDRL_EPS_DATA * 1e11);
387 hdrl_image_delete(expect);
388
389 hdrl_image * coef1 = hdrl_imagelist_get(out_coef, 1);
390 expect = hdrl_image_new(10, 10);
391 hdrl_image_add_scalar(expect, (hdrl_value){0.49, 0.0351317});
392 hdrl_image_set_pixel(expect, 3, 4, (hdrl_value){0.49, 0.0399607});
393 hdrl_test_image_abs(coef1, expect, HDRL_EPS_DATA * 1e11);
394 hdrl_image_delete(expect);
395
396 cpl_image * cexpect = cpl_image_new(10, 10, HDRL_TYPE_DATA);
397 /* dof*scaling between weights and measured error */
398 cpl_image_add_scalar(cexpect, 1.831e-29);
399 cpl_test_image_abs(out_chi2, cexpect, DBL_EPSILON * 1e9);
400 cpl_image_delete(cexpect);
401
402 cexpect = cpl_image_new(10, 10, HDRL_TYPE_DATA);
403 cpl_image_add_scalar(cexpect, 3);
404 cpl_image_set(cexpect, 3, 4, 2);
405 cpl_test_image_abs(out_dof, cexpect, 0);
406 cpl_image_delete(cexpect);
407
408 hdrl_imagelist_delete(out_coef);
409 cpl_image_delete(out_chi2);
410 cpl_image_delete(out_dof);
411
412
413 /* smoke test less good pixels than fit degree */
414 hdrl_image_reject(hdrl_imagelist_get(input, 0), 2, 2);
415 hdrl_image_reject(hdrl_imagelist_get(input, 1), 2, 2);
416 hdrl_image_reject(hdrl_imagelist_get(input, 2), 2, 2);
417 hdrl_image_reject(hdrl_imagelist_get(input, 3), 2, 2);
418 hdrl_fit_polynomial_imagelist2(input, samp, 3, &out_coef, &out_chi2, &out_dof);
419
421 cpl_imagelist_delete(samp);
422 hdrl_imagelist_delete(out_coef);
423 cpl_image_delete(out_chi2);
424 cpl_image_delete(out_dof);
425 }
426}
427
428void test_real_data(void)
429{
430 /* pixel from vcam ramp data with poisson model error (gain 2.4, ron 10) */
431 double x[] = { 3., 3., 5., 5., 7., 7., 10., 10., 12., 12.,
432 15., 15., 20., 20. };
433 double y[] = { 3441, 3420, 5606, 5586, 7814, 7815, 11003, 10970,
434 13292, 13198, 16347, 16175, 21267, 21318 };
435 double e[] = { 39.16312027, 39.05124664, 49.35416031, 49.26966476,
436 57.92955399, 57.93315125, 68.4440155 , 68.34349823,
437 75.08883667, 74.82757568, 83.13392639, 82.7017746 ,
438 94.66387939, 94.77605438 };
439 cpl_vector * vx = cpl_vector_wrap(ARRAY_LEN(x), x);
440 cpl_vector * vy = cpl_vector_wrap(ARRAY_LEN(y), y);
441 cpl_vector * ve = cpl_vector_wrap(ARRAY_LEN(e), e);
442
443 double exp_c[] = {296.10245659, 1063.12005477};
444 cpl_matrix * design = polyvander1d(vx, 1);
445 hdrl_ls_fit_result * res = polyfit1d(vx, vy, ve, 1);
446 matrix_eq(res->coef, exp_c, DBL_EPSILON * 2e10);
447
448 cpl_vector_unwrap(vx);
449 cpl_vector_unwrap(vy);
450 cpl_vector_unwrap(ve);
451 cpl_matrix_delete(design);
452 hdrl_ls_fit_result_delete(res);
453}
454
455/*----------------------------------------------------------------------------*/
459/*----------------------------------------------------------------------------*/
460int main(void)
461{
462 cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_WARNING);
463
464 test_vander1d();
465 test_fit();
466 test_poisson();
467 test_imglistfit();
468 test_imglistfit2();
469 test_real_data();
470
471 return cpl_test_end(0);
472}
cpl_error_code hdrl_fit_polynomial_imagelist2(const hdrl_imagelist *list, const cpl_imagelist *samplepos, const int degree, hdrl_imagelist **coef, cpl_image **chi2, cpl_image **dof)
weighted least squares polynomial fit of each pixel of a imagelist
Definition hdrl_fit.c:540
cpl_error_code hdrl_fit_polynomial_imagelist(const hdrl_imagelist *list, const cpl_vector *samplepos, const int degree, hdrl_imagelist **coef, cpl_image **chi2, cpl_image **dof)
weighted least squares polynomial fit of each pixel of a imagelist
Definition hdrl_fit.c:367
cpl_error_code hdrl_image_set_pixel(hdrl_image *self, cpl_size xpos, cpl_size ypos, hdrl_value value)
set pixel values of hdrl_image
Definition hdrl_image.c:594
cpl_error_code hdrl_image_add_scalar(hdrl_image *self, hdrl_value value)
Elementwise addition of a scalar to an image.
cpl_error_code hdrl_image_reject(hdrl_image *self, cpl_size xpos, cpl_size ypos)
mark pixel as bad
Definition hdrl_image.c:427
hdrl_image * hdrl_image_new(cpl_size nx, cpl_size ny)
create new zero filled hdrl image
Definition hdrl_image.c:311
void hdrl_image_delete(hdrl_image *himg)
delete hdrl_image
Definition hdrl_image.c:379
cpl_error_code hdrl_imagelist_set(hdrl_imagelist *himlist, hdrl_image *himg, cpl_size pos)
Insert an image into an imagelist.
void hdrl_imagelist_delete(hdrl_imagelist *himlist)
Free all memory used by a hdrl_imagelist object including the images.
cpl_size hdrl_imagelist_get_size(const hdrl_imagelist *himlist)
Get the number of images in the imagelist.
hdrl_imagelist * hdrl_imagelist_new(void)
Create an empty imagelist.
hdrl_image * hdrl_imagelist_get(const hdrl_imagelist *himlist, cpl_size inum)
Get an image from a list of images.