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
00029
00030
00031
00032
00033
00034 #include "utl_cube2spectrum.h"
00035 #include "eclipse.h"
00036 #include "fitshead.h"
00037 #include "spectrum_ops.h"
00038 #include "sinfoni_key_names.h"
00039 #include "utilities_scired.h"
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052 void
00053 change_header(fits_header * header, char * name) {
00054 int cenpix=0;
00055 float cenLambda=0;
00056 float dispersion=0;
00057 char blank[FILENAMESZ] ;
00058
00059
00060 cenpix = spiffi_get_cenpix(name);
00061 cenLambda = spiffi_get_cenLambda(name);
00062 dispersion = spiffi_get_dispersion(name);
00063
00064 fits_header_mod(header, "CTYPE2", "WAVE" , "wavelength axis in microns");
00065
00066 sprintf(blank, "%d", cenpix) ;
00067 fits_header_mod(header, "CRPIX2", blank , "Reference pixel");
00068 sprintf(blank, "%f", cenLambda) ;
00069 fits_header_mod(header, "CRVAL2", blank , "central wavelength");
00070 sprintf(blank, "%f", dispersion) ;
00071 fits_header_mod(header, "CDELT2", blank , "microns per pixel");
00072
00073 fits_header_mod(header, "CUNIT2", "MICRON" , "spectral unit");
00074 fits_header_mod(header, "CTYPE1", "ONESPEC" , "one spectrum in y-direction");
00075 fits_header_del(header, "CRPIX1");
00076 fits_header_del(header, "CRVAL1");
00077 fits_header_del(header, "CDELT1");
00078 fits_header_del(header, "CUNIT1");
00079 fits_header_del(header, "CTYPE3");
00080 fits_header_del(header, "CRPIX3");
00081 fits_header_del(header, "CRVAL3");
00082 fits_header_del(header, "CDELT3");
00083 fits_header_del(header, "CUNIT3");
00084 fits_header_del(header, "CD1_1");
00085 fits_header_del(header, "CD1_2");
00086 fits_header_del(header, "CD2_1");
00087 fits_header_del(header, "CD2_2");
00088
00089 }
00090
00091
00098
00099 int si_utl_cube2spectrum(
00100 cpl_parameterlist * parlist,
00101 cpl_frameset * framelist)
00102 {
00103 const char * fctid = "si_utl_cube2spectrum" ;
00104 cpl_parameter * param =NULL;
00105 const char * operation =NULL;
00106 const char * aperture =NULL;
00107 const char * name_i=NULL;
00108 int llx=0;
00109 int urx=0;
00110 int lly=0;
00111 int ury=0;
00112 int lo_rej=0;
00113 int hi_rej=0;
00114 int centerx=0;
00115 int centery=0;
00116 int radius=0;
00117
00118 double disp=0;
00119 double cpix=0;
00120 double clam=0;
00121
00122 cpl_frame * frm_cub=NULL;
00123
00124 const char * name_o=NULL ;
00125 cpl_propertylist * plist =NULL;
00126 cpl_image * image =NULL;
00127 cpl_image * imagew =NULL;
00128 cpl_frame * product_frame=NULL;
00129
00130 OneImage * im_spec=NULL;
00131 OneCube * cube=NULL;
00132 Vector * spec=NULL;
00133 fits_header* head=NULL;
00134
00135
00136
00137 param = cpl_parameterlist_find(parlist, "sinfoni.si_utl_cube2spectrum.op");
00138 operation = cpl_parameter_get_string(param);
00139
00140
00141 param = cpl_parameterlist_find(parlist, "sinfoni.si_utl_cube2spectrum.ap");
00142 aperture = cpl_parameter_get_string(param);
00143
00144
00145 name_o = "out_spec.fits";
00146
00147
00148 param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.llx");
00149 llx = cpl_parameter_get_int(param) ;
00150
00151
00152 param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.lly");
00153 lly = cpl_parameter_get_int(param) ;
00154
00155
00156 param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.urx");
00157 urx = cpl_parameter_get_int(param) ;
00158
00159
00160 param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.ury");
00161 ury = cpl_parameter_get_int(param) ;
00162
00163
00164 param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.lo_rej");
00165 lo_rej = cpl_parameter_get_int(param) ;
00166
00167
00168 param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.hi_rej");
00169 hi_rej = cpl_parameter_get_int(param) ;
00170
00171
00172 param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.centerx");
00173 centerx = cpl_parameter_get_int(param) ;
00174
00175
00176 param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.centery");
00177 centery = cpl_parameter_get_int(param) ;
00178
00179
00180 param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.radius");
00181 radius = cpl_parameter_get_int(param) ;
00182
00183
00184
00185 if (sinfoni_dfs_set_groups(framelist)) {
00186 cpl_msg_error(fctid, "Cannot identify RAW and CALIB frames") ;
00187 return -1 ;
00188 }
00189
00190
00191 if ((frm_cub = cpl_frameset_find(framelist, SI_UTL_CUBE2SPECTRUM_CUBE))==NULL) {
00192 cpl_msg_error(fctid, "SOF does not have a file tagged as %s",SI_UTL_CUBE2SPECTRUM_CUBE);
00193 return -1 ;
00194 }
00195
00196
00197 if ((plist=cpl_propertylist_load(cpl_frame_get_filename(frm_cub),
00198 0)) == NULL) {
00199 cpl_msg_error(fctid, "Cannot read the FITS header") ;
00200 return -1 ;
00201 }
00202
00203
00204 name_i = cpl_frame_get_filename(frm_cub);
00205 cube = load_cube((char*)name_i);
00206 head = fits_read_header((char*)name_i);
00207
00208 cpix = (double) spiffi_get_cenpix((char*)name_i);
00209 clam = (double) spiffi_get_cenLambda((char*)name_i);
00210 disp = (double) spiffi_get_dispersion((char*)name_i);
00211
00212
00213
00214 if(strcmp(aperture,"rectangle") ==0) {
00215 if (llx<0) {
00216 cpl_msg_warning(fctid,"llx=%d too low set it to 0",llx);
00217 llx=0;
00218 }
00219 if (lly<0) {
00220 cpl_msg_warning(fctid,"lly=%d too low set it to 0",lly);
00221 lly=0;
00222 }
00223 if (urx>cube->lx-1) {
00224 cpl_msg_warning(fctid,"urx=%d too large set it to %d",urx,cube->lx-1);
00225 urx=cube->lx-1;
00226 }
00227 if (ury>cube->ly-1) {
00228 cpl_msg_warning(fctid,"ury=%d too large set it to %d",ury,cube->ly-1);
00229 ury=cube->ly-1;
00230 }
00231 if(llx>=urx) {
00232 cpl_msg_error(fctid,"llx>=urx!");
00233 cpl_propertylist_delete(plist) ;
00234 cpl_frame_delete(product_frame) ;
00235 cpl_image_delete(image) ;
00236 destroy_cube(cube);
00237 destroy_image(im_spec);
00238 fits_header_destroy(head);
00239 return -1;
00240 }
00241 if(lly>=ury) {
00242 cpl_msg_error(fctid,"lly>=ury!");
00243 cpl_propertylist_delete(plist) ;
00244 cpl_frame_delete(product_frame) ;
00245 cpl_image_delete(image) ;
00246 destroy_cube(cube);
00247 destroy_image(im_spec);
00248 fits_header_destroy(head);
00249 return -1;
00250 }
00251
00252 }
00253
00254 if(strcmp(aperture,"circle") ==0) {
00255 if((centerx-radius) < 0) {
00256 cpl_msg_error(fctid,"It is not possible to set centerx-radius<0.");
00257 cpl_propertylist_delete(plist) ;
00258 cpl_frame_delete(product_frame) ;
00259 cpl_image_delete(image) ;
00260 destroy_cube(cube);
00261 destroy_image(im_spec);
00262 fits_header_destroy(head);
00263 return -1;
00264 }
00265
00266 if((centery-radius) < 0) {
00267 cpl_msg_error(fctid,"It is not possible to set centery-radius<0.");
00268 cpl_propertylist_delete(plist) ;
00269 cpl_frame_delete(product_frame) ;
00270 cpl_image_delete(image) ;
00271 destroy_cube(cube);
00272 destroy_image(im_spec);
00273 fits_header_destroy(head);
00274 return -1;
00275 }
00276
00277 if((centerx+radius) >= cube->lx) {
00278 cpl_msg_error(fctid,"It is not possible to set centerx+radius >= cube x sixe");
00279 cpl_propertylist_delete(plist) ;
00280 cpl_frame_delete(product_frame) ;
00281 cpl_image_delete(image) ;
00282 destroy_cube(cube);
00283 destroy_image(im_spec);
00284 fits_header_destroy(head);
00285 return -1;
00286 }
00287
00288 if((centery+radius) >= cube->ly) {
00289 cpl_msg_error(fctid,"It is not possible to set centery+radius >= cube y size.");
00290 cpl_propertylist_delete(plist) ;
00291 cpl_frame_delete(product_frame) ;
00292 cpl_image_delete(image) ;
00293 destroy_cube(cube);
00294 destroy_image(im_spec);
00295 fits_header_destroy(head);
00296 return -1;
00297 }
00298
00299 if((radius) < 0) {
00300 cpl_msg_error(fctid,"It is not possible to set radius<0.");
00301 cpl_propertylist_delete(plist) ;
00302 cpl_frame_delete(product_frame) ;
00303 cpl_image_delete(image) ;
00304 destroy_cube(cube);
00305 destroy_image(im_spec);
00306 fits_header_destroy(head);
00307 return -1;
00308 }
00309
00310 }
00311 if(strcmp(operation,"average") ==0) {
00312 if (strcmp(aperture,"rectangle") ==0) {
00313 spec = meanRectangleOfCubeSpectra(cube,llx,lly,urx,ury);
00314 } else if (strcmp(aperture,"circle")==0) {
00315 spec = meanCircleOfCubeSpectra(cube,centerx,centery,radius);
00316 } else {
00317 cpl_msg_error(fctid,"Aperture not supported");
00318 cpl_msg_info(fctid,"Supported apertures are only:");
00319 cpl_msg_info(fctid,"rectangle, circle:");
00320 cpl_propertylist_delete(plist) ;
00321 cpl_frame_delete(product_frame) ;
00322 cpl_image_delete(image) ;
00323 destroy_cube(cube);
00324 destroy_image(im_spec);
00325 fits_header_destroy(head);
00326 return -1;
00327 }
00328 } else if (strcmp(operation,"clean_mean") ==0) {
00329 if (strcmp(aperture,"rectangle") == 0) {
00330 spec = cleanmeanRectangleOfCubeSpectra(cube,llx,lly,urx,ury,lo_rej,hi_rej);
00331 } else if (strcmp(aperture,"circle")==0) {
00332 spec = cleanmeanCircleOfCubeSpectra(cube,centerx,centery,radius,lo_rej,hi_rej);
00333 } else {
00334 cpl_msg_error(fctid,"Aperture not supported");
00335 cpl_msg_info(fctid,"Supported apertures are only:");
00336 cpl_msg_info(fctid,"rectangle, circle:");
00337 cpl_propertylist_delete(plist) ;
00338 cpl_frame_delete(product_frame) ;
00339 cpl_image_delete(image) ;
00340 destroy_cube(cube);
00341 destroy_image(im_spec);
00342 fits_header_destroy(head);
00343 return -1;
00344 }
00345 } else if (strcmp(operation,"median") ==0) {
00346 if (strcmp(aperture,"rectangle")==0) {
00347 spec = medianRectangleOfCubeSpectra(cube,llx,lly,urx,ury);
00348 } else if (strcmp(aperture,"circle")==0) {
00349 spec = medianCircleOfCubeSpectra(cube,centerx,centery,radius);
00350 } else {
00351 cpl_msg_error(fctid,"Aperture not supported");
00352 cpl_msg_info(fctid,"Supported apertures are only:");
00353 cpl_msg_info(fctid,"rectangle, circle:");
00354 cpl_propertylist_delete(plist) ;
00355 cpl_frame_delete(product_frame) ;
00356 cpl_image_delete(image) ;
00357 destroy_cube(cube);
00358 destroy_image(im_spec);
00359 fits_header_destroy(head);
00360 return -1;
00361 }
00362 } else if (strcmp(operation,"sum") ==0) {
00363 if (strcmp(aperture,"rectangle")==0) {
00364 spec = sumRectangleOfCubeSpectra(cube,llx,lly,urx,ury);
00365 } else if (strcmp(aperture,"circle")==0) {
00366 spec = sumCircleOfCubeSpectra(cube,centerx,centery,radius);
00367 } else {
00368 cpl_msg_error(fctid,"Aperture not supported");
00369 cpl_msg_info(fctid,"Supported apertures are only:");
00370 cpl_msg_info(fctid,"rectangle, circle:");
00371 cpl_propertylist_delete(plist) ;
00372 cpl_frame_delete(product_frame) ;
00373 cpl_image_delete(image) ;
00374 destroy_cube(cube);
00375 destroy_image(im_spec);
00376 fits_header_destroy(head);
00377 return -1;
00378 }
00379 } else if (strcmp(operation,"extract") == 0) {
00380 if (strcmp(aperture,"rectangle")==0) {
00381 spec = medianRectangleOfCubeSpectra(cube,llx,lly,urx,ury);
00382 } else if (strcmp(aperture,"circle")==0) {
00383 spec = medianCircleOfCubeSpectra(cube,centerx,centery,radius);
00384 } else {
00385 cpl_msg_error(fctid,"Aperture not supported");
00386 cpl_msg_info(fctid,"Supported apertures are only:");
00387 cpl_msg_info(fctid,"rectangle, circle:");
00388 cpl_propertylist_delete(plist) ;
00389 cpl_frame_delete(product_frame) ;
00390 cpl_image_delete(image) ;
00391 destroy_cube(cube);
00392 destroy_image(im_spec);
00393 fits_header_destroy(head);
00394 return -1;
00395 }
00396 } else {
00397 cpl_msg_error(fctid,"Operation not supported");
00398 cpl_msg_info(fctid,"Supported operations are only:");
00399 cpl_msg_info(fctid,"average, clean_mean, median, sum, extract :");
00400 cpl_propertylist_delete(plist) ;
00401 cpl_frame_delete(product_frame) ;
00402 cpl_image_delete(image) ;
00403 destroy_cube(cube);
00404 destroy_image(im_spec);
00405 fits_header_destroy(head);
00406 return -1;
00407 }
00408
00409 im_spec = vectorToImage(spec);
00410 change_header(head, (char*)name_i);
00411
00412
00413 imagew=cpl_image_wrap_float(im_spec->lx,im_spec->ly,im_spec->data);
00414 image=cpl_image_duplicate(imagew);
00415 cpl_image_unwrap(imagew);
00416
00417
00418
00419 name_o = "out_ima.fits" ;
00420
00421
00422 product_frame = cpl_frame_new();
00423 cpl_frame_set_filename(product_frame, name_o) ;
00424 cpl_frame_set_tag(product_frame, SI_UTL_CUBE2SPECTRUM_PROIMA) ;
00425 cpl_frame_set_type(product_frame, CPL_FRAME_TYPE_IMAGE) ;
00426 cpl_frame_set_group(product_frame, CPL_FRAME_GROUP_PRODUCT) ;
00427 cpl_frame_set_level(product_frame, CPL_FRAME_LEVEL_FINAL) ;
00428 if (cpl_error_get_code()) {
00429 cpl_msg_error(fctid, "Error while initialising the product frame") ;
00430 cpl_propertylist_delete(plist) ;
00431 cpl_frame_delete(product_frame) ;
00432 cpl_image_delete(image) ;
00433 return -1 ;
00434 }
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449 if (cpl_image_save(image, name_o, CPL_BPP_DEFAULT, plist,
00450 CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
00451 cpl_msg_error(fctid, "Could not save product");
00452 cpl_propertylist_delete(plist) ;
00453 cpl_frame_delete(product_frame) ;
00454 cpl_image_delete(image) ;
00455 return -1 ;
00456 }
00457 cpl_propertylist_delete(plist) ;
00458 set_wcs_spectrum (image, (char*)name_o,clam, disp, cpix);
00459
00460
00461 cpl_frameset_insert(framelist, product_frame) ;
00462
00463 destroy_cube(cube);
00464 destroy_image(im_spec);
00465 fits_header_destroy(head);
00466 cpl_image_delete(image) ;
00467
00468
00469 if (cpl_error_get_code())
00470 return -1 ;
00471 else
00472 return 0 ;
00473 }