irplib_strehl.c

00001 /* $Id: irplib_strehl.c,v 1.26 2006/11/21 09:49:01 llundin Exp $
00002  *
00003  * This file is part of the irplib package
00004  * Copyright (C) 2002,2003 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00019  */
00020 
00021 /*
00022  * $Author: llundin $
00023  * $Date: 2006/11/21 09:49:01 $
00024  * $Revision: 1.26 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 /*-----------------------------------------------------------------------------
00033                                    Includes
00034  -----------------------------------------------------------------------------*/
00035 
00036 #include <string.h>
00037 #include <assert.h>
00038 #include <math.h>
00039 #include <float.h>
00040 
00041 #include <cpl.h>
00042 
00043 #include "irplib_utils.h"
00044 #include "irplib_strehl.h"
00045 
00046 /*----------------------------------------------------------------------------*/
00050 /*----------------------------------------------------------------------------*/
00051 
00052 /*-----------------------------------------------------------------------------
00053                                    Define
00054  -----------------------------------------------------------------------------*/
00055 
00056 #define IRPLIB_PI                   3.1415926535897932384626433
00057 
00058 /* Determined empirically by C. Lidman for Strehl error computation */
00059 #define STREHL_ERROR_COEFFICIENT    (IRPLIB_PI * 0.007 / 0.0271)
00060 
00061 #ifndef IRPLIB_STREHL_RAD_CENTRAL
00062 #define IRPLIB_STREHL_RAD_CENTRAL 5
00063 #endif
00064 
00065 #define IRPLIB_DISK_BG_MIN_PIX_NB    30
00066 #define IRPLIB_DISK_BG_REJ_LOW       0.1
00067 #define IRPLIB_DISK_BG_REJ_HIGH      0.1
00068 
00069 /*-----------------------------------------------------------------------------
00070                                    Functions prototypes
00071  -----------------------------------------------------------------------------*/
00072 
00073 static cpl_image * irplib_strehl_generate_otf(double, double, double, double, 
00074                                               int, double);
00075 static double PSF_H1(double, double, double);
00076 static double PSF_H2(double, double);
00077 static double PSF_G(double, double);
00078 static double PSF_sinc(double);
00079 static double PSF_TelOTF(double, double);
00080 
00081 static cpl_error_code irplib_strehl_disk_max(const cpl_image *, double, double,
00082                                              double, double *);
00083 
00084 
00085 /*-----------------------------------------------------------------------------
00086                                    Functions code
00087  -----------------------------------------------------------------------------*/
00090 /*----------------------------------------------------------------------------*/
00119 /*----------------------------------------------------------------------------*/
00120 cpl_error_code irplib_strehl_compute(const cpl_image *   im,
00121                                      double              m1,
00122                                      double              m2,
00123                                      double              lam,
00124                                      double              dlam,
00125                                      double              pscale,
00126                                      int                 size,
00127                                      double              xpos,
00128                                      double              ypos,
00129                                      double              r1,
00130                                      double              r2,
00131                                      double              r3,
00132                                      int                 noise_box_sz,
00133                                      int                 noise_nsamples,
00134                                      double          *   strehl,
00135                                      double          *   strehl_err,
00136                                      double          *   star_bg,
00137                                      double          *   star_peak,
00138                                      double          *   star_flux,
00139                                      double          *   psf_peak,
00140                                      double          *   psf_flux,
00141                                      double          *   bg_noise)
00142 {
00143     cpl_image  * psf;
00144     double       star_radius, max_radius;
00145     /* FIXME: Arbitrary choice of image border */
00146     const double window_size = (double)(IRPLIB_STREHL_RAD_CENTRAL);
00147     int          ring[4];
00148 
00149 
00150     assert(window_size > 0.0);
00151     assert(STREHL_ERROR_COEFFICIENT > 0.0);
00152 
00153     /* Test inputs */
00154     cpl_ensure_code(im != NULL,         CPL_ERROR_NULL_INPUT);
00155     cpl_ensure_code(strehl != NULL,     CPL_ERROR_NULL_INPUT);
00156     cpl_ensure_code(strehl_err != NULL, CPL_ERROR_NULL_INPUT);
00157     cpl_ensure_code(star_bg != NULL,    CPL_ERROR_NULL_INPUT);
00158     cpl_ensure_code(star_peak != NULL,  CPL_ERROR_NULL_INPUT);
00159     cpl_ensure_code(star_flux != NULL,  CPL_ERROR_NULL_INPUT);
00160     cpl_ensure_code(psf_peak != NULL,   CPL_ERROR_NULL_INPUT);
00161     cpl_ensure_code(psf_flux != NULL,   CPL_ERROR_NULL_INPUT);
00162 
00163     cpl_ensure_code(pscale > 0.0,      CPL_ERROR_ILLEGAL_INPUT);
00164 
00165     cpl_ensure_code(r1 > 0.0,      CPL_ERROR_ILLEGAL_INPUT);
00166     cpl_ensure_code(r2 > 0.0,      CPL_ERROR_ILLEGAL_INPUT);
00167     cpl_ensure_code(r3 > r2,       CPL_ERROR_ILLEGAL_INPUT);
00168 
00169     /* Computing a Strehl ratio is a story between an ideal PSF */
00170     /* and a candidate image supposed to approximate this ideal PSF. */
00171 
00172     /* Generate first appropriate PSF to find max peak */
00173     psf = irplib_strehl_generate_psf(m1, m2, lam, dlam, pscale, size);
00174     cpl_ensure_code(psf != NULL,      CPL_ERROR_ILLEGAL_OUTPUT);
00175 
00176     /* Compute flux in PSF and find max peak */
00177     *psf_peak = cpl_image_get_max(psf);
00178     cpl_image_delete(psf);
00179 
00180     assert( *psf_peak > 0.0); /* The ideal PSF has a positive maximum */
00181     *psf_flux = 1.0; /* The psf flux, cpl_image_get_flux(psf), is always 1 */
00182 
00183     /* Measure the background in the candidate image */
00184     *star_bg = irplib_strehl_ring_background(im, xpos, ypos, r2/pscale, r3/pscale,
00185                                              IRPLIB_BG_METHOD_AVER_REJ);
00186 
00187     /* Compute star_radius in pixels */
00188     star_radius = r1/pscale;
00189 
00190     /* Measure the flux on the candidate image */
00191     *star_flux = irplib_strehl_disk_flux(im, xpos, ypos, star_radius, *star_bg);
00192 
00193     cpl_ensure_code(*star_flux > 0.0,      CPL_ERROR_ILLEGAL_OUTPUT);
00194 
00195     /* Find the peak value on the central part of the candidate image */
00196     max_radius = window_size < star_radius ? window_size : star_radius;
00197     cpl_ensure_code(!irplib_strehl_disk_max(im, xpos, ypos, max_radius,
00198                                             star_peak), cpl_error_get_code());
00199     *star_peak -= *star_bg;
00200 
00201     cpl_ensure_code(*star_peak > 0.0,      CPL_ERROR_ILLEGAL_OUTPUT);
00202 
00203     /* Compute Strehl */
00204     /* (StarPeak / StarFlux) / (PsfPeak / PsfFlux) */
00205     *strehl = (*star_peak * *psf_flux ) / ( *star_flux * *psf_peak);
00206 
00207     if (*strehl > 1)
00208         cpl_msg_warning(cpl_func, "Extreme Strehl-ratio=%g, star_peak=%g, "
00209                         "star_flux=%g, psf_peak=%g, psf_flux=%g", *strehl,
00210                         *star_peak, *star_flux, *psf_peak, *psf_flux);
00211 
00212     /* Compute Strehl error */
00213     /* FIXME: cpl_flux_get_noise_ring() should take a ring specification
00214        using doubles */
00215     ring[0] = (int)(0.5 + xpos);
00216     ring[1] = (int)(0.5 + ypos);
00217     ring[2] = (int)(0.5 + r2/pscale);
00218     ring[3] = (int)(0.5 + r3/pscale);
00219 
00220     cpl_ensure_code(!cpl_flux_get_noise_ring(im, ring, noise_box_sz,
00221                                              noise_nsamples, bg_noise, NULL),
00222                     cpl_error_get_code());
00223 
00224     *strehl_err = STREHL_ERROR_COEFFICIENT * (*bg_noise) * pscale * 
00225         star_radius * star_radius / *star_flux;
00226 
00227     /* This check should not be able to fail, but just to be sure */
00228     cpl_ensure_code(*strehl_err >= 0.0,       CPL_ERROR_ILLEGAL_OUTPUT);
00229 
00230     return CPL_ERROR_NONE;
00231 }
00232 
00233 /*----------------------------------------------------------------------------*/
00249 /*----------------------------------------------------------------------------*/
00250 double irplib_strehl_disk_flux(const cpl_image * im,
00251                                double            xpos,
00252                                double            ypos,
00253                                double            rad,
00254                                double            bg)
00255 {
00256     const double    sqr = rad * rad;
00257     double          sqrest;
00258     const float *   pim;
00259     double          flux = 0.0;
00260     double          yj, xi; 
00261     int             nx, ny;
00262     int             lx, ly, ux, uy;
00263     int             i, j;
00264 
00265 
00266     /* Check entries */
00267     cpl_ensure(im != NULL, CPL_ERROR_NULL_INPUT, 0.0);
00268     cpl_ensure(cpl_image_get_type(im) == CPL_TYPE_FLOAT,
00269                CPL_ERROR_UNSUPPORTED_MODE, 0.0);
00270     cpl_ensure(rad > 0.0, CPL_ERROR_ILLEGAL_INPUT, 0.0);
00271  
00272     nx = cpl_image_get_size_x(im);
00273     ny = cpl_image_get_size_y(im);
00274 
00275     /* Round down */
00276     lx = (int)(xpos - rad);
00277     ly = (int)(ypos - rad);
00278     if (lx < 0) lx = 0;
00279     if (ly < 0) ly = 0;
00280 
00281     /* Round up */
00282     ux = (int)(xpos + rad) + 1;
00283     uy = (int)(ypos + rad) + 1;
00284     if (ux > (nx-1)) ux = nx-1;
00285     if (uy > (ny-1)) uy = ny-1;
00286 
00287     /* im is actually not modified */
00288     pim = cpl_image_get_data_float((cpl_image*)im);
00289     for (j=ly ; j<uy ; j++) {
00290         yj = (double)j - ypos;
00291         sqrest = sqr - yj * yj;
00292         for (i=lx; i<ux ; i++) {
00293             xi = (double)i - xpos;
00294             if (sqrest >= xi * xi && isnan(pim[i+j*nx]) == 0) {
00295                 flux += (double)pim[i+j*nx] - bg;
00296             }
00297         }
00298     }
00299     return flux;
00300 }
00301 
00302 /*----------------------------------------------------------------------------*/
00316 /*----------------------------------------------------------------------------*/
00317 double irplib_strehl_ring_background(const cpl_image *       im,
00318                                      double                  xpos,
00319                                      double                  ypos,
00320                                      double                  rad_int,
00321                                      double                  rad_ext,
00322                                      irplib_strehl_bg_method mode)
00323 {
00324     int             npix;
00325     const double    sqr_int = rad_int * rad_int;
00326     const double    sqr_ext = rad_ext * rad_ext;
00327     double          dist;
00328     cpl_vector  *   pix_arr;
00329     const float *   pim;
00330     double          flux = 0.0;
00331     double          yj, xi; 
00332     int             lx, ly, ux, uy;
00333     int             nx, ny;
00334     int             i, j;
00335 
00336     /* Check entries */
00337     cpl_ensure(im != NULL, CPL_ERROR_NULL_INPUT, 0.0);
00338     cpl_ensure(cpl_image_get_type(im) == CPL_TYPE_FLOAT,
00339                CPL_ERROR_UNSUPPORTED_MODE, 0.0);
00340     cpl_ensure(rad_int > 0.0, CPL_ERROR_ILLEGAL_INPUT, 0.0);
00341     cpl_ensure(rad_ext > rad_int, CPL_ERROR_ILLEGAL_INPUT, 0.0);
00342 
00343     cpl_ensure(mode == IRPLIB_BG_METHOD_AVER_REJ ||
00344                mode == IRPLIB_BG_METHOD_MEDIAN,
00345                CPL_ERROR_UNSUPPORTED_MODE, 0.0);
00346 
00347     nx = cpl_image_get_size_x(im);
00348     ny = cpl_image_get_size_y(im);
00349 
00350     /* Round down */
00351     lx = (int)(xpos - rad_ext);
00352     ly = (int)(ypos - rad_ext);
00353     if (lx < 0) lx = 0;
00354     if (ly < 0) ly = 0;
00355 
00356     /* Round up */
00357     ux = (int)(xpos + rad_ext) + 1;
00358     uy = (int)(ypos + rad_ext) + 1;
00359     if (ux > (nx-1)) ux = nx-1;
00360     if (uy > (ny-1)) uy = ny-1;
00361 
00362     npix = (ux - lx + 1) * (uy - ly + 1);
00363     cpl_ensure(npix >= IRPLIB_DISK_BG_MIN_PIX_NB, CPL_ERROR_DATA_NOT_FOUND, 0.0);
00364 
00365     /* Allocate pixel array to hold values in the ring */
00366     pix_arr = cpl_vector_new(npix);
00367 
00368     /* Count number of pixels in the ring */
00369     /* Retrieve all pixels which belong to the ring */
00370     /* im is actually _not_ modified */
00371     pim = cpl_image_get_data_float((cpl_image*)im);
00372     npix = 0;
00373     for (j=ly ; j<uy ; j++) {
00374         yj = (double)j - ypos;
00375         for (i=lx ; i<ux; i++) {
00376             xi = (double)i - xpos;
00377             dist = yj * yj + xi * xi;
00378             if (sqr_int <= dist && dist <= sqr_ext && isnan(pim[i+j*nx]) == 0) {
00379                 cpl_vector_set(pix_arr, npix, (double)pim[i+j*nx]);
00380                 npix++;
00381             }
00382         }
00383     }
00384 
00385     if (npix < IRPLIB_DISK_BG_MIN_PIX_NB) {
00386         cpl_vector_delete(pix_arr);
00387         cpl_ensure(0, CPL_ERROR_DATA_NOT_FOUND, 0.0);
00388     }
00389 
00390     /* Should not be able to fail now */
00391 
00392     /* Resize pixel array to actual number of values within the ring */
00393     cpl_vector_set_size(pix_arr, npix);
00394 
00395     if (mode == IRPLIB_BG_METHOD_AVER_REJ) {
00396         const int low_ind  = (int)((double)npix * IRPLIB_DISK_BG_REJ_LOW);
00397         const int high_ind = (int)((double)npix
00398                                    * (1.0 - IRPLIB_DISK_BG_REJ_HIGH));
00399 
00400         /* Sort the array */
00401         cpl_vector_sort(pix_arr, 1);
00402 
00403         for (i=low_ind ; i<high_ind ; i++) {
00404             flux += cpl_vector_get(pix_arr, i);
00405         }
00406         if (high_ind - low_ind > 1) flux /= (double)(high_ind - low_ind);
00407     } else /* if (mode == IRPLIB_BG_METHOD_MEDIAN) */ {
00408         flux = cpl_vector_get_median(pix_arr);
00409     }
00410 
00411     cpl_vector_delete(pix_arr);
00412 
00413     return flux;
00414 }
00415 
00416 /*----------------------------------------------------------------------------*/
00436 /*----------------------------------------------------------------------------*/
00437 cpl_image * irplib_strehl_generate_psf(
00438         double   m1,
00439         double   m2,
00440         double   lam,
00441         double   dlam,
00442         double   pscale,
00443         int      size)
00444 {
00445     cpl_image * otf_image = irplib_strehl_generate_otf(m1, m2, lam, dlam, 
00446             size, pscale);
00447 
00448     if (otf_image == NULL) return NULL;
00449 
00450     /* Transform back to real space
00451        - Normalization is unnecessary, due to the subsequent normalisation.
00452        - An OTF is point symmetric about its center, i.e. it is even,
00453          i.e. the real space image is real.
00454        - Because of this a forward FFT works as well.
00455        - If the PSF ever needs to have its images halves swapped add
00456          CPL_FFT_SWAP_HALVES to the FFT call.
00457      */
00458 
00459     if (cpl_image_fft(otf_image, NULL, CPL_FFT_UNNORMALIZED) ||
00460 
00461         /* Compute absolute values of PSF */
00462         cpl_image_abs(otf_image) ||
00463 
00464         /* Normalize PSF to get flux=1  */
00465         cpl_image_normalise(otf_image, CPL_NORM_FLUX)) {
00466 
00467         cpl_image_delete(otf_image);
00468         return NULL;
00469     }
00470     
00471     return otf_image;
00472 }
00473 
00476 /*----------------------------------------------------------------------------*/
00492 /*----------------------------------------------------------------------------*/
00493 static cpl_image * irplib_strehl_generate_otf(
00494         double  m1,
00495         double  m2,
00496         double  lam,
00497         double  dlam,
00498         int     size,
00499         double  pscale)
00500 {
00501     cpl_image   *   otf_image;
00502     double      *   otf_data;
00503     double          obs_ratio ;  /* m1 / m2    */
00504     double          f_max ;      /* cut-off frequency        */
00505     int             pix0 ;      /* Pixel corresponding to the zero frequency */
00506     double          a, x, y;
00507     double          f, rsq, fc, invfc, lambda;
00508     double          sincy;
00509     double          invsize;
00510     register int    pos;
00511     int             i, j, k;
00512 
00513 
00514     cpl_ensure(m2   > 0.0,      CPL_ERROR_ILLEGAL_INPUT, NULL);
00515     cpl_ensure(m1   > m2,       CPL_ERROR_ILLEGAL_INPUT, NULL);
00516     cpl_ensure(lam  > 0.0,      CPL_ERROR_ILLEGAL_INPUT, NULL);
00517     cpl_ensure(dlam > 0.0,      CPL_ERROR_ILLEGAL_INPUT, NULL);
00518     cpl_ensure(size > 0,        CPL_ERROR_ILLEGAL_INPUT, NULL);
00519     cpl_ensure(pscale > 0.0,    CPL_ERROR_ILLEGAL_INPUT, NULL);
00520 
00521     /* Convert pixel scale from sec to radians, microns in meters    */
00522     pscale /= (double)206265;
00523     lam /= (double)1.0e6;
00524     dlam /= (double)1.0e6;
00525 
00526     /* Obscuration ratio    */
00527     obs_ratio = m2 / m1;
00528 
00529     /* Pixel corresponding to the zero frequency    */
00530     pix0 = size/2;
00531     invsize = (double)1.0 / (double)size;
00532 
00533     /* Cut-off frequency in pixels  */
00534     f_max = m1 * pscale * (double)size / lam;
00535 
00536     /* Allocate for output image    */
00537     otf_image = cpl_image_new(size, size, CPL_TYPE_DOUBLE);
00538     if (otf_image==NULL) return NULL;
00539     otf_data = cpl_image_get_data_double(otf_image);
00540 
00541     /* Now compute the OTF  */
00542     /* OPTIMIZED CODE !!! LIMITED READABILITY !!!   */
00543 
00544     for (k=1 ; k<=9 ; k++) {    /* iteration on the wavelength  */
00545         /* Compute intermediate cut-off frequency   */
00546         lambda = (double)(lam - dlam*(double)(k-5)/8.0);
00547         fc = (double)f_max * (double)lam / lambda;
00548         invfc = 1.0 / fc;
00549 
00550         /* Convolution with the detector pixels */
00551         pos = 0;
00552         for (j=0 ; j<size ; j++) {
00553             y = (double)(j-pix0);
00554             sincy = PSF_sinc(IRPLIB_PI * y * invsize);
00555             for (i=0 ; i<size ; i++) {
00556                 x = (double)(i-pix0);
00557                 rsq = x*x + y*y;
00558                 if (rsq < fc*fc) {
00559                     if (rsq < 0.01)
00560                         a = 1.0;
00561                     else {
00562                         f = sqrt(rsq) * invfc;
00563                         a = PSF_TelOTF(f,obs_ratio) *
00564                             PSF_sinc(IRPLIB_PI * x * invsize) * sincy;
00565                     }
00566                 } else {
00567                     a = 0.0;
00568                 }
00569                 otf_data[pos++] += a / 9.0;
00570             }
00571         }
00572     }
00573     return otf_image;
00574 }
00575 
00576 /*----------------------------------------------------------------------------*
00577  * H1 function
00578  *----------------------------------------------------------------------------*/
00579 static double PSF_H1(
00580         double  f,
00581         double  u,
00582         double  v)
00583 {
00584     const double e = fabs(1.0-v) > 0.0 ? -1.0 : 1.0; /* e = 1.0 iff v = 1.0 */
00585 
00586     return((v*v/IRPLIB_PI)*acos((f/v)*(1.0+e*(1.0-u*u)/(4.0*f*f))));
00587 }
00588 
00589 /*----------------------------------------------------------------------------*
00590  * H2 function
00591  *----------------------------------------------------------------------------*/
00592 static double PSF_H2(double  f,
00593                      double  u)
00594 {
00595     const double tmp1 = (2.0 * f) / (1.0 + u);
00596     const double tmp2 = (1.0 - u) / (2.0 * f);
00597 
00598     return -1.0 * (f/IRPLIB_PI) * (1.0+u)
00599         * sqrt((1.0-tmp1*tmp1)*(1.0-tmp2*tmp2));
00600 }   
00601 
00602 /*----------------------------------------------------------------------------*
00603  * G function
00604  *----------------------------------------------------------------------------*/
00605 static double PSF_G(double  f, 
00606                     double  u)
00607 {
00608     if (f <= (1.0-u)/2.0) return(u*u);
00609     if (f >= (1.0+u)/2.0) return(0.0);
00610     else return(PSF_H1(f,u,1.0) + PSF_H1(f,u,u) + PSF_H2(f,u));
00611 }   
00612 
00613 /*----------------------------------------------------------------------------*
00614  * sinc function
00615  *----------------------------------------------------------------------------*/
00616 static double PSF_sinc(double x)
00617 {
00618   return fabs(x) > fabs(sin(x)) ? sin(x)/x : 1.0;
00619 }
00620 
00621 /*----------------------------------------------------------------------------*
00622  * Telescope OTF function 
00623  *----------------------------------------------------------------------------*/
00624 static double PSF_TelOTF(double  f,  
00625                          double  u)
00626 {
00627     return((PSF_G(f,1.0)+u*u*PSF_G(f/u,1.0)-2.0*PSF_G(f,u))/(1.0-u*u));
00628 }
00629 
00630 /*----------------------------------------------------------------------------*/
00642 /*----------------------------------------------------------------------------*/
00643 static cpl_error_code irplib_strehl_disk_max(const cpl_image * self,
00644                                              double            xpos,
00645                                              double            ypos,
00646                                              double            radius,
00647                                              double          * ppeak)
00648 {
00649 
00650     const double    sqr = radius * radius;
00651     double          sqrest;
00652     const float *   pself;
00653     float           peak = FLT_MAX;  /* Avoid (false) uninit warning */
00654     double          yj, xi; 
00655     int             nx, ny;
00656     int             lx, ly, ux, uy;
00657     int             i, j;
00658     cpl_boolean     first = CPL_TRUE;
00659 
00660 
00661     /* Check entries */
00662     cpl_ensure_code(ppeak != NULL, CPL_ERROR_NULL_INPUT);
00663     cpl_ensure_code(self  != NULL, CPL_ERROR_NULL_INPUT);
00664     cpl_ensure_code(cpl_image_get_type(self) == CPL_TYPE_FLOAT,
00665                     CPL_ERROR_UNSUPPORTED_MODE);
00666     cpl_ensure_code(radius > 0.0, CPL_ERROR_ILLEGAL_INPUT);
00667  
00668     nx = cpl_image_get_size_x(self);
00669     ny = cpl_image_get_size_y(self);
00670 
00671     /* Round down */
00672     lx = (int)(xpos - radius);
00673     ly = (int)(ypos - radius);
00674     if (lx < 0) lx = 0;
00675     if (ly < 0) ly = 0;
00676 
00677     /* Round up */
00678     ux = (int)(xpos + radius) + 1;
00679     uy = (int)(ypos + radius) + 1;
00680     if (ux > (nx-1)) ux = nx-1;
00681     if (uy > (ny-1)) uy = ny-1;
00682 
00683     /* im is actually not modified */
00684     pself = cpl_image_get_data_float((cpl_image*)self);
00685     for (j=ly ; j<uy ; j++) {
00686         yj = (double)j - ypos;
00687         sqrest = sqr - yj * yj;
00688         for (i=lx; i<ux ; i++) {
00689             xi = (double)i - xpos;
00690             if (sqrest >= xi * xi && isnan(pself[i+j*nx]) == 0) {
00691                 if (first || pself[i+j*nx] > peak) {
00692                     first = CPL_FALSE;
00693                     peak = pself[i+j*nx];
00694                 }
00695             }
00696         }
00697     }
00698 
00699     cpl_ensure_code(!first, CPL_ERROR_DATA_NOT_FOUND);
00700 
00701     *ppeak = (double)peak;
00702 
00703     return CPL_ERROR_NONE;
00704 }

Generated on Wed Jan 17 08:33:41 2007 for SINFONI Pipeline Reference Manual by  doxygen 1.4.4