CR2RE Pipeline Reference Manual 1.7.0
irplib_oddeven.c
1/* $Id: irplib_oddeven.c,v 1.9 2012-01-12 11:50: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: 2012-01-12 11:50:41 $
24 * $Revision: 1.9 $
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 <math.h>
37#include <cpl.h>
38
39#include "irplib_oddeven.h"
40
41/*-----------------------------------------------------------------------------
42 Functions prototypes
43 -----------------------------------------------------------------------------*/
44
45static cpl_imagelist *irplib_oddeven_cube_conv_xy_rtheta(cpl_imagelist *);
46static cpl_imagelist *irplib_oddeven_cube_conv_rtheta_xy(cpl_imagelist *);
47
48/*----------------------------------------------------------------------------*/
52/*----------------------------------------------------------------------------*/
53
55
56/*----------------------------------------------------------------------------*/
64/*----------------------------------------------------------------------------*/
65int
66irplib_oddeven_monitor(const cpl_image *in, int iquad, double *r_even)
67{
68 cpl_image *extracted;
69 cpl_image *labels;
70 int *plabels;
71 int llx, lly, urx, ury;
72 int nx, ny;
73 double f_even, f_tot;
74 cpl_apertures *aperts;
75 int i, j;
76
77 /* Test entries */
78 if (in == NULL || r_even == NULL)
79 return -1;
80 nx = cpl_image_get_size_x(in);
81 ny = cpl_image_get_size_y(in);
82
83 switch (iquad) {
84 case 1:
85 llx = 1;
86 lly = 1;
87 urx = nx / 2;
88 ury = ny / 2;
89 break;
90 case 2:
91 llx = (nx / 2) + 1;
92 lly = 1;
93 urx = nx;
94 ury = ny / 2;
95 break;
96 case 3:
97 llx = 1;
98 lly = (ny / 2) + 1;
99 urx = nx / 2;
100 ury = ny;
101 break;
102 case 4:
103 llx = (nx / 2) + 1;
104 lly = (ny / 2) + 1;
105 urx = nx;
106 ury = ny;
107 break;
108 case 0:
109 llx = 1;
110 lly = 1;
111 urx = nx;
112 ury = ny;
113 break;
114 default:
115 cpl_msg_error(cpl_func, "Unsupported mode");
116 *r_even = 0.0;
117 return -1;
118 }
119
120 /* Extract quadrant */
121 if ((extracted = cpl_image_extract(in, llx, lly, urx, ury)) == NULL) {
122 cpl_msg_error(cpl_func, "Cannot extract quadrant");
123 *r_even = 0.0;
124 return -1;
125 }
126 nx = cpl_image_get_size_x(extracted);
127 ny = cpl_image_get_size_y(extracted);
128
129 /* Get f_tot */
130 f_tot = cpl_image_get_median(extracted);
131 if (fabs(f_tot) < 1e-6) {
132 cpl_msg_warning(cpl_func, "Quadrant median is 0.0");
133 cpl_image_delete(extracted);
134 *r_even = 0.0;
135 return -1;
136 }
137
138 /* Create the label image to define the even columns */
139 labels = cpl_image_new(nx, ny, CPL_TYPE_INT);
140 plabels = cpl_image_get_data_int(labels);
141 for (i = 0; i < nx; i++) {
142 if (i % 2)
143 for (j = 0; j < ny; j++)
144 plabels[i + j * nx] = 0;
145 else
146 for (j = 0; j < ny; j++)
147 plabels[i + j * nx] = 1;
148 }
149
150 /* Get the median of even columns */
151 if ((aperts = cpl_apertures_new_from_image(extracted, labels)) == NULL) {
152 cpl_msg_error(cpl_func, "Cannot compute the even columns median");
153 cpl_image_delete(extracted);
154 cpl_image_delete(labels);
155 *r_even = 0.0;
156 return -1;
157 }
158 cpl_image_delete(extracted);
159 cpl_image_delete(labels);
160 f_even = cpl_apertures_get_median(aperts, 1);
161 cpl_apertures_delete(aperts);
162
163 /* Compute the even rate and return */
164 *r_even = f_even / f_tot;
165 return 0;
166}
167
168/*----------------------------------------------------------------------------*/
174/*----------------------------------------------------------------------------*/
175cpl_image *
176irplib_oddeven_correct(const cpl_image *in)
177{
178 cpl_image *in_real;
179 cpl_image *in_imag;
180 cpl_imagelist *freq_i;
181 cpl_imagelist *freq_i_amp;
182 cpl_image *cur_im;
183 double *pcur_im;
184 cpl_image *cleaned;
185 int nx;
186 cpl_vector *hf_med;
187
188 /* Test entries */
189 if (in == NULL)
190 return NULL;
191
192 nx = cpl_image_get_size_x(in);
193
194 /* Local copy of the input image in DOUBLE */
195 in_real = cpl_image_cast(in, CPL_TYPE_DOUBLE);
196 in_imag = cpl_image_duplicate(in_real);
197 cpl_image_multiply_scalar(in_imag, 0.0);
198
199 /* Apply FFT to input image */
200 cpl_image_fft(in_real, in_imag, CPL_FFT_DEFAULT);
201
202 /* Put the result in an image list */
203 freq_i = cpl_imagelist_new();
204 cpl_imagelist_set(freq_i, in_real, 0);
205 cpl_imagelist_set(freq_i, in_imag, 1);
206
207 /* Convert to amplitude/phase */
208 freq_i_amp = irplib_oddeven_cube_conv_xy_rtheta(freq_i);
209 cpl_imagelist_delete(freq_i);
210
211 /* Correct the odd-even frequency */
212 cur_im = cpl_imagelist_get(freq_i_amp, 0);
213 pcur_im = cpl_image_get_data_double(cur_im);
214 /* Odd-even frequency will be replaced by
215 the median of the 5 values around */
216 hf_med = cpl_vector_new(5);
217
218 cpl_vector_set(hf_med, 0, pcur_im[nx / 2 + 1]);
219 cpl_vector_set(hf_med, 1, pcur_im[nx / 2 + 2]);
220 cpl_vector_set(hf_med, 2, pcur_im[nx / 2 + 3]);
221 cpl_vector_set(hf_med, 3, pcur_im[nx / 2]);
222 cpl_vector_set(hf_med, 4, pcur_im[nx / 2 - 1]);
223
224 pcur_im[nx / 2 + 1] = cpl_vector_get_median(hf_med);
225 cpl_vector_delete(hf_med);
226
227 /* Convert to X/Y */
228 freq_i = irplib_oddeven_cube_conv_rtheta_xy(freq_i_amp);
229 cpl_imagelist_delete(freq_i_amp);
230 /* FFT back to image space */
231 cpl_image_fft(cpl_imagelist_get(freq_i, 0), cpl_imagelist_get(freq_i, 1),
232 CPL_FFT_INVERSE);
233 cleaned = cpl_image_cast(cpl_imagelist_get(freq_i, 0), CPL_TYPE_FLOAT);
234 cpl_imagelist_delete(freq_i);
235 return cleaned;
236}
237
239
240/*----------------------------------------------------------------------------*/
251/*----------------------------------------------------------------------------*/
252static cpl_imagelist *
253irplib_oddeven_cube_conv_xy_rtheta(cpl_imagelist *cube_in)
254{
255 cpl_imagelist *cube_out;
256 double re, im;
257 double mod, phase;
258 int nx, ny, np;
259 cpl_image *tmp_im;
260 double *pim1;
261 double *pim2;
262 double *pim3;
263 double *pim4;
264 int i, j;
265
266 /* Error handling : test entries */
267 if (cube_in == NULL)
268 return NULL;
269 np = cpl_imagelist_get_size(cube_in);
270 if (np != 2)
271 return NULL;
272
273 /* Initialise */
274 tmp_im = cpl_imagelist_get(cube_in, 0);
275 pim1 = cpl_image_get_data_double(tmp_im);
276 nx = cpl_image_get_size_x(tmp_im);
277 ny = cpl_image_get_size_y(tmp_im);
278 tmp_im = cpl_imagelist_get(cube_in, 1);
279 pim2 = cpl_image_get_data_double(tmp_im);
280
281 /* Allocate cube_out */
282 cube_out = cpl_imagelist_duplicate(cube_in);
283
284 tmp_im = cpl_imagelist_get(cube_out, 0);
285 pim3 = cpl_image_get_data_double(tmp_im);
286 tmp_im = cpl_imagelist_get(cube_out, 1);
287 pim4 = cpl_image_get_data_double(tmp_im);
288 /* Convert */
289 for (j = 0; j < ny; j++) {
290 for (i = 0; i < nx; i++) {
291 re = (double)pim1[i + j * nx];
292 im = (double)pim2[i + j * nx];
293 mod = (double)(sqrt(re * re + im * im));
294 if (re != 0.0)
295 phase = (double)atan2(im, re);
296 else
297 phase = 0.0;
298 pim3[i + j * nx] = mod;
299 pim4[i + j * nx] = phase;
300 }
301 }
302 return cube_out;
303}
304
305/*----------------------------------------------------------------------------*/
318/*----------------------------------------------------------------------------*/
319static cpl_imagelist *
320irplib_oddeven_cube_conv_rtheta_xy(cpl_imagelist *cube_in)
321{
322 cpl_imagelist *cube_out;
323 double re, im;
324 double mod, phase;
325 int nx, ny, np;
326 cpl_image *tmp_im;
327 double *pim1;
328 double *pim2;
329 double *pim3;
330 double *pim4;
331 int i, j;
332
333 /* Error handling : test entries */
334 if (cube_in == NULL)
335 return NULL;
336 np = cpl_imagelist_get_size(cube_in);
337 if (np != 2)
338 return NULL;
339
340 /* Initialise */
341 tmp_im = cpl_imagelist_get(cube_in, 0);
342 pim1 = cpl_image_get_data_double(tmp_im);
343 nx = cpl_image_get_size_x(tmp_im);
344 ny = cpl_image_get_size_y(tmp_im);
345 tmp_im = cpl_imagelist_get(cube_in, 1);
346 pim2 = cpl_image_get_data_double(tmp_im);
347
348 /* Allocate cube_out */
349 cube_out = cpl_imagelist_duplicate(cube_in);
350
351 tmp_im = cpl_imagelist_get(cube_out, 0);
352 pim3 = cpl_image_get_data_double(tmp_im);
353 tmp_im = cpl_imagelist_get(cube_out, 1);
354 pim4 = cpl_image_get_data_double(tmp_im);
355 /* Convert */
356 for (j = 0; j < ny; j++) {
357 for (i = 0; i < nx; i++) {
358 mod = (double)pim1[i + j * nx];
359 phase = (double)pim2[i + j * nx];
360 re = (double)(mod * cos(phase));
361 im = (double)(mod * sin(phase));
362 pim3[i + j * nx] = re;
363 pim4[i + j * nx] = im;
364 }
365 }
366 return cube_out;
367}
irplib_oddeven_monitor
int irplib_oddeven_monitor(const cpl_image *in, int iquad, double *r_even)
Estimate the odd/even rate in an image quadrant.
Definition irplib_oddeven.c:66
irplib_oddeven_correct
cpl_image * irplib_oddeven_correct(const cpl_image *in)
Correct the odd/even in an image.
Definition irplib_oddeven.c:176
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2