29#include "hdrl_imagelist.h"
37#define ARRAY_LEN(a) sizeof((a))/sizeof((a)[0])
46#define matrix_eq(m, exp, eps) \
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); \
56static void test_vander1d(
void)
59 double p[] = {1, 2, 3, 4};
60 double exp[] = {1., 1., 1.,
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);
70 v = polyvander1d(pv, 2);
71 cpl_test_error(CPL_ERROR_NONE);
72 matrix_eq(v, exp, DBL_EPSILON * 5);
75 cpl_vector_unwrap(pv);
80 double p[] = {1, 2, 3, 4};
81 double e[] = {1, 1, 4, 1};
82 double exp[] = {1., 1., 1.,
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);
91 matrix_eq(v, exp, DBL_EPSILON * 5);
92 cpl_vector_unwrap(pv);
93 cpl_vector_unwrap(ev);
99static void test_fit(
void)
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);
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]));
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);
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);
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);
155 hdrl_ls_fit_result_delete(r);
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);
162 cpl_vector_unwrap(pv);
163 cpl_vector_unwrap(ev);
164 cpl_vector_unwrap(sv);
183static inline cpl_vector * hdrl_ls_fit_result_get_fit_interval(
184 const hdrl_ls_fit_result * r,
185 const cpl_vector * data,
189 double mse = hdrl_ls_fit_result_get_chi2(r, data, errors) /
190 hdrl_ls_fit_result_get_residual_dof(r);
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);
204void test_poisson(
void)
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. };
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);
215 cpl_vector * ve_model = cpl_vector_duplicate(vx);
216 cpl_vector_power(ve_model, 0.5);
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);
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);
228 cpl_vector * pred_e =
229 hdrl_ls_fit_result_get_fit_interval(res, vy, ve_model);
231 cpl_test_vector_abs(pred_e, ve_real, 0.7);
232 cpl_vector_delete(pred_e);
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);
242void test_imglistfit(
void)
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;
252 cpl_test_error(CPL_ERROR_NULL_INPUT);
255 cpl_test_error(CPL_ERROR_NULL_INPUT);
258 cpl_test_error(CPL_ERROR_NULL_INPUT);
261 cpl_test_error(CPL_ERROR_NULL_INPUT);
264 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
267 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
269 for (
size_t i = 0; i < n; i++) {
270 double t = (i + 1) * 100;
274 (hdrl_value){0.5 * t + 100 - i, sqrt(0.5 * t)});
279 cpl_vector_set(samp, i, t);
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);
292 hdrl_test_image_abs(coef0, expect, HDRL_EPS_DATA * 1e11);
299 hdrl_test_image_abs(coef1, expect, HDRL_EPS_DATA * 1e11);
302 cpl_image * cexpect = cpl_image_new(10, 10, HDRL_TYPE_DATA);
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);
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);
315 cpl_image_delete(out_chi2);
316 cpl_image_delete(out_dof);
327 cpl_vector_delete(samp);
329 cpl_image_delete(out_chi2);
330 cpl_image_delete(out_dof);
334void test_imglistfit2(
void)
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;
344 cpl_test_error(CPL_ERROR_NULL_INPUT);
347 cpl_test_error(CPL_ERROR_NULL_INPUT);
350 cpl_test_error(CPL_ERROR_NULL_INPUT);
353 cpl_test_error(CPL_ERROR_NULL_INPUT);
356 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
359 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
361 for (
size_t i = 0; i < n; i++) {
362 double t = (i + 1) * 100;
366 (hdrl_value){0.5 * t + 100 - i, sqrt(0.5 * t)});
368 cpl_image * sampi = cpl_image_new(10, 10, HDRL_TYPE_DATA);
369 cpl_image_add_scalar(sampi, t);
371 cpl_image_reject(sampi, 3, 4);
373 cpl_imagelist_set(samp, sampi, i);
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);
386 hdrl_test_image_abs(coef0, expect, HDRL_EPS_DATA * 1e11);
393 hdrl_test_image_abs(coef1, expect, HDRL_EPS_DATA * 1e11);
396 cpl_image * cexpect = cpl_image_new(10, 10, HDRL_TYPE_DATA);
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);
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);
409 cpl_image_delete(out_chi2);
410 cpl_image_delete(out_dof);
421 cpl_imagelist_delete(samp);
423 cpl_image_delete(out_chi2);
424 cpl_image_delete(out_dof);
428void test_real_data(
void)
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);
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);
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);
462 cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_WARNING);
471 return cpl_test_end(0);
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
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
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
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
hdrl_image * hdrl_image_new(cpl_size nx, cpl_size ny)
create new zero filled hdrl image
void hdrl_image_delete(hdrl_image *himg)
delete hdrl_image
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.