24#ifndef HDRL_USE_PRIVATE
26#define HDRL_USE_PRIVATE
29#include "hdrl_fringe.h"
30#include "hdrl_prototyping.h"
31#include "hdrl_random.h"
45void test_hermite_sum(
void);
46void hdrl_fringe_compute_test_input(
void);
47void hdrl_fringe_correct_test_input(
void);
48void hdrl_fringe_hermite_test(
void);
49void hdrl_mime_fringe_amplitudes_test(
void);
50void hdrl_mime_fringe_amplitudes_ls_test(
void);
52static cpl_matrix * hdrl_mime_hermite_functions_create(
int n,
double center,
53 double scale,
const cpl_matrix * x);
56void hdrl_fringe_hermite_test(
void)
88 w = (b - a) / (n - 1);
91 f = cpl_matrix_new(n, 1);
93 for (i = 0; i < n; i++)
95 double xi = (cpl_matrix_get_data(x))[i];
96 (cpl_matrix_get_data(f))[i] = (1.0 + xi) * exp(-0.5 * xi * xi);
100 funs = hdrl_mime_hermite_functions_create(nfun, center, scale, NULL);
102 cpl_test_error(CPL_ERROR_NULL_INPUT);
103 funs = hdrl_mime_hermite_functions_create(0, center, scale, x);
105 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
106 funs = hdrl_mime_hermite_functions_create(nfun, center, 0.0, x);
108 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
109 funs = hdrl_mime_hermite_functions_create(nfun, center, -1.0, x);
111 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
114 funs = hdrl_mime_hermite_functions_create(nfun, center, scale, x);
115 cpl_test_nonnull(funs);
116 cpl_test_error(CPL_ERROR_NONE);
120 cpl_test_nonnull(coeffs);
121 cpl_test_error(CPL_ERROR_NONE);
122 cpl_test_eq(cpl_matrix_multiply_scalar(coeffs, w), CPL_ERROR_NONE);
126 hseries = hdrl_mime_hermite_series_create(nfun, center, scale, NULL, x);
127 cpl_test_null(hseries);
128 cpl_test_error(CPL_ERROR_NULL_INPUT);
129 hseries = hdrl_mime_hermite_series_create(nfun, center, scale, coeffs, NULL);
130 cpl_test_null(hseries);
131 cpl_test_error(CPL_ERROR_NULL_INPUT);
132 hseries = hdrl_mime_hermite_series_create(nfun, center, scale, NULL, NULL);
133 cpl_test_null(hseries);
134 cpl_test_error(CPL_ERROR_NULL_INPUT);
135 hseries = hdrl_mime_hermite_series_create(0, center, scale, coeffs, x);
136 cpl_test_null(hseries);
137 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
138 hseries = hdrl_mime_hermite_series_create(nfun, center, -1.0, coeffs, x);
139 cpl_test_null(hseries);
140 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
141 hseries = hdrl_mime_hermite_series_create(nfun, center, 0.0, coeffs, x);
142 cpl_test_null(hseries);
143 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
146 hseries = hdrl_mime_hermite_series_create(nfun, center, scale, coeffs, x);
147 cpl_test_nonnull(hseries);
148 cpl_test_error(CPL_ERROR_NONE);
149 cpl_test_eq(cpl_matrix_subtract(hseries, f), CPL_ERROR_NONE);
152 cpl_matrix_delete(x);
153 cpl_matrix_delete(f);
154 cpl_matrix_delete(funs);
155 cpl_matrix_delete(coeffs);
156 cpl_matrix_delete(hseries);
160void test_fringe_mime_gmx1(
void)
170 params[2] = sqrt(0.5);
174 params[5] = sqrt(0.5);
178 hdrl_mime_gmix1(x, params, &y);
179 cpl_test_leq(fabs(y - 3.678794411714423e-01), 1e-15);
182 hdrl_mime_gmix1(x, params, &y);
183 cpl_test_leq(fabs(y - 1.0), 1e-15);
186 hdrl_mime_gmix1(x, params, &y);
187 cpl_test_leq(fabs(y - 3.678794411714423e-01), 1e-15);
189 hdrl_mime_gmix_derivs1(x, params, result);
190 cpl_test_leq(fabs(result[0] - 3.678794411714423e-01), 1e-15);
191 cpl_test_leq(fabs(result[1] - 7.357588823428847e-01), 1e-15);
192 cpl_test_leq(fabs(result[2] - 1.040520190045778e+00), 1e-15);
193 cpl_test_leq(fabs(result[3] - 3.678794411714423e-01), 1e-15);
194 cpl_test_zero(result[4]);
195 cpl_test_zero(result[5]);
200 params[2] = sqrt(0.5);
204 params[5] = sqrt(0.5);
208 hdrl_mime_gmix1(x, params, &y);
209 cpl_test_leq(fabs(y - 3.678794411714423e-01), 1e-15);
212 hdrl_mime_gmix1(x, params, &y);
213 cpl_test_leq(fabs(y - 1.0), 1e-15);
216 hdrl_mime_gmix1(x, params, &y);
217 cpl_test_leq(fabs(y - 3.678794411714423e-01), 1e-15);
219 hdrl_mime_gmix_derivs1(x, params, result);
220 cpl_test_leq(fabs(result[0] - 3.678794411714423e-01), 1e-15);
221 cpl_test_zero(result[1]);
222 cpl_test_zero(result[2]);
223 cpl_test_leq(fabs(result[3] - 3.678794411714423e-01), 1e-15);
224 cpl_test_leq(fabs(result[4] - 7.357588823428847e-01), 1e-15);
225 cpl_test_leq(fabs(result[5] - 1.040520190045778e+00), 1e-15);
248static cpl_matrix * hdrl_mime_hermite_functions_create(
int n,
double center,
249 double scale,
const cpl_matrix * x)
259 cpl_ensure(x, CPL_ERROR_NULL_INPUT, NULL);
261 if (n < 1 || scale <= 0.0)
263 cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT);
268 nr = cpl_matrix_get_nrow(x) * cpl_matrix_get_ncol(x);
271 funs = cpl_matrix_new(nr, n);
274 rt = 1.0 / sqrt(sqrt(CPL_MATH_PI));
277 m = cpl_matrix_get_data(funs);
278 mx = cpl_matrix_get_data_const(x);
279 for (i = 0; i < nr; i++, m += n)
281 xi = (mx[i] - center) / scale;
282 m[0] = rt * exp(-0.5 * xi * xi);
286 m = cpl_matrix_get_data(funs);
287 mx = cpl_matrix_get_data_const(x);
289 for (i = 0; i < nr; i++, m += n)
291 xi = (mx[i] - center) / scale;
292 m[1] = rt * sqrt(2.0) * xi * exp(-0.5 * xi * xi);
296 m = cpl_matrix_get_data(funs);
297 for (i = 0; i < nr; i++, m += n)
299 xi = (mx[i] - center) / scale;
300 for (j = 2; j < n; j++)
302 m[j] = sqrt(2.0) * xi * m[j - 1] - sqrt(j - 1) * m[j - 2];
303 m[j] = m[j] / sqrt(j);
308 cpl_matrix_multiply_scalar(funs, 1 / sqrt(scale));
314void test_hermite_sum(
void)
320 double center, scale;
343 sums = hdrl_mime_hermite_functions_sums_create(nfun, center, scale, NULL);
345 cpl_test_error(CPL_ERROR_NULL_INPUT);
346 sums = hdrl_mime_hermite_functions_sums_create(0, center, scale, x);
348 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
349 sums = hdrl_mime_hermite_functions_sums_create(nfun, center, 0.0, x);
351 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
352 sums = hdrl_mime_hermite_functions_sums_create(nfun, center, -1.0, x);
354 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
357 sums = hdrl_mime_hermite_functions_sums_create(nfun, center, scale, x);
358 cpl_test_nonnull(sums);
359 cpl_test_error(CPL_ERROR_NONE);
361 cpl_test_leq(fabs(cpl_matrix_get(sums, 0, 0) - 1.685081590066050e+00),
363 cpl_test_leq(fabs(cpl_matrix_get(sums, 1, 0) - 9.093843925414908e-01),
365 cpl_test_leq(fabs(cpl_matrix_get(sums, 2, 0) - 4.521636055506448e-01),
367 cpl_test_leq(fabs(cpl_matrix_get(sums, 3, 0) - 8.130124958110769e-01),
369 cpl_test_leq(fabs(cpl_matrix_get(sums, 4, 0) - 8.013868548156017e-01),
373 cpl_matrix_delete(x);
374 cpl_matrix_delete(sums);
380void hdrl_fringe_compute_test_input(
void)
384 hdrl_parameter * collapse_params;
388 cpl_imagelist * ilist_obj = cpl_imagelist_new();
389 cpl_mask * stat_mask = cpl_mask_new(nx, ny);
392 cpl_image * contrib_map;
399 cpl_image * cima1 = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
400 cpl_image * cima2 = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
404 cpl_image_add_scalar(cima1, 1.);
405 cpl_image_add_scalar(cima2, 10.);
409 &master, &contrib_map, NULL);
410 cpl_test_error(CPL_ERROR_NULL_INPUT);
411 cpl_test_null(master);
412 cpl_test_null(contrib_map);
414 cpl_table *qctable = (cpl_table*)0xbadbeef;
416 &master, &contrib_map, &qctable);
417 cpl_test_error(CPL_ERROR_NULL_INPUT);
418 cpl_test_null(master);
419 cpl_test_null(contrib_map);
422 &master, &contrib_map, NULL);
423 cpl_test_error(CPL_ERROR_NULL_INPUT);
424 cpl_test_null(master);
425 cpl_test_null(contrib_map);
429 cpl_imagelist_set(ilist_obj, cima1, 0);
431 &master, &contrib_map, NULL);
432 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
433 cpl_test_null(master);
434 cpl_test_null(contrib_map);
435 cpl_imagelist_set(ilist_obj, cima2, 1);
441 &master, &contrib_map, NULL);
442 cpl_test_error(CPL_ERROR_NULL_INPUT);
443 cpl_test_null(master);
444 cpl_test_null(contrib_map);
451 &master, &contrib_map, NULL);
452 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
453 cpl_test_null(master);
454 cpl_test_null(contrib_map);
480 &master, &contrib_map, NULL);
481 cpl_test_error(CPL_ERROR_NONE);
482 cpl_test_nonnull(master);
483 cpl_test_nonnull(contrib_map);
485 cpl_image_delete(contrib_map);
491 cpl_imagelist_delete(ilist_obj);
493 cpl_mask_delete(stat_mask);
497void hdrl_fringe_correct_test_input(
void)
503 cpl_imagelist * ilist_obj = cpl_imagelist_new();
504 cpl_mask * stat_mask = cpl_mask_new(nx, ny);
515 cpl_image * cima1 = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
516 cpl_image * cima2 = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
520 cpl_image_add_scalar(cima1, 1.);
521 cpl_image_add_scalar(cima2, 10.);
525 cpl_test_error(CPL_ERROR_NULL_INPUT);
529 cpl_imagelist_set(ilist_obj, cima1, 0);
531 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
532 cpl_imagelist_set(ilist_obj, cima2, 1);
538 cpl_test_error(CPL_ERROR_NULL_INPUT);
545 cpl_test_error(CPL_ERROR_INCOMPATIBLE_INPUT);
554 cpl_test_error(CPL_ERROR_NONE);
557 cpl_test_error(CPL_ERROR_NONE);
560 cpl_test_error(CPL_ERROR_NONE);
566 hdrl_random_state * rng = hdrl_random_state_new(1, NULL);
571 for(cpl_size i = 1; i <= sizex; i++) {
572 for(cpl_size j = 1; j <= sizey; j++) {
573 cpl_image_set(cplima1, i, j, (
double)hdrl_random_poisson(rng, 100.));
574 cpl_image_set(cplima2, i, j, (
double)hdrl_random_poisson(rng, 200.));
578 hdrl_random_state_delete(rng);
579 float* pobj = cpl_image_get_data_float(cima1);
580 int sx = cpl_image_get_size_x(cima1);
581 int sy = cpl_image_get_size_y(cima1);
582 int i_min = 0.25 * sx;
583 int i_max = 0.75 * sx;
584 int j_min = 0.25 * sy;
585 int j_max = 0.75 * sy;
587 for(
int j = j_min; j < j_max; j++) {
588 for(
int i = i_min; i < i_max; i++) {
593 pobj = cpl_image_get_data_float(cima2);
594 for(
int j = j_min; j < j_max; j++) {
595 for(
int i = i_min; i < i_max; i++) {
604 cpl_binary* pmask = cpl_mask_get_data(stat_mask);
606 for(
int j = j_min; j < j_max; j++) {
607 for(
int i = i_min; i < i_max; i++) {
613 cpl_image* contrib_map=NULL;
614 hdrl_parameter * collapse_params;
618 &masterfringe, &contrib_map, NULL);
620 cpl_test_error(CPL_ERROR_NONE);
626 cpl_imagelist_delete(ilist_obj);
628 cpl_image_delete(contrib_map);
630 cpl_mask_delete(stat_mask);
636void hdrl_mime_fringe_amplitudes_test(
void)
639 cpl_image * image_double = cpl_image_new(10, 10, CPL_TYPE_DOUBLE);
640 cpl_image * image_float = cpl_image_new(10, 10, CPL_TYPE_FLOAT);
641 cpl_mask * mask = cpl_mask_new(10, 10);
645 matrix = hdrl_mime_fringe_amplitudes(NULL, NULL);
646 cpl_test_error(CPL_ERROR_NULL_INPUT);
647 cpl_test_null(matrix);
649 matrix = hdrl_mime_fringe_amplitudes(image_double, NULL);
650 cpl_test_error(CPL_ERROR_NULL_INPUT);
651 cpl_test_null(matrix);
653 matrix = hdrl_mime_fringe_amplitudes(NULL, mask);
654 cpl_test_error(CPL_ERROR_NULL_INPUT);
655 cpl_test_null(matrix);
658 matrix = hdrl_mime_fringe_amplitudes(image_float, mask);
659 cpl_test_error(CPL_ERROR_INVALID_TYPE);
660 cpl_test_null(matrix);
664 cpl_test_error(CPL_ERROR_NONE);
665 matrix = hdrl_mime_fringe_amplitudes(image_double, mask);
666 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
667 cpl_test_null(matrix);
669 cpl_mask_delete(mask);
670 cpl_image_delete(image_float);
671 cpl_image_delete(image_double);
675void hdrl_mime_fringe_amplitudes_ls_test(
void)
678 cpl_image * image_double1 = cpl_image_new(10, 10, CPL_TYPE_DOUBLE);
679 cpl_image * image_double2 = cpl_image_new(10, 10, CPL_TYPE_DOUBLE);
680 cpl_image * image_float = cpl_image_new(10, 10, CPL_TYPE_FLOAT);
681 cpl_mask * mask = cpl_mask_new(10, 10);
685 matrix = hdrl_mime_fringe_amplitudes_ls(NULL, NULL, NULL);
686 cpl_test_error(CPL_ERROR_NULL_INPUT);
687 cpl_test_null(matrix);
689 matrix = hdrl_mime_fringe_amplitudes_ls(NULL, mask, image_double2);
690 cpl_test_error(CPL_ERROR_NULL_INPUT);
691 cpl_test_null(matrix);
693 matrix = hdrl_mime_fringe_amplitudes_ls(image_double1, NULL, image_double2);
694 cpl_test_error(CPL_ERROR_NULL_INPUT);
695 cpl_test_null(matrix);
697 matrix = hdrl_mime_fringe_amplitudes_ls(image_double1, mask, NULL);
698 cpl_test_error(CPL_ERROR_NULL_INPUT);
699 cpl_test_null(matrix);
702 matrix = hdrl_mime_fringe_amplitudes_ls(image_float, mask, image_double2);
703 cpl_test_error(CPL_ERROR_INVALID_TYPE);
704 cpl_test_null(matrix);
706 matrix = hdrl_mime_fringe_amplitudes_ls(image_double1, mask, image_float);
707 cpl_test_error(CPL_ERROR_INVALID_TYPE);
708 cpl_test_null(matrix);
712 cpl_test_error(CPL_ERROR_NONE);
713 matrix = hdrl_mime_fringe_amplitudes_ls(image_double1, mask, image_double2);
714 cpl_test_error(CPL_ERROR_ILLEGAL_INPUT);
715 cpl_test_null(matrix);
717 cpl_mask_delete(mask);
718 cpl_image_delete(image_float);
719 cpl_image_delete(image_double2);
720 cpl_image_delete(image_double1);
732 cpl_test_init(PACKAGE_BUGREPORT, CPL_MSG_WARNING);
734 hdrl_fringe_hermite_test();
735 test_fringe_mime_gmx1();
737 hdrl_fringe_compute_test_input();
738 hdrl_fringe_correct_test_input();
739 hdrl_mime_fringe_amplitudes_test();
740 hdrl_mime_fringe_amplitudes_ls_test();
743 return cpl_test_end(0);
hdrl_parameter * hdrl_collapse_mean_parameter_create(void)
create a parameter object for mean
cpl_error_code hdrl_fringe_correct(hdrl_imagelist *ilist_fringe, const cpl_imagelist *ilist_obj, const cpl_mask *stat_mask, const hdrl_image *masterfringe, cpl_table **qctable)
Scales and subtracts the master fringe from the images.
cpl_error_code hdrl_fringe_compute(hdrl_imagelist *ilist_fringe, const cpl_imagelist *ilist_obj, const cpl_mask *stat_mask, const hdrl_parameter *collapse_params, hdrl_image **master, cpl_image **contrib_map, cpl_table **qctable)
Calculates the master fringe and contribution map based on the.
cpl_error_code hdrl_image_add_scalar(hdrl_image *self, hdrl_value value)
Elementwise addition of a scalar to an image.
cpl_size hdrl_image_get_size_y(const hdrl_image *self)
return size of Y dimension of image
cpl_size hdrl_image_get_size_x(const hdrl_image *self)
return size of X dimension of image
cpl_image * hdrl_image_get_image(hdrl_image *himg)
get data as cpl image
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
hdrl_image * hdrl_imagelist_unset(hdrl_imagelist *himlist, cpl_size pos)
Remove an image from an imagelist.
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.
hdrl_imagelist * hdrl_imagelist_new(void)
Create an empty imagelist.
void hdrl_parameter_delete(hdrl_parameter *obj)
shallow delete of a parameter
cpl_matrix * hdrl_mime_matrix_product_left_transpose_create(const cpl_matrix *mat1, const cpl_matrix *mat2)
Create the product of the transpose of a matrix with another matrix.
cpl_matrix * hdrl_mime_matrix_linspace_create(int n, double a, double b)
Create equally spaced nodes.