utl_cube_combine.c

00001 /* $Id: utl_cube_combine.c,v 1.7 2005/10/15 10:19:33 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/10/15 10:19:33 $
00024  * $Revision: 1.7 $
00025  * $Name:  $
00026  */
00027 
00028 
00029 /*-----------------------------------------------------------------------------
00030                                 Includes
00031  -----------------------------------------------------------------------------*/
00032 
00033 
00034 #include "utl_cube_combine.h"
00035 #include "eclipse.h"
00036 #include "utilities.h"
00037 #include <stdio.h>
00038 #include "sinfoni_key_names.h"
00039 
00040 /*-----------------------------------------------------------------------------
00041                             Functions prototypes
00042  -----------------------------------------------------------------------------*/
00043 
00044 
00045 /*-----------------------------------------------------------------------------
00046                             Static variables
00047  -----------------------------------------------------------------------------*/
00048 
00049 
00050 /*-----------------------------------------------------------------------------
00051                                 Functions code
00052  -----------------------------------------------------------------------------*/
00053 /*----------------------------------------------------------------------------*/
00060 /*----------------------------------------------------------------------------*/
00061 int si_utl_cube_combine(
00062         cpl_parameterlist   *   parlist, 
00063         cpl_frameset        *   framelist)
00064 {
00065     const char          *   fctid = "si_utl_cube_combine" ;
00066     cpl_parameter       *   param=NULL ;
00067     const char          *   operation=NULL ;
00068     const char          *   method=NULL ;
00069     double                  sigma=0 ;
00070      cpl_image           *   ima1=NULL ;
00071     cpl_image           *   ima2=NULL ;
00072 
00073     const char                *   name_o=NULL ;
00074     const char                *   name_i=NULL ;
00075     const char                *   name_m=NULL ;
00076    
00077     cpl_propertylist    *   plist=NULL ;
00078     cpl_image           *   image=NULL ;
00079     cpl_frame           *   product_frame=NULL;
00080     cpl_frame           *   frame=NULL;
00081   
00082  
00083     int xsize=0;
00084     int ysize=0;
00085 
00086     int i=0;
00087     int n=0;
00088   
00089     OneCube ** cube_object;
00090     OneCube * cube_jitter=NULL;
00091     OneCube * cube_mask=NULL;
00092     OneCube ** cube_out;
00093     char**    files;
00094     int nframes=0;
00095     float time=0;
00096     float * times=NULL;
00097     float * offsetx=NULL;
00098     float * offsety=NULL;
00099     float * offs_x=NULL;
00100     float * offs_y=NULL;
00101 
00102     char ** fileoffs;
00103     char * filename=NULL;
00104     char * tmpname=NULL;
00105     char * tmptag=NULL;
00106     float tmpoffx=0;
00107     float tmpoffy=0;
00108  
00109     cpl_imagelist * cub_jit=NULL;
00110     cpl_imagelist * cub_msk=NULL;
00111     OneImage* eimg=NULL;
00112     cpl_image * cimg=NULL;
00113 
00114     FILE* file_list=NULL;
00115     int cnt=0;
00116     int np=0;
00117 
00118     /* HOW TO RETRIEVE INPUT PARAMETERS */
00119     /* --stropt */
00120     param = cpl_parameterlist_find(parlist, "sinfoni.si_utl_cube_combine.op");
00121     operation = cpl_parameter_get_string(param);
00122 
00123     param = cpl_parameterlist_find(parlist, "sinfoni.si_utl_cube_combine.method");
00124     method = cpl_parameter_get_string(param);
00125 
00126 
00127     param = cpl_parameterlist_find(parlist, "sinfoni.si_utl_cube_combine.name_o");
00128     name_o = cpl_parameter_get_string(param);
00129 
00130 
00131     param = cpl_parameterlist_find(parlist, "sinfoni.si_utl_cube_combine.name_i");
00132     name_i = cpl_parameter_get_string(param);
00133 
00134 
00135     /* --doubleopt */
00136     param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube_combine.xsize");
00137     xsize = cpl_parameter_get_int(param) ;
00138 
00139 
00140     param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube_combine.ysize");
00141     ysize = cpl_parameter_get_int(param) ;
00142 
00143 
00144     param = cpl_parameterlist_find(parlist,"sinfoni.si_utl_cube_combine.sigma");
00145     sigma = cpl_parameter_get_double(param) ;
00146 
00147 
00148     /* Now performing the data reduction */
00149     /* Let's generate one image for the example */
00150 
00151     if ( NULL == (file_list = fopen (name_i, "r" ) ) )
00152     {
00153       cpl_msg_error(fctid,"cannot open %s\n", name_i) ;
00154       return -1 ;
00155     }
00156     cnt = 0 ;
00157     tmpname = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char));
00158     tmptag = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char));
00159     while ( fscanf( file_list, "%f %f",&tmpoffx, &tmpoffy) != EOF )
00160     {
00161       cnt ++ ;
00162     }
00163 
00164 
00165     fclose(file_list);
00166     n= cnt ;
00167     files = (char**) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00168     offs_x = (float*) cpl_calloc(cnt,sizeof(float*));
00169     offs_y = (float*) cpl_calloc(cnt,sizeof(float*));
00170     cnt=0;   
00171     if ( NULL == (file_list = fopen (name_i, "r" ) ) )
00172     {
00173       cpl_msg_error(fctid,"cannot open %s\n", name_i) ;
00174       return -1 ;
00175     }
00176 
00177     while ( fscanf( file_list, "%f %f",&tmpoffx,&tmpoffy ) != EOF ) {
00178       frame=cpl_frameset_get_frame(framelist,cnt);
00179       files[cnt]=cpl_frame_get_filename(frame);
00180       offs_x[cnt]=tmpoffx;
00181       offs_x[cnt]=tmpoffy;
00182       printf("filename=%s\n",files[cnt]);
00183       cnt ++ ;
00184     }
00185     nframes=cnt;
00186     fclose(file_list);
00187 
00188  
00189     times   = new_float_array(nframes);
00190     offsetx = new_float_array(nframes);
00191     offsety = new_float_array(nframes);
00192 
00193     fileoffs = (char**) cpl_calloc(nframes,sizeof(char*));
00194     filename = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char));
00195     for (n=0;n<nframes;n++){
00196       array_set_value(offsetx, offs_x[n], n);
00197       array_set_value(offsety, offs_y[n], n);
00198     }
00199 
00200 
00201 
00202 
00203     printf("ok1\n");
00204 
00205     cube_object = new_cube_list(nframes);
00206     if(cube_object == NULL) {
00207          cpl_msg_error(fctid,"could not allocate memory") ;
00208      cpl_free(fileoffs);
00209      cpl_free(filename);
00210      destroy_array(times);
00211      destroy_array(offsetx);
00212      destroy_array(offsety);
00213      cpl_free(files);
00214      cpl_free(offs_x);
00215      cpl_free(offs_y);
00216      cpl_free(tmpname);
00217      cpl_free(tmptag);
00218      return -1 ;
00219     }
00220 
00221     cube_out = new_cube_list(nframes);
00222     for (i=0;i<nframes;i++) {
00223       printf("file=%s\n",files[i]);
00224        time = spiffi_get_exptime(files[i]);
00225       array_set_value(times, time, i);
00226       cube_out[i]=load_cube(files[i]); 
00227       insert_cube(cube_object, cube_out[i],  i);
00228     }
00229     np = cube_out[0]->np;
00230 
00231 
00232     for (i=0;i< nframes;i++){
00233       destroy_cube(cube_out[i]);
00234     }
00235     destroy_cube_list(cube_out);
00236         printf("ok2\n");
00237 
00238 
00239     cube_jitter = newCube(xsize, ysize, np);       
00240     if (cube_jitter == NULL){
00241       cpl_msg_error (fctid,"could allocate new data cube");
00242 
00243 
00244       destroy_cube_list(cube_object);
00245       destroy_cube(cube_jitter);
00246       cpl_free(fileoffs);
00247       cpl_free(filename);
00248       destroy_array(times);
00249       destroy_array(offsetx);
00250       destroy_array(offsety);
00251       cpl_free(files);
00252       cpl_free(offs_x);
00253       cpl_free(offs_y);
00254       cpl_free(tmpname);
00255       cpl_free(tmptag);
00256       return -1;
00257     } 
00258        printf("ok3\n");
00259        printf("file0=%s\n",files[0]);
00260     if ((plist=cpl_propertylist_load(files[0],0)) == NULL) {
00261         cpl_msg_error(fctid, "Cannot read the FITS header") ;
00262         return -1 ;
00263     }
00264     if(strcmp(method,"clean_mean")==0) {
00265        cube_mask = combineJitteredCubes(cube_object, cube_jitter, nframes, offsetx, 
00266                                 offsety, times, (char*) "tanh");
00267     } else {
00268        cube_mask = combineCubes(cube_object, cube_jitter, nframes, 
00269                          offsetx, offsety, sigma, (char*) "tanh");
00270     }
00271        printf("ok4\n");
00272     if (cube_mask == NULL){
00273       cpl_msg_error (fctid,"could not merge the jittered data cubes\n");
00274       return -1;
00275     }
00276     name_m="out_cube_mask.fits";
00277        printf("ok5\n");
00278     save_cube_to_fits_hdr_dump( cube_jitter, (char*)name_o, NULL );
00279     save_cube_to_fits_hdr_dump( cube_mask, (char*) name_m, NULL);
00280        printf("ok6\n");
00281 
00282     /* HOW TO SAVE A PRODUCT ON DISK  */
00283 
00284     /* Create product frame */
00285 
00286     product_frame = cpl_frame_new();
00287     cpl_frame_set_filename(product_frame, name_o) ;
00288     cpl_frame_set_tag(product_frame, SI_UTL_CUBE_COMBINE_PROCUBE) ;
00289     cpl_frame_set_type(product_frame, CPL_FRAME_TYPE_IMAGE) ;
00290     cpl_frame_set_group(product_frame, CPL_FRAME_GROUP_PRODUCT) ;
00291     cpl_frame_set_level(product_frame, CPL_FRAME_LEVEL_FINAL) ;
00292     if (cpl_error_get_code()) {
00293         cpl_msg_error(fctid, "Error while initialising the product frame") ;
00294         cpl_propertylist_delete(plist) ;
00295         cpl_frame_delete(product_frame) ;
00296         cpl_image_delete(image) ;
00297         return -1 ;
00298     }
00299 
00300        printf("ok7\n");
00301 
00302     
00303     /* Add DataFlow keywords */
00304     /*
00305     if (cpl_dfs_setup_product_header(plist, product_frame, framelist, parlist,
00306             "si_utl_cube_combine", "SINFONI", KEY_VALUE_HPRO_DID) != CPL_ERROR_NONE) {
00307         cpl_msg_error(fctid, "Problem in the product DFS-compliance") ;
00308         cpl_propertylist_delete(plist) ;
00309         cpl_frame_delete(product_frame) ;
00310         cpl_image_delete(image) ;
00311         return -1 ;
00312     }
00313     */
00314     /* Save the file */
00315 
00316 
00317 
00318     cub_jit=cpl_imagelist_new(cube_jitter->lx,cube_jitter->ly,cube_jitter->np,CPL_TYPE_FLOAT);
00319     for(i=0;i<cube_jitter->np;i++) {
00320       eimg=cube_jitter->plane[i];
00321       cimg=cpl_image_wrap_float(eimg->lx,eimg->ly,eimg->data);
00322       cpl_imagelist_set(cub_jit,cimg,i);
00323     }
00324 
00325 
00326     if (cpl_imagelist_save(cub_jit, name_o, CPL_BPP_DEFAULT, plist,
00327                        CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
00328         cpl_msg_error(fctid, "Could not save product");
00329         cpl_propertylist_delete(plist) ;
00330         cpl_frame_delete(product_frame) ;
00331         cpl_imagelist_delete(cub_jit) ;
00332         return -1 ;
00333     }
00334     cpl_frameset_insert(framelist, product_frame) ;
00335 
00336 
00337 
00338 
00339     product_frame = cpl_frame_new();
00340     cpl_frame_set_filename(product_frame, name_m) ;
00341     cpl_frame_set_tag(product_frame, SI_UTL_CUBE_COMBINE_PROMASK) ;
00342     cpl_frame_set_type(product_frame, CPL_FRAME_TYPE_IMAGE) ;
00343     cpl_frame_set_group(product_frame, CPL_FRAME_GROUP_PRODUCT) ;
00344     cpl_frame_set_level(product_frame, CPL_FRAME_LEVEL_FINAL) ;
00345     if (cpl_error_get_code()) {
00346         cpl_msg_error(fctid, "Error while initialising the product frame") ;
00347         cpl_propertylist_delete(plist) ;
00348         cpl_frame_delete(product_frame) ;
00349         cpl_image_delete(image) ;
00350         return -1 ;
00351     }
00352 
00353     cub_msk=cpl_imagelist_new(cube_mask->lx,cube_mask->ly,cube_mask->np,CPL_TYPE_FLOAT);
00354     for(i=0;i<cube_mask->np;i++) {
00355       eimg=cube_mask->plane[i];
00356       cimg=cpl_image_wrap_float(eimg->lx,eimg->ly,eimg->data);
00357       cpl_imagelist_set(cub_msk,cimg,i);
00358     }
00359 
00360 
00361     if (cpl_imagelist_save(cub_msk, name_m, CPL_BPP_DEFAULT, plist,
00362                        CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
00363         cpl_msg_error(fctid, "Could not save product");
00364         cpl_propertylist_delete(plist) ;
00365         cpl_frame_delete(product_frame) ;
00366         cpl_imagelist_delete(cub_msk) ;
00367         return -1 ;
00368     }
00369     cpl_frameset_insert(framelist, product_frame) ;
00370 
00371 
00372 
00373     cpl_propertylist_delete(plist) ;
00374     cpl_image_delete(ima1);
00375     cpl_image_delete(ima2);
00376     cpl_image_delete(image) ;
00377 
00378 
00379 
00380     destroy_array(times);
00381     destroy_array(offsetx);
00382     destroy_array(offsety);
00383     destroy_cube(cube_mask);
00384     destroy_cube(cube_jitter);
00385     destroy_cube_list(cube_object);
00386 
00387     
00388 
00389 
00390     cpl_free(fileoffs);
00391     cpl_free(filename);
00392     cpl_free(offs_x);
00393     cpl_free(offs_y);
00394     /* Log the saved file in the input frameset */
00395     
00396     /* Return */
00397     if (cpl_error_get_code()) {
00398         return -1 ;
00399     }
00400     else { 
00401         return 0 ;
00402     }
00403 }

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