utl_cube2spectrum.c

00001 /* $Id: utl_cube2spectrum.c,v 1.7 2005/09/05 14:25:31 amodigli Exp $
00002  *
00003  * This file is part of the IIINSTRUMENT Pipeline
00004  * Copyright (C) 2002,2003 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00019  */
00020 
00021 /*
00022  * $Author: amodigli $
00023  * $Date: 2005/09/05 14:25:31 $
00024  * $Revision: 1.7 $
00025  * $Name:  $
00026  */
00027 
00028 
00029 /*-----------------------------------------------------------------------------
00030                                 Includes
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                             Functions prototypes
00043  -----------------------------------------------------------------------------*/
00044 
00045 /*-----------------------------------------------------------------------------
00046                             Static variables
00047  -----------------------------------------------------------------------------*/
00048 
00049 /*-----------------------------------------------------------------------------
00050                                 Functions code
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     /* HOW TO RETRIEVE INPUT PARAMETERS */
00136     /* --stropt */
00137     param = cpl_parameterlist_find(parlist, "sinfoni.si_utl_cube2spectrum.op");
00138     operation = cpl_parameter_get_string(param);
00139  
00140     /* --stropt */
00141     param = cpl_parameterlist_find(parlist, "sinfoni.si_utl_cube2spectrum.ap");
00142     aperture = cpl_parameter_get_string(param);
00143  
00144     /* --stropt */
00145     name_o = "out_spec.fits";
00146  
00147     /* --intopt */
00148     param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.llx");
00149     llx = cpl_parameter_get_int(param) ;
00150 
00151     /* --intopt */
00152     param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.lly");
00153     lly = cpl_parameter_get_int(param) ;
00154 
00155     /* --intopt */
00156     param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.urx");
00157     urx = cpl_parameter_get_int(param) ;
00158 
00159     /* --intopt */
00160     param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.ury");
00161     ury = cpl_parameter_get_int(param) ;
00162 
00163     /* --intopt */
00164     param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.lo_rej");
00165     lo_rej = cpl_parameter_get_int(param) ;
00166 
00167     /* --intopt */
00168     param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.hi_rej");
00169     hi_rej = cpl_parameter_get_int(param) ;
00170 
00171     /* --intopt */
00172     param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.centerx");
00173     centerx = cpl_parameter_get_int(param) ;
00174 
00175     /* --intopt */
00176     param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.centery");
00177     centery = cpl_parameter_get_int(param) ;
00178 
00179     /* --intopt */
00180     param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube2spectrum.radius");
00181     radius = cpl_parameter_get_int(param) ;
00182   
00183 
00184     /* Identify the RAW and CALIB frames in the input frameset */
00185     if (sinfoni_dfs_set_groups(framelist)) {
00186         cpl_msg_error(fctid, "Cannot identify RAW and CALIB frames") ;
00187         return -1 ;
00188     }
00189  
00190     /* HOW TO ACCESS INPUT DATA */
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     /* Now performing the data reduction */
00203     /* Let's generate one image for the example */
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     /* save_image_to_fits_hdr_dump ( im_spec, (char*)name_o, head, -32 ); */ 
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     /* HOW TO SAVE A PRODUCT ON DISK  */
00418     /* Set the file name */
00419     name_o = "out_ima.fits" ;
00420 
00421     /* Create product frame */
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     /* Add DataFlow keywords */
00437     /*
00438     if (cpl_dfs_setup_product_header(plist, product_frame, framelist, parlist,
00439             "si_utl_cube2spectrum", "SINFONI", KEY_VALUE_HPRO_DID) != CPL_ERROR_NONE) {
00440         cpl_msg_error(fctid, "Problem in the product DFS-compliance") ;
00441         cpl_propertylist_delete(plist) ;
00442         cpl_frame_delete(product_frame) ;
00443         cpl_image_delete(image) ;
00444         return -1 ;
00445     }
00446     */
00447 
00448     /* Save the file */
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     /* Log the saved file in the input frameset */
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     /* Return */
00469     if (cpl_error_get_code()) 
00470         return -1 ;
00471     else 
00472         return 0 ;
00473 }

Generated on Wed Oct 26 13:08:56 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001