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 cpl_propertylist_delete(header); \
00135 } while (0)
00136
00137
00138
00139
00140
00147 void fors_img_science(cpl_frameset *frames, const cpl_parameterlist *parameters)
00148 {
00149
00150 cpl_frameset *sci_frame = NULL;
00151 fors_image *sci = NULL;
00152
00153
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
00162
00163
00164
00165
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
00175 extract_method *em = NULL;
00176
00177
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
00184 em = fors_extract_method_new(parameters, context);
00185 assure( !cpl_error_get_code(), return,
00186 "Could not get extraction parameters" );
00187
00188
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
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
00207
00208
00209
00210
00211
00212
00213
00214
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
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
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
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
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
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
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
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
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419 phot = fors_create_sources_table(stars);
00420 assure( !cpl_error_get_code(), return,
00421 "Failed to create extracted sources table");
00422
00423
00424
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
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
00478
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
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