#include "i3_mcmc.h" #include "i3_load_data.h" #include "i3_math.h" #include "i3_model.h" #include "i3_minimizer.h" #include "i3_fisher.h" #include "i3_mcmc.h" #include "i3_utils.h" #include "i3_great.h" #include "i3_image_fits.h" #include "i3_analyze.h" #include "i3_version.h" #include "i3_logging.h" #include int main_analyze(int argc, char * argv[]); int main_generate(int argc, char * argv[]); int parse_generate_parameters(int argc, char * argv[], int first_arg, i3_model * model, i3_parameter_set * params, i3_moffat_psf * psf, i3_flt *noise_sigma); void usage_message_and_exit() { fprintf(stderr, "Standard usage:\n"); fprintf(stderr, "parameter_filename image_file object_catalogue_file psf_file output_filename \ [first_image_to_process] [last_image_to_process] [additional ini file options]\n"); fprintf(stderr, "\nFor help on usage to generate test-images run im3shape --generate\n\n"); exit(1); } int main(int argc, char * argv[]){ // Totalhis driver now gives very similar output to the great08 driver on LowNoise_Known - except for x, y flip if (argc==1) usage_message_and_exit(); if (argc==2 && (strcmp(argv[1], "-v")==0 || strcmp(argv[1],"--version")==0 )){ printf("im3shape-version: %s\n",HG_VERSION_STRING); exit(0); } if (0==strcmp(argv[1], "--generate") || 0==strcmp(argv[1],"-g")){ return main_generate(argc, argv); } return main_analyze(argc, argv); } int main_analyze(int argc, char * argv[]) { int n_args=6; if(argcverbosity); // If args 6 and 7 are numbers then they must be image ranges, so use them. //Otherwise they must be command line options, so read them. if (argc>=7){ if (isdigit(argv[6][0])) { first_image = atoi(argv[6]); if (argc>=8){ if (isdigit(argv[7][0])) { last_image = atoi(argv[7]); if (argc>=9) i3_options_read_command_line(options, argc, 8, argv); } else i3_options_read_command_line(options, argc, 7, argv); } } else i3_options_read_command_line(options, argc, 6, argv); } //Read in the full image. // i3_image * full_image = i3_read_fits_image(input_image_filename); // if (!full_image) I3_FATAL("Could not load image from fits file.",1); int status=0; fitsfile * image_file = NULL; fits_open_file(&image_file, input_image_filename, READONLY, &status); long full_image_dimensions[2]; fits_get_img_size(image_file, 2, full_image_dimensions, &status); int full_image_nx = full_image_dimensions[0]; int full_image_ny = full_image_dimensions[1]; if (status) I3_FATAL("Could not open fits file", 1); //Load the object catalogue i3_catalog * galaxy_cat = i3_catalog_read(input_obj_filename); // Currently just quits silently if no galaxy catalogue file is found - ideally put error message and quit //Read in the segmentation mask image if any should be used. fitsfile * segmentation_mask_file = NULL; if (options->use_segmentation_mask){ fits_open_file(&segmentation_mask_file, options->segmentation_mask_filename, READONLY, &status); if (status) I3_FATAL("Could not load segmentation mask image from fits file.",1); fits_get_img_size(image_file, 2, full_image_dimensions, &status); if ((full_image_dimensions[0]!=full_image_nx) || (full_image_dimensions[1]!=full_image_ny)){ I3_FATAL("Segementation map different size to image", 4); } } // Identify the type of the psf and create it // In general we should do something smarter with switching the PSF catalog type. For now let's just define both pointers, one of them won't be used i3_moffat_psf * cat_psf_moffat; int psf_type; // only needed when we are using moffat_catalog i3_flt truncation = options->psf_truncation_pixels; // only needed when we are using moffat_catalog i3_moffat_psf * psf_func; // pointer for current psf function i3_image * psf_image; // pointer for current psf image int npsf; if(!strcmp(options->psf_input,"psf_image_cube")){ // check the number of PSFs in the cube npsf = i3_fits_cube_number_images(input_psf_filename, 1); i3_print(i3_verb_standard,"using %s %s with %d psfs \n",input_psf_filename,options->psf_input,npsf); // read the first psf image from the cube int tmp = 0; i3_image ** temp_psf_image_ptr = i3_read_fits_cube(input_psf_filename, 1, 0, 1, &tmp); if (temp_psf_image_ptr==NULL){ I3_FATAL("Unable to open PSF cube",tmp); } i3_image * temp_psf_image = temp_psf_image_ptr[0]; // check the size of the images in the cube int n_hires_pixels = (options->stamp_size + options->padding)*options->upsampling; int n_psf_pixels = temp_psf_image->nx; if(n_psf_pixels != n_hires_pixels){ printf("psf image size = %d\nhires image size =%d\n",n_psf_pixels,n_hires_pixels); I3_FATAL("provided PSF image resolution doesn't match the settings provided in options file",2); } i3_image_destroy( temp_psf_image ); }else if(!strcmp(options->psf_input,"moffat_catalog")){ cat_psf_moffat = i3_psf_catalog_read(input_psf_filename, &npsf); // assign first from the catalog psf_func = cat_psf_moffat; //Select the type of PSF model from the options file. if (options->airy_psf) { psf_type = GREAT10_PSF_TYPE_AIRY; // Ideally would remove the "GREAT10" part of this name } else { psf_type = GREAT10_PSF_TYPE_MOFFAT; // Ditto re GREAT10 } i3_print(i3_verb_standard,"using %s with %d psfs \n",options->psf_input,npsf); //printf("%f %f\n",cat_psf_moffat->beta,cat_psf_moffat->e1); }else if(!strcmp(options->psf_input,"psf_image_single")){ // read the first psf image i3_print(i3_verb_standard,"using single psf for all galaxies: %s \n",options->psf_input); psf_image = i3_read_fits_image(input_psf_filename); // check the size of the images in the cube int n_hires_pixels = (options->stamp_size + options->padding)*options->upsampling; int n_psf_pixels = psf_image->nx; if(n_psf_pixels != n_hires_pixels){ printf("psf image size = %d\nhires image size =%d\n",n_psf_pixels,n_hires_pixels); I3_FATAL("provided PSF image resolution doesn't match the settings provided in options file",3); } //i3_image_destroy(this_psf_image); }else I3_FATAL("Unknown PSF input method",4); //Create basic information and space for the chosen model (e.g. sersics) i3_model * model = i3_model_create(options->model_name,options); if (model==NULL) I3_FATAL("Unknown model type specified in parameter file. Either typo or model not compiled in Makefile",5); //Make a dummy weight image. Assume noise level is constant across image bool individual_weights = false; i3_image * weight_stamp = i3_image_create(options->stamp_size,options->stamp_size); if(options->noise_sigma>0) i3_image_fill(weight_stamp,pow(options->noise_sigma,-2)); else i3_image_fill(weight_stamp,1.); fitsfile * image_mask_file=NULL; // If specified, we have a mask which removes image sections. // Load it as a full-size image if (strcmp(options->image_mask, "")){ individual_weights = true; fits_open_file(&image_mask_file, options->image_mask, READONLY, &status); if (status) I3_FATAL("Could not read image mask from specified mask file",5); fits_get_img_size(image_file, 2, full_image_dimensions, &status); if ((full_image_dimensions[0]!=full_image_nx) || (full_image_dimensions[1]!=full_image_ny)){ I3_FATAL("Segementation map different size to image", 4); } } else{ } //Set the minimizer method as levmar i3_minimizer_method method; method = (i3_minimizer_method) options->minimizer_method; if (method != i3_minimizer_method_levmar){ I3_WARNING("Note: Not using LevMar - other methods not well tested."); } // is this an OK place to put this?? I guess we don't want to set in ini file? //Set up space for the best likelihood, best parameters and best image i3_flt bestL = 0; i3_parameter_set * params_ml; // help, this is specific to sersics!!! i3_image * bestI = i3_image_create(options->stamp_size,options->stamp_size); //Prepare the output file FILE * output = fopen(output_filename,"w"); if(!output) I3_FATAL("Cannot open output file for writing",7); if(!strcmp(model->name,"sersics") || !strcmp(model->name,"logsersics")) fprintf(output,"%s",i3_analyze_sersics_header_line()); // this is not true for now // fprintf(output,"# ID like %s flux_ratio S/N min_residuals max_residuals ....lots of other junk...\n",i3_model_header_line(model)); //The concept of a different image and stamp size is now obselete (given that we have the catalogue) int image_size = options->stamp_size; time_t time_start, time_end; i3_flt cpu_time_used; int gals_counter = 0; if(options->timeit){ time_start = time(NULL); }; i3_print(i3_verb_quiet, "Analyzing %d galaxies", galaxy_cat->n); // Loop over rows in the galaxy catalogue file for (int i=0;in;i++){ // Get the postage stamp to analyze if (ilast_image) continue; // Count number of processed galaxies gals_counter += 1; clock_t start_time = clock(); i3_catalog_row * row = galaxy_cat->row + i; gal_id identifier = row->id; int ix = (int) (row->x - image_size/2.0); int iy = (int) (row->y - image_size/2.0); //printf("row = %d;\n",row->catalog_row); //printf("id = %d; x = %d; y = %d\n",identifier,ix,iy); // Need to catch objects with ix, iy <0 and ix+image_size > fits file. // This is a very naff hack by Sarah - just set image to zero if roughly off cut-out (to 1 or 2 pixels accuracy) //JAZ The i3_image_copy_part function already does a check for images off the edge // of the full image, and returns NULL if so. So we can just rely on that. // i3_image *galaxy = i3_image_copy_part(full_image,ix,iy,image_size,image_size); i3_image * galaxy = i3_fits_image_copy_part(image_file, ix, iy, image_size, image_size); if (galaxy==NULL){ i3_print(i3_verb_quiet,"Skipping galaxy %ld since off edge: ix=%d iy=%d size=%d, fits_image_size=%d,%d",identifier,ix,iy,image_size,full_image_nx,full_image_ny); continue; } // Load segmentation mask which is needed for background estimation i3_image * mask_stamp=NULL; if (options->use_segmentation_mask){ // Extract the part of the mask to be used for this galaxy i3_image * mask_chunk = i3_fits_image_copy_part(segmentation_mask_file,ix,iy,image_size,image_size); // Now cut the stamp out of the chunk. if (mask_chunk) mask_stamp=i3_image_cut_out_stamp(mask_chunk,options->stamp_size); i3_image_destroy(mask_chunk); if (mask_stamp==NULL){ i3_print(i3_verb_quiet,"Failed to cut out segmentation mask %ld stamp at: ix=%d iy=%d size=%d, fits_image_size=%d,%d",identifier,ix,iy,image_size,full_image_nx,full_image_ny); i3_image_destroy(mask_stamp); i3_image_destroy(galaxy); continue; } } // Arghh, we need to fit a background level for ClusterSTEP, DES etc // Bodge this in really messily for now.... Sorry... SLB. // This is not a bad fix for ClusterSTEP but will bias GREAT runs which have no background // Ideally the options file would include a switch for whether or not we estimate the background // This "background_estimate" step should go into a function really... if (options->background_subtract) { i3_flt background_estimate = 0.0; if (options->use_segmentation_mask) I3_FATAL("Code missing for seg mask and bg subtraction in C",7); else background_estimate = i3_image_subtract_background(galaxy); i3_print(i3_verb_noisy,"background_estimate = %20.10f",background_estimate); } // Not convinced the below line is currently doing this (seems to run with ix<0) if (!galaxy) I3_FATAL("Failed to get galaxy from full image",8); i3_image * stamp = i3_image_cut_out_stamp(galaxy,options->stamp_size); if (!stamp) I3_FATAL("Failed to cut stamp out from galaxy",9); //Scale the image and noise to unit maximum if (options->rescale_stamp) { i3_flt scaling = 1.0/i3_image_sum(stamp); i3_image_scale(stamp,scaling); if(options->noise_sigma>0) i3_image_fill(weight_stamp,pow(options->noise_sigma*scaling,-2)); } // If there is a weighting to apply, do that now. if (individual_weights){ i3_image *mask_chunk = i3_fits_image_copy_part(image_mask_file,ix,iy,image_size,image_size); i3_image *mask_stamp2 = i3_image_cut_out_stamp(mask_chunk,options->stamp_size); i3_image_multiply_mask(weight_stamp, mask_stamp); i3_print(i3_verb_debug, "Stamp mask fraction = %lf\n", 1.0 - i3_image_sum(mask_stamp)/mask_stamp->n); i3_image_destroy(mask_chunk); i3_image_destroy(mask_stamp2); } // Modify weighting map by setting those pixel to zero that are assigned // to objects other than the one currently processed if (options->use_segmentation_mask){ i3_modify_weight_map_by_segmentation_mask(identifier, weight_stamp, mask_stamp); // If asked, save the mask stamp here so that we do not need to save if for later. if (options->save_images){ char filename[256]; snprintf(filename,256,"%s/mask.%ld.fits",options->output_directory,identifier); i3_image_save_fits(mask_stamp,filename); } i3_image_destroy(mask_stamp); } //Get the relevant PSF parameters for this postage stamp // Again, PSF switching is not great here, we should be doing something smarter i3_data_set * data_set; if(!strcmp(options->psf_input,"psf_image_cube")){ //printf("npsf = %d\n",npsf); int tmp = 0; i3_image ** temp_psf_image_ptr; if(npsf > 1) temp_psf_image_ptr = i3_read_fits_cube(input_psf_filename, 1, i,1, &tmp); //printf("n_pad = %d\nn_sub =%d\nn_pix =%d\n",options->padding,options->stamp_size,options->upsampling); else temp_psf_image_ptr = i3_read_fits_cube(input_psf_filename, 1, 0,1, &tmp); if (temp_psf_image_ptr==NULL){ I3_FATAL("Not enough images in PSF cube",i+1); } i3_image * temp_psf_image = temp_psf_image_ptr[0]; data_set = i3_build_dataset_with_psf_image( options, identifier, stamp, weight_stamp, temp_psf_image ); i3_image_destroy(temp_psf_image); //printf("nx = %d\n",this_psf_image->nx); } if(!strcmp(options->psf_input,"psf_image_single")){ //i3_image * this_psf_image = i3_read_fits_image(input_psf_filename); data_set = i3_build_dataset_with_psf_image(options, identifier, stamp, weight_stamp, psf_image); //i3_image_destroy(this_psf_image); } if(!strcmp(options->psf_input,"moffat_catalog")){ if (i>=npsf) I3_FATAL("There were not enough lines in the PSF catalog for all the objects",i); psf_func = cat_psf_moffat + i; i3_print(i3_verb_standard,"i3_moffat_psf parameters used for galaxy %d: %f %f %f %f %d %d\n", i, psf_func->beta, psf_func->fwhm, psf_func->e1, psf_func->e2, psf_func->x, psf_func->y); data_set = i3_build_dataset_truncated_psf(options,identifier,stamp,weight_stamp,psf_func,psf_type,truncation); } //Modify the PSF FWHM definition to match the GREAT08 (and therefore also GREAT10?) simulation code //this_psf->fwhm = this_psf->fwhm * pow(1 - this_psf->e1*this_psf->e1 - this_psf->e2*this_psf->e2 , 0.5); //printf("%2.2f %2.2f %2.2f %2.2f \n",this_psf->beta,this_psf->fwhm,this_psf->e1,this_psf->e2); //Now run the main image analysis on this postage stamp - see i3_analyze.c //i3_flt mincheck = i3_image_min(stamp); //printf("mincheck = %20.10f\n",mincheck); //i3_data_set * data_set = i3_build_dataset(options,stamp,weight_stamp,this_psf,psf_type); // char optname[256]; // snprintf(optname,256,"%s/c_options.txt",options->output_directory); // i3_options_save(data_set->options,optname); /* for debugging (Michael) i3_image_print_first_element_of_each_row(stamp); */ // argh haven't thought about the value of weight_stamp for non-GREAT10 images... params_ml = i3_analyze_dataset_method(data_set,model,options,&bestL,bestI,method); //Calculate minimum and maximum of residuals map i3_image * residuals = i3_image_weighted_add_image(stamp,1.0,bestI,-1.0); i3_flt min_residuals = i3_image_min_masked(residuals, data_set->weight); i3_flt max_residuals = i3_image_max_masked(residuals, data_set->weight); // get additional image information i3_flt image_info = 0.0; #if !NO_SERSICS i3_flt image_args = 0.5; if(0==strcmp(model->name,"sersics") && 0==strcmp(options->psf_input,"moffat_catalog")){ i3_sersics_get_image_info(params_ml, data_set, &image_args, &image_info); } else{ image_info = i3_image_standard_fwhm(bestI); } #endif i3_flt model_min = i3_image_min(bestI); i3_flt model_max = i3_image_max(bestI); clock_t end_time = clock(); double time_taken = (double)(end_time-start_time) / CLOCKS_PER_SEC; //Save the results save_result_im3shape(options,model,params_ml,bestL,output,row,data_set,min_residuals,max_residuals, model_min, model_max, time_taken, &image_info); if(options->verbosity>=i3_verb_noisy){ printf("im3shape ML parameters: \n"); } if(options->verbosity>=i3_verb_standard){ save_result_im3shape(options,model,params_ml,bestL,stdout,row,data_set,min_residuals,max_residuals, model_min, model_max, time_taken, &image_info); } if (options->save_images){ char filename[256]; snprintf(filename,256,"%s/model.%ld.fits",options->output_directory,identifier); i3_image_save_fits(bestI,filename); snprintf(filename,256,"%s/image.%ld.fits",options->output_directory,identifier); i3_image_save_fits(stamp,filename); snprintf(filename,256,"%s/weight.%ld.fits",options->output_directory,identifier); i3_image_save_fits(weight_stamp,filename); //snprintf(filename,256,"%s/psf.%ld.fits",options->output_directory,identifier); //i3_image_save_fits(psf_image,filename); #if !NO_SERSICS if(options->save_sersics && !strcmp(model->name,"sersics")) { i3_sersics_model_image_save_components((i3_sersics_parameter_set*) params_ml, data_set, bestI, true); } #endif } //Clean up fflush(output); free(params_ml); i3_dataset_destroy(data_set); i3_image_destroy(stamp); i3_image_destroy(residuals); i3_image_destroy(galaxy); } // End loop over rows in galaxy catalogue file // } if(options->timeit){ time_end = time(NULL); cpu_time_used = (double)(time_end - time_start); printf("Total time for %d galaxies: %4.2fm\ttime per galaxy: %2.2fs\n",gals_counter, cpu_time_used/60.,cpu_time_used/(i3_flt)gals_counter); }; fclose(output); //free(cat_psf_moffat); i3_image_destroy(bestI); i3_model_destroy(model); i3_image_destroy(weight_stamp); fits_close_file(image_file, &status); i3_catalog_destroy(galaxy_cat); free(options); return 0; } int main_generate(int argc, char * argv[]) { if (argc<3){ fprintf(stderr, "Syntax for generating images:\nim3shape --generate params.ini out.fits param_name_1=param_value_1 param_name_2=param_value_2 ... \n"); return 1; } char * parameter_file = argv[2]; char * output_filename = argv[3]; i3_options * options = i3_options_default(); i3_options_read(options, parameter_file); i3_set_verbosity(options->verbosity); i3_model * model = i3_model_create(options->model_name, options); if(!model){ fprintf(stderr,"Could not find a model with name %s.\n",options->model_name); return 2; } i3_parameter_set * params = i3_model_option_starts(model, options); i3_moffat_psf * psf = i3_small_round_psf(); i3_flt noise_sigma=0.0; int error = parse_generate_parameters(argc, argv, 4, model, params, psf, &noise_sigma); if (error) return error; // Create empty image i3_image * stamp = i3_image_create(options->stamp_size, options->stamp_size); i3_image_fill(stamp, 0.0); // and constant weight. Should not matter too much i3_image * weight = i3_image_create(options->stamp_size, options->stamp_size); i3_image_fill(weight, 1.0); // Collect together into dataset i3_data_set * data_set = i3_build_dataset_truncated_psf( options, 1, stamp, weight, psf, options->airy_psf ? GREAT10_PSF_TYPE_AIRY : GREAT10_PSF_TYPE_MOFFAT, options->psf_truncation_pixels); // Generate image in likelihood function model->likelihood(stamp,params,data_set); // Add some white noise, if desired if (noise_sigma != 0.0) i3_image_add_white_noise(stamp, noise_sigma); // Save result error = i3_image_save_fits(stamp, output_filename); return error; } int parse_generate_parameters(int argc, char * argv[], int first_arg, i3_model * model, i3_parameter_set * params, i3_moffat_psf * psf, i3_flt *noise_sigma){ for (int i=first_arg; iname); return 4; } } } return 0; }