fors_img_science_impl.c

00001 /* $Id: fors_img_science_impl.c,v 1.39 2008/08/04 12:05:54 cizzo Exp $
00002  *
00003  * This file is part of the FORS Library
00004  * Copyright (C) 2002-2006 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  02110-1301 USA
00019  */
00020 
00021 /*
00022  * $Author: cizzo $
00023  * $Date: 2008/08/04 12:05:54 $
00024  * $Revision: 1.39 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 #include <fors_img_science_impl.h>
00033 
00034 #include <fors_extract.h>
00035 #include <fors_tools.h>
00036 #include <fors_setting.h>
00037 #include <fors_data.h>
00038 #include <fors_image.h>
00039 #include <fors_qc.h>
00040 #include <fors_dfs.h>
00041 #include <fors_utils.h>
00042 
00043 #include <cpl.h>
00044 #include <math.h>
00045 #include <stdbool.h>
00046 
00053 const char *const fors_img_science_name = "fors_img_science";
00054 const char *const fors_img_science_description_short = "Reduce scientific exposure";
00055 const char *const fors_img_science_author = "Jonas M. Larsen";
00056 const char *const fors_img_science_email = PACKAGE_BUGREPORT;
00057 const char *const fors_img_science_description = 
00058 "Input files:\n"
00059 "  DO category:               Type:       Explanation:             Number:\n"
00060 "  SCIENCE_IMG                Raw         Science image               1\n"
00061 "  MASTER_BIAS                FITS image  Master bias                 1\n"
00062 "  MASTER_SKY_FLAT_IMG        FITS image  Master sky flat field       1\n"
00063 "\n"
00064 "Output files:\n"
00065 "  DO category:               Data type:  Explanation:\n"
00066 "  SCIENCE_REDUCED_IMG        FITS image  Reduced science image\n"
00067 "  PHOT_BACKGROUND_SCI_IMG    FITS image  Reduced science image background\n"
00068 "  SOURCES_SCI_IMG            FITS image  Unfiltered SExtractor output\n"
00069 "  OBJECT_TABLE_SCI_IMG       FITS table  Extracted sources properties\n";
00070 
00071 
00072 static double
00073 get_image_quality(const fors_star_list *sources, double *image_quality_err,
00074                   double *stellarity,
00075                   double *ellipticity,
00076                   double *ellipticity_rms);
00077 
00078 #undef cleanup
00079 #define cleanup \
00080 do { \
00081     cpl_free((void *)full_name); \
00082 } while (0)
00083 
00089 void fors_img_science_define_parameters(cpl_parameterlist *parameters)
00090 {
00091     cpl_parameter *p;
00092     const char *context = cpl_sprintf("fors.%s", fors_img_science_name);
00093     const char *full_name = NULL;
00094     const char *name;
00095 
00096     name = "cr_remove";
00097     full_name = cpl_sprintf("%s.%s", context, name);
00098     p = cpl_parameter_new_value(full_name,
00099                                 CPL_TYPE_BOOL,
00100                                 "Cosmic ray removal",
00101                                 context,
00102                                 false);
00103     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00104     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00105     cpl_parameterlist_append(parameters, p);
00106     cpl_free((void *)full_name); full_name = NULL;
00107 
00108     fors_extract_define_parameters(parameters, context);
00109     
00110     cpl_free((void *)context);
00111 
00112     return;
00113 }
00114 
00115 
00116 #undef cleanup
00117 #define cleanup \
00118 do { \
00119     cpl_frameset_delete(sci_frame); \
00120     cpl_frameset_delete(master_bias_frame); \
00121     cpl_frameset_delete(master_flat_frame); \
00122     fors_image_delete(&sci); \
00123     fors_image_delete_const(&master_bias); \
00124     fors_image_delete(&master_flat); \
00125     cpl_table_delete(phot); \
00126     cpl_table_delete(sources); \
00127     cpl_image_delete(background); \
00128     fors_extract_method_delete(&em); \
00129     fors_star_list_delete(&stars, fors_star_delete); \
00130     cpl_free((void *)context); \
00131     fors_setting_delete(&setting); \
00132     cpl_propertylist_delete(qc); \
00133     cpl_propertylist_delete(product_header); \
00134 } while (0)
00135 
00136 /* %%% Removed from cleanup
00137     cpl_frameset_delete(phot_table); \
00138 */
00139 
00146 void fors_img_science(cpl_frameset *frames, const cpl_parameterlist *parameters)
00147 {
00148     /* Raw */
00149     cpl_frameset *sci_frame      = NULL;
00150     fors_image *sci              = NULL;
00151 
00152     /* Calibration */
00153     cpl_frameset *master_bias_frame = NULL;
00154     const fors_image *master_bias   = NULL; 
00155 
00156     cpl_frameset *master_flat_frame = NULL;
00157     fors_image *master_flat         = NULL; 
00158 
00159 /* %%%
00160     cpl_frameset *phot_table        = NULL;
00161     double ext_coeff, dext_coeff;
00162 */
00163 
00164     /* Products */
00165     cpl_propertylist *qc = cpl_propertylist_new();
00166     cpl_propertylist *product_header = cpl_propertylist_new();
00167     cpl_table *phot = NULL;
00168     fors_extract_sky_stats sky_stats;
00169     cpl_image *background = NULL;
00170     cpl_table *sources = NULL;
00171 
00172     /* Parameters */
00173     extract_method  *em = NULL;
00174 
00175     /* Other */
00176     const char *context   = cpl_sprintf("fors.%s", fors_img_science_name);
00177     fors_star_list *stars = NULL;
00178     fors_setting *setting = NULL;
00179 
00180     /* Get parameters */
00181     em = fors_extract_method_new(parameters, context);
00182     assure( !cpl_error_get_code(), return, 
00183             "Could not get extraction parameters" );
00184     
00185     /* Find raw */
00186     sci_frame = fors_frameset_extract(frames, SCIENCE_IMG);
00187     assure( cpl_frameset_get_size(sci_frame) == 1, return, 
00188             "Exactly 1 %s required. %d found", 
00189             SCIENCE_IMG, cpl_frameset_get_size(sci_frame) );
00190 
00191     /* Find calibration */
00192     master_bias_frame = fors_frameset_extract(frames, MASTER_BIAS);
00193     assure( cpl_frameset_get_size(master_bias_frame) == 1, return, 
00194             "One %s required. %d found", 
00195             MASTER_BIAS, cpl_frameset_get_size(master_bias_frame) );
00196 
00197     master_flat_frame = fors_frameset_extract(frames, MASTER_SKY_FLAT_IMG);
00198     assure( cpl_frameset_get_size(master_flat_frame) == 1, return, 
00199             "One %s required. %d found", 
00200             MASTER_SKY_FLAT_IMG, cpl_frameset_get_size(master_flat_frame) );
00201 
00202 /* %%%
00203     phot_table = fors_frameset_extract(frames, PHOT_TABLE);
00204     assure( cpl_frameset_get_size(phot_table) == 1, return, 
00205             "One %s required. %d found",
00206             PHOT_TABLE, cpl_frameset_get_size(phot_table));
00207 */
00208 
00209     /* Done finding frames */
00210 
00211     /* Get instrument setting */
00212     setting = fors_setting_new(cpl_frameset_get_first(sci_frame));
00213     assure( !cpl_error_get_code(), return, "Could not get instrument setting" );
00214 
00215 
00216     master_bias = fors_image_load(cpl_frameset_get_first(master_bias_frame), 
00217                                   NULL, setting, NULL);
00218     assure( !cpl_error_get_code(), return, 
00219             "Could not load master bias");
00220 
00221     /* Load raw frames, subtract bias */
00222     sci = fors_image_load(cpl_frameset_get_first(sci_frame), master_bias, 
00223                           setting, NULL);
00224     assure( !cpl_error_get_code(), return, "Could not load standard image");
00225     fors_image_delete_const(&master_bias);
00226 
00227     /* Load master flat */
00228     master_flat = fors_image_load(cpl_frameset_get_first(master_flat_frame), 
00229                                   NULL, setting, NULL);
00230     assure( !cpl_error_get_code(), return, "Could not load master flat");
00231     
00232     /* Divide by normalized flat */
00233     fors_image_divide_scalar(master_flat,
00234                              fors_image_get_median(master_flat, NULL), -1.0);
00235 
00236     fors_image_divide(sci, master_flat);
00237     assure( !cpl_error_get_code(), return, "Could not divide by master flat");
00238     fors_image_delete(&master_flat);
00239 
00240     /* Extract sources */
00241     stars = fors_extract(sci, setting, em,
00242              &sky_stats, &background, &sources);
00243     assure( !cpl_error_get_code(), return, "Could not extract objects");  
00244 
00245     /* QC */
00246     fors_qc_start_group(qc, fors_qc_dic_version, setting->instrument);
00247     
00248     fors_qc_write_group_heading(cpl_frameset_get_first(sci_frame),
00249                                 PHOTOMETRY_TABLE,
00250                                 setting->instrument);
00251     assure( !cpl_error_get_code(), return, "Could not write %s QC parameters", 
00252             PHOTOMETRY_TABLE);
00253 
00254 
00255     double sky_mag;
00256     if (sky_stats.mean > 0) {
00257         sky_mag = -2.5*log(sky_stats.mean /
00258                            (setting->pixel_scale*setting->pixel_scale))/log(10);
00259     }
00260     else {
00261         cpl_msg_warning(cpl_func, 
00262                         "Average sky background is negative (%f ADU), "
00263                         "cannot compute magnitude, setting QC.SKYAVG to zero",
00264                         sky_stats.mean);
00265         sky_mag = 0;
00266     }
00267     fors_qc_write_qc_double(qc,
00268                             sky_mag,
00269                             "QC.SKYAVG",
00270                             "mag/arcsec^2",
00271                             "Mean of sky background",
00272                             setting->instrument);
00273     
00274     if (sky_stats.median > 0) {
00275         sky_mag = -2.5*log(sky_stats.median /
00276                            (setting->pixel_scale*setting->pixel_scale))/log(10);
00277     }
00278     else {
00279         cpl_msg_warning(cpl_func, 
00280                         "Median sky background is negative (%f ADU), "
00281                         "cannot compute magnitude, setting QC.SKYMED to zero",
00282                         sky_mag);
00283         sky_mag = 0;
00284     }
00285     fors_qc_write_qc_double(qc,
00286                             -2.5*log(sky_stats.median /
00287                                      (setting->pixel_scale*setting->pixel_scale)
00288                                 )/log(10),
00289                             "QC.SKYMED",
00290                             "mag/arcsec^2",
00291                             "Median of sky background",
00292                             setting->instrument);
00293 
00294     /* deltaM = -2.5*log10(e)*deltaF/F */
00295     fors_qc_write_qc_double(qc,
00296                             fabs(-2.5 * (1.0/log(10)) *
00297                                  sky_stats.rms/sky_stats.median),
00298                             "QC.SKYRMS",
00299                             "mag/arcsec^2",
00300                             "Standard deviation of sky background",
00301                             setting->instrument);
00302 
00303     double image_quality_error;
00304     double stellarity;
00305     double ellipticity, ellipticity_rms;
00306     double image_quality = get_image_quality(stars, 
00307                                              &image_quality_error,
00308                                              &stellarity,
00309                                              &ellipticity,
00310                                              &ellipticity_rms);
00311 
00312     if (image_quality > 0.) {
00313         image_quality *= TWOSQRT2LN2 * setting->pixel_scale;
00314         image_quality_error *= TWOSQRT2LN2 * setting->pixel_scale;
00315     }
00316 
00317     fors_qc_write_qc_double(qc,
00318                             image_quality,
00319                             "QC.IMGQU",
00320                             "arcsec",
00321                             "Image quality of scientific exposure",
00322                             setting->instrument);
00323 
00324     fors_qc_write_qc_double(qc,
00325                             image_quality_error,
00326                             "QC.IMGQUERR",
00327                             "arcsec",
00328                             "Uncertainty of image quality",
00329                             setting->instrument);
00330 
00331     fors_qc_write_qc_double(qc,
00332                             stellarity,
00333                             "QC.STELLAVG",
00334                             NULL,
00335                             "Mean stellarity index",
00336                             setting->instrument);
00337 
00338     fors_qc_write_qc_double(qc,
00339                             ellipticity,
00340                             "QC.IMGQUELL",
00341                             NULL,
00342                             "Mean star ellipticity",
00343                             setting->instrument);
00344 
00345     fors_qc_write_qc_double(qc,
00346                             ellipticity_rms,
00347                             "QC.IMGQUELLERR",
00348                             NULL,
00349                             "Standard deviation of star ellipticities",
00350                             setting->instrument);
00351 
00352     fors_qc_end_group();
00353 
00354     /* Save SCIENCE_REDUCED, PHOT_BACKGROUND_SCI_IMG */
00355 
00356     fors_dfs_add_wcs(qc, cpl_frameset_get_first(sci_frame), setting);
00357     fors_dfs_add_exptime(qc, cpl_frameset_get_first(sci_frame), 0.);
00358     fors_dfs_add_wcs(product_header, cpl_frameset_get_first(sci_frame), 
00359                      setting);
00360     fors_dfs_add_exptime(product_header, cpl_frameset_get_first(sci_frame), 0.);
00361 
00362     fors_dfs_save_image(frames, sci, SCIENCE_REDUCED_IMG,
00363                         qc, parameters, fors_img_science_name, 
00364                         cpl_frameset_get_first(sci_frame));
00365     assure( !cpl_error_get_code(), return, "Saving %s failed",
00366             SCIENCE_REDUCED_IMG);
00367 
00368     fors_image_delete(&sci);
00369     
00370     dfs_save_image(frames, background, PHOT_BACKGROUND_SCI_IMG,
00371                    product_header, parameters, fors_img_science_name, 
00372                    setting->version);
00373     assure( !cpl_error_get_code(), return, "Saving %s failed",
00374             PHOT_BACKGROUND_SCI_IMG);
00375 
00376     cpl_image_delete(background); background = NULL;
00377     
00378     /* Load filter coefficients */
00379 
00380 /* %%%
00381     fors_phot_table_load(cpl_frameset_get_first(phot_table), setting,
00382                          NULL, NULL, 
00383              &ext_coeff, &dext_coeff,
00384              NULL, NULL);
00385     assure( !cpl_error_get_code(), return, "Could not load photometry table" );
00386 */
00387 
00388     /* Correct for atmospheric extinction */
00389 /* %%%
00390     fors_star_ext_corr(stars, setting, ext_coeff, dext_coeff,
00391                        cpl_frameset_get_first(sci_frame));
00392     assure( !cpl_error_get_code(), return, 
00393             "Extinction correction failed");
00394 */
00395 
00396     /* Create, save FITS product */
00397     phot = fors_create_sources_table(stars);
00398     assure( !cpl_error_get_code(), return,
00399             "Failed to create extracted sources table");
00400 
00401     /*
00402      * Eliminate unused columns from photometry table
00403      */
00404 
00405     cpl_table_erase_column(phot, "INSTR_CMAG");
00406     cpl_table_erase_column(phot, "DINSTR_CMAG");
00407     cpl_table_erase_column(phot, "OBJECT");
00408     cpl_table_erase_column(phot, "MAG");
00409     cpl_table_erase_column(phot, "DMAG");
00410     cpl_table_erase_column(phot, "CAT_MAG");
00411     cpl_table_erase_column(phot, "DCAT_MAG");
00412     cpl_table_erase_column(phot, "COLOR");
00413     cpl_table_erase_column(phot, "USE_CAT");
00414     cpl_table_erase_column(phot, "SHIFT_X");
00415     cpl_table_erase_column(phot, "SHIFT_Y");
00416     cpl_table_erase_column(phot, "ZEROPOINT");
00417     cpl_table_erase_column(phot, "DZEROPOINT");
00418     cpl_table_erase_column(phot, "WEIGHT");
00419 
00420     fors_dfs_save_table(frames, sources, SOURCES_SCI,
00421                         NULL, parameters, fors_img_science_name, 
00422                         cpl_frameset_get_first(sci_frame));
00423     assure( !cpl_error_get_code(), return, "Saving %s failed",
00424             SOURCES_SCI);
00425 
00426     fors_dfs_save_table(frames, phot, PHOTOMETRY_TABLE,
00427                         NULL, parameters, fors_img_science_name, 
00428                         cpl_frameset_get_first(sci_frame));
00429     assure( !cpl_error_get_code(), return, "Saving %s failed",
00430             PHOTOMETRY_TABLE);
00431 
00432     cleanup;
00433     return;
00434 }
00435 
00436 #undef cleanup
00437 #define cleanup
00438 
00444 static bool
00445 is_star(const fors_star *s, void *data)
00446 {
00447     data = data;
00448     assure( s != NULL, return false, NULL );
00449 
00450 /*FIXME
00451   All stars for the moment... */
00452 
00453     return s->stellarity_index >= 0.0;
00454 }
00455 
00456 #undef cleanup
00457 #define cleanup \
00458 do { \
00459     fors_star_list_delete(&stars, fors_star_delete); \
00460 } while(0)
00461 
00474 static double
00475 get_image_quality(const fors_star_list *sources, double *image_quality_err,
00476                   double *stellarity,
00477                   double *ellipticity,
00478                   double *ellipticity_rms)
00479 {
00480     fors_star_list *stars = fors_star_list_extract(sources,
00481                                                    fors_star_duplicate,
00482                                                    is_star, NULL);
00483 
00484     double fwhm;
00485     if (fors_star_list_size(stars) >= 1) {
00486         *image_quality_err = fors_star_list_mad(stars, fors_star_extension, NULL) 
00487             * STDEV_PR_MAD;
00488         
00489         fwhm = fors_star_list_median(stars, fors_star_extension , NULL);
00490         
00491         *stellarity      = fors_star_list_mean(stars, fors_star_stellarity, NULL);
00492         *ellipticity     = fors_star_list_mean(stars, fors_star_ellipticity, NULL);
00493         *ellipticity_rms = fors_star_list_mad(stars, fors_star_ellipticity, NULL)
00494             * STDEV_PR_MAD;
00495     }
00496     else {
00497         cpl_msg_warning(cpl_func, "No stars found! Cannot compute image quality, "
00498                         "setting QC parameters to -1");
00499 
00500         /* -1 is not a valid value for any of these */
00501         *image_quality_err = -1;
00502         fwhm = -1;
00503         *stellarity = -1;
00504         *ellipticity = -1;
00505         *ellipticity_rms = -1;
00506     }
00507 
00508     cleanup;
00509     return fwhm;
00510 }
00511 

Generated on Wed Sep 10 07:31:51 2008 for FORS Pipeline Reference Manual by  doxygen 1.4.6