CR2RE Pipeline Reference Manual 1.7.0
irplib_flat.c
1/* $Id: irplib_flat.c,v 1.15 2007-08-07 12:15:41 llundin 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: llundin $
23 * $Date: 2007-08-07 12:15:41 $
24 * $Revision: 1.15 $
25 * $Name: not supported by cvs2svn $
26 */
27
28#ifdef HAVE_CONFIG_H
29#include <config.h>
30#endif
31
32/*-----------------------------------------------------------------------------
33 Includes
34 -----------------------------------------------------------------------------*/
35
36#include "irplib_flat.h"
37#include <math.h>
38
39/*-----------------------------------------------------------------------------
40 Functions prototypes
41 -----------------------------------------------------------------------------*/
42
43static double *irplib_flat_fit_proportional(double *, double *, int);
44
45/*----------------------------------------------------------------------------*/
49/*----------------------------------------------------------------------------*/
50
52
53/*----------------------------------------------------------------------------*/
84/*----------------------------------------------------------------------------*/
85cpl_imagelist *
86irplib_flat_fit_set(cpl_imagelist *raw, int mode)
87{
88 double *plane_med = NULL;
89 double *slope = NULL;
90 cpl_image *gain = NULL;
91 double *pgain = NULL;
92 cpl_image *intercept = NULL;
93 double *pintercept = NULL;
94 cpl_image *sq_err = NULL;
95 double *psq_err = NULL;
96 double *timeline = NULL;
97 float *raw_im_data = NULL;
98 cpl_imagelist *result = NULL;
99 const int nx = cpl_image_get_size_x(cpl_imagelist_get(raw, 0));
100 const int ny = cpl_image_get_size_y(cpl_imagelist_get(raw, 0));
101 const int ni = cpl_imagelist_get_size(raw);
102 int i, j;
103
104 /* Check entries */
105 if (raw == NULL)
106 return NULL;
107 if ((mode != 0) && (mode != 1))
108 return NULL;
109 if (cpl_image_get_type(cpl_imagelist_get(raw, 0)) != CPL_TYPE_FLOAT)
110 return NULL;
111 if (cpl_imagelist_get_size(raw) <= 1)
112 return NULL;
113
114 /* Compute median for all planes */
115 plane_med = cpl_malloc(ni * sizeof(double));
116 for (i = 0; i < ni; i++)
117 plane_med[i] = cpl_image_get_median(cpl_imagelist_get(raw, i));
118
119 /* Create result images */
120 gain = cpl_image_new(nx, ny, CPL_TYPE_DOUBLE);
121 pgain = cpl_image_get_data_double(gain);
122 if (mode == 1) {
123 intercept = cpl_image_new(nx, ny, CPL_TYPE_DOUBLE);
124 pintercept = cpl_image_get_data_double(intercept);
125 }
126 sq_err = cpl_image_new(nx, ny, CPL_TYPE_DOUBLE);
127 psq_err = cpl_image_get_data_double(sq_err);
128 timeline = cpl_malloc(ni * sizeof(double));
129
130 /* Loop on all pixel positions */
131 cpl_msg_info(cpl_func, "Computing gains for all positions (long)...");
132 for (i = 0; i < nx * ny; i++) {
133 /* extract time line */
134 for (j = 0; j < ni; j++) {
135 raw_im_data = cpl_image_get_data_float(cpl_imagelist_get(raw, j));
136 timeline[j] = (double)raw_im_data[i];
137 }
138 /* Fit slope to this time line */
139 if (mode == 1) {
140 slope = irplib_flat_fit_slope_robust(plane_med, timeline, ni);
141 pintercept[i] = slope[0];
142 pgain[i] = slope[1];
143 psq_err[i] = slope[2];
144 /* Set results in output images */
145 }
146 else {
147 slope = irplib_flat_fit_proportional(plane_med, timeline, ni);
148 /* Set results in output images */
149 pgain[i] = slope[0];
150 psq_err[i] = slope[1];
151 }
152 cpl_free(slope);
153 }
154 cpl_free(plane_med);
155 cpl_free(timeline);
156
157 /* Return */
158 result = cpl_imagelist_new();
159 if (mode == 1) {
160 cpl_imagelist_set(result, gain, 0);
161 cpl_imagelist_set(result, intercept, 1);
162 cpl_imagelist_set(result, sq_err, 2);
163 }
164 else {
165 cpl_imagelist_set(result, gain, 0);
166 cpl_imagelist_set(result, sq_err, 1);
167 }
168 return result;
169}
170
171/* @cond */
172#define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
173#define MAX_ITERATE 30
174/* @endcond */
175/*----------------------------------------------------------------------------*/
192/*----------------------------------------------------------------------------*/
193double *
194irplib_flat_fit_slope_robust(double *x, double *y, int np)
195{
196 double *c;
197 double aa, bb, bcomp, b1, b2, del, abdevt, f1, f2, sigb, d, sum;
198 double sx, sy, sxy, sxx, chisq;
199 cpl_vector *arr;
200 double *parr;
201 double aa_ls, bb_ls;
202 int iter;
203 int i;
204
205 /* Check entries */
206 if (x == NULL || y == NULL)
207 return NULL;
208
209 c = cpl_malloc(3 * sizeof(double));
210
211 sx = sy = sxx = sxy = 0.00;
212 for (i = 0; i < np; i++) {
213 sx += x[i];
214 sy += y[i];
215 sxy += x[i] * y[i];
216 sxx += x[i] * x[i];
217 }
218
219 del = np * sxx - sx * sx;
220 aa_ls = aa = (sxx * sy - sx * sxy) / del;
221 bb_ls = bb = (np * sxy - sx * sy) / del;
222
223 chisq = 0.00;
224 for (i = 0; i < np; i++) {
225 double temp = y[i] - (aa + bb * x[i]);
226 temp *= temp;
227 chisq += temp;
228 }
229
230 arr = cpl_vector_new(np);
231 parr = cpl_vector_get_data(arr);
232 sigb = sqrt(chisq / del);
233 b1 = bb;
234
235 bcomp = b1;
236 sum = 0.00;
237 for (i = 0; i < np; i++) {
238 parr[i] = y[i] - bcomp * x[i];
239 }
240 aa = cpl_vector_get_median(arr); /* arr permuted */
241 abdevt = 0.0;
242 for (i = 0; i < np; i++) {
243 d = y[i] - (bcomp * x[i] + aa);
244 abdevt += fabs(d);
245 if (fabs(y[i]) > 1e-7)
246 d /= fabs(y[i]);
247 if (fabs(d) > 1e-7)
248 sum += (d >= 0.0 ? x[i] : -x[i]);
249 }
250 f1 = sum;
251 b2 = bb + SIGN(3.0 * sigb, f1);
252 bcomp = b2;
253 sum = 0.00;
254 for (i = 0; i < np; i++)
255 parr[i] = y[i] - bcomp * x[i];
256 aa = cpl_vector_get_median(arr); /* arr permuted */
257 abdevt = 0.0;
258 for (i = 0; i < np; i++) {
259 d = y[i] - (bcomp * x[i] + aa);
260 abdevt += fabs(d);
261 if (fabs(y[i]) > 1e-7)
262 d /= fabs(y[i]);
263 if (fabs(d) > 1e-7)
264 sum += (d >= 0.0 ? x[i] : -x[i]);
265 }
266 f2 = sum;
267
268 if (fabs(b2 - b1) < 1e-7) {
269 c[0] = aa;
270 c[1] = bb;
271 c[2] = abdevt / (double)np;
272 cpl_vector_delete(arr);
273 return c;
274 }
275
276 iter = 0;
277 while (f1 * f2 > 0.0) {
278 bb = 2.0 * b2 - b1;
279 b1 = b2;
280 f1 = f2;
281 b2 = bb;
282
283 bcomp = b2;
284 sum = 0.00;
285 for (i = 0; i < np; i++)
286 parr[i] = y[i] - bcomp * x[i];
287 aa = cpl_vector_get_median(arr); /* arr permuted */
288 abdevt = 0.0;
289 for (i = 0; i < np; i++) {
290 d = y[i] - (bcomp * x[i] + aa);
291 abdevt += fabs(d);
292 if (fabs(y[i]) > 1e-7)
293 d /= fabs(y[i]);
294 if (fabs(d) > 1e-7)
295 sum += (d >= 0.0 ? x[i] : -x[i]);
296 }
297 f2 = sum;
298 iter++;
299 if (iter >= MAX_ITERATE)
300 break;
301 }
302 if (iter >= MAX_ITERATE) {
303 c[0] = aa_ls;
304 c[1] = bb_ls;
305 c[2] = -1.0;
306 cpl_vector_delete(arr);
307 return c;
308 }
309
310 sigb = 0.01 * sigb;
311 while (fabs(b2 - b1) > sigb) {
312 double f;
313 bb = 0.5 * (b1 + b2);
314 if ((fabs(bb - b1) < 1e-7) || (fabs(bb - b2) < 1e-7))
315 break;
316 bcomp = bb;
317 sum = 0.00;
318 for (i = 0; i < np; i++)
319 parr[i] = y[i] - bcomp * x[i];
320 aa = cpl_vector_get_median(arr); /* arr permuted */
321 abdevt = 0.0;
322 for (i = 0; i < np; i++) {
323 d = y[i] - (bcomp * x[i] + aa);
324 abdevt += fabs(d);
325 if (fabs(y[i]) > 1e-7)
326 d /= fabs(y[i]);
327 if (fabs(d) > 1e-7)
328 sum += (d >= 0.0 ? x[i] : -x[i]);
329 }
330 f = sum;
331
332 if (f * f1 >= 0.0) {
333 f1 = f;
334 b1 = bb;
335 }
336 else {
337 /*f2=f;*/
338 b2 = bb;
339 }
340 }
341 cpl_vector_delete(arr);
342 c[0] = aa;
343 c[1] = bb;
344 c[2] = abdevt / np;
345 return c;
346}
347#undef MAX_ITERATE
348#undef SIGN
349
350
352
353/*----------------------------------------------------------------------------*/
377/*----------------------------------------------------------------------------*/
378#define FITPROP_BIG_SLOPE 1e30
379static double *
380irplib_flat_fit_proportional(double *x, double *y, int np)
381{
382 cpl_vector *slopes;
383 double *pslopes;
384 double *med_slope;
385 double sq_err;
386 int i;
387
388 /* Check entries */
389 if (x == NULL || y == NULL)
390 return NULL;
391
392 slopes = cpl_vector_new(np);
393 pslopes = cpl_vector_get_data(slopes);
394 for (i = 0; i < np; i++) {
395 if (fabs(x[i]) > 1e-30)
396 pslopes[i] = y[i] / x[i];
397 else
398 pslopes[i] = FITPROP_BIG_SLOPE;
399 }
400 med_slope = cpl_malloc(2 * sizeof(double));
401 med_slope[0] = cpl_vector_get_median(slopes); /* slopes permuted */
402 cpl_vector_delete(slopes);
403
404 sq_err = 0.00;
405 for (i = 0; i < np; i++) {
406 double val = med_slope[0] * x[i];
407 sq_err += (val - y[i]) * (val - y[i]);
408 }
409 sq_err /= (double)np;
410 med_slope[1] = sq_err;
411
412 return med_slope;
413#undef FITPROP_BIG_SLOPE
414}
irplib_flat_fit_set
cpl_imagelist * irplib_flat_fit_set(cpl_imagelist *raw, int mode)
Compute a flat-field out of a set of exposures.
Definition irplib_flat.c:86
irplib_flat_fit_slope_robust
double * irplib_flat_fit_slope_robust(double *x, double *y, int np)
Fit a slope to a list of points (robust fit).
Definition irplib_flat.c:194
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2