#include "i3_image_fits.h" #include "i3_sersics.h" #include "i3_mcmc.h" #include "math.h" #include "i3_math.h" #include "float.h" #include "i3_model_tools.h" #include "i3_sersics_definition.h" #include "i3_options.h" #include #include #include #include #include "i3_catalog.h" #include "i3_analyze.h" #define BEERMAT_MAX_E 0.95 // THIS VALUE SHOULD BE TESTED/OPTIMIZED #define BEERMAT_R_PARAM 0.2 // THIS VALUE SHOULD BE TESTED/OPTIMIZED /* * For a discussion of the meaning of the beermat R param see... TODO PUT SOME DESCRIPTION AND REF */ /* void i3_sersics_init_interp(gsl_vector * x1,gsl_vector * x2,gsl_matrix * X,dataset){ // // // int n_sub = dataset->upsampling; // int n_pix = dataset->stamp_size; // int n_pad = dataset->padding; // padding is adding // int n_all = n_pix+n_pad; // int n_win = (n_pad+1) * n_sub; // t // // static gsl_vector * i3_sersics_interp_x1 = gsl_vector_alloc(n_win*n_win); // static gsl_vector * i3_sersics_interp_x2 = gsl_vector_alloc(n_win*n_win); // static gsl_matrix * i3_sersics_interp_X = get_polynomial_basis( i3_sersics_interp_x1, i3_sersics_interp_x2 ); // // x1 = i3_sersics_interp_x1; // x2 = i3_sersics_interp_x1; // X = i3_sersics_interp_X; // // // } */ /* gsl_matrix * get_polynomial_basis(gsl_vector * vec_x1, gsl_vector * vec_x2){ int nx = vec_x1->size; int nf = 36; // number of features - 36 fro 2D polynomial expansion up to 7th order gsl_matrix * X = gsl_matrix_alloc(nx, nf); // double * x1 = vec_x1->data; // double * x2 = vec_x2->data; for(int i=0;isize1; // number of data points int nf = X->size2; // number of features in polynomial expansion // prepare the output matrices gsl_matrix * cov = gsl_matrix_alloc(nf,nf); gsl_matrix_set_all(cov,1.); gsl_vector * c = gsl_vector_alloc(nf); double chisq = 0; double alpha = 1.; double beta = 0.; // do least squares gsl_multifit_linear_workspace * lls_workspace = gsl_multifit_linear_alloc(np, nf); gsl_multifit_linear(X, y, c, cov, &chisq, lls_workspace); // prepare the data for the grid search on the interpolated space - this should be done somewhere outside the optimiser loop // this bit probably can be faster // double dxy = 1./(double)n_sub; int index_max_l = gsl_vector_max_index(y); //printf("index_max_l %d\n",index_max_l); double min_g1 = gsl_vector_get(x1,index_max_l) - dxy ; double max_g1 = gsl_vector_get(x1,index_max_l) + dxy ; double min_g2 = gsl_vector_get(x2,index_max_l) - dxy ; double max_g2 = gsl_vector_get(x2,index_max_l) + dxy ; //int ni = 200; // number of points for interpolation int ni2 = ni*ni; double dxyp = dxy/ni; //printf("creating interp grid %f %f %f %f %f\n",min_g1+n_pix/2,max_g1+n_pix/2,min_g2+n_pix/2,max_g2+n_pix/2,dxyp); int nij = 0; for(int i=0;ix0,paramset->y0,paramset->e1,paramset->e2,paramset->radius,paramset->radius_ratio,paramset->flux_ratio); // calculate how many likelihood evaluations are taken dataset->n_logL_evals++; // get the image parameters int n_sub = dataset->upsampling; int n_pix = dataset->stamp_size; int n_pad = dataset->padding; // padding is adding int n_all = n_pix+n_pad; int n_win = (n_pad+1) * n_sub; // t // init the model images i3_image * image_dvc = i3_image_create(n_all*n_sub,n_all*n_sub); i3_image_zero(image_dvc);; i3_image * image_exp = i3_image_create(n_all*n_sub,n_all*n_sub); i3_image_zero(image_exp); // for this model we always generate a model in the middle of the postage stamp paramset->x0 = n_pix/2; paramset->y0 = n_pix/2; // use the eta paratersitation i3_sersics_components( paramset,dataset,image_dvc,image_exp ); // set up the temp lo res image i3_image * img_tmp1 = i3_image_create(n_pix,n_pix); i3_image * img_tmp2 = i3_image_create(n_pix,n_pix); // set up the likelihood grid // starting to use linear regression notation phi([x1,x2]) = X, Xc=y gsl_vector * x1 = gsl_vector_alloc(n_win*n_win); gsl_vector * x2 = gsl_vector_alloc(n_win*n_win); gsl_vector * y = gsl_vector_alloc(n_win*n_win); double tx1, tx2, ty; // temporaty x1 x2 likelihood //printf("starting the loop - we should get %d values of xy_logL, n_win %d\n",n_win*n_win,n_win); int ns = 0; for(int na1 = 0; na1image, dataset->weight, &(dataset->amplitude1), &(dataset->amplitude2) ); //if(isnan(ty)) printf("logL max ampl 2 comp returned nan\n"); gsl_vector_set(x1,ns,tx1); gsl_vector_set(x2,ns,tx2); gsl_vector_set(y,ns,ty); //printf("n_win %d x1 %f y %f\n",n_win,gsl_vector_get(x1,ns),gsl_vector_get(x2,ns)); ns++; } } int ni = 200; // number of points to interpolate double dxy = 1./(double)n_sub; // delta xy gsl_vector * g1 = gsl_vector_alloc(ni*ni); gsl_vector * g2 = gsl_vector_alloc(ni*ni); gsl_vector * p = gsl_vector_alloc(ni*ni); gsl_vector_set_all(p,0.); get_interpolated_likelihood_in_xy0(x1,x2,y,g1,g2,p,ni,dxy); // find min and max of the interpolated values; int index_max_g = gsl_vector_max_index(p); i3_flt best_L = (i3_flt)gsl_vector_get(p,index_max_g); i3_flt best_x0 = (i3_flt)gsl_vector_get(g1,index_max_g) + n_pix/2; i3_flt best_y0 = (i3_flt)gsl_vector_get(g2,index_max_g) + n_pix/2; //if(isnan(best_L)) printf("interpolated likelihood is nan\n"); // Get the chi^2 value // if (!(dataset->image) || !(dataset->weight)) I3_FATAL("Image or weight not set for run of milestone model",1); // if (low_resolution_image_dvc->nx != dataset->image->nx || model_image->ny != dataset->image->ny) I3_FATAL("Model image and data image different sizes",1); // if (low_resolution_image_dvc->nx != dataset->weight->nx || model_image->ny != dataset->weight->ny) I3_FATAL("Model image and weight image different sizes. Possibly you have not set the weight image at all?",1); // // //JZ The function i3_logL_weight_map_maximise_amplitude_2components as a by-product puts the full combined model image into its first image argument, over-writing it. // i3_flt like = i3_logL_weight_map_maximise_amplitude_2components(low_resolution_image_dvc, low_resolution_image_exp, dataset->image, dataset->weight, &(dataset->amplitude1), &(dataset->amplitude2) ); // // i3_image_copy_into(low_resolution_image_dvc, model_image); // clean up i3_image_destroy( img_tmp1 ); i3_image_destroy( img_tmp2 ); i3_image_destroy( image_dvc ); i3_image_destroy( image_exp ); gsl_vector_free(x1); gsl_vector_free(x2); gsl_vector_free(y); gsl_vector_free(g1); gsl_vector_free(g2); gsl_vector_free(p); // some verb //paramset->x0 = best_x0; //paramset->y0 = best_y0; //printf("x0 y0 e1 e2 %f %f %f %f %f\n",paramset->x0,paramset->y0,paramset->e1, paramset->e2, paramset->radius); //i3_ellipticity_eta_to_e(paramset->e1,paramset->e2,&e1t,&e2t);paramset->e1=e1t; paramset->e2=e2t; //e1=paramset->e1;e2=paramset->e2; // printf("% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\n",best_L,best_x0,best_y0,paramset->e,paramset->theta,paramset->radius,paramset->radius_ratio,paramset->flux_ratio); //printf("% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\n",best_L,paramset->x0,paramset->y0,paramset->e1,paramset->e2,paramset->radius,paramset->radius_ratio,paramset->flux_ratio); return best_L; } */ void i3_sersics_model_image_with_moffat(i3_sersics_parameter_set * p, i3_options * options, i3_moffat_psf * psf_func, i3_image * model_image){ // This function is a convenience function for the python interface // where we don't have a i3_data_set object yet! // The below could be done without creating a data_set structure, // as all information is contained in the options struct. However, // for convenience and reducing code replication we do it like in // the im3shape driver int psf_type; // only needed when we are using moffat_catalog i3_flt truncation = options->psf_truncation_pixels; 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 } // Choose dummy identifier gal_id identifier = 0; //Make a dummy weight image. Assume noise level is constant across image i3_image * weight_stamp = i3_image_create(options->stamp_size,options->stamp_size); i3_image * stamp = i3_image_create(options->stamp_size,options->stamp_size); i3_data_set * data_set = i3_build_dataset_truncated_psf(options,identifier,stamp,weight_stamp,psf_func,psf_type,truncation); i3_sersics_model_image(p, data_set, model_image); i3_dataset_destroy(data_set); i3_image_destroy(weight_stamp); i3_image_destroy(stamp); } void i3_sersics_model_image_with_psf_image(i3_sersics_parameter_set * p, i3_options * options, i3_image * psf_image, i3_image * model_image){ // This function is a convenience function for the python interface // where we don't have a i3_data_set object yet! // The below could be done without creating a data_set structure, // as all information is contained in the options struct. However, // for convenience and reducing code replication we do it like in // the im3shape driver // The difference to i3_sersics_model_image_with_moffat is // that here we hand in a psf image rather than a moffat parameter struct // Choose dummy identifier gal_id identifier = 0; //Make a dummy weight image. Assume noise level is constant across image i3_image * weight_stamp = i3_image_create(options->stamp_size,options->stamp_size); i3_image * stamp = i3_image_create(options->stamp_size,options->stamp_size); i3_image_zero(stamp); i3_data_set * data_set = i3_build_dataset_with_psf_image(options, identifier, stamp, weight_stamp, psf_image); i3_sersics_model_image(p, data_set, model_image); i3_dataset_destroy(data_set); i3_image_destroy(weight_stamp); i3_image_destroy(stamp); } void i3_sersics_model_image(i3_sersics_parameter_set * p, i3_data_set * dataset, i3_image * model_image) { if ( dataset->options->use_hankel_transform ) //i3_sersics_model_image_using_hankel_transform(p, dataset, model_image, false); i3_sersics_model_image_save_components(p, dataset, model_image, false); else i3_sersics_model_image_save_components(p, dataset, model_image, false); } void i3_sersics_model_image_using_hankel_transform(i3_sersics_parameter_set * p, i3_data_set * dataset, i3_image * model_image, bool save_components){ // get the image parameters int n_sub = dataset->upsampling; int n_pix = dataset->stamp_size; int n_pad = dataset->padding; int n_all = n_pix+n_pad; // allocate memory for fourier image i3_fourier * fourier_image = i3_fourier_create(n_all, n_all); i3_fourier_zero(fourier_image); // Prepare the Sersic parameters for the bulge and the disc. i3_flt sersic_dvc = p->bulge_index; i3_flt sersic_exp = p->disc_index; i3_flt ab_dvc = (p->radius*n_sub)*(p->radius*n_sub); i3_flt bulge_disc_ratio = p->radius_ratio; i3_flt ab_exp = ab_dvc/bulge_disc_ratio/bulge_disc_ratio; i3_flt A_dvc = p->bulge_A; i3_flt A_exp = p->disc_A; // These are pixel center coordinates i3_flt x0 = p->x0 * n_sub + n_sub*n_pad/2; i3_flt y0 = p->y0 * n_sub + n_sub*n_pad/2; i3_flt e1_bulge = p->e1; i3_flt e2_bulge = p->e2; i3_flt e1_disc = p->e1; i3_flt e2_disc = p->e2; // These values are zero by default i3_flt delta_e_bulge = p->delta_e_bulge; i3_flt delta_theta_bulge = p->delta_theta_bulge; // The default value of these options is zero, meaning no truncation // i3_flt bulge_truncation = dataset->options->sersics_bulge_truncation * n_sub; // i3_flt disc_truncation = dataset->options->sersics_disc_truncation * n_sub; // make the components //If used, add to the bulge and disc components //This will break without warning if delta_e_bulge pushes us over the circle edge! if ((delta_e_bulge!=0)||(delta_theta_bulge!=0)){ i3_flt e = i3_sqrt(e1_bulge*e1_bulge+e2_bulge*e2_bulge); i3_flt theta = 0.5*i3_atan2(e2_bulge,e1_bulge); e += delta_e_bulge; theta += delta_theta_bulge; e1_bulge = e * i3_cos(2*theta); e2_bulge = e * i3_sin(2*theta); } if (A_dvc!=0.0) i3_add_fourier_space_sersic(ab_dvc, e1_bulge, e2_bulge, A_dvc, x0, y0, sersic_dvc, fourier_image); if (A_exp!=0.0) i3_add_fourier_space_sersic(ab_exp, e1_disc, e2_disc, A_exp, x0, y0, sersic_exp, fourier_image); /* Now get psf fourier image (at same resolution as model fourier image) * from dataset. Convolve (i.e. multiply) with model fourier image. * Then do the reverse fft. Won't doing this at low res cause aliasing? */ } void i3_sersics_model_image_save_components(i3_sersics_parameter_set * p, i3_data_set * dataset, i3_image * model_image, bool save_components){ // get the image parameters int n_sub = dataset->upsampling; int n_pix = dataset->stamp_size; int n_pad = dataset->padding; int n_all = n_pix+n_pad; // allocate memory for high res image i3_image * image_hires_pad = i3_image_create(n_all*n_sub,n_all*n_sub); i3_image_zero(image_hires_pad); // Prepare the Sersic parameters for the bulge and the disc. i3_flt sersic_dvc = p->bulge_index; i3_flt sersic_exp = p->disc_index; i3_flt ab_dvc = (p->radius*n_sub)*(p->radius*n_sub); i3_flt bulge_disc_ratio = p->radius_ratio; i3_flt ab_exp = ab_dvc/bulge_disc_ratio/bulge_disc_ratio; i3_flt A_dvc = p->bulge_A; i3_flt A_exp = p->disc_A; // These are pixel center coordinates i3_flt x0 = p->x0 * n_sub + n_sub*n_pad/2; i3_flt y0 = p->y0 * n_sub + n_sub*n_pad/2; i3_flt e1_bulge = p->e1; i3_flt e2_bulge = p->e2; i3_flt e1_disc = p->e1; i3_flt e2_disc = p->e2; // These values are zero by default i3_flt delta_e_bulge = p->delta_e_bulge; i3_flt delta_theta_bulge = p->delta_theta_bulge; // The default value of these options is zero, meaning no truncation i3_flt bulge_truncation = dataset->options->sersics_bulge_truncation * n_sub; i3_flt disc_truncation = dataset->options->sersics_disc_truncation * n_sub; // make the components //If used, add to the bulge and disc components //This will break without warning if delta_e_bulge pushes us over the circle edge! if ((delta_e_bulge!=0)||(delta_theta_bulge!=0)){ i3_flt e = i3_sqrt(e1_bulge*e1_bulge+e2_bulge*e2_bulge); i3_flt theta = 0.5*i3_atan2(e2_bulge,e1_bulge); e += delta_e_bulge; theta += delta_theta_bulge; e1_bulge = e * i3_cos(2*theta); e2_bulge = e * i3_sin(2*theta); } if(dataset->central_pixel_upsampling){ int n_central_upsampling = dataset->n_central_pixel_upsampling; int n_central_pixels_to_upsample = dataset->n_central_pixels_to_upsample; if(save_components){ // If you are saving the components then you cannot be too bothered about speed. // In that case we make the images twice, saving them once if (A_dvc!=0.0) i3_add_real_space_sersic_truncated_radius_upsample_central(image_hires_pad, ab_dvc, e1_bulge, e2_bulge, A_dvc, x0, y0, sersic_dvc, bulge_truncation, n_central_pixels_to_upsample, n_central_upsampling); char filename[256]; snprintf(filename, 256, "%s/bulge_%ld.fits", dataset->options->output_directory, dataset->identifier); i3_image_save_fits(image_hires_pad,filename); i3_image_zero(image_hires_pad); if (A_exp!=0.0) i3_add_real_space_sersic_truncated_radius_upsample_central(image_hires_pad, ab_exp, e1_disc, e2_disc, A_exp, x0, y0, sersic_exp, disc_truncation, n_central_pixels_to_upsample, n_central_upsampling); snprintf(filename, 256, "%s/disc_%ld.fits", dataset->options->output_directory, dataset->identifier); i3_image_save_fits(image_hires_pad,filename); i3_image_zero(image_hires_pad); } //i3_flt bulge_flux = 0.0; //i3_flt disc_flux = 0.0; //i3_flt total_flux = 0.0; if (A_dvc!=0.0) i3_add_real_space_sersic_truncated_radius_upsample_central(image_hires_pad, ab_dvc, e1_bulge, e2_bulge, A_dvc, x0, y0, sersic_dvc, bulge_truncation, n_central_pixels_to_upsample, n_central_upsampling); //for (int i=0;in;i++) bulge_flux+=image_hires_pad->data[i]; if (A_exp!=0.0) i3_add_real_space_sersic_truncated_radius_upsample_central(image_hires_pad, ab_exp, e1_disc, e2_disc, A_exp, x0, y0, sersic_exp, disc_truncation, n_central_pixels_to_upsample, n_central_upsampling); //for (int i=0;in;i++) total_flux+=image_hires_pad->data[i]; //disc_flux=total_flux-bulge_flux; //dataset->flux_ratio = bulge_flux/disc_flux; //printf("bulge flux %2.2f \n",bulge_flux); //printf("disc flux %2.2f \n",disc_flux); }else{ if(save_components){ if (A_dvc!=0.0) i3_add_real_space_sersic_truncated_radius(image_hires_pad, ab_dvc, e1_bulge, e2_bulge, A_dvc, x0, y0, sersic_dvc, bulge_truncation); char filename[256]; snprintf(filename, 256, "bulge_%ld.fits", dataset->identifier); i3_image_save_fits(image_hires_pad,filename); i3_image_zero(image_hires_pad); if (A_exp!=0.0) i3_add_real_space_sersic_truncated_radius(image_hires_pad, ab_exp, e1_disc, e2_disc, A_exp, x0, y0, sersic_exp, disc_truncation); snprintf(filename, 256, "disc_%ld.fits", dataset->identifier); i3_image_save_fits(image_hires_pad,filename); i3_image_zero(image_hires_pad); } if (A_dvc!=0.0) i3_add_real_space_sersic_truncated_radius(image_hires_pad, ab_dvc, e1_bulge, e2_bulge, A_dvc, x0, y0, sersic_dvc, bulge_truncation); if (A_exp!=0.0) i3_add_real_space_sersic_truncated_radius(image_hires_pad, ab_exp, e1_disc, e2_disc, A_exp, x0, y0, sersic_exp, disc_truncation); } // Convolve with the pre-prepared point spread function and Sinc kernel if( dataset->psf_downsampler_kernel ){ i3_image_convolve_fourier( image_hires_pad, dataset->psf_downsampler_kernel ); } // downsample int cut_start = 0; int cut_step = 0; //printf("tk3 %d %d %d %d \n",model_image->nx,n_all*n_sub,n_pix*n_sub,n_pix); if(model_image->nx == n_pix){ cut_start = (n_pad/2)*n_sub + n_sub/2; cut_step = n_sub; i3_image_dsample_cut_into( image_hires_pad, model_image, cut_start, cut_start, cut_step); }else if(model_image->nx == n_pix*n_sub){ cut_start = n_sub/2; cut_step = 1; i3_image_dsample_cut_into( image_hires_pad, model_image, cut_start, cut_start, cut_step); }else if(model_image->nx == n_all*n_sub){ cut_start = 0; cut_step = 1; i3_image_dsample_cut_into( image_hires_pad, model_image, cut_start, cut_start, cut_step); //i3_image_copy_into( image_hires_pad, model_image ); }else { printf("model_image %d image_hires_pad %d",(int)model_image->nx,(int)image_hires_pad->nx); I3_FATAL("model_image size is not correct\n",1); }; // clean up i3_image_destroy( image_hires_pad ); } void i3_sersics_model_jacobian(i3_sersics_parameter_set * p, i3_data_set * dataset, i3_flt * jac){ // get the image parameters int n_sub = dataset->upsampling; int n_pix = dataset->image->nx; int n_pad = dataset->padding; int n_all = n_pix+n_pad; int n = n_pix*n_pix; int m = 7; register int i, k; // allocate memory for high res image i3_image * jacobian[m]; i3_image * jacobian_hires[m]; i3_image * jacobian_hires_exp[m]; i3_image * jacobian_hires_dvc[m]; for (k=0; kbulge_index; i3_flt sersic_exp = p->disc_index; i3_flt ab_dvc = (p->radius*n_sub)*(p->radius*n_sub); i3_flt bulge_disc_ratio = p->radius_ratio; i3_flt ab_exp = ab_dvc/bulge_disc_ratio/bulge_disc_ratio; i3_flt truncation_factor = 40.0*n_sub; i3_flt A_dvc = p->bulge_A; i3_flt A_exp = p->disc_A; i3_flt x0 = p->x0*n_sub + n_pad*n_sub/2; i3_flt y0 = p->y0*n_sub + n_pad*n_sub/2; i3_flt e1_bulge = p->e1; i3_flt e2_bulge = p->e2; i3_flt e1_disc = p->e1; i3_flt e2_disc = p->e2; i3_flt delta_e_bulge = p->delta_e_bulge; i3_flt delta_theta_bulge = p->delta_theta_bulge; // make the components //If used, add to the bulge and disc components //This will break without warning if delta_e_bulge pushes us over the circle edge! if ((dataset->options->sersics_delta_e_bulge_fixed == 0)||(dataset->options->sersics_delta_e_bulge_fixed == 0)){ if ((delta_e_bulge!=0)||(delta_theta_bulge!=0)){ printf("Never reach this case\n"); i3_flt e = i3_sqrt(e1_bulge*e1_bulge+e2_bulge*e2_bulge); i3_flt theta = 0.5*i3_atan2(e2_bulge,e1_bulge); e += delta_e_bulge; theta += delta_theta_bulge; e1_bulge = e * i3_cos(2*theta); e2_bulge = e * i3_sin(2*theta); } } // ===================================================================================== // central pixel upsampling // ===================================================================================== if(dataset->central_pixel_upsampling){ int n_central_upsampling = dataset->n_central_pixel_upsampling; int n_central_pixels_to_upsample = dataset->n_central_pixels_to_upsample; // == Compute jacobian images ============================ i3_add_real_space_sersic_truncated_radius_upsample_central_jac(ab_dvc, e1_bulge, e2_bulge, 1.0, x0, y0, sersic_dvc, truncation_factor, n_central_pixels_to_upsample, n_central_upsampling, jacobian_hires_dvc); i3_add_real_space_sersic_truncated_radius_upsample_central_jac(ab_exp, e1_disc, e2_disc, 1.0, x0, y0, sersic_exp, truncation_factor, n_central_pixels_to_upsample, n_central_upsampling, jacobian_hires_exp); // == Compute jacobian image wrt x0,y0,e1,e2 ============================ for (k=0; k<4; k++){ i3_image_weighted_add_images_into(jacobian_hires_dvc[k], A_dvc, jacobian_hires_exp[k], A_exp, jacobian_hires[k] ); } // == Compute jacobian image wrt radius ============================ // Note the multiplicative scaling factors, ab_dvc and ab_exp are functions // of the radius itself, hence need to be multiplied with the corresponding // scaling factor according to the chain rule i3_image_weighted_add_images_into(jacobian_hires_dvc[4], 2*p->radius*n_sub*n_sub*A_dvc, jacobian_hires_exp[4], 2*p->radius*n_sub*n_sub/bulge_disc_ratio/bulge_disc_ratio*A_exp, jacobian_hires[4] ); // == Compute jacobian image wrt A_dvc,A_exp ============================ i3_image_add_image_into(jacobian_hires[5], jacobian_hires_dvc[5]); i3_image_add_image_into(jacobian_hires[6], jacobian_hires_exp[6]); } // ===================================================================================== // no central pixel upsampling // ===================================================================================== else{ /* // == Compute jacobian images wrt bulge and disc amplitude ============================ i3_add_real_space_sersic_truncated_radius(ab_dvc, e1_bulge, e2_bulge, 1.0, x0, y0, sersic_dvc, truncation_factor, jac_hires_pad_bulge_A); i3_add_real_space_sersic_truncated_radius(ab_exp, e1_disc, e2_disc, 1.0, x0, y0, sersic_exp, truncation_factor, jac_hires_pad_disc_A); */ } // ===================================================================================== // Convolve with the pre-prepared point spread function and Sinc kernel // ========================================= ============================================ if( dataset->psf_downsampler_kernel ){ for (k=0; kpsf_downsampler_kernel ); } } // ===================================================================================== // downsample // ===================================================================================== int cut_start = 0; int cut_step = 0; //printf("tk3 %d %d %d %d \n",model_image->nx,n_all*n_sub,n_pix*n_sub,n_pix); if(jacobian[0]->nx == n_pix){ cut_start = (n_pad/2)*n_sub + n_sub/2; cut_step = n_sub; for (k=0; knx == n_pix*n_sub){ cut_start = n_sub/2; cut_step = 1; for (k=0; knx == n_all*n_sub){ cut_start = 0; cut_step = 1; for (k=0; kdata[i] = jacobian[2]->data[i] - jacobian[7]->data[i]; jacobian[3]->data[i] = jacobian[3]->data[i] - jacobian[7]->data[i]; } */ // ===================================================================================== // Copy into jac // ===================================================================================== for (k=0; kdata[i]*n_sub; // Mulitply by n_sub to get scale right } } else{ for(i=0; idata[i]; } } } // ===================================================================================== // clean up // ===================================================================================== //printf("Here now %d\n\n",m); for (k=0; kupsampling; int n_pix = dataset->image->nx; int n_pad = dataset->padding; int n_all = n_pix+n_pad; int n = n_pix*n_pix; int m = 7; register int i, k; // allocate memory for high res image i3_image * jacobian[m]; i3_image * jacobian_hires[m]; i3_image * jacobian_hires_exp[m]; i3_image * jacobian_hires_dvc[m]; for (k=0; kbulge_index; i3_flt sersic_exp = p->disc_index; i3_flt ab_dvc = (p->radius*n_sub)*(p->radius*n_sub); i3_flt bulge_disc_ratio = p->radius_ratio; i3_flt ab_exp = ab_dvc/bulge_disc_ratio/bulge_disc_ratio; i3_flt truncation_factor = 40.0*n_sub; i3_flt A_dvc = p->bulge_A; i3_flt A_exp = p->disc_A; i3_flt x0 = p->x0*n_sub + n_pad*n_sub/2; i3_flt y0 = p->y0*n_sub + n_pad*n_sub/2; i3_flt e1_bulge = p->e1; i3_flt e2_bulge = p->e2; i3_flt e1_disc = p->e1; i3_flt e2_disc = p->e2; i3_flt delta_e_bulge = p->delta_e_bulge; i3_flt delta_theta_bulge = p->delta_theta_bulge; // make the components //If used, add to the bulge and disc components //This will break without warning if delta_e_bulge pushes us over the circle edge! if ((dataset->options->sersics_delta_e_bulge_fixed == 0)||(dataset->options->sersics_delta_e_bulge_fixed == 0)){ if ((delta_e_bulge!=0)||(delta_theta_bulge!=0)){ printf("Never reach this case\n"); i3_flt e = i3_sqrt(e1_bulge*e1_bulge+e2_bulge*e2_bulge); i3_flt theta = 0.5*i3_atan2(e2_bulge,e1_bulge); e += delta_e_bulge; theta += delta_theta_bulge; e1_bulge = e * i3_cos(2*theta); e2_bulge = e * i3_sin(2*theta); } } // ===================================================================================== // central pixel upsampling // ===================================================================================== if(dataset->central_pixel_upsampling){ int n_central_upsampling = dataset->n_central_pixel_upsampling; int n_central_pixels_to_upsample = dataset->n_central_pixels_to_upsample; // == Compute jacobian images ============================ i3_add_real_space_sersic_truncated_radius_upsample_central_jac_exact(ab_dvc, e1_bulge, e2_bulge, 1.0, x0, y0, sersic_dvc, truncation_factor, n_central_pixels_to_upsample, n_central_upsampling, jacobian_hires_dvc); i3_add_real_space_sersic_truncated_radius_upsample_central_jac_exact(ab_exp, e1_disc, e2_disc, 1.0, x0, y0, sersic_exp, truncation_factor, n_central_pixels_to_upsample, n_central_upsampling, jacobian_hires_exp); // == Compute jacobian image wrt x0,y0,e1,e2 ============================ for (k=0; k<4; k++){ i3_image_weighted_add_images_into(jacobian_hires_dvc[k], A_dvc, jacobian_hires_exp[k], A_exp, jacobian_hires[k] ); } // == Compute jacobian image wrt radius ============================ // Note the multiplicative scaling factors, ab_dvc and ab_exp are functions // of the radius itself, hence need to be multiplied with the corresponding // scaling factor according to the chain rule i3_image_weighted_add_images_into(jacobian_hires_dvc[4], 2*p->radius*n_sub*n_sub*A_dvc, jacobian_hires_exp[4], 2*p->radius*n_sub*n_sub/bulge_disc_ratio/bulge_disc_ratio*A_exp, jacobian_hires[4] ); // == Compute jacobian image wrt A_dvc,A_exp ============================ i3_image_add_image_into(jacobian_hires[5], jacobian_hires_dvc[5]); i3_image_add_image_into(jacobian_hires[6], jacobian_hires_exp[6]); } // ===================================================================================== // no central pixel upsampling // ===================================================================================== else{ /* // == Compute jacobian images wrt bulge and disc amplitude ============================ i3_add_real_space_sersic_truncated_radius(ab_dvc, e1_bulge, e2_bulge, 1.0, x0, y0, sersic_dvc, truncation_factor, jac_hires_pad_bulge_A); i3_add_real_space_sersic_truncated_radius(ab_exp, e1_disc, e2_disc, 1.0, x0, y0, sersic_exp, truncation_factor, jac_hires_pad_disc_A); */ } // ===================================================================================== // Convolve with the pre-prepared point spread function and Sinc kernel // ========================================= ============================================ if( dataset->psf_downsampler_kernel ){ for (k=0; kpsf_downsampler_kernel ); } } // ===================================================================================== // downsample // ===================================================================================== int cut_start = 0; int cut_step = 0; //printf("tk3 %d %d %d %d \n",model_image->nx,n_all*n_sub,n_pix*n_sub,n_pix); if(jacobian[0]->nx == n_pix){ cut_start = (n_pad/2)*n_sub + n_sub/2; cut_step = n_sub; for (k=0; knx == n_pix*n_sub){ cut_start = n_sub/2; cut_step = 1; for (k=0; knx == n_all*n_sub){ cut_start = 0; cut_step = 1; for (k=0; kdata[i] = jacobian[2]->data[i] - jacobian[7]->data[i]; jacobian[3]->data[i] = jacobian[3]->data[i] - jacobian[7]->data[i]; } */ // ===================================================================================== // Copy into jac // ===================================================================================== for (k=0; kdata[i]*n_sub; // Mulitply by n_sub to get scale right } } else{ for(i=0; idata[i]; } } } // ===================================================================================== // clean up // ===================================================================================== //printf("Here now %d\n\n",m); for (k=0; kupsampling; int n_pix = dataset->image->nx; int n_pad = dataset->padding; int n_all = n_pix+n_pad; int n = n_pix*n_pix; int m = 7; register int i, k; // allocate memory for high res image i3_image * jacobian[m]; i3_image * jacobian_hires[m]; i3_image * jacobian_hires_exp[m]; i3_image * jacobian_hires_dvc[m]; for (k=0; kbulge_index; i3_flt sersic_exp = p->disc_index; i3_flt ab_dvc = (p->radius*n_sub)*(p->radius*n_sub); i3_flt bulge_disc_ratio = p->radius_ratio; i3_flt ab_exp = ab_dvc/bulge_disc_ratio/bulge_disc_ratio; i3_flt truncation_factor = 40.0*n_sub; i3_flt A_dvc = p->bulge_A; i3_flt A_exp = p->disc_A; i3_flt x0 = p->x0*n_sub + n_pad*n_sub/2; i3_flt y0 = p->y0*n_sub + n_pad*n_sub/2; i3_flt e1_bulge = p->e1; i3_flt e2_bulge = p->e2; i3_flt e1_disc = p->e1; i3_flt e2_disc = p->e2; i3_flt delta_e_bulge = p->delta_e_bulge; i3_flt delta_theta_bulge = p->delta_theta_bulge; // make the components //If used, add to the bulge and disc components //This will break without warning if delta_e_bulge pushes us over the circle edge! if ((dataset->options->sersics_delta_e_bulge_fixed == 0)||(dataset->options->sersics_delta_e_bulge_fixed == 0)){ if ((delta_e_bulge!=0)||(delta_theta_bulge!=0)){ printf("Never reach this case\n"); i3_flt e = i3_sqrt(e1_bulge*e1_bulge+e2_bulge*e2_bulge); i3_flt theta = 0.5*i3_atan2(e2_bulge,e1_bulge); e += delta_e_bulge; theta += delta_theta_bulge; e1_bulge = e * i3_cos(2*theta); e2_bulge = e * i3_sin(2*theta); } } // ===================================================================================== // central pixel upsampling // ===================================================================================== if(dataset->central_pixel_upsampling){ int n_central_upsampling = dataset->n_central_pixel_upsampling; int n_central_pixels_to_upsample = dataset->n_central_pixels_to_upsample; // == Compute jacobian images ============================ i3_add_real_space_sersic_truncated_radius_upsample_central_jac_exact_vec(ab_dvc, e1_bulge, e2_bulge, 1.0, x0, y0, sersic_dvc, truncation_factor, n_central_pixels_to_upsample, n_central_upsampling, jacobian_hires_dvc); i3_add_real_space_sersic_truncated_radius_upsample_central_jac_exact_vec(ab_exp, e1_disc, e2_disc, 1.0, x0, y0, sersic_exp, truncation_factor, n_central_pixels_to_upsample, n_central_upsampling, jacobian_hires_exp); // == Compute jacobian image wrt x0,y0,e1,e2 ============================ for (k=0; k<4; k++){ i3_image_weighted_add_images_into(jacobian_hires_dvc[k], A_dvc, jacobian_hires_exp[k], A_exp, jacobian_hires[k] ); } // == Compute jacobian image wrt radius ============================ // Note the multiplicative scaling factors, ab_dvc and ab_exp are functions // of the radius itself, hence need to be multiplied with the corresponding // scaling factor according to the chain rule i3_image_weighted_add_images_into(jacobian_hires_dvc[4], 2*p->radius*n_sub*n_sub*A_dvc, jacobian_hires_exp[4], 2*p->radius*n_sub*n_sub/bulge_disc_ratio/bulge_disc_ratio*A_exp, jacobian_hires[4] ); // == Compute jacobian image wrt A_dvc,A_exp ============================ i3_image_add_image_into(jacobian_hires[5], jacobian_hires_dvc[5]); i3_image_add_image_into(jacobian_hires[6], jacobian_hires_exp[6]); } // ===================================================================================== // no central pixel upsampling // ===================================================================================== else{ /* // == Compute jacobian images wrt bulge and disc amplitude ============================ i3_add_real_space_sersic_truncated_radius(ab_dvc, e1_bulge, e2_bulge, 1.0, x0, y0, sersic_dvc, truncation_factor, jac_hires_pad_bulge_A); i3_add_real_space_sersic_truncated_radius(ab_exp, e1_disc, e2_disc, 1.0, x0, y0, sersic_exp, truncation_factor, jac_hires_pad_disc_A); */ } // ===================================================================================== // Convolve with the pre-prepared point spread function and Sinc kernel // ========================================= ============================================ if( dataset->psf_downsampler_kernel ){ for (k=0; kpsf_downsampler_kernel ); } } // ===================================================================================== // downsample // ===================================================================================== int cut_start = 0; int cut_step = 0; //printf("tk3 %d %d %d %d \n",model_image->nx,n_all*n_sub,n_pix*n_sub,n_pix); if(jacobian[0]->nx == n_pix){ cut_start = (n_pad/2)*n_sub + n_sub/2; cut_step = n_sub; for (k=0; knx == n_pix*n_sub){ cut_start = n_sub/2; cut_step = 1; for (k=0; knx == n_all*n_sub){ cut_start = 0; cut_step = 1; for (k=0; kdata[i] = jacobian[2]->data[i] - jacobian[7]->data[i]; jacobian[3]->data[i] = jacobian[3]->data[i] - jacobian[7]->data[i]; } */ // ===================================================================================== // Copy into jac // ===================================================================================== for (k=0; kdata[i]*n_sub; // Mulitply by n_sub to get scale right } } else{ for(i=0; idata[i]; } } } // ===================================================================================== // clean up // ===================================================================================== //printf("Here now %d\n\n",m); for (k=0; kimage; i3_image * weighted_image = i3_image_copy(image); // Find centroids from weighted quadrupole moments using narrow weight function i3_flt weight_radius = 3.0; int weight_iterations = 10; i3_image_compute_weighted_moments(image,weight_radius,weight_iterations,&Q); if (Q.x0 > options->sersics_x0_min && Q.x0 < options->sersics_x0_max) start->x0 = Q.x0; if (Q.y0 > options->sersics_y0_min && Q.y0 < options->sersics_y0_max) start->y0 = Q.y0; //printf("start->x0=%5.3f\t start->y0=%5.3f\n",start->x0,start->y0); // Estimate starting size. Do the correct procedure for Gaussians. weight_radius = 10.0; for (int i=0;iy0,start->x0,weight_radius); // Need to run compute_center to get the sum used by compute_quadruople i3_image_compute_center(weighted_image,&Q); i3_image_compute_quadrupole(weighted_image,&Q); if ((Q.qxx>0)&(Q.qyy>0)){ // We do times sqrt(2) to compensate for weight by itself weight_radius=pow(2.0,0.5) * pow(Q.qxx*Q.qyy,0.25); } else { // If here, must be picking up a lot of noise - so need smaller weight radius weight_radius=weight_radius/2; } } // fwhm = (2*(sqrt(2*log(2)))) * sigma for a Gaussian i3_flt convolved_fwhm_estimate=2.3548*weight_radius; // Turned that off because it was crushin when we were using image of PSF, no functional form // TODO: make estimation of unconvolved fwhm work with images of PSFs //i3_flt unconvolved_fwhm_estimate = 0.5*pow(convolved_fwhm_estimate*convolved_fwhm_estimate-data_set->psf_form->fwhm*data_set->psf_form->fwhm,0.5); // Modify the default start value if the quad moms output is within the minmax i3_flt unconvolved_fwhm_estimate = 0.5*convolved_fwhm_estimate; if (unconvolved_fwhm_estimate > options->sersics_radius_min && unconvolved_fwhm_estimate < options->sersics_radius_max){ start->radius=unconvolved_fwhm_estimate; } //printf("unconvolved_fwhm_estimate=%5.3f\t start->radius=%5.3f\n",unconvolved_fwhm_estimate,start->radius); //Estimate e and theta using weighted (uncorrected) quadrupole moments //Not yet totally convinced this is better than using e=0.2 and random theta... //start->e=0.2; //start->theta=1; // This is my random number for theta (in radians) NB. 0 is bad //printf("weight_radius=%5.3f\n",weight_radius); //i3_weight_image(image,weighted_image,start->y0,start->x0,weight_radius); // Need to run compute_center to get the sum used by compute_quadruople //i3_image_compute_center(weighted_image,&Q); //i3_imag // Clean up i3_image_destroy(weighted_image); } i3_flt i3_sersics_likelihood(i3_image * model_image, i3_sersics_parameter_set * paramset, i3_data_set * dataset){ // calculate how many likelihood evaluations are taken dataset->n_logL_evals++; // prepare the model image i3_image_zero( model_image ); // prepare the beermatted param set i3_sersics_parameter_set * paramset_beermat = i3_sersics_beermat_params(paramset,dataset->options); // get the model image i3_sersics_model_image( paramset_beermat, dataset, model_image); // Get the chi^2 value if (!(dataset->image) || !(dataset->weight)) I3_FATAL("Image or weight not set for run of sersics model",1); if (model_image->nx != dataset->image->nx || model_image->ny != dataset->image->ny ) I3_FATAL("Model image and data image different sizes",1); if (model_image->nx != dataset->weight->nx || model_image->ny != dataset->weight->ny) I3_FATAL("Model image and weight image different sizes. Possibly you have not set the weight image at all?",1); #ifdef SERSICS_TRUE_LIKE i3_flt chi2 = i3_chi2_weight_map(model_image, dataset->image, dataset->weight ); #else i3_flt chi2=666.; #endif // clean up //i3_flt flux_ratio = paramset->bulge_A/(paramset->bulge_A+paramset->disc_A); //if(isnan(logL)) fprintf(stderr,"nanerror \t % 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\n",logL,paramset->x0,paramset->y0,paramset->e1,paramset->e2,paramset->radius,paramset->radius_ratio,paramset->bulge_A,paramset->disc_A,paramset->bulge_index,paramset->disc_index,paramset->delta_e_bulge,paramset->delta_theta_bulge); if (dataset->options->verbosity>10) { fprintf(stderr,"% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e\t% 2.4e \n",chi2,paramset_beermat->x0,paramset_beermat->y0,paramset_beermat->e1,paramset_beermat->e2,paramset_beermat->radius,paramset_beermat->bulge_A,paramset_beermat->disc_A); } if (dataset->options->verbosity>11) { // Save the model image that was tried char filename[256]; snprintf( filename, 256, "%s/trial_model_image_%04d.fits", dataset->options->output_directory, dataset->n_logL_evals); i3_image_save_fits(model_image,filename); //snprintf( filename, 256, "%s/dataset_image_%04d.fits", dataset->options->output_directory, dataset->n_logL_evals); //i3_image_save_fits(dataset->image,filename); //snprintf( filename, 256, "%s/weight_image_%04d.fits", dataset->options->output_directory, dataset->n_logL_evals); //i3_image_save_fits(dataset->weight,filename); } free(paramset_beermat); return -chi2/2.0; } void i3_sersics_beermat_mapping(i3_sersics_parameter_set * input, i3_sersics_parameter_set * output, i3_options * options) { // Most parametersare the same in the two parameterizations memcpy(output, input, sizeof(i3_sersics_parameter_set)); // BARNEY: DO WE NEED TO REMEMBER TO FREE HP LATER ON USING FREE()? // NOW include the beermat changes, mapping only these internal variables i3_e12_to_e12beermat(input->e1, input->e2, BEERMAT_MAX_E, &(output->e1), &(output->e2)); i3_r_to_rbeermat(input->radius, BEERMAT_R_PARAM, &(output->radius)); if(options->sersics_beermat_amplitudes_positive){ // output->bulge_A = i3_fabs(input->bulge_A); // output->disc_A = i3_fabs(output->disc_A); i3_flt eps = 1e-1; output->bulge_A = sqrt(input->bulge_A* input->bulge_A + eps); output->disc_A = sqrt(input->disc_A * input->disc_A + eps); } } i3_sersics_parameter_set * i3_sersics_beermat_params(i3_sersics_parameter_set * p , i3_options * options){ i3_sersics_parameter_set * hp = malloc(sizeof(i3_sersics_parameter_set)); memcpy(hp, p, sizeof(i3_sersics_parameter_set)); // BARNEY: DO WE NEED TO REMEMBER TO FREE HP LATER ON USING FREE()? i3_flt e1 = p->e1; i3_flt e2 = p->e2; i3_flt r = p->radius; // NOW include the beermat changes, mapping only these internal variables i3_flt e1_beermat; i3_flt e2_beermat; i3_flt r_beermat; i3_e12_to_e12beermat(e1, e2, BEERMAT_MAX_E, &e1_beermat, &e2_beermat); i3_r_to_rbeermat(r, BEERMAT_R_PARAM, &r_beermat); hp->e1 = e1_beermat; hp->e2 = e2_beermat; hp->radius = r_beermat; if(options->sersics_beermat_amplitudes_positive){ // hp->bulge_A = i3_fabs(p->bulge_A); // hp->disc_A = i3_fabs(p->disc_A); // hp->bulge_A = p->bulge_A * p->bulge_A; // hp->disc_A = p->disc_A * p->disc_A; // fprintf(stdout,"aaa\n"); i3_flt eps = 1e-4; hp->bulge_A = sqrt(p->bulge_A* p->bulge_A + eps); hp->disc_A = sqrt(p->disc_A * p->disc_A + eps); } return hp; } void i3_sersics_get_image_info(i3_sersics_parameter_set * params, i3_data_set * dataset, i3_flt * image_args, i3_flt * image_info){ // get the image parameters int original_dataset_upsampling=dataset->upsampling ; dataset->upsampling = 9; int n_sub = dataset->upsampling; int n_pix = dataset->image->nx; int n_pad = dataset->padding; int n_all = n_pix+n_pad; // prepare the circular psf i3_moffat_psf psf_form_circular; psf_form_circular.beta = dataset->psf_form->beta; psf_form_circular.fwhm = dataset->psf_form->fwhm; psf_form_circular.e1 = 0.; psf_form_circular.e2 = 0.; // prepare the parameter set i3_sersics_parameter_set * params_circular = malloc(sizeof(i3_sersics_parameter_set)); memcpy(params_circular,params,sizeof(i3_sersics_parameter_set)); params_circular->x0 = n_pix/2.; params_circular->y0 = n_pix/2.; params_circular->e1 = 0; params_circular->e2 = 0; // create image space //i3_image * image_npsf = i3_image_create( n_all*n_sub, n_all*n_sub ); i3_image_zero( image_npsf ); i3_image * image_wpsf = i3_image_create( n_all*n_sub, n_all*n_sub ); i3_image_zero( image_wpsf ); i3_fourier * fourier_ker = dataset->psf_downsampler_kernel; i3_fourier * kernel_circular = i3_fourier_conv_kernel_moffat(n_pix, n_sub, n_pad, psf_form_circular.beta, psf_form_circular.fwhm, psf_form_circular.e1, psf_form_circular.e2, dataset->options->psf_truncation_pixels, GREAT10_PSF_TYPE_MOFFAT); //i3_build_great10_kernel(&psf_form_circular, n_pix, n_sub, n_pad, GREAT10_PSF_TYPE_MOFFAT); //i3_fourier * kernel_elliptic = i3_build_great10_kernel( dataset->psf_form, n_pix, n_sub, n_pad, GREAT10_PSF_TYPE_MOFFAT); // get the unconvolved image //dataset->psf_downsampler_kernel = NULL; //i3_sersics_model_image(params_circular,dataset,image_npsf); //if(dataset->options->save_images) i3_image_save_fits( image_npsf, "image_test_npsf.fits" ); // get the convolved image dataset->psf_downsampler_kernel = kernel_circular; i3_sersics_model_image(params_circular,dataset,image_wpsf); // if(dataset->options->save_images) i3_image_save_fits( image_wpsf, "image_test_wpsf.fits" ); // measure what needs to be measured image_info[0] = i3_image_get_fwhm(image_wpsf, 0.5, params_circular->x0/2.)/(i3_flt)n_sub ; //image_info[1] = measure_fwxm(image_wpsf, image_args[1])/(i3_flt)n_sub ; //image_info[2] = measure_fwxm(image_wpsf, image_args[2])/(i3_flt)n_sub ; //image_info[3] = measure_entropy(image_npsf); //image_info[4] = measure_entropy(image_wpsf); //image_info[5] = measure_kldivergence(image_npsf,image_wpsf); //image_info[6] = measure_kldivergence(image_wpsf,image_npsf); //image_info[3] = 66.; //image_info[4] = 66.; //image_info[5] = 66.; //image_info[6] = 66.; // do the same with elliptical profiles //dataset->psf_downsampler_kernel = NULL; //i3_sersics_model_image(params,dataset,image_npsf); //dataset->psf_downsampler_kernel = kernel_elliptic; //i3_sersics_model_image(params,dataset,image_wpsf); //image_info[7] = measure_kldivergence(image_npsf,image_wpsf); //image_info[8] = measure_kldivergence(image_wpsf,image_npsf); //image_info[7] = 66.; //image_info[8] = 66.; // clean up dataset->psf_downsampler_kernel = fourier_ker; dataset->upsampling = original_dataset_upsampling; i3_image_destroy( image_wpsf ); //i3_image_destroy( image_npsf ); i3_fourier_destroy( kernel_circular ); //i3_fourier_destroy( kernel_elliptic ); free(params_circular); }