CR2RE Pipeline Reference Manual 1.7.0
irplib_ksigma_clip.c
1/* $Id: irplib_ksigma_clip.c,v 1.1 2011-11-02 13:18:28 amodigli Exp $
2 *
3 * This file is part of the irplib package
4 * Copyright (C) 2002, 2003 European Southern Observatory
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02111-1307 USA
19 */
20
21/*
22 * $Author: amodigli $
23 * $Date: 2011-11-02 13:18:28 $
24 * $Revision: 1.1 $
25 * $Name: not supported by cvs2svn $
26 */
27
28#ifdef HAVE_CONFIG_H
29#include <config.h>
30#endif
31
32#include <complex.h>
33
34/*---------------------------------------------------------------------------
35 Includes
36 ---------------------------------------------------------------------------*/
37
38#include <math.h>
39#include <string.h>
40#include <assert.h>
41#include <float.h>
42
43#include <cpl.h>
44
45#include "irplib_ksigma_clip.h"
46
47#include "irplib_hist.h"
48#include "irplib_utils.h"
49
50/*--------------------------------------------------------------------------*/
51
52/*
53 * @defgroup ksigmaclip kappa sigma clip functions
54 */
55
56/*--------------------------------------------------------------------------*/
57
58/*---------------------------------------------------------------------------
59 Defines
60 ---------------------------------------------------------------------------*/
61
62/*---------------------------------------------------------------------------
63 Private function prototypes
64 ---------------------------------------------------------------------------*/
65
66static cpl_error_code irplib_ksigma_clip_double(const double *pi,
67 int llx,
68 int lly,
69 int urx,
70 int ury,
71 int nx,
72 double var_sum,
73 int npixs,
74 double kappa,
75 int nclip,
76 double tolerance,
77 double *mean,
78 double *stdev);
79
80static cpl_error_code irplib_ksigma_clip_float(const float *pi,
81 int llx,
82 int lly,
83 int urx,
84 int ury,
85 int nx,
86 double var_sum,
87 int npixs,
88 double kappa,
89 int nclip,
90 double tolerance,
91 double *mean,
92 double *stdev);
93
94static cpl_error_code irplib_ksigma_clip_int(const int *pi,
95 int llx,
96 int lly,
97 int urx,
98 int ury,
99 int nx,
100 double var_sum,
101 int npixs,
102 double kappa,
103 int nclip,
104 double tolerance,
105 double *mean,
106 double *stdev);
107
108
109/*---------------------------------------------------------------------------*/
161/*---------------------------------------------------------------------------*/
162cpl_error_code
163irplib_ksigma_clip(const cpl_image *img,
164 int llx,
165 int lly,
166 int urx,
167 int ury,
168 double kappa,
169 int nclip,
170 double tolerance,
171 double *kmean,
172 double *kstdev)
173{
174 cpl_errorstate inistate = cpl_errorstate_get();
175
176 int nx, ny;
177
178 cpl_stats *stats;
179 double mean, stdev, var_sum;
180 int npixs;
181
182 cpl_ensure_code(img != NULL, CPL_ERROR_NULL_INPUT);
183
184 nx = cpl_image_get_size_x(img);
185 ny = cpl_image_get_size_y(img);
186
187 cpl_ensure_code(llx > 0 && urx > llx && urx <= nx && lly > 0 && ury > lly &&
188 ury <= ny,
189 CPL_ERROR_ILLEGAL_INPUT);
190
191 cpl_ensure_code(tolerance >= 0.0, CPL_ERROR_ILLEGAL_INPUT);
192 cpl_ensure_code(kappa > 1.0, CPL_ERROR_ILLEGAL_INPUT);
193 cpl_ensure_code(nclip > 0, CPL_ERROR_ILLEGAL_INPUT);
194
195 stats =
196 cpl_stats_new_from_image_window(img, CPL_STATS_MEAN | CPL_STATS_STDEV,
197 llx, lly, urx, ury);
198
199 npixs = cpl_stats_get_npix(stats); /* Non-bad pixels in window */
200 mean = cpl_stats_get_mean(stats);
201 stdev = cpl_stats_get_stdev(stats);
202 var_sum = stdev * stdev * (npixs - 1);
203
204 cpl_stats_delete(stats);
205
206 /* img, llx etc. may cause errors: Check and propagate */
207 cpl_ensure_code(cpl_errorstate_is_equal(inistate), cpl_error_get_code());
208
209 switch (cpl_image_get_type(img)) {
210 case CPL_TYPE_DOUBLE:
211 skip_if(irplib_ksigma_clip_double(
212 cpl_image_get_data_double_const(img), llx, lly, urx, ury, nx,
213 var_sum, npixs, kappa, nclip, tolerance, &mean, &stdev));
214 break;
215 case CPL_TYPE_FLOAT:
216 skip_if(irplib_ksigma_clip_float(
217 cpl_image_get_data_float_const(img), llx, lly, urx, ury, nx,
218 var_sum, npixs, kappa, nclip, tolerance, &mean, &stdev));
219 break;
220 case CPL_TYPE_INT:
221 skip_if(irplib_ksigma_clip_int(cpl_image_get_data_int_const(img),
222 llx, lly, urx, ury, nx, var_sum,
223 npixs, kappa, nclip, tolerance,
224 &mean, &stdev));
225 break;
226 default:
227 /* It is an error in CPL to reach this point */
228 assert(0);
229 }
230
231 *kmean = mean;
232 if (kstdev != NULL)
233 *kstdev = stdev; /* Optional */
234
235 end_skip;
236
237 return cpl_error_get_code();
238}
239
240#define CONCAT(a, b) a##_##b
241#define CONCAT2X(a, b) CONCAT(a, b)
242
243#define CPL_TYPE double
244#include "irplib_ksigma_clip_body.h"
245#undef CPL_TYPE
246
247#define CPL_TYPE float
248#include "irplib_ksigma_clip_body.h"
249#undef CPL_TYPE
250
251#define CPL_TYPE int
252#include "irplib_ksigma_clip_body.h"
253#undef CPL_TYPE
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2