/* * This file is part of the ESO SINFONI Pipeline * Copyright (C) 2004,2005 European Southern Observatory * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, 51 Franklin St, Fifth Floor, Boston, MA 02111-1307 USA */ /*---------------------------------------------------------------------------- File name : sinfo_bp_norm.c Author : J. Schreiber Created on : May 5, 2003 Description : Different methods for searching for bad pixels used in the recipe sinfo_rec_mflat ---------------------------------------------------------------------------*/ #ifdef HAVE_CONFIG_H # include #endif /*---------------------------------------------------------------------------- Includes ---------------------------------------------------------------------------*/ #include "sinfo_bp_norm.h" #include "sinfo_image_ops.h" #include "sinfo_detlin.h" #include "sinfo_badnorm_ini_by_cpl.h" #include "sinfo_baddist_ini_by_cpl.h" #include "sinfo_pro_save.h" #include "sinfo_functions.h" #include "sinfo_pro_types.h" #include "sinfo_hidden.h" #include "sinfo_error.h" #include "sinfo_utils_wrappers.h" /**@{*/ /** * @addtogroup sinfo_bad_pix_search Bad Pixel Search * * TBD */ /*---------------------------------------------------------------------------- Defines ---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------- Function Definitions ---------------------------------------------------------------------------*/ /** @name sinfo_new_bp_search_norm @memo finds pixels whose intensity differs from the median more than a given threshold @param plugin_id recipe name @param config input parameters list @param sof input set of frames @param procatg PRO.CATG of product. @return integer (0 if it worked, -1 if it doesn't) @doc this function searches for static bad pixels in stacks of flatfield frames, that means it stacks the flatfield in a data cube and takes the sinfo_median along the z-axis to be sure to have no cosmics, then the intensity tilt of each column is removed. The standard deviation and mean of the pixel values within a defined rectangular zone is determined in a way that the extreme low and high values are cut off. The next step is to indicate each pixel as bad that has a deviation from the mean greater than a user defined factor times the standard deviation. This step can be leaped over if wished. The next step is to calculate the median of the 8 nearest neighbors for each pixel and to determine the deviation of each pixel value from this sinfo_median. If this deviation is greater than a defined factor times the standard deviation the pixel value is replaced by the sinfo_median. The last step is repeated to be able to consider also clusters of bad pixels. Then the resulting image is compared with the input image with the column intensity tilt removed and each changed pixel is indicated as bad. Finally, a bad pixel mask is produced that means each good pixel is marked with 1 and each bad pixel with 0. */ int sinfo_new_bp_search_normal(const char* plugin_id, cpl_parameterlist* config, cpl_frameset* sof, cpl_frameset* ref_set, const char* procatg) { bad_config * cfg = NULL; cpl_imagelist * image_list = NULL; cpl_image ** med = NULL; cpl_image * medImage = NULL; cpl_image * medIm = NULL; cpl_image * colImage = NULL; cpl_image * compImage = NULL; cpl_image * maskImage = NULL; cpl_image * threshIm = NULL; Stats * stats = NULL; cpl_parameter *p = NULL; int no = 0; float lo_cut = 0.; float hi_cut = 0.; int i = 0; int n = 0; int half_box_size = 0; cpl_frameset* raw = NULL; cpl_table* qclog_tbl = NULL; char key_value[FILE_NAME_SZ]; /* parse the file names and parameters to the bad_config data structure cfg */ sinfo_check_rec_status(0); check_nomsg(raw = cpl_frameset_new()); sinfo_check_rec_status(1); if (strcmp(procatg, PRO_BP_MAP_NO) == 0) { cknull(cfg = sinfo_parse_cpl_input_badnorm(config, sof, procatg, &raw), "could not parse cpl input!"); } else if (strcmp(procatg, PRO_BP_MAP_DI) == 0) { cknull(cfg = sinfo_parse_cpl_input_baddist(config, sof, procatg, &raw), "could not parse cpl input!"); } else if (strcmp(procatg, PRO_DEFAULT) == 0) { cknull(cfg = sinfo_parse_cpl_input_badnorm(config, sof, procatg, &raw), "could not parse cpl input!"); } else { sinfo_msg_error("Error: PRO.CATG %s, not supported!", procatg); goto cleanup; } sinfo_check_rec_status(2); /* take a clean mean of the frames */ sinfo_msg("Takes a clean mean of the frames"); check_nomsg(image_list = cpl_imagelist_new()); sinfo_check_rec_status(3); for (i = 0; i < cfg->nframes; i++) { if (sinfo_is_fits_file(cfg->framelist[i]) != 1) { sinfo_msg_error("Input file %s is not FITS", cfg->framelist[i]); goto cleanup; } check_nomsg( cpl_imagelist_set(image_list, cpl_image_load(cfg->framelist[i], CPL_TYPE_FLOAT, 0, 0), i)); } /* finally take the average image of the cube by rejecting the extreme values */ check_nomsg(no = cpl_imagelist_get_size(image_list)); lo_cut = (floor)(cfg->loReject * no + 0.5); hi_cut = (floor)(cfg->hiReject * no + 0.5); cknull(medImage=cpl_imagelist_collapse_minmax_create(image_list, lo_cut, hi_cut), "error in sinfo_average_with_rejection"); /* free memory */ sinfo_free_imagelist(&image_list); /*---------------------------------------------- * remove the intensity tilt from every column * and compute the standard deviation on a rectangular zone */ cknull(medIm = sinfo_new_thresh_image(medImage, cfg->mincut, cfg->maxcut), "error sinfo_new_thresh_image"); cknull(colImage = sinfo_new_col_tilt( medIm, cfg->sigmaFactor ), "sinfo_colTilt failed" ); cknull(stats = sinfo_new_image_stats_on_rectangle(colImage, cfg->loReject, cfg->hiReject, cfg->llx, cfg->lly, cfg->urx, cfg->ury), " sinfo_get_image_stats_on_vig failed"); if (stats != NULL ) { sinfo_msg("Clean stdev: %f\n", stats->cleanstdev ); sinfo_msg("Clean mean: %f\n", stats->cleanmean ); } /* cknull_nomsg(qclog_tbl = sinfo_qclog_init()); ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC ICTILT STDEV", stats->cleanstdev, "Intensity column clean stdev","%g")); ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC ICTILT MEAN", stats->cleanmean, "Intensity column clean mean","%g")); ck0(sinfo_pro_save_ima(colImage,raw,sof, (char*) BP_NORM_INT_COL_TILT_CORR_OUT_FILENAME, PRO_INT_COL_TILT_COR,qclog_tbl,plugin_id,config), "cannot save ima %s", BP_NORM_INT_COL_TILT_CORR_OUT_FILENAME); sinfo_free_table(&qclog_tbl); */ /* indicate pixels with great deviations from the clean mean as bad */ if (cfg->threshInd == 1) { cknull(threshIm = sinfo_new_thresh_image(colImage, stats->cleanmean-cfg->meanfactor*stats->cleanstdev, stats->cleanmean+cfg->meanfactor*stats->cleanstdev), " sinfo_threshImage failed" ); } if (cfg->threshInd == 0) { threshIm = colImage; } /* AMO here invalid fread? */ med = (cpl_image**) cpl_calloc(cfg->iterations, sizeof(cpl_image*)); /* filter iteratively the images by a sinfo_median filter of the nearest neighbors under the condition of a deviation greater than a factor times the standard deviation */ sinfo_msg("Apply sinfo_median filter on pixel nearest neighbors"); if (cfg->methodInd == 1) { if (cfg->factor > 0) { cknull(med[0]=sinfo_new_median_image(threshIm, -cfg->factor*stats->cleanstdev), " sinfo_medianImage failed (1)" ); for (i = 0; i < cfg->iterations - 1; i++) { cknull(med[i + 1] = sinfo_new_median_image(med[i], -cfg->factor* stats->cleanstdev), "sinfo_medianImage failed (2)"); } } else { cknull(med[0] = sinfo_new_median_image(threshIm, -cfg->factor), " sinfo_medianImage failed (1)" ); for (i = 0; i < cfg->iterations - 1; i++) { cknull(med[i+1] = sinfo_new_median_image(med[i], -cfg->factor), " sinfo_medianImage failed (%d)",i ); } } } else if (cfg->methodInd == 2) { cknull(med[0] = sinfo_new_abs_dist_image(threshIm, -cfg->factor), " sinfo_absDistImage failed (1)" ); for (i = 0; i < cfg->iterations - 1; i++) { cknull(med[i+1] = sinfo_new_abs_dist_image(med[i], -cfg->factor), " sinfo_absDistImage failed (2)" ); } } else if (cfg->methodInd == 3) { cknull( med[0] = sinfo_new_mean_image_in_spec(threshIm, -cfg->factor * stats->cleanstdev), "sinfo_meanImageInSpec failed (1)"); for (i = 0; i < cfg->iterations - 1; i++) { cknull( med[i + 1] = sinfo_new_mean_image_in_spec(med[i], -cfg->factor * stats->cleanstdev), " sinfo_meanImageInSpec failed (2)"); } } else if (cfg->methodInd == 4) { half_box_size = (cfg->urx - cfg->llx) / 2; cknull(med[0] = sinfo_new_local_median_image(threshIm, -cfg->factor, cfg->loReject, cfg->hiReject, half_box_size), " sinfo_localMedianImage failed (1)" ); for (i = 0; i < cfg->iterations - 1; i++) { cknull(med[i+1] = sinfo_new_local_median_image(med[i], -cfg->factor, cfg->loReject, cfg->hiReject, half_box_size), " sinfo_localMedianImage failed (2)" ); } } else { sinfo_msg_error (" wrong indicator methodInd !" ); goto cleanup; } /* compare the filtered image with the input image */ cknull(compImage = sinfo_new_compare_images(threshIm, med[cfg->iterations - 1], medImage), " sinfo_compareImages failed" ); /* generate the bad pixel mask */ sinfo_msg("Generates bad pixel map"); cknull(maskImage = sinfo_new_promote_image_to_mask( compImage, &n ), " error in sinfo_promoteImageToMask" ); sinfo_msg("No of bad pixels: %d\n", n ); cknull_nomsg(qclog_tbl = sinfo_qclog_init()); check_nomsg(p = cpl_parameterlist_find(config, "sinfoni.bp.method")); snprintf(key_value, MAX_NAME_SIZE - 1, "%s", cpl_parameter_get_string(p)); ck0_nomsg( sinfo_qclog_add_string(qclog_tbl, "QC BP-MAP METHOD", key_value, "BP search method", "%s")); ck0_nomsg( sinfo_qclog_add_int(qclog_tbl, "QC BP-MAP NBADPIX", n, "No of bad pixels", "%d")); ck0( sinfo_pro_save_ima(maskImage, ref_set, sof, cfg->outName, procatg, qclog_tbl, plugin_id, config), "cannot save ima %s", cfg->outName); sinfo_free_table(&qclog_tbl); sinfo_free_image(&maskImage); sinfo_free_image(&compImage); sinfo_free_image_array(&med, cfg->iterations); /* if (med != NULL) { for ( i = 0 ; i < cfg->iterations ; i++ ) { if(med[i] != NULL) { cpl_image_delete(med[i]) ; med[i]=NULL; } } cpl_free(med) ; med=NULL; } */ if (stats != NULL ) { cpl_free(stats); stats = NULL; } sinfo_free_image(&medIm); sinfo_free_image(&medImage); sinfo_free_image(&colImage); if (cfg->threshInd == 1) { sinfo_free_image(&threshIm); } sinfo_badnorm_free(&cfg); sinfo_free_frameset(&raw); return 0; cleanup: sinfo_free_image_array(&med, cfg->iterations); /* if (med != NULL) { for ( i = 0 ; i < cfg->iterations ; i++ ) { if(med[i] != NULL) { cpl_image_delete(med[i]) ; med[i]=NULL; } } cpl_free(med) ; med=NULL; } */ sinfo_free_image(&compImage); sinfo_free_image(&maskImage); sinfo_free_image(&threshIm); sinfo_free_table(&qclog_tbl); sinfo_free_image(&threshIm); if (stats != NULL ) { cpl_free(stats); stats = NULL; } sinfo_free_image(&medIm); sinfo_free_image(&medImage); sinfo_free_image(&colImage); sinfo_free_imagelist(&image_list); sinfo_free_frameset(&raw); sinfo_badnorm_free(&cfg); return -1; } /**@}*/