00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031
00032
00033
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
00054
00055
00056 #define IRPLIB_PI 3.1415926535897932384626433
00057
00058
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
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
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
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
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
00170
00171
00172
00173 psf = irplib_strehl_generate_psf(m1, m2, lam, dlam, pscale, size);
00174 cpl_ensure_code(psf != NULL, CPL_ERROR_ILLEGAL_OUTPUT);
00175
00176
00177 *psf_peak = cpl_image_get_max(psf);
00178 cpl_image_delete(psf);
00179
00180 assert( *psf_peak > 0.0);
00181 *psf_flux = 1.0;
00182
00183
00184 *star_bg = irplib_strehl_ring_background(im, xpos, ypos, r2/pscale, r3/pscale,
00185 IRPLIB_BG_METHOD_AVER_REJ);
00186
00187
00188 star_radius = r1/pscale;
00189
00190
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
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
00204
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
00213
00214
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
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
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
00276 lx = (int)(xpos - rad);
00277 ly = (int)(ypos - rad);
00278 if (lx < 0) lx = 0;
00279 if (ly < 0) ly = 0;
00280
00281
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
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
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
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
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
00366 pix_arr = cpl_vector_new(npix);
00367
00368
00369
00370
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
00391
00392
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
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 {
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
00451
00452
00453
00454
00455
00456
00457
00458
00459 if (cpl_image_fft(otf_image, NULL, CPL_FFT_UNNORMALIZED) ||
00460
00461
00462 cpl_image_abs(otf_image) ||
00463
00464
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 ;
00504 double f_max ;
00505 int pix0 ;
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
00522 pscale /= (double)206265;
00523 lam /= (double)1.0e6;
00524 dlam /= (double)1.0e6;
00525
00526
00527 obs_ratio = m2 / m1;
00528
00529
00530 pix0 = size/2;
00531 invsize = (double)1.0 / (double)size;
00532
00533
00534 f_max = m1 * pscale * (double)size / lam;
00535
00536
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
00542
00543
00544 for (k=1 ; k<=9 ; k++) {
00545
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
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
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;
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
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
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
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
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;
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
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
00672 lx = (int)(xpos - radius);
00673 ly = (int)(ypos - radius);
00674 if (lx < 0) lx = 0;
00675 if (ly < 0) ly = 0;
00676
00677
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
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 }