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 #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
00137
00138
00139
00146 void fors_img_science(cpl_frameset *frames, const cpl_parameterlist *parameters)
00147 {
00148
00149 cpl_frameset *sci_frame = NULL;
00150 fors_image *sci = NULL;
00151
00152
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
00161
00162
00163
00164
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
00173 extract_method *em = NULL;
00174
00175
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
00181 em = fors_extract_method_new(parameters, context);
00182 assure( !cpl_error_get_code(), return,
00183 "Could not get extraction parameters" );
00184
00185
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
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
00204
00205
00206
00207
00208
00209
00210
00211
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
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
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
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
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
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
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
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
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397 phot = fors_create_sources_table(stars);
00398 assure( !cpl_error_get_code(), return,
00399 "Failed to create extracted sources table");
00400
00401
00402
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
00451
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
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