CR2RE Pipeline Reference Manual 1.7.0
hdrl_collapse.c
1/*
2 * This file is part of the HDRL
3 * Copyright (C) 2013 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_collapse.h"
29#include "hdrl_parameter.h"
30#include "hdrl_parameter_defs.h"
31#include "hdrl_sigclip.h"
32#include "hdrl_mode.h"
33#include "hdrl_utils.h"
34#include "hdrl_types.h"
35#include <cpl.h>
36
37#include <string.h>
38#include <math.h>
39/*-----------------------------------------------------------------------------
40 Types
41-----------------------------------------------------------------------------*/
42
43/* function performing imglist -> image reduction */
44typedef cpl_error_code (hdrl_collapse_imagelist_to_image_f)(
45 const cpl_imagelist * data,
46 const cpl_imagelist * errors,
47 cpl_image ** out,
48 cpl_image ** err,
49 cpl_image ** contrib,
50 void * parameters,
51 void * extra_out);
52
53struct hdrl_collapse_imagelist_to_image_s {
54 /* function performing imglist -> image reduction */
55 hdrl_collapse_imagelist_to_image_f * func;
56 /* create extra out storage */
57 void * (*create_eout)(const cpl_image *);
58 /* move extra out storage */
59 cpl_error_code (*move_eout)(void *, void *, const cpl_size);
60 /* unwrap extra out storage */
61 hdrl_free * unwrap_eout;
62 /* delete extra out storage */
63 void (*delete_eout)(void *);
64 /* parameters the reduction function requires */
65 hdrl_parameter * parameters;
66};
67
68/* function performing imglist -> image reduction */
69typedef cpl_error_code (hdrl_collapse_imagelist_to_vector_f)(
70 const cpl_imagelist * data,
71 const cpl_imagelist * errors,
72 cpl_vector ** out,
73 cpl_vector ** err,
74 cpl_array ** contrib,
75 void * parameters,
76 void * extra_out);
77
78struct hdrl_collapse_imagelist_to_vector_s {
79 /* function performing imglist -> vector reduction */
80 hdrl_collapse_imagelist_to_vector_f * func;
81 /* create extra out storage */
82 void * (*create_eout)(cpl_size);
83 /* move extra out storage */
84 cpl_error_code (*move_eout)(void *, void *, const cpl_size);
85 /* unwrap extra out storage */
86 hdrl_free * unwrap_eout;
87 /* delete extra out storage */
88 void (*delete_eout)(void *);
89 /* parameters the reduction function requires */
90 hdrl_parameter * parameters;
91};
92
93static cpl_error_code
94hdrl_collapse_mean(const cpl_imagelist * data,
95 const cpl_imagelist * errors,
96 cpl_image ** out, cpl_image ** err,
97 cpl_image ** contrib, void *, void *);
98static cpl_error_code
99hdrl_collapse_median(const cpl_imagelist * data,
100 const cpl_imagelist * errors,
101 cpl_image ** out, cpl_image ** err,
102 cpl_image ** contrib, void *, void *);
103static cpl_error_code
104hdrl_collapse_sigclip(const cpl_imagelist * data,
105 const cpl_imagelist * errors,
106 cpl_image ** out, cpl_image ** err,
107 cpl_image ** contrib,
108 void * parameters, void * extra_out);
109
110static cpl_error_code
111hdrl_collapse_minmax(const cpl_imagelist * data,
112 const cpl_imagelist * errors,
113 cpl_image ** out, cpl_image ** err,
114 cpl_image ** contrib,
115 void * parameters, void * extra_out);
116static cpl_error_code
117hdrl_collapse_mode(const cpl_imagelist * data,
118 const cpl_imagelist * errors,
119 cpl_image ** out, cpl_image ** err,
120 cpl_image ** contrib, void *, void *);
121
122static cpl_error_code
123hdrl_collapse_weighted_mean(const cpl_imagelist * data_,
124 const cpl_imagelist * errors_,
125 cpl_image ** out, cpl_image ** err,
126 cpl_image ** contrib, void *, void *);
127
128/*-----------------------------------------------------------------------------
129 Static
130 -----------------------------------------------------------------------------*/
131
132/*----------------------------------------------------------------------------*/
139/*----------------------------------------------------------------------------*/
140
142
144
145
146/*-----------------------------------------------------------------------------
147 Collapse Parameters
148 -----------------------------------------------------------------------------*/
149static void hdrl_nop_free(void * HDRL_UNUSED(x))
150{
151 return;
152}
153
154static void * hdrl_mean_alloc(size_t n);
155static void * hdrl_median_alloc(size_t n);
156static void * hdrl_weighted_mean_alloc(size_t n);
157
158/* Mean COLLAPSE */
159typedef hdrl_parameter_empty hdrl_collapse_mean_parameter;
160static hdrl_parameter_typeobj hdrl_collapse_mean_parameter_type = {
161 HDRL_PARAMETER_COLLAPSE_MEAN, /* type */
162 (hdrl_alloc *)&hdrl_mean_alloc, /* fp_alloc */
163 (hdrl_free *)&hdrl_nop_free, /* fp_free */
164 NULL, /* fp_destroy */
165 sizeof(hdrl_collapse_mean_parameter), /* obj_size */
166};
167
168HDRL_PARAMETER_SINGLETON(HDRL_COLLAPSE_MEAN,
169 hdrl_collapse_mean_parameter_type,
170 hdrl_mean_alloc)
171
172/* Median COLLAPSE */
173typedef hdrl_parameter_empty hdrl_collapse_median_parameter;
174static hdrl_parameter_typeobj hdrl_collapse_median_parameter_type = {
175 HDRL_PARAMETER_COLLAPSE_MEDIAN, /* type */
176 (hdrl_alloc *)&hdrl_median_alloc, /* fp_alloc */
177 (hdrl_free *)&hdrl_nop_free, /* fp_free */
178 NULL, /* fp_destroy */
179 sizeof(hdrl_collapse_median_parameter), /* obj_size */
180};
181
182HDRL_PARAMETER_SINGLETON(HDRL_COLLAPSE_MEDIAN,
183 hdrl_collapse_median_parameter_type,
184 hdrl_median_alloc)
185
186/* Weighted Mean COLLAPSE */
187typedef hdrl_parameter_empty hdrl_collapse_weighted_mean_parameter;
188static hdrl_parameter_typeobj hdrl_collapse_weighted_mean_parameter_type = {
189 HDRL_PARAMETER_COLLAPSE_WEIGHTED_MEAN, /* type */
190 (hdrl_alloc *)&hdrl_weighted_mean_alloc, /* fp_alloc */
191 (hdrl_free *)&hdrl_nop_free, /* fp_free */
192 NULL, /* fp_destroy */
193 sizeof(hdrl_collapse_weighted_mean_parameter), /* obj_size */
194};
195
196HDRL_PARAMETER_SINGLETON(HDRL_COLLAPSE_WEIGHTED_MEAN,
197 hdrl_collapse_weighted_mean_parameter_type,
198 hdrl_weighted_mean_alloc)
199
200/* Sigma-Clipping COLLAPSE */
201typedef struct {
202 HDRL_PARAMETER_HEAD;
203 double kappa_low;
204 double kappa_high;
205 int niter;
206} hdrl_collapse_sigclip_parameter;
207static hdrl_parameter_typeobj hdrl_collapse_sigclip_parameter_type = {
208 HDRL_PARAMETER_COLLAPSE_SIGCLIP, /* type */
209 (hdrl_alloc *)&cpl_malloc, /* fp_alloc */
210 (hdrl_free *)&cpl_free, /* fp_free */
211 NULL, /* fp_destroy */
212 sizeof(hdrl_collapse_sigclip_parameter), /* obj_size */
213};
214
215
216/* Minmax-Clipping COLLAPSE */
217typedef struct {
218 HDRL_PARAMETER_HEAD;
219 double nlow;
220 double nhigh;
221} hdrl_collapse_minmax_parameter;
222static hdrl_parameter_typeobj hdrl_collapse_minmax_parameter_type = {
223 HDRL_PARAMETER_COLLAPSE_MINMAX, /* type */
224 (hdrl_alloc *)&cpl_malloc, /* fp_alloc */
225 (hdrl_free *)&cpl_free, /* fp_free */
226 NULL, /* fp_destroy */
227 sizeof(hdrl_collapse_minmax_parameter), /* obj_size */
228};
229
230/* Mode COLLAPSE */
231typedef struct {
232 HDRL_PARAMETER_HEAD;
233 double histo_min;
234 double histo_max;
235 double bin_size;
236 hdrl_mode_type method;
237 cpl_size error_niter;
238} hdrl_collapse_mode_parameter;
239static hdrl_parameter_typeobj hdrl_collapse_mode_parameter_type = {
240 HDRL_PARAMETER_COLLAPSE_MODE, /* type */
241 (hdrl_alloc *)&cpl_malloc, /* fp_alloc */
242 (hdrl_free *)&cpl_free, /* fp_free */
243 NULL, /* fp_destroy */
244 sizeof(hdrl_collapse_mode_parameter), /* obj_size */
245};
246
248
249/* ---------------------------------------------------------------------------*/
254/* ---------------------------------------------------------------------------*/
255hdrl_parameter * hdrl_collapse_mean_parameter_create(void)
256{
257 hdrl_parameter_empty * p = (hdrl_parameter_empty *)
258 hdrl_parameter_new(&hdrl_collapse_mean_parameter_type);
259 return (hdrl_parameter *)p;
260}
261
262/* ---------------------------------------------------------------------------*/
267/* ---------------------------------------------------------------------------*/
268cpl_boolean hdrl_collapse_parameter_is_mean(const hdrl_parameter * self)
269{
270 return hdrl_parameter_check_type(self, &hdrl_collapse_mean_parameter_type);
271}
272
273/* ---------------------------------------------------------------------------*/
278/* ---------------------------------------------------------------------------*/
279hdrl_parameter * hdrl_collapse_median_parameter_create(void)
280{
281 hdrl_parameter_empty * p = (hdrl_parameter_empty *)
282 hdrl_parameter_new(&hdrl_collapse_median_parameter_type);
283 return (hdrl_parameter *)p;
284}
285
286/* ---------------------------------------------------------------------------*/
291/* ---------------------------------------------------------------------------*/
292cpl_boolean hdrl_collapse_parameter_is_median(const hdrl_parameter * self)
293{
294 return hdrl_parameter_check_type(self, &hdrl_collapse_median_parameter_type);
295}
296
297/* ---------------------------------------------------------------------------*/
302/* ---------------------------------------------------------------------------*/
303hdrl_parameter * hdrl_collapse_weighted_mean_parameter_create(void)
304{
305 hdrl_parameter_empty * p = (hdrl_parameter_empty *)
306 hdrl_parameter_new(&hdrl_collapse_weighted_mean_parameter_type);
307 return (hdrl_parameter *)p;
308}
309
310/* ---------------------------------------------------------------------------*/
315/* ---------------------------------------------------------------------------*/
316cpl_boolean hdrl_collapse_parameter_is_weighted_mean(
317 const hdrl_parameter * self)
318{
319 return hdrl_parameter_check_type(self,
320 &hdrl_collapse_weighted_mean_parameter_type);
321}
322
323/* ---------------------------------------------------------------------------*/
331/* ---------------------------------------------------------------------------*/
332hdrl_parameter *
333hdrl_collapse_sigclip_parameter_create(double kappa_low, double kappa_high, int niter)
334{
335 hdrl_collapse_sigclip_parameter * p = (hdrl_collapse_sigclip_parameter *)
336 hdrl_parameter_new(&hdrl_collapse_sigclip_parameter_type);
337 p->kappa_low = kappa_low;
338 p->kappa_high = kappa_high;
339 p->niter = niter;
340
341 if (hdrl_collapse_sigclip_parameter_verify((hdrl_parameter*)p) !=
342 CPL_ERROR_NONE) {
343 hdrl_parameter_delete((hdrl_parameter*)p);
344 return NULL;
345 }
346 return (hdrl_parameter *)p;
347}
348
349/* ---------------------------------------------------------------------------*/
354/* ---------------------------------------------------------------------------*/
355cpl_boolean hdrl_collapse_parameter_is_sigclip(const hdrl_parameter * self)
356{
357 return hdrl_parameter_check_type(self,
358 &hdrl_collapse_sigclip_parameter_type);
359}
360
362
363/*----------------------------------------------------------------------------*/
370/*----------------------------------------------------------------------------*/
371cpl_error_code hdrl_collapse_sigclip_parameter_verify(
372 const hdrl_parameter * param)
373{
374 cpl_error_ensure(param != NULL, CPL_ERROR_NULL_INPUT,
375 return CPL_ERROR_NULL_INPUT,
376 "NULL Collapse Sigclip Parameters");
377
378 cpl_error_ensure(hdrl_collapse_parameter_is_sigclip(param),
379 CPL_ERROR_INCOMPATIBLE_INPUT, return
380 CPL_ERROR_INCOMPATIBLE_INPUT,
381 "Not a Sigclip parameter");
382
383 const hdrl_collapse_sigclip_parameter * param_loc =
384 (const hdrl_collapse_sigclip_parameter *)param ;
385
386 cpl_error_ensure(param_loc->niter > 0, CPL_ERROR_ILLEGAL_INPUT,
387 return CPL_ERROR_ILLEGAL_INPUT,
388 "sigma-clipping iter (%d) value must be > 0",
389 param_loc->niter);
390
391 return CPL_ERROR_NONE ;
392}
393
395
396/* ---------------------------------------------------------------------------*/
402/* ---------------------------------------------------------------------------*/
403double hdrl_collapse_sigclip_parameter_get_kappa_high(const hdrl_parameter * p)
404{
405 cpl_ensure(p, CPL_ERROR_NULL_INPUT, -1);
406 cpl_ensure(hdrl_collapse_parameter_is_sigclip(p),
407 CPL_ERROR_INCOMPATIBLE_INPUT, -1);
408
409 return ((const hdrl_collapse_sigclip_parameter *)p)->kappa_high;
410}
411
412/* ---------------------------------------------------------------------------*/
418/* ---------------------------------------------------------------------------*/
419double hdrl_collapse_sigclip_parameter_get_kappa_low(const hdrl_parameter * p)
420{
421 cpl_ensure(p, CPL_ERROR_NULL_INPUT, -1);
422 cpl_ensure(hdrl_collapse_parameter_is_sigclip(p),
423 CPL_ERROR_INCOMPATIBLE_INPUT, -1);
424
425 return ((const hdrl_collapse_sigclip_parameter *)p)->kappa_low;
426}
427
428/* ---------------------------------------------------------------------------*/
434/* ---------------------------------------------------------------------------*/
435int hdrl_collapse_sigclip_parameter_get_niter(const hdrl_parameter * p)
436{
437 cpl_ensure(p, CPL_ERROR_NULL_INPUT, -1);
438 cpl_ensure(hdrl_collapse_parameter_is_sigclip(p),
439 CPL_ERROR_INCOMPATIBLE_INPUT, -1);
440
441 return ((const hdrl_collapse_sigclip_parameter *)p)->niter;
442}
443
444/* ---------------------------------------------------------------------------*/
452/* ---------------------------------------------------------------------------*/
453hdrl_parameter *
454hdrl_collapse_minmax_parameter_create(double nlow, double nhigh)
455{
456 hdrl_collapse_minmax_parameter * p = (hdrl_collapse_minmax_parameter *)
457 hdrl_parameter_new(&hdrl_collapse_minmax_parameter_type);
458 p->nlow = nlow;
459 p->nhigh = nhigh;
460
461 if (hdrl_collapse_minmax_parameter_verify((hdrl_parameter*)p) !=
462 CPL_ERROR_NONE) {
463 hdrl_parameter_delete((hdrl_parameter*)p);
464 return NULL;
465 }
466 return (hdrl_parameter *)p;
467}
468
469/* ---------------------------------------------------------------------------*/
474/* ---------------------------------------------------------------------------*/
475cpl_boolean hdrl_collapse_parameter_is_minmax(const hdrl_parameter * self)
476{
477 return hdrl_parameter_check_type(self,
478 &hdrl_collapse_minmax_parameter_type);
479}
480
482
483/*----------------------------------------------------------------------------*/
490/*----------------------------------------------------------------------------*/
491cpl_error_code hdrl_collapse_minmax_parameter_verify(
492 const hdrl_parameter * param)
493{
494
495 cpl_error_ensure(param != NULL, CPL_ERROR_NULL_INPUT,
496 return CPL_ERROR_NULL_INPUT,
497 "NULL Collapse Minmax Parameters");
498
499 cpl_error_ensure(hdrl_collapse_parameter_is_minmax(param),
500 CPL_ERROR_INCOMPATIBLE_INPUT, return
501 CPL_ERROR_INCOMPATIBLE_INPUT,
502 "Not a minmax parameter");
503
504 const hdrl_collapse_minmax_parameter * param_loc =
505 (const hdrl_collapse_minmax_parameter *)param ;
506
507 cpl_error_ensure(param_loc->nlow >= 0, CPL_ERROR_ILLEGAL_INPUT,
508 return CPL_ERROR_ILLEGAL_INPUT,
509 "nlow value (%g) must be >= 0",
510 param_loc->nlow);
511
512 cpl_error_ensure(param_loc->nhigh >= 0, CPL_ERROR_ILLEGAL_INPUT,
513 return CPL_ERROR_ILLEGAL_INPUT,
514 "nhigh value (%g) must be >= 0",
515 param_loc->nlow);
516
517 return CPL_ERROR_NONE ;
518}
519
521
522/* ---------------------------------------------------------------------------*/
528/* ---------------------------------------------------------------------------*/
529double hdrl_collapse_minmax_parameter_get_nhigh(const hdrl_parameter * p)
530{
531 cpl_ensure(p, CPL_ERROR_NULL_INPUT, -1);
532 cpl_ensure(hdrl_collapse_parameter_is_minmax(p),
533 CPL_ERROR_INCOMPATIBLE_INPUT, -1);
534
535 return ((const hdrl_collapse_minmax_parameter *)p)->nhigh;
536}
537
538/* ---------------------------------------------------------------------------*/
544/* ---------------------------------------------------------------------------*/
545double hdrl_collapse_minmax_parameter_get_nlow(const hdrl_parameter * p)
546{
547 cpl_ensure(p, CPL_ERROR_NULL_INPUT, -1);
548 cpl_ensure(hdrl_collapse_parameter_is_minmax(p),
549 CPL_ERROR_INCOMPATIBLE_INPUT, -1);
550
551 return ((const hdrl_collapse_minmax_parameter *)p)->nlow;
552}
553
554/* ---------------------------------------------------------------------------*/
565/* ---------------------------------------------------------------------------*/
566hdrl_parameter *
567hdrl_collapse_mode_parameter_create(double histo_min,
568 double histo_max,
569 double bin_size,
570 hdrl_mode_type mode_method,
571 cpl_size error_niter)
572{
573 hdrl_collapse_mode_parameter * p = (hdrl_collapse_mode_parameter *)
574 hdrl_parameter_new(&hdrl_collapse_mode_parameter_type);
575 p->histo_min = histo_min;
576 p->histo_max = histo_max;
577 p->bin_size = bin_size;
578 p->method = mode_method;
579 p->error_niter = error_niter;
580 if (hdrl_collapse_mode_parameter_verify((hdrl_parameter*)p) !=
581 CPL_ERROR_NONE) {
582 hdrl_parameter_delete((hdrl_parameter*)p);
583 return NULL;
584 }
585 return (hdrl_parameter *)p;
586}
587
588/* ---------------------------------------------------------------------------*/
593/* ---------------------------------------------------------------------------*/
594cpl_boolean hdrl_collapse_parameter_is_mode(const hdrl_parameter * self)
595{
596 return hdrl_parameter_check_type(self,
597 &hdrl_collapse_mode_parameter_type);
598}
599
601
602/*----------------------------------------------------------------------------*/
609/*----------------------------------------------------------------------------*/
610cpl_error_code hdrl_collapse_mode_parameter_verify(
611 const hdrl_parameter * param)
612{
613
614 cpl_error_ensure(param != NULL, CPL_ERROR_NULL_INPUT,
615 return CPL_ERROR_NULL_INPUT,
616 "NULL Collapse Mode Parameters");
617
618 cpl_error_ensure(hdrl_collapse_parameter_is_mode(param),
619 CPL_ERROR_INCOMPATIBLE_INPUT, return
620 CPL_ERROR_INCOMPATIBLE_INPUT,
621 "Not a mode parameter");
622
623 const hdrl_collapse_mode_parameter * param_loc =
624 (const hdrl_collapse_mode_parameter *)param ;
625
626
627 cpl_error_ensure(param_loc->bin_size >= 0, CPL_ERROR_ILLEGAL_INPUT,
628 return CPL_ERROR_ILLEGAL_INPUT,
629 "bin_size value (%g) must be >= 0",
630 param_loc->bin_size);
631
632 cpl_error_ensure(param_loc->error_niter >= 0, CPL_ERROR_ILLEGAL_INPUT,
633 return CPL_ERROR_ILLEGAL_INPUT,
634 "error_niter value (%lld) must be >= 0",
635 param_loc->error_niter);
636
637 cpl_error_ensure(param_loc->method == HDRL_MODE_MEDIAN ||
638 param_loc->method == HDRL_MODE_FIT ||
639 param_loc->method == HDRL_MODE_WEIGHTED ,
640 CPL_ERROR_ILLEGAL_INPUT,
641 return CPL_ERROR_ILLEGAL_INPUT,
642 "Please check the computation method of the mode. It has to be "
643 "%d, or %d, or %d", HDRL_MODE_MEDIAN, HDRL_MODE_WEIGHTED,
644 HDRL_MODE_FIT);
645
646 return CPL_ERROR_NONE ;
647}
648
650
651/* ---------------------------------------------------------------------------*/
657/* ---------------------------------------------------------------------------*/
658double hdrl_collapse_mode_parameter_get_histo_min(const hdrl_parameter * p)
659{
660 cpl_ensure(p, CPL_ERROR_NULL_INPUT, -1);
661 cpl_ensure(hdrl_collapse_parameter_is_mode(p),
662 CPL_ERROR_INCOMPATIBLE_INPUT, -1);
663
664 return ((const hdrl_collapse_mode_parameter *)p)->histo_min;
665}
666
667/* ---------------------------------------------------------------------------*/
673/* ---------------------------------------------------------------------------*/
674double hdrl_collapse_mode_parameter_get_histo_max(const hdrl_parameter * p)
675{
676 cpl_ensure(p, CPL_ERROR_NULL_INPUT, -1);
677 cpl_ensure(hdrl_collapse_parameter_is_mode(p),
678 CPL_ERROR_INCOMPATIBLE_INPUT, -1);
679
680 return ((const hdrl_collapse_mode_parameter *)p)->histo_max;
681}
682
683/* ---------------------------------------------------------------------------*/
689/* ---------------------------------------------------------------------------*/
690double hdrl_collapse_mode_parameter_get_bin_size(const hdrl_parameter * p)
691{
692 cpl_ensure(p, CPL_ERROR_NULL_INPUT, -1);
693 cpl_ensure(hdrl_collapse_parameter_is_mode(p),
694 CPL_ERROR_INCOMPATIBLE_INPUT, -1);
695
696 return ((const hdrl_collapse_mode_parameter *)p)->bin_size;
697}
698
699/* ---------------------------------------------------------------------------*/
705/* ---------------------------------------------------------------------------*/
706hdrl_mode_type hdrl_collapse_mode_parameter_get_method(const hdrl_parameter * p)
707{
708 cpl_ensure(p, CPL_ERROR_NULL_INPUT, -1);
709 cpl_ensure(hdrl_collapse_parameter_is_mode(p),
710 CPL_ERROR_INCOMPATIBLE_INPUT, -1);
711
712 return ((const hdrl_collapse_mode_parameter *)p)->method;
713}
714
715/* ---------------------------------------------------------------------------*/
721/* ---------------------------------------------------------------------------*/
722cpl_size hdrl_collapse_mode_parameter_get_error_niter(const hdrl_parameter * p)
723{
724 cpl_ensure(p, CPL_ERROR_NULL_INPUT, -1);
725 cpl_ensure(hdrl_collapse_parameter_is_mode(p),
726 CPL_ERROR_INCOMPATIBLE_INPUT, -1);
727
728 return ((const hdrl_collapse_mode_parameter *)p)->error_niter;
729}
730
731
732/* ---------------------------------------------------------------------------*/
747/* ---------------------------------------------------------------------------*/
748cpl_parameterlist * hdrl_collapse_parameter_create_parlist(
749 const char *base_context,
750 const char *prefix,
751 const char *method_def,
752 hdrl_parameter *sigclip_def,
753 hdrl_parameter *minmax_def,
754 hdrl_parameter *mode_def)
755{
756 cpl_ensure(base_context && prefix, CPL_ERROR_NULL_INPUT, NULL);
757
758 cpl_ensure(hdrl_collapse_parameter_is_sigclip(sigclip_def)
759 && hdrl_collapse_parameter_is_minmax( minmax_def )
760 && hdrl_collapse_parameter_is_mode( mode_def ),
761 CPL_ERROR_INCOMPATIBLE_INPUT, NULL);
762
763 char *name;
764 cpl_parameterlist *parlist = cpl_parameterlist_new();
765 cpl_parameter *p;
766 char *context = hdrl_join_string(".", 2, base_context, prefix);
767
768 /* --prefix.method */
769 name = hdrl_join_string(".", 2, context, "method");
770 p = cpl_parameter_new_enum(name, CPL_TYPE_STRING,
771 "Method used for collapsing the data", context,
772 method_def, 6, "MEAN", "WEIGHTED_MEAN", "MEDIAN", "SIGCLIP",
773 "MINMAX", "MODE");
774 cpl_free(name) ;
775 name = hdrl_join_string(".", 2, prefix, "method");
776 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
777 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
778 cpl_free(name) ;
779 cpl_parameterlist_append(parlist, p);
780
781 /* --prefix.sigclip.xxx */
782 name = hdrl_join_string(".", 2, prefix, "sigclip");
783 cpl_parameterlist * psigclip = hdrl_sigclip_parameter_create_parlist(
784 base_context, name, sigclip_def);
785 cpl_free(name) ;
786
787 for (cpl_parameter * par = cpl_parameterlist_get_first(psigclip);
788 par != NULL;
789 par = cpl_parameterlist_get_next(psigclip)) {
790 cpl_parameterlist_append(parlist, cpl_parameter_duplicate(par));
791 }
792 cpl_parameterlist_delete(psigclip);
793
794 /* --prefix.minmax.xxx */
795 name = hdrl_join_string(".", 2, prefix, "minmax");
796 cpl_parameterlist * pminmax = hdrl_minmax_parameter_create_parlist(
797 base_context, name, minmax_def);
798 cpl_free(name) ;
799
800
801 for (cpl_parameter * par = cpl_parameterlist_get_first(pminmax);
802 par != NULL;
803 par = cpl_parameterlist_get_next(pminmax)) {
804 cpl_parameterlist_append(parlist, cpl_parameter_duplicate(par));
805 }
806 cpl_parameterlist_delete(pminmax);
807
808
809 /* --prefix.mode.xxx */
810 name = hdrl_join_string(".", 2, prefix, "mode");
811 cpl_parameterlist * pmode = hdrl_mode_parameter_create_parlist(
812 base_context, name, mode_def);
813 cpl_free(name) ;
814
815
816 for (cpl_parameter * par = cpl_parameterlist_get_first(pmode);
817 par != NULL;
818 par = cpl_parameterlist_get_next(pmode)) {
819 cpl_parameterlist_append(parlist, cpl_parameter_duplicate(par));
820 }
821 cpl_parameterlist_delete(pmode);
822
823 cpl_free(context);
824
825 if (cpl_error_get_code()) {
826 cpl_parameterlist_delete(parlist);
827 return NULL;
828 }
829
830 return parlist;
831}
832
833/* ---------------------------------------------------------------------------*/
844/* ---------------------------------------------------------------------------*/
845hdrl_parameter * hdrl_collapse_parameter_parse_parlist(
846 const cpl_parameterlist * parlist,
847 const char * prefix)
848{
849 cpl_ensure(prefix && parlist, CPL_ERROR_NULL_INPUT, NULL);
850 hdrl_parameter * p = NULL;
851
852 /* Get the Method parameter */
853 char * name = hdrl_join_string(".", 2, prefix, "method");
854 const cpl_parameter * par = cpl_parameterlist_find_const(parlist, name);
855 const char * value = cpl_parameter_get_string(par);
856 if (value == NULL) {
857 cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
858 "Parameter %s not found", name);
859 cpl_free(name);
860 return NULL;
861 }
862
863 /* Switch on the methods */
864 if(!strcmp(value, "MEDIAN")) {
865 p = hdrl_collapse_median_parameter_create();
866 } else if (!strcmp(value, "WEIGHTED_MEAN")) {
867 p = hdrl_collapse_weighted_mean_parameter_create();
868 } else if (!strcmp(value, "MEAN")) {
869 p = hdrl_collapse_mean_parameter_create();
870 } else if (!strcmp(value, "SIGCLIP")) {
871 double kappa_low, kappa_high;
872 int niter;
873 hdrl_sigclip_parameter_parse_parlist(parlist, prefix, &kappa_low,
874 &kappa_high, &niter);
875 p = hdrl_collapse_sigclip_parameter_create(kappa_low,kappa_high,niter);
876 } else if (!strcmp(value, "MINMAX")) {
877 double nlow, nhigh;
878 hdrl_minmax_parameter_parse_parlist(parlist, prefix, &nlow,
879 &nhigh);
880 p = hdrl_collapse_minmax_parameter_create(nlow,nhigh);
881
882 } else if (!strcmp(value, "MODE")) {
883 double histo_min, histo_max, bin_size;
884 cpl_size error_niter;
885 hdrl_mode_type method;
886 hdrl_mode_parameter_parse_parlist(parlist, prefix, &histo_min,
887 &histo_max,
888 &bin_size,
889 &method,
890 &error_niter);
891 p = hdrl_collapse_mode_parameter_create(histo_min, histo_max, bin_size,
892 method, error_niter);
893 } else {
894 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
895 "%s not a valid method for %s", value, name);
896 cpl_free(name);
897 return NULL;
898 }
899 cpl_free(name);
900 return p;
901}
902
904
905/* ---------------------------------------------------------------------------*/
918/* ---------------------------------------------------------------------------*/
919static cpl_imagelist *
920wrap_synced_errlist(const cpl_imagelist * data, const cpl_imagelist * errors)
921{
922 cpl_imagelist * nerrors = cpl_imagelist_new();
923 for (size_t i = 0; i < (size_t)cpl_imagelist_get_size(errors); i++) {
924 const cpl_image * img = cpl_imagelist_get_const(data, i);
925 const cpl_image * err = cpl_imagelist_get_const(errors, i);
926
927 CPL_DIAG_PRAGMA_PUSH_IGN(-Wcast-qual);
928
929 cpl_image * nerr = cpl_image_wrap(cpl_image_get_size_x(err),
930 cpl_image_get_size_y(err),
931 cpl_image_get_type(err),
932 (void*)cpl_image_get_data_const(err));
933 cpl_mask_delete(hcpl_image_set_bpm(nerr,
934 (cpl_mask*)cpl_image_get_bpm_const(img)));
935
936 CPL_DIAG_PRAGMA_POP;
937
938 cpl_imagelist_set(nerrors, nerr, i);
939 }
940 return nerrors;
941}
942
943/* ---------------------------------------------------------------------------*/
949/* ---------------------------------------------------------------------------*/
950static void
951unwrap_synced_errlist(cpl_imagelist * nerrors)
952{
953 for (size_t i = 0; i < (size_t)cpl_imagelist_get_size(nerrors); i++) {
954
955 CPL_DIAG_PRAGMA_PUSH_IGN(-Wcast-qual);
956
957 cpl_image * err = (cpl_image*)cpl_imagelist_get_const(nerrors, i);
958
959 CPL_DIAG_PRAGMA_POP;
960
961 cpl_image_unset_bpm(err);
962 cpl_image_unwrap(err);
963 }
964 cpl_imagelist_unwrap(nerrors);
965}
966
967/* ---------------------------------------------------------------------------*/
979/* ---------------------------------------------------------------------------*/
980static cpl_image * imagelist_sqsum(
981 const cpl_imagelist * data,
982 cpl_image ** pcontrib)
983{
984 cpl_image * contrib = cpl_image_new_from_accepted(data);
985 cpl_image * res = NULL;
986
987 for (cpl_size i = 0; i < cpl_imagelist_get_size(data); i++) {
988 const cpl_image * img = cpl_imagelist_get_const(data, i);
989 cpl_image * sqerr = cpl_image_multiply_create(img, img);
990 if (cpl_image_get_bpm_const(sqerr)) {
991 cpl_image_fill_rejected(sqerr, 0.0);
992 cpl_image_accept_all(sqerr);
993 }
994
995 if (i == 0) {
996 res = sqerr;
997 } else {
998 cpl_image_add(res, sqerr);
999 cpl_image_delete(sqerr);
1000 }
1001 }
1002 cpl_mask * allbad = cpl_mask_threshold_image_create(contrib, 0, 0);
1003 cpl_image_reject_from_mask(res, allbad);
1004 cpl_mask_delete(allbad);
1005
1006 if (pcontrib)
1007 *pcontrib = contrib;
1008 else
1009 cpl_image_delete(contrib);
1010
1011 return res;
1012}
1013
1014/* ---------------------------------------------------------------------------*/
1040/* ---------------------------------------------------------------------------*/
1041static cpl_error_code
1042hdrl_collapse_mean(const cpl_imagelist * data,
1043 const cpl_imagelist * errors,
1044 cpl_image ** out, cpl_image ** err,
1045 cpl_image ** contrib,
1046 void * HDRL_UNUSED(parameters),
1047 void * HDRL_UNUSED(extra_out))
1048{
1049 /* (\Sum_i^n x_i) / n */
1050 /* \sqrt(\Sum_i^n x_i^2) / n */
1051 cpl_errorstate prestate = cpl_errorstate_get();
1052 *out = cpl_imagelist_collapse_create(data);
1053 /* ignore division by 0 on all pixels zero error */
1054 if (*out == NULL) {
1055 cpl_errorstate_set(prestate);
1056 *out = cpl_image_duplicate(cpl_imagelist_get_const(data, 0));
1057 cpl_image_accept_all(*out);
1058 cpl_mask_not(cpl_image_get_bpm(*out));
1059 *err = cpl_image_duplicate(cpl_imagelist_get_const(errors, 0));
1060 cpl_image_accept_all(*err);
1061 cpl_mask_not(cpl_image_get_bpm(*err));
1062 *contrib = cpl_image_new(cpl_image_get_size_x(*err),
1063 cpl_image_get_size_y(*err), CPL_TYPE_INT);
1064 cpl_image_fill_rejected(*out, NAN);
1065 cpl_image_fill_rejected(*err, NAN);
1066 }
1067 else {
1068 *err = imagelist_sqsum(errors, contrib);
1069 cpl_image_power(*err, 0.5);
1070 cpl_image_divide(*err, *contrib);
1071 cpl_image_fill_rejected(*out, NAN);
1072 cpl_image_fill_rejected(*err, NAN);
1073 }
1074
1075 return cpl_error_get_code();
1076}
1077
1078/* ---------------------------------------------------------------------------*/
1103/* ---------------------------------------------------------------------------*/
1104static cpl_error_code
1105hdrl_collapse_weighted_mean(const cpl_imagelist * data_,
1106 const cpl_imagelist * errors_,
1107 cpl_image ** out, cpl_image ** err,
1108 cpl_image ** contrib,
1109 void * HDRL_UNUSED(parameters),
1110 void * HDRL_UNUSED(extra_out))
1111{
1112 /* (\Sum_i^n w_i * x_i) / (\Sum_i^n w_i) */
1113 /* 1 / \sqrt(\Sum_i^n w_i^2) */
1114 cpl_errorstate prestate = cpl_errorstate_get();
1115 cpl_imagelist * data = cpl_imagelist_duplicate(data_);
1116 cpl_imagelist * errors = cpl_imagelist_new();
1117 cpl_image * tmperr;
1118 cpl_imagelist_cast(errors, errors_,
1119 cpl_image_get_type(cpl_imagelist_get(data, 0)));
1120 cpl_imagelist_power(errors, -2);
1121 cpl_imagelist_multiply(data, errors);
1122 *contrib = cpl_image_new_from_accepted(data);
1123 *out = cpl_imagelist_collapse_create(data);
1124 if (*out == NULL) {
1125 cpl_errorstate_set(prestate);
1126 *out = cpl_image_duplicate(cpl_imagelist_get_const(data, 0));
1127 cpl_image_accept_all(*out);
1128 cpl_mask_not(cpl_image_get_bpm(*out));
1129 *err = cpl_image_duplicate(cpl_imagelist_get_const(errors, 0));
1130 cpl_image_accept_all(*err);
1131 cpl_mask_not(cpl_image_get_bpm(*err));
1132 cpl_image_fill_rejected(*out, NAN);
1133 cpl_image_fill_rejected(*err, NAN);
1134 cpl_imagelist_delete(errors);
1135 cpl_imagelist_delete(data);
1136 return cpl_error_get_code();
1137 }
1138 cpl_imagelist_delete(data);
1139 tmperr = cpl_imagelist_collapse_create(errors);
1140 cpl_imagelist_delete(errors);
1141 cpl_image_multiply(*out, *contrib);
1142 cpl_image_multiply(tmperr, *contrib);
1143 cpl_image_divide(*out, tmperr);
1144 cpl_image_power(tmperr, -0.5);
1145 if (cpl_image_get_type(cpl_imagelist_get_const(errors_, 0)) ==
1146 cpl_image_get_type(cpl_imagelist_get_const(data_, 0))) {
1147 *err = tmperr;
1148 }
1149 else {
1150 *err = cpl_image_cast(tmperr,
1151 cpl_image_get_type(cpl_imagelist_get_const(errors_, 0)));
1152 cpl_image_delete(tmperr);
1153 }
1154
1155 cpl_image_fill_rejected(*out, NAN);
1156 cpl_image_fill_rejected(*err, NAN);
1157
1158 return cpl_error_get_code();
1159}
1160
1161/* ---------------------------------------------------------------------------*/
1178/* ---------------------------------------------------------------------------*/
1179static cpl_error_code
1180hdrl_collapse_median(const cpl_imagelist * data,
1181 const cpl_imagelist * errors,
1182 cpl_image ** out, cpl_image ** err,
1183 cpl_image ** contrib,
1184 void * HDRL_UNUSED(parameters),
1185 void * HDRL_UNUSED(extra_out))
1186{
1187 cpl_errorstate prestate = cpl_errorstate_get();
1188 /* same as mean with scaling by \sqrt(\pi / 2) */
1189 *out = cpl_imagelist_collapse_median_create(data);
1190 *err = imagelist_sqsum(errors, contrib);
1191 cpl_image_power(*err, 0.5);
1192 cpl_image_divide(*err, *contrib);
1193 if (cpl_error_get_code() == CPL_ERROR_DIVISION_BY_ZERO) {
1194 cpl_errorstate_set(prestate);
1195 cpl_image_accept_all(*out);
1196 cpl_mask_not(cpl_image_get_bpm(*out));
1197 cpl_image_accept_all(*err);
1198 cpl_mask_not(cpl_image_get_bpm(*err));
1199 cpl_image_fill_rejected(*out, NAN);
1200 cpl_image_fill_rejected(*err, NAN);
1201 return cpl_error_get_code();
1202 }
1203 /* scale error so it estimates stdev of normal distribution */
1204 cpl_image_multiply_scalar(*err, sqrt(CPL_MATH_PI_2));
1205 /* revert scaling for contrib <= 2 as median == mean in this case */
1206 cpl_image * tmp = cpl_image_cast(*contrib, CPL_TYPE_DOUBLE);
1207 cpl_image_threshold(tmp, 2.1, 2.1, 1. / sqrt(CPL_MATH_PI_2), 1.);
1208 cpl_image_multiply(*err, tmp);
1209 cpl_image_delete(tmp);
1210 cpl_image_fill_rejected(*out, NAN);
1211 cpl_image_fill_rejected(*err, NAN);
1212
1213 return cpl_error_get_code();
1214}
1215
1216/* ---------------------------------------------------------------------------*/
1233/* ---------------------------------------------------------------------------*/
1234static cpl_error_code
1235hdrl_collapse_sigclip(const cpl_imagelist * data,
1236 const cpl_imagelist * errors,
1237 cpl_image ** out, cpl_image ** err,
1238 cpl_image ** contrib,
1239 void * parameters, void * extra_out)
1240{
1241 /* same as mean but working on the not-clipped values only */
1242 hdrl_collapse_sigclip_parameter * par = parameters;
1243 hdrl_sigclip_image_output * eout =
1244 (hdrl_sigclip_image_output *)extra_out;
1245 cpl_ensure_code(par, CPL_ERROR_NULL_INPUT);
1246 const cpl_image * img = cpl_imagelist_get_const(data, 0);
1247 size_t nx = cpl_image_get_size_x(img);
1248 size_t ny = cpl_image_get_size_y(img);
1249 *out = cpl_image_new(nx, ny, HDRL_TYPE_DATA);
1250 *err = cpl_image_new(nx, ny, HDRL_TYPE_ERROR);
1251 *contrib = cpl_image_new(nx, ny, CPL_TYPE_INT);
1252 hdrl_vector_cache * cache =
1253 hdrl_vector_cache_new(cpl_imagelist_get_size(data), nx * 2);
1254
1255 /* sigmaclip along imagelist axis */
1256 for (size_t y = 1; y < ny + 1; y++) {
1257 cpl_vector * vdv[nx];
1258 cpl_vector * vev[nx];
1259 hdrl_imagelist_to_vector_row(data, y, vdv, cache);
1260 hdrl_imagelist_to_vector_row(errors, y, vev, cache);
1261 for (size_t x = 1; x < nx + 1; x++) {
1262 cpl_vector * vd = vdv[x - 1];
1263 cpl_vector * ve = vev[x - 1];
1264 if (vd && ve) {
1265 double m, e, rej_low, rej_high;
1266 cpl_size naccepted;
1267 hdrl_kappa_sigma_clip(vd, ve,
1268 par->kappa_low, par->kappa_high,
1269 par->niter, CPL_TRUE,
1270 &m, &e, &naccepted,
1271 &rej_low, &rej_high);
1272 cpl_image_set(*out, x, y, m);
1273 cpl_image_set(*err, x, y, e);
1274 cpl_image_set(*contrib, x, y, naccepted);
1275 if (eout) {
1276 cpl_image_set(eout->reject_low, x, y, rej_low);
1277 cpl_image_set(eout->reject_high, x, y, rej_high);
1278 }
1279 }
1280 else {
1281 cpl_image_set(*out, x, y, NAN);
1282 cpl_image_set(*err, x, y, NAN);
1283 cpl_image_reject(*out, x, y);
1284 cpl_image_reject(*err, x, y);
1285 cpl_image_set(*contrib, x, y, 0);
1286 if (eout) {
1287 cpl_image_set(eout->reject_low, x, y, 0.);
1288 cpl_image_set(eout->reject_high, x, y, 0.);
1289 }
1290 }
1291 hdrl_cplvector_delete_to_cache(cache, vd);
1292 hdrl_cplvector_delete_to_cache(cache, ve);
1293 }
1294 }
1295 hdrl_vector_cache_delete(cache);
1296
1297 return cpl_error_get_code();
1298}
1299
1300
1301
1302
1303/* ---------------------------------------------------------------------------*/
1321/* ---------------------------------------------------------------------------*/
1322static cpl_error_code
1323hdrl_collapse_minmax(const cpl_imagelist * data,
1324 const cpl_imagelist * errors,
1325 cpl_image ** out, cpl_image ** err,
1326 cpl_image ** contrib,
1327 void * parameters, void * extra_out)
1328{
1329 /* same as mean but working on the not-clipped values only */
1330 hdrl_collapse_minmax_parameter * par = parameters;
1331 cpl_ensure_code(par, CPL_ERROR_NULL_INPUT);
1332 hdrl_minmax_image_output * eout =
1333 (hdrl_minmax_image_output *)extra_out;
1334 const cpl_image * img = cpl_imagelist_get_const(data, 0);
1335 size_t nx = cpl_image_get_size_x(img);
1336 size_t ny = cpl_image_get_size_y(img);
1337 *out = cpl_image_new(nx, ny, HDRL_TYPE_DATA);
1338 *err = cpl_image_new(nx, ny, HDRL_TYPE_ERROR);
1339 *contrib = cpl_image_new(nx, ny, CPL_TYPE_INT);
1340 hdrl_vector_cache * cache =
1341 hdrl_vector_cache_new(cpl_imagelist_get_size(data), nx * 2);
1342
1343 /* minmaxclip along imagelist axis */
1344 for (size_t y = 1; y < ny + 1; y++) {
1345 cpl_vector * vdv[nx];
1346 cpl_vector * vev[nx];
1347 hdrl_imagelist_to_vector_row(data, y, vdv, cache);
1348 hdrl_imagelist_to_vector_row(errors, y, vev, cache);
1349 for (size_t x = 1; x < nx + 1; x++) {
1350 cpl_vector * vd = vdv[x - 1];
1351 cpl_vector * ve = vev[x - 1];
1352 if (vd && ve) {
1353 double m, e, rej_low, rej_high;
1354 cpl_size naccepted;
1355 hdrl_minmax_clip(vd, ve,
1356 par->nlow, par->nhigh, CPL_TRUE,
1357 &m, &e, &naccepted,
1358 &rej_low, &rej_high);
1359 cpl_image_set(*out, x, y, m);
1360 cpl_image_set(*err, x, y, e);
1361 cpl_image_set(*contrib, x, y, naccepted);
1362 if (eout) {
1363 cpl_image_set(eout->reject_low, x, y, rej_low);
1364 cpl_image_set(eout->reject_high, x, y, rej_high);
1365 }
1366 }
1367 else {
1368 cpl_image_set(*out, x, y, NAN);
1369 cpl_image_set(*err, x, y, NAN);
1370 cpl_image_reject(*out, x, y);
1371 cpl_image_reject(*err, x, y);
1372 cpl_image_set(*contrib, x, y, 0);
1373 if (eout) {
1374 cpl_image_set(eout->reject_low, x, y, 0.);
1375 cpl_image_set(eout->reject_high, x, y, 0.);
1376 }
1377 }
1378 hdrl_cplvector_delete_to_cache(cache, vd);
1379 hdrl_cplvector_delete_to_cache(cache, ve);
1380 }
1381 }
1382 hdrl_vector_cache_delete(cache);
1383
1384 return cpl_error_get_code();
1385}
1386
1387
1388
1389/* ---------------------------------------------------------------------------*/
1406/* ---------------------------------------------------------------------------*/
1407static cpl_error_code
1408hdrl_collapse_mode(const cpl_imagelist * data,
1409 const cpl_imagelist * errors,
1410 cpl_image ** out,
1411 cpl_image ** err,
1412 cpl_image ** contrib,
1413 void * parameters,
1414 void * HDRL_UNUSED(extra_out))
1415{
1416
1417 hdrl_collapse_mode_parameter * par = parameters;
1418 cpl_ensure_code(par, CPL_ERROR_NULL_INPUT);
1419 const cpl_image * img = cpl_imagelist_get_const(data, 0);
1420 size_t nx = cpl_image_get_size_x(img);
1421 size_t ny = cpl_image_get_size_y(img);
1422 *out = cpl_image_new(nx, ny, HDRL_TYPE_DATA);
1423 *err = cpl_image_new(nx, ny, HDRL_TYPE_ERROR);
1424 *contrib = cpl_image_new(nx, ny, CPL_TYPE_INT);
1425 hdrl_vector_cache * cache =
1426 hdrl_vector_cache_new(cpl_imagelist_get_size(data), nx * 2);
1427
1428 /* mode along imagelist axis */
1429 for (size_t y = 1; y < ny + 1; y++) {
1430 cpl_vector * vdv[nx];
1431 cpl_vector * vev[nx];
1432 hdrl_imagelist_to_vector_row(data, y, vdv, cache);
1433 hdrl_imagelist_to_vector_row(errors, y, vev, cache);
1434 for (size_t x = 1; x < nx + 1; x++) {
1435 cpl_vector * vd = vdv[x - 1];
1436 cpl_vector * ve = vev[x - 1];
1437 cpl_errorstate prestate = cpl_errorstate_get();
1438 if (vd && ve) {
1439 double m, e;
1440 cpl_size naccepted;
1441 cpl_errorstate prestate1 = cpl_errorstate_get();
1442 if (hdrl_mode_clip(vd,
1443 par->histo_min,
1444 par->histo_max,
1445 par->bin_size,
1446 par->method,
1447 par->error_niter,
1448 &m, &e, &naccepted) == CPL_ERROR_NONE) {
1449 cpl_image_set(*out, x, y, m);
1450 cpl_image_set(*err, x, y, e);
1451 cpl_image_set(*contrib, x, y, naccepted);
1452 } else {
1453 cpl_image_set(*out, x, y, NAN);
1454 cpl_image_set(*err, x, y, NAN);
1455 cpl_image_reject(*out, x, y);
1456 cpl_image_reject(*err, x, y);
1457 cpl_image_set(*contrib, x, y, 0);
1458 cpl_errorstate_set(prestate1);
1459 }
1460 }
1461 else {
1462 cpl_image_set(*out, x, y, NAN);
1463 cpl_image_set(*err, x, y, NAN);
1464 cpl_image_reject(*out, x, y);
1465 cpl_image_reject(*err, x, y);
1466 cpl_image_set(*contrib, x, y, 0);
1467 cpl_errorstate_set(prestate);
1468 }
1469 hdrl_cplvector_delete_to_cache(cache, vd);
1470 hdrl_cplvector_delete_to_cache(cache, ve);
1471 }
1472 }
1473 hdrl_vector_cache_delete(cache);
1474
1475 return cpl_error_get_code();}
1476
1477
1478static void *
1479hdrl_nop_create_eout_vec(cpl_size HDRL_UNUSED(size))
1480{
1481 return NULL;
1482}
1483
1484static void *
1485hdrl_nop_create_eout_img(const cpl_image * HDRL_UNUSED(img))
1486{
1487 return NULL;
1488}
1489
1490static cpl_error_code
1491hdrl_nop_move_eout(void * HDRL_UNUSED(dst_),
1492 void * HDRL_UNUSED(src_), const cpl_size HDRL_UNUSED(y))
1493{
1494 return CPL_ERROR_NONE;
1495}
1496
1497static void
1498hdrl_nop_unwrap_eout(void * HDRL_UNUSED(dst))
1499{
1500 return;
1501}
1502
1503/* ---------------------------------------------------------------------------*/
1508/* ---------------------------------------------------------------------------*/
1509hdrl_collapse_imagelist_to_image_t *
1510 hdrl_collapse_imagelist_to_image_mean(void)
1511{
1512 hdrl_collapse_imagelist_to_image_t * s = cpl_calloc(1, sizeof(*s));
1513 s->func = &hdrl_collapse_mean;
1514 s->create_eout = &hdrl_nop_create_eout_img;
1515 s->move_eout = &hdrl_nop_move_eout;
1516 s->unwrap_eout = &hdrl_nop_unwrap_eout;
1517 s->delete_eout = &hdrl_nop_unwrap_eout;
1518 return s;
1519}
1520
1521/* ---------------------------------------------------------------------------*/
1526/* ---------------------------------------------------------------------------*/
1527hdrl_collapse_imagelist_to_image_t *
1528 hdrl_collapse_imagelist_to_image_weighted_mean(void)
1529{
1530 hdrl_collapse_imagelist_to_image_t * s = cpl_calloc(1, sizeof(*s));
1531 s->func = &hdrl_collapse_weighted_mean;
1532 s->create_eout = &hdrl_nop_create_eout_img;
1533 s->move_eout = &hdrl_nop_move_eout;
1534 s->unwrap_eout = &hdrl_nop_unwrap_eout;
1535 s->delete_eout = &hdrl_nop_unwrap_eout;
1536 return s;
1537}
1538
1539/* ---------------------------------------------------------------------------*/
1544/* ---------------------------------------------------------------------------*/
1545hdrl_collapse_imagelist_to_image_t *
1546 hdrl_collapse_imagelist_to_image_median(void)
1547{
1548 hdrl_collapse_imagelist_to_image_t * s = cpl_calloc(1, sizeof(*s));
1549 s->func = &hdrl_collapse_median;
1550 s->create_eout = &hdrl_nop_create_eout_img;
1551 s->move_eout = &hdrl_nop_move_eout;
1552 s->unwrap_eout = &hdrl_nop_unwrap_eout;
1553 s->delete_eout = &hdrl_nop_unwrap_eout;
1554 return s;
1555}
1556
1557static void *
1558hdrl_sigclip_create_eout_img(const cpl_image * img)
1559{
1560 cpl_ensure(img, CPL_ERROR_NULL_INPUT, NULL);
1561 hdrl_sigclip_image_output * eout = cpl_calloc(1, sizeof(*eout));
1562 eout->reject_low = cpl_image_new(cpl_image_get_size_x(img),
1563 cpl_image_get_size_y(img),
1564 cpl_image_get_type(img));
1565 eout->reject_high = cpl_image_new(cpl_image_get_size_x(img),
1566 cpl_image_get_size_y(img),
1567 cpl_image_get_type(img));
1568 /* add masks for thread safety on move */
1569 cpl_image_get_bpm(eout->reject_low);
1570 cpl_image_get_bpm(eout->reject_high);
1571 return eout;
1572}
1573
1574static void
1575hdrl_sigclip_delete_eout_img(void * eout_)
1576{
1577 if (eout_ == NULL) {
1578 return;
1579 }
1580 hdrl_sigclip_image_output * eout = (hdrl_sigclip_image_output*)eout_;
1581 cpl_image_delete(eout->reject_low);
1582 cpl_image_delete(eout->reject_high);
1583 cpl_free(eout);
1584}
1585
1586static cpl_error_code
1587hdrl_sigclip_move_eout_img(void * dst_, void * src_, const cpl_size y)
1588{
1589 hdrl_sigclip_image_output * dst = dst_;
1590 hdrl_sigclip_image_output * src = src_;
1591 cpl_ensure_code(dst, CPL_ERROR_NULL_INPUT);
1592 cpl_ensure_code(src, CPL_ERROR_NULL_INPUT);
1593 cpl_ensure_code(y > 0, CPL_ERROR_ACCESS_OUT_OF_RANGE);
1594 cpl_ensure_code(y <= cpl_image_get_size_y(dst->reject_low),
1595 CPL_ERROR_ACCESS_OUT_OF_RANGE);
1596
1597 cpl_image_copy(dst->reject_low, src->reject_low, 1, y);
1598 cpl_image_copy(dst->reject_high, src->reject_high, 1, y);
1599 cpl_image_delete(src->reject_low);
1600 cpl_image_delete(src->reject_high);
1601 cpl_free(src);
1602 return cpl_error_get_code();
1603}
1604
1605
1606/* ---------------------------------------------------------------------------*/
1619/* ---------------------------------------------------------------------------*/
1620hdrl_collapse_imagelist_to_image_t *
1621hdrl_collapse_imagelist_to_image_sigclip(double kappa_low,
1622 double kappa_high,
1623 int niter)
1624{
1625 hdrl_collapse_imagelist_to_image_t * s = cpl_calloc(1, sizeof(*s));
1626 hdrl_parameter * sp =
1627 hdrl_collapse_sigclip_parameter_create(kappa_low, kappa_high, niter);
1628 s->func = &hdrl_collapse_sigclip;
1629 s->create_eout = &hdrl_sigclip_create_eout_img;
1630 s->move_eout = &hdrl_sigclip_move_eout_img;
1631 s->unwrap_eout = &cpl_free;
1632 s->delete_eout = &hdrl_sigclip_delete_eout_img;
1633 s->parameters = sp;
1634 return s;
1635}
1636
1637/* ---------------------------------------------------------------------------*/
1649/* ---------------------------------------------------------------------------*/
1650hdrl_collapse_imagelist_to_image_t *
1651hdrl_collapse_imagelist_to_image_minmax(double nlow,
1652 double nhigh)
1653{
1654 hdrl_collapse_imagelist_to_image_t * s = cpl_calloc(1, sizeof(*s));
1655 hdrl_parameter * sp =
1656 hdrl_collapse_minmax_parameter_create(nlow, nhigh);
1657 s->func = &hdrl_collapse_minmax;
1658 s->create_eout = &hdrl_sigclip_create_eout_img;
1659 s->move_eout = &hdrl_sigclip_move_eout_img;
1660 s->unwrap_eout = &cpl_free;
1661 s->delete_eout = &hdrl_sigclip_delete_eout_img;
1662 s->parameters = sp;
1663 return s;
1664}
1665/* ---------------------------------------------------------------------------*/
1680/* ---------------------------------------------------------------------------*/
1681hdrl_collapse_imagelist_to_image_t *
1682hdrl_collapse_imagelist_to_image_mode( double histo_min,
1683 double histo_max,
1684 double bin_size,
1685 hdrl_mode_type method,
1686 cpl_size error_niter)
1687{
1688 hdrl_collapse_imagelist_to_image_t * s = cpl_calloc(1, sizeof(*s));
1689 hdrl_parameter * sp =
1690 hdrl_collapse_mode_parameter_create(histo_min, histo_max, bin_size,
1691 method, error_niter);
1692 s->func = &hdrl_collapse_mode;
1693 s->create_eout = &hdrl_nop_create_eout_img;
1694 s->move_eout = &hdrl_nop_move_eout;
1695 s->unwrap_eout = &hdrl_nop_unwrap_eout;
1696 s->delete_eout = &hdrl_nop_unwrap_eout;
1697 s->parameters = sp;
1698 return s;
1699}
1700
1701/* ---------------------------------------------------------------------------*/
1715/* ---------------------------------------------------------------------------*/
1716cpl_error_code
1717hdrl_collapse_imagelist_to_image_call(hdrl_collapse_imagelist_to_image_t * f,
1718 const cpl_imagelist * data,
1719 const cpl_imagelist * errors,
1720 cpl_image ** out,
1721 cpl_image ** err,
1722 cpl_image ** contrib,
1723 void ** eout)
1724{
1725 cpl_ensure_code(f, CPL_ERROR_NULL_INPUT);
1726 cpl_ensure_code(data, CPL_ERROR_NULL_INPUT);
1727 cpl_ensure_code(errors, CPL_ERROR_NULL_INPUT);
1728 cpl_ensure_code(out, CPL_ERROR_NULL_INPUT);
1729 cpl_ensure_code(err, CPL_ERROR_NULL_INPUT);
1730 cpl_ensure_code(contrib, CPL_ERROR_NULL_INPUT);
1731 cpl_ensure_code(cpl_imagelist_get_size(data) ==
1732 cpl_imagelist_get_size(errors),
1733 CPL_ERROR_INCOMPATIBLE_INPUT);
1734 if (eout) {
1735 *eout = f->create_eout(cpl_imagelist_get_const(data, 0));
1736 }
1737
1738 cpl_imagelist * nerrors = wrap_synced_errlist(data, errors);
1739 if (nerrors == NULL) {
1740 return cpl_error_get_code();
1741 }
1742
1743 cpl_error_code errcode = f->func(data, nerrors, out, err, contrib,
1744 f->parameters, eout ? *eout : NULL);
1745
1746 unwrap_synced_errlist(nerrors);
1747
1748 return errcode;
1749}
1750
1751/* ---------------------------------------------------------------------------*/
1762/* ---------------------------------------------------------------------------*/
1763void *
1764hdrl_collapse_imagelist_to_image_create_eout(
1765 hdrl_collapse_imagelist_to_image_t * f,
1766 const cpl_image * data)
1767{
1768 cpl_ensure(f, CPL_ERROR_NULL_INPUT, NULL);
1769 cpl_ensure(data, CPL_ERROR_NULL_INPUT, NULL);
1770
1771 return f->create_eout(data);
1772}
1773
1774/* ---------------------------------------------------------------------------*/
1783/* ---------------------------------------------------------------------------*/
1784void
1785hdrl_collapse_imagelist_to_image_unwrap_eout(
1786 hdrl_collapse_imagelist_to_image_t * f,
1787 void * eout)
1788{
1789 if (f != NULL) {
1790 f->unwrap_eout(eout);
1791 }
1792}
1793
1794/* ---------------------------------------------------------------------------*/
1803/* ---------------------------------------------------------------------------*/
1804void
1805hdrl_collapse_imagelist_to_image_delete_eout(
1806 hdrl_collapse_imagelist_to_image_t * f,
1807 void * eout)
1808{
1809 if (f != NULL) {
1810 f->delete_eout(eout);
1811 }
1812}
1813
1814/* ---------------------------------------------------------------------------*/
1827/* ---------------------------------------------------------------------------*/
1828cpl_error_code
1829hdrl_collapse_imagelist_to_image_move_eout(
1830 hdrl_collapse_imagelist_to_image_t * f,
1831 void * dst,
1832 void * src,
1833 cpl_size y)
1834{
1835 cpl_ensure_code(f, CPL_ERROR_NULL_INPUT);
1836
1837 return f->move_eout(dst, src, y);
1838}
1839
1840/* ---------------------------------------------------------------------------*/
1846/* ---------------------------------------------------------------------------*/
1847void hdrl_collapse_imagelist_to_image_delete(hdrl_collapse_imagelist_to_image_t * p)
1848{
1849 if (p) {
1850 hdrl_parameter_delete(p->parameters);
1851 }
1852 cpl_free(p);
1853}
1854
1855/* ---------------------------------------------------------------------------*/
1866/* ---------------------------------------------------------------------------*/
1867void * hdrl_collapse_imagelist_to_vector_create_eout(
1868 hdrl_collapse_imagelist_to_vector_t * f,
1869 const cpl_size size)
1870{
1871 cpl_ensure(f, CPL_ERROR_NULL_INPUT, NULL);
1872 cpl_ensure(size > 0, CPL_ERROR_ILLEGAL_INPUT, NULL);
1873
1874 return f->create_eout(size);
1875}
1876
1877/* ---------------------------------------------------------------------------*/
1886/* ---------------------------------------------------------------------------*/
1887void
1888hdrl_collapse_imagelist_to_vector_unwrap_eout(
1889 hdrl_collapse_imagelist_to_vector_t * f,
1890 void * eout)
1891{
1892 if (f != NULL) {
1893 f->unwrap_eout(eout);
1894 }
1895}
1896
1897/* ---------------------------------------------------------------------------*/
1906/* ---------------------------------------------------------------------------*/
1907void hdrl_collapse_imagelist_to_vector_delete_eout(
1908 hdrl_collapse_imagelist_to_vector_t * f,
1909 void * eout)
1910{
1911 if (f != NULL) {
1912 f->delete_eout(eout);
1913 }
1914}
1915
1916/* ---------------------------------------------------------------------------*/
1929/* ---------------------------------------------------------------------------*/
1930cpl_error_code hdrl_collapse_imagelist_to_vector_move_eout(
1931 hdrl_collapse_imagelist_to_vector_t * f,
1932 void * dst,
1933 void * src,
1934 cpl_size y)
1935{
1936 cpl_ensure_code(f, CPL_ERROR_NULL_INPUT);
1937
1938 return f->move_eout(dst, src, y);
1939}
1940
1941/* ---------------------------------------------------------------------------*/
1962/* ---------------------------------------------------------------------------*/
1963static cpl_error_code
1964reduce_imagelist_to_vector_mean(const cpl_imagelist * data,
1965 const cpl_imagelist * errors,
1966 cpl_vector ** out, cpl_vector ** err,
1967 cpl_array ** contrib,
1968 void * HDRL_UNUSED(parameters),
1969 void * HDRL_UNUSED(extra_out))
1970{
1971 size_t nz = cpl_imagelist_get_size(data);
1972 *out = cpl_vector_new(nz);
1973 *err = cpl_vector_new(nz);
1974 *contrib = cpl_array_new(nz, CPL_TYPE_INT);
1975
1976 for (size_t i = 0; i < nz; i++) {
1977 const cpl_image * img = cpl_imagelist_get_const(data, i);
1978 const cpl_image * ierr = cpl_imagelist_get_const(errors, i);
1979 size_t naccepted = hdrl_get_image_good_npix(img);
1980
1981 if (naccepted != 0) {
1982 double error = sqrt(cpl_image_get_sqflux(ierr)) / naccepted;
1983
1984 cpl_vector_set(*out, i, cpl_image_get_mean(img));
1985 cpl_vector_set(*err, i, error);
1986 }
1987 else {
1988 cpl_vector_set(*out, i, NAN);
1989 cpl_vector_set(*err, i, NAN);
1990 }
1991 cpl_array_set_int(*contrib, i, naccepted);
1992 }
1993
1994 return cpl_error_get_code();
1995}
1996
1997/* ---------------------------------------------------------------------------*/
2002/* ---------------------------------------------------------------------------*/
2003hdrl_collapse_imagelist_to_vector_t *
2004 hdrl_collapse_imagelist_to_vector_mean(void)
2005{
2006 hdrl_collapse_imagelist_to_vector_t * s = cpl_calloc(1, sizeof(*s));
2007 s->create_eout = &hdrl_nop_create_eout_vec;
2008 s->move_eout = &hdrl_nop_move_eout;
2009 s->unwrap_eout = &hdrl_nop_unwrap_eout;
2010 s->delete_eout = &hdrl_nop_unwrap_eout;
2011 s->func = &reduce_imagelist_to_vector_mean;
2012 return s;
2013}
2014
2015/* ---------------------------------------------------------------------------*/
2040/* ---------------------------------------------------------------------------*/
2041static cpl_error_code
2042reduce_imagelist_to_vector_weighted_mean(const cpl_imagelist * data,
2043 const cpl_imagelist * errors,
2044 cpl_vector ** out, cpl_vector ** err,
2045 cpl_array ** contrib,
2046 void * HDRL_UNUSED(parameters),
2047 void * HDRL_UNUSED(extra_out))
2048{
2049 size_t nz = cpl_imagelist_get_size(data);
2050 *out = cpl_vector_new(nz);
2051 *err = cpl_vector_new(nz);
2052 *contrib = cpl_array_new(nz, CPL_TYPE_INT);
2053
2054 for (size_t i = 0; i < nz; i++) {
2055 cpl_image * img =
2056 cpl_image_duplicate(cpl_imagelist_get_const(data, i));
2057 cpl_image * ierr =
2058 cpl_image_duplicate(cpl_imagelist_get_const(errors, i));
2059 size_t naccepted = hdrl_get_image_good_npix(img);
2060
2061 if (naccepted != 0) {
2062 /* (\Sum_i^n w_i * x_i) / (\Sum_i^n w_i) */
2063 /* 1 / \sqrt(\Sum_i^n w_i^2) */
2064 cpl_image_power(ierr, -2);
2065 /* ierr = weights now */
2066 cpl_image_multiply(img, ierr);
2067 double sum_v = cpl_image_get_mean(img) * naccepted;
2068 double sum_w = cpl_image_get_mean(ierr) * naccepted;
2069 double wmean = sum_v / sum_w;
2070 double error = 1. / sqrt(sum_w);
2071
2072 cpl_vector_set(*out, i, wmean);
2073 cpl_vector_set(*err, i, error);
2074 }
2075 else {
2076 cpl_vector_set(*out, i, NAN);
2077 cpl_vector_set(*err, i, NAN);
2078 }
2079 cpl_array_set_int(*contrib, i, naccepted);
2080 cpl_image_delete(img);
2081 cpl_image_delete(ierr);
2082 }
2083
2084 return cpl_error_get_code();
2085}
2086
2087/* ---------------------------------------------------------------------------*/
2092/* ---------------------------------------------------------------------------*/
2093hdrl_collapse_imagelist_to_vector_t *
2094 hdrl_collapse_imagelist_to_vector_weighted_mean(void)
2095{
2096 hdrl_collapse_imagelist_to_vector_t * s = cpl_calloc(1, sizeof(*s));
2097 s->create_eout = &hdrl_nop_create_eout_vec;
2098 s->move_eout = &hdrl_nop_move_eout;
2099 s->unwrap_eout = &hdrl_nop_unwrap_eout;
2100 s->delete_eout = &hdrl_nop_unwrap_eout;
2101 s->func = &reduce_imagelist_to_vector_weighted_mean;
2102 return s;
2103}
2104
2105/* ---------------------------------------------------------------------------*/
2127/* ---------------------------------------------------------------------------*/
2128static cpl_error_code
2129reduce_imagelist_to_vector_median(const cpl_imagelist * data,
2130 const cpl_imagelist * errors,
2131 cpl_vector ** out, cpl_vector ** err,
2132 cpl_array ** contrib,
2133 void * HDRL_UNUSED(parameters),
2134 void * HDRL_UNUSED(extra_out))
2135{
2136 size_t nz = cpl_imagelist_get_size(data);
2137 *out = cpl_vector_new(nz);
2138 *err = cpl_vector_new(nz);
2139 *contrib = cpl_array_new(nz, CPL_TYPE_INT);
2140
2141 for (size_t i = 0; i < nz; i++) {
2142 const cpl_image * img = cpl_imagelist_get_const(data, i);
2143 const cpl_image * ierr = cpl_imagelist_get_const(errors, i);
2144 size_t naccepted = hdrl_get_image_good_npix(img);
2145
2146 if (naccepted != 0) {
2147 double error = sqrt(cpl_image_get_sqflux(ierr)) / naccepted;
2148 /* sqrt(statistical efficiency on normal data)*/
2149 if (naccepted > 2) {
2150 error *= sqrt(CPL_MATH_PI_2);
2151 }
2152
2153 cpl_vector_set(*out, i, cpl_image_get_median(img));
2154 cpl_vector_set(*err, i, error);
2155 }
2156 else {
2157 cpl_vector_set(*out, i, NAN);
2158 cpl_vector_set(*err, i, NAN);
2159 }
2160 cpl_array_set_int(*contrib, i, naccepted);
2161 }
2162
2163 return cpl_error_get_code();
2164}
2165
2166/* ---------------------------------------------------------------------------*/
2171/* ---------------------------------------------------------------------------*/
2172hdrl_collapse_imagelist_to_vector_t *
2173 hdrl_collapse_imagelist_to_vector_median(void)
2174{
2175 hdrl_collapse_imagelist_to_vector_t * s = cpl_calloc(1, sizeof(*s));
2176 s->create_eout = &hdrl_nop_create_eout_vec;
2177 s->move_eout = &hdrl_nop_move_eout;
2178 s->unwrap_eout = &hdrl_nop_unwrap_eout;
2179 s->delete_eout = &hdrl_nop_unwrap_eout;
2180 s->func = &reduce_imagelist_to_vector_median;
2181 return s;
2182}
2183
2184static void *
2185hdrl_sigclip_create_eout_vec(cpl_size size)
2186{
2187 hdrl_sigclip_vector_output * eout = cpl_calloc(1, sizeof(*eout));
2188 eout->reject_low = cpl_vector_new(size);
2189 eout->reject_high = cpl_vector_new(size);
2190 return eout;
2191}
2192
2193static void
2194hdrl_sigclip_delete_eout_vec(void * eout_)
2195{
2196 if (eout_ == NULL) {
2197 return;
2198 }
2199 hdrl_sigclip_vector_output * eout = (hdrl_sigclip_vector_output*)eout_;
2200 cpl_vector_delete(eout->reject_low);
2201 cpl_vector_delete(eout->reject_high);
2202 cpl_free(eout);
2203}
2204
2205static cpl_error_code
2206hdrl_sigclip_move_eout_vec(void * dst_, void * src_, const cpl_size y)
2207{
2208 hdrl_sigclip_vector_output * dst = dst_;
2209 hdrl_sigclip_vector_output * src = src_;
2210 cpl_ensure_code(dst, CPL_ERROR_NULL_INPUT);
2211 cpl_ensure_code(src, CPL_ERROR_NULL_INPUT);
2212 cpl_ensure_code(y >= 0, CPL_ERROR_ACCESS_OUT_OF_RANGE);
2213 cpl_ensure_code(y < cpl_vector_get_size(dst->reject_low),
2214 CPL_ERROR_ACCESS_OUT_OF_RANGE);
2215
2216 double * ddst = cpl_vector_get_data(dst->reject_low);
2217 double * dsrc = cpl_vector_get_data(src->reject_low);
2218 memcpy(ddst + y, dsrc, cpl_vector_get_size(src->reject_low));
2219 ddst = cpl_vector_get_data(dst->reject_high);
2220 dsrc = cpl_vector_get_data(src->reject_high);
2221 memcpy(ddst + y, dsrc, cpl_vector_get_size(src->reject_high));
2222 cpl_vector_delete(src->reject_low);
2223 cpl_vector_delete(src->reject_high);
2224 cpl_free(src);
2225 return cpl_error_get_code();
2226}
2227
2228/* ---------------------------------------------------------------------------*/
2245/* ---------------------------------------------------------------------------*/
2246static cpl_error_code
2247reduce_imagelist_to_vector_sigclip(const cpl_imagelist * data,
2248 const cpl_imagelist * errors,
2249 cpl_vector ** out, cpl_vector ** err,
2250 cpl_array ** contrib, void * parameters,
2251 void * extra_out)
2252{
2253 hdrl_collapse_sigclip_parameter * par = parameters;
2254 hdrl_minmax_vector_output * eout =
2255 (hdrl_minmax_vector_output *)extra_out;
2256
2257 cpl_size nz = cpl_imagelist_get_size(data);
2258 *out = cpl_vector_new(nz);
2259 *err = cpl_vector_new(nz);
2260 *contrib = cpl_array_new(nz, CPL_TYPE_INT);
2261
2262 /* sigmaclip on each image of the imagelist */
2263 for (cpl_size z = 0; z < nz ; z++) {
2264 double corr, error, low, high;
2265 cpl_size contribution;
2266 if (hdrl_kappa_sigma_clip_image(cpl_imagelist_get_const(data, z),
2267 cpl_imagelist_get_const(errors, z),
2268 par->kappa_low,
2269 par->kappa_high,
2270 par->niter,
2271 &corr,
2272 &error,
2273 &contribution,
2274 &low,
2275 &high) != CPL_ERROR_NONE) {
2276 break;
2277 }
2278 cpl_vector_set(*out, z, corr);
2279 cpl_vector_set(*err, z, error);
2280 cpl_array_set_int(*contrib, z, contribution);
2281
2282 if (eout) {
2283 cpl_vector_set(eout->reject_low, z, low);
2284 cpl_vector_set(eout->reject_high, z, high);
2285 }
2286 }
2287
2288 return cpl_error_get_code();
2289}
2290
2291/* ---------------------------------------------------------------------------*/
2305/* ---------------------------------------------------------------------------*/
2306hdrl_collapse_imagelist_to_vector_t *
2307hdrl_collapse_imagelist_to_vector_sigclip(double kappa_low, double kappa_high,
2308 int niter)
2309{
2310 hdrl_collapse_imagelist_to_vector_t * s = cpl_calloc(1, sizeof(*s));
2311 hdrl_parameter * sp =
2312 hdrl_collapse_sigclip_parameter_create(kappa_low, kappa_high, niter);
2313 s->func = &reduce_imagelist_to_vector_sigclip;
2314 s->create_eout = &hdrl_sigclip_create_eout_vec;
2315 s->move_eout = &hdrl_sigclip_move_eout_vec;
2316 s->unwrap_eout = &cpl_free;
2317 s->delete_eout = &hdrl_sigclip_delete_eout_vec;
2318 s->parameters = sp;
2319 return s;
2320}
2321
2322/* ---------------------------------------------------------------------------*/
2338/* ---------------------------------------------------------------------------*/
2339static cpl_error_code
2340reduce_imagelist_to_vector_minmax(const cpl_imagelist * data,
2341 const cpl_imagelist * errors,
2342 cpl_vector ** out, cpl_vector ** err,
2343 cpl_array ** contrib, void * parameters,
2344 void * extra_out)
2345{
2346 hdrl_collapse_minmax_parameter * par = parameters;
2347 hdrl_sigclip_vector_output * eout =
2348 (hdrl_sigclip_vector_output *)extra_out;
2349 cpl_size nz = cpl_imagelist_get_size(data);
2350 *out = cpl_vector_new(nz);
2351 *err = cpl_vector_new(nz);
2352 *contrib = cpl_array_new(nz, CPL_TYPE_INT);
2353
2354 /* minmax on each image of the imagelist */
2355 for (cpl_size z = 0; z < nz ; z++) {
2356 double corr, error, low, high;
2357 cpl_size contribution;
2358 if (hdrl_minmax_clip_image(cpl_imagelist_get_const(data, z),
2359 cpl_imagelist_get_const(errors, z),
2360 par->nlow,
2361 par->nhigh,
2362 &corr,
2363 &error,
2364 &contribution,
2365 &low, &high) != CPL_ERROR_NONE) {
2366 break;
2367 }
2368 cpl_vector_set(*out, z, corr);
2369 cpl_vector_set(*err, z, error);
2370 cpl_array_set_int(*contrib, z, contribution);
2371
2372 if (eout) {
2373 cpl_vector_set(eout->reject_low, z, low);
2374 cpl_vector_set(eout->reject_high, z, high);
2375 }
2376 }
2377
2378 return cpl_error_get_code();
2379}
2380
2381/* ---------------------------------------------------------------------------*/
2392/* ---------------------------------------------------------------------------*/
2393hdrl_collapse_imagelist_to_vector_t *
2394hdrl_collapse_imagelist_to_vector_minmax(double nlow, double nhigh)
2395{
2396 hdrl_collapse_imagelist_to_vector_t * s = cpl_calloc(1, sizeof(*s));
2397 hdrl_parameter * sp =
2398 hdrl_collapse_minmax_parameter_create(nlow, nhigh);
2399 s->func = &reduce_imagelist_to_vector_minmax;
2400 s->create_eout = &hdrl_sigclip_create_eout_vec;
2401 s->move_eout = &hdrl_sigclip_move_eout_vec;
2402 s->unwrap_eout = &cpl_free;
2403 s->delete_eout = &hdrl_sigclip_delete_eout_vec;
2404 s->parameters = sp;
2405 return s;
2406}
2407
2408/* ---------------------------------------------------------------------------*/
2428/* ---------------------------------------------------------------------------*/
2429static cpl_error_code
2430reduce_imagelist_to_vector_mode(const cpl_imagelist * data,
2431 const cpl_imagelist * HDRL_UNUSED(errors),
2432 cpl_vector ** out, cpl_vector ** err,
2433 cpl_array ** contrib, void * parameters,
2434 void * HDRL_UNUSED(extra_out))
2435{
2436
2437 hdrl_collapse_mode_parameter * par = parameters;
2438 cpl_size nz = cpl_imagelist_get_size(data);
2439 *out = cpl_vector_new(nz);
2440 *err = cpl_vector_new(nz);
2441 *contrib = cpl_array_new(nz, CPL_TYPE_INT);
2442
2443 /* mode on each image of the imagelist */
2444 for (cpl_size z = 0; z < nz ; z++) {
2445 double corr, error;
2446 cpl_size contribution;
2447 if (hdrl_mode_clip_image(cpl_imagelist_get_const(data, z),
2448 par->histo_min,
2449 par->histo_max,
2450 par->bin_size,
2451 par->method,
2452 par->error_niter,
2453 &corr,
2454 &error,
2455 &contribution) != CPL_ERROR_NONE) {
2456 break;
2457 }
2458 cpl_vector_set(*out, z, corr);
2459 cpl_vector_set(*err, z, error);
2460 cpl_array_set_int(*contrib, z, contribution);
2461
2462 }
2463
2464 return cpl_error_get_code();
2465}
2466
2467/* ---------------------------------------------------------------------------*/
2472/* ---------------------------------------------------------------------------*/
2473hdrl_collapse_imagelist_to_vector_t *
2474 hdrl_collapse_imagelist_to_vector_mode(double histo_min, double histo_max,
2475 double bin_size, hdrl_mode_type method,
2476 cpl_size error_niter)
2477{
2478
2479 hdrl_collapse_imagelist_to_vector_t * s = cpl_calloc(1, sizeof(*s));
2480 hdrl_parameter * sp =
2481 hdrl_collapse_mode_parameter_create(histo_min, histo_max, bin_size, method, error_niter);
2482
2483 s->func = &reduce_imagelist_to_vector_mode;
2484 s->create_eout = &hdrl_nop_create_eout_vec;
2485 s->move_eout = &hdrl_nop_move_eout;
2486 s->unwrap_eout = &hdrl_nop_unwrap_eout;
2487 s->delete_eout = &hdrl_nop_unwrap_eout;
2488 s->parameters = sp;
2489 return s;
2490}
2491
2492/* ---------------------------------------------------------------------------*/
2506/* ---------------------------------------------------------------------------*/
2507cpl_error_code
2508hdrl_collapse_imagelist_to_vector_call(hdrl_collapse_imagelist_to_vector_t * f,
2509 const cpl_imagelist * data,
2510 const cpl_imagelist * errors,
2511 cpl_vector ** out,
2512 cpl_vector ** err,
2513 cpl_array ** contrib,
2514 void ** eout)
2515{
2516 cpl_ensure_code(f, CPL_ERROR_NULL_INPUT);
2517 cpl_ensure_code(data, CPL_ERROR_NULL_INPUT);
2518 cpl_ensure_code(errors, CPL_ERROR_NULL_INPUT);
2519 cpl_ensure_code(out, CPL_ERROR_NULL_INPUT);
2520 cpl_ensure_code(err, CPL_ERROR_NULL_INPUT);
2521 cpl_ensure_code(contrib, CPL_ERROR_NULL_INPUT);
2522 if (eout) {
2523 *eout = f->create_eout(cpl_imagelist_get_size(data));
2524 }
2525
2526 cpl_imagelist * nerrors = wrap_synced_errlist(data, errors);
2527 if (nerrors == NULL) {
2528 return cpl_error_get_code();
2529 }
2530
2531 cpl_error_code errcode = f->func(data, errors, out, err, contrib,
2532 f->parameters, eout ? *eout : NULL);
2533
2534 unwrap_synced_errlist(nerrors);
2535
2536 return errcode;
2537}
2538
2539/* ---------------------------------------------------------------------------*/
2545/* ---------------------------------------------------------------------------*/
2546void
2547hdrl_collapse_imagelist_to_vector_delete(hdrl_collapse_imagelist_to_vector_t * p)
2548{
2549 if (p) {
2550 cpl_free(p->parameters);
2551 }
2552 cpl_free(p);
2553}
2554
2556
hdrl_collapse_mode_parameter_get_bin_size
double hdrl_collapse_mode_parameter_get_bin_size(const hdrl_parameter *p)
get size of the histogram bins
Definition hdrl_collapse.c:690
hdrl_collapse_mode_parameter_get_histo_min
double hdrl_collapse_mode_parameter_get_histo_min(const hdrl_parameter *p)
get min value
Definition hdrl_collapse.c:658
hdrl_collapse_parameter_is_weighted_mean
cpl_boolean hdrl_collapse_parameter_is_weighted_mean(const hdrl_parameter *self)
check if parameter is a weighted mean parameter
Definition hdrl_collapse.c:316
hdrl_collapse_sigclip_parameter_get_kappa_low
double hdrl_collapse_sigclip_parameter_get_kappa_low(const hdrl_parameter *p)
get low kappa
Definition hdrl_collapse.c:419
hdrl_collapse_mean_parameter_create
hdrl_parameter * hdrl_collapse_mean_parameter_create(void)
create a parameter object for mean
Definition hdrl_collapse.c:255
hdrl_collapse_sigclip_parameter_create
hdrl_parameter * hdrl_collapse_sigclip_parameter_create(double kappa_low, double kappa_high, int niter)
create a parameter object for sigclipped mean
Definition hdrl_collapse.c:333
hdrl_collapse_weighted_mean_parameter_create
hdrl_parameter * hdrl_collapse_weighted_mean_parameter_create(void)
create a parameter object for weighted mean
Definition hdrl_collapse.c:303
hdrl_collapse_parameter_is_mean
cpl_boolean hdrl_collapse_parameter_is_mean(const hdrl_parameter *self)
check if parameter is a mean parameter
Definition hdrl_collapse.c:268
hdrl_collapse_parameter_is_median
cpl_boolean hdrl_collapse_parameter_is_median(const hdrl_parameter *self)
check if parameter is a median parameter
Definition hdrl_collapse.c:292
hdrl_collapse_mode_parameter_get_histo_max
double hdrl_collapse_mode_parameter_get_histo_max(const hdrl_parameter *p)
get high value
Definition hdrl_collapse.c:674
hdrl_collapse_sigclip_parameter_get_niter
int hdrl_collapse_sigclip_parameter_get_niter(const hdrl_parameter *p)
get maximum number of clipping iterations
Definition hdrl_collapse.c:435
hdrl_collapse_mode_parameter_get_error_niter
cpl_size hdrl_collapse_mode_parameter_get_error_niter(const hdrl_parameter *p)
get the error type of the mode
Definition hdrl_collapse.c:722
hdrl_collapse_parameter_is_minmax
cpl_boolean hdrl_collapse_parameter_is_minmax(const hdrl_parameter *self)
check if parameter is a minmax mean parameter
Definition hdrl_collapse.c:475
hdrl_collapse_parameter_is_mode
cpl_boolean hdrl_collapse_parameter_is_mode(const hdrl_parameter *self)
check if parameter is a mode parameter
Definition hdrl_collapse.c:594
hdrl_collapse_parameter_is_sigclip
cpl_boolean hdrl_collapse_parameter_is_sigclip(const hdrl_parameter *self)
check if parameter is a sigclip mean parameter
Definition hdrl_collapse.c:355
hdrl_collapse_sigclip_parameter_get_kappa_high
double hdrl_collapse_sigclip_parameter_get_kappa_high(const hdrl_parameter *p)
get high kappa
Definition hdrl_collapse.c:403
hdrl_collapse_minmax_parameter_get_nlow
double hdrl_collapse_minmax_parameter_get_nlow(const hdrl_parameter *p)
get low value
Definition hdrl_collapse.c:545
hdrl_collapse_median_parameter_create
hdrl_parameter * hdrl_collapse_median_parameter_create(void)
create a parameter object for median
Definition hdrl_collapse.c:279
hdrl_collapse_mode_parameter_create
hdrl_parameter * hdrl_collapse_mode_parameter_create(double histo_min, double histo_max, double bin_size, hdrl_mode_type mode_method, cpl_size error_niter)
create a parameter object for the mode
Definition hdrl_collapse.c:567
hdrl_collapse_parameter_create_parlist
cpl_parameterlist * hdrl_collapse_parameter_create_parlist(const char *base_context, const char *prefix, const char *method_def, hdrl_parameter *sigclip_def, hdrl_parameter *minmax_def, hdrl_parameter *mode_def)
Create parameters for the collapse.
Definition hdrl_collapse.c:748
hdrl_collapse_mode_parameter_get_method
hdrl_mode_type hdrl_collapse_mode_parameter_get_method(const hdrl_parameter *p)
get the mode determination method
Definition hdrl_collapse.c:706
hdrl_collapse_minmax_parameter_create
hdrl_parameter * hdrl_collapse_minmax_parameter_create(double nlow, double nhigh)
create a parameter object for min-max rejected mean
Definition hdrl_collapse.c:454
hdrl_collapse_parameter_parse_parlist
hdrl_parameter * hdrl_collapse_parameter_parse_parlist(const cpl_parameterlist *parlist, const char *prefix)
parse parameterlist for imagelist reduction method
Definition hdrl_collapse.c:845
hdrl_collapse_minmax_parameter_get_nhigh
double hdrl_collapse_minmax_parameter_get_nhigh(const hdrl_parameter *p)
get high value
Definition hdrl_collapse.c:529
hdrl_parameter_delete
void hdrl_parameter_delete(hdrl_parameter *obj)
shallow delete of a parameter
Definition hdrl_parameter.c:148
hdrl_join_string
char * hdrl_join_string(const char *sep_, int n,...)
join strings together
Definition hdrl_utils.c:810
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2