sinfo_new_cube_ops.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 /*************************************************************************
00020 * E.S.O. - VLT project
00021 *
00022 *
00023 *
00024 * who       when      what
00025 * --------  --------  ----------------------------------------------
00026 * schreib  17/05/00  created
00027 */
00028 
00029 /************************************************************************
00030 *   NAME
00031 *       sinfo_new_cube_ops.c -
00032 *       cube arithmetic routines
00033 *
00034 *   SYNOPSIS
00035 *    #include "sinfo_new_cube_ops.h"
00036 *
00037 *
00038 *
00039 *    2) cpl_imagelist *
00040 *       sinfo_new_cube_ops( cpl_imagelist   *   cube1,
00041 *                cpl_imagelist  *   cube2,
00042 *                int        operation) 
00043 *
00044 *    3) cpl_imagelist *
00045 *       sinfo_new_cube_const_ops(
00046 *                     cpl_imagelist    * cube1,
00047 *                     double    constant,
00048 *                     int       operation)
00049 *
00050 *    4) cpl_imagelist *
00051 *       sinfo_new_cube_sub(
00052 *               cpl_imagelist   *   c1,
00053 *               cpl_imagelist   *   c2 )
00054 *
00055 *    5) cpl_imagelist *
00056 *       sinfo_new_cube_add(
00057 *               cpl_imagelist   *   c1,
00058 *               cpl_imagelist   *   c2  )
00059 *    6) cpl_imagelist *
00060 *       sinfo_new_cube_mul(
00061 *               cpl_imagelist   *   c1,
00062 *               cpl_imagelist   *   c2 )
00063 *
00064 *    7) cpl_imagelist *
00065 *       sinfo_new_cube_div(
00066 *               cpl_imagelist   *   c1,
00067 *               cpl_imagelist   *   c2 )
00068 *
00069 *    8) cpl_imagelist * sinfo_new_add_image_to_cube(cpl_imagelist * cu, 
00070                                                     cpl_image * im)
00071 *
00072 *    9) cpl_imagelist * sinfo_new_sub_image_from_cube (cpl_imagelist * cu, 
00073                                                     cpl_image * im)
00074 *
00075 *    10) cpl_imagelist * sinfo_new_mul_image_to_cube(cpl_imagelist * cu, 
00076                                                     cpl_image * im)
00077 *
00078 *    11) cpl_imagelist * sinfo_new_div_cube_by_image(cpl_imagelist * cu, 
00079                                                     cpl_image * im)
00080 * 
00081 *    12) cpl_imagelist * sinfo_new_add_spectrum_to_cube(cpl_imagelist *cu, 
00082                                                     Vector *spec)
00083 * 
00084 *    13) cpl_imagelist * sinfo_new_sub_spectrum_from_cube(cpl_imagelist *cu, 
00085                                                     Vector *spec)
00086 *
00087 *    14) cpl_imagelist * sinfo_new_mul_spectrum_to_cube(cpl_imagelist *cu, 
00088                                                     Vector *spec)
00089 *
00090 *    15) cpl_imagelist * sinfo_new_div_cube_by_spectrum(cpl_imagelist *cu, 
00091                                                     Vector *spec)
00092 *
00093 *    16) Vector * sinfo_new_clean_mean_of_spectra(cpl_imagelist * cube, 
00094 *                                    int llx, 
00095 *                                    int lly, 
00096 *                                    int urx, 
00097 *                                    int ury,
00098 *                                    double lo_reject,
00099 *                                    double hi_reject)
00100 *
00101 *    17) cpl_image * sinfo_new_median_cube(cpl_imagelist * cube) 
00102 *
00103 *    18) cpl_image * sinfo_new_average_cube_to_image(cpl_imagelist * cube) 
00104 *
00105 *    19) cpl_image * sinfo_new_sum_cube_to_image(cpl_imagelist * cube)
00106 *
00107 *    20) cpl_image * 
00108          sinfo_new_average_cube_to_image_between_waves (cpl_imagelist * cube, 
00109 *                                                   float     dispersion,
00110 *                                                   float     centralWave,
00111 *                                                   float     initialLambda,
00112 *                                                   float     finalLambda)  
00113 *
00114 *    21) cpl_image * sinfo_new_extract_image_from_cube(cpl_imagelist * cube, 
00115                                                     int plane_index) 
00116 *
00117 *    22) Vector * sinfo_new_extract_spectrum_from_cube( cpl_imagelist * cube, 
00118                                                     int x_pos, int y_pos )
00119 *    23) cpl_imagelist * 
00120          sinfo_new_combine_jittered_cubes ( cpl_imagelist ** cubes,
00121 *                                         cpl_imagelist  * mergedCube,
00122 *                                         int        n_cubes,
00123 *                                         float    * cumoffsetx,
00124 *                                         float    * cumoffsety,
00125 *                                         float    * exptimes,
00126 *                                         char     * kernel_type )
00127 *    24) cpl_imagelist * sinfo_new_interpol_cube_simple( cpl_imagelist * cube, 
00128 *                                      cpl_imagelist * badcube,
00129 *                                      int       maxdist )
00130 *
00131 *
00132 *    25) cpl_imagelist * sinfo_cube_zshift(const cpl_imagelist * cube, 
00133 *                                          const double shift,
00134 *                                          double* rest)
00135 *
00136 *    26) cpl_imagelist * sinfo_cube_zshift_poly(const cpl_imagelist * cube, 
00137 *                                               const double shift,
00138 *                                               const int    order)
00139 *
00140 *    27) cpl_imagelist * sinfo_cube_zshift_spline3(const cpl_imagelist * cube, 
00141 *                                                  const double shift)
00142 *
00143 *
00144 *
00145 *
00146 *   DESCRIPTION
00147 *    2) 4 operations between 2 cubes
00148 *    3) 4 operations between a cube and a constant
00149 *    4) subtract one cube from another
00150 *    5) add a cube to another
00151 *    6) multiply two cubes
00152 *    7) divide two cubes
00153 *    8) add an image to all planes in the cube
00154 *    9) subtract an image from all planes in the cube
00155 *    10) multiply an image to all planes in the cube
00156 *    11) divide all planes in the cube by an image
00157 *    12) adds a spectrum (in z-direction) to all data 
00158 *                        points in a cube
00159 *    13) subtracts a spectrum (in z-direction) from all 
00160 *                        data points in a cube
00161 *    14) multiplies a spectrum (in z-direction) to all data 
00162 *                        points in a cube 
00163 *    15) divides all data points of a cube by a spectrum 
00164 *                        (in z-direction) 
00165 *    16) averaging routine to get a better spectral S/N, sorts
00166 *        the values of the same z-position, cuts the lowest and 
00167 *        highest values according to given thresholds and then
00168 *        takes the average within the x-y plane , cannot have 
00169 *        a sum of low and high rejected values greater than 90%
00170 *        of all values
00171 *    17) determines the sinfo_new_median value in every pixel position
00172 *        by considering all pixels along the third axis.
00173 *        ZERO pixels in a plane are not considered. If all
00174 *        pixels at a position are not valid the result will
00175 *        be 'ZERO'.
00176 *    18) determines the average value in every pixel position
00177 *        by considering all pixels along the third axis.
00178 *        ZERO pixels in a plane are not considered. If all
00179 *        pixels at a position are not valid the result will
00180 *        be 'ZERO'.
00181 *    19) determines the sum value in every pixel position
00182 *        by considering all pixels along the third axis.
00183 *        ZERO pixels in a plane are not considered. If all
00184 *        pixels at a position are not valid the result will
00185 *        be 'ZERO'.
00186 *    20) determines the average value in every pixel position
00187 *        by considering only the pixels along the third axis
00188 *        which lie between the given wavelength values.
00189 *        These values are first recalculated to plane indices
00190 *        by using the given dispersion and minimum wavelength in 
00191 *        the cube.
00192 *        ZERO pixels in a plane are not considered. If all
00193 *        pixels at a position are not valid the result will
00194 *        be 'ZERO'.
00195 *    21) returns the wanted image plane of the cube
00196 *    22) returns the wanted single spectrum of the cube
00197 *    23) merges jittered data cubes to one bigger cube
00198 *        by averaging the overlap regions weighted by
00199 *        the integration times. The x, y size of the final data
00200 *        cube is user given, and should be between 32 and 64
00201 *        pixels, while the relative pixel-offset (sub-pixel
00202 *        accuracy) of the single cubes with respect to the
00203 *        first cube in the list is read from the SEQ CUMOFFSETX,Y
00204 *        fits header keyword.
00205 *   24)  interpolates bad pixel of an object cube if a bad pixel
00206 *        mask cube is available by using the nearest neighbors
00207 *        in 3 dimensions.
00208 *
00209 *   25)  shifts an imagelist by a given amount to integer pixel accuracy
00210 *   26)  shifts an imagelist by a given amount to sub-pixel accuracy
00211 *   27)  shifts an imagelist by a given amount to sub-pixel accuracy
00212 *   FILES
00213 *
00214 *   ENVIRONMENT
00215 *
00216 *   RETURN VALUES
00217 *
00218 *   CAUTIONS
00219 *
00220 *   EXAMPLES
00221 *
00222 *   SEE ALSO
00223 *
00224 *   BUGS
00225 *
00226 *------------------------------------------------------------------------
00227 */
00228 #ifdef HAVE_CONFIG_H
00229 #  include <config.h>
00230 #endif
00231 
00232 #include "sinfo_vltPort.h"
00233 
00234 /*
00235  * System Headers
00236  */
00237 
00238 /*
00239  * Local Headers
00240  */
00241 #include "sinfo_dfs.h"
00242 #include "sinfo_new_cube_ops.h"
00243 #include <sys/types.h>
00244 #include <sys/times.h>
00245 #include <math.h>
00246 #include "sinfo_resampling.h"
00247 #include "sinfo_function_1d.h"
00248 #include "sinfo_error.h"
00249 #include "sinfo_globals.h"
00250 #include "sinfo_utils_wrappers.h"
00251 /*----------------------------------------------------------------------------
00252  *                          Function codes
00253  *--------------------------------------------------------------------------*/
00277 cpl_imagelist *
00278 sinfo_new_cube_ops(
00279     cpl_imagelist   *   cube1,
00280     cpl_imagelist   *   cube2,
00281     int     operation) 
00282 {
00283 
00284     if (cube1==NULL || cube2==NULL) 
00285     {
00286         sinfo_msg_error("null cubes");
00287         return NULL ;
00288     }
00289 
00290     switch(operation) 
00291     {
00292     case '+':
00293     return sinfo_new_cube_add(cube1, cube2) ;
00294     break ;
00295     case '-':
00296     return sinfo_new_cube_sub(cube1, cube2) ;
00297     break ;
00298 
00299     case '*':
00300     return sinfo_new_cube_mul(cube1, cube2) ;
00301     break ;
00302 
00303     case '/':
00304     return sinfo_new_cube_div(cube1, cube2) ;
00305     break ;
00306 
00307     default:
00308     sinfo_msg_error("illegal requested operation: aborting cube arithmetic") ;
00309     return NULL ;
00310     }
00311 }
00312 
00313  
00314 /*----------------------------------------------------------------------------
00315    Function :   sinfo_new_cube_const_ops()
00316    In       :   1 cube, 1 constant, operation to perform
00317    Out      :   result cube
00318    Job      :   4 operations between a cube and a constant
00319    Notice   :   possible operations are:
00320                 Addition    '+'
00321                 Subtraction     '-'
00322                 Multiplication  '*'
00323                 Division    '/'
00324                 Logarithm   'l'
00325                 Power       '^'
00326                 Exponentiation  'e'
00327 
00328  ---------------------------------------------------------------------------*/
00329 
00330 cpl_imagelist *
00331 sinfo_new_cube_const_ops(
00332     cpl_imagelist   *   c1,
00333     double      constant,
00334     int     operation)
00335 {
00336     int ilx1=0;
00337     int ily1=0;
00338     int inp1=0;
00339     cpl_imagelist* c2=NULL;
00340     cpl_image* img1=NULL;
00341  
00342 
00343 
00344     if (c1 == NULL) 
00345     {
00346          sinfo_msg_error("null cube") ;
00347          return NULL ;
00348     }
00349     inp1=cpl_imagelist_get_size(c1);
00350     img1=cpl_imagelist_get(c1,0);
00351     ilx1=cpl_image_get_size_x(img1);
00352     ily1=cpl_image_get_size_y(img1);
00353 
00354 
00355 
00356 
00357 
00358     if ((constant == 0.0) && (operation == '/')) 
00359     {
00360         sinfo_msg_error("division by zero requested "
00361                         "in cube/constant operation") ;
00362         return NULL ;
00363     }
00364 
00365     if ( NULL == (c2 = cpl_imagelist_new()) ) 
00366     {
00367         sinfo_msg_error ("cannot allocate new cube" ) ;
00368         return NULL ;
00369     }
00370 
00371     c2=cpl_imagelist_duplicate(c1);
00372     if(operation == '+') {
00373       cpl_imagelist_add_scalar(c2,constant);
00374     } else if (operation == '-') {
00375       cpl_imagelist_subtract_scalar(c2,constant);
00376     } else if (operation == '*') {
00377       cpl_imagelist_multiply_scalar(c2,constant);
00378     } else if (operation == '/') {
00379       cpl_imagelist_divide_scalar(c2,constant);
00380 
00381     } else {
00382       sinfo_msg_error("operation not supported");
00383       return NULL;
00384     }
00385     return c2 ;
00386 }
00387 
00388 
00389 /*----------------------------------------------------------------------------
00390  * Function :   sinfo_new_cube_sub()
00391  * In       :   two cubes
00392  * Out      :   result cube 
00393  * Job      :   subtract one cube from another
00394  *--------------------------------------------------------------------------*/
00395 
00396 cpl_imagelist *
00397 sinfo_new_cube_sub(
00398     cpl_imagelist   *   c1,
00399     cpl_imagelist   *   c2
00400 )
00401 {
00402     cpl_imagelist   *                 c3 ;
00403     ulong32         i ;
00404     int                 np ;
00405     int ilx1=0;
00406     int ily1=0;
00407     int inp1=0;
00408     int ilx2=0;
00409     int ily2=0;
00410     int inp2=0;
00411 
00412 
00413     cpl_image* i_img=NULL;
00414     cpl_image* img1=NULL;
00415     cpl_image* img2=NULL;
00416     cpl_image* img3=NULL;
00417     float* p1data=NULL;
00418     float* p2data=NULL;
00419     float* p3data=NULL;
00420 
00421 
00422 
00423     inp1=cpl_imagelist_get_size(c1);
00424     i_img=cpl_imagelist_get(c1,0);
00425     ilx1=cpl_image_get_size_x(i_img);
00426     ily1=cpl_image_get_size_y(i_img);
00427 
00428 
00429     inp2=cpl_imagelist_get_size(c2);
00430     i_img=cpl_imagelist_get(c2,0);
00431     ilx2=cpl_image_get_size_x(i_img);
00432     ily2=cpl_image_get_size_y(i_img);
00433 
00434     if ((ilx1 != ilx2) ||
00435     (ily1 != ily2)) 
00436     {
00437     sinfo_msg_error("incompatible size: cannot subtract") ;
00438     return NULL ;
00439     }
00440 
00441     if ((inp2 != inp1) &&
00442     (inp2 != 1)) 
00443     {
00444     sinfo_msg_error("cannot compute with these number of planes") ;
00445     return NULL ;
00446     }
00447 
00448     if ( NULL == (c3 = cpl_imagelist_new()) ) 
00449     {
00450         sinfo_msg_error ("cannot allocate new cube" ) ;
00451         return NULL ;
00452     }
00453     
00454     for (np=0 ; np < inp1 ; np++) 
00455     {
00456       img3=cpl_image_new(ilx1,ily1,CPL_TYPE_FLOAT);
00457       cpl_imagelist_set(c3,img3,np);
00458     }       
00459 
00460 
00461     for (np=0 ; np < inp1 ; np++) 
00462     {
00463       img1=cpl_imagelist_get(c1,np);
00464       p1data=cpl_image_get_data_float(img1);
00465       img2=cpl_imagelist_get(c2,np);
00466       p2data=cpl_image_get_data_float(img2);
00467       img3=cpl_imagelist_get(c3,np);
00468       p3data=cpl_image_get_data_float(img3);
00469 
00470         for (i=0 ; i< (ulong32)ilx1*ily1 ; i++) 
00471         {
00472             p3data[i] = p1data[i] - p2data[i] ;
00473     } 
00474     }
00475 
00476     return c3 ;
00477 }
00478 
00479 
00480 /*----------------------------------------------------------------------------
00481  * Function :   sinfo_new_cube_add()
00482  * In       :   two cubes
00483  * Out      :   result cube 
00484  * Job      :   add a cube to another
00485  *--------------------------------------------------------------------------*/
00486 
00487 cpl_imagelist *
00488 sinfo_new_cube_add(
00489     cpl_imagelist   *   c1,
00490     cpl_imagelist   *   c2
00491 )
00492 {
00493     cpl_imagelist  *          c3 ;
00494     ulong32     i ;
00495     int         np ;
00496     int ilx1=0;
00497     int ily1=0;
00498     int inp1=0;
00499     int ilx2=0;
00500     int ily2=0;
00501     int inp2=0;
00502 
00503 
00504     cpl_image* i_img=NULL;
00505     cpl_image* img1=NULL;
00506     cpl_image* img2=NULL;
00507     cpl_image* img3=NULL;
00508     float* p1data=NULL;
00509     float* p2data=NULL;
00510     float* p3data=NULL;
00511 
00512 
00513 
00514     inp1=cpl_imagelist_get_size(c1);
00515     i_img=cpl_imagelist_get(c1,0);
00516     ilx1=cpl_image_get_size_x(i_img);
00517     ily1=cpl_image_get_size_y(i_img);
00518 
00519 
00520     inp2=cpl_imagelist_get_size(c2);
00521     i_img=cpl_imagelist_get(c2,0);
00522     ilx2=cpl_image_get_size_x(i_img);
00523     ily2=cpl_image_get_size_y(i_img);
00524     if ((ilx1 != ilx2) || (ily1 != ily2)) 
00525     {
00526     sinfo_msg_error("incompatible size: cannot add") ;
00527     return NULL ;
00528     }
00529     if ((inp2 != inp1) && (inp2 != 1)) 
00530     {
00531     sinfo_msg_error("cannot compute with these number of planes") ;
00532     return NULL ;
00533     }
00534 
00535     if (NULL == (c3 = cpl_imagelist_new()) )
00536     {
00537         sinfo_msg_error ("cannot allocate new cube") ; 
00538         return NULL ;
00539     }
00540 
00541     for (np=0 ; np < inp1 ; np++) 
00542     {
00543       img3=cpl_image_new(ilx1,ily1,CPL_TYPE_FLOAT);
00544       cpl_imagelist_set(c3,img3,np);
00545     }       
00546 
00547     for (np=0 ; np < inp1 ; np++) 
00548     {
00549       img1=cpl_imagelist_get(c1,np);
00550       p1data=cpl_image_get_data_float(img1);
00551       img2=cpl_imagelist_get(c2,np);
00552       p2data=cpl_image_get_data_float(img2);
00553       img3=cpl_imagelist_get(c3,np);
00554       p3data=cpl_image_get_data_float(img3);
00555         for (i=0 ; i< (ulong32)ilx1*ily1 ; i++) 
00556         {
00557         p3data[i] = p1data[i] + p2data[i] ;
00558         }
00559     }
00560 
00561     return c3 ;
00562 }
00563 
00564 /*----------------------------------------------------------------------------
00565  * Function :   sinfo_new_cube_mul()
00566  * In       :   two cubes
00567  * Out      :   result cube
00568  * Job      :   multiply 2 cubes    
00569  *--------------------------------------------------------------------------*/
00570 
00571 cpl_imagelist *
00572 sinfo_new_cube_mul(
00573     cpl_imagelist   *   c1,
00574     cpl_imagelist   *   c2
00575 )
00576 {
00577     cpl_imagelist             *c3 ;
00578     ulong32     i ;
00579     int             np ;
00580     int ilx1=0;
00581     int ily1=0;
00582     int inp1=0;
00583     int ilx2=0;
00584     int ily2=0;
00585     int inp2=0;
00586 
00587 
00588     cpl_image* i_img=NULL;
00589     cpl_image* img1=NULL;
00590     cpl_image* img2=NULL;
00591     cpl_image* img3=NULL;
00592     float* p1data=NULL;
00593     float* p2data=NULL;
00594     float* p3data=NULL;
00595 
00596 
00597 
00598 
00599     inp1=cpl_imagelist_get_size(c1);
00600     i_img=cpl_imagelist_get(c1,0);
00601     ilx1=cpl_image_get_size_x(i_img);
00602     ily1=cpl_image_get_size_y(i_img);
00603 
00604 
00605     inp2=cpl_imagelist_get_size(c2);
00606     i_img=cpl_imagelist_get(c2,0);
00607     ilx2=cpl_image_get_size_x(i_img);
00608     ily2=cpl_image_get_size_y(i_img);
00609 
00610     if ((ilx1 != ilx2) || (ily1 != ily2)) 
00611     {
00612     sinfo_msg_error("incompatible size: cannot multiply") ;
00613     return NULL ;
00614     }
00615 
00616     if ((inp2 != inp1) && (inp2 != 1)) 
00617     {
00618     sinfo_msg_error("cannot compute with these number of planes") ;
00619     return NULL ;
00620     }
00621 
00622     if ( NULL == (c3 = cpl_imagelist_new()) )
00623     {
00624         sinfo_msg_error ("cannot allocate new cube" ) ;
00625         return NULL ;
00626     }
00627 
00628 
00629     for (np=0 ; np < inp1 ; np++) 
00630     {
00631       img3=cpl_image_new(ilx1,ily1,CPL_TYPE_FLOAT);
00632       cpl_imagelist_set(c3,img3,np);
00633     }       
00634 
00635     for (np=0 ; np < inp1 ; np++) 
00636     {
00637       img1=cpl_imagelist_get(c1,np);
00638       p1data=cpl_image_get_data_float(img1);
00639       img2=cpl_imagelist_get(c2,np);
00640       p2data=cpl_image_get_data_float(img2);
00641       img3=cpl_imagelist_get(c3,np);
00642       p3data=cpl_image_get_data_float(img3);
00643         for (i=0 ; i< (ulong32)ilx1*ilx2 ; i++) 
00644         {
00645             p3data[i] = p1data[i] * p2data[i] ;
00646         }
00647     }
00648 
00649     return c3 ;
00650 }
00651 
00652 
00653 /*----------------------------------------------------------------------------
00654  * Function :   sinfo_new_cube_div()
00655  * In       :   two cubes
00656  * Out      :   result cube 
00657  * Job      :   divide 2 cubes  
00658  *--------------------------------------------------------------------------*/
00659 
00660 cpl_imagelist *
00661 sinfo_new_cube_div(
00662     cpl_imagelist   *   c1,
00663     cpl_imagelist   *   c2
00664 )
00665 {
00666     cpl_imagelist *           c3 ;
00667     ulong32     i ;
00668     int         np ;
00669     int ilx1=0;
00670     int ily1=0;
00671     int inp1=0;
00672     int ilx2=0;
00673     int ily2=0;
00674     int inp2=0;
00675 
00676 
00677     cpl_image* i_img=NULL;
00678     cpl_image* img1=NULL;
00679     cpl_image* img2=NULL;
00680     cpl_image* img3=NULL;
00681     float* p1data=NULL;
00682     float* p2data=NULL;
00683     float* p3data=NULL;
00684 
00685 
00686     inp1=cpl_imagelist_get_size(c1);
00687     i_img=cpl_imagelist_get(c1,0);
00688     ilx1=cpl_image_get_size_x(i_img);
00689     ily1=cpl_image_get_size_y(i_img);
00690 
00691 
00692     inp2=cpl_imagelist_get_size(c2);
00693     i_img=cpl_imagelist_get(c2,0);
00694     ilx2=cpl_image_get_size_x(i_img);
00695     ily2=cpl_image_get_size_y(i_img);
00696 
00697     if ((ilx1 != ilx2) ||
00698     (ily1 != ily2)) 
00699     {
00700     sinfo_msg_error("incompatible size: cannot divide") ;
00701     return NULL ;
00702     }
00703 
00704     if ((inp2 != inp1) && (inp2 != 1)) 
00705     {
00706     sinfo_msg_error("cannot compute with these number of planes") ;
00707     return NULL ;
00708     }
00709 
00710     if (NULL == (c3 = cpl_imagelist_new()) ) 
00711     {
00712         sinfo_msg_error ("cannot allocate a new cube") ;
00713         return NULL ;
00714     }
00715 
00716     for (np=0 ; np < inp1 ; np++) 
00717     {
00718       img3=cpl_image_new(ilx1,ily1,CPL_TYPE_FLOAT);
00719       cpl_imagelist_set(c3,img3,np);
00720     }       
00721         
00722     for (np=0 ; np < inp1 ; np++) 
00723     {
00724       img1=cpl_imagelist_get(c1,np);
00725       p1data=cpl_image_get_data_float(img1);
00726       img2=cpl_imagelist_get(c2,np);
00727       p2data=cpl_image_get_data_float(img2);
00728       img3=cpl_imagelist_get(c3,np);
00729       p3data=cpl_image_get_data_float(img3);
00730 
00731 
00732         for (i=0 ; i< (ulong32) ilx1*ily1 ; i++) 
00733         {
00734             if (fabs((double)p2data[i]) < 1e-10)
00735             {
00736             p3data[i] = 0.0 ;
00737             }
00738             else
00739             {
00740                 p3data[i] = p1data[i] / p2data[i] ;
00741             }
00742         }
00743     }
00744 
00745     return c3 ;
00746 }
00747 
00748 
00749 
00750 /*---------------------------------------------------------------------------
00751    Function :   sinfo_new_add_image_to_cube()
00752    In       :   1 allocated cube, 1 allocated image
00753    Out      :   result cube
00754    Job      :   add an image to all planes in the cube
00755  ---------------------------------------------------------------------------*/
00756 
00757 cpl_imagelist * 
00758 sinfo_new_add_image_to_cube(cpl_imagelist * cu, cpl_image * im)
00759 {
00760     cpl_imagelist *   cube ;
00761     int            i ;
00762     int clx=0;
00763     int cly=0;
00764     int cnp=0;
00765     int ilx=0;
00766     int ily=0;
00767 
00768 
00769     cpl_image* i_img=NULL;
00770 
00771     if (cu==NULL || im==NULL) 
00772     {
00773        sinfo_msg_error ("null cube or null image") ;
00774        return NULL ;
00775     }
00776     cnp=cpl_imagelist_get_size(cu);
00777     i_img=cpl_imagelist_get(cu,0);
00778     clx=cpl_image_get_size_x(i_img);
00779     cly=cpl_image_get_size_y(i_img);
00780 
00781     ilx=cpl_image_get_size_x(im);
00782     ily=cpl_image_get_size_y(im);
00783 
00784     if ((clx != ilx) || (cly != ily)) 
00785     {
00786         sinfo_msg_error("incompatible size: cannot add image to cube") ;
00787     return NULL ;
00788     }
00789 
00790     cube = cpl_imagelist_duplicate (cu) ;
00791 
00792     for (i=0 ; i<cnp ; i++)
00793     {
00794       /* AMO 
00795         here may be we have to use cpl_image_add_create and cpl_imagelist_set 
00796        */
00797     cpl_image_add(cpl_imagelist_get(cube,i), im) ;
00798     }
00799 
00800     return cube ;
00801 }
00802 
00803 /*---------------------------------------------------------------------------
00804    Function :   sinfo_new_sub_image_from_cube()
00805    In       :   1 allocated cube, 1 allocated image
00806    Out      :       result cube
00807    Job      :   subtract an image from all planes in the cube
00808  ---------------------------------------------------------------------------*/
00809 
00810 cpl_imagelist * 
00811 sinfo_new_sub_image_from_cube (cpl_imagelist * cu, cpl_image * im)
00812 {
00813     cpl_imagelist   * cube ;
00814     int            i ;
00815     int clx=0;
00816     int cly=0;
00817     int cnp=0;
00818     int ilx=0;
00819     int ily=0;
00820 
00821 
00822     cpl_image* i_img=NULL;
00823 
00824     if (cu==NULL || im==NULL) 
00825     {
00826         sinfo_msg_error ("null cube or null image") ;
00827         return NULL ;
00828     }
00829     cnp=cpl_imagelist_get_size(cu);
00830     i_img=cpl_imagelist_get(cu,0);
00831     clx=cpl_image_get_size_x(i_img);
00832     cly=cpl_image_get_size_y(i_img);
00833 
00834     ilx=cpl_image_get_size_x(im);
00835     ily=cpl_image_get_size_y(im);
00836 
00837     if ((clx != ilx) || (cly != ily)) 
00838     {
00839 
00840     sinfo_msg_error("incompatible size: cannot subtract image from cube") ;
00841         return NULL ;
00842     }
00843 
00844     cube = cpl_imagelist_duplicate (cu) ;
00845 
00846     for (i=0 ; i<cnp ; i++) 
00847     {
00848       /* AMO 
00849         here may be we have to use cpl_image_add_create and cpl_imagelist_set 
00850        */
00851     cpl_image_subtract(cpl_imagelist_get(cube,i), im) ;
00852     }
00853     return cube ;
00854 }
00855 
00856 /*---------------------------------------------------------------------------
00857    Function :   sinfo_new_mul_image_to_cube()
00858    In       :   1 allocated cube, 1 allocated image
00859    Out      :   result cube
00860    Job      :   multiply an image to all planes in the cube
00861  ---------------------------------------------------------------------------*/
00862 
00863 cpl_imagelist * 
00864 sinfo_new_mul_image_to_cube(cpl_imagelist * cu, cpl_image * im)
00865 {
00866     cpl_imagelist   * cube ;  
00867     int            i ;
00868     int clx=0;
00869     int cly=0;
00870     int cnp=0;
00871     int ilx=0;
00872     int ily=0;
00873 
00874 
00875     cpl_image* i_img=NULL;
00876  
00877     if (cu==NULL || im==NULL) 
00878     {
00879         sinfo_msg_error("null cube or null image") ;
00880         return NULL ;
00881     }
00882     cnp=cpl_imagelist_get_size(cu);
00883     i_img=cpl_imagelist_get(cu,0);
00884     clx=cpl_image_get_size_x(i_img);
00885     cly=cpl_image_get_size_y(i_img);
00886 
00887     ilx=cpl_image_get_size_x(im);
00888     ily=cpl_image_get_size_y(im);
00889 
00890     if ((clx != ilx) || (cly != ily)) 
00891     {
00892     sinfo_msg_error("incompatible size: cannot multiply image by cube") ;
00893     return NULL ;
00894     }
00895 
00896     cube = cpl_imagelist_duplicate (cu) ;
00897 
00898     for (i=0 ; i<cnp ; i++) 
00899     { 
00900       /* AMO 
00901         here may be we have to use cpl_image_add_create and cpl_imagelist_set 
00902        */
00903     cpl_image_multiply(cpl_imagelist_get(cube,i), im) ;
00904     }
00905 
00906     return cube ;
00907 }
00908 
00909 /*---------------------------------------------------------------------------
00910    Function :   sinfo_new_div_cube_by_image()
00911    In       :   1 allocated cube, 1 allocated image
00912    Out      :   result cube
00913    Job      :   divide all planes in the cube by an image
00914  ---------------------------------------------------------------------------*/
00915 
00916 cpl_imagelist * 
00917 sinfo_new_div_cube_by_image(cpl_imagelist * cu, cpl_image * im)
00918 {
00919     cpl_imagelist   * cube ;
00920     int            i ;
00921     int clx=0;
00922     int cly=0;
00923     int cnp=0;
00924     int ilx=0;
00925     int ily=0;
00926 
00927 
00928     cpl_image* i_img=NULL;
00929 
00930     if (cu==NULL || im==NULL) 
00931     {
00932         sinfo_msg_error ("null cube or null image") ;
00933         return NULL ;
00934     }
00935     cnp=cpl_imagelist_get_size(cu);
00936     i_img=cpl_imagelist_get(cu,0);
00937     clx=cpl_image_get_size_x(i_img);
00938     cly=cpl_image_get_size_y(i_img);
00939 
00940     ilx=cpl_image_get_size_x(im);
00941     ily=cpl_image_get_size_y(im);
00942 
00943     if ((clx != ilx) || (cly != ily)) 
00944     {
00945     sinfo_msg_error("incompatible size: cannot divide cube by image") ;
00946     return NULL ;
00947     }
00948 
00949     cube = cpl_imagelist_duplicate (cu) ;
00950 
00951     for (i=0 ; i<cnp ; i++) 
00952     {
00953       /* AMO 
00954         here may be we have to use cpl_image_add_create and cpl_imagelist_set 
00955        */
00956     cpl_image_divide(cpl_imagelist_get(cube,i), im) ;
00957     }
00958 
00959     return cube ;
00960 }
00961 
00962 
00963 /*---------------------------------------------------------------------------
00964    Function :   sinfo_new_add_spectrum_to_cube()
00965    In       :   1 allocated cube, 1 allocated spectrum sinfo_vector
00966    Out      :   result cube
00967    Job      :   adds a spectrum (in z-direction) to all data 
00968                         points in a cube
00969  ---------------------------------------------------------------------------*/
00970 
00971 cpl_imagelist * 
00972 sinfo_new_add_spectrum_to_cube(cpl_imagelist *cu, Vector *spec)
00973 {
00974     cpl_imagelist *   cube ;
00975     int         i ,j ;
00976     int ilx=0;
00977     int ily=0;
00978     int inp=0;
00979     float* pidata=NULL;
00980     float* podata=NULL;
00981     cpl_image* i_img=NULL;
00982     cpl_image* o_img=NULL;
00983     
00984     if (cu == NULL || spec == NULL)
00985     {
00986         sinfo_msg_error ("null cube or null spectrum") ;
00987         return NULL ;
00988     }
00989     inp=cpl_imagelist_get_size(cu);
00990     i_img=cpl_imagelist_get(cu,0);
00991     ilx=cpl_image_get_size_x(i_img);
00992     ily=cpl_image_get_size_y(i_img);
00993 
00994     if ( inp != spec -> n_elements )
00995     {
00996         sinfo_msg_error("cube length and spectrum length are not compatible") ;
00997         return NULL ;
00998     }
00999  
01000     if ( NULL == (cube = cpl_imagelist_new ()) ) 
01001     {
01002         sinfo_msg_error ("cannot allocate new cube" ) ;
01003         return NULL ;
01004     }    
01005     for ( i = 0; i < inp; i++)
01006     {
01007       o_img=cpl_image_new(ilx,ily,CPL_TYPE_FLOAT);
01008       cpl_imagelist_set(cube,o_img,i);
01009     } 
01010 
01011 
01012     for ( i = 0; i < inp; i++)
01013     {
01014       i_img=cpl_imagelist_get(cu,i);
01015       pidata=cpl_image_get_data_float(i_img);
01016       o_img=cpl_imagelist_get(cube,i);
01017       podata=cpl_image_get_data_float(o_img);
01018         for ( j = 0; j < (int) ilx*ily; j++)
01019         {       
01020             podata[j] = pidata[j] + spec -> data[i] ;
01021         }
01022     }    
01023     
01024     return cube ;
01025 }
01026 
01027 /*---------------------------------------------------------------------------
01028    Function :   sinfo_new_sub_spectrum_from_cube()
01029    In       :   1 allocated cube, 1 allocated spectrum sinfo_vector
01030    Out      :   result cube
01031    Job      :   subtracts a spectrum (in z-direction) from all 
01032                         data points in a cube
01033  ---------------------------------------------------------------------------*/
01034 
01035 cpl_imagelist * 
01036 sinfo_new_sub_spectrum_from_cube(cpl_imagelist *cu, Vector *spec)
01037 {
01038     cpl_imagelist *   cube ;
01039     int         i ,j ;
01040     int ilx=0;
01041     int ily=0;
01042     int inp=0;
01043     float* pidata=NULL;
01044     float* podata=NULL;
01045     cpl_image* i_img=NULL;
01046     cpl_image* o_img=NULL;
01047     
01048     if (cu == NULL || spec == NULL)
01049     {
01050         sinfo_msg_error ("null cube or null spectrum") ;
01051         return NULL ;
01052     }
01053     inp=cpl_imagelist_get_size(cu);
01054     i_img=cpl_imagelist_get(cu,0);
01055     ilx=cpl_image_get_size_x(i_img);
01056     ily=cpl_image_get_size_y(i_img);
01057  
01058     if ( inp != spec -> n_elements )
01059     {
01060         sinfo_msg_error("cube length and spectrum length are not compatible") ;
01061         return NULL ;
01062     }
01063 
01064     if ( NULL == (cube = cpl_imagelist_new()) ) 
01065     {
01066         sinfo_msg_error ("cannot allocate new cube" ) ;
01067         return NULL ;
01068     }
01069     for ( i = 0; i < inp; i++)
01070     {
01071       o_img=cpl_image_new(ilx,ily,CPL_TYPE_FLOAT);
01072       cpl_imagelist_set(cube,o_img,i);
01073     } 
01074 
01075     for ( i = 0; i < inp; i++)
01076     {
01077       i_img=cpl_imagelist_get(cu,i);
01078       pidata=cpl_image_get_data_float(i_img);
01079       o_img=cpl_imagelist_get(cube,i);
01080       podata=cpl_image_get_data_float(o_img);
01081         for ( j = 0; j < (int) ilx*ily; j++)
01082         {       
01083             if ( isnan(pidata[j]) || isnan(spec -> data[i]) )
01084             {
01085                 podata[j] = ZERO ;
01086             }
01087             else
01088             {
01089                 podata[j] = pidata[j] - spec -> data[i] ;
01090             }
01091         }
01092     }    
01093     
01094     return cube ;
01095 }
01096 
01097 
01098 /*---------------------------------------------------------------------------
01099    Function :   sinfo_new_mul_spectrum_to_cube()
01100    In       :   1 allocated cube, 1 allocated spectrum sinfo_vector
01101    Out      :   result cube
01102    Job      :   multiplies a spectrum (in z-direction) to all data 
01103                         points in a cube 
01104  ---------------------------------------------------------------------------*/
01105 
01106 cpl_imagelist * 
01107 sinfo_new_mul_spectrum_to_cube(cpl_imagelist *cu, Vector *spec)
01108 {
01109     cpl_imagelist *   cube ;
01110     int         i ,j ;
01111     int ilx=0;
01112     int ily=0;
01113     int inp=0;
01114     float* pidata=NULL;
01115     float* podata=NULL;
01116     cpl_image* i_img=NULL;
01117     cpl_image* o_img=NULL;
01118     
01119     if (cu == NULL || spec == NULL)
01120     {
01121         sinfo_msg_error ("null cube or null spectrum") ;
01122         return NULL ;
01123     }
01124     inp=cpl_imagelist_get_size(cu);
01125     i_img=cpl_imagelist_get(cu,0);
01126     ilx=cpl_image_get_size_x(i_img);
01127     ily=cpl_image_get_size_y(i_img);
01128 
01129     if ( inp != spec -> n_elements )
01130     {
01131         sinfo_msg_error("cube length and spectrum length are not compatible") ;
01132         return NULL ;
01133     }
01134 
01135     if ( NULL == (cube = cpl_imagelist_new ()) ) 
01136     {
01137         sinfo_msg_error ("cannot allocate new cube" ) ;
01138         return NULL ;
01139     }
01140 
01141     for ( i = 0; i < inp; i++)
01142     {
01143       o_img=cpl_image_new(ilx,ily,CPL_TYPE_FLOAT);
01144       cpl_imagelist_set(cube,o_img,i);
01145     } 
01146 
01147     for ( i = 0; i < inp; i++)
01148     {
01149       i_img=cpl_imagelist_get(cu,i);
01150       pidata=cpl_image_get_data_float(i_img);
01151       o_img=cpl_imagelist_get(cube,i);
01152       podata=cpl_image_get_data_float(o_img);
01153         for ( j = 0; j < (int) ilx*ily; j++)
01154         {       
01155             if ( isnan(pidata[j]) || isnan(spec->data[i]) )
01156             {
01157                 podata[j] = ZERO ;
01158             }
01159             else
01160             {
01161                 podata[j] = pidata[j] * spec -> data[i] ;
01162             }
01163         }
01164     }    
01165     
01166     return cube ;
01167 }
01168 
01169 
01170 /*---------------------------------------------------------------------------
01171    Function :   sinfo_new_div_cube_by_spectrum()
01172    In       :   1 allocated cube, 1 allocated spectrum sinfo_vector
01173    Out      :   result cube
01174    Job      :   divides all data points of a cube by a spectrum 
01175                         (in z-direction) 
01176  ---------------------------------------------------------------------------*/
01177 
01178 cpl_imagelist * 
01179 sinfo_new_div_cube_by_spectrum(cpl_imagelist *cu, Vector *spec)
01180 {
01181     cpl_imagelist *   cube ;
01182     float       help ;
01183     int         i ,j ;
01184     int ilx=0;
01185     int ily=0;
01186     int inp=0;
01187     float* pidata=NULL;
01188     float* podata=NULL;
01189     cpl_image* i_img=NULL;
01190     cpl_image* o_img=NULL;
01191     
01192     if (cu == NULL || spec == NULL)
01193     {
01194         sinfo_msg_error ("null cube or null spectrum") ;
01195         return NULL ;
01196     }
01197     inp=cpl_imagelist_get_size(cu);
01198     i_img=cpl_imagelist_get(cu,0);
01199     ilx=cpl_image_get_size_x(i_img);
01200     ily=cpl_image_get_size_y(i_img);
01201  
01202     if ( inp != spec -> n_elements )
01203     {
01204         sinfo_msg_error("cube length and spectrum length are not compatible") ;
01205         return NULL ;
01206     }
01207 
01208     if (NULL == (cube = cpl_imagelist_new ()) ) 
01209     {
01210         sinfo_msg_error ("cannot allocate new cube") ;
01211         return NULL ;
01212     }
01213 
01214     for ( i = 0; i < inp; i++)
01215     {
01216       o_img=cpl_image_new(ilx,ily,CPL_TYPE_FLOAT);
01217       cpl_imagelist_set(cube,o_img,i);
01218     } 
01219 
01220     for ( i = 0; i < inp; i++)
01221     {
01222 
01223       i_img=cpl_imagelist_get(cu,i);
01224       pidata=cpl_image_get_data_float(i_img);
01225       o_img=cpl_imagelist_get(cube,i);
01226       podata=cpl_image_get_data_float(o_img);
01227         for ( j = 0; j < (int) ilx*ily; j++)
01228         {       
01229             if (!isnan(spec->data[i]) && spec->data[i] != 0.)
01230             {
01231                 help = 1/spec->data[i] ;
01232                 if ( help > THRESH )
01233                 {
01234                     help = 1. ;
01235                 }
01236             }
01237             else
01238             {
01239                 help = ZERO ;
01240             }
01241 
01242             if ( isnan(help) || isnan(pidata[j]) )
01243             {
01244                 podata[j] = ZERO ;
01245             }
01246             else
01247             {
01248                 podata[j] = pidata[j] * help ;
01249             }
01250         }
01251     }    
01252     return cube ;
01253 }
01254 
01255 
01256 /*---------------------------------------------------------------------------
01257    Function :   sinfo_new_clean_mean_of_spectra()
01258    In       :   1 allocated cube, position of rectangle in x-y plane ,
01259                         low and high cut threshold
01260    Out      :   result spectrum sinfo_vector
01261    Job      :   averaging routine to get a better spectral S/N, sorts
01262                         the values of the same z-position, cuts the lowest and 
01263                         highest values according to given thresholds and then
01264                         takes the average within the x-y plane , cannot have 
01265                         a sum of low and high rejected values greater than 90%
01266                         of all values
01267  ---------------------------------------------------------------------------*/
01268 
01269 Vector * 
01270 sinfo_new_clean_mean_of_spectra(cpl_imagelist * cube, 
01271                              int llx, 
01272                              int lly, 
01273                              int urx, 
01274                              int ury,
01275                              double lo_reject,
01276                              double hi_reject)
01277 {
01278     Vector                           * mean ;
01279     pixelvalue                   *local_rectangle ;
01280     int                    i, j, k, l, m ;
01281     int             recsize, lo_n, hi_n, nv ;
01282     int ilx=0;
01283     int ily=0;
01284     int inp=0;
01285     float* pidata=NULL;
01286     cpl_image* i_img=NULL;
01287 
01288     if ( cube == NULL || cpl_imagelist_get_size(cube) < 1 )
01289     {
01290         sinfo_msg_error ("no cube to take the mean of his spectra") ;
01291         return NullVector ;
01292     }
01293     inp=cpl_imagelist_get_size(cube);
01294     i_img=cpl_imagelist_get(cube,0);
01295     ilx=cpl_image_get_size_x(i_img);
01296     ily=cpl_image_get_size_y(i_img);
01297  
01298     if ((llx<1) || (llx>ilx) ||
01299         (urx<1) || (urx>ilx) ||
01300         (lly<1) || (lly>ily) ||
01301         (ury<1) || (ury>ily) ||
01302         (llx>=urx) || (lly>=ury)) 
01303     {
01304         sinfo_msg_error("invalid rectangle coordinates:") ;
01305         sinfo_msg_error("lower left is [%d %d] upper right is [%d %d]", 
01306                         llx, lly, urx, ury) ;
01307         return NullVector ;
01308     }
01309 
01310     if ((lo_reject + hi_reject) > 0.9) 
01311     {
01312         sinfo_msg_error("illegal rejection thresholds: [%f] and [%f]", 
01313                         lo_reject, hi_reject) ;
01314         sinfo_msg_error("threshold sum should not be over 0.9"
01315                         " aborting average") ;
01316         return NullVector ;
01317     }
01318 
01319     /* shift from FITS coordinates to C coordinates */
01320     llx -- ; 
01321     lly -- ;
01322     urx -- ; 
01323     ury -- ;
01324     
01325     recsize = (urx - llx + 1) * (ury - lly + 1) ;
01326 
01327     lo_n = (int) (recsize * lo_reject + 0.5) ;
01328     hi_n = (int) (recsize * hi_reject + 0.5) ;
01329 
01330     if (lo_n + hi_n >= recsize)
01331     { 
01332         sinfo_msg_error ("everything would be rejected") ;
01333         return NullVector;
01334     }
01335  
01336     /* allocate a new sinfo_vector to store the average spectral values */
01337     if (NULL == (mean = sinfo_new_vector (inp)) ) 
01338     {
01339         sinfo_msg_error ("cannot allocate a new sinfo_vector") ;
01340         return NullVector ;
01341     }
01342 
01343     /*------------------------------------------------------------------------  
01344      *  loop through the cube planes, through the x axis and the y-axis of the 
01345      *  plane rectangle and store pixel values in a buffer.
01346      */
01347     for ( i = 0 ; i < inp ; i++ )
01348     {
01349       i_img=cpl_imagelist_get(cube,i);
01350       pidata=cpl_image_get_data_float(i_img);
01351       m = 0 ;
01352       local_rectangle=(pixelvalue *)cpl_calloc(recsize, sizeof (pixelvalue*));
01353 
01354         for ( j = lly ; j <= ury ; j++ )
01355         {
01356             for ( k = llx ; k <= urx ; k++ )
01357             {
01358                 local_rectangle[m] = pidata[k + j * ilx] ;
01359                 m ++ ;
01360             }             
01361         }
01362         /*sorts the pixelvalues in the buffer*/
01363         sinfo_pixel_qsort (local_rectangle, recsize) ; 
01364       
01365         nv = 0 ;
01366         for ( l = lo_n ; l < (recsize - hi_n) ; l++ )
01367         {
01368             mean -> data[i] += local_rectangle[l] ;
01369             nv ++;
01370         }
01371         mean -> data[i] /= nv ;
01372 
01373         cpl_free ( local_rectangle ) ;
01374     }
01375     return mean ;
01376 }
01377 
01378 
01379 /*---------------------------------------------------------------------------
01380    Function :sinfo_new_median_cube()
01381    In       :1 allocated cube
01382    Out      :result image
01383    Job      :determines the sinfo_new_median value in every pixel position
01384                  by considering all pixels along the third axis.
01385                  ZERO pixels in a plane are not considered. If all
01386                  pixels at a position are not valid the result will
01387                  be 'ZERO'.
01388  ---------------------------------------------------------------------------*/
01389 cpl_image * 
01390 sinfo_new_median_cube(cpl_imagelist * cube) 
01391 {
01392     cpl_image  *         im ;
01393     pixelvalue *    buffer ;
01394     int        i, j, k, nz ;
01395     int ilx=0;
01396     int ily=0;
01397     int inp=0;
01398     float* pidata=NULL;
01399     float* podata=NULL;
01400     cpl_image* i_img=NULL;
01401 
01402     if ( cube == NULL )
01403     {
01404         sinfo_msg_error ("null cube") ;
01405         return NULL ;
01406     }
01407     inp=cpl_imagelist_get_size(cube);
01408     i_img=cpl_imagelist_get(cube,0);
01409     ilx=cpl_image_get_size_x(i_img);
01410     ily=cpl_image_get_size_y(i_img);
01411 
01412     /* allocate memory */
01413     if (NULL == (im = cpl_image_new (ilx, ily, CPL_TYPE_FLOAT )) ) 
01414     {
01415         sinfo_msg_error ("cannot allocate new image") ;
01416         return NULL ;
01417     }
01418 
01419     /*------------------------------------------------------------------------ 
01420      * transfer each sinfo_vector in z direction in a buffer and collect 
01421        only non-blank data.
01422      */
01423      
01424     podata=cpl_image_get_data_float(im);
01425     for ( i = 0 ; i < (int) ilx*ily ; i++ )
01426     {
01427         buffer = (pixelvalue *) cpl_calloc (inp, sizeof (pixelvalue *));
01428         k = 0 ;
01429         for ( j = 0 ; j < inp ; j ++ )  
01430         {
01431           i_img=cpl_imagelist_get(cube,j);
01432       pidata=cpl_image_get_data_float(i_img);
01433             if ( !isnan(pidata[i]) )
01434             {
01435                 buffer[k] = pidata[i] ;
01436                 k ++ ;
01437             }
01438         } 
01439         nz = k ;
01440 
01441         /* proceed depending on the number of valid pixels */
01442         if ( nz > 2 )
01443         {
01444             podata[i] = sinfo_new_median ( buffer, nz ) ;
01445         }
01446         else if (nz == 2)
01447         {
01448             podata[i] = (buffer[0] + buffer[1]) / 2. ;
01449         }
01450         else if (nz == 1)
01451         {
01452             podata[i] = buffer[0] ;
01453         }
01454         else if (nz == 0)
01455         {       
01456             podata[i] = ZERO ;
01457         }
01458 
01459         cpl_free ( buffer ) ;
01460     }
01461 
01462     return im ;
01463 }
01464 
01465 
01466 /*---------------------------------------------------------------------------
01467    Function :   sinfo_new_average_cube_to_image()
01468    In       :   1 allocated cube
01469    Out      :   result image
01470    Job      :   determines the average value in every pixel position
01471                         by considering all pixels along the third axis.
01472                         ZERO pixels in a plane are not considered. If all
01473                         pixels at a position are not valid the result will
01474                         be 'ZERO'.
01475  ---------------------------------------------------------------------------*/
01476 cpl_image * 
01477 sinfo_new_average_cube_to_image(cpl_imagelist * cube) 
01478 {
01479     cpl_image  *      im ;
01480     int        i, j, nz ;
01481     int ilx=0;
01482     int ily=0;
01483     int inp=0;
01484     float* pidata=NULL;
01485     float* podata=NULL;
01486     cpl_image* i_img=NULL;
01487 
01488     if ( cube == NULL )
01489     {
01490         sinfo_msg_error ("null cube") ;
01491         return NULL ;
01492     }
01493     inp=cpl_imagelist_get_size(cube);
01494     i_img=cpl_imagelist_get(cube,0);
01495     ilx=cpl_image_get_size_x(i_img);
01496     ily=cpl_image_get_size_y(i_img);
01497 
01498     /* allocate memory */
01499     if (NULL == (im = cpl_image_new (ilx, ily,CPL_TYPE_FLOAT )) ) 
01500     {
01501         sinfo_msg_error ("cannot allocate new image") ;
01502         return NULL ;
01503     }
01504 
01505     /*------------------------------------------------------------------------ 
01506      * transfer each vector in z direction in a buffer and collect 
01507        only non-blank data.
01508      */
01509     
01510     podata=cpl_image_get_data_float(im);
01511     for ( i = 0 ; i < (int) ilx*ily ; i++ )
01512     {
01513         nz = 0 ;
01514         for ( j = 0 ; j < inp ; j ++ )  
01515         {
01516           i_img=cpl_imagelist_get(cube,j);
01517       pidata=cpl_image_get_data_float(i_img);
01518             if ( !isnan(pidata[i]) )
01519             {
01520                 nz ++ ;
01521                 podata[i] += pidata[i] ; 
01522             }
01523         } 
01524         
01525         /* proceed depending on the number of valid pixels */
01526         if ( nz >= 1 )
01527         {
01528             podata[i] /= nz ;
01529         }
01530         else if (nz == 0)
01531         {       
01532             podata[i] = ZERO ;
01533         }
01534     }
01535 
01536     return im ;
01537 }
01538 
01539 /*---------------------------------------------------------------------------
01540    Function     :       sinfo_new_sum_cube_to_image()
01541    In           :       1 allocated cube
01542    Out          :       result image
01543    Job          :       determines the sum value in every pixel position
01544                         by considering all pixels along the third axis.
01545                         ZERO pixels in a plane are not considered. If all
01546                         pixels at a position are not valid the result will
01547                         be 'ZERO'.
01548  ---------------------------------------------------------------------------*/
01549 cpl_image * 
01550 sinfo_new_sum_cube_to_image(cpl_imagelist * cube)
01551 {
01552     cpl_image  *      im ;
01553     int        i, j, nz ;
01554     int ilx=0;
01555     int ily=0;
01556     int inp=0;
01557     float* pidata=NULL;
01558     float* podata=NULL;
01559     cpl_image* i_img=NULL;
01560 
01561     if ( cube == NULL )
01562     {
01563         sinfo_msg_error ("null cube") ;
01564         return NULL ;
01565     }
01566     inp=cpl_imagelist_get_size(cube);
01567     i_img=cpl_imagelist_get(cube,0);
01568     ilx=cpl_image_get_size_x(i_img);
01569     ily=cpl_image_get_size_y(i_img);
01570 
01571     /* allocate memory */
01572     if (NULL == (im = cpl_image_new (ilx, ily, CPL_TYPE_FLOAT )) )
01573     {
01574         sinfo_msg_error ("cannot allocate new image") ;
01575         return NULL ;
01576     }
01577 
01578     /*-------------------------------------------------------------------------
01579      * transfer each vector in z direction in a buffer and collect only 
01580        non-blank data.
01581      */
01582 
01583     podata=cpl_image_get_data_float(im);
01584     for ( i = 0 ; i < (int) ilx*ily ; i++ )
01585     {
01586         nz = 0 ;
01587         for ( j = 0 ; j < inp ; j ++ )
01588         {
01589           i_img=cpl_imagelist_get(cube,j);
01590       pidata=cpl_image_get_data_float(i_img);
01591             if ( !isnan(pidata[i]) )
01592             {
01593                 nz++ ;
01594                 podata[i] += pidata[i] ;
01595             }
01596         }
01597 
01598         /* proceed depending on the number of valid pixels */
01599         if (nz == 0)
01600         {
01601             podata[i] = ZERO ;
01602         }
01603     }
01604 
01605     return im ;
01606 }
01607 
01608 /*---------------------------------------------------------------------------
01609    Function sinfo_new_average_cube_to_image_between_waves()
01610    In       cube: data cube to collapse
01611                 dispersion: dispersion per pixel in microns/pixel 
01612                 (derived from fits header information)
01613                 centralWave: central wavelength in the cube in microns
01614                                        (derived from fits header information)
01615                 initialLambda, finalLambda: wavelength values in microns
01616                                             within which the cube is averaged
01617    Out      :resulting averaged image
01618    Job      :determines the average value in every pixel position
01619                  by considering only the pixels along the third axis
01620                  which lie between the given wavelength values.
01621                  These values are first recalculated to plane indices
01622                  by using the given dispersion and minimum wavelength in 
01623                  the cube.
01624                  ZERO pixels in a plane are not considered. If all
01625                  pixels at a position are not valid the result will
01626                  be 'ZERO'.
01627  ---------------------------------------------------------------------------*/
01628 cpl_image * 
01629 sinfo_new_average_cube_to_image_between_waves (cpl_imagelist * cube, 
01630                                            float     dispersion,
01631                                            float     centralWave,
01632                                            float     initialLambda,
01633                                            float     finalLambda)  
01634 {
01635     cpl_image  *      im ;
01636     int        firstPlane ;
01637     int        lastPlane ;
01638     int        i, j, nz ;
01639     float      minWave ;
01640     int ilx=0;
01641     int ily=0;
01642     int inp=0;
01643     float* pidata=NULL;
01644     float* podata=NULL;
01645     cpl_image* i_img=NULL;
01646 
01647     if ( cube == NULL )
01648     {
01649         sinfo_msg_error ("null cube") ;
01650         return NULL ;
01651     }
01652     i_img=cpl_imagelist_get(cube,0);
01653     ilx=cpl_image_get_size_x(i_img);
01654     ily=cpl_image_get_size_y(i_img);
01655 
01656     inp=cpl_imagelist_get_size(cube);
01657 
01658     minWave = centralWave - (float) (inp / 2)*dispersion ;
01659 
01660     if ( dispersion <= 0. || minWave <= 0. )
01661     {
01662         sinfo_msg_error ("wrong dispersion or minimum wavelength given") ;
01663         return NULL ;
01664     }
01665 
01666     if ( initialLambda < minWave || 
01667         (initialLambda >= minWave + dispersion * inp) )
01668     {
01669         sinfo_msg_error ("wrong initial wavelength given") ;
01670         return NULL ;
01671     }
01672 
01673     if ( finalLambda <= minWave || 
01674         (finalLambda > minWave + dispersion * inp) )
01675     {
01676         sinfo_msg_error ("wrong final wavelength given") ;
01677         return NULL ;
01678     }
01679 
01680     /* allocate memory */
01681     if (NULL == (im = cpl_image_new (ilx, ily, CPL_TYPE_FLOAT )) ) 
01682     {
01683         sinfo_msg_error ("cannot allocate new image") ;
01684         return NULL ;
01685     }
01686 
01687     /* transfer the wavelength range to image plane indices */
01688     firstPlane = sinfo_new_nint ((double) ((initialLambda - minWave) / 
01689                                           dispersion)) ; 
01690     lastPlane  = sinfo_new_nint ((double) ((finalLambda - minWave) / 
01691                                           dispersion)) ; 
01692  
01693     if ( firstPlane < 0 || firstPlane >= inp ||
01694          lastPlane  < 0 || lastPlane  >  inp )
01695     {
01696         sinfo_msg_error ("wrong values given!") ;
01697         return NULL ;
01698     }
01699 
01700     /*------------------------------------------------------------------------ 
01701      * transfer each vector in z direction in a buffer and collect only 
01702        non-blank data.
01703      */
01704  
01705 
01706 
01707     podata=cpl_image_get_data_float(im);   
01708     for ( i = 0 ; i < (int) ilx*ily ; i++ )
01709     {
01710         nz = 0 ;
01711         
01712         for ( j = firstPlane ; j <= lastPlane ; j ++ )  
01713         {
01714           i_img=cpl_imagelist_get(cube,j);
01715       pidata=cpl_image_get_data_float(i_img);
01716             if ( !isnan(pidata[i]) )
01717             {
01718                 nz ++ ;
01719                 podata[i] += pidata[i] ; 
01720             }
01721         } 
01722         
01723         /* proceed depending on the number of valid pixels */
01724         if ( nz >= 1 )
01725         {
01726             podata[i] /= nz ;
01727         }
01728         else if (nz == 0)
01729         {       
01730             podata[i] = ZERO ;
01731         }
01732     }
01733 
01734     return im ;
01735 }
01736 
01737 /*---------------------------------------------------------------------------
01738    Function :   sinfo_new_extract_image_from_cube()
01739    In       :   1 allocated cube
01740                         index of cube plane
01741    Out      :   extracted image
01742    Job      :   returns the wanted image plane of the cube
01743  ---------------------------------------------------------------------------*/
01744 cpl_image * 
01745 sinfo_new_extract_image_from_cube(cpl_imagelist * cube, int plane_index) 
01746 {
01747     if ( cube == NULL )
01748     {
01749         sinfo_msg_error ("null cube") ;
01750         return NULL ;
01751     }
01752 
01753     if ( plane_index < 0 || plane_index >= cpl_imagelist_get_size(cube) )
01754     {
01755         sinfo_msg_error ("wrong plane index for image to be extracted") ;
01756         return NULL ;
01757     }
01758 
01759     return cpl_imagelist_get(cube,plane_index) ;
01760 }
01761 
01762 /*---------------------------------------------------------------------------
01763    Function :sinfo_new_extract_spectrum_from_cube()
01764    In       :cube: 1 allocated cube
01765                  x_pos, y_pos: x, y pixel position of the 
01766                                spectrum counted from 0
01767    Out      :extracted spectral sinfo_vector object
01768    Job      :returns the wanted single spectrum of the cube
01769  ---------------------------------------------------------------------------*/
01770 Vector * 
01771 sinfo_new_extract_spectrum_from_cube(cpl_imagelist * cube, 
01772                                      int x_pos, int y_pos)
01773 {
01774     Vector * returnedSpectrum ;
01775     int i ;
01776     int ilx=0;
01777     int ily=0;
01778     int inp=0;
01779     float* pidata=NULL;
01780     cpl_image* i_img=NULL;
01781 
01782     if ( cube == NULL )
01783     {
01784         sinfo_msg_error ("no cube given!") ;
01785         return NullVector ;
01786     }
01787     i_img=cpl_imagelist_get(cube,0);
01788     ilx=cpl_image_get_size_x(i_img);
01789     ily=cpl_image_get_size_y(i_img);
01790     inp=cpl_imagelist_get_size(cube);
01791  
01792     if ( x_pos < 0 || x_pos >= ilx )
01793     {
01794         sinfo_msg_error ("wrong x-positon of spectrum given!") ;
01795         return NullVector ;
01796     }
01797 
01798     if ( y_pos < 0 || y_pos >= ily )
01799     {
01800         sinfo_msg_error ("wrong y-positon of spectrum given!") ;
01801         return NullVector ;
01802     }
01803 
01804     /* allocate memory */
01805     if ( NULL == (returnedSpectrum = sinfo_new_vector ( inp )) )
01806     {
01807         sinfo_msg_error ("cannot allocate new spectrum!") ;
01808         return NullVector ;
01809     }
01810 
01811     for ( i = 0 ; i < inp ; i++ )
01812     {
01813       i_img=cpl_imagelist_get(cube,i);
01814       pidata=cpl_image_get_data_float(i_img);
01815       returnedSpectrum -> data[i] = pidata[x_pos + ilx*y_pos] ;
01816     }
01817 
01818     return returnedSpectrum ;
01819 }
01820 
01821 /*---------------------------------------------------------------------------
01822    Function     :       sinfo_new_combine_jittered_cubes()
01823    In           :       cubes: list of jittered cubes to mosaic
01824                         mergedCube: resulting merged cube containing the 
01825                                       jittered cubes
01826                         n_cubes: number of cubes in the list to merge
01827                         cumoffsetx,y: array of relative x, y pixel offsets
01828                                       with respect to the first frame in the 
01829                                       same sequence as the cube list.
01830                         exptimes: exposure times array giving the time
01831                                   in the same sequence as the cube list
01832                         kernel_type: the name of the interpolation kernel
01833                                      that you want to generate using the 
01834                                      eclipse routine 
01835                                      sinfo_generate_interpolation_kernel()
01836                                      Supported kernels are:
01837                                      NULL:      default kernel, currently tanh
01838                                      "default": dito
01839                                      "tanh":    Hyperbolic tangent
01840                                      "sinc2":   Square sinc
01841                                      "lanczos": Lanczos2 kernel
01842                                      "hamming": Hamming kernel
01843                                      "hann":    Hann kernel
01844    Out          :       mask: cube of the same size as combinedCube 
01845                               containing 0 for blank (ZERO pixels) and
01846                               the summed integration times for 
01847                               overlapping regions
01848                         mergedCube: final data cube containing the 
01849                                     jittered cubes
01850    Job          :       merges jittered data cubes to one bigger cube
01851                         by averaging the overlap regions weighted by
01852                         the integration times. The x, y size of the final data
01853                         cube is user given, and should be between 32 and 64 
01854                         pixels, while the relative pixel-offset (sub-pixel
01855                         accuracy) of the single cubes with respect to the
01856                         first cube in the list is read from the 
01857                         SEQ CUMOFFSETX,Y
01858                         fits header keyword. 
01859  ---------------------------------------------------------------------------*/
01860 cpl_imagelist * 
01861 sinfo_new_combine_jittered_cubes ( cpl_imagelist ** cubes,
01862                                  cpl_imagelist  * mergedCube,
01863                                  int        n_cubes,
01864                                  float    * cumoffsetx,
01865                                  float    * cumoffsety,
01866                                  float    * exptimes,
01867                                  char     * kernel_type )
01868 {
01869 
01870     int i=0 ;
01871     int x=0;
01872     int y=0;
01873     int z=0; 
01874     int llx0=0;
01875     int lly0=0; 
01876     int posx=0;
01877     int posy=0;
01878     float weight=0;
01879     cpl_imagelist * mask=NULL;
01880     double * kernel=NULL;
01881     /*cpl_image * shiftedImage ;*/
01882 
01883     int* llx=NULL ;
01884     int* lly=NULL ;
01885 
01886     float* sub_offsetx=NULL ;
01887     float* sub_offsety=NULL ;
01888 
01889     cpl_imagelist ** tmpcubes=NULL ; 
01890     pixelvalue * tmpspace=NULL;
01891     
01892 
01893   int ilx=0;
01894   int ily=0;
01895   int olx=0;
01896   int oly=0;
01897   int mlx=0;
01898   int onp=0;
01899   int inp=0;
01900 
01901 
01902 
01903   float* podata=NULL;
01904   float* pmdata=NULL;
01905   float* ptdata=NULL;
01906 
01907   cpl_image* i_img=NULL;
01908   cpl_image* o_img=NULL;
01909   cpl_image* m_img=NULL;
01910   cpl_image* t_img=NULL;
01911 
01912 
01913     if ( cubes == NULL )
01914     {
01915         sinfo_msg_error ("no cube list given!") ;
01916         return NULL ;
01917     }
01918     if ( n_cubes <= 0 )
01919     {
01920         sinfo_msg_error ("wrong number of data cubes in list!") ;
01921         return NULL ;
01922     }
01923     if ( cumoffsetx == NULL || cumoffsety == NULL )
01924     {
01925         sinfo_msg_error ("no cumoffsetx/y given!") ;
01926         return NULL ;
01927     }
01928     if ( exptimes == NULL )
01929     {
01930         sinfo_msg_error ("no exposure time array given!") ;
01931         return NULL ;
01932     }
01933     
01934     o_img=cpl_imagelist_get(mergedCube,0);
01935     olx=cpl_image_get_size_x(o_img);
01936     oly=cpl_image_get_size_y(o_img);
01937     onp=cpl_imagelist_get_size(mergedCube);
01938     if ( NULL == (mask = cpl_imagelist_new()) )
01939     {
01940         sinfo_msg_error ("could not allocate cube!") ;
01941         return NULL ;
01942     }
01943     for(i=0;i<onp;i++){
01944       o_img=cpl_image_new(olx,oly,CPL_TYPE_FLOAT);
01945       cpl_imagelist_set(mergedCube,o_img,i);
01946     }
01947 
01948     i_img=cpl_imagelist_get(cubes[0],0);
01949     ilx=cpl_image_get_size_x(i_img);
01950     ily=cpl_image_get_size_y(i_img);
01951       
01952     inp=cpl_imagelist_get_size(cubes[0]);
01953 
01954     /*-------------------------------------------------------------------- 
01955      * center the cubes within the allocated big cube
01956      * that means define the (0,0) positions of the cubes in the image planes
01957      * to sub-pixel accuracy by using cumoffsetx,y and the reference cube
01958      */
01959     /* position of first reference frame, centered in big cube */
01960     llx0 = olx/2 - ilx/2 ;
01961     lly0 = oly/2 - ily/2 ;
01962 
01963     /*--------------------------------------------------------------------
01964      * go through the frame list and determine the lower left edge position
01965      * of the shifted cubes. Additionnally, the sub-pixel offsets are 
01966      * determined.
01967      */
01968 
01969     llx=cpl_calloc(n_cubes,sizeof(int)); ;
01970     lly=cpl_calloc(n_cubes,sizeof(int)) ;
01971 
01972     sub_offsetx=cpl_calloc(n_cubes,sizeof(float)) ;
01973     sub_offsety=cpl_calloc(n_cubes,sizeof(float)) ;
01974 
01975     for ( i = 0 ; i < n_cubes ; i++ )
01976     {
01977         llx[i] = llx0 - sinfo_new_nint(cumoffsetx[i]) ;
01978         sub_offsetx[i] = (float)sinfo_new_nint(cumoffsetx[i]) - cumoffsetx[i] ;
01979         lly[i] = lly0 - sinfo_new_nint(cumoffsety[i]) ;      
01980         sub_offsety[i] = (float)sinfo_new_nint(cumoffsety[i]) - cumoffsety[i] ;
01981     }
01982 
01983  
01984     /* -------------------------------------------------------------
01985      * shift the cubes according to the computed sub-pixel offsets 
01986      * that means shift the single image planes of each cube
01987      * first determine an interpolation kernel
01988      */
01989     if ( NULL == (kernel = sinfo_generate_interpolation_kernel(kernel_type)))
01990     {
01991         sinfo_msg_warning ("could not generate desired interpolation kernel"
01992                            " or no kernel_typ was given, the default kernel"
01993                            " is used now!") ;
01994     }
01995     /* go through the frame list */
01996     
01997 
01998     tmpcubes=(cpl_imagelist**)cpl_calloc(n_cubes,sizeof(cpl_imagelist*)) ; 
01999 
02000     for ( i = 0 ; i < n_cubes ; i++ )
02001     {
02002         tmpspace = cpl_calloc(ilx, ily*sizeof(pixelvalue)) ;
02003         tmpcubes[i] = cpl_imagelist_new();
02004         
02005         for ( z = 0 ; z < inp ; z++ )
02006         {
02007 
02008             
02009             t_img=sinfo_new_shift_image(cpl_imagelist_get(cubes[i],z), 
02010                                   sub_offsetx[i], sub_offsety[i], kernel);
02011 
02012         if (t_img==NULL)
02013             {
02014                 sinfo_msg_error ("could not shift image plane no %d"
02015                                  " in cube no %d!", z, i) ;
02016                 cpl_imagelist_delete(mergedCube) ;
02017                 cpl_imagelist_delete(mask) ;
02018                 cpl_free(kernel) ;
02019                 return NULL ;
02020             }
02021             cpl_imagelist_set(tmpcubes[i],t_img,z);
02022         }
02023     cpl_free(tmpspace);
02024     }
02025   
02026     /*------------------------------------------------------------------------- 
02027      * Build the mask data cube.
02028      * The mask is 0 where no data is available, otherwise the integration 
02029        time of one frame, respectively the summed integration
02030      * times in the overlapping regions are inserted
02031      */
02032     /* go through the frame list */
02033     for ( i = 0 ; i < n_cubes ; i++ )
02034     {
02035 
02036         /* go through the first image plane of the big data cube */
02037         for ( y = 0 ; y < oly ; y++ )
02038         {
02039             for ( x = 0 ; x < olx ; x++ )
02040             {
02041                 /* find the position of the present cube and 
02042                    go through the single spectra */
02043                 if ( y >= lly[i] && y < lly[i]+ily &&
02044                      x >= llx[i] && x < llx[i]+ilx )
02045                 {
02046                     posx = x - llx[i] ;
02047                     posy = y - lly[i] ;
02048                     for ( z = 0 ; z < onp ; z++ )
02049                     {
02050               t_img=cpl_imagelist_get(tmpcubes[i],z);
02051                       ptdata=cpl_image_get_data_float(t_img);
02052               m_img=cpl_imagelist_get(mask,z);
02053                       pmdata=cpl_image_get_data_float(m_img);
02054                         if (!isnan(ptdata[posx+posy*ilx]) &&
02055                                          ptdata[posx+posy*ilx] != 0.)
02056                         {
02057                             pmdata[x+y*mlx] += exptimes[i] ;
02058                         }
02059                     }
02060                 }
02061             }
02062         }
02063     }
02064 
02065 
02066 
02067 
02068 
02069  
02070     /* calculate a weighted average using the 
02071        exposure time of the single frames 
02072        of the overlapping regions of the cubes */
02073     for ( i = 0 ; i < n_cubes ; i++ )
02074     {
02075 
02076         /* go through the first image plane of the big data cube */
02077         for ( y = 0 ; y < oly ; y++ )
02078         {
02079 
02080             for ( x = 0 ; x < olx ; x++ )
02081             {
02082 
02083                 /* find the position of the present cube 
02084                    and go through the single spectra */
02085                 if ( y >= lly[i] && y < lly[i]+ily &&
02086                      x >= llx[i] && x < llx[i]+ilx )
02087                 {
02088 
02089                     posx = x - llx[i] ;
02090                     posy = y - lly[i] ;
02091                     for ( z = 0 ; z < onp ; z++ )
02092                     {
02093 
02094               t_img=cpl_imagelist_get(tmpcubes[i],z);
02095                       ptdata=cpl_image_get_data_float(t_img);
02096               m_img=cpl_imagelist_get(mask,z);
02097                       pmdata=cpl_image_get_data_float(m_img);
02098                       mlx=cpl_image_get_size_x(m_img);
02099 
02100               o_img=cpl_imagelist_get(mergedCube,z);
02101                       podata=cpl_image_get_data_float(o_img);
02102                       podata[x+y*olx]=0;
02103                         if (!isnan(ptdata[posx+posy*ilx]))
02104                         {
02105                             if (pmdata[x+y*mlx] != 0.)
02106                             {
02107                 /* adjust the intensities to 
02108                                    the first reference cube */
02109                                 weight = exptimes[0] / pmdata[x+y*mlx] ;
02110                             }
02111                             else
02112                             {
02113                                 weight = 0. ;
02114                             }
02115                             podata[x+y*olx] += 
02116                                weight*ptdata[posx+posy*ilx] ;
02117                         }
02118                     }
02119                 }
02120             }
02121         }
02122     }
02123 
02124 
02125     
02126   
02127     /* convert the "free space" in the cube to blank pixels */
02128     /* convert_0_to_ZERO_for_cubes(mergedCube) ; */ 
02129     cpl_free(kernel) ; /* originated by eclise-malloc */
02130     for( i = 0 ; i < n_cubes ; i++ )
02131     {
02132         cpl_imagelist_delete (tmpcubes[i]) ;
02133     }
02134 
02135     cpl_free(tmpcubes); ;
02136     cpl_free(llx); ;
02137     cpl_free(lly) ;
02138 
02139     cpl_free(sub_offsetx) ;
02140     cpl_free(sub_offsety) ;
02141 
02142     return mask ;
02143 }
02144 
02145 
02146 
02147 
02148 
02149 int 
02150 sinfo_new_combine_jittered_cubes_range ( cpl_imagelist ** cubes,
02151                                  cpl_imagelist  * mergedCube,
02152                                  cpl_imagelist  * mask,
02153                                  int        n_cubes,
02154                                  float    * cumoffsetx,
02155                                  float    * cumoffsety,
02156                                  double    * exptimes,
02157                                  char     * kernel_type,
02158                                  const int z_min, const int z_max )
02159 {
02160 
02161     int i,m ;
02162     int x, y, z ; 
02163     int llx0, lly0 ; 
02164     int posx, posy ;
02165     float weight ;
02166     double * kernel ;
02167     /* const int z_siz=z_max-z_min; */
02168     /*cpl_image * shiftedImage ;*/
02169 
02170     cpl_imagelist ** tmpcubes=NULL ;
02171     int* llx=NULL ;
02172     int* lly=NULL ;
02173     float* sub_offsetx=NULL ;
02174     float* sub_offsety=NULL ;
02175 
02176     pixelvalue * tmpspace;
02177 
02178 
02179 
02180   int ilx=0;
02181   int ily=0;
02182   int olx=0;
02183   int oly=0;
02184   int mlx=0;
02185   int mly=0;
02186 
02187 
02188 
02189 
02190   float* podata=NULL;
02191   float* pmdata=NULL;
02192   float* ptdata=NULL;
02193 
02194   cpl_image* i_img=NULL;
02195   cpl_image* o_img=NULL;
02196   cpl_image* m_img=NULL;
02197   cpl_image* t_img=NULL;
02198 
02199      
02200    if ( cubes == NULL )
02201     {
02202         sinfo_msg_error ("no cube list given!") ;
02203         return -1 ;
02204     }
02205     if ( n_cubes <= 0 )
02206     {
02207         sinfo_msg_error ("wrong number of data cubes in list!") ;
02208         return -1 ;
02209     }
02210     if ( cumoffsetx == NULL || cumoffsety == NULL )
02211     {
02212         sinfo_msg_error ("no cumoffsetx/y given!") ;
02213         return -1;
02214     }
02215     if ( exptimes == NULL )
02216     {
02217         sinfo_msg_error ("no exposure time array given!") ;
02218         return -1 ;
02219     }
02220     
02221     /*
02222     for ( i = 0 ; i < n_cubes ; i++ )
02223     {
02224     tmpcubes[i] = sinfo_copy_cube_range(cubes[i],z_min,z_max);
02225     }
02226     */
02227     
02228     o_img=cpl_imagelist_get(mergedCube,z_min);
02229     olx=cpl_image_get_size_x(o_img);
02230     oly=cpl_image_get_size_y(o_img);
02231     i_img=cpl_imagelist_get(cubes[0],0);
02232     ilx=cpl_image_get_size_x(i_img);
02233     ily=cpl_image_get_size_y(i_img);
02234     mlx=olx;
02235     mly=oly;
02236 
02237 
02238     /*-------------------------------------------------------------------- 
02239      * center the cubes within the allocated big cube
02240      * that means define the (0,0) positions of the cubes in the image planes
02241      * to sub-pixel accuracy by using cumoffsetx,y and the reference cube
02242      */
02243     /* position of first reference frame, centered in big cube */
02244     llx0 = olx/2 - ilx/2 ;
02245     lly0 = oly/2 - ily/2 ;
02246 
02247     /*--------------------------------------------------------------------
02248      * go through the frame list and determine the lower left edge position
02249      * of the shifted cubes. Additionnally, the sub-pixel offsets are 
02250      * determined.
02251      */
02252 
02253 
02254     llx=cpl_calloc(n_cubes,sizeof(int)) ;
02255     lly=cpl_calloc(n_cubes,sizeof(int)) ;
02256     sub_offsetx=cpl_calloc(n_cubes,sizeof(float)) ;
02257     sub_offsety=cpl_calloc(n_cubes,sizeof(float)) ;
02258 
02259     for ( i = 0 ; i < n_cubes ; i++ )
02260     {
02261         llx[i] = llx0 - sinfo_new_nint(cumoffsetx[i]) ;
02262         sub_offsetx[i] = (float)sinfo_new_nint(cumoffsetx[i]) - cumoffsetx[i] ;
02263         lly[i] = lly0 - sinfo_new_nint(cumoffsety[i]) ;      
02264         sub_offsety[i] = (float)sinfo_new_nint(cumoffsety[i]) - cumoffsety[i] ;
02265      }
02266 
02267 
02268     /* -------------------------------------------------------------
02269      * shift the cubes according to the computed sub-pixel offsets 
02270      * that means shift the single image planes of each cube
02271      * first determine an interpolation kernel
02272      */
02273     if ( NULL == (kernel = sinfo_generate_interpolation_kernel( kernel_type )) )
02274     {
02275         sinfo_msg_warning ("could not generate desired interpolation kernel"
02276                            " or no kernel_typ was given, the default kernel"
02277                            " is used now!") ;
02278     }
02279     /* go through the frame list */
02280     
02281     tmpcubes=(cpl_imagelist**)cpl_calloc(n_cubes,sizeof(cpl_imagelist*)) ;
02282 
02283     for ( i = 0 ; i < n_cubes ; i++ )
02284     {
02285         tmpspace = cpl_calloc(ilx, ily*sizeof(pixelvalue)) ;
02286     /* 
02287         tmpcubes[i] = sinfo_new_cube(cubes[i]->lx,cubes[i]->ly,z_siz);
02288     */
02289         tmpcubes[i] = cpl_imagelist_new();
02290         for ( z = z_min, m=0 ; z< z_max ; z++, m++ )
02291         {
02292       /*
02293             sinfo_shiftImageinCube( cubes[i]->plane[z], 
02294                                     sub_offsetx[i], 
02295                                     sub_offsety[i], 
02296                                     kernel, 
02297                                     tmpcubes[i]->plane[m], 
02298                                     tmpspace ); 
02299       */
02300 
02301             t_img=sinfo_new_shift_image(cpl_imagelist_get(cubes[i],z), 
02302                                   sub_offsetx[i], sub_offsety[i], kernel);
02303 
02304 
02305         if (t_img==NULL)
02306             {
02307                 sinfo_msg_error ("could not shift image plane no %d"
02308                                  "in cube no %d!", z, i) ;
02309                 cpl_imagelist_delete(mergedCube) ;
02310                 cpl_imagelist_delete(mask) ;
02311                 cpl_free(kernel) ;
02312                 return -1 ;
02313             }
02314             cpl_imagelist_set(tmpcubes[i],t_img,m);
02315             m_img=cpl_image_new(mlx,mly,CPL_TYPE_FLOAT);
02316             cpl_imagelist_set(mask,m_img,z);
02317         }
02318     cpl_free(tmpspace);
02319     }
02320   
02321     /*------------------------------------------------------------------------- 
02322      * Build the mask data cube.
02323      * The mask is 0 where no data is available, otherwise the 
02324        integration time of 
02325      * one frame, respectively the summed integration
02326      * times in the overlapping regions are inserted
02327      */
02328     /* go through the frame list */
02329     for ( z = z_min, m=0 ; z < z_max ; z++, m++ )
02330       {
02331 
02332         /* go through the first image plane of the big data cube */
02333         for ( y = 0 ; y < oly ; y++ )
02334         {
02335              for ( x = 0 ; x < olx ; x++ )
02336             {
02337 
02338 
02339           for ( i = 0 ; i < n_cubes ; i++ )
02340         {
02341 
02342           i_img=cpl_imagelist_get(cubes[i],0);
02343           ilx=cpl_image_get_size_x(i_img);
02344           ily=cpl_image_get_size_y(i_img);
02345 
02346 
02347                 /* find the position of the present cube and go 
02348                    through the single spectra */
02349                 if ( y >= lly[i] && y < lly[i]+ily &&
02350                      x >= llx[i] && x < llx[i]+ilx )
02351                 {
02352                     posx = x - llx[i] ;
02353                     posy = y - lly[i] ;
02354 
02355 
02356               t_img=cpl_imagelist_get(tmpcubes[i],m);
02357               ptdata=cpl_image_get_data_float(t_img);
02358               m_img=cpl_imagelist_get(mask,z);
02359               pmdata=cpl_image_get_data_float(m_img);
02360                       mlx=cpl_image_get_size_x(m_img);
02361 
02362                          if (!isnan(ptdata[posx+posy*ilx]) &&
02363                                     ptdata[posx+posy*ilx] != 0.)
02364                         {
02365                              pmdata[x+y*mlx] += (float)exptimes[i] ;
02366             } else if (isnan(ptdata[posx+posy*ilx])) {
02367               sinfo_msg_debug("ptdata %d, %d, %d is NAN\t",x,y,z); 
02368             } else if (ptdata[posx+posy*ilx] == 0.) {
02369               sinfo_msg_debug("ptdata %d, %d, %d is 0\t",x,y,z); 
02370             }
02371     
02372                 } else {
02373           sinfo_msg_debug("point %d, %d, %d outside range\n",x,y,z); 
02374         }
02375         }
02376         }
02377     }
02378       }
02379 
02380     /* calculate a weighted average using the exposure time of the 
02381        single frames of the overlapping regions of the cubes */
02382     for ( z = z_min, m = 0 ; z < z_max ; z++, m++ )
02383       {
02384     o_img=cpl_imagelist_get(mergedCube,z);
02385     podata=cpl_image_get_data_float(o_img);
02386     m_img=cpl_imagelist_get(mask,z);
02387     pmdata=cpl_image_get_data_float(m_img);
02388     mlx=cpl_image_get_size_x(m_img);
02389 
02390         /* go through the first image plane of the big data cube */
02391         for ( y = 0 ; y < oly ; y++ )
02392         {
02393             for ( x = 0 ; x < olx ; x++ )
02394             {
02395 
02396                 /* find the position of the present cube and 
02397                    go through the single spectra */
02398 
02399           for ( i = 0 ; i < n_cubes ; i++ )
02400         {
02401 
02402                 if ( y >= lly[i] && y < lly[i]+ily &&
02403                      x >= llx[i] && x < llx[i]+ilx )
02404                 {
02405                     posx = x - llx[i] ;
02406                     posy = y - lly[i] ;
02407 
02408 
02409 
02410 
02411               t_img=cpl_imagelist_get(tmpcubes[i],m);
02412               ptdata=cpl_image_get_data_float(t_img);
02413               /* To prevent black regions in peculiar batterfly cases
02414               podata[x+y*olx]=0;
02415               */
02416                         if (!isnan(ptdata[posx+posy*ilx]))
02417                         {
02418                             if (pmdata[x+y*mlx] != 0.)
02419                             {
02420                 /* adjust the intensities to the 
02421                                    first reference cube */
02422                                 weight = exptimes[0] / pmdata[x+y*mlx] ;
02423                             }
02424                             else
02425                             {
02426                                 weight = 0. ;
02427                             }
02428                             podata[x+y*olx] += weight*ptdata[posx+posy*ilx] ;
02429 
02430                         }
02431                     }
02432                 }
02433             }
02434         }
02435     }
02436  
02437 
02438     /* convert the "free space" in the cube to blank pixels */
02439      /* convert_0_to_ZERO_for_cubes(mergedCube) ; */
02440     cpl_free(kernel) ; /* originated by eclise-malloc */
02441     for( i = 0 ; i < n_cubes ; i++ )
02442     {
02443         cpl_imagelist_delete (tmpcubes[i]) ;
02444     }
02445 
02446  
02447     cpl_free(tmpcubes) ;
02448     cpl_free(llx) ;
02449     cpl_free(lly) ;
02450     cpl_free(sub_offsetx) ;
02451     cpl_free(sub_offsety) ;
02452 
02453      return 0 ;
02454 }
02455 
02456 
02457 
02458 int 
02459 sinfo_new_combine_jittered_cubes_thomas_range(cpl_imagelist ** cubes,
02460                     cpl_imagelist  * mergedCube,
02461                     cpl_imagelist  * mask,
02462                     int        n_cubes,
02463                     float    * cumoffsetx,
02464                     float    * cumoffsety,
02465                     double    * exptimes,
02466                     char     * kernel_type,
02467                     const int z_min,
02468                     const int z_max,
02469                                     const double kappa )
02470 {
02471 
02472     int i ;
02473     int x, y, z ; 
02474     int llx0, lly0 ; 
02475     int posx, posy ;
02476     /* float weight ; */
02477 
02478     double med=0;
02479     double sig=0;
02480     double avg=0;
02481     int nc=0;
02482 
02483 
02484     double * kernel ;
02485 
02486     cpl_vector* val=NULL;
02487     cpl_vector* msk=NULL;
02488     int nclip=0;
02489     int ks=0;
02490     int ovr=0;
02491 
02492     double msk_sum=0;
02493     double val_msk_sum=0;
02494     /*cpl_image * shiftedImage ;*/
02495 
02496 
02497     int* llx=NULL;
02498     int* lly=NULL ;
02499     float* sub_offsetx=NULL ;
02500     float* sub_offsety=NULL ;
02501     cpl_imagelist ** tmpcubes=NULL ;
02502 
02503     pixelvalue * tmpspace;
02504     const int z_siz=z_max-z_min;   
02505     int m=0;
02506 
02507 
02508 
02509 
02510 
02511 
02512   int ilx=0;
02513   int ily=0;
02514   int olx=0;
02515   int oly=0;
02516   int mlx=0;
02517   int mly=0;
02518 
02519 
02520   int inp=0;
02521   int onp=0;
02522 
02523 
02524   float* podata=NULL;
02525   float* pmdata=NULL;
02526   float* ptdata=NULL;
02527   float* pvdata=NULL;
02528 
02529   cpl_image* i_img=NULL;
02530   cpl_image* o_img=NULL;
02531   cpl_image* m_img=NULL;
02532   cpl_image* t_img=NULL;
02533   cpl_image* v_img=NULL;
02534 
02535 
02536   
02537     if ( cubes == NULL )
02538     {
02539         sinfo_msg_error ("no cube list given!") ;
02540         return -1 ;
02541     }
02542     if ( n_cubes <= 0 )
02543     {
02544         sinfo_msg_error ("wrong number of data cubes in list!") ;
02545         return -1 ;
02546     }
02547     if ( cumoffsetx == NULL || cumoffsety == NULL )
02548     {
02549         sinfo_msg_error ("no cumoffsetx/y given!") ;
02550         return -1 ;
02551     }
02552     if ( exptimes == NULL )
02553     {
02554         sinfo_msg_error ("no exposure time array given!") ;
02555         return -1 ;
02556     }
02557     if (z_siz <= 0 )
02558     {
02559         sinfo_msg_error ("z_max <= z_min given!") ;
02560         return -1 ;
02561     }
02562     
02563     i_img=cpl_imagelist_get(cubes[0],0);
02564     o_img=cpl_imagelist_get(mergedCube,0);
02565     ilx=cpl_image_get_size_x(i_img);
02566     ily=cpl_image_get_size_y(i_img);
02567     olx=cpl_image_get_size_x(o_img);
02568     oly=cpl_image_get_size_y(o_img);
02569     mlx=olx;
02570     mly=oly;
02571 
02572   
02573     /*-------------------------------------------------------------------- 
02574      * center the cubes within the allocated big cube
02575      * that means define the (0,0) positions of the cubes in the image planes
02576      * to sub-pixel accuracy by using cumoffsetx,y and the reference cube
02577      */
02578     /* position of first reference frame, centered in big cube */
02579     llx0 = olx/2 - ilx/2 ;
02580     lly0 = oly/2 - ily/2 ;
02581 
02582     /*--------------------------------------------------------------------
02583      * go through the frame list and determine the lower left edge position
02584      * of the shifted cubes. Additionnally, the sub-pixel offsets are 
02585      * determined.
02586      */
02587 
02588     llx=cpl_calloc(n_cubes,sizeof(int));
02589     lly=cpl_calloc(n_cubes,sizeof(int)) ;
02590     sub_offsetx=cpl_calloc(n_cubes,sizeof(float)) ;
02591     sub_offsety=cpl_calloc(n_cubes,sizeof(float)) ;
02592 
02593 
02594 
02595     for ( i = 0 ; i < n_cubes ; i++ )
02596     {
02597         llx[i] = llx0 - sinfo_new_nint(cumoffsetx[i]) ;
02598         sub_offsetx[i] = (float)sinfo_new_nint(cumoffsetx[i]) - cumoffsetx[i] ;
02599         lly[i] = lly0 - sinfo_new_nint(cumoffsety[i]) ;      
02600         sub_offsety[i] = (float)sinfo_new_nint(cumoffsety[i]) - cumoffsety[i] ;
02601     }
02602 
02603   
02604     /* -------------------------------------------------------------
02605      * shift the cubes according to the computed sub-pixel offsets 
02606      * that means shift the single image planes of each cube
02607      * first determine an interpolation kernel
02608      */
02609     if ( NULL == (kernel = sinfo_generate_interpolation_kernel(kernel_type)) )
02610     {
02611         sinfo_msg_warning ("could not generate desired interpolation kernel"
02612                            "or no kernel_typ was given, the default kernel"
02613                            "is used now!") ;
02614     }
02615     /* go through the frame list */
02616     
02617     tmpcubes=(cpl_imagelist**)cpl_calloc(n_cubes,sizeof(cpl_imagelist*)) ;
02618 
02619     for ( i = 0 ; i < n_cubes ; i++ )
02620     {
02621       i_img=cpl_imagelist_get(cubes[i],0);
02622       ilx=cpl_image_get_size_x(i_img);
02623       ily=cpl_image_get_size_y(i_img);
02624       inp=cpl_imagelist_get_size(cubes[i]);
02625       tmpspace = cpl_calloc(ilx, ily*sizeof(pixelvalue)) ;
02626       tmpcubes[i]=cpl_imagelist_new();
02627       /* was
02628     tmpcubes[i]=sinfo_new_cube(cubes[i]->lx,cubes[i]->ly,cubes[i]->np);
02629       */
02630         for ( z = z_min, m=0 ; z < z_max ; z++, m++ )
02631         {
02632               t_img=sinfo_new_shift_image(cpl_imagelist_get(cubes[i],z), 
02633                                           sub_offsetx[i], 
02634                                           sub_offsety[i], 
02635                                           kernel);
02636 
02637         if (t_img==NULL)
02638             {
02639                 sinfo_msg_error("could not shift image plane no %d "
02640                                 "in cube no %d!", z, i) ;
02641                 cpl_imagelist_delete(mergedCube) ;
02642                 cpl_imagelist_delete(mask) ;
02643                 cpl_free(kernel) ;
02644                 return -1 ;
02645             }
02646 
02647             cpl_imagelist_set(tmpcubes[i],t_img,m);
02648             m_img=cpl_image_new(mlx,mly,CPL_TYPE_FLOAT);
02649             cpl_imagelist_set(mask,m_img,z);
02650         }
02651     cpl_free(tmpspace);
02652     /*tmpcubes[i]=sinfo_copy_cube(cubes[i]);*/
02653     }
02654 
02655     /*------------------------------------------------------------------------- 
02656      * Build the mask data cube.
02657      * The mask is 0 where no data is available, otherwise the integration 
02658        time of one frame, respectively the summed integration
02659      * times in the overlapping regions are inserted
02660      */
02661     /* go through the frame list */
02662 
02663 
02664     o_img=cpl_imagelist_get(mergedCube,0);
02665     olx=cpl_image_get_size_x(o_img);
02666     oly=cpl_image_get_size_y(o_img);
02667     onp=cpl_imagelist_get_size(mergedCube);
02668  
02669     for ( i = 0 ; i < n_cubes ; i++ )
02670     {
02671 
02672 
02673       i_img=cpl_imagelist_get(cubes[i],0);
02674       ilx=cpl_image_get_size_x(i_img);
02675       ily=cpl_image_get_size_y(i_img);
02676       inp=cpl_imagelist_get_size(cubes[i]);
02677 
02678         /* go through the first image plane of the big data cube */
02679         for ( y = 0 ; y < oly ; y++ )
02680         {
02681 
02682             for ( x = 0 ; x < olx ; x++ )
02683             {
02684 
02685                 /* find the position of the present cube and go 
02686                    through the single spectra */
02687                 if ( y >= lly[i] && y < lly[i]+ily &&
02688                      x >= llx[i] && x < llx[i]+ilx )
02689                 {
02690                     posx = x - llx[i] ;
02691                     posy = y - lly[i] ;
02692 
02693                     for ( z = z_min,m=0 ; z < z_max ; z++,m++ )
02694                     {
02695               t_img=cpl_imagelist_get(tmpcubes[i],m);
02696               ptdata=cpl_image_get_data_float(t_img);
02697               m_img=cpl_imagelist_get(mask,z);
02698               pmdata=cpl_image_get_data_float(m_img);
02699                       mlx=cpl_image_get_size_x(m_img);
02700 
02701                         if (!isnan(ptdata[posx+posy*ilx]) &&
02702                                          ptdata[posx+posy*ilx] != 0.)
02703                         {
02704                            pmdata[x+y*mlx] += (float)exptimes[i]  ;
02705                         }
02706                     }
02707                 }
02708             }
02709         }
02710     }
02711  
02712     m=0;
02713     for ( z = z_min; z < z_max ; z++ ) {
02714       m_img=cpl_imagelist_get(mask,z);
02715       pmdata=cpl_image_get_data_float(m_img);
02716       o_img=cpl_imagelist_get(mergedCube,z);
02717       podata=cpl_image_get_data_float(o_img);
02718       mlx=cpl_image_get_size_x(m_img);
02719       mly=cpl_image_get_size_y(m_img);
02720       /*     
02721       podata[x+y*olx]=0;
02722       */
02723       /* go through the first image plane of the big data cube */
02724       for ( y = 0 ; y < oly ; y++ ) {
02725     for ( x = 0 ; x < olx ; x++ ) {
02726       avg=0;
02727       nc=0;
02728           /* computes nc */
02729       for ( i = 0 ; i < n_cubes ; i++ ) {
02730             t_img=cpl_imagelist_get(tmpcubes[i],m);
02731             ptdata=cpl_image_get_data_float(t_img);
02732         if ( y >= lly[i] && y < lly[i]+ily &&
02733          x >= llx[i] && x < llx[i]+ilx ) {
02734           posx = x - llx[i] ;
02735           posy = y - lly[i] ;
02736           if (!isnan(ptdata[posx+posy*ilx]) &&
02737                  ptdata[posx+posy*ilx] != 0.) {
02738         nc++;
02739           }
02740         }
02741       }
02742           if( nc > 0 ) {
02743         msk=cpl_vector_new(n_cubes);
02744         for (i=0;i<n_cubes;i++) {
02745           cpl_vector_set(msk,i,1);
02746         }
02747 
02748         /* k-s clipping */
02749         nclip=0;
02750 
02751 
02752         for (ks=0;ks<nc;ks++) {
02753           sig=0;
02754           med=0;
02755           ovr=0;
02756           if(nc-nclip >0) {
02757         val=cpl_vector_new(nc-nclip);
02758           }
02759           /* fill val */
02760           for ( i = 0 ; i < n_cubes ; i++ ) {
02761         t_img=cpl_imagelist_get(tmpcubes[i],m);
02762         ptdata=cpl_image_get_data_float(t_img);
02763             if ( y >= lly[i] && y < lly[i]+ily &&
02764              x >= llx[i] && x < llx[i]+ilx ) {
02765           posx = x - llx[i] ;
02766           posy = y - lly[i] ;
02767               if (!isnan(ptdata[posx+posy*ilx]) &&
02768                            ptdata[posx+posy*ilx] != 0. &&
02769                                   (cpl_vector_get(msk,i) != 0)) {
02770             cpl_vector_set(val,ovr,
02771                            (double)ptdata[posx+posy*ilx]);
02772             ovr++;
02773           }
02774         }
02775           }
02776 
02777           /* get avg, med, sig */
02778           if(ovr>0) {
02779         avg=cpl_vector_get_mean(val);
02780         med=cpl_vector_get_median(val);
02781         if(ovr>1) {
02782           sig=cpl_vector_get_stdev(val);
02783         } else {
02784           sig=0;
02785         }
02786         cpl_vector_delete(val);
02787           }
02788 
02789           for ( i = 0 ; i < n_cubes ; i++ ) {
02790         t_img=cpl_imagelist_get(tmpcubes[i],m);
02791         ptdata=cpl_image_get_data_float(t_img);
02792         /* Do k-s clipping at each pixel */
02793             if ( y >= lly[i] && y < lly[i]+ily &&
02794              x >= llx[i] && x < llx[i]+ilx ) {
02795           posx = x - llx[i] ;
02796           posy = y - lly[i] ;
02797               if (!isnan(ptdata[posx+posy*ilx]) &&
02798                      ptdata[posx+posy*ilx] != 0. &&
02799                             (cpl_vector_get(msk,i) != 0)) {  
02800             if(abs((ptdata[posx+posy*ilx]-med))> kappa*sig) {
02801                     ptdata[posx+posy*ilx]=0;
02802                       pmdata[x+y*mlx] -= exptimes[i]  ;
02803               cpl_vector_set(msk,i,0);
02804               nclip++;
02805             }
02806           }
02807         }
02808           }
02809         } /* end of k-s clipping */
02810 
02811         msk_sum=0;
02812         val_msk_sum=0;
02813         for ( i = 0 ; i < n_cubes ; i++ ) {
02814         v_img=cpl_imagelist_get(tmpcubes[i],m);
02815         pvdata=cpl_image_get_data_float(v_img);
02816           /* computes sky at each point */
02817           if ( y >= lly[i] && y < lly[i]+ily &&
02818            x >= llx[i] && x < llx[i]+ilx ) {
02819             posx = x - llx[i] ;
02820             posy = y - lly[i] ;
02821             if (!isnan(pvdata[posx+posy*ilx]) &&
02822                    pvdata[posx+posy*ilx] != 0. &&
02823                           (cpl_vector_get(msk,i) != 0)) {  
02824                   msk_sum+=pmdata[x+y*mlx];
02825                   val_msk_sum+=pvdata[posx+posy*ilx]*
02826                        pmdata[x+y*mlx];
02827         }
02828           }
02829         }
02830         podata[x+y*olx]=val_msk_sum/msk_sum;
02831         cpl_vector_delete(msk);
02832       } /* end check if overlap nc >0  */
02833     } /* end loop over x */
02834       } /* end loop over y */
02835       m++;
02836     } /* end loop over z */
02837 
02838 
02839     /* calculate a weighted average using the exposure time of the 
02840        single frames of the overlapping regions of the cubes */
02841     /*
02842     for ( i = 0 ; i < n_cubes ; i++ )
02843     {
02844 
02845       // go through the first image plane of the big data cube 
02846         for ( y = 0 ; y < oly ; y++ )
02847         {
02848             for ( x = 0 ; x < olx ; x++ )
02849             {
02850           // find the position of the present cube and 
02851               // go through the single spectra 
02852                 if ( y >= lly[i] && y < lly[i]+ily &&
02853                      x >= llx[i] && x < llx[i]+ilx )
02854                 {
02855                     posx = x - llx[i] ;
02856                     posy = y - lly[i] ;
02857                     for ( z = z_min ; z < z_max ; z++ )
02858                     {
02859                 t_img=cpl_imagelist_get(tmpcubes[i],z);
02860                 ptdata=cpl_image_get_data_float(t_img);
02861                         m_img=cpl_imagelist_get(mask,z);
02862                         pmdata=cpl_image_get_data_float(m_img);
02863                         o_img=cpl_imagelist_get(mergedCube,z);
02864                         podata=cpl_image_get_data_float(o_img);
02865                         mlx=cpl_image_get_size_x(m_img);
02866                         mly=cpl_image_get_size_y(m_img);
02867                         if (!isnan(ptdata[posx+posy*ilx]))
02868                         {
02869                             if (pmdata[x+y*mask->lx] != 0.)
02870                             {
02871                   // adjust the intensities to the first reference cube 
02872                                 weight = exptimes[0] / pmdata[x+y*mlx] ;
02873                             }
02874                             else
02875                             {
02876                                 weight = 0. ;
02877                             }
02878                             
02879                             podata[x+y*olx] += weight*ptdata[posx+posy*ilx] ;
02880 
02881                         }
02882                     }
02883                 }
02884             }
02885         }
02886     }
02887     */
02888 
02889     /* convert the "free space" in the cube to blank pixels */
02890     /* convert_0_to_ZERO_for_cubes(mergedCube) ; */
02891       /* convert_0_to_ZERO_for_cubes(mergedSky) ; */
02892     cpl_free(kernel) ; /* originated by eclise-malloc */
02893     for( i = 0 ; i < n_cubes ; i++ )
02894     {
02895         cpl_imagelist_delete (tmpcubes[i]) ;
02896     }
02897 
02898     cpl_free(tmpcubes);
02899     cpl_free(llx);
02900     cpl_free(lly) ;
02901     cpl_free(sub_offsetx) ;
02902     cpl_free(sub_offsety) ;
02903 
02904 
02905     return 0 ;
02906 }
02907 
02908 /*---------------------------------------------------------------------------
02909    Function :   sinfo_new_interpol_cube_simple()
02910    In       :   cube: 1 allocated cube
02911                         badcube: bad pixel mask cube (0: bad, 1: good pixel)
02912                         maxdist: maximal pixel distance from bad pixel 
02913                                  to search for good pixels, don't make this
02914                                  value too big!
02915    Out      :   interpolated cube, and corrected bad pixel mask cube
02916    Job      :   interpolates bad pixel of an object cube if a bad pixel
02917                         mask cube is available by using the nearest neighbors
02918                         in 3 dimensions.
02919  ---------------------------------------------------------------------------*/
02920 cpl_imagelist * 
02921 sinfo_new_interpol_cube_simple( cpl_imagelist * cube, 
02922                       cpl_imagelist * badcube,
02923                       int       maxdist )
02924 {
02925   cpl_imagelist * intercube ;
02926   float*     goodNeighbors=NULL ;
02927   int z, row, col ;
02928   int nx, ny, nz ;
02929   int llx, lly, llz ;
02930   int zi, coli, rowi ;
02931   int n ;
02932 
02933 
02934 
02935 
02936   int clx=0;
02937   int cly=0;
02938   int blx=0;
02939   int bly=0;
02940 
02941   int cnp=0;
02942 
02943 
02944   float* pbdata=NULL;
02945   float* pidata=NULL;
02946   float* pbzidata=NULL;
02947   float* pczidata=NULL;
02948 
02949   cpl_image* c_img=NULL;
02950   cpl_image* b_img=NULL;
02951   cpl_image* i_img=NULL;
02952 
02953   cpl_image* bzi_img=NULL;
02954   cpl_image* czi_img=NULL;
02955 
02956 
02957 
02958     if ( cube == NULL || badcube == NULL )
02959     {
02960         sinfo_msg_error("no cube given!") ;
02961         return NULL ;
02962     }
02963     if ( maxdist < 1 )
02964     {
02965         sinfo_msg_error("wrong maxrad given!") ;
02966         return NULL ;
02967     }
02968     intercube = cpl_imagelist_duplicate(cube) ;
02969 
02970     goodNeighbors=cpl_calloc((2*maxdist+1)*(2*maxdist+1)*(2*maxdist+1) -1,
02971                              sizeof(float)) ;
02972 
02973     cnp=cpl_imagelist_get_size(cube);
02974     for ( z = 0 ; z < cnp ; z++ )
02975     {
02976       b_img=cpl_imagelist_get(badcube,z);
02977       i_img=cpl_imagelist_get(intercube,z);
02978       pbdata=cpl_image_get_data_float(b_img);
02979       pidata=cpl_image_get_data_float(i_img);
02980       blx=cpl_image_get_size_x(b_img);
02981       bly=cpl_image_get_size_y(b_img);
02982       
02983       c_img=cpl_imagelist_get(cube,z);
02984       clx=cpl_image_get_size_x(c_img);
02985       cly=cpl_image_get_size_y(c_img);
02986 
02987         for ( row = 0 ; row < cly ; row++ )
02988         {
02989             for ( col = 0 ; col < clx ; col++ )
02990             {
02991                 if ( pbdata[col+row*clx] == 0 )
02992                 {
02993                     /* determine the lower left sinfo_edge of the cube */
02994                     llx = col - maxdist ;
02995                     nx = 2*maxdist +1 ;
02996                     if (llx < 0) 
02997                     {
02998                        nx += llx ;
02999                        llx = 0 ; 
03000                     }
03001                     if ( llx + nx > clx )
03002                     {
03003                         nx -= (llx + nx - clx) ;
03004                     }
03005  
03006                     lly = row - maxdist ;
03007                     ny = 2*maxdist +1 ;
03008                     if (lly < 0) 
03009                     {
03010                        ny += lly ;
03011                        lly = 0 ; 
03012                     }
03013                     if ( lly + ny > cly )
03014                     {
03015                         ny -= (lly + ny - cly) ;
03016                     }
03017                     
03018                     llz = z - maxdist ;
03019                     nz = 2*maxdist +1 ;
03020                     if (llz < 0) 
03021                     {
03022                        nz += llz ;
03023                        llz = 0 ; 
03024                     }
03025                     if ( llz + nz > cnp )
03026                     {
03027                         nz -= (llz + nz - cnp) ;
03028                     }
03029                     n = 0 ;
03030                     for ( zi = llz ; zi < llz+nz ; zi++ )
03031                     {
03032               bzi_img=cpl_imagelist_get(badcube,zi);
03033               czi_img=cpl_imagelist_get(cube,zi);
03034                       pbzidata=cpl_image_get_data_float(bzi_img);
03035                       pczidata=cpl_image_get_data_float(czi_img);
03036 
03037                         for ( rowi = lly ; rowi < lly+ny ; rowi++ )
03038                         {
03039                             for ( coli = llx ; coli < llx+nx ; coli++ )
03040                             {
03041                                 if ( pbzidata[coli+rowi*blx] == 1 )
03042                                 {
03043                   goodNeighbors[n] = pczidata[coli+rowi*clx] ;
03044                   n++ ;
03045                                 }
03046                             }
03047                         }
03048                     }
03049                     if ( n > 0 )
03050                     {
03051                         pidata[col+row*clx]=sinfo_new_median(goodNeighbors,n);
03052                         pbdata[col+row*clx]=1 ;
03053                     }
03054                     else
03055                     {
03056                         continue ;
03057                     }
03058                 }
03059             }
03060         }
03061     }
03062     cpl_free(goodNeighbors) ;
03063     return intercube ;
03064 }
03065 
03066 /*---------------------------------------------------------------------------
03067    Function     :       sinfo_new_combine_cubes()
03068    In           :       cubes: list of jittered cubes to mosaic
03069                         mergedCube: resulting merged cube containing the 
03070                                       jittered cubes
03071                         n_cubes: number of cubes in the list to merge
03072                         cumoffsetx,y: array of relative x, y pixel offsets
03073                                       with respect to the first frame in the 
03074                                       same sequence as the cube list.
03075                         factor:      sigma factor beyond which pixels are
03076                                      thrown away.
03077                         kernel_type: the name of the interpolation kernel
03078                                      that you want to generate using the 
03079                                      eclipse routine 
03080                                      sinfo_generate_interpolation_kernel()
03081                                      Supported kernels are:
03082                                      NULL:      default kernel, currently tanh
03083                                      "default": dito
03084                                      "tanh":    Hyperbolic tangent
03085                                      "sinc2":   Square sinc
03086                                      "lanczos": Lanczos2 kernel
03087                                      "hamming": Hamming kernel
03088                                      "hann":    Hann kernel
03089    Out          :       mask: cube of the same size as combinedCube 
03090                               containing 0 for blank (ZERO pixels) and
03091                               n used pixels for overlapping regions
03092                         mergedCube: final data cube containing the jittered 
03093                                     cubes
03094    Job          :merges jittered data cubes to one bigger cube
03095                  by taking the sinfo_new_median in each pixel and throw 
03096                  away the high deviation pixels (bigger than factor * sigma)
03097                  The x, y size of the final data
03098                  cube is user given, and should be between 32 and 64 
03099                  pixels, while the relative pixel-offset (sub-pixel
03100                  accuracy) of the single cubes with respect to the
03101                  first cube in the list is read from the SEQ CUMOFFSETX,Y
03102                  fits header keyword. 
03103  ---------------------------------------------------------------------------*/
03104 cpl_imagelist * 
03105 sinfo_new_combine_cubes ( cpl_imagelist ** cubes,
03106                          cpl_imagelist  * mergedCube,
03107                          int        n_cubes,
03108                          float    * cumoffsetx,
03109                          float    * cumoffsety,
03110                          float      factor,
03111                          char     * kernel_type )
03112 {
03113   int i=0 ;
03114   int x=0;
03115   int y=0;
03116   int z=0; 
03117   int llx0=0;
03118   int lly0=0; 
03119   int posx=0;
03120   int posy=0;
03121   cpl_imagelist * mask=NULL ;
03122   double * kernel=NULL ;
03123   cpl_image * shiftedImage=NULL ;
03124   int n=0;
03125   int ns=0;
03126   double sum=0; 
03127   double sum2=0;
03128   double mean=0;
03129   double sigma=0;
03130 
03131   cpl_imagelist ** tmpcubes=NULL ;
03132 
03133   int* llx=NULL ;
03134   int* lly=NULL ;
03135 
03136   float* sub_offsetx=NULL ;
03137   float* sub_offsety=NULL ;
03138   float* cubedata=NULL ;
03139 
03140   int mlx=0;
03141   int mly=0;
03142   int clx=0;
03143   int cly=0;
03144   int mnp=0;
03145   int cnp=0;
03146 
03147 
03148   float* ptdata=NULL;
03149   float* podata=NULL;
03150   float* pmdata=NULL;
03151 
03152   cpl_image* tmp_img=NULL;
03153   cpl_image* o_img=NULL;
03154   cpl_image* m_img=NULL;
03155   cpl_image* c_img=NULL;
03156   cpl_image* t_img=NULL;
03157 
03158 
03159 
03160 
03161 
03162   if ( cubes == NULL )
03163     {
03164         sinfo_msg_error ("no cube list given!") ;
03165         return NULL ;
03166     }
03167 
03168   if ( n_cubes <= 0 )
03169     {
03170         sinfo_msg_error ("wrong number of data cubes in list!") ;
03171         return NULL ;
03172     }
03173 
03174   if ( cumoffsetx == NULL || cumoffsety == NULL )
03175     {
03176         sinfo_msg_error ("no cumoffsetx/y given!") ;
03177         return NULL ;
03178     }
03179 
03180   if ( factor <= 0. )
03181     {
03182         sinfo_msg_error ("wrong factor given!") ;
03183         return NULL ;
03184     }
03185   
03186   m_img=cpl_imagelist_get(mergedCube,0);
03187   mlx=cpl_image_get_size_x(m_img);
03188   mly=cpl_image_get_size_y(m_img);
03189   cnp=cpl_imagelist_get_size(cubes[0]);
03190   c_img=cpl_imagelist_get(cubes[0],0);
03191   clx=cpl_image_get_size_x(c_img);
03192   cly=cpl_image_get_size_y(c_img);
03193 
03194 
03195   tmpcubes=(cpl_imagelist**)cpl_calloc(n_cubes,sizeof(cpl_imagelist*)) ;
03196 
03197         /* allocation for a cube structure without the image planes  */
03198     /*
03199    for ( i = 0 ; i < n_cubes ; i++ )
03200     {
03201          tmpcubes[i] = (cpl_imagelist*)cpl_malloc(sizeof(cpl_imagelist)) ;
03202         tmpcubes[i]->plane = (cpl_image**)cpl_calloc(cubes[0]->np , 
03203                               sizeof(cpl_image*)) ;
03204 
03205         tmpcubes[i]->lx = cubes[0]->lx ;
03206         tmpcubes[i]->ly = cubes[0]->ly ;
03207         tmpcubes[i]->np = cubes[0]->np ;
03208         tmpcubes[i]->nbpix = (ulong32)cubes[0]->lx * 
03209                              (ulong32)cubes[0]->ly * 
03210                              (ulong32)cubes[0]->np ;
03211         tmpcubes[i]->history = (char*)NULL ;
03212         tmpcubes[i]->n_comments = 0 ;
03213         tmpcubes[i]->orig_ptype = BPP_DEFAULT ;
03214         tmpcubes[i]->filename = NULL ;
03215     }
03216     */
03217     tmpcubes[0]=cpl_imagelist_duplicate(cubes[0]);
03218    
03219     /*-------------------------------------------------------------------- 
03220      * center the cubes within the allocated big cube
03221      * that means define the (0,0) positions of the cubes in the image planes
03222      * to sub-pixel accuracy by using cumoffsetx,y and the reference cube
03223      */
03224     /* position of first reference frame, centered in big cube */
03225     llx0 = mlx/2 - clx/2 ;
03226     lly0 = mly/2 - cly/2 ;
03227 
03228     /*--------------------------------------------------------------------
03229      * go through the frame list and determine the lower left edge position
03230      * of the shifted cubes. Additionnally, the sub-pixel offsets are 
03231      * determined.
03232      */
03233 
03234 
03235     llx=cpl_calloc(n_cubes,sizeof(int)) ;
03236     lly=cpl_calloc(n_cubes,sizeof(int)) ;
03237     sub_offsetx=cpl_calloc(n_cubes,sizeof(float)) ;
03238     sub_offsety=cpl_calloc(n_cubes,sizeof(float)) ;
03239 
03240     for ( i = 0 ; i < n_cubes ; i++ )
03241     {
03242         llx[i] = llx0 - sinfo_new_nint(cumoffsetx[i]) ;
03243         sub_offsetx[i] = (float)sinfo_new_nint(cumoffsetx[i]) - cumoffsetx[i] ;
03244         lly[i] = lly0 - sinfo_new_nint(cumoffsety[i]) ;      
03245         sub_offsety[i] = (float)sinfo_new_nint(cumoffsety[i]) - cumoffsety[i] ;
03246     }
03247 
03248     /* -------------------------------------------------------------
03249      * shift the cubes according to the computed sub-pixel offsets 
03250      * that means shift the single image planes of each cube
03251      * first determine an interpolation kernel
03252      */
03253     if ( NULL == (kernel = sinfo_generate_interpolation_kernel(kernel_type)) )
03254     {
03255         sinfo_msg_warning ("could not generate desired interpolation kernel"
03256                            " or no kernel_typ was given, the default kernel"
03257                            " is used now!") ;
03258     }
03259     /* go through the frame list */
03260     for ( i = 0 ; i < n_cubes ; i++ )
03261     {
03262       /* go through the image planes and shift each plane by a 
03263          sub-pixel value */
03264       for ( z = 0 ; z < cnp ; z++ )
03265         {
03266       tmp_img=cpl_imagelist_get(cubes[i],z);
03267       if ( NULL == (shiftedImage = sinfo_new_shift_image(tmp_img,
03268                                sub_offsetx[i],
03269                                sub_offsety[i], 
03270                                kernel ) ) )
03271             {
03272           sinfo_msg_error ("could not shift image plane no %d "
03273                                "in cube no %d!", z, i) ;
03274           cpl_imagelist_delete(mergedCube) ;
03275           cpl_imagelist_delete(mask) ;
03276           cpl_free(kernel) ;
03277           return NULL ;
03278             }
03279       cpl_imagelist_set(tmpcubes[i],shiftedImage,z);
03280         }
03281     }
03282   
03283     cubedata=cpl_calloc(n_cubes,sizeof(float)) ;
03284 
03285     for ( y = 0 ; y < mly ; y++ )
03286     {
03287         for ( x = 0 ; x < mlx ; x++ )
03288         {
03289             for ( z = 0 ; z < mnp ; z++ )
03290             {
03291                 sum = 0. ;
03292                 sum2 = 0. ;
03293                 n = 0 ;
03294                 for ( i = 0 ; i < n_cubes ; i++ )
03295                 {
03296           c_img=cpl_imagelist_get(cubes[i],z);
03297 
03298           clx=cpl_image_get_size_x(c_img);
03299           cly=cpl_image_get_size_y(c_img);
03300 
03301           t_img=cpl_imagelist_get(tmpcubes[i],z);
03302                   ptdata=cpl_image_get_data_float(t_img);
03303 
03304           m_img=cpl_imagelist_get(mergedCube,z);
03305                   pmdata=cpl_image_get_data_float(m_img);
03306           o_img=cpl_imagelist_get(mask,z);
03307                   podata=cpl_image_get_data_float(o_img);
03308                   /* 
03309                     find the position of the present cube and go 
03310                     through the single spectra 
03311                    */
03312                     if ( y >= lly[i] && y < lly[i]+cly &&
03313                          x >= llx[i] && x < llx[i]+clx )
03314                     {
03315                         posx = x - llx[i] ;
03316                         posy = y - lly[i] ;
03317                         if (!isnan(ptdata[posx+posy*clx]))
03318                         {
03319                             sum += ptdata[posx+posy*clx] ;
03320                             sum2 += (ptdata[posx+posy*clx] *
03321                                      ptdata[posx+posy*clx]) ;
03322                             cubedata[n] = ptdata[posx+posy*clx] ;
03323                             n++ ;
03324                         }
03325                     }
03326                 }
03327                 
03328                 if ( n == 0 ) 
03329                 {
03330                     mean = 0. ;
03331                     sigma = 0. ;
03332                     pmdata[x+y*mlx] = 0. ;
03333                     podata[x+y*mlx] = 0 ;
03334                 }
03335                 else if ( n == 1 ) 
03336                 {
03337                     mean = sum ;
03338                     sigma = 0. ;
03339                     pmdata[x+y*mlx] = mean ;
03340                     podata[x+y*mlx] = 1 ;
03341                 }
03342                 else
03343                 {
03344                     mean = sum/(double)n ;
03345                     sigma = sqrt( (sum2 - sum*mean) / (double)(n - 1) ) ;
03346                     ns = 0 ; 
03347                     for ( i = 0 ; i < n ; i++ )
03348                     {
03349                         if ( cubedata[i] > mean+factor*sigma || 
03350                              cubedata[i] < mean-factor*sigma )
03351                         {
03352                             continue ;
03353                         }
03354                         else
03355                         {
03356                             pmdata[x+y*mlx] += cubedata[i] ;
03357                             ns++ ;
03358                         }
03359                     }
03360                     if ( ns == 0 )
03361                     {
03362                         pmdata[x+y*mlx] = 0. ;
03363                     }
03364                     else
03365                     {
03366                         pmdata[x+y*mlx] /= (float)ns ;
03367                     }
03368                     podata[x+y*mlx] = (float)ns ;
03369                 }
03370             }
03371         }
03372     }
03373 
03374     for( i = 0 ; i < n_cubes ; i++ )
03375     {
03376         cpl_imagelist_delete (tmpcubes[i]) ;
03377     }
03378     cpl_free(tmpcubes);
03379     cpl_free(llx);
03380     cpl_free(lly);
03381     cpl_free(sub_offsetx);
03382     cpl_free(sub_offsety);
03383     cpl_free(cubedata);
03384 
03385     /* convert the "free space" in the cube to blank pixels */
03386     sinfo_new_convert_0_to_ZERO_for_cubes(mergedCube) ; 
03387     cpl_free(kernel) ;
03388     return mask ;
03389 }
03390 
03391 cpl_imagelist * 
03392 sinfo_new_bin_cube(cpl_imagelist *cu, 
03393                              int xscale, 
03394                              int yscale, 
03395                              int xmin, 
03396                              int xmax, 
03397                              int ymin, 
03398                              int ymax)
03399 {
03400   int i,j,k;
03401   cpl_imagelist * cube;
03402   int ilx=0;
03403   int ily=0;
03404   int olx=0;
03405   int oly=0;
03406   int inp=0;
03407 
03408   float* pidata=NULL;
03409   float* podata=NULL;
03410   cpl_image* i_img=NULL;
03411   cpl_image* o_img=NULL;
03412 
03413 
03414   /* old code
03415   if (NULL == (cube = sinfo_newCube (xmax-xmin+1,ymax-ymin+1, cu->np)) ) 
03416   {
03417       sinfo_msg_error ("cannot allocate new cube") ;
03418       return NULL ;
03419   }
03420   */
03421   inp=cpl_imagelist_get_size(cu);
03422   i_img=cpl_imagelist_get(cu,0);
03423   ilx=cpl_image_get_size_x(i_img);
03424   ily=cpl_image_get_size_y(i_img);
03425   olx=xmax-xmin+1;
03426   oly=ymax-ymin+1;
03427 
03428     
03429   cube=cpl_imagelist_new();
03430   for ( i = 0 ; i < inp ; i++ ) {
03431     o_img = cpl_image_new(olx,oly,CPL_TYPE_FLOAT);
03432     cpl_imagelist_set(cube,o_img,i);
03433   }
03434 
03435 
03436   for (i=0;i<inp;i++){
03437       i_img=cpl_imagelist_get(cu,i);
03438       pidata=cpl_image_get_data_float(i_img);
03439       o_img=cpl_imagelist_get(cube,i);
03440       podata=cpl_image_get_data_float(o_img);
03441       for (j=0 ; j < olx ; j++) {
03442           for (k=0 ; k< oly ; k++) {
03443           podata[j+k*olx]=pidata[((int) (j+xmin)/xscale)+
03444                                      ((int) (k+ymin)/yscale)*ilx]/
03445                                        (xscale*yscale);
03446       }
03447       }
03448   }
03449 
03450   return cube;
03451 }
03452 
03453 
03454 cpl_imagelist * 
03455 sinfo_new_scale_cube(cpl_imagelist *cu, 
03456                                float xscale, 
03457                                float yscale, 
03458                                char * kernel_type)
03459 {
03460     cpl_imagelist   *   cube ;
03461     int             i, j, k, l ;
03462     int             lx_out, ly_out ;
03463     double          cur ;
03464     double      *   invert_transform ;
03465     double          neighbors[16] ;
03466     double          rsc[8],
03467                     sumrs ;
03468     double      param[6];
03469     double          x, y ;
03470     int             px, py ;
03471     int             pos ;
03472     int             tabx, taby ;
03473     double      *   kernel ;
03474     int             leaps[16] ;
03475     int                 ilx=0;
03476     int                 ily=0;
03477     int                 tlx=0;
03478     int                 tly=0;
03479     int                 inp;
03480     float*              podata=0;
03481     cpl_image*          in_img=NULL;
03482     cpl_image*          ou_img=NULL;
03483 
03484 
03485     if (cu == NULL)
03486     {
03487         sinfo_msg_error ("null cube") ;
03488         return NULL ;
03489     }
03490     
03491     param[0]=xscale;
03492     param[1]=0;
03493     param[2]=0;
03494     param[3]=0;
03495     param[4]=yscale;
03496     param[5]=0;
03497     
03498 
03499     invert_transform = sinfo_invert_linear_transform(param) ;
03500     if (invert_transform == NULL) {
03501         sinfo_msg_error("cannot compute sinfo_invert transform: "
03502                         "aborting warping") ;
03503         return NULL ;
03504     }
03505 
03506     /* Generate default interpolation kernel */
03507     kernel = sinfo_generate_interpolation_kernel(kernel_type) ;
03508     if (kernel == NULL) {
03509         sinfo_msg_error("cannot generate kernel: aborting resampling") ;
03510         return NULL ;
03511     }
03512 
03513     /* Compute new image size   */
03514     /* Compute new image size   */
03515     ilx=cpl_image_get_size_x(cpl_imagelist_get(cu,0));
03516     ily=cpl_image_get_size_y(cpl_imagelist_get(cu,0));
03517     inp=cpl_imagelist_get_size(cu);
03518 
03519     lx_out = (int) ilx*xscale ;
03520     ly_out = (int) ily*yscale ;
03521     
03522     cube=cpl_imagelist_new();
03523     for ( l = 0 ; l < inp ; i++ ) {
03524      in_img = cpl_image_new(ilx,ily,CPL_TYPE_FLOAT);
03525      cpl_imagelist_set(cube,in_img,l);
03526     }
03527 
03528     /* old code
03529     if (NULL == (cube = sinfo_newCube (lx_out, ly_out, cu->np)) ) 
03530     {
03531         sinfo_msg_error (" cannot allocate new cube") ;
03532         return NULL ;
03533     }
03534     */
03535 
03536     for (l=0;l<inp;l++){
03537       in_img=cpl_imagelist_get(cu,l);
03538       ou_img=cpl_imagelist_get(cube,l);
03539       tlx=cpl_image_get_size_x(in_img);
03540       tly=cpl_image_get_size_y(in_img);
03541       podata=cpl_image_get_data_float(ou_img);
03542         /* Pre compute leaps for 16 closest neighbors positions */
03543         leaps[0] = -1 - tlx ;
03544         leaps[1] =    - tlx ;
03545         leaps[2] =  1 - tlx ;
03546         leaps[3] =  2 - tlx ;
03547 
03548         leaps[4] = -1 ;
03549         leaps[5] =  0 ;
03550         leaps[6] =  1 ;
03551         leaps[7] =  2 ;
03552 
03553         leaps[8] = -1 + tlx ;
03554         leaps[9] =      tlx ;
03555         leaps[10]=  1 + tlx ;
03556         leaps[11]=  2 + tlx ;
03557 
03558         leaps[12]= -1 + 2*tlx ;
03559         leaps[13]=      2*tlx ;
03560         leaps[14]=  1 + 2*tlx ;
03561         leaps[15]=  2 + 2*tlx ;
03562 
03563         /* Double loop on the output image  */
03564         for (j=0 ; j < ly_out ; j++) {
03565             for (i=0 ; i< lx_out ; i++) {
03566                 /* Compute the original source for this pixel   */
03567 
03568                 x = invert_transform[0] * (double)i +
03569                     invert_transform[1] * (double)j +
03570                 invert_transform[2] ; 
03571 
03572                 y = invert_transform[3] * (double)i +
03573                 invert_transform[4] * (double)j +
03574                 invert_transform[5] ; 
03575 
03576             /* Which is the closest integer positioned neighbor?    */
03577                 px = (int)x ;
03578         py = (int)y ;
03579 
03580                 if ((px < 1) ||
03581                     (px > (tlx-2)) ||
03582                     (py < 1) ||
03583                     (py > (tly-2)))
03584                     podata[i+j*lx_out] = (pixelvalue)0.0 ;
03585                 else {
03586                     /* Now feed the positions for the closest 16 neighbors  */
03587                     pos = px + py * tlx ;
03588                     for (k=0 ; k<16 ; k++){
03589                         if(!isnan(podata[(int)(pos+leaps[k])])) neighbors[k] = 
03590                          (double)(podata[(int)(pos+leaps[k])]) ;
03591                 else neighbors[k]=0;
03592             }
03593 
03594                     /* Which tabulated value index shall we use?    */
03595                     tabx = (x - (double)px) * (double)(TABSPERPIX) ; 
03596                     taby = (y - (double)py) * (double)(TABSPERPIX) ; 
03597 
03598                     /* Compute resampling coefficients  */
03599                     /* rsc[0..3] in x, rsc[4..7] in y   */
03600 
03601                     rsc[0] = kernel[TABSPERPIX + tabx] ;
03602                     rsc[1] = kernel[tabx] ;
03603                     rsc[2] = kernel[TABSPERPIX - tabx] ;
03604                     rsc[3] = kernel[2 * TABSPERPIX - tabx] ;
03605                     rsc[4] = kernel[TABSPERPIX + taby] ;
03606                     rsc[5] = kernel[taby] ;
03607                     rsc[6] = kernel[TABSPERPIX - taby] ;
03608                     rsc[7] = kernel[2 * TABSPERPIX - taby] ;
03609 
03610                     sumrs = (rsc[0]+rsc[1]+rsc[2]+rsc[3]) *
03611                         (rsc[4]+rsc[5]+rsc[6]+rsc[7]) ;
03612 
03613                     /* Compute interpolated pixel now   */
03614                     cur =   rsc[4] * (  rsc[0]*neighbors[0] +
03615                                     rsc[1]*neighbors[1] +
03616                                     rsc[2]*neighbors[2] +
03617                                     rsc[3]*neighbors[3] ) +
03618                         rsc[5] * (  rsc[0]*neighbors[4] +
03619                                     rsc[1]*neighbors[5] +
03620                                     rsc[2]*neighbors[6] +
03621                                     rsc[3]*neighbors[7] ) +
03622                         rsc[6] * (  rsc[0]*neighbors[8] +
03623                                     rsc[1]*neighbors[9] +
03624                                     rsc[2]*neighbors[10] +
03625                                     rsc[3]*neighbors[11] ) +
03626                         rsc[7] * (  rsc[0]*neighbors[12] +
03627                                     rsc[1]*neighbors[13] +
03628                                     rsc[2]*neighbors[14] +
03629                                     rsc[3]*neighbors[15] ) ; 
03630 
03631                     /* Affect the value to the output image */
03632                     podata[i+j*lx_out] = (pixelvalue)(cur/sumrs) ;
03633                     /* done ! */
03634                 }       
03635             }
03636         }
03637     }
03638     cpl_free(kernel) ;
03639     cpl_free(invert_transform) ;
03640     return cube ;
03641 }
03642 
03643 
03653 cpl_imagelist * 
03654 sinfo_cube_zshift(const cpl_imagelist * cube_inp, 
03655                   const double shift,
03656                   double* sub_shift)
03657 {
03658 
03659     cpl_imagelist * cube_out=NULL ;
03660     cpl_image* img_inp=NULL;
03661     cpl_image* img_out=NULL;
03662     int        col, row,z ;
03663     int        int_shift ;
03664     int ilx=0;
03665     int ily=0;
03666     int ilz=0;
03667 
03668     int olx=0;
03669     int oly=0;
03670     int olz=0;
03671     int i=0;
03672     float* pidata=NULL;
03673     float* podata=NULL;
03674 
03675     cknull(cube_inp,"no input cube given!") ;
03676     check_nomsg(img_inp=cpl_imagelist_get(cube_inp,0));
03677     check_nomsg(ilx=cpl_image_get_size_x(img_inp));
03678     check_nomsg(ily=cpl_image_get_size_y(img_inp));
03679     check_nomsg(ilz=cpl_imagelist_get_size(cube_inp));
03680 
03681     olx=ilx;
03682     oly=ily;
03683     olz=ilz;
03684 
03685     int_shift = sinfo_new_nint(shift) ;
03686     *sub_shift = shift - (double) int_shift ;
03687     if ( int_shift == 0 )
03688     {
03689         cube_out =cpl_imagelist_duplicate(cube_inp) ;
03690         return cube_out ;
03691     }
03692     else
03693     {
03694       /* allocate memory */
03695       cknull(cube_out = cpl_imagelist_new(),"could not allocate memory!") ;
03696       for ( i = 0 ; i < olz ; i++ ) {
03697         check_nomsg(img_out=cpl_image_new(olx,oly,CPL_TYPE_FLOAT));
03698         check_nomsg(cpl_imagelist_set(cube_out,img_out,i));
03699      }
03700     }
03701 
03702     for(z=0; z< ilz; z++) {
03703       if ( (z-int_shift >= 0 ) && (z - int_shift < olz) ) {
03704         check_nomsg(img_inp=cpl_imagelist_get(cube_inp,z));
03705         check_nomsg(img_out=cpl_imagelist_get(cube_out,z-int_shift));
03706     check_nomsg(pidata=cpl_image_get_data_float(img_inp));
03707     check_nomsg(podata=cpl_image_get_data_float(img_out));
03708     for ( col = 0 ; col < ilx ; col++ ) {
03709       for ( row = 0 ; row < ily ; row++ ) {
03710         podata[col+row*olx] = pidata[col+row*olx] ;
03711       }
03712     }
03713       }
03714     }
03715     return cube_out ;
03716 
03717  cleanup:
03718     sinfo_free_imagelist(&cube_out);
03719     return NULL ;
03720 }
03721 
03731 cpl_imagelist * 
03732 sinfo_cube_zshift_poly(const cpl_imagelist * cube_inp, 
03733                        const double sub_shift,
03734                        const int    order)
03735 {
03736   cpl_imagelist * cube_out ;
03737 
03738   float* spec=NULL ;
03739   float* corrected_spec=NULL ;
03740   float* xnum=NULL ;
03741 
03742   float sum=0;
03743   float new_sum=0 ;
03744   float eval=0 ;
03745   float * imageptr=NULL ;
03746   int row=0;
03747   int col=0 ;
03748   int firstpos=0 ;
03749   int n_points=0 ;
03750   int i=0 ;
03751   int flag=0;
03752   int ilx=0;
03753   int ily=0;
03754   int ilz=0;
03755 
03756   int olx=0;
03757   int oly=0;
03758   int olz=0;
03759   int z=0;
03760 
03761   float* pidata=NULL;
03762   float* podata=NULL;
03763   cpl_image* img_inp=NULL;
03764   cpl_image* img_out=NULL;
03765 
03766   if ( cube_inp == NULL ) {
03767     sinfo_msg_error("no imagelist given!") ;
03768     return NULL ;
03769   }
03770 
03771   img_inp=cpl_imagelist_get(cube_inp,0);
03772 
03773   ilx=cpl_image_get_size_x(img_inp);
03774   ily=cpl_image_get_size_y(img_inp);
03775   ilz=cpl_imagelist_get_size(cube_inp);
03776 
03777   if ( order <= 0 ) {
03778     sinfo_msg_error("wrong order of interpolation polynom given!") ;
03779     return NULL ;
03780   }
03781 
03782 
03783   olx=ilx;
03784   oly=ily;
03785   olz=ilz;
03786   /* allocate memory */
03787 
03788   if ( NULL == (cube_out = cpl_imagelist_new()) ) {
03789     sinfo_msg_error ("could not allocate memory!") ;
03790     return NULL ;
03791   } else {
03792     for ( i = 0 ; i < ilz ; i++ ) {
03793       img_out=cpl_image_new(olx,oly,CPL_TYPE_FLOAT);
03794       cpl_imagelist_set(cube_out,img_out,i);
03795     }
03796   }
03797   
03798 
03799   n_points = order + 1 ;
03800   if ( n_points % 2 == 0 ) {
03801     firstpos = (int)(n_points/2) - 1 ;
03802   } else {
03803     firstpos = (int)(n_points/2) ;
03804   }
03805 
03806   spec=cpl_calloc(ilz,sizeof(float)) ;
03807   corrected_spec=cpl_calloc(ilz,sizeof(float)) ;
03808   xnum=cpl_calloc(order+1,sizeof(float)) ;
03809 
03810 
03811   /* fill the xa[] array for the polint function */
03812   for ( i = 0 ; i < n_points ; i++ ) {
03813     xnum[i] = i ;
03814   }
03815 
03816   for ( col = 0 ; col < ilx ; col++ ) {
03817     for ( row = 0 ; row < ily ; row++ ) {
03818       for( z=0; z< ilz; z++) {
03819         corrected_spec[z] = 0. ;
03820       }
03821       sum = 0. ;
03822       for ( z = 0 ; z < ilz ; z++ ) {
03823     img_inp=cpl_imagelist_get(cube_inp,z);
03824         pidata=cpl_image_get_data_float(img_inp);
03825     spec[z] = pidata[col + row*ilx] ;
03826     if (isnan(spec[z]) ) {
03827           spec[z] = 0. ;
03828                   
03829       for ( i = z - firstpos ; i < z-firstpos+n_points ; i++ ) {
03830         if ( i < 0 ) continue ;
03831             if ( i >= ilz) continue  ;
03832             corrected_spec[i] = ZERO ;
03833           }
03834         }
03835         if ( z != 0 && z != ilz - 1 ) {
03836           sum += spec[z] ;
03837         }
03838 
03839       }
03840        
03841       new_sum = 0. ;
03842       for ( z = 0 ; z < ilz ; z++ ) {
03843    
03844         /* ---------------------------------------------------------------
03845          * now determine the arrays of size n_points with which the
03846          * polynom is determined and determine the position eval
03847          * where the polynom is evaluated in polynomial interpolation.
03848          * Take care of the points near the row edges!
03849          */
03850         if (isnan(corrected_spec[z])) continue ;
03851         if ( z - firstpos < 0 ) {
03852           imageptr = &spec[0] ;
03853           eval     = sub_shift + z ;
03854         } else if ( z - firstpos + n_points >= ilz ) {
03855           imageptr = &spec[ilz - n_points] ;
03856           eval     = sub_shift + z + n_points - ilz ;
03857         } else {
03858       imageptr = &spec[z-firstpos] ;
03859           eval     = sub_shift + firstpos ;
03860         }
03861 
03862         flag=0;
03863         corrected_spec[z]=sinfo_new_nev_ille(xnum,imageptr,order,eval,&flag);
03864         if ( z != 0 && z != ilz - 1 ) {
03865           new_sum += corrected_spec[z] ;
03866         }
03867       }
03868 
03869       /* fill the output spectrum */
03870       for (z = 0 ; z < ilz ; z++ )
03871       {
03872         img_out=cpl_imagelist_get(cube_out,z);
03873         podata=cpl_image_get_data_float(img_out);
03874         if ( new_sum == 0. ) {
03875           new_sum = 1. ;
03876         }
03877         if ( z == 0 ) {
03878           podata[col+row*olx] = ZERO ;
03879         } else if ( z == ilz - 1 ) {
03880           podata[col+row*olx] = ZERO ;
03881         } else if ( isnan(corrected_spec[z]) ) {
03882           podata[col+row*olx] = ZERO ;
03883         } else {
03884           corrected_spec[z] *= sum / new_sum ;
03885           podata[col+row*olx] = corrected_spec[z] ;
03886     }
03887       }
03888 
03889     }
03890   }
03891 
03892   cpl_free(spec) ;
03893   cpl_free(corrected_spec) ;
03894   cpl_free(xnum) ;
03895   return cube_out ;
03896 
03897 
03898 }
03899 
03908 cpl_imagelist * 
03909 sinfo_cube_zshift_spline3(const cpl_imagelist * cube_inp, 
03910                           const double sub_shift)
03911 {
03912 
03913   cpl_imagelist * cube_out=NULL ;
03914   float* spec=NULL ;
03915   float* corrected_spec=NULL ;
03916   float* xnum=NULL ;
03917   float* eval=NULL ;
03918   float sum=0;
03919   float new_sum=0 ;
03920   int row=0;
03921   int col=0;
03922   int i=0;
03923   int z=0;
03924 
03925   int ilx=0;
03926   int ily=0;
03927   int ilz=0;
03928   int olx=0;
03929   int oly=0;
03930   int olz=0;
03931 
03932   float* pidata=NULL;
03933   float* podata=NULL;
03934   cpl_image* img_inp=NULL;
03935   cpl_image* img_out=NULL;
03936 
03937   if ( cube_inp == NULL ) {
03938     sinfo_msg_error("no imagelist given!") ;
03939     return NULL ;
03940   }
03941 
03942   img_inp=cpl_imagelist_get(cube_inp,0);
03943   ilx=cpl_image_get_size_x(img_inp);
03944   ily=cpl_image_get_size_y(img_inp);
03945   ilz=cpl_imagelist_get_size(cube_inp);
03946 
03947 
03948   olx=ilx;
03949   oly=ily;
03950   olz=ilz;
03951   /* allocate memory */
03952   if ( NULL == (cube_out = cpl_imagelist_new()) ) {
03953     sinfo_msg_error ("could not allocate memory!") ;
03954     return NULL ;
03955   } else {
03956     for ( i = 0 ; i < ilz ; i++ ) {
03957       img_out=cpl_image_new(olx,oly,CPL_TYPE_FLOAT);
03958       cpl_imagelist_set(cube_out,img_out,i);
03959     }
03960   }
03961 
03962   xnum=cpl_calloc(ilz,sizeof(float)) ;
03963   /* fill the xa[] array for the spline function */
03964   for ( i = 0 ; i < ilz ; i++ ) {
03965     xnum[i] = i ;
03966   }
03967 
03968   spec=cpl_calloc(ilz,sizeof(float)) ;
03969   corrected_spec=cpl_calloc(ilz,sizeof(float)) ;
03970   eval=cpl_calloc(ilz,sizeof(float)) ;
03971 
03972   for ( col = 0 ; col < ilx ; col++ ) {
03973     for ( row = 0 ; row < ily ; row++ ) {
03974       sum = 0. ;
03975       for ( z = 0 ; z < ilz ; z++ ) {
03976     img_inp=cpl_imagelist_get(cube_inp,z);
03977         pidata=cpl_image_get_data_float(img_inp);
03978     spec[z] = pidata[col + row*ilx] ;
03979     if (isnan(spec[z]) ) {
03980       for ( i = z-1 ; i <= z+1 ; i++ ) {
03981         if ( i < 0 ) continue ;
03982         if ( i >= ilz) continue ;
03983         corrected_spec[i] = ZERO ;
03984       }
03985       spec[z] = 0. ;
03986     }
03987     sum += spec[z] ;
03988     eval[z] = (float)sub_shift+(float)z ;
03989       }
03990       /* now we do the spline interpolation*/
03991       if ( -1 == sinfo_function1d_natural_spline( xnum, spec, ilz, eval, 
03992                                               corrected_spec, ilz ) )
03993         {
03994       sinfo_msg_error("error in spline interpolation!") ;
03995       return NULL ;
03996         }
03997         
03998       new_sum = 0. ;
03999       for ( z = 0 ; z < ilz ; z++ ) {
04000     if ( isnan(corrected_spec[z]) ) {
04001       continue ;
04002     }
04003     new_sum += corrected_spec[z] ;
04004       }
04005       /* fill output imagelist */
04006       for ( z = 0 ; z < ilz ; z++ ) {
04007         img_out=cpl_imagelist_get(cube_out,z);
04008         podata=cpl_image_get_data_float(img_out);
04009     if ( new_sum == 0. ) new_sum =1. ;
04010     {
04011       if ( isnan(corrected_spec[z]) ) {
04012         podata[col+row*olx] = ZERO ;
04013       } else {
04014         corrected_spec[z] *= sum / new_sum ;
04015         podata[col+row*olx] = corrected_spec[z] ;
04016       }
04017     }
04018       }
04019     }
04020   }
04021   cpl_free(xnum);
04022   cpl_free(spec) ;
04023   cpl_free(corrected_spec) ;
04024   cpl_free(eval) ;
04025 
04026   return cube_out ;
04027 }

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