fors_img_science_impl.c

00001 /* $Id: fors_img_science_impl.c,v 1.42 2009/06/12 13:09:55 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: 2009/06/12 13:09:55 $
00024  * $Revision: 1.42 $
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     cpl_propertylist_delete(header); \
00135 } while (0)
00136 
00137 /* %%% Removed from cleanup
00138     cpl_frameset_delete(phot_table); \
00139 */
00140 
00147 void fors_img_science(cpl_frameset *frames, const cpl_parameterlist *parameters)
00148 {
00149     /* Raw */
00150     cpl_frameset *sci_frame      = NULL;
00151     fors_image *sci              = NULL;
00152 
00153     /* Calibration */
00154     cpl_frameset *master_bias_frame = NULL;
00155     const fors_image *master_bias   = NULL; 
00156 
00157     cpl_frameset *master_flat_frame = NULL;
00158     fors_image *master_flat         = NULL; 
00159 
00160 /* %%%
00161     cpl_frameset *phot_table        = NULL;
00162     double ext_coeff, dext_coeff;
00163 */
00164 
00165     /* Products */
00166     cpl_propertylist *qc = cpl_propertylist_new();
00167     cpl_propertylist *product_header = cpl_propertylist_new();
00168     cpl_propertylist *header = NULL;
00169     cpl_table *phot = NULL;
00170     fors_extract_sky_stats sky_stats;
00171     cpl_image *background = NULL;
00172     cpl_table *sources = NULL;
00173 
00174     /* Parameters */
00175     extract_method  *em = NULL;
00176 
00177     /* Other */
00178     const char *context   = cpl_sprintf("fors.%s", fors_img_science_name);
00179     fors_star_list *stars = NULL;
00180     fors_setting *setting = NULL;
00181     double avg_airmass = 0.0;
00182 
00183     /* Get parameters */
00184     em = fors_extract_method_new(parameters, context);
00185     assure( !cpl_error_get_code(), return, 
00186             "Could not get extraction parameters" );
00187     
00188     /* Find raw */
00189     sci_frame = fors_frameset_extract(frames, SCIENCE_IMG);
00190     assure( cpl_frameset_get_size(sci_frame) == 1, return, 
00191             "Exactly 1 %s required. %d found", 
00192             SCIENCE_IMG, cpl_frameset_get_size(sci_frame) );
00193 
00194     /* Find calibration */
00195     master_bias_frame = fors_frameset_extract(frames, MASTER_BIAS);
00196     assure( cpl_frameset_get_size(master_bias_frame) == 1, return, 
00197             "One %s required. %d found", 
00198             MASTER_BIAS, cpl_frameset_get_size(master_bias_frame) );
00199 
00200     master_flat_frame = fors_frameset_extract(frames, MASTER_SKY_FLAT_IMG);
00201     assure( cpl_frameset_get_size(master_flat_frame) == 1, return, 
00202             "One %s required. %d found", 
00203             MASTER_SKY_FLAT_IMG, cpl_frameset_get_size(master_flat_frame) );
00204 
00205 /* %%%
00206     phot_table = fors_frameset_extract(frames, PHOT_TABLE);
00207     assure( cpl_frameset_get_size(phot_table) == 1, return, 
00208             "One %s required. %d found",
00209             PHOT_TABLE, cpl_frameset_get_size(phot_table));
00210 */
00211 
00212     /* Done finding frames */
00213 
00214     /* Get instrument setting */
00215     setting = fors_setting_new(cpl_frameset_get_first(sci_frame));
00216     assure( !cpl_error_get_code(), return, "Could not get instrument setting" );
00217 
00218 
00219     master_bias = fors_image_load(cpl_frameset_get_first(master_bias_frame), 
00220                                   NULL, setting, NULL);
00221     assure( !cpl_error_get_code(), return, 
00222             "Could not load master bias");
00223 
00224     /* Load raw frames, subtract bias */
00225     sci = fors_image_load(cpl_frameset_get_first(sci_frame), master_bias, 
00226                           setting, NULL);
00227     assure( !cpl_error_get_code(), return, "Could not load standard image");
00228     fors_image_delete_const(&master_bias);
00229 
00230     /* Load master flat */
00231     master_flat = fors_image_load(cpl_frameset_get_first(master_flat_frame), 
00232                                   NULL, setting, NULL);
00233     assure( !cpl_error_get_code(), return, "Could not load master flat");
00234     
00235     /* Divide by normalized flat */
00236     fors_image_divide_scalar(master_flat,
00237                              fors_image_get_median(master_flat, NULL), -1.0);
00238 
00239     fors_image_divide(sci, master_flat);
00240     assure( !cpl_error_get_code(), return, "Could not divide by master flat");
00241     fors_image_delete(&master_flat);
00242 
00243     /* Extract sources */
00244     stars = fors_extract(sci, setting, em,
00245              &sky_stats, &background, &sources);
00246     assure( !cpl_error_get_code(), return, "Could not extract objects");  
00247 
00248     /* QC */
00249     fors_qc_start_group(qc, fors_qc_dic_version, setting->instrument);
00250     
00251     fors_qc_write_group_heading(cpl_frameset_get_first(sci_frame),
00252                                 PHOTOMETRY_TABLE,
00253                                 setting->instrument);
00254     assure( !cpl_error_get_code(), return, "Could not write %s QC parameters", 
00255             PHOTOMETRY_TABLE);
00256 
00257 
00258     double sky_mag;
00259     if (sky_stats.mean > 0) {
00260         sky_mag = -2.5*log(sky_stats.mean /
00261                            (setting->pixel_scale*setting->pixel_scale))/log(10);
00262     }
00263     else {
00264         cpl_msg_warning(cpl_func, 
00265                         "Average sky background is negative (%f ADU), "
00266                         "cannot compute magnitude, setting QC.SKYAVG to zero",
00267                         sky_stats.mean);
00268         sky_mag = 0;
00269     }
00270     fors_qc_write_qc_double(qc,
00271                             sky_mag,
00272                             "QC.SKYAVG",
00273                             "mag/arcsec^2",
00274                             "Mean of sky background",
00275                             setting->instrument);
00276     
00277     if (sky_stats.median > 0) {
00278         sky_mag = -2.5*log(sky_stats.median /
00279                            (setting->pixel_scale*setting->pixel_scale))/log(10);
00280     }
00281     else {
00282         cpl_msg_warning(cpl_func, 
00283                         "Median sky background is negative (%f ADU), "
00284                         "cannot compute magnitude, setting QC.SKYMED to zero",
00285                         sky_mag);
00286         sky_mag = 0;
00287     }
00288     fors_qc_write_qc_double(qc,
00289                             -2.5*log(sky_stats.median /
00290                                      (setting->pixel_scale*setting->pixel_scale)
00291                                 )/log(10),
00292                             "QC.SKYMED",
00293                             "mag/arcsec^2",
00294                             "Median of sky background",
00295                             setting->instrument);
00296 
00297     /* deltaM = -2.5*log10(e)*deltaF/F */
00298     fors_qc_write_qc_double(qc,
00299                             fabs(-2.5 * (1.0/log(10)) *
00300                                  sky_stats.rms/sky_stats.median),
00301                             "QC.SKYRMS",
00302                             "mag/arcsec^2",
00303                             "Standard deviation of sky background",
00304                             setting->instrument);
00305 
00306     double image_quality_error;
00307     double stellarity;
00308     double ellipticity, ellipticity_rms;
00309     double image_quality = get_image_quality(stars, 
00310                                              &image_quality_error,
00311                                              &stellarity,
00312                                              &ellipticity,
00313                                              &ellipticity_rms);
00314 
00315     if (image_quality > 0.) {
00316         image_quality *= TWOSQRT2LN2 * setting->pixel_scale;
00317         image_quality_error *= TWOSQRT2LN2 * setting->pixel_scale;
00318     }
00319 
00320     fors_qc_write_qc_double(qc,
00321                             image_quality,
00322                             "QC.IMGQU",
00323                             "arcsec",
00324                             "Image quality of scientific exposure",
00325                             setting->instrument);
00326 
00327     fors_qc_write_qc_double(qc,
00328                             image_quality_error,
00329                             "QC.IMGQUERR",
00330                             "arcsec",
00331                             "Uncertainty of image quality",
00332                             setting->instrument);
00333 
00334     fors_qc_write_qc_double(qc,
00335                             stellarity,
00336                             "QC.STELLAVG",
00337                             NULL,
00338                             "Mean stellarity index",
00339                             setting->instrument);
00340 
00341     fors_qc_write_qc_double(qc,
00342                             ellipticity,
00343                             "QC.IMGQUELL",
00344                             NULL,
00345                             "Mean star ellipticity",
00346                             setting->instrument);
00347 
00348     fors_qc_write_qc_double(qc,
00349                             ellipticity_rms,
00350                             "QC.IMGQUELLERR",
00351                             NULL,
00352                             "Standard deviation of star ellipticities",
00353                             setting->instrument);
00354 
00355     fors_qc_end_group();
00356 
00357     /* Save SCIENCE_REDUCED, PHOT_BACKGROUND_SCI_IMG */
00358 
00359 /* %%% */
00360 
00361     header = cpl_propertylist_load(
00362                         cpl_frame_get_filename(
00363                            cpl_frameset_get_first(sci_frame)), 0);
00364 
00365     if (header == NULL) {
00366         cpl_msg_error(cpl_func, "Failed to load raw header");
00367         cleanup;
00368         return;
00369     }
00370 
00371     avg_airmass = fors_get_airmass(header);
00372 
00373     cpl_propertylist_update_double(qc, "AIRMASS", avg_airmass);
00374     cpl_propertylist_update_double(product_header, "AIRMASS", avg_airmass);
00375 
00376 /* %%% */
00377 
00378     fors_dfs_add_wcs(qc, cpl_frameset_get_first(sci_frame), setting);
00379     fors_dfs_add_exptime(qc, cpl_frameset_get_first(sci_frame), 0.);
00380     fors_dfs_add_wcs(product_header, cpl_frameset_get_first(sci_frame), 
00381                      setting);
00382     fors_dfs_add_exptime(product_header, cpl_frameset_get_first(sci_frame), 0.);
00383 
00384     fors_dfs_save_image(frames, sci, SCIENCE_REDUCED_IMG,
00385                         qc, parameters, fors_img_science_name, 
00386                         cpl_frameset_get_first(sci_frame));
00387     assure( !cpl_error_get_code(), return, "Saving %s failed",
00388             SCIENCE_REDUCED_IMG);
00389 
00390     fors_image_delete(&sci);
00391     
00392     dfs_save_image(frames, background, PHOT_BACKGROUND_SCI_IMG,
00393                    product_header, parameters, fors_img_science_name, 
00394                    setting->version);
00395     assure( !cpl_error_get_code(), return, "Saving %s failed",
00396             PHOT_BACKGROUND_SCI_IMG);
00397 
00398     cpl_image_delete(background); background = NULL;
00399     
00400     /* Load filter coefficients */
00401 
00402 /* %%%
00403     fors_phot_table_load(cpl_frameset_get_first(phot_table), setting,
00404                          NULL, NULL, 
00405              &ext_coeff, &dext_coeff,
00406              NULL, NULL);
00407     assure( !cpl_error_get_code(), return, "Could not load photometry table" );
00408 */
00409 
00410     /* Correct for atmospheric extinction */
00411 /* %%%
00412     fors_star_ext_corr(stars, setting, ext_coeff, dext_coeff,
00413                        cpl_frameset_get_first(sci_frame));
00414     assure( !cpl_error_get_code(), return, 
00415             "Extinction correction failed");
00416 */
00417 
00418     /* Create, save FITS product */
00419     phot = fors_create_sources_table(stars);
00420     assure( !cpl_error_get_code(), return,
00421             "Failed to create extracted sources table");
00422 
00423     /*
00424      * Eliminate unused columns from photometry table
00425      */
00426 
00427     cpl_table_erase_column(phot, "INSTR_CMAG");
00428     cpl_table_erase_column(phot, "DINSTR_CMAG");
00429     cpl_table_erase_column(phot, "OBJECT");
00430     cpl_table_erase_column(phot, "MAG");
00431     cpl_table_erase_column(phot, "DMAG");
00432     cpl_table_erase_column(phot, "CAT_MAG");
00433     cpl_table_erase_column(phot, "DCAT_MAG");
00434     cpl_table_erase_column(phot, "COLOR");
00435     /* new columns since 4.4.10 */
00436     if (cpl_table_has_column(phot, "DCOLOR"))
00437         cpl_table_erase_column(phot, "DCOLOR");
00438     if (cpl_table_has_column(phot, "COV_CATM_COL"))
00439         cpl_table_erase_column(phot, "COV_CATM_COL");
00440     cpl_table_erase_column(phot, "USE_CAT");
00441     cpl_table_erase_column(phot, "SHIFT_X");
00442     cpl_table_erase_column(phot, "SHIFT_Y");
00443     cpl_table_erase_column(phot, "ZEROPOINT");
00444     cpl_table_erase_column(phot, "DZEROPOINT");
00445     cpl_table_erase_column(phot, "WEIGHT");
00446 
00447     fors_dfs_save_table(frames, sources, SOURCES_SCI,
00448                         NULL, parameters, fors_img_science_name, 
00449                         cpl_frameset_get_first(sci_frame));
00450     assure( !cpl_error_get_code(), return, "Saving %s failed",
00451             SOURCES_SCI);
00452 
00453     fors_dfs_save_table(frames, phot, PHOTOMETRY_TABLE,
00454                         NULL, parameters, fors_img_science_name, 
00455                         cpl_frameset_get_first(sci_frame));
00456     assure( !cpl_error_get_code(), return, "Saving %s failed",
00457             PHOTOMETRY_TABLE);
00458 
00459     cleanup;
00460     return;
00461 }
00462 
00463 #undef cleanup
00464 #define cleanup
00465 
00471 static bool
00472 is_star(const fors_star *s, void *data)
00473 {
00474     data = data;
00475     assure( s != NULL, return false, NULL );
00476 
00477 /*FIXME
00478   All stars for the moment... */
00479 
00480     return s->stellarity_index >= 0.7;
00481 }
00482 
00483 #undef cleanup
00484 #define cleanup \
00485 do { \
00486     fors_star_list_delete(&stars, fors_star_delete); \
00487 } while(0)
00488 
00501 static double
00502 get_image_quality(const fors_star_list *sources, double *image_quality_err,
00503                   double *stellarity,
00504                   double *ellipticity,
00505                   double *ellipticity_rms)
00506 {
00507     fors_star_list *stars = fors_star_list_extract(sources,
00508                                                    fors_star_duplicate,
00509                                                    is_star, NULL);
00510 
00511     double fwhm;
00512     if (fors_star_list_size(stars) >= 1) {
00513         *image_quality_err = fors_star_list_mad(stars, fors_star_extension, NULL) 
00514             * STDEV_PR_MAD;
00515         
00516         fwhm = fors_star_list_median(stars, fors_star_extension , NULL);
00517         
00518         *stellarity      = fors_star_list_mean(stars, fors_star_stellarity, NULL);
00519         *ellipticity     = fors_star_list_mean(stars, fors_star_ellipticity, NULL);
00520         *ellipticity_rms = fors_star_list_mad(stars, fors_star_ellipticity, NULL)
00521             * STDEV_PR_MAD;
00522     }
00523     else {
00524         cpl_msg_warning(cpl_func, "No stars found! Cannot compute image quality, "
00525                         "setting QC parameters to -1");
00526 
00527         /* -1 is not a valid value for any of these */
00528         *image_quality_err = -1;
00529         fwhm = -1;
00530         *stellarity = -1;
00531         *ellipticity = -1;
00532         *ellipticity_rms = -1;
00533     }
00534 
00535     cleanup;
00536     return fwhm;
00537 }
00538 

Generated on Wed Jun 24 14:11:17 2009 for FORS Pipeline Reference Manual by  doxygen 1.4.7