CR2RE Pipeline Reference Manual 1.7.0
irplib_strehl.c
1/* $Id: irplib_strehl.c,v 1.43 2009-11-18 21:37:48 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: 2009-11-18 21:37:48 $
24 * $Revision: 1.43 $
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_strehl.h"
37#include "irplib_utils.h"
38
39#include <assert.h>
40#include <stdint.h>
41#include <math.h>
42
43/*----------------------------------------------------------------------------*/
47/*----------------------------------------------------------------------------*/
48
49/*-----------------------------------------------------------------------------
50 Define
51 -----------------------------------------------------------------------------*/
52
53#ifndef IRPLIB_STREHL_RAD_CENTRAL
54#define IRPLIB_STREHL_RAD_CENTRAL 5
55#endif
56
57#ifndef IRPLIB_STREHL_DETECT_LEVEL
58#define IRPLIB_STREHL_DETECT_LEVEL 5.0
59#endif
60
61#define IRPLIB_DISK_BG_MIN_PIX_NB 30
62#define IRPLIB_DISK_BG_REJ_LOW 0.1
63#define IRPLIB_DISK_BG_REJ_HIGH 0.1
64
65#ifdef CPL_MIN
66#define IRPLIB_MIN CPL_MIN
67#else
68#define IRPLIB_MIN(A, B) (((A) < (B)) ? (A) : (B))
69#endif
70
71#ifdef CPL_MAX
72#define IRPLIB_MAX CPL_MAX
73#else
74#define IRPLIB_MAX(A, B) (((A) > (B)) ? (A) : (B))
75#endif
76
77/*-----------------------------------------------------------------------------
78 Functions prototypes
79 -----------------------------------------------------------------------------*/
80
81static cpl_image *
82irplib_strehl_generate_otf(double, double, double, double, int, double);
83static double PSF_H1(double, double, double);
84static double PSF_H2(double, double);
85static double PSF_G(double, double);
86static double PSF_sinc_norm(double);
87static double PSF_TelOTF(double, double);
88
89#ifndef IRPLIB_NO_FIT_GAUSSIAN
90#ifdef IRPLIB_STREHL_USE_CPL_IMAGE_FIT_GAUSSIAN
91static double irplib_gaussian_2d(double, double, double, double, double);
92#endif
93
94#if defined CPL_VERSION_CODE && CPL_VERSION_CODE >= CPL_VERSION(6, 9, 1)
95#define irplib_gaussian_eval_2d cpl_gaussian_eval_2d
96#else
97static double irplib_gaussian_eval_2d(const cpl_array *, double, double);
98#endif
99
100static uint32_t irplib_roundup_power2(uint32_t v) CPL_ATTR_CONST;
101
102static cpl_error_code irplib_gaussian_maxpos(const cpl_image *,
103 double,
104 double,
105 double,
106 double *,
107 double *,
108 double *);
109
110static cpl_error_code irplib_closeset_aperture(const cpl_apertures *self,
111 const double x,
112 const double y,
113 int *ind);
114#endif
115
116/*-----------------------------------------------------------------------------
117 Functions code
118 -----------------------------------------------------------------------------*/
120
121/*----------------------------------------------------------------------------*/
151/*----------------------------------------------------------------------------*/
152cpl_error_code
153irplib_strehl_compute(const cpl_image *im,
154 double m1,
155 double m2,
156 double lam,
157 double dlam,
158 double pscale,
159 int size,
160 double xpos,
161 double ypos,
162 double r1,
163 double r2,
164 double r3,
165 int noise_box_sz,
166 int noise_nsamples,
167 double *strehl,
168 double *strehl_err,
169 double *star_bg,
170 double *star_peak,
171 double *star_flux,
172 double *psf_peak,
173 double *psf_flux,
174 double *bg_noise)
175{
176 cpl_image *psf;
177 double star_radius, max_radius;
178
179 /* FIXME: Arbitrary choice of image border */
180 const double window_size = (double)(IRPLIB_STREHL_RAD_CENTRAL);
181
182 /* Determined empirically by C. Lidman for Strehl error computation */
183 const double strehl_error_coefficient = CPL_MATH_PI * 0.007 / 0.0271;
184 double ring[4];
185 /* cpl_flux_get_noise_ring() must succeed with this many tries */
186 int ring_tries = 3;
187#ifndef IRPLIB_NO_FIT_GAUSSIAN
188 double xposfit = 0.0, yposfit = 0.0, peak = 0.0;
189 cpl_error_code code;
190#endif
191 cpl_errorstate prestate = cpl_errorstate_get();
192
193 /* Check compile-time constant */
194 cpl_ensure_code(window_size > 0.0, CPL_ERROR_ILLEGAL_INPUT);
195
196 /* Test inputs */
197 cpl_ensure_code(im != NULL, CPL_ERROR_NULL_INPUT);
198 cpl_ensure_code(strehl != NULL, CPL_ERROR_NULL_INPUT);
199 cpl_ensure_code(strehl_err != NULL, CPL_ERROR_NULL_INPUT);
200 cpl_ensure_code(star_bg != NULL, CPL_ERROR_NULL_INPUT);
201 cpl_ensure_code(star_peak != NULL, CPL_ERROR_NULL_INPUT);
202 cpl_ensure_code(star_flux != NULL, CPL_ERROR_NULL_INPUT);
203 cpl_ensure_code(psf_peak != NULL, CPL_ERROR_NULL_INPUT);
204 cpl_ensure_code(psf_flux != NULL, CPL_ERROR_NULL_INPUT);
205
206 cpl_ensure_code(pscale > 0.0, CPL_ERROR_ILLEGAL_INPUT);
207
208 cpl_ensure_code(r1 > 0.0, CPL_ERROR_ILLEGAL_INPUT);
209 cpl_ensure_code(r2 > 0.0, CPL_ERROR_ILLEGAL_INPUT);
210 cpl_ensure_code(r3 > r2, CPL_ERROR_ILLEGAL_INPUT);
211
212 /* Computing a Strehl ratio is a story between an ideal PSF */
213 /* and a candidate image supposed to approximate this ideal PSF. */
214
215 /* Generate first appropriate PSF to find max peak */
216 psf = irplib_strehl_generate_psf(m1, m2, lam, dlam, pscale, size);
217 if (psf == NULL) {
218 return cpl_error_set_where(cpl_func);
219 }
220
221 /* Compute flux in PSF and find max peak */
222 *psf_peak = cpl_image_get_max(psf);
223 cpl_image_delete(psf);
224
225 assert(*psf_peak > 0.0); /* The ideal PSF has a positive maximum */
226 *psf_flux = 1.0; /* The psf flux, cpl_image_get_flux(psf), is always 1 */
227
228#ifndef IRPLIB_NO_FIT_GAUSSIAN
229 code = irplib_gaussian_maxpos(im, IRPLIB_STREHL_DETECT_LEVEL, xpos, ypos,
230 &xposfit, &yposfit, &peak);
231 if (code) {
232 cpl_errorstate_set(prestate);
233 }
234 else {
235 xpos = xposfit;
236 ypos = yposfit;
237 }
238#endif
239
240 /* Measure the background in the candidate image */
241 *star_bg =
242 irplib_strehl_ring_background(im, xpos, ypos, r2 / pscale, r3 / pscale,
243 IRPLIB_BG_METHOD_AVER_REJ);
244 if (!cpl_errorstate_is_equal(prestate)) {
245 return cpl_error_set_where(cpl_func);
246 }
247
248 /* Compute star_radius in pixels */
249 star_radius = r1 / pscale;
250
251 /* Measure the flux on the candidate image */
252 *star_flux = irplib_strehl_disk_flux(im, xpos, ypos, star_radius, *star_bg);
253
254 if (*star_flux <= 0.0) {
255 return cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT,
256 "Non-positive star flux=%g (Star "
257 "background=%g)",
258 *star_flux, *star_bg);
259 }
260
261 /* Find the peak value on the central part of the candidate image */
262 max_radius = window_size < star_radius ? window_size : star_radius;
263 cpl_ensure_code(!irplib_strehl_disk_max(im, xpos, ypos, max_radius,
264 star_peak),
265 cpl_error_get_code());
266 *star_peak -= *star_bg;
267
268 if (*star_flux <= 0.0) {
269 return cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT,
270 "Non-positive star peak=%g (Star "
271 "background=%g, Star flux=%g)",
272 *star_flux, *star_bg, *star_flux);
273 }
274
275 /* Compute Strehl */
276 /* (StarPeak / StarFlux) / (PsfPeak / PsfFlux) */
277 *strehl = (*star_peak * *psf_flux) / (*star_flux * *psf_peak);
278
279#ifndef IRPLIB_NO_FIT_GAUSSIAN
280 if (code == CPL_ERROR_NONE && peak > *star_peak && *star_peak > 0.0 &&
281 *strehl * peak / *star_peak <= 1.0) {
282 cpl_msg_debug(cpl_func, "Increasing Strehl from %g: %g (%g)", *strehl,
283 *strehl * peak / *star_peak, peak / *star_peak);
284 *strehl *= peak / *star_peak;
285 *star_peak = peak;
286 }
287#endif
288
289 /* Compute Strehl error */
290 ring[0] = xpos;
291 ring[1] = ypos;
292 ring[2] = r2 / pscale;
293 ring[3] = r3 / pscale;
294
295 while (cpl_flux_get_noise_ring(im, ring, noise_box_sz, noise_nsamples,
296 bg_noise, NULL) &&
297 --ring_tries > 0)
298 ;
299 if (ring_tries > 0) {
300 cpl_errorstate_set(prestate); /* Recover, if an error happened */
301 }
302 else {
303 return cpl_error_set_where(cpl_func);
304 }
305
306 *strehl_err = strehl_error_coefficient * (*bg_noise) * pscale *
307 star_radius * star_radius / *star_flux;
308
309 if (*strehl > 1.0) {
310 cpl_msg_warning(cpl_func,
311 "Extreme Strehl-ratio=%g (strehl-error=%g, "
312 "star_peak=%g, star_flux=%g, psf_peak=%g, psf_flux=%g)",
313 *strehl, *strehl_err, *star_peak, *star_flux, *psf_peak,
314 *psf_flux);
315 }
316
317 /* This check should not be able to fail, but just to be sure */
318 return *strehl_err >= 0.0
319 ? CPL_ERROR_NONE
320 : cpl_error_set_message(
321 cpl_func, CPL_ERROR_ILLEGAL_OUTPUT,
322 "Negative strehl-error=%g (Strehl-ratio=%g, "
323 "star_peak=%g, star_flux=%g, psf_peak=%g, "
324 "psf_flux=%g",
325 *strehl_err, *strehl, *star_peak, *star_flux, *psf_peak,
326 *psf_flux);
327}
328
329/*----------------------------------------------------------------------------*/
342/*----------------------------------------------------------------------------*/
343double
344irplib_strehl_disk_flux(const cpl_image *im,
345 double xpos,
346 double ypos,
347 double rad,
348 double bg)
349{
350 const int nx = cpl_image_get_size_x(im);
351 const int ny = cpl_image_get_size_y(im);
352 /* Round down */
353 const int lx = (int)(xpos - rad);
354 const int ly = (int)(ypos - rad);
355 /* Round up */
356 const int ux = (int)(xpos + rad) + 1;
357 const int uy = (int)(ypos + rad) + 1;
358
359 const double sqr = rad * rad;
360 double flux = 0.0;
361 int i, j;
362
363
364 /* Check entries */
365 cpl_ensure(im != NULL, CPL_ERROR_NULL_INPUT, 0.0);
366 cpl_ensure(rad > 0.0, CPL_ERROR_ILLEGAL_INPUT, 0.0);
367
368 for (j = IRPLIB_MAX(ly, 0); j < IRPLIB_MIN(uy, ny); j++) {
369 const double yj = (double)j - ypos;
370 for (i = IRPLIB_MAX(lx, 0); i < IRPLIB_MIN(ux, nx); i++) {
371 const double xi = (double)i - xpos;
372 const double dist = yj * yj + xi * xi;
373 if (dist <= sqr) {
374 int isbad;
375 const double value = cpl_image_get(im, i + 1, j + 1, &isbad);
376
377 if (!isbad) {
378 flux += value - bg;
379 }
380 }
381 }
382 }
383
384 return flux;
385}
386
387/*----------------------------------------------------------------------------*/
399/*----------------------------------------------------------------------------*/
400double
401irplib_strehl_ring_background(const cpl_image *im,
402 double xpos,
403 double ypos,
404 double rad_int,
405 double rad_ext,
406 irplib_strehl_bg_method mode)
407{
408 const int nx = cpl_image_get_size_x(im);
409 const int ny = cpl_image_get_size_y(im);
410 /* Round down */
411 const int lx = (int)(xpos - rad_ext);
412 const int ly = (int)(ypos - rad_ext);
413 /* Round up */
414 const int ux = (int)(xpos + rad_ext) + 1;
415 const int uy = (int)(ypos + rad_ext) + 1;
416 int mpix, npix;
417 const double sqr_int = rad_int * rad_int;
418 const double sqr_ext = rad_ext * rad_ext;
419 cpl_vector *pix_arr;
420 double flux = 0.0;
421 int i, j;
422
423 /* Check entries */
424 cpl_ensure(im != NULL, CPL_ERROR_NULL_INPUT, 0.0);
425 cpl_ensure(rad_int > 0.0, CPL_ERROR_ILLEGAL_INPUT, 0.0);
426 cpl_ensure(rad_ext > rad_int, CPL_ERROR_ILLEGAL_INPUT, 0.0);
427
428 cpl_ensure(mode == IRPLIB_BG_METHOD_AVER_REJ ||
429 mode == IRPLIB_BG_METHOD_MEDIAN,
430 CPL_ERROR_UNSUPPORTED_MODE, 0.0);
431
432 mpix = (int)((2.0 * rad_ext + 1.0) * (2.0 * rad_ext + 1.0));
433
434 /* Allocate pixel array to hold values in the ring */
435 pix_arr = cpl_vector_new(mpix);
436
437 /* Count number of pixels in the ring */
438 /* Retrieve all pixels which belong to the ring */
439 npix = 0;
440 for (j = IRPLIB_MAX(ly, 0); j < IRPLIB_MIN(uy, ny); j++) {
441 const double yj = (double)j - ypos;
442 for (i = IRPLIB_MAX(lx, 0); i < IRPLIB_MIN(ux, nx); i++) {
443 const double xi = (double)i - xpos;
444 const double dist = yj * yj + xi * xi;
445 if (sqr_int <= dist && dist <= sqr_ext) {
446 int isbad;
447 const double value = cpl_image_get(im, i + 1, j + 1, &isbad);
448
449 if (!isbad) {
450 cpl_vector_set(pix_arr, npix, value);
451 npix++;
452 }
453 }
454 }
455 }
456
457 assert(npix <= mpix);
458
459 if (npix < IRPLIB_DISK_BG_MIN_PIX_NB) {
460 cpl_vector_delete(pix_arr);
461 (void)cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND,
462 "Need "
463 "at least %d (not %d <= %d) samples to "
464 "compute noise",
465 IRPLIB_DISK_BG_MIN_PIX_NB, npix, mpix);
466 return 0.0;
467 }
468
469 /* Should not be able to fail now */
470
471 /* Resize pixel array to actual number of values within the ring */
472 pix_arr = cpl_vector_wrap(npix, (double *)cpl_vector_unwrap(pix_arr));
473
474 if (mode == IRPLIB_BG_METHOD_AVER_REJ) {
475 const int low_ind = (int)((double)npix * IRPLIB_DISK_BG_REJ_LOW);
476 const int high_ind =
477 (int)((double)npix * (1.0 - IRPLIB_DISK_BG_REJ_HIGH));
478
479 /* Sort the array */
480 cpl_vector_sort(pix_arr, CPL_SORT_ASCENDING);
481
482 for (i = low_ind; i < high_ind; i++) {
483 flux += cpl_vector_get(pix_arr, i);
484 }
485 if (high_ind - low_ind > 1)
486 flux /= (double)(high_ind - low_ind);
487 }
488 else /* if (mode == IRPLIB_BG_METHOD_MEDIAN) */ {
489 flux = cpl_vector_get_median(pix_arr);
490 }
491
492 cpl_vector_delete(pix_arr);
493
494 return flux;
495}
496
497/*----------------------------------------------------------------------------*/
517/*----------------------------------------------------------------------------*/
518cpl_image *
519irplib_strehl_generate_psf(double m1,
520 double m2,
521 double lam,
522 double dlam,
523 double pscale,
524 int size)
525{
526 cpl_image *otf_image =
527 irplib_strehl_generate_otf(m1, m2, lam, dlam, size, pscale);
528
529 if (otf_image == NULL ||
530
531 /* Transform back to real space
532 - Normalization is unnecessary, due to the subsequent normalisation.
533 - An OTF is point symmetric about its center, i.e. it is even,
534 i.e. the real space image is real.
535 - Because of this a forward FFT works as well.
536 - If the PSF ever needs to have its images halves swapped add
537 CPL_FFT_SWAP_HALVES to the FFT call.
538 */
539
540 cpl_image_fft(otf_image, NULL, CPL_FFT_UNNORMALIZED) ||
541
542 /* Compute absolute values of PSF */
543 cpl_image_abs(otf_image) ||
544
545 /* Normalize PSF to get flux=1 */
546 cpl_image_normalise(otf_image, CPL_NORM_FLUX)) {
547 (void)cpl_error_set_where(cpl_func);
548 cpl_image_delete(otf_image);
549 otf_image = NULL;
550 }
551
552 return otf_image;
553}
554
556
557/*----------------------------------------------------------------------------*/
573/*----------------------------------------------------------------------------*/
574static cpl_image *
575irplib_strehl_generate_otf(double m1,
576 double m2,
577 double lam,
578 double dlam,
579 int size,
580 double pscale)
581{
582 double *otf_data;
583 /* Obscuration ratio, m1 / m2 */
584 const double obs_ratio = m1 != 0.0 ? m2 / m1 : 0.0;
585 /* pixel scale converted from Arsecond to radian */
586 const double rpscale = pscale * CPL_MATH_2PI / (double)(360 * 60 * 60);
587 /* Cut-off frequency in pixels per central wavelength (in m) */
588 const double f_max = m1 * rpscale * (double)size;
589
590 /* Pixel corresponding to the zero frequency */
591 const int pix0 = size / 2;
592 int i, j;
593
594
595 cpl_ensure(m2 > 0.0, CPL_ERROR_ILLEGAL_INPUT, NULL);
596 cpl_ensure(m1 > m2, CPL_ERROR_ILLEGAL_INPUT, NULL);
597 cpl_ensure(dlam > 0.0, CPL_ERROR_ILLEGAL_INPUT, NULL);
598 cpl_ensure(pscale > 0.0, CPL_ERROR_ILLEGAL_INPUT, NULL);
599 cpl_ensure(size > 0, CPL_ERROR_ILLEGAL_INPUT, NULL);
600 /* Due the the FFT, size is actually required to be a power of two */
601 cpl_ensure(size % 2 == 0, CPL_ERROR_ILLEGAL_INPUT, NULL);
602
603 /* Ensure positive lambda */
604 cpl_ensure(2.0 * lam > dlam, CPL_ERROR_ILLEGAL_INPUT, NULL);
605
606 /* Convert wavelengths from micron to meter */
607 lam /= 1.0e6;
608 dlam /= 1.0e6;
609
610 /* Allocate the output pixel buffer */
611 otf_data = (double *)cpl_malloc(size * size * sizeof(*otf_data));
612
613 /* Convolution with the detector pixels */
614 /* The OTF is point symmetric so the whole image can be computed from the
615 values of a single octant. */
616 /* The image could be created with calloc() and j limited by
617 f_max / (mlam - mdlam * 0.5) but this is not faster */
618 for (j = 0; j <= pix0; j++) {
619 double sinc_y_9 = 0.0; /* Avoid uninit warning */
620 for (i = 0; i <= j; i++) {
621 if (i == 0 && j == 0) {
622 otf_data[size * pix0 + pix0] = 1.0;
623 }
624 else {
625 const double x = (double)i;
626 const double y = (double)j;
627 const double sqdist = x * x + y * y;
628 double f_lambda = 0.0;
629 double sinc_xy_9 = 0.0; /* Zero if OTF is zero */
630 double otfxy = 0.0;
631 int k;
632
633 assert(j > 0);
634
635 /* 9 iterations on the wavelength */
636 /* Unrolling the loop is not faster (due to the break?) */
637 for (k = 4; k >= -4; k--) {
638 /* Compute intermediate cut-off frequency */
639 const double lambda = lam - dlam * (double)k / 8.0;
640
641 /* A decreasing k ensures that we either enter on the first
642 iteration or not at all */
643 if (sqdist * lambda * lambda >= f_max * f_max)
644 break;
645
646 if (k == 4) {
647 f_lambda = sqrt(sqdist) / f_max;
648 if (i == 0) {
649 /* Sinc(x = 0) == 1 */
650 sinc_xy_9 = sinc_y_9 =
651 PSF_sinc_norm(y / (double)size) / 9.0;
652 }
653 else {
654 sinc_xy_9 =
655 sinc_y_9 * PSF_sinc_norm(x / (double)size);
656 }
657 }
658
659 otfxy += PSF_TelOTF(f_lambda * lambda, obs_ratio);
660 }
661 otfxy *= sinc_xy_9;
662
663 /* When i == j the same value is written to the same
664 position twice. That's probably faster than a guard */
665 otf_data[size * (pix0 - j) + pix0 - i] = otfxy;
666 otf_data[size * (pix0 - i) + pix0 - j] = otfxy;
667 if (i < pix0) {
668 otf_data[size * (pix0 - j) + pix0 + i] = otfxy;
669 otf_data[size * (pix0 + i) + pix0 - j] = otfxy;
670 if (j < pix0) {
671 otf_data[size * (pix0 + j) + pix0 - i] = otfxy;
672 otf_data[size * (pix0 - i) + pix0 + j] = otfxy;
673 otf_data[size * (pix0 + j) + pix0 + i] = otfxy;
674 otf_data[size * (pix0 + i) + pix0 + j] = otfxy;
675 }
676 }
677 }
678 }
679 }
680
681 return cpl_image_wrap_double(size, size, otf_data);
682}
683
684/*----------------------------------------------------------------------------*
685 * H1 function
686 *----------------------------------------------------------------------------*/
687static double
688PSF_H1(double f, double u, double v)
689{
690 const double e = fabs(1.0 - v) > 0.0 ? -1.0 : 1.0; /* e = 1.0 iff v = 1.0 */
691
692 return ((v * v / CPL_MATH_PI) *
693 acos((f / v) * (1.0 + e * (1.0 - u * u) / (4.0 * f * f))));
694}
695
696/*----------------------------------------------------------------------------*
697 * H2 function
698 *----------------------------------------------------------------------------*/
699static double
700PSF_H2(double f, double u)
701{
702 const double tmp1 = (2.0 * f) / (1.0 + u);
703 const double tmp2 = (1.0 - u) / (2.0 * f);
704
705 return -1.0 * (f / CPL_MATH_PI) * (1.0 + u) *
706 sqrt((1.0 - tmp1 * tmp1) * (1.0 - tmp2 * tmp2));
707}
708
709/*----------------------------------------------------------------------------*
710 * G function
711 *----------------------------------------------------------------------------*/
712static double
713PSF_G(double f, double u)
714{
715 if (f <= (1.0 - u) / 2.0)
716 return (u * u);
717 if (f >= (1.0 + u) / 2.0)
718 return (0.0);
719 else
720 return (PSF_H1(f, u, 1.0) + PSF_H1(f, u, u) + PSF_H2(f, u));
721}
722
723/*----------------------------------------------------------------------------*/
731/*----------------------------------------------------------------------------*/
732static double
733PSF_sinc_norm(double x)
734{
735 /* This static function should not be called with zero, but handle it
736 anyway. */
737 return x != 0.0 ? sin(x * CPL_MATH_PI) / (x * CPL_MATH_PI) : 1.0;
738}
739
740/*----------------------------------------------------------------------------*
741 * Telescope OTF function
742 *----------------------------------------------------------------------------*/
743static double
744PSF_TelOTF(double f, double u)
745{
746 return ((PSF_G(f, 1.0) + u * u * PSF_G(f / u, 1.0) - 2.0 * PSF_G(f, u)) /
747 (1.0 - u * u));
748}
749
750/*----------------------------------------------------------------------------*/
761/*----------------------------------------------------------------------------*/
762cpl_error_code
763irplib_strehl_disk_max(const cpl_image *self,
764 double xpos,
765 double ypos,
766 double radius,
767 double *ppeak)
768{
769 const int nx = cpl_image_get_size_x(self);
770 const int ny = cpl_image_get_size_y(self);
771 /* Round down */
772 const int lx = (int)(xpos - radius);
773 const int ly = (int)(ypos - radius);
774 /* Round up */
775 const int ux = (int)(xpos + radius) + 1;
776 const int uy = (int)(ypos + radius) + 1;
777
778 const double sqr = radius * radius;
779 cpl_boolean first = CPL_TRUE;
780 int i, j;
781
782
783 /* Check entries */
784 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
785 cpl_ensure_code(ppeak != NULL, CPL_ERROR_NULL_INPUT);
786 cpl_ensure_code(radius > 0.0, CPL_ERROR_ILLEGAL_INPUT);
787
788
789 for (j = IRPLIB_MAX(ly, 0); j < IRPLIB_MIN(uy, ny); j++) {
790 const double yj = (double)j - ypos;
791 for (i = IRPLIB_MAX(lx, 0); i < IRPLIB_MIN(ux, nx); i++) {
792 const double xi = (double)i - xpos;
793 const double dist = yj * yj + xi * xi;
794 if (dist <= sqr) {
795 int isbad;
796 const double value = cpl_image_get(self, i + 1, j + 1, &isbad);
797
798 if (!isbad && (first || value > *ppeak)) {
799 first = CPL_FALSE;
800 *ppeak = value;
801 }
802 }
803 }
804 }
805
806 return first ? cpl_error_set(cpl_func, CPL_ERROR_DATA_NOT_FOUND)
807 : CPL_ERROR_NONE;
808}
809
810#ifndef IRPLIB_NO_FIT_GAUSSIAN
811#ifdef IRPLIB_STREHL_USE_CPL_IMAGE_FIT_GAUSSIAN
812/*----------------------------------------------------------------------------*/
828/*----------------------------------------------------------------------------*/
829static double
830irplib_gaussian_2d(double x, double y, double norm, double sig_x, double sig_y)
831{
832 /* Copied from CPL */
833 return norm /
834 (sig_x * sig_y * CPL_MATH_2PI *
835 exp(x * x / (2.0 * sig_x * sig_x) + y * y / (2.0 * sig_y * sig_y)));
836}
837#endif
838
839#if defined CPL_VERSION_CODE && CPL_VERSION_CODE >= CPL_VERSION(6, 9, 1)
840#else
841/*----------------------------------------------------------------------------*/
860/*----------------------------------------------------------------------------*/
861static double
862irplib_gaussian_eval_2d(const cpl_array *self, double x, double y)
863{
864 cpl_errorstate prestate = cpl_errorstate_get();
865 const double B = cpl_array_get_double(self, 0, NULL);
866 const double A = cpl_array_get_double(self, 1, NULL);
867 const double R = cpl_array_get_double(self, 2, NULL);
868 const double M_x = cpl_array_get_double(self, 3, NULL);
869 const double M_y = cpl_array_get_double(self, 4, NULL);
870 const double S_x = cpl_array_get_double(self, 5, NULL);
871 const double S_y = cpl_array_get_double(self, 6, NULL);
872
873 double value = 0.0;
874
875 if (!cpl_errorstate_is_equal(prestate)) {
876 (void)cpl_error_set_where(cpl_func);
877 }
878 else if (cpl_array_get_size(self) != 7) {
879 (void)cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_INPUT);
880 }
881 else if (fabs(R) < 1.0 && S_x != 0.0 && S_y != 0.0) {
882 const double x_n = (x - M_x) / S_x;
883 const double y_n = (y - M_y) / S_y;
884
885 value = B + A / (CPL_MATH_2PI * S_x * S_y * sqrt(1 - R * R)) *
886 exp(-0.5 / (1 - R * R) *
887 (x_n * x_n + y_n * y_n - 2.0 * R * x_n * y_n));
888 }
889 else if (fabs(R) > 1.0) {
890 (void)cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT,
891 "fabs(R=%g) > 1", R);
892 }
893 else {
894 (void)cpl_error_set_message(cpl_func, CPL_ERROR_DIVISION_BY_ZERO,
895 "R=%g. Sigma=(%g, %g)", R, S_x, S_y);
896 }
897
898 return value;
899}
900#endif
901
902/*----------------------------------------------------------------------------*/
909/*----------------------------------------------------------------------------*/
910static uint32_t
911irplib_roundup_power2(uint32_t v)
912{
913 v |= v >> 1;
914 v |= v >> 2;
915 v |= v >> 4;
916 v |= v >> 8;
917 v |= v >> 16;
918
919 return v + 1;
920}
921
922
923static cpl_error_code
924irplib_closeset_aperture(const cpl_apertures *self,
925 const double x,
926 const double y,
927 int *ind)
928{
929 double dist = INFINITY;
930
931 const int nsize = cpl_apertures_get_size(self);
932 int i;
933
934 cpl_ensure_code(nsize > 0, cpl_error_get_code());
935 cpl_ensure_code(ind, CPL_ERROR_NULL_INPUT);
936
937 *ind = -1;
938
939 for (i = 1; i <= nsize; ++i) {
940 const double this_x = cpl_apertures_get_centroid_x(self, i);
941 const double this_y = cpl_apertures_get_centroid_y(self, i);
942 const double this_d_sq = pow(this_x - x, 2.0) + pow(this_y - y, 2.0);
943
944 if (this_d_sq < dist) {
945 dist = this_d_sq;
946 *ind = i;
947 }
948 }
949 return CPL_ERROR_NONE;
950}
951
952/*----------------------------------------------------------------------------*/
963/*----------------------------------------------------------------------------*/
964static cpl_error_code
965irplib_gaussian_maxpos(const cpl_image *self,
966 double sigma,
967 double x_initial,
968 double y_initial,
969 double *pxpos,
970 double *pypos,
971 double *ppeak)
972{
973 const cpl_size nx = cpl_image_get_size_x(self);
974 const cpl_size ny = cpl_image_get_size_y(self);
975 int iretry = 3; /* Number retries with decreasing sigma */
976 int ifluxapert = 0;
977 double med_dist;
978 const double median = cpl_image_get_median_dev(self, &med_dist);
979 cpl_mask *selection;
980 cpl_size nlabels = 0;
981 cpl_image *labels = NULL;
982 cpl_apertures *aperts;
983 cpl_size npixobj;
984 double objradius;
985 cpl_size winsize;
986 cpl_size xposmax, yposmax;
987 double xposcen, yposcen;
988 double valmax, valfit = -1.0;
989#ifdef IRPLIB_STREHL_USE_CPL_IMAGE_FIT_GAUSSIAN
990 double norm, xcen, ycen, sig_x, sig_y, fwhm_x, fwhm_y;
991#endif
992 cpl_array *gauss_parameters = NULL;
993 cpl_errorstate prestate = cpl_errorstate_get();
994 cpl_error_code code = CPL_ERROR_NONE;
995
996
997 cpl_ensure_code(sigma > 0.0, CPL_ERROR_ILLEGAL_INPUT);
998
999 selection = cpl_mask_new(nx, ny);
1000
1001 for (; iretry > 0 && nlabels == 0; iretry--, sigma *= 0.5) {
1002 /* Compute the threshold */
1003 const double threshold = median + sigma * med_dist;
1004
1005
1006 /* Select the pixel above the threshold */
1007 code = cpl_mask_threshold_image(selection, self, threshold, DBL_MAX,
1008 CPL_BINARY_1);
1009
1010 if (code)
1011 break;
1012
1013 /* Labelise the thresholded selection */
1014 cpl_image_delete(labels);
1015 labels = cpl_image_labelise_mask_create(selection, &nlabels);
1016 }
1017 sigma *= 2.0; /* FIXME: unelegant */
1018
1019 cpl_mask_delete(selection);
1020
1021 if (code) {
1022 cpl_image_delete(labels);
1023 return cpl_error_set_where(cpl_func);
1024 }
1025 else if (nlabels == 0) {
1026 cpl_image_delete(labels);
1027 return cpl_error_set(cpl_func, CPL_ERROR_DATA_NOT_FOUND);
1028 }
1029
1030 aperts = cpl_apertures_new_from_image(self, labels);
1031
1032 /* Find the aperture closest to the provided coordinates */
1033 code = irplib_closeset_aperture(aperts, x_initial, y_initial, &ifluxapert);
1034
1035 if (code) {
1036 cpl_apertures_delete(aperts);
1037 cpl_image_delete(labels);
1038 return cpl_error_set(cpl_func, CPL_ERROR_DATA_NOT_FOUND);
1039 }
1040
1041 npixobj = cpl_apertures_get_npix(aperts, ifluxapert);
1042 objradius = sqrt((double)npixobj * CPL_MATH_1_PI);
1043 /* Size is power of two for future noise filtering w. fft */
1044 winsize =
1045 IRPLIB_MIN(IRPLIB_MIN(nx, ny),
1046 irplib_roundup_power2((uint32_t)(3.0 * objradius + 0.5)));
1047
1048 xposmax = cpl_apertures_get_maxpos_x(aperts, ifluxapert);
1049 yposmax = cpl_apertures_get_maxpos_y(aperts, ifluxapert);
1050 xposcen = cpl_apertures_get_centroid_x(aperts, ifluxapert);
1051 yposcen = cpl_apertures_get_centroid_y(aperts, ifluxapert);
1052 valmax = cpl_apertures_get_max(aperts, ifluxapert);
1053
1054 cpl_apertures_delete(aperts);
1055 cpl_image_delete(labels);
1056
1057 cpl_msg_debug(cpl_func, "Object radius at S/R=%g: %g (window-size=%u)",
1058 sigma, objradius, (unsigned)winsize);
1059 cpl_msg_debug(cpl_func, "Object-peak @ (%d, %d) = %g", (int)xposmax,
1060 (int)yposmax, valmax);
1061
1062 gauss_parameters = cpl_array_new(7, CPL_TYPE_DOUBLE);
1063 cpl_array_set_double(gauss_parameters, 0, median);
1064
1065 code = cpl_fit_image_gaussian(self, NULL, xposcen, yposcen, winsize,
1066 winsize, gauss_parameters, NULL, NULL, NULL,
1067 NULL, NULL, NULL, NULL, NULL, NULL);
1068 if (!code) {
1069 const double M_x = cpl_array_get_double(gauss_parameters, 3, NULL);
1070 const double M_y = cpl_array_get_double(gauss_parameters, 4, NULL);
1071
1072 valfit = irplib_gaussian_eval_2d(gauss_parameters, M_x, M_y);
1073
1074 if (!cpl_errorstate_is_equal(prestate)) {
1075 code = cpl_error_get_code();
1076 }
1077 else {
1078 *pxpos = M_x;
1079 *pypos = M_y;
1080 *ppeak = valfit;
1081
1082 cpl_msg_debug(cpl_func, "Gauss-fit @ (%g, %g) = %g", M_x, M_y,
1083 valfit);
1084 }
1085 }
1086 cpl_array_delete(gauss_parameters);
1087
1088#ifdef IRPLIB_STREHL_USE_CPL_IMAGE_FIT_GAUSSIAN
1089 if (code || valfit < valmax) {
1090 cpl_errorstate_set(prestate);
1091
1092 code = cpl_image_fit_gaussian(self, xposcen, yposcen,
1093 (int)(2.0 * objradius), &norm, &xcen,
1094 &ycen, &sig_x, &sig_y, &fwhm_x, &fwhm_y);
1095
1096 if (!code) {
1097 valfit = irplib_gaussian_2d(0.0, 0.0, norm, sig_x, sig_y);
1098
1099 cpl_msg_debug(cpl_func,
1100 "Gauss-Fit @ (%g, %g) = %g. norm=%g, "
1101 "sigma=(%g, %g)",
1102 xcen, ycen, valfit, norm, sig_x, sig_y);
1103
1104 if (valfit > valmax) {
1105 *pxpos = xcen;
1106 *pypos = ycen;
1107 *ppeak = valfit;
1108 }
1109 }
1110 }
1111#endif
1112
1113 if (code || valfit < valmax) {
1114 cpl_errorstate_set(prestate);
1115 *pxpos = xposcen;
1116 *pypos = yposcen;
1117 *ppeak = valmax;
1118 }
1119
1120 return code ? cpl_error_set_where(cpl_func) : CPL_ERROR_NONE;
1121}
1122#endif
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2