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
00033
00034
00035
00036
00037 #include "sinfo_utl_cube2spectrum.h"
00038 #include "sinfo_functions.h"
00039 #include "sinfo_spectrum_ops.h"
00040 #include "sinfo_key_names.h"
00041 #include "sinfo_pro_save.h"
00042 #include "sinfo_utilities_scired.h"
00043 #include "sinfo_globals.h"
00044 #include "sinfo_error.h"
00045 #include "sinfo_utils_wrappers.h"
00046
00047
00048
00049
00050
00051
00052
00053
00054 static void
00055 sinfo_change_header(const char * name);
00063
00064
00065
00066
00067 static void
00068 sinfo_change_header(const char * name) {
00069 int cenpix=0;
00070 double cenLambda=0;
00071 double dispersion=0;
00072 cpl_propertylist* plist=NULL;
00073
00074 plist=cpl_propertylist_load(name,0);
00075 cenpix = sinfo_pfits_get_crpix3(plist);
00076 cenLambda = sinfo_pfits_get_crval3(plist);
00077 dispersion = sinfo_pfits_get_cdelt3(plist);
00078
00079 cpl_propertylist_set_string(plist,"CTYPE2","WAVE");
00080 cpl_propertylist_set_comment(plist,"CTYPE2","wavelength axis in microns");
00081
00082 cpl_propertylist_set_double(plist,"CRPIX2",(double)cenpix);
00083 cpl_propertylist_set_comment(plist,"CRPIX2","Reference pixel");
00084
00085 cpl_propertylist_set_double(plist,"CRVAL2",cenLambda);
00086 cpl_propertylist_set_comment(plist,"CRVAL2","central wavelength");
00087
00088 cpl_propertylist_set_double(plist,"CDELT2",dispersion);
00089 cpl_propertylist_set_comment(plist,"CDELT2","microns per pixel");
00090
00091 cpl_propertylist_set_string(plist,"CUNIT2","mum");
00092 cpl_propertylist_set_comment(plist,"CUNIT2","spectral unit");
00093
00094 cpl_propertylist_set_string(plist,"CTYPE1","ONESPEC");
00095 cpl_propertylist_set_comment(plist,"CTYPE1","one spectrum in y-direction");
00096
00097 cpl_propertylist_erase_regexp(plist, "CRPIX1",0);
00098 cpl_propertylist_erase_regexp(plist, "CRVAL1",0);
00099 cpl_propertylist_erase_regexp(plist, "CDELT1",0);
00100 cpl_propertylist_erase_regexp(plist, "CUNIT1",0);
00101
00102 cpl_propertylist_erase_regexp(plist, "CTYPE3",0);
00103 cpl_propertylist_erase_regexp(plist, "CRPIX3",0);
00104 cpl_propertylist_erase_regexp(plist, "CRVAL3",0);
00105 cpl_propertylist_erase_regexp(plist, "CDELT3",0);
00106 cpl_propertylist_erase_regexp(plist, "CUNIT3",0);
00107
00108 cpl_propertylist_erase_regexp(plist, "CD1_1",0);
00109 cpl_propertylist_erase_regexp(plist, "CD1_2",0);
00110 cpl_propertylist_erase_regexp(plist, "CD2_1",0);
00111 cpl_propertylist_erase_regexp(plist, "CD2_2",0);
00112 sinfo_free_propertylist(&plist);
00113
00114 }
00115
00116
00124
00125 int sinfo_utl_cube2spectrum(
00126 cpl_parameterlist * parlist,
00127 cpl_frameset * framelist,
00128 const char* tag)
00129 {
00130 cpl_parameter * param =NULL;
00131 const char * operation =NULL;
00132 const char * aperture =NULL;
00133 char name_i [MAX_NAME_SIZE];
00134 int llx=0;
00135 int urx=0;
00136 int lly=0;
00137 int ury=0;
00138 int lo_rej=0;
00139 int hi_rej=0;
00140 int centerx=0;
00141 int centery=0;
00142 int radius=0;
00143 int clx=0;
00144 int cly=0;
00145
00146 double disp=0;
00147 double cpix=0;
00148 double clam=0;
00149
00150 cpl_frame * frm_cub=NULL;
00151
00152 char ima_o[MAX_NAME_SIZE];
00153 char tbl_o[MAX_NAME_SIZE];
00154 char tag_i[MAX_NAME_SIZE];
00155 char tag_o[MAX_NAME_SIZE];
00156 cpl_propertylist * plist =NULL;
00157 cpl_image * image =NULL;
00158 cpl_frame * product_frame=NULL;
00159
00160 cpl_image * im_spec=NULL;
00161 cpl_table * tbl_spec=NULL;
00162 cpl_image * img=NULL;
00163 cpl_imagelist * cube=NULL;
00164 Vector * spec=NULL;
00165
00166
00167
00168
00169 check_nomsg(param=cpl_parameterlist_find(parlist,
00170 "sinfoni.sinfo_utl_cube2spectrum.op"));
00171 check_nomsg(operation=cpl_parameter_get_string(param));
00172
00173
00174 check_nomsg(param=cpl_parameterlist_find(parlist,
00175 "sinfoni.sinfo_utl_cube2spectrum.ap"));
00176 check_nomsg(aperture=cpl_parameter_get_string(param));
00177
00178
00179
00180 if (tag == NULL) {
00181 strcpy(ima_o,"out_spec_ima.fits");
00182 strcpy(tbl_o,"out_spec_tbl.fits");
00183 strcpy(tag_i,SI_UTL_CUBE2SPECTRUM_CUBE);
00184 strcpy(tag_o,SI_UTL_CUBE2SPECTRUM_PROIMA);
00185 } else {
00186 snprintf(ima_o,MAX_NAME_SIZE-1,"%s%s",tag,"_spec_ima.fits");
00187 snprintf(tbl_o,MAX_NAME_SIZE-1,"%s%s",tag,"_spec_tbl.fits");
00188 snprintf(tag_o,MAX_NAME_SIZE-1,"%s%s",tag,"_SPCT");
00189 strcpy(tag_i, tag);
00190 }
00191
00192
00193
00194 check_nomsg(param=cpl_parameterlist_find(parlist,
00195 "sinfoni.sinfo_utl_cube2spectrum.llx"));
00196 check_nomsg(llx=cpl_parameter_get_int(param));
00197
00198
00199 check_nomsg(param=cpl_parameterlist_find(parlist,
00200 "sinfoni.sinfo_utl_cube2spectrum.lly"));
00201 check_nomsg(lly=cpl_parameter_get_int(param));
00202
00203
00204 check_nomsg(param=cpl_parameterlist_find(parlist,
00205 "sinfoni.sinfo_utl_cube2spectrum.urx"));
00206 check_nomsg(urx=cpl_parameter_get_int(param));
00207
00208
00209 check_nomsg(param=cpl_parameterlist_find(parlist,
00210 "sinfoni.sinfo_utl_cube2spectrum.ury"));
00211 check_nomsg(ury=cpl_parameter_get_int(param));
00212
00213
00214 check_nomsg(param=cpl_parameterlist_find(parlist,
00215 "sinfoni.sinfo_utl_cube2spectrum.lo_rej"));
00216 check_nomsg(lo_rej=cpl_parameter_get_int(param));
00217
00218
00219 check_nomsg(param=cpl_parameterlist_find(parlist,
00220 "sinfoni.sinfo_utl_cube2spectrum.hi_rej"));
00221 check_nomsg(hi_rej=cpl_parameter_get_int(param));
00222
00223
00224 check_nomsg(param=cpl_parameterlist_find(parlist,
00225 "sinfoni.sinfo_utl_cube2spectrum.centerx"));
00226 check_nomsg(centerx=cpl_parameter_get_int(param));
00227
00228
00229 check_nomsg(param=cpl_parameterlist_find(parlist,
00230 "sinfoni.sinfo_utl_cube2spectrum.centery"));
00231 check_nomsg(centery=cpl_parameter_get_int(param));
00232
00233
00234 check_nomsg(param=cpl_parameterlist_find(parlist,
00235 "sinfoni.sinfo_utl_cube2spectrum.radius"));
00236 check_nomsg(radius=cpl_parameter_get_int(param));
00237
00238
00239 if (sinfo_dfs_set_groups(framelist)) {
00240 sinfo_msg_error( "Cannot identify RAW and CALIB frames") ;
00241 goto cleanup;
00242 }
00243
00244
00245 cknull(frm_cub=cpl_frameset_find(framelist,tag_i),
00246 "SOF does not have a file tagged as %s",tag_i);
00247
00248
00249
00250 check_nomsg(strcpy(name_i,cpl_frame_get_filename(frm_cub)));
00251 check_nomsg(cube = cpl_imagelist_load((char*)name_i,CPL_TYPE_FLOAT,0));
00252
00253 check_nomsg(img=cpl_imagelist_get(cube,0));
00254 check_nomsg(clx=cpl_image_get_size_x(img));
00255 check_nomsg(cly=cpl_image_get_size_y(img));
00256 check(plist=cpl_propertylist_load(name_i,0),
00257 "Cannot read the FITS header") ;
00258
00259 cpix = (double) sinfo_pfits_get_crpix3(plist);
00260 clam = (double) sinfo_pfits_get_crval3(plist);
00261 disp = (double) sinfo_pfits_get_cdelt3(plist);
00262 sinfo_free_propertylist(&plist);
00263 if(strcmp(aperture,"rectangle") ==0) {
00264 if (llx<0) {
00265 sinfo_msg_warning("llx=%d too low set it to 0",llx);
00266 llx=0;
00267 }
00268 if (lly<0) {
00269 sinfo_msg_warning("lly=%d too low set it to 0",lly);
00270 lly=0;
00271 }
00272 if (urx>clx-1) {
00273 sinfo_msg_warning("urx=%d too large set it to %d",urx,clx-1);
00274 urx=clx-1;
00275 }
00276 if (ury>cly-1) {
00277 sinfo_msg_warning("ury=%d too large set it to %d",ury,cly-1);
00278 ury=cly-1;
00279 }
00280 if(llx>=urx) {
00281 sinfo_msg_error("llx>=urx!");
00282 goto cleanup;
00283 }
00284 if(lly>=ury) {
00285 sinfo_msg_error("lly>=ury!");
00286 goto cleanup;
00287 }
00288
00289 }
00290
00291 if(strcmp(aperture,"circle") ==0) {
00292 if((centerx-radius) < 0) {
00293 sinfo_msg_error("It is not possible to set centerx-radius<0.");
00294 goto cleanup;
00295 }
00296
00297 if((centery-radius) < 0) {
00298 sinfo_msg_error("It is not possible to set centery-radius<0.");
00299 goto cleanup;
00300 }
00301
00302 if((centerx+radius) >= clx) {
00303 sinfo_msg_error("It is not possible to set centerx+radius >= cube x sixe");
00304 goto cleanup;
00305 }
00306
00307 if((centery+radius) >= cly) {
00308 sinfo_msg_error("It is not possible to set centery+radius >= cube y size.");
00309 goto cleanup;
00310 }
00311
00312 if((radius) < 0) {
00313 sinfo_msg_error("It is not possible to set radius<0.");
00314 goto cleanup;
00315 }
00316
00317 }
00318
00319 if(strcmp(operation,"average") ==0) {
00320 if (strcmp(aperture,"rectangle") ==0) {
00321 spec = sinfo_new_mean_rectangle_of_cube_spectra(cube,llx,lly,urx,ury);
00322 } else if (strcmp(aperture,"circle")==0) {
00323 spec = sinfo_new_mean_circle_of_cube_spectra(cube,centerx,
00324 centery,radius);
00325 } else {
00326 sinfo_msg_error("Aperture not supported");
00327 sinfo_msg("Supported apertures are only:");
00328 sinfo_msg("rectangle, circle:");
00329 goto cleanup;
00330 }
00331 } else if (strcmp(operation,"clean_mean") ==0) {
00332 if (strcmp(aperture,"rectangle") == 0) {
00333 spec = sinfo_new_clean_mean_rectangle_of_cube_spectra(cube,llx,lly,
00334 urx,ury,
00335 lo_rej,hi_rej);
00336 } else if (strcmp(aperture,"circle")==0) {
00337 spec = sinfo_new_clean_mean_circle_of_cube_spectra(cube,centerx,
00338 centery,radius,
00339 lo_rej,hi_rej);
00340 } else {
00341 sinfo_msg_error("Aperture not supported");
00342 sinfo_msg("Supported apertures are only:");
00343 sinfo_msg("rectangle, circle:");
00344 goto cleanup;
00345 }
00346 } else if (strcmp(operation,"median") ==0) {
00347 if (strcmp(aperture,"rectangle")==0) {
00348 spec = sinfo_new_median_rectangle_of_cube_spectra(cube,llx,lly,
00349 urx,ury);
00350 } else if (strcmp(aperture,"circle")==0) {
00351 spec = sinfo_new_median_circle_of_cube_spectra(cube,centerx,
00352 centery,radius);
00353 } else {
00354 sinfo_msg_error("Aperture not supported");
00355 sinfo_msg("Supported apertures are only:");
00356 sinfo_msg("rectangle, circle:");
00357 goto cleanup;
00358 }
00359 } else if (strcmp(operation,"sum") ==0) {
00360 if (strcmp(aperture,"rectangle")==0) {
00361 spec = sinfo_new_sum_rectangle_of_cube_spectra(cube,llx,lly,urx,ury);
00362 } else if (strcmp(aperture,"circle")==0) {
00363 spec = sinfo_new_sum_circle_of_cube_spectra(cube,centerx,
00364 centery,radius);
00365 } else {
00366 sinfo_msg_error("Aperture not supported");
00367 sinfo_msg("Supported apertures are only:");
00368 sinfo_msg("rectangle, circle:");
00369 goto cleanup;
00370 }
00371 } else if (strcmp(operation,"extract") == 0) {
00372 if (strcmp(aperture,"rectangle")==0) {
00373 spec = sinfo_new_median_rectangle_of_cube_spectra(cube,llx,lly,
00374 urx,ury);
00375 } else if (strcmp(aperture,"circle")==0) {
00376 spec = sinfo_new_median_circle_of_cube_spectra(cube,centerx,
00377 centery,radius);
00378 } else {
00379 sinfo_msg_error("Aperture not supported");
00380 sinfo_msg("Supported apertures are only:");
00381 sinfo_msg("rectangle, circle:");
00382 goto cleanup;
00383 }
00384 } else {
00385 sinfo_msg_error("Operation not supported");
00386 sinfo_msg("Supported operations are only:");
00387 sinfo_msg("average, clean_mean, median, sum, extract :");
00388 goto cleanup;
00389 }
00390 im_spec = sinfo_new_vector_to_image(spec);
00391
00392 check_nomsg(sinfo_change_header(name_i));
00393
00394
00395
00396
00397
00398 ck0(sinfo_pro_save_ima(im_spec,framelist,framelist,ima_o,
00399 tag_o,NULL,cpl_func,parlist),"failed to save ima");
00400
00401
00402 sinfo_new_set_wcs_spectrum (im_spec,ima_o,clam, disp, cpix);
00403
00404 sinfo_stectrum_ima2table(im_spec,ima_o,&tbl_spec);
00405 ck0(sinfo_pro_save_tbl(tbl_spec,framelist,framelist,tbl_o,
00406 tag_o,NULL,cpl_func,parlist),
00407 "failed to save spectrum");
00408
00409 cleanup:
00410 sinfo_free_propertylist(&plist) ;
00411 sinfo_free_frame(&product_frame) ;
00412 sinfo_free_image(&image) ;
00413 sinfo_free_imagelist(&cube);
00414 sinfo_free_image(&im_spec);
00415
00416
00417
00418 sinfo_free_table(&tbl_spec);
00419
00420
00421
00422 if (cpl_error_get_code())
00423 return -1 ;
00424 else
00425 return 0 ;
00426
00427
00428 }