sinfoni/sinfo_utl_cube2spectrum.c

00001 /* $Id: sinfo_utl_cube2spectrum.c,v 1.6 2006/12/06 08:46:17 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: 2006/12/06 08:46:17 $
00024  * $Revision: 1.6 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #  include <config.h>
00030 #endif
00031 
00032 /*----------------------------------------------------------------------------
00033                                 Includes
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                             Functions prototypes
00048  ----------------------------------------------------------------------------*/
00049 
00050 /*----------------------------------------------------------------------------
00051                             Static variables
00052  ----------------------------------------------------------------------------*/
00053 
00054 static void 
00055 sinfo_change_header(const char * name);
00063 /*-----------------------------------------------------------------------------
00064                             Function implementation
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     /*fits_header* head=NULL; */
00166 
00167      /* HOW TO RETRIEVE INPUT PARAMETERS */
00168     /* --stropt */
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     /* --stropt */
00174     check_nomsg(param=cpl_parameterlist_find(parlist, 
00175                                          "sinfoni.sinfo_utl_cube2spectrum.ap"));
00176     check_nomsg(aperture=cpl_parameter_get_string(param));
00177      /* --stropt */
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     /* --intopt */
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     /* --intopt */
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     /* --intopt */
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     /* --intopt */
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     /* --intopt */
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     /* --intopt */
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     /* --intopt */
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     /* --intopt */
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     /* --intopt */
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     /* Identify the RAW and CALIB frames in the input frameset */
00239     if (sinfo_dfs_set_groups(framelist)) {
00240         sinfo_msg_error( "Cannot identify RAW and CALIB frames") ;
00241         goto cleanup;
00242     }
00243 
00244     /* HOW TO ACCESS INPUT DATA */
00245     cknull(frm_cub=cpl_frameset_find(framelist,tag_i),
00246       "SOF does not have a file tagged as %s",tag_i);
00247    
00248     /* Now performing the data reduction */
00249     /* Let's generate one image for the example */
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     /* head = sinfo_fits_read_header((char*)name_i); */
00392     check_nomsg(sinfo_change_header(name_i));
00393 
00394     /* HOW TO SAVE A PRODUCT ON DISK  */
00395     /* Set the file name */
00396  
00397     /* Create product frame */
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     /* This generate seg fault
00416     if(spec != NULL) sinfo_free_svector(&spec);
00417     */
00418     sinfo_free_table(&tbl_spec);
00419     /*if(head != NULL) sinfo_fits_header_destroy(head); */
00420  
00421     /* Return */
00422     if (cpl_error_get_code()) 
00423         return -1 ;
00424     else 
00425         return 0 ;
00426 
00427 
00428 }

Generated on Wed Jan 17 08:33:44 2007 for SINFONI Pipeline Reference Manual by  doxygen 1.4.4