cube_ops.c

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

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