sinfo_utilities.c

00001 /*
00002  * This file is part of the ESO SINFONI Pipeline
00003  * Copyright (C) 2004,2005 European Southern Observatory
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or
00008  * (at your option) any later version.
00009  *
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00018  */
00019 #ifdef HAVE_CONFIG_H
00020 #  include <config.h>
00021 #endif
00022 #include "sinfo_utilities.h"
00023 #include "sinfo_error.h"
00024 #include "sinfo_utils_wrappers.h"
00025 #include "sinfo_function_1d.h"
00037 int
00038 sinfo_table_column_dump(cpl_table* t, const char* name, cpl_type type)
00039 {
00040   int nrow=0;
00041   int i=0;
00042   int* pi=NULL;
00043   float* pf=NULL;
00044   double* pd=NULL;
00045   char** ps=NULL;
00046   int index=0;
00047 
00048   nrow=cpl_table_get_nrow(t);
00049 
00050   switch(type) {
00051 
00052   case CPL_TYPE_INT:
00053     pi=cpl_table_get_data_int(t,name);
00054     for(i=0;i<nrow;i++) {
00055       sinfo_msg("val=%d",pi[i]);
00056     }
00057     break;
00058   case CPL_TYPE_FLOAT:
00059     pf=cpl_table_get_data_float(t,name);
00060     for(i=0;i<nrow;i++) {
00061       sinfo_msg("val=%g",pf[i]);
00062     }
00063     break;
00064   case CPL_TYPE_DOUBLE:
00065     pd=cpl_table_get_data_double(t,name);
00066     for(i=0;i<nrow;i++) {
00067       sinfo_msg("val=%g",pd[i]);
00068     }
00069     break;
00070   case CPL_TYPE_STRING:
00071     ps=cpl_table_get_data_string(t,name);
00072     for(i=0;i<nrow;i++) {
00073       sinfo_msg("val=%s",ps[i]);
00074     }
00075     break;
00076   default:
00077     sinfo_msg_error("Wrong column type");
00078     cpl_error_set(cpl_func, CPL_ERROR_TYPE_MISMATCH);
00079     return index;
00080 
00081   }
00082   return 0;
00083 }
00084 
00085 
00086 
00096 cpl_table*
00097 sinfo_table_shift_column_spline3(cpl_table* t, 
00098                                  const char* col, 
00099                                  const double shift)
00100 {
00101   cpl_table* out=NULL;
00102   int nrow=0;
00103   int i=0;
00104   int z=0;
00105 
00106   float sum=0;
00107   float new_sum=0;
00108 
00109   float* pi=NULL;
00110   float* po=NULL;
00111   float* eval=NULL;
00112   float* xnum=NULL;
00113   float* spec=NULL;
00114   float* corrected_spec=NULL;
00115 
00116   cknull(t,"null input table");
00117   out=cpl_table_duplicate(t);
00118 
00119   nrow=cpl_table_get_nrow(t);
00120   check_nomsg(cpl_table_cast_column(t,col,"FINT",CPL_TYPE_FLOAT));
00121   check_nomsg(cpl_table_cast_column(out,col,"FINT",CPL_TYPE_FLOAT));
00122   pi=cpl_table_get_data_float(t,"FINT");
00123   po=cpl_table_get_data_float(out,"FINT");
00124   
00125 
00126 
00127   xnum=cpl_calloc(nrow,sizeof(float)) ;
00128   /* fill the xa[] array for the spline function */
00129   for ( i = 0 ; i < nrow ; i++ ) {
00130     xnum[i] = i ;
00131   }
00132 
00133   spec=cpl_calloc(nrow,sizeof(float)) ;
00134   corrected_spec=cpl_calloc(nrow,sizeof(float)) ;
00135   eval=cpl_calloc(nrow,sizeof(float)) ;
00136 
00137   sum = 0. ;
00138   for ( z = 0 ; z < nrow ; z++ ) {
00139     spec[z] = pi[z] ;
00140     if (isnan(spec[z]) ) {
00141       for ( i = z-1 ; i <= z+1 ; i++ ) {
00142     if ( i < 0 ) continue ;
00143     if ( i >= nrow) continue ;
00144     corrected_spec[i] = ZERO ;
00145       }
00146       spec[z] = 0. ;
00147     }
00148     sum += spec[z] ;
00149     eval[z] = (float)shift+(float)z ;
00150   }
00151   /* now we do the spline interpolation*/
00152   if ( -1 == sinfo_function1d_natural_spline(xnum,spec, nrow, 
00153                                              eval,corrected_spec, nrow))
00154     {
00155       sinfo_msg_error("error in spline interpolation!") ;
00156       goto cleanup;
00157     }
00158         
00159   new_sum = 0. ;
00160   for ( z = 0 ; z < nrow ; z++ ) {
00161     if ( isnan(corrected_spec[z]) ) {
00162       continue ;
00163     }
00164     new_sum += corrected_spec[z] ;
00165   }
00166   /* fill output imagelist */
00167   for ( z = 0 ; z < nrow ; z++ ) {
00168     if ( new_sum == 0. ) new_sum =1. ; 
00169     {
00170       if ( isnan(corrected_spec[z]) ) {
00171     po[z] = ZERO ;
00172       } else {
00173     corrected_spec[z] *= sum / new_sum ;
00174     po[z] = corrected_spec[z] ;
00175       }
00176     }
00177   }
00178 
00179   sinfo_free_float(&xnum);
00180   sinfo_free_float(&spec) ;
00181   sinfo_free_float(&corrected_spec) ;
00182   sinfo_free_float(&eval) ;
00183 
00184   check_nomsg(cpl_table_erase_column(t,"FINT"));
00185   check_nomsg(cpl_table_erase_column(out,col));
00186   check_nomsg(cpl_table_cast_column(out,"FINT",col,CPL_TYPE_DOUBLE));
00187   check_nomsg(cpl_table_erase_column(out,"FINT"));
00188 
00189   return out;
00190  cleanup:
00191 
00192   sinfo_free_float(&xnum);
00193   sinfo_free_float(&spec) ;
00194   sinfo_free_float(&corrected_spec) ;
00195   sinfo_free_float(&eval) ;
00196   sinfo_free_table(&out);
00197   return NULL;
00198 
00199 
00200 }
00201 
00202 
00212 cpl_table*
00213 sinfo_table_shift_column_int(const cpl_table* t, 
00214                              const char* col, 
00215                              const double s, 
00216                                    double* r)
00217 {
00218   cpl_table* out=NULL;
00219   int is=(int)s;
00220   int nrow=0;
00221   int i=0;
00222 
00223   double* pi=NULL;
00224   double* po=NULL;
00225 
00226   cknull(t,"null input table");
00227   out=cpl_table_duplicate(t);
00228   *r=s-is;
00229   nrow=cpl_table_get_nrow(t);
00230   pi=cpl_table_get_data_double(t,col);
00231   po=cpl_table_get_data_double(out,col);
00232   for(i=0;i<nrow;i++) {
00233     if( ((i-is) >=0) && ((i-is) < nrow)) {
00234       po[i-is]=pi[i];
00235     }
00236   }
00237   return out;
00238  cleanup:
00239   sinfo_free_table(&out);
00240   return NULL;
00241 
00242 }
00243 
00244 
00255 cpl_table*
00256 sinfo_table_shift_column_poly(cpl_table* t, 
00257                               const char* col, 
00258                               const double shift,
00259                               const int order)
00260 {
00261   cpl_table* out=NULL;
00262   int nrow=0;
00263   int i=0;
00264   int flag=0;
00265   int n_points=0;
00266   int firstpos=0;
00267   int z=0;
00268   float eval=0;
00269   float sum=0;
00270   float new_sum=0;
00271   float* pi=NULL;
00272   float* po=NULL;
00273   float* spec=NULL ;
00274   float* corrected_spec=NULL ;
00275   float* xnum=NULL ;
00276   float* tableptr=NULL;
00277 
00278   cknull(t,"null input table");
00279   if ( order <= 0 ) {
00280     sinfo_msg_error("wrong order of interpolation polynom given!") ;
00281     goto cleanup;
00282   }
00283 
00284   out=cpl_table_duplicate(t);
00285 
00286   nrow=cpl_table_get_nrow(t);
00287   cpl_table_cast_column(t,col,"FINT",CPL_TYPE_FLOAT);
00288   cpl_table_cast_column(out,col,"FINT",CPL_TYPE_FLOAT);
00289   pi=cpl_table_get_data_float(t,"FINT");
00290   po=cpl_table_get_data_float(out,"FINT");
00291 
00292   n_points = order + 1 ;
00293   if ( n_points % 2 == 0 ) {
00294     firstpos = (int)(n_points/2) - 1 ;
00295   } else {
00296     firstpos = (int)(n_points/2) ;
00297   }
00298   spec=cpl_calloc(nrow,sizeof(float)) ;
00299   corrected_spec=cpl_calloc(nrow,sizeof(float)) ;
00300   xnum=cpl_calloc(order+1,sizeof(float)) ;
00301   /* fill the xa[] array for the polint function */
00302   for ( i = 0 ; i < n_points ; i++ ) {
00303     xnum[i] = i ;
00304   }
00305 
00306 
00307   for(i=0;i<nrow;i++) {
00308     corrected_spec[i] = 0. ;
00309   }
00310 
00311   sum = 0. ;
00312   for ( z = 0 ; z < nrow ; z++ ) {
00313     spec[z] = pi[z] ;
00314     if (isnan(spec[z]) ) {
00315       spec[z] = 0. ;
00316                   
00317       for ( i = z - firstpos ; i < z-firstpos+n_points ; i++ ) {
00318     if ( i < 0 ) continue ;
00319     if ( i >= nrow) continue  ;
00320     corrected_spec[i] = ZERO ;
00321       }
00322     }
00323     if ( z != 0 && z != nrow - 1 ) {
00324       sum += spec[z] ;
00325     }
00326   }
00327 
00328   new_sum = 0. ;
00329   for ( z = 0 ; z < nrow ; z++ ) {
00330     /* ---------------------------------------------------------------
00331      * now determine the arrays of size n_points with which the
00332      * polynom is determined and determine the position eval
00333      * where the polynom is evaluated in polynomial interpolation.
00334      * Take care of the points near the row edges!
00335      */
00336     if (isnan(corrected_spec[z])) continue ;
00337     if ( z - firstpos < 0 ) {
00338       tableptr = &spec[0] ;
00339       eval     = shift + z ;
00340     } else if ( z - firstpos + n_points >= nrow ) {
00341       tableptr = &spec[nrow - n_points] ;
00342       eval     = shift + z + n_points - nrow ;
00343     } else {
00344       tableptr = &spec[z-firstpos] ;
00345       eval     = shift + firstpos ;
00346     }
00347 
00348     flag=0;
00349     corrected_spec[z]=sinfo_new_nev_ille(xnum,tableptr,order,eval,&flag);
00350     if ( z != 0 && z != nrow - 1 ) {
00351       new_sum += corrected_spec[z] ;
00352     }
00353   }
00354 
00355   /* fill the output spectrum */
00356   for (z = 0 ; z < nrow ; z++ ) {
00357     if ( new_sum == 0. ) {
00358       new_sum = 1. ;
00359     }
00360     if ( z == 0 ) {
00361       po[z] = ZERO ;
00362     } else if ( z == nrow - 1 ) {
00363       po[z] = ZERO ;
00364     } else if ( isnan(corrected_spec[z]) ) {
00365       po[z] = ZERO ;
00366     } else {
00367       corrected_spec[z] *= sum / new_sum ;
00368       po[z] = corrected_spec[z] ;
00369     }
00370   }
00371   check_nomsg(cpl_table_erase_column(t,"FINT"));
00372   check_nomsg(cpl_table_erase_column(out,col));
00373   check_nomsg(cpl_table_cast_column(out,"FINT",col,CPL_TYPE_DOUBLE));
00374   check_nomsg(cpl_table_erase_column(out,"FINT"));
00375 
00376   sinfo_free_float(&spec) ;
00377   sinfo_free_float(&corrected_spec) ;
00378   sinfo_free_float(&xnum) ;
00379 
00380   return out;
00381  cleanup:
00382 
00383 
00384   sinfo_free_float(&spec) ;
00385   sinfo_free_float(&corrected_spec) ;
00386   sinfo_free_float(&xnum) ;
00387   sinfo_free_table(&out);
00388   return NULL;
00389 
00390 
00391 }
00392 
00393 
00394 
00395 
00396     void sinfo_new_array_set_value( float * array, float value, int i )
00397     {
00398         array[i] = value ;
00399     }
00400     float sinfo_new_array_get_value( float * array, int i )
00401     {
00402         return array[i] ;
00403     }
00404 
00405  
00406 
00407     void sinfo_new_destroy_array(float ** array)
00408     {
00409       if(*array != NULL) {
00410         cpl_free( *array ) ;
00411         *array = NULL;
00412       }
00413     }
00414 
00415     void sinfo_new_destroy_stringarray(char ** array, int size_x)
00416     {
00417         int i ;
00418       
00419         for ( i = 0 ; i < size_x ; i++ )
00420         {
00421             cpl_free( array[i] ) ;
00422         }
00423         cpl_free( array ) ;
00424     }
00425 
00426     void sinfo_new_destroy_2Dintarray(int *** array, int size_x)
00427     {
00428         int i ;
00429    
00430         if((*array) != NULL) {
00431       for ( i = 0 ; i < size_x ; i++ ) {
00432         if((*array)[i] != NULL) {
00433           cpl_free( (*array)[i] );
00434           (*array)[i]=NULL;
00435         }
00436       }
00437       cpl_free( *array ) ;
00438       *array=NULL;
00439     }
00440 
00441     }
00442 
00443 
00444     void sinfo_new_intarray_set_value( int * array, int value, int i )
00445     {
00446         array[i] = value ;
00447     }
00448     float sinfo_new_array2D_get_value( float ** array, int x, int y )
00449     {
00450         return array[x][y] ;
00451     }
00452     int ** sinfo_new_2Dintarray( int size_x, int size_y)
00453     {
00454         int ** retVal ;
00455         int i ;
00456 
00457         retVal = (int **) cpl_calloc( size_x, sizeof (int*) ) ;
00458         for ( i = 0 ; i < size_x ; i++ )
00459         {
00460             retVal[i] = (int *) cpl_calloc( size_y, sizeof (int)) ;
00461         }
00462         return retVal ;
00463     }
00464 
00465     float * sinfo_new_floatarray( int size)
00466     {
00467         return (float *) cpl_calloc( size, sizeof (float) ) ;
00468     }
00469 
00470 
00471     void sinfo_new_destroy_2Dfloatarray(float *** array, int size_x)
00472     {
00473         int i ;
00474         if((*array) != NULL) {
00475       for ( i = 0 ; i < size_x ; i++ ) {
00476         if((*array)[i] != NULL) {
00477           cpl_free( (*array)[i] );
00478           (*array)[i]=NULL;
00479         }
00480       }
00481       cpl_free( *array ) ;
00482       *array=NULL;
00483     }
00484     }
00485 
00486     void sinfo_new_destroy_2Ddoublearray(double *** array, int size_x)
00487     {
00488         int i ;
00489 
00490         if((*array) != NULL) {
00491       for ( i = 0 ; i < size_x ; i++ ) {
00492         if((*array)[i] != NULL) {
00493           cpl_free( (*array)[i] );
00494           (*array)[i]=NULL;
00495         }
00496       }
00497       cpl_free( *array ) ;
00498       *array=NULL;
00499     }
00500 
00501     }
00502 
00503 
00504     void sinfo_new_array2D_set_value( float ** array, float value, int x, int y )
00505     {
00506         array[x][y] = value ;
00507     }
00508 
00509     double sinfo_new_doublearray_get_value( double * array, int i )
00510     {
00511         return array[i] ;
00512     }
00513     void sinfo_new_doublearray_set_value( double * array, double value, int i )
00514     {
00515         array[i] = value ;
00516     }
00517 
00518    void sinfo_new_destroy_doublearray(double * array)
00519     {
00520         cpl_free( array ) ;
00521     }
00522     double * sinfo_new_doublearray( int size)
00523     {
00524         return (double *) cpl_calloc( size, sizeof (double) ) ;
00525     }
00526 
00527     double ** sinfo_new_2Ddoublearray( int size_x, int size_y)
00528     {
00529         double ** retVal ;
00530         int i ;
00531 
00532         retVal = (double **) cpl_calloc( size_x, sizeof (double*) ) ;
00533         for ( i = 0 ; i < size_x ; i++ )
00534         {
00535             retVal[i] = (double *) cpl_calloc( size_y, sizeof (double)) ;
00536         }
00537         return retVal ;
00538     }
00539 
00540     float ** sinfo_new_2Dfloatarray( int size_x, int size_y)
00541     {
00542         float ** retVal ;
00543         int i ;
00544 
00545         retVal = (float **) cpl_calloc( size_x, sizeof (float*) ) ;
00546         for ( i = 0 ; i < size_x ; i++ )
00547         {
00548             retVal[i] = (float *) cpl_calloc( size_y, sizeof (float)) ;
00549         }
00550         return retVal ;
00551     }
00552 
00553 
00554     int * sinfo_new_intarray( int size)
00555     {
00556         return (int *) cpl_calloc( size, sizeof (int) ) ;
00557     }
00558     void sinfo_new_destroy_intarray(int ** array)
00559     {
00560         cpl_free( *array ) ;
00561         *array=NULL;
00562     }
00563 
00564     int sinfo_new_intarray_get_value( int * array, int i )
00565     {
00566         return array[i] ;
00567     }
00568 
00569     float sinfo_new_Stats_get_cleanstdev(Stats * stats)
00570     {
00571         return stats -> cleanstdev ;
00572     }
00573     float sinfo_new_Stats_get_cleanmean(Stats * stats)
00574     {
00575         return stats -> cleanmean ;
00576     }
00577 
00578 char * sinfo_new_get_basename(const char *filename)
00579 {
00580   char *p ;
00581   p = strrchr (filename, '/');
00582   return p ? p + 1 : (char *) filename;
00583 }
00584 
00585 
00586 
00587 char * sinfo_new_get_rootname(const char * filename)
00588 {
00589     static char path[MAX_NAME_SIZE+1];
00590     char * lastdot ;
00591 
00592     if (strlen(filename)>MAX_NAME_SIZE) return NULL ;
00593     memset(path, MAX_NAME_SIZE, 0);
00594     strcpy(path, filename);
00595     lastdot = strrchr(path, '.');
00596     if (lastdot == NULL) return path ;
00597     if ((!strcmp(lastdot, ".fits")) || (!strcmp(lastdot, ".FITS")) ||
00598         (!strcmp(lastdot, ".paf")) || (!strcmp(lastdot, ".PAF")) ||
00599         (!strcmp(lastdot, ".dat")) || (!strcmp(lastdot, ".DAT")) ||
00600         (!strcmp(lastdot, ".fits")) || (!strcmp(lastdot, ".TFITS")) ||
00601         (!strcmp(lastdot, ".ascii")) || (!strcmp(lastdot, ".ASCII")))
00602     {
00603         lastdot[0] = (char)0;
00604     }
00605     return path ;
00606 }
00607 
00608 
00609 
00610 /*----------------------------------------------------------------------------*/
00618 /*----------------------------------------------------------------------------*/
00619 cpl_imagelist * sinfo_new_frameset_to_iset(cpl_frameset * fset)
00620 {
00621     cpl_imagelist   *   iset=NULL ;
00622     char        **  filenames ;
00623     int             nfiles=0 ;
00624 
00625     /* Test entries */
00626     if (fset == NULL) return NULL ;
00627 
00628     /* Get the filenames */
00629     if ((filenames = sinfo_new_frameset_to_filenames(fset, &nfiles)) == NULL) {
00630         sinfo_msg_error( "Cannot get the files names") ;
00631         return NULL ;
00632     }
00633     /* Load image set */
00634     if ((iset = sinfo_new_imagelist_load_frameset(fset, 
00635                                          CPL_TYPE_FLOAT, 0, 0)) == NULL) {
00636         sinfo_msg_error( "Cannot load *** the image set") ;
00637         sinfo_msg_error((char* ) cpl_error_get_message());
00638 
00639         cpl_free(filenames) ;
00640         return NULL ;
00641     }
00642 
00643     /* Free and Return  */
00644     cpl_free(filenames) ;
00645     return iset ;
00646 }
00647 #include "cpl_imagelist_io.h"
00648 
00649 /*----------------------------------------------------------------------------*/
00660 /*----------------------------------------------------------------------------*/
00661 cpl_imagelist * 
00662 sinfo_new_imagelist_load_frameset(const cpl_frameset * frameset,
00663                                   cpl_type type, 
00664                                   int pnum,
00665                                   int extnum)
00666 {
00667     cpl_image  * image     = NULL;
00668     cpl_imagelist  * imagelist = NULL;
00669     cpl_frame  * frame     = cpl_frameset_get_first(frameset);
00670     const int nz = cpl_frameset_get_size(frameset);
00671     int       i;
00672 
00673     /* Require imagelist to contain at least one image */
00674     cpl_ensure(nz > 0, CPL_ERROR_DATA_NOT_FOUND, NULL);
00675 
00676     for (i = 0; frame != NULL; i++, frame = cpl_frameset_get_next(frameset)) {
00677 
00678         const char * name = cpl_frame_get_filename(frame);
00679         if (name == NULL) break; /* Error check */
00680 
00681 
00682         image = cpl_image_load(name, type, pnum, extnum);
00683 
00684         if (image == NULL) break; /* Error check */
00685 
00686         if (i == 0) {
00687             const int nx = cpl_image_get_size_x(image);
00688             const int ny = cpl_image_get_size_y(image);
00689 
00690             if (nx < 1 || ny < 1) break; /* Error check */
00691             imagelist = cpl_imagelist_new();
00692         if (imagelist == NULL) break; /* Error check */
00693         }
00694 
00695         if (cpl_imagelist_set(imagelist, image, i)) break;
00696         image = NULL; /* Image is now part of the imagelist */
00697 
00698     }
00699 
00700     if (i != nz) {
00701         /* Error handling */
00702         cpl_image_delete(image);
00703         cpl_imagelist_delete(imagelist);
00704         imagelist = NULL;
00705     }
00706     return imagelist;
00707 
00708 }
00709 
00720 char ** sinfo_new_frameset_to_filenames(cpl_frameset *set, int *nfiles)
00721 {
00722     char **filenames=NULL;
00723 
00724     int nbframes=0;
00725     int i=0;
00726 
00727     cpl_frame *curr_frame;
00728 
00729     if (set == NULL) {
00730         return NULL;
00731     }
00732 
00733     nbframes = cpl_frameset_get_size(set);
00734 
00735     if (nbframes < 1) {
00736         return NULL;
00737     }
00738     /*
00739      * Create the list of filenames and fill it
00740      */
00741     filenames = cpl_malloc(nbframes * sizeof(char *));
00742 
00743     curr_frame = cpl_frameset_get_first(set);
00744     for (i = 0; i < nbframes; i++) {
00745        filenames[i]=(char*) cpl_frame_get_filename(curr_frame);
00746         curr_frame = cpl_frameset_get_next(set);
00747     }
00748 
00749 
00750     /*
00751      * Set the number of files found
00752      */
00753     *nfiles = nbframes;
00754 
00755     return filenames;
00756 
00757 }
00758 
00759 

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