sinfo_cube_construct.c

00001 /*
00002  * This file is part of the ESO SINFONI Pipeline
00003  * Copyright (C) 2004,2005 European Southern Observatory
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or
00008  * (at your option) any later version.
00009  *
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00018  */
00019 /************************************************************************
00020 * E.S.O. - VLT project
00021 *
00022 *
00023 *
00024 * who       when      what
00025 * --------  --------  ----------------------------------------------
00026 * schreib  30/08/00  created
00027 */
00028 
00029 /************************************************************************
00030 *   NAME
00031 *        sinfo_cube_construct.c -
00032 *        some procedures to construct a data cube
00033 *
00034 *   SYNOPSIS
00035 *
00036 *   1) cpl_image * sinfo_new_convolve_ns_image_by_gauss( cpl_image * lineImage,
00037 *                                         int        hw )
00038 *
00039 *   2) float * sinfo_north_south_test( cpl_image * ns_image,
00040 *                                int        n_slitlets,
00041 *                                int        halfWidth,
00042 *                                float      fwhm,
00043 *                                float      minDiff,
00044 *                                float      estimated_dist,
00045 *                                float      devtol )
00046 *
00047 *   3) cpl_imagelist * sinfo_new_make_cube ( cpl_image * calibImage,
00048 *                           float    * distances,
00049 *                           float    * correct_diff_dist )
00050 *
00051 *   4) cpl_imagelist * sinfo_new_make_cube_spi ( cpl_image *  calibImage,
00052 *                              float    ** slit_edges,
00053 *                              float    *  shift )
00054 *
00055 *   5) cpl_imagelist * sinfo_new_make_cube_dist ( cpl_image * calibImage,
00056 *                               float      firstCol,
00057 *                               float    * distances,
00058 *                               float    * shift )
00059 *
00060 *   6) cpl_imagelist * sinfo_new_make_3D_cube_dist ( cpl_image * calibImage,
00061 *                                 float      firstCol,
00062 *                                 float    * distances,
00063 *                                 float    * shift )
00064 *
00065 *   7) cpl_imagelist * sinfo_new_make_3D_cube ( cpl_image * calibImage,
00066 *                             int      * kpixshift,  
00067 *                             int        kpixfirst )
00068 *
00069 *   8) cpl_imagelist * 
00070        sinfo_new_determine_mask_cube(cpl_imagelist * sourceMaskCube,
00071 *                                    float     lowLimit,
00072 *                                    float     highLimit )
00073 *
00074 *   9) cpl_imagelist * sinfo_new_interpol_cube ( cpl_imagelist * sourceCube,
00075 *                               cpl_imagelist * maskCube,
00076 *                               int       n_neighbors, 
00077 *                               int       max_radius )
00078 *
00079 *   10) cpl_imagelist * sinfo_new_fine_tune_cube( cpl_imagelist * cube,
00080 *                              float   * correct_diff_dist )
00081 *
00082 *   11) cpl_imagelist * sinfo_new_fine_tune_cube_by_FFT( cpl_imagelist * cube,
00083 *                                   float   * correct_diff_dist )
00084 *
00085 *   12) cpl_imagelist * sinfo_new_fine_tune_cube_by_spline(cpl_imagelist * cube,
00086 *                                        float   * correct_diff_dist )
00087 *
00088 *   DESCRIPTION
00089 *
00090 *   1) convolves a north-south-test image with a sinfo_gaussian
00091 *      with user given integer half width by using the eclipse
00092 *      routine sinfo_function1d_filter_lowpass().
00093 *   2) determines the distances of the slitlets
00094 *   3) makes a data cube out of a resampled source image
00095 *      this SPIFFI specific routine takes into account the
00096 *      Spiffi slitlet order on the detector.
00097 *      Also shifts the resulting image rows by one pixel if
00098 *      necessary according to the distances array gained from
00099 *      the north-south test routine.
00100 *      Can do the same with the bad pixel map image to generate a
00101 *      bad pixel mask cube.
00102 *   4) makes a data cube out of a resampled source image
00103 *      this SPIFFI specific routine takes into account the
00104 *      Spiffi slitlet order on the detector.
00105 *      This routine takes fitted slitlet positions into account.
00106 *      Can do the same with the bad pixel map image to generate a
00107 *      bad pixel mask cube.
00108 *   5) makes a data cube out of a resampled source image
00109 *      this SPIFFI specific routine takes into account the
00110 *      Spiffi slitlet order on the detector.
00111 *      Also shifts the resulting image rows by one pixel if
00112 *      necessary according to the distances array gained from
00113 *      the north-south test routine.
00114 *      Can do the same with the bad pixel map image to generate a
00115 *      bad pixel mask cube.
00116 *   6) makes a data cube out of a resampled source image
00117 *      this 3D specific routine takes into account the
00118 *      3D slitlet order on the detector.
00119 *      Also shifts the resulting image rows by one pixel if
00120 *      necessary according to the distances array gained from
00121 *      the north-south test routine.
00122 *      Can do the same with the bad pixel map image to generate a
00123 *      bad pixel mask cube.
00124 *   7) makes a data cube out of a resampled source image
00125 *      this MPE 3D specific routine takes into account the
00126 *      3D slitlet order on the detector.
00127 *      Also shifts the resulting image row by an integer pixel shift if
00128 *      necessary according to the input kpixshift array 
00129 *      Can do the same with the bad pixel map image to generate a
00130 *      bad pixel mask cube.
00131 *   8) converts resampled bad pixels to real bad pixels in data cubes.
00132 *   9) Bad pixel interpolation 3D like (saturated pixels exist):
00133 *      interpolates the bad pixels of the source cube by
00134 *      using the nearest neighbors. 
00135 *      first it is checked if the bad pixel is interpolatable:
00136 *      it is only interpolatable if the number of good pixels 
00137 *      in its spectrum of length 2*n_neighbors+1 exceeds 3 and
00138 *      if there is at least one good pixel on either side of the
00139 *      central pixel.
00140 *      Afterwards good neighboring pixels are searched within the 
00141 *      image plane of the bad pixel by using an increasing pixel radius. 
00142 *      Good pixels mean, the corresponding spectral pixels of the 
00143 *      bad pixel and its spatial neighboring pixel must have
00144 *      at least 2 valid pixel pairs to be able to be used for
00145 *      the interpolation. The search is stopped if 9 valid neighboring
00146 *      pixels are found. 
00147 *      Now normalize the found spectral values, collect the valid pixels 
00148 *      (there must be at least 18) and take the sinfo_median of the valid 
00149 *      pixels with which the bad pixel is replaced.
00150 *   10) fine tunes each row in the right position according 
00151 *      to the distances of the slitlets to each other
00152 *      (output of the north-south test).
00153 *      This means that the rows must be realigned by a 
00154 *      fraction of a pixel to accomodate non-integer slit 
00155 *      length. The fractional realignment is done by using
00156 *      the polynomial interpolation algorithm of N.R. 
00157 *      Each row is rescaled so that the total flux is
00158 *      conserved.
00159 *  11) fine tunes each row in the right position according 
00160 *      to the distances of the slitlets to each other
00161 *      (output of the north-south test).
00162 *      This means that the rows must be realigned by a 
00163 *      fraction of a pixel to accomodate non-integer slit 
00164 *      length. The fractional realignment is done by using
00165 *      the FFT algorithm four1() of N.R. 
00166 *  12) fine tunes each row in the right position according 
00167 *      to the distances of the slitlets to each other
00168 *      (output of the north-south test).
00169 *      This means that the rows must be realigned by a 
00170 *      fraction of a pixel to accomodate non-integer slit 
00171 *      length. The fractional realignment is done by using
00172 *      the spline interpolation algorithm splint in connection
00173 *      with the algorithm spline of N.R. 
00174 *      This algorithms assume that each row is a tabulated
00175 *      function. The first derivatives of the interpolating
00176 *      function at the first and last point must be given.
00177 *      These are set higher than 1xe^30, so the routine
00178 *      sets the corresponding boundary condition for a natural
00179 *      spline, with zero second derivative on that boundary.
00180 *      Each row is rescaled so that the total flux is
00181 *      conserved.
00182 *
00183 *   FILES
00184 *
00185 *   ENVIRONMENT
00186 *
00187 *   RETURN VALUES
00188 *
00189 *   CAUTIONS
00190 *
00191 *   EXAMPLES
00192 *
00193 *   SEE ALSO
00194 *
00195 *   BUGS
00196 *
00197 *------------------------------------------------------------------------
00198 */
00199 
00200 #ifdef HAVE_CONFIG_H
00201 #  include <config.h>
00202 #endif
00203 #define POSIX_SOURCE 1
00204 #include "sinfo_vltPort.h"
00205 
00206 /*
00207  * System Headers
00208  */
00209 
00210 /*
00211  * Local Headers
00212  */
00213 #include "sinfo_function_1d.h"
00214 #include "sinfo_cube_construct.h"
00215 #include "sinfo_spectrum_ops.h"
00216 #include "sinfo_wave_calibration.h"
00217 #include "sinfo_utilities.h"
00218 #include "sinfo_local_types.h"
00219 #include "sinfo_fft_base.h"
00228 /*----------------------------------------------------------------------------
00229  *                          Function codes
00230  *--------------------------------------------------------------------------*/
00241 cpl_image * 
00242 sinfo_new_convolve_ns_image_by_gauss( cpl_image * lineImage,
00243                                    int        hw )
00244 {
00245     cpl_image * returnImage ;
00246     float* row_buffer=NULL ;
00247     float * filter ;
00248     int col, row ;
00249     int ilx=0;
00250     int ily=0;
00251  
00252     float* pidata=NULL;
00253     float* podata=NULL;
00254 
00255     if ( lineImage == NULL )
00256     {
00257         sinfo_msg_error("no input image given!\n") ;
00258         return NULL ;
00259     }
00260     ilx=cpl_image_get_size_x(lineImage);
00261     ily=cpl_image_get_size_y(lineImage);
00262     pidata=cpl_image_get_data_float(lineImage);
00263     if ( hw < 1 )
00264     {
00265         sinfo_msg_error(" wrong half width given!\n") ;
00266         return NULL ;
00267     }
00268 
00269     /* allocate memory for returned image */
00270     if ( NULL == ( returnImage = cpl_image_new(ilx,ily,CPL_TYPE_FLOAT ) ))
00271     {
00272         sinfo_msg_error("cannot allocate a new image\n");
00273         return NULL ;
00274     }
00275     podata=cpl_image_get_data_float(returnImage);
00276 
00277     /* go through the image rows and save them in a buffer */
00278     row_buffer=cpl_calloc(ily,sizeof(float)) ;
00279 
00280     for ( row = 0 ; row < ily ; row++ )
00281     {
00282         for ( col = 0 ; col < ilx ; col++ )
00283         {
00284             if ( isnan(pidata[col+row*ilx]) )
00285             {
00286                 row_buffer[col] = 0. ;
00287             }
00288             else
00289             {
00290                 row_buffer[col] = pidata[col + row*ilx] ;
00291             }
00292         }
00293 
00294         /*--------------------------------------------------------------------
00295          * now low pass filter the rows by the gaussian and fill the return
00296          * image.
00297          */
00298         filter = sinfo_function1d_filter_lowpass( row_buffer,
00299                                             ilx,
00300                                             LOW_PASS_GAUSSIAN,
00301                                             hw ) ;
00302         for ( col = 0 ; col < ily ; col++ )
00303         {
00304             podata[col + row*ilx] = filter[col] ;
00305         }
00306         /* deallocate memory */
00307         sinfo_function1d_del (filter) ;
00308     }
00309     cpl_free(row_buffer);    
00310     return returnImage ;
00311 }
00312 
00329 float * 
00330 sinfo_north_south_test( cpl_image * ns_image,
00331                           int        n_slitlets,
00332                           int        halfWidth,
00333                           float      fwhm,
00334                           float      minDiff,
00335                           float      estimated_dist,
00336                           float      devtol,
00337               int        bottom,
00338               int        top )
00339 {
00340     int i, j, k, m, row, col, n, ni, na ;
00341     int position, counter, iters ;
00342     int xdim, ndat, its, numpar ;
00343     pixelvalue row_buf[cpl_image_get_size_x(ns_image)] ;
00344     float sum, mean, maxval ;
00345     float tol, lab ;
00346     float * distances ;
00347     float distances_buf[cpl_image_get_size_y(ns_image)][n_slitlets-1] ;
00348     float x_position[n_slitlets] ;
00349     float * xdat, * wdat ;
00350     int * mpar ;
00351     int found[3*n_slitlets], found_clean[3*n_slitlets] ;
00352     int found_cleanit[3*n_slitlets] ;
00353     Vector * line ;
00354     FitParams ** par ;
00355     int foundit, begin, end ;
00356     int zeroindicator ;
00357     int ilx=0;
00358     int ily=0;
00359  
00360     float* pidata=NULL;
00361  
00362     if ( ns_image == NULL )
00363     {
00364         sinfo_msg_error("sorry, no image given\n") ;
00365         return NULL ;
00366     }
00367     ilx=cpl_image_get_size_x(ns_image);
00368     ily=cpl_image_get_size_y(ns_image);
00369     pidata=cpl_image_get_data_float(ns_image);
00370 
00371 
00372     if ( n_slitlets < 1 )
00373     {
00374         sinfo_msg_error("wrong number of slitlets given\n") ;
00375         return NULL ;
00376     }
00377     if ( halfWidth < 0 || halfWidth >= estimated_dist )
00378     {
00379         sinfo_msg_error("wrong half width given\n") ;
00380         return NULL ;
00381     }
00382     if ( fwhm <= 0. )
00383     {
00384         sinfo_msg_error("wrong fwhm given\n") ;
00385         return NULL ;
00386     }
00387     if ( minDiff < 1. )
00388     {
00389         sinfo_msg_error("wrong minDiff given\n") ;
00390         return NULL ;
00391     }
00392 
00393     /* allocate memory for output array */
00394     if (NULL == (distances = (float *) cpl_calloc ( n_slitlets - 1 , 
00395                                                     sizeof (float) ))) 
00396     {
00397         sinfo_msg_error("could not allocate memory\n") ;
00398         return NULL ;
00399     }
00400 
00401     /* go through the image rows */
00402     for ( row = bottom ; row < top ; row++ )
00403     {
00404         zeroindicator = 0 ;
00405 
00406         /* initialize the distance buffer */
00407         for ( i = 0 ; i < n_slitlets-1 ; i++ )
00408         {
00409             distances_buf[row][i] = ZERO ;
00410         }
00411 
00412         /* fill the row buffer array with image data */
00413         for ( col = 0 ; col < ilx ; col++ )
00414         {
00415             row_buf[col] = pidata[col + row*ilx] ;
00416         }
00417 
00418         /* determine the mean of the row data */
00419         sum = 0. ;
00420         n = 0 ;
00421         for ( i = 0 ; i < ilx ; i++ )
00422         {
00423             if ( isnan(row_buf[i]) )
00424             {
00425                 continue ;
00426             }
00427             sum += row_buf[i] ;
00428             n++ ;
00429         }
00430         mean = sum / (float)n ;
00431 
00432 
00433         /* store the positions of image values greater than the mean */
00434         n = 0 ;
00435         for ( i = 0 ; i < ilx ; i++ )
00436         {
00437             if (isnan(row_buf[i]))
00438             {
00439                 continue ;
00440             }
00441             if ( row_buf[i] >  sqrt(mean*mean*9) )
00442             {
00443                 found[n] = i ;
00444                 n++ ;
00445             } 
00446         }
00447        
00448         if ( n < n_slitlets )
00449         {
00450             sinfo_msg_warning("t1 wrong number of intensity columns found "
00451                               "in row: %d, found number: %d, mean: %g",
00452                               row, n, mean) ;
00453             continue ;
00454         }
00455         else
00456         { 
00457             /* find the maximum value position around the found columns */
00458             na = 0 ;
00459             for ( i = 1 ; i < n ; i ++ )
00460             {
00461                 if ( found[i] - found[i-1] < halfWidth )
00462                 {
00463                     begin = found[i] - halfWidth ;
00464                     if ( begin < 0 )
00465                     {
00466                         begin = 0 ;
00467                     }
00468                     end = found[i] + halfWidth ;
00469                     if ( end >= ilx )
00470                     {
00471                         end = ilx - 1 ;
00472                     }
00473                     /* find the maximum value inside the box 
00474                        around the found positions*/
00475                     maxval = -FLT_MAX ;
00476                     foundit = 0 ;
00477                     for ( j = begin ; j <= end ; j++ )
00478                     {
00479                         /* do not consider boxes that contain bad pixels */
00480                         if (isnan(row_buf[j]))
00481                         {
00482                             continue ;
00483                         }
00484                         if (row_buf[j] >= maxval )
00485                         {
00486                             maxval = row_buf[j] ;
00487                             foundit = j ;
00488                         }
00489                     }
00490                     if (maxval == -FLT_MAX)
00491                     {
00492                         continue ;
00493                     }
00494                     for ( k = 0 ; k < na ; k++ )
00495                     {
00496                         if ( found_cleanit[k] >= begin && 
00497                              found_cleanit[k] < foundit )
00498                         {
00499                             na-- ;
00500                         }
00501                     }
00502                     for ( k = 0 ; k < n ; k++ )
00503                     {
00504                         if ( found[k] == foundit)
00505                         {
00506                  if (na>0){
00507                             if ( found_cleanit[na-1] != found[k] )
00508                             {
00509                                 found_cleanit[na] = found[k] ;
00510                                 na++ ;
00511                             }
00512              }
00513              else{
00514                 found_cleanit[na] = found[k] ;  
00515                             na++ ;
00516              } 
00517                       }
00518                     }
00519                 }
00520                 else
00521                 {
00522                     if ( i == 1 )
00523                     {
00524                         found_cleanit[na] = found[0] ;
00525                         na++ ;
00526                         found_cleanit[na] = found[1] ;
00527                         na++ ;
00528                     }
00529                     else
00530                     {   
00531                 if (na>0){
00532                             if ( found_cleanit[na-1] != found[i-1])
00533                             {
00534                                 found_cleanit[na] = found[i-1] ;
00535                                 na++ ;
00536                             }
00537                             if ( found_cleanit[na-1] != found[i])
00538                             {
00539                                 found_cleanit[na] = found[i] ;
00540                                 na++ ;
00541                             }
00542                         }
00543             else
00544                         {
00545                             found_cleanit[na] = found[i] ;
00546                             na++ ;
00547                         }
00548             }  
00549                 }
00550             }
00551             /* determine only one pixel position for each slitlet intensity */
00552             j = 1 ;
00553             for ( i = 1 ; i < na ; i++ )
00554             {
00555                 if ( (float)(found_cleanit[i] - found_cleanit[i-1]) < 
00556                             (estimated_dist - devtol) ||
00557                      (float)(found_cleanit[i] - found_cleanit[i-1]) > 
00558                             (estimated_dist + devtol) )
00559                 {
00560                     continue ;
00561                 }
00562                 else
00563                 {
00564                     found_clean[j-1] = found_cleanit[i-1] ;
00565                     found_clean[j]   = found_cleanit[i] ;
00566                     j++ ;
00567                 }
00568             }
00569         }
00570         if ( j > n_slitlets )
00571         {
00572             /* check the distance again */
00573             ni = 1 ;
00574             for ( i = 1 ; i < j ; i++ )
00575             {
00576                 if ( (float)(found_clean[i] - found_clean[i-1]) < 
00577                             (estimated_dist - devtol ) ||
00578                      (float)(found_clean[i] - found_clean[i-1]) > 
00579                             (estimated_dist + devtol ) )
00580                 { 
00581                     continue ;
00582                 }
00583                 else
00584                 {
00585 
00586                     found_clean[ni-1] = found_clean[i-1] ;
00587                     found_clean[ni]   = found_clean[i] ;
00588                     ni++ ;
00589                 }
00590             }
00591             if ( ni != n_slitlets )
00592             {
00593                 sinfo_msg_warning("t2 wrong number of intensity columns"
00594                                   " found in row: %d,  found number: %d",
00595                                   row, ni) ;
00596                 continue ;
00597             }
00598             else 
00599             {
00600                 j = ni ;
00601             }
00602         }
00603         else if ( j < n_slitlets )
00604         {
00605             cpl_msg_debug ("north_south_test3:",
00606                             "t3 wrong number of intensity columns "
00607                             "found in row: %d , found number: %d, mean: %g\n", 
00608                             row, j, mean) ;
00609             continue ;
00610         }
00611         counter = 0 ;
00612         /* go through the found intensity pixels in one row */
00613         for ( i = 0 ; i < j ; i++ )
00614         {
00615             /* allocate memory for the array where the line is fitted in */
00616             if ( NULL == (line = sinfo_new_vector (2*halfWidth + 1)) )
00617             {
00618                 sinfo_msg_error ("cannot allocate new Vector \n") ;
00619                 cpl_free(distances) ;
00620                 return NULL ;
00621             }
00622 
00623             /* allocate memory */
00624             xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00625             wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00626             mpar = (int *)   cpl_calloc( MAXPAR, sizeof (int) ) ;
00627             par = sinfo_new_fit_params(1) ;
00628 
00629             m = 0 ;
00630             for ( k = found_clean[i]-halfWidth ; 
00631                   k <= found_clean[i]+halfWidth ; k++ )
00632             {
00633                 if ( k < 0 )
00634                 {
00635                     k = 0. ;
00636                 }
00637                 else if ( k >= ilx )
00638                 {
00639                     k = ilx - 1 ;
00640                 }
00641                 else if ( isnan(row_buf[k]) )
00642                 {
00643                     zeroindicator = 1 ;
00644                     break ;
00645                 }
00646                 else
00647                 {
00648                     line -> data[m] = row_buf[k] ;
00649                     m++ ;
00650                 }
00651             }
00652             if ( zeroindicator == 1 )
00653             {
00654                 sinfo_new_destroy_vector(line) ;
00655                 cpl_free(xdat) ;
00656                 cpl_free(wdat) ;
00657                 cpl_free(mpar) ;
00658                 sinfo_new_destroy_fit_params(&par) ;
00659                 break ;
00660             }
00661 
00662             /*----------------------------------------------------------------
00663              * go through the spectral sinfo_vector
00664              * determine the maximum pixel value in the spectral sinfo_vector
00665              */
00666             maxval = -FLT_MAX ;
00667             position = -INT32_MAX ;
00668             for ( k = 0 ; k < m ; k++ )
00669             {
00670                 xdat[k] = k ;
00671                 wdat[k] = 1.0 ;
00672                 if ( line -> data[k] >= maxval )
00673                 {
00674                     maxval = line -> data[k] ;
00675                     position = k ;
00676                 }
00677             }
00678 
00679             /* set initial values for the fitting routine */
00680             xdim     = XDIM ;
00681             ndat     = line -> n_elements ;
00682             numpar   = MAXPAR ;
00683             tol      = TOL ;
00684             lab      = LAB ;
00685             its      = ITS ;
00686             (*par) -> fit_par[1] = fwhm ;
00687             (*par) -> fit_par[2] = (float) position ;
00688             (*par) -> fit_par[3] = (float) (line -> data[0] + 
00689                                    line -> data[line->n_elements - 1]) / 2.0 ;
00690             (*par) -> fit_par[0]  = maxval - ((*par) -> fit_par[3]) ;
00691 
00692 
00693             /* exclude negative peaks and low signal cases */
00694             if ( (*par) -> fit_par[0] < minDiff )
00695             {
00696                 sinfo_msg_warning ("sorry, signal of line too low to fit "
00697                                    "in row: %d in slitlet %d\n", row, i) ;
00698                 sinfo_new_destroy_vector(line) ;
00699                 cpl_free(xdat) ;
00700                 cpl_free(wdat) ;
00701                 cpl_free(mpar) ;
00702                 sinfo_new_destroy_fit_params(&par) ;
00703                 continue ;
00704             }
00705 
00706             for ( k = 0 ; k < MAXPAR ; k++ )
00707             {
00708                 (*par) -> derv_par[k] = 0.0 ;
00709                 mpar[k] = 1 ;
00710             }
00711             /* finally, do the least square fit using a Gaussian */
00712             if ( 0 > ( iters = sinfo_new_lsqfit_c( xdat, &xdim, 
00713                                                    line -> data, wdat, &ndat, 
00714                                                    (*par) -> fit_par,
00715                                                    (*par) -> derv_par, mpar, 
00716                                                    &numpar, &tol, &its, &lab)) )
00717             {
00718           /*
00719                 cpl_msg_debug ("north_south_test:",
00720                                "sinfo_lsqfit_c: least squares fit failed,"
00721                                " error no.: %d in row: %d in slitlet %d\n",
00722                                 iters, row, i) ;
00723           */
00724                 sinfo_new_destroy_vector(line) ;
00725                 cpl_free(xdat) ;
00726                 cpl_free(wdat) ;
00727                 cpl_free(mpar) ;
00728                 sinfo_new_destroy_fit_params(&par) ;
00729                 continue ;
00730             }
00731 
00732             /* check for negative fit results */
00733             if ( (*par) -> fit_par[0] <= 0. || 
00734                  (*par) -> fit_par[1] <= 0. ||
00735                  (*par) -> fit_par[2] < 0. )
00736             {
00737                 sinfo_msg_warning ("negative parameters as fit result, "
00738                                    "not used! in row %d in slitlet %d", 
00739                                    row, i) ;
00740                 sinfo_new_destroy_vector(line) ;
00741                 cpl_free(xdat) ;
00742                 cpl_free(wdat) ;
00743                 cpl_free(mpar) ;
00744                 sinfo_new_destroy_fit_params(&par) ;
00745                 continue ;
00746             }
00747 
00748             /* correct the fitted position for the given row of the line 
00749                in image coordinates */
00750             (*par) -> fit_par[2] =  (float) (found_clean[i] - halfWidth) + 
00751                                             (*par) -> fit_par[2] ;
00752             x_position[counter] = (*par) -> fit_par[2] ;
00753             counter ++ ;
00754 
00755             /* free memory */
00756             sinfo_new_destroy_fit_params(&par) ;
00757             sinfo_new_destroy_vector ( line ) ;
00758             cpl_free ( xdat ) ;
00759             cpl_free ( wdat ) ;
00760             cpl_free ( mpar ) ;
00761         }
00762         if (zeroindicator == 1)
00763         {
00764             sinfo_msg_debug ("bad pixel in fitting box in row: %d\n", row) ;
00765             continue ;
00766         }
00767 
00768         if ( counter != n_slitlets )
00769         {
00770             continue ;
00771             sinfo_msg_warning("wrong number of slitlets found in row: %d",row);
00772         }
00773         /* store the distances between the sources in a buffer */
00774         for ( i = 1 ; i < n_slitlets ; i++ )
00775         {
00776             distances_buf[row][i-1] = x_position[i] - x_position[i-1] ;
00777         }
00778     }
00779 
00780     /* ----------------------------------------------------------------
00781      * go through the rows again and take the mean of the distances, 
00782      * throw away the runaways 
00783      */
00784     for ( i = 0 ; i < n_slitlets-1 ; i++ )
00785     {
00786         n   = 0 ;
00787         sum = 0. ;
00788         for ( row = bottom ; row < top ; row++ )
00789         {
00790             if ( fabs( distances_buf[row][i] - estimated_dist ) > devtol || 
00791                  isnan(distances_buf[row][i]) )
00792             {
00793             /*
00794           sinfo_msg("dist=%g devtol=%g isan=%d", 
00795             distances_buf[row][i],
00796             devtol,
00797             isnan(distances_buf[row][i]));
00798             */
00799                 continue ;
00800             }
00801             sum += distances_buf[row][i] ;
00802             n++ ;
00803         }
00804         if ( n < 2 )
00805         {
00806             sinfo_msg_error("distances array could not be determined "
00807                             "completely!, deviations of distances from number "
00808                             "of slitlets too big\n" ) ;
00809             cpl_free(distances) ;
00810             return NULL ;
00811         }
00812         else
00813         {
00814             distances[i] = sum / (float)n ;
00815         }
00816     }
00817     return distances ; 
00818 }
00819 
00843 cpl_imagelist * 
00844 sinfo_new_make_cube ( cpl_image * calibImage,
00845                      float    * distances,
00846                      float    * correct_diff_dist )
00847 {
00848     cpl_imagelist * returnCube ;
00849     int imsize, kslit, kpix ;
00850     int slit_index ;
00851     int z, col, recol ;
00852     int ilx=0;
00853     int ily=0;
00854 
00855     float* podata=NULL;
00856     float* pidata=NULL;
00857     cpl_image* o_img;
00858 
00859     if ( NULL == calibImage )
00860     {
00861         sinfo_msg_error("no resampled image given!\n") ;
00862         return NULL ;
00863     }
00864     ilx=cpl_image_get_size_x(calibImage);
00865     ily=cpl_image_get_size_y(calibImage);
00866     pidata=cpl_image_get_data_float(calibImage);
00867 
00868     if ( NULL == distances )
00869     {
00870         sinfo_msg_error("no distances array from ns_test given!/n") ;
00871         return NULL ;
00872     }
00873 
00874     if ( NULL == correct_diff_dist )
00875     {
00876         sinfo_msg_error("correct_diff_dist array is not allocated!/n") ;
00877         return NULL ;
00878     }
00879        
00880     if ( N_SLITLETS != 32 )
00881     {
00882         sinfo_msg_error ("wrong number of slitlets given \n" ) ;
00883         return NULL ;
00884     }
00885     imsize = ilx / N_SLITLETS ;
00886 
00887     /* allocate memory */  
00888     if ( NULL == (returnCube = cpl_imagelist_new()) )
00889     {
00890         sinfo_msg_error ("cannot allocate new cube \n" ) ;
00891         return NULL ;
00892     }
00893 
00894     /* now build the data cube out of the resampled image */
00895     for ( z = 0 ; z < ily ; z++ ) /* go through the z-axis */
00896     {
00897 
00898       o_img=cpl_image_new(imsize,N_SLITLETS,CPL_TYPE_FLOAT);
00899       podata=cpl_image_get_data_float(o_img);
00900         kpix       = 0 ;
00901         kslit      = 0 ;
00902         slit_index = -1 ;
00903         recol      = -1 ;
00904         for ( col = 0 ; col < ilx ; col++ ) /* go through the image columns */
00905         {
00906             if ( col % imsize == 0 )
00907             {
00908                 recol = 0 ;
00909                 kslit = col/imsize ;
00910                 /* sort the slitlets in the right spiffi specific way */
00911                 switch (kslit)
00912                 {
00913                     case 0:
00914                         slit_index = 8 ;
00915                         break ;
00916                     case 1:
00917                         slit_index = 7 ;
00918                         break ;
00919                     case 2:
00920                         slit_index = 9 ;
00921                         break ;
00922                     case 3:
00923                         slit_index = 6 ;
00924                         break ;
00925                     case 4:
00926                         slit_index = 10 ;
00927                         break ;
00928                     case 5:
00929                         slit_index = 5 ;
00930                         break ;
00931                     case 6:
00932                         slit_index = 11 ;
00933                         break ;
00934                     case 7:
00935                         slit_index = 4 ;
00936                         break ;
00937                     case 8:
00938                         slit_index = 12 ;
00939                         break ;
00940                     case 9:
00941                         slit_index = 3 ;
00942                         break ;
00943                     case 10:
00944                         slit_index = 13 ;
00945                         break ;
00946                     case 11:
00947                         slit_index = 2 ;
00948                         break ;
00949                     case 12:
00950                         slit_index = 14 ;
00951                         break ;
00952                     case 13:
00953                         slit_index = 1 ;
00954                         break ;
00955                     case 14:
00956                         slit_index = 15 ;
00957                         break ;
00958                     case 15:
00959                         slit_index = 0 ;
00960                         break ;
00961                     case 16:
00962                         slit_index = 31 ;
00963                         break ;
00964                     case 17:
00965                         slit_index = 16 ;
00966                         break ;
00967                     case 18:
00968                         slit_index = 30 ;
00969                         break ;
00970                     case 19:
00971                         slit_index = 17 ;
00972                         break ;
00973                     case 20:
00974                         slit_index = 29 ;
00975                         break ;
00976                     case 21:
00977                         slit_index = 18 ;
00978                         break ;
00979                     case 22:
00980                         slit_index = 28 ;
00981                         break ;
00982                     case 23:
00983                         slit_index = 19 ;
00984                         break ;
00985                     case 24:
00986                         slit_index = 27 ;
00987                         break ;
00988                     case 25:
00989                         slit_index = 20 ;
00990                         break ;
00991                     case 26:
00992                         slit_index = 26 ;
00993                         break ;
00994                     case 27:
00995                         slit_index = 21 ;
00996                         break ;
00997                     case 28:
00998                         slit_index = 25 ;
00999                         break ;
01000                     case 29:
01001                         slit_index = 22 ;
01002                         break ;
01003                     case 30:
01004                         slit_index = 24 ;
01005                         break ;
01006                     case 31:
01007                         slit_index = 23 ;
01008                         break ;
01009                     default:
01010                         sinfo_msg_error("wrong slitlet index: couldn't be a "
01011                                "spiffi image,  there must be 32 slitlets!") ;
01012                         cpl_imagelist_delete(returnCube) ;
01013                         return NULL ;
01014                         break ;
01015                 }
01016 
01017                 if ( kslit != 0 )
01018                 {
01019                     /*-------------------------------------------------------- 
01020                      * shift the first pixel by an integer if the absolute 
01021                      * amount of distances[]
01022                      * is bigger than 0.5 
01023                      */
01024                     kpix = sinfo_new_nint(distances[kslit-1]) ;
01025 
01026                     /*----------------------------------------------- 
01027                      * now sort the distances array according to the row order 
01028                      * and add a 0 value for the first (reference) slitlet 
01029                      * that means row 8 
01030                      */
01031                     correct_diff_dist[slit_index] = distances[kslit-1] - 
01032                                                     (float)kpix ;
01033                 }
01034                 /* refer all distances to the first slitlet */
01035                 else
01036                 {
01037                     correct_diff_dist[slit_index] = 0. ;
01038                 }
01039             }
01040 
01041             /* fill each cube plane with one image row */
01042             podata[recol+slit_index*imsize] = pidata[col+kpix+z*ilx];
01043             recol++ ;
01044 
01045             if ( recol > imsize )
01046             {
01047                 sinfo_msg_error("wrong column of reconstructed "
01048                                 "image, shouldn't happen!\n") ; 
01049                 cpl_imagelist_delete(returnCube) ;
01050                 return NULL ;
01051             }
01052         }
01053     }
01054     return returnCube ;
01055 }
01072 cpl_imagelist * 
01073 sinfo_new_make_cube_spi ( cpl_image *  calibImage,
01074                         float    ** slit_edges,
01075                         float    *  shift )
01076 {
01077     cpl_imagelist * returnCube ;
01078     float diff, start ;
01079     float * center ;
01080     int * row_index ;
01081     int slit ;
01082     int col, z ;
01083     int imsize ;
01084     int * beginCol ;
01085     int col_counter ;
01086     int ilx=0;
01087     int ily=0;
01088 
01089     float* podata=NULL;
01090     float* pidata=NULL;
01091     cpl_image* o_img;
01092 
01093 
01094     if ( NULL == calibImage )
01095     {
01096         sinfo_msg_error("no resampled image given!\n") ;
01097         return NULL ;
01098     }
01099     ilx=cpl_image_get_size_x(calibImage);
01100     ily=cpl_image_get_size_y(calibImage);
01101     pidata=cpl_image_get_data_float(calibImage);
01102 
01103     if ( NULL == slit_edges )
01104     {
01105         sinfo_msg_error("no slit_edges array given from sinfo_fitSlits()!/n") ;
01106         return NULL ;
01107     }
01108 
01109     if ( N_SLITLETS != 32 )
01110     {
01111         sinfo_msg_error ("wrong number of slitlets given \n" ) ;
01112         return NULL ;
01113     }
01114     imsize = ilx / N_SLITLETS ;
01115 
01116     /* allocate memory */  
01117     if ( NULL == (row_index = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01118     {
01119         sinfo_msg_error ("cannot allocate memory \n" ) ;
01120         return NULL ;
01121     }
01122     if ( NULL == (beginCol = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01123     {
01124         sinfo_msg_error ("cannot allocate memory \n" ) ;
01125         cpl_free(row_index) ;
01126         return NULL ;
01127     }
01128     if ( NULL == (center = (float*) cpl_calloc(N_SLITLETS, sizeof(float)) ) )
01129     {
01130         sinfo_msg_error ("cannot allocate memory \n" ) ;
01131         cpl_free (row_index) ;
01132         cpl_free (beginCol) ;
01133         return NULL ;
01134     }
01135     if ( NULL == (returnCube = cpl_imagelist_new()) )
01136     {
01137         sinfo_msg_error ("cannot allocate new cube \n" ) ;
01138         cpl_free (row_index) ;
01139         cpl_free (beginCol) ;
01140         cpl_free (center) ;
01141         return NULL ;
01142     }
01143     /* determine the absolute center of the slitlets and the distances 
01144        inside the image*/
01145     for ( slit = 0 ; slit < N_SLITLETS ; slit++ ) 
01146     /* go through the slitlets of each row of the resampled image */
01147     {
01148         center[slit] = (slit_edges[slit][1] + slit_edges[slit][0]) / 2. ;
01149         /* -------------------------------------------------------------
01150          * sort the slitlets in the right spiffi specific way
01151          * the row_index describes the row index of the current slitlet 
01152          * in the resulting cube images.
01153          */
01154         switch (slit)
01155         {
01156             case 0:
01157                 row_index[0] = 8 ;
01158                 break ;
01159             case 1:
01160                 row_index[1] = 7 ;
01161                 break ;
01162             case 2:
01163                 row_index[2] = 9 ;
01164                 break ;
01165             case 3:
01166                 row_index[3] = 6 ;
01167                 break ;
01168             case 4:
01169                 row_index[4] = 10 ;
01170                 break ;
01171             case 5:
01172                 row_index[5] = 5 ;
01173                 break ;
01174             case 6:
01175                 row_index[6] = 11 ;
01176                 break ;
01177             case 7:
01178                 row_index[7] = 4 ;
01179                 break ;
01180             case 8:
01181                 row_index[8] = 12 ;
01182                 break ;
01183             case 9:
01184                 row_index[9] = 3 ;
01185                 break ;
01186             case 10:
01187                 row_index[10] = 13 ;
01188                 break ;
01189             case 11:
01190                 row_index[11] = 2 ;
01191                 break ;
01192             case 12:
01193                 row_index[12] = 14 ;
01194                 break ;
01195             case 13:
01196                 row_index[13] = 1 ;
01197                 break ;
01198             case 14:
01199                 row_index[14] = 15 ;
01200                 break ;
01201             case 15:
01202                 row_index[15] = 0 ;
01203                 break ;
01204             case 16:
01205                 row_index[16] = 31 ;
01206                 break ;
01207             case 17:
01208                 row_index[17] = 16 ;
01209                 break ;
01210             case 18:
01211                 row_index[18] = 30 ;
01212                 break ;
01213             case 19:
01214                 row_index[19] = 17 ;
01215                 break ;
01216             case 20:
01217                 row_index[20] = 29 ;
01218                 break ;
01219             case 21:
01220                 row_index[21] = 18 ;
01221                 break ;
01222             case 22:
01223                 row_index[22] = 28 ;
01224                 break ;
01225             case 23:
01226                 row_index[23] = 19 ;
01227                 break ;
01228             case 24:
01229                 row_index[24] = 27 ;
01230                 break ;
01231             case 25:
01232                 row_index[25] = 20 ;
01233                 break ;
01234             case 26:
01235                 row_index[26] = 26 ;
01236                 break ;
01237             case 27:
01238                 row_index[27] = 21 ;
01239                 break ;
01240             case 28:
01241                 row_index[28] = 25 ;
01242                 break ;
01243             case 29:
01244                 row_index[29] = 22 ;
01245                 break ;
01246             case 30:
01247                 row_index[30] = 24 ;
01248                 break ;
01249             case 31:
01250                 row_index[31] = 23 ;
01251                 break ;
01252             default:
01253                 sinfo_msg_error("wrong slitlet index: couldn't be a spiffi "
01254                                 "image,  there must be 32 slitlets!\n") ;
01255                 cpl_imagelist_delete(returnCube) ;
01256                 cpl_free (row_index) ;
01257                 cpl_free (beginCol) ;
01258                 cpl_free (center) ;
01259                 return NULL ;
01260         }
01261         /* determine the integer column on which the slitlet starts, center the
01262            slitlet on the image row */
01263         start = center[slit] - (float) (imsize - 1)/2. ;
01264         beginCol[slit] = sinfo_new_nint (start) ;
01265         /* determine the error of using integer pixels */
01266         diff = start - (float)beginCol[slit] ;
01267 
01268         /*-------------------------------------------------------------------- 
01269          * determine the output shift values by which the rows are finally 
01270            shifted, consider the integer pixel errors  
01271          * resort shift array to get the row index 
01272          */
01273         shift[row_index[slit]] = diff ;
01274     }   
01275 
01276     /* now build the data cube out of the resampled image */
01277     for ( z = 0 ; z < ily ; z++ ) /* go through the z-axis */
01278     {
01279       o_img=cpl_image_new(imsize,N_SLITLETS,CPL_TYPE_FLOAT);
01280       podata=cpl_image_get_data_float(o_img);
01281       for ( slit = 0 ; slit < N_SLITLETS ; slit++ )
01282         {
01283       col_counter = beginCol[slit] ;
01284       /* each slitlet is centered on the final image row */
01285       for ( col = 0 ; col < imsize ; col++ )
01286             {
01287           if ( col_counter > ilx-1 )
01288                 {
01289           col_counter-- ;
01290                 }
01291           if ( col_counter + z*ilx < 0 )
01292                 {
01293           podata[col+row_index[slit]*imsize] = pidata[0] ;
01294                 }
01295               else
01296                 {   
01297                   podata[col+row_index[slit]*imsize]=pidata[col_counter+z*ilx];
01298                 }
01299 
01300                 col_counter++ ;
01301             }
01302         }
01303       cpl_imagelist_set(returnCube,o_img,z);    
01304     }
01305     cpl_free (row_index) ;
01306     cpl_free (beginCol) ;
01307     cpl_free (center) ;
01308 
01309     return returnCube ;
01310 }
01337 cpl_imagelist * 
01338 sinfo_new_make_cube_dist ( cpl_image * calibImage,
01339                          float      firstCol,
01340                          float    * distances,
01341                          float    * shift )
01342 {
01343     cpl_imagelist * returnCube ;
01344     float di ;
01345     float diff, start ;
01346     int * row_index ;
01347     int slit ;
01348     int col, z ;
01349     int imsize ;
01350     int * beginCol ;
01351     int col_counter ;
01352     int ilx=0;
01353     int ily=0;
01354 
01355     float* podata=NULL;
01356     float* pidata=NULL;
01357     cpl_image* o_img;
01358 
01359     if ( NULL == calibImage )
01360     {
01361         sinfo_msg_error(" no resampled image given!\n") ;
01362         return NULL ;
01363     }
01364     ilx=cpl_image_get_size_x(calibImage);
01365     ily=cpl_image_get_size_y(calibImage);
01366     pidata=cpl_image_get_data_float(calibImage);
01367 
01368     if ( NULL == distances )
01369     {
01370         sinfo_msg_error("no distances array given from north_south_test()!") ;
01371         return NULL ;
01372     }
01373 
01374     if ( N_SLITLETS != 32 )
01375     {
01376         sinfo_msg_error ("wrong number of slitlets given \n" ) ;
01377         return NULL ;
01378     }
01379     imsize = ilx / N_SLITLETS ;
01380 
01381     /* allocate memory */  
01382     if ( NULL == (row_index = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01383     {
01384         sinfo_msg_error ("cannot allocate memory \n" ) ;
01385         return NULL ;
01386     }
01387     if ( NULL == (beginCol = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01388     {
01389         sinfo_msg_error ("cannot allocate memory \n" ) ;
01390         cpl_free(row_index) ;
01391         return NULL ;
01392     }
01393     if ( NULL == (returnCube = cpl_imagelist_new()) )
01394     {
01395         sinfo_msg_error ("cannot allocate new cube \n" ) ;
01396         cpl_free(row_index) ;
01397         cpl_free(beginCol) ;
01398         return NULL ;
01399     }
01400 
01401     di = 0. ;
01402     /* determine the absolute beginning of the slitlets and the distances 
01403        inside the image*/
01404     for ( slit = 0 ; slit < N_SLITLETS ; slit++ ) 
01405     /* go through the slitlets of each row of the resampled image */
01406     {
01407 
01408         /* -------------------------------------------------------------
01409          * sort the slitlets in the right spiffi specific way
01410          * the row_index describes the row index of the current slitlet 
01411          * in the resulting cube images.
01412          */
01413         switch (slit)
01414         {
01415             case 0:
01416                 row_index[0] = 8 ;
01417                 break ;
01418             case 1:
01419                 row_index[1] = 7 ;
01420                 break ;
01421             case 2:
01422                 row_index[2] = 9 ;
01423                 break ;
01424             case 3:
01425                 row_index[3] = 6 ;
01426                 break ;
01427             case 4:
01428                 row_index[4] = 10 ;
01429                 break ;
01430             case 5:
01431                 row_index[5] = 5 ;
01432                 break ;
01433             case 6:
01434                 row_index[6] = 11 ;
01435                 break ;
01436             case 7:
01437                 row_index[7] = 4 ;
01438                 break ;
01439             case 8:
01440                 row_index[8] = 12 ;
01441                 break ;
01442             case 9:
01443                 row_index[9] = 3 ;
01444                 break ;
01445             case 10:
01446                 row_index[10] = 13 ;
01447                 break ;
01448             case 11:
01449                 row_index[11] = 2 ;
01450                 break ;
01451             case 12:
01452                 row_index[12] = 14 ;
01453                 break ;
01454             case 13:
01455                 row_index[13] = 1 ;
01456                 break ;
01457             case 14:
01458                 row_index[14] = 15 ;
01459                 break ;
01460             case 15:
01461                 row_index[15] = 0 ;
01462                 break ;
01463             case 16:
01464                 row_index[16] = 31 ;
01465                 break ;
01466             case 17:
01467                 row_index[17] = 16 ;
01468                 break ;
01469             case 18:
01470                 row_index[18] = 30 ;
01471                 break ;
01472             case 19:
01473                 row_index[19] = 17 ;
01474                 break ;
01475             case 20:
01476                 row_index[20] = 29 ;
01477                 break ;
01478             case 21:
01479                 row_index[21] = 18 ;
01480                 break ;
01481             case 22:
01482                 row_index[22] = 28 ;
01483                 break ;
01484             case 23:
01485                 row_index[23] = 19 ;
01486                 break ;
01487             case 24:
01488                 row_index[24] = 27 ;
01489                 break ;
01490             case 25:
01491                 row_index[25] = 20 ;
01492                 break ;
01493             case 26:
01494                 row_index[26] = 26 ;
01495                 break ;
01496             case 27:
01497                 row_index[27] = 21 ;
01498                 break ;
01499             case 28:
01500                 row_index[28] = 25 ;
01501                 break ;
01502             case 29:
01503                 row_index[29] = 22 ;
01504                 break ;
01505             case 30:
01506                 row_index[30] = 24 ;
01507                 break ;
01508             case 31:
01509                 row_index[31] = 23 ;
01510                 break ;
01511             default:
01512                 sinfo_msg_error("wrong slitlet index: couldn't be a "
01513                                 "spiffi image,  there must be 32 slitlets!") ;
01514                 cpl_imagelist_delete(returnCube) ;
01515                 cpl_free(row_index) ;
01516                 cpl_free(beginCol) ;
01517                 return NULL ;
01518         }
01519 
01520         /* determine the integer column on which the slitlet starts */
01521         if ( slit == 0 )
01522         {
01523             start = firstCol ;
01524         }
01525         else
01526         {
01527             di += distances[slit-1] ;
01528             start = firstCol + di ;
01529         }
01530         beginCol[slit] = sinfo_new_nint(start) ;
01531 
01532         /* determine the error of using integer pixels, its always smaller 
01533            than 1 */
01534         diff = start - (float)beginCol[slit] ;
01535 
01536         /*---------------------------------------------------------------- 
01537          * determine the output shift values by which the rows are finally 
01538          * shifted, consider the integer pixel errors and resort shift array 
01539          * to get the row index 
01540          */
01541         shift[row_index[slit]] = diff ;
01542     }   
01543 
01544     /* now build the data cube out of the resampled image */
01545     for ( z = 0 ; z < ily ; z++ ) /* go through the z-axis */
01546     {
01547       o_img=cpl_image_new(imsize,N_SLITLETS,CPL_TYPE_FLOAT);
01548       podata=cpl_image_get_data_float(o_img);
01549       for ( slit = 0 ; slit < N_SLITLETS ; slit++ )
01550         {
01551       col_counter = beginCol[slit] ;
01552       /* each slitlet is centered on the final image row */
01553       for ( col = 0 ; col < imsize ; col++ )
01554             {
01555           if ( col_counter > ilx-1 )
01556                 {
01557           col_counter-- ;
01558                 }
01559                 if ( col_counter + z*ilx < 0 )
01560                 {
01561           podata[col+row_index[slit]*imsize] = podata[0] ;
01562                 }
01563                 else
01564                 {   
01565                   podata[col+row_index[slit]*imsize]=pidata[col_counter+z*ilx];
01566                 }
01567 
01568                 col_counter++ ;
01569             }
01570         }  
01571       cpl_imagelist_set(returnCube,o_img,z);  
01572     }
01573     cpl_free (row_index) ;
01574     cpl_free (beginCol) ;
01575 
01576     return returnCube ;
01577 }
01604 cpl_imagelist * 
01605 sinfo_new_make_3D_cube_dist ( cpl_image * calibImage,
01606                            float      firstCol,
01607                            float    * distances,
01608                            float    * shift )
01609 {
01610     cpl_imagelist * returnCube ;
01611     float di ;
01612     float diff, start ;
01613     int * row_index ;
01614     int slit ;
01615     int col, z ;
01616     int imsize ;
01617     int * beginCol ;
01618     int col_counter ;
01619     int ilx=0;
01620     int ily=0;
01621 
01622     float* podata=NULL;
01623     float* pidata=NULL;
01624     cpl_image* o_img;
01625 
01626     if ( NULL == calibImage )
01627     {
01628         sinfo_msg_error(" no resampled image given!\n") ;
01629         return NULL ;
01630     }
01631     ilx=cpl_image_get_size_x(calibImage);
01632     ily=cpl_image_get_size_y(calibImage);
01633     pidata=cpl_image_get_data_float(calibImage);
01634 
01635     if ( NULL == distances )
01636     {
01637         sinfo_msg_error("no distances array given from north_south_test()!") ;
01638         return NULL ;
01639     }
01640 
01641     if ( N_SLITLETS != 16 )
01642     {
01643         sinfo_msg_error ("wrong number of slitlets given \n" ) ;
01644         return NULL ;
01645     }
01646     imsize = ilx / N_SLITLETS ;
01647 
01648     /* allocate memory */  
01649     if ( NULL == (row_index = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01650     {
01651         sinfo_msg_error ("cannot allocate memory \n" ) ;
01652         return NULL ;
01653     }
01654     if ( NULL == (beginCol = (int*) cpl_calloc(N_SLITLETS, sizeof(int)) ) )
01655     {
01656         sinfo_msg_error ("cannot allocate memory \n" ) ;
01657         cpl_free(row_index) ;
01658         return NULL ;
01659     }
01660     if ( NULL == (returnCube = cpl_imagelist_new()) )
01661     {
01662         sinfo_msg_error ("cannot allocate new cube \n" ) ;
01663         cpl_free(row_index) ;
01664         cpl_free(beginCol) ;
01665         return NULL ;
01666     }
01667 
01668     di = 0. ;
01669     /* determine the absolute beginning of the slitlets and the distances 
01670        inside the image*/
01671     for ( slit = 0 ; slit < N_SLITLETS ; slit++ ) 
01672     /* go through the slitlets of each row of the resampled image */
01673     {
01674 
01675         /* --------------------------------------------------------------
01676          * sort the slitlets in the right 3D specific way
01677          * the row_index describes the row index of the current slitlet 
01678          * in the resulting cube images.
01679          */
01680         row_index[slit] = slit ;
01681 
01682         /* determine the integer column on which the slitlet starts */
01683         if ( slit == 0 )
01684         {
01685             start = firstCol ;
01686         }
01687         else
01688         {
01689             di += distances[slit-1] ;
01690             start = firstCol + di ;
01691         }
01692         beginCol[slit] = sinfo_new_nint(start) ;
01693 
01694         /* determine the error of using integer pixels, 
01695            `its always smaller than 1 */
01696         diff = start - (float)beginCol[slit] ;
01697 
01698         /*---------------------------------------------------------------- 
01699          * determine the output shift values by which the rows are finally 
01700            shifted, consider the integer pixel errors and resort shift array 
01701            to get the row index 
01702          */
01703         shift[row_index[slit]] = diff ;
01704     }   
01705 
01706     /* now build the data cube out of the resampled image */
01707     for ( z = 0 ; z < ily ; z++ ) /* go through the z-axis */
01708     {
01709       o_img=cpl_image_new(imsize,N_SLITLETS,CPL_TYPE_FLOAT);
01710       podata=cpl_image_get_data_float(o_img);
01711         for ( slit = 0 ; slit < N_SLITLETS ; slit++ )
01712         {
01713             col_counter = beginCol[slit] ;
01714             /* each slitlet is centered on the final image row */
01715             for ( col = 0 ; col < imsize ; col++ )
01716             {
01717                 if ( col_counter > ilx-1 )
01718                 {
01719                     col_counter-- ;
01720                 }
01721                 podata[col+row_index[slit]*imsize]=pidata[col_counter+z*ilx];
01722                 col_counter++ ;
01723             }
01724         } 
01725         cpl_imagelist_set(returnCube,o_img,z);   
01726     }
01727     cpl_free (row_index) ;
01728     cpl_free (beginCol) ;
01729 
01730     return returnCube ;
01731 }
01732 
01750 cpl_imagelist * 
01751 sinfo_new_make_3D_cube ( cpl_image * calibImage,
01752                        int      * kpixshift, 
01753                        int        kpixfirst )
01754 {
01755     cpl_imagelist * returnCube ;
01756     int imsize, kslit, kpix ;
01757     int z, col, recol ;
01758     int ilx=0;
01759     int ily=0;
01760 
01761     float* podata=NULL;
01762     float* pidata=NULL;
01763     cpl_image* o_img;
01764 
01765     if ( NULL == calibImage )
01766     {
01767         sinfo_msg_error("no resampled image given!\n") ;
01768         return NULL ;
01769     }
01770     ilx=cpl_image_get_size_x(calibImage);
01771     ily=cpl_image_get_size_y(calibImage);
01772     pidata=cpl_image_get_data_float(calibImage);
01773 
01774     if ( NULL == kpixshift )
01775     {
01776         sinfo_msg_error("no shift array given!/n") ;
01777         return NULL ;
01778     }
01779 
01780     if ( kpixfirst < 0 )
01781     {
01782         sinfo_msg_error("wrong first valid pixel given!/n") ;
01783         return NULL ;
01784     }
01785 
01786     if ( N_SLITLETS != 16 )
01787     {
01788         sinfo_msg_error ("wrong number of slitlets given \n" ) ;
01789         return NULL ;
01790     }
01791     imsize = ilx / N_SLITLETS ;
01792 
01793     if ( NULL == (returnCube = cpl_imagelist_new()) )
01794     {
01795         sinfo_msg_error ("cannot allocate new cube \n" ) ;
01796         return NULL ;
01797     }
01798 
01799     /* now build the data cube out of the resampled image */
01800     for ( z = 0 ; z < ily ; z++ ) /* go through the z-axis */
01801     {
01802       o_img=cpl_image_new(imsize,N_SLITLETS,CPL_TYPE_FLOAT);
01803       podata=cpl_image_get_data_float(o_img);
01804         kpix       = 0 ;
01805         kslit      = 0 ;
01806         recol      = -1 ;
01807         for ( col = 0 ; col < ilx ; col++ ) /* go through the image columns */
01808         {
01809             if ( col % imsize == 0 ) 
01810             {
01811                 recol = 0 ;
01812                 kslit = col/imsize ;
01813                 kpix  = kpixfirst + kpixshift[kslit] ;
01814             }
01815 
01816             /* fill each cube plane with one image row */
01817             podata[recol+kslit*imsize] = pidata[col+kpix+z*ilx] ;
01818             recol++ ;
01819             if ( recol > imsize )
01820             {
01821                 sinfo_msg_error("wrong column of reconstructed image, i"
01822                                 "shouldn't happen!\n") ; 
01823                 cpl_imagelist_delete(returnCube) ;
01824                 return NULL ;
01825             }
01826         }
01827         cpl_imagelist_set(returnCube,o_img,z);
01828     }
01829     return returnCube ;
01830 }
01831 
01844 cpl_imagelist * 
01845 sinfo_new_determine_mask_cube ( cpl_imagelist * sourceMaskCube,
01846                               float     lowLimit,
01847                               float     highLimit )
01848 {
01849     cpl_imagelist * retCube ; 
01850     int z, n ;
01851     int ilx=0;
01852     int ily=0;
01853     int inp=0;
01854     int olx=0;
01855     int oly=0;
01856     int onp=0;
01857     float* podata=NULL;
01858     cpl_image* o_img;
01859 
01860     if ( sourceMaskCube == NULL )
01861     {
01862         sinfo_msg_error("no cube given!\n") ;
01863         return NULL ;
01864     }
01865     ilx=cpl_image_get_size_x(cpl_imagelist_get(sourceMaskCube,0));
01866     ily=cpl_image_get_size_y(cpl_imagelist_get(sourceMaskCube,0));
01867     inp=cpl_imagelist_get_size(sourceMaskCube);
01868 
01869 
01870     if ( lowLimit > 0. )
01871     {
01872         sinfo_msg_error("lowLimit wrong!\n") ;
01873         return NULL ;
01874     }
01875     if ( highLimit >= 1. || highLimit < 0. )
01876     {
01877         sinfo_msg_error("highLimit wrong!\n") ;
01878         return NULL ;
01879     }
01880 
01881     retCube = cpl_imagelist_duplicate (sourceMaskCube) ;
01882     onp=inp;
01883     olx=ilx;
01884     oly=ily;
01885 
01886     for ( z = 0 ; z < onp ; z++ )
01887     {
01888       o_img=cpl_imagelist_get(retCube,0);
01889       podata=cpl_image_get_data_float(o_img);
01890         for ( n = 0 ; n < (int) olx*oly; n++ )
01891         {
01892             if ( podata[n] == 0. )
01893             {
01894                continue ;
01895             }
01896             if ( podata[n] == 1. )
01897             {
01898                continue ;
01899             }
01900             if ( podata[n] >= lowLimit && 
01901                  podata[n] <= highLimit )
01902             {
01903                 podata[n] = 0. ;
01904             }
01905             else 
01906             {
01907                 podata[n] = 1. ;
01908             }
01909         }
01910     }
01911     return retCube ;
01912 }
01951 cpl_imagelist * 
01952 sinfo_new_interpol_cube ( cpl_imagelist * sourceCube,
01953                          cpl_imagelist * maskCube,
01954                          int       n_neighbors, /* 7 */
01955                          int       max_radius ) /* 5 */
01956 {
01957     cpl_imagelist  * returnCube ;
01958     float** spec=NULL ;
01959     float* spec1=NULL ;
01960     int n_im, n_bad, n_bad1, n_bad2 ;
01961     int n_planes, specn, nspec1 ;
01962     int i, m, n, z, ni, kk, p ;
01963     int dis, dismin, dismax ;
01964     int agreed ;
01965     int xcordi, ycordi, xcordm, ycordm ;
01966 
01967 
01968 
01969     int ilx=0;
01970     int ily=0;
01971     int inp=0;
01972    
01973     float* pidata=NULL;
01974     float* pmdata=NULL;
01975     float* podata=NULL;
01976     cpl_image* i_img=NULL;
01977     cpl_image* m_img=NULL;
01978     cpl_image* o_img=NULL;
01979 
01980     if ( NULL == sourceCube )
01981     {
01982         sinfo_msg_error(" no source cube given!\n") ;
01983         return NULL ;
01984     }
01985 
01986 
01987     ilx=cpl_image_get_size_x(cpl_imagelist_get(sourceCube,0));
01988     ily=cpl_image_get_size_y(cpl_imagelist_get(sourceCube,0));
01989     inp=cpl_imagelist_get_size(sourceCube);
01990 
01991     if ( NULL == maskCube )
01992     {
01993         sinfo_msg_error("no bad pixel mask cube given!\n") ;
01994         return NULL ;
01995     }
01996 
01997     if ( n_neighbors <= 0 )
01998     {
01999         sinfo_msg_error("wrong number of neighbors in the spectral "
02000                         "direction given!") ;
02001         return NULL ;
02002     }
02003    
02004     if ( max_radius <= 0 )
02005     {
02006         sinfo_msg_error("wrong maximal radius for interpolation inside "
02007                         "an image plane given!") ;
02008         return NULL ;
02009     }
02010 
02011     returnCube = cpl_imagelist_duplicate(sourceCube) ;
02012     
02013     n_im     = ilx * ily ;
02014     n_planes = inp ;
02015 
02016     spec1=cpl_calloc(300,sizeof(float)) ;
02017     spec=sinfo_new_2Dfloatarray(100,2*n_neighbors+1) ;
02018 
02019     /* loop over the image planes and look for bad pixels and correct them */
02020     for ( z = 0 ; z < n_planes ; z++ ) /* go through image planes */
02021     {
02022       m_img=cpl_imagelist_get(maskCube,z);
02023       pmdata=cpl_image_get_data_float(m_img);
02024       o_img=cpl_imagelist_get(returnCube,z);
02025       podata=cpl_image_get_data_float(o_img);
02026 
02027         /*-------------------------------------------------------------------
02028          * determine n, the length of one wing in one spectrum with which the 
02029          * bad pixel will be interpolated. The length of one wing is 
02030            n_neighbors but less at the edges of the cube. 
02031          */
02032         if ( z < n_neighbors )
02033         {
02034             n = z ;
02035         }
02036         else if ( n_planes - z <= n_neighbors)
02037         {
02038             n = n_planes - z -1 ;
02039         }
02040         else
02041         { 
02042             n = n_neighbors ;
02043         }
02044 
02045         for ( i = 0 ; i < n_im ; i ++ ) /* go through one image */
02046         {
02047             /* continue if the pixel is a good one */
02048             if ( pmdata[i] != 0. )
02049             {
02050                 continue ;
02051             }
02052         
02053             /*-------------------------------------------------------------
02054              * exclude pixels with too many bad neighbors in the spectrum.
02055              * exit if: too few good pixels in the neighboring spectrum or 
02056              * good pixels are only on one side of the spectrum.
02057              */
02058             n_bad  = 0 ;
02059             n_bad1 = 0 ;
02060             n_bad2 = 0 ;
02061             /* go through the neighbor spectral pixels */
02062             for ( ni = z-n ; ni <= z+n ; ni++ ) 
02063             {
02064                 if ( pmdata[i] == 0. )
02065                 {
02066                     n_bad++ ;
02067                     /* count bad pixels on either spectral side of 
02068                        the bad pixel to be interpolated */
02069                     if ( ni < z )
02070                     {
02071                         n_bad1++ ;
02072                     }
02073                     if ( ni > z )
02074                     {
02075                         n_bad2++ ;
02076                     }
02077                 }
02078             }
02079          
02080             /*--------------------------------------------------------------- 
02081              * now the criteria are checked which the neighborhood in the 
02082                spectral dimension has to match if the pixel is interpolatable.
02083              * The total number of the good pixel in the spectrum must be more 
02084                than 3 and there must be at least one good pixel on either side 
02085                of the central pixel.
02086              */
02087             if ( (2*n+1 - n_bad) < 3 || (n - n_bad1) < 1 || (n - n_bad2) < 1 )
02088             {
02089                 continue ;
02090             }
02091             
02092             /* read the master spectrum into the first row of the array spec */
02093             kk = 0 ;
02094             for ( ni = z-n ; ni <= z+n ; ni++ )
02095             {
02096           i_img=cpl_imagelist_get(sourceCube,ni);
02097               pidata=cpl_image_get_data_float(i_img);
02098                 spec[1][kk] = pmdata[i] != 0. ? pidata[i] : ZERO ;
02099                 kk++ ; /* length of spectrum */
02100             }
02101             
02102             /* look for appropriate neighbors in the x-y neighborhood */
02103             agreed = 1 ; /* loop guard */
02104             specn  = 2 ; /* number of spectra in spec. 
02105                             First is master spectrum */
02106             dismin = 0 ; /* x+y minimal distance to bad pixel */
02107             dismax = 1 ; /* x+y maximal distance to bad pixel */
02108             do
02109             {
02110                 for ( m = 0 ; m < n_im ; m++ )
02111                 {
02112                     if ( pmdata[m] == 0. )
02113                     {
02114                         continue ;
02115                     }
02116 
02117                     /* --------------------------------------------------------
02118                      * determine the x and y coordinates of the bad pixel (i)
02119                      * and the pixels used to interpolate (m) 
02120                      */
02121                     xcordi = i % ilx ;
02122                     xcordm = m % ilx ;
02123                     ycordi = i / ilx ;
02124                     ycordm = m / ilx ;
02125                     /*----------------------------------------------------- 
02126                      * check the distance: take only close pixels
02127                      * extension 'i' is coordinate of the bad pixel to be 
02128                        interpolated
02129                      */
02130                     dis = abs(xcordi-xcordm) + abs(ycordi-ycordm) ;
02131                     if ( dis <= dismin || dis > dismax )
02132                     {
02133                         continue ;
02134                     }
02135                     /*--------------------------------------------------------
02136                      * check on number of bad pixels in the spectrum of a 
02137                      * neighbor pixel; reject it if it contains less than 2 
02138                      * usable pixel pairs. a bit more explanation:
02139                      * let this be a 15 pixel spectrum with the pixel to be 
02140                      * interpolated denoted by '0' and other bad pixels marked 
02141                      * with 'b'. Good pixels are marked with '1'. Below a 
02142                      * neighbor spectrum is drawn containing bad pixels as 
02143                      * well. The third line shows the position of the usable 
02144                      * pixel pairs, spectral
02145                      * positions, where both spectra have valid pixels.
02146                      *
02147                      *   1 1 1 b b 1 1 0 b 1 b b 1 b b
02148                      *   b 1 1 1 b b 1 1 1 1 1 1 b b 1
02149                      *     ^ ^       ^     ^             4 good pixel pairs
02150                      */
02151       
02152                     n_bad = 0 ;
02153                     for ( ni = z-n ; ni <= z+n ; ni++ )
02154                     {
02155                         if ( pmdata[i] == 0. || pmdata[m] == 0. )
02156                         {
02157                             n_bad++ ;
02158                         }
02159                     }
02160                     if ( n_bad > 2*n-1 ) 
02161                     /* we need at least 2 usable pixel pairs */
02162                     {
02163                         continue ;
02164                     }
02165                     
02166                     /* transfer the spectrum to the next position 
02167                        of array spec */
02168                     kk = 0 ;
02169                     for ( ni = z-n ; ni <= z+n ; ni++ )
02170                     {
02171               i_img=cpl_imagelist_get(sourceCube,ni);
02172               pidata=cpl_image_get_data_float(i_img);
02173                         spec[specn][kk] = pmdata[m] != 0. ? pidata[m] : ZERO ;
02174                         kk++ ;
02175                     }
02176                     specn++ ;
02177                     if ( specn > 10 ) /* if we have 9 neighbors then break */
02178                     {
02179                         agreed = 0 ;
02180                         break ;
02181                     }
02182                 }
02183                 /* if no break, increase search radius and continue */
02184                 dismin++ ;        
02185                 dismax++ ;
02186                 /* if search radius is too big, exit with fewer good neighbors */
02187                 if ( dismax > max_radius )
02188                 {
02189                     agreed = 0 ;
02190                 }
02191             }   while(agreed) ;    
02192                        
02193             specn-- ;
02194             dismax -= 2 ;
02195 
02196             /* ---------------------------------------------------------------
02197              * Take the master spectrum with the bad pixel in the middle and 
02198                divide it by each of the neighbor spectra and normalize the 
02199                division with the value in  the center position.
02200              */
02201             for ( kk = 0 ; kk < 2*n+1 ; kk++ )
02202             {
02203                 if ( kk == n )    /* do not divide the master bad pixel */
02204                 {
02205                     continue ;
02206                 }
02207 
02208                 /* do not divide bad pixels in the master spectrum */
02209                 if ( isnan(spec[1][kk]) ) 
02210                 {
02211                     for ( p = 2 ; p <= specn ; p++ )
02212                     {
02213                         spec[p][kk] = ZERO ;
02214                     }
02215                 }    
02216                 else       /* all is well, now divide */
02217                 {
02218                     for ( p = 2 ; p <= specn ; p++ )
02219                     {
02220                         if ( !isnan(spec[p][kk]) && spec[p][kk] != 0. &&
02221                              !isnan(spec[p][n]) )
02222                         {
02223                             spec[p][kk] = spec[1][kk] / 
02224                                           spec[p][kk] * spec[p][n] ;
02225                         }
02226                         else
02227                         {
02228                             spec[p][kk] = ZERO ;
02229                         }
02230                     }
02231                 }
02232             }
02233  
02234             /*-----------------------------------------------------------------
02235              * determine the sinfo_median of all values. With 9 good neighbors 
02236              * and at least 2 good values per neighbor we have between 18 and 
02237              * 9*14 values for the statistics. If there are not enough good 
02238              * neighbors available, only continue if we have collected at 
02239              * least 18 values.
02240              */
02241             nspec1 = 0 ;  
02242             /* collect the good values in the array spec1 */ 
02243             for ( p = 2 ; p <= specn ; p++ )
02244             {
02245                 for ( kk = 0 ; kk < 2*n+1 ; kk++ )
02246                 {
02247                     if ( !isnan(spec[p][kk]) && kk != n )
02248                     {
02249                         spec1[nspec1] = spec[p][kk] ;
02250                         nspec1++ ;
02251                     }
02252                 }
02253             }
02254             
02255             /* now test if we have at least 18 values */
02256             if ( nspec1 < 18 )
02257             {
02258                 continue ;
02259             }
02260  
02261             /* interpolate the bad pixel by the sinfo_median of spec1 */
02262             podata[i] = sinfo_new_median(spec1, nspec1) ;
02263             pmdata[i] = 1 ;
02264         }
02265     }           
02266     sinfo_new_destroy_2Dfloatarray(&spec,2*n_neighbors+1) ;
02267     cpl_free(spec1);
02268     return returnCube ;
02269 }
02290 cpl_imagelist * 
02291 sinfo_new_fine_tune_cube( cpl_imagelist * cube,
02292                                      float   * correct_diff_dist,
02293                                      int       n_order )
02294 {
02295     cpl_imagelist * returnCube ;
02296     float* row_data=NULL ;
02297     float* corrected_row_data=NULL ;
02298     float* xnum=NULL ;
02299     float sum, new_sum ;
02300     float eval/*, dy*/ ;
02301     float * imageptr ;
02302     int row, col ;
02303     int i, z ;
02304     int imsize, n_points ;
02305     int firstpos ;
02306     int  flag;
02307     int ilx=0;
02308     int ily=0;
02309     int inp=0;
02310    
02311     float* pidata=NULL;
02312     float* podata=NULL;
02313     cpl_image* i_img=NULL;
02314     cpl_image* o_img=NULL;
02315 
02316 
02317     if ( NULL == cube )
02318     {
02319         sinfo_msg_error("no input cube given!\n") ;
02320         return NULL ;
02321     }
02322     ilx=cpl_image_get_size_x(cpl_imagelist_get(cube,0));
02323     ily=cpl_image_get_size_y(cpl_imagelist_get(cube,0));
02324     inp=cpl_imagelist_get_size(cube);
02325 
02326     if ( NULL == correct_diff_dist )
02327     {
02328         sinfo_msg_error("no distances array from ns_test given!n") ;
02329         return NULL ;
02330     }
02331   
02332     if ( n_order <= 0 )
02333     {
02334         sinfo_msg_error("wrong order of interpolation polynom given!") ;
02335     returnCube = cpl_imagelist_duplicate(cube);
02336         return returnCube ;
02337     }
02338 
02339     returnCube = cpl_imagelist_duplicate(cube);
02340     
02341     imsize = ily ;
02342     if ( imsize != N_SLITLETS )
02343     {
02344         sinfo_msg_error ("wrong image size\n" ) ;
02345         return NULL ;
02346     }
02347 
02348     n_points = n_order + 1 ;
02349     if ( n_points % 2 == 0 )
02350     {
02351         firstpos = (int)(n_points/2) - 1 ;
02352     }
02353     else
02354     {
02355         firstpos = (int)(n_points/2) ;
02356     }
02357     xnum=cpl_calloc(n_order+1,sizeof(float)) ;
02358 
02359     for ( i = 0 ; i < n_points ; i++ )
02360     {
02361         xnum[i] = i ;
02362     }    
02363 
02364     row_data=cpl_calloc(ilx,sizeof(float)) ;
02365     corrected_row_data=cpl_calloc(ilx,sizeof(float)) ;
02366 
02367     for ( z = 0 ; z < inp ; z++ )
02368     {
02369       i_img=cpl_imagelist_get(cube,z);
02370       pidata=cpl_image_get_data_float(i_img);
02371       o_img=cpl_imagelist_get(returnCube,z);
02372       podata=cpl_image_get_data_float(o_img);
02373 
02374 
02375         for ( row = 0 ; row < imsize ; row++ )
02376         {
02377             for ( col = 0 ; col < ilx ; col++ )
02378             {
02379                 corrected_row_data[col] = 0. ;
02380             }
02381             sum = 0. ; 
02382             for ( col = 0 ; col < ilx ; col++ )
02383             {
02384                 row_data[col] = pidata[col+row*ilx] ;
02385                 if ( isnan(row_data[col]) )
02386                 {
02387                     row_data[col] = 0. ;
02388                     for ( i = col - firstpos ; 
02389                           i < col -firstpos+n_points ; i++ )
02390                     {
02391                         if ( i < 0 ) continue ;
02392                         if ( i >= ilx) continue ; 
02393                         corrected_row_data[i] = ZERO ;
02394                     }
02395                 }
02396                 if ( col != 0 && col != ilx - 1 )
02397                 {
02398                     sum += row_data[col] ;
02399                 }
02400             }
02401 
02402            
02403             new_sum = 0. ;
02404             for ( col = 0 ; col < ilx ; col++ )
02405             {
02406                 
02407                 if ( isnan(corrected_row_data[col]) )
02408                 {
02409                     continue ;
02410                 }
02411                 if ( col - firstpos < 0 )
02412                 {
02413                     imageptr = &row_data[0] ;
02414                     eval     = correct_diff_dist[row] + col ;
02415                 }
02416                 else if ( col - firstpos + n_points >= ilx )
02417                 {
02418                     imageptr = &row_data[ilx - n_points] ;
02419                     eval     = correct_diff_dist[row] + col + n_points - ilx ;
02420                 }
02421                 else
02422                 {
02423                     imageptr = &row_data[col-firstpos] ;
02424                     eval     = correct_diff_dist[row] + firstpos ;
02425                 }
02426 
02427         
02428         flag = 0;
02429         corrected_row_data[col]=sinfo_new_nev_ille(xnum, imageptr, 
02430                                                        n_order, eval, &flag);
02431 
02432                
02433                 if ( col != 0 && col != ilx - 1 )
02434                 {
02435                     new_sum += corrected_row_data[col] ;
02436                 }
02437             }
02438             for ( col = 0 ; col < ilx ; col++ )
02439             {
02440                 
02441                 if ( col == 0 )
02442                 {
02443                     podata[col+row*ilx] = ZERO ;
02444                 }
02445                 else if ( col == ilx - 1 )
02446                 {
02447                     podata[col+row*ilx] = ZERO ;
02448                 }
02449                 else
02450                 {
02451                     if ( isnan(corrected_row_data[col]) ) 
02452                     {
02453                         podata[col+row*ilx] = ZERO ;
02454                     }
02455                     else
02456                     {
02457                         if ( new_sum == 0. ) new_sum = 1. ;
02458                      
02459                         podata[col+row*ilx] = corrected_row_data[col] ;
02460                     }
02461                 }
02462             }
02463         }
02464     }       
02465 
02466     cpl_free(xnum) ;
02467     cpl_free(row_data) ;
02468     cpl_free(corrected_row_data) ;
02469 
02470     return returnCube ;
02471 }
02472 
02492 cpl_imagelist * 
02493 sinfo_new_fine_tune_cube_by_FFT( cpl_imagelist * cube,
02494                                            float   * correct_diff_dist )
02495 {
02496     cpl_imagelist * returnCube ;
02497 
02498     float* row_data=NULL ;
02499     dcomplex* data=NULL ;
02500     dcomplex* corrected_data=NULL ;
02501 
02502     unsigned nn[1];
02503     /*float corrected_row_data[cube->lx] ;*/
02504     float phi, pphi ;
02505     float coph, siph ;
02506     int row, col ;
02507     int i, z ;
02508     int imsize ;
02509     int blank_indicator ;
02510 
02511 
02512     int ilx=0;
02513     int ily=0;
02514     int inp=0;
02515    
02516     float* pidata=NULL;
02517     float* podata=NULL;
02518     cpl_image* i_img=NULL;
02519     cpl_image* o_img=NULL;
02520 
02521 
02522 
02523     if ( NULL == cube )
02524     {
02525         sinfo_msg_error(" no input cube given!\n") ;
02526         return NULL ;
02527     }
02528     ilx=cpl_image_get_size_x(cpl_imagelist_get(cube,0));
02529     ily=cpl_image_get_size_y(cpl_imagelist_get(cube,0));
02530     inp=cpl_imagelist_get_size(cube);
02531 
02532     nn[1] = ilx ;
02533     if ( NULL == correct_diff_dist )
02534     {
02535         sinfo_msg_error("no distances array from ns_test given!") ;
02536         return NULL ;
02537     }
02538 
02539     returnCube = cpl_imagelist_duplicate( cube ) ;
02540 
02541     imsize = ily ;
02542     if ( imsize != N_SLITLETS )
02543     {
02544         sinfo_msg_error ("wrong image size\n" ) ;
02545         return NULL ;
02546     }
02547 
02548     data=cpl_calloc(ilx,sizeof(dcomplex)) ;
02549     corrected_data=cpl_calloc(ilx,sizeof(dcomplex)) ;
02550 
02551     row_data=cpl_calloc(ilx,sizeof(float)) ;
02552     /* loop over the image planes */
02553     for ( z = 0 ; z < inp ; z++ )
02554     {
02555       i_img=cpl_imagelist_get(cube,z);
02556       pidata=cpl_image_get_data_float(i_img);
02557       o_img=cpl_imagelist_get(returnCube,z);
02558       podata=cpl_image_get_data_float(o_img);
02559         /* consider one row at a time */
02560         for ( row = 0 ; row < imsize ; row++ )
02561         {
02562             blank_indicator = 1 ;
02563             for ( col = 0 ; col < ilx ; col++ )
02564             {
02565                 /* transfer the row data to a double sized array */
02566                 row_data[col] = pidata[col+row*ilx] ;
02567         data[col].x = row_data[col] ;
02568         data[col].y = 0. ;
02569                 /* if row contains a blank pixel proceed */
02570                 if ( isnan(row_data[col]) )
02571                 {
02572                     blank_indicator = 0 ; 
02573                 }
02574             }
02575 
02576             /* if row contains a blank don't apply FFT but proceed */
02577             if ( blank_indicator == 0 )
02578             {
02579                 for ( col = 0 ; col < ilx ; col++ )
02580                 {
02581                     podata[col+row*ilx] = ZERO ;
02582                 }
02583                 continue ;
02584             }
02585            
02586             /* FFT algorithm of eclipse */
02587             sinfo_fftn( data, nn, 1, 1 ) ;
02588  
02589             /* calculate the corrected phase shift for each frequency */
02590             phi = 2*PI_NUMB/(float)ilx * correct_diff_dist[row] ;
02591             for ( i = 0 ; i < ilx ; i++ )
02592             {
02593                 /* positive frequencies */
02594                 if ( i <= ilx/2 )
02595                 {
02596                     /* phase shift */
02597                     pphi = phi * (float)(i) ;
02598                     /* Euler factor */
02599                     coph = cos ( pphi ) ;
02600                     siph = sin ( pphi ) ;
02601                 }
02602                 else /* negative frequencies */
02603                 {
02604                     /* phase shift */
02605                     pphi = phi * (float)(i - ilx/2) ;
02606                     /* Euler factor */
02607                     coph = cos ( pphi ) ;
02608                     siph = sin ( pphi ) ;
02609                 }
02610 
02611                 /* ------------------------------------------------------------
02612                  * now calculate the shift in the pixel space by multiplying
02613                  * the fourier transform by the Euler factor of the phase shift
02614                  * and inverse fourier transforming.
02615                  * used Fourier pair: h(x-x0) <==> H(k)*exp(2*pi*i*k*x0) 
02616                  */
02617                 /* calculate real part */
02618                 corrected_data[i].x   = data[i].x * coph - data[i].y * siph ; 
02619                 /* calculate imaginary part */
02620                 corrected_data[i].y = data[i].x * siph + data[i].y * coph ;
02621             }
02622  
02623             /* transform back: inverse FFT */
02624             sinfo_fftn( corrected_data, nn, 1, -1 ) ;
02625 
02626             /* normalize */ 
02627             for ( i = 0 ; i < ilx ; i++ )
02628             {
02629                 corrected_data[i].x /= ilx ;
02630         corrected_data[i].y /= ilx ;
02631             } 
02632 
02633             /* now transfer row to output, leave the left-most 
02634                and right-most pixel column */
02635             for ( col = 0 ; col < ilx ; col++ )
02636             {
02637                 if ( col == 0 )
02638                 {
02639                     podata[col+row*ilx] = ZERO ;
02640                 }
02641                 else if ( col == ilx - 1 )
02642                 {
02643                     podata[col+row*ilx] = ZERO ;
02644                 }
02645                 else
02646                 {
02647                     podata[col+row*ilx] = corrected_data[col].x ; 
02648                 }
02649             }
02650         }
02651     }
02652 
02653     cpl_free(data) ;
02654     cpl_free(corrected_data) ;
02655 
02656 
02657     cpl_free(row_data);
02658     return returnCube ;
02659 }
02687 cpl_imagelist * sinfo_new_fine_tune_cube_by_spline ( cpl_imagelist * cube,
02688                                                float   * correct_diff_dist )
02689 {
02690     cpl_imagelist * returnCube ;
02691 
02692     float* row_data=NULL ;
02693     float* corrected_row_data=NULL ;
02694     float* xnum=NULL ;
02695     float* eval=NULL ;
02696 
02697     float sum, new_sum ;
02698     int row, col ;
02699     int i, z ;
02700     int imsize ;
02701     int ilx=0;
02702     int ily=0;
02703     int inp=0;
02704    
02705     float* pidata=NULL;
02706     float* podata=NULL;
02707     cpl_image* i_img=NULL;
02708     cpl_image* o_img=NULL;
02709 
02710  
02711     if ( NULL == cube )
02712     {
02713         sinfo_msg_error("no input cube given!\n") ;
02714         return NULL ;
02715     }
02716     ilx=cpl_image_get_size_x(cpl_imagelist_get(cube,0));
02717     ily=cpl_image_get_size_y(cpl_imagelist_get(cube,0));
02718     inp=cpl_imagelist_get_size(cube);
02719 
02720     if ( NULL == correct_diff_dist )
02721     {
02722         sinfo_msg_error("no distances array from ns_test given!/n") ;
02723         return NULL ;
02724     }
02725   
02726     imsize = ily ;
02727     if ( imsize != N_SLITLETS )
02728     {
02729         sinfo_msg_error ("wrong image size\n" ) ;
02730         return NULL ;
02731     }
02732   
02733     returnCube = cpl_imagelist_duplicate( cube ) ;
02734 
02735     row_data=cpl_calloc(ilx,sizeof(float)) ;
02736     corrected_row_data=cpl_calloc(ilx,sizeof(float)) ;
02737     xnum=cpl_calloc(ilx,sizeof(float)) ;
02738     eval=cpl_calloc(ilx,sizeof(float)) ;
02739   
02740     /* fill the xa[] array for a polynomial interpolation */
02741     for ( i = 0 ; i < ilx ; i++ )
02742     {
02743         xnum[i] = i ;
02744     }    
02745 
02746     /* loop over the image planes */
02747     for ( z = 0 ; z < inp ; z++ )
02748     {
02749       i_img=cpl_imagelist_get(cube,z);
02750       pidata=cpl_image_get_data_float(i_img);
02751       o_img=cpl_imagelist_get(returnCube,z);
02752       podata=cpl_image_get_data_float(o_img);
02753         /* consider 1 row at a time */
02754         for ( row = 0 ; row < imsize ; row++ )
02755         {
02756             for ( col = 0 ; col < ilx ; col++ )
02757             {
02758                 corrected_row_data[col] = 0. ;
02759             }
02760         sum = 0. ; /* initialize flux for later rescaling */
02761             /* go through the columns and compute the flux for each 
02762                row (leave the sinfo_edge points) */
02763             for ( col = 0 ; col < ilx ; col++ )
02764             {   
02765             eval[col] = correct_diff_dist[row] + (float)col ;
02766                 row_data[col] = pidata[col+row*ilx] ;
02767                 if (col != 0 && col != ilx - 1 && !isnan(row_data[col]) )
02768                 {
02769                     sum += row_data[col] ;
02770                 }
02771                 if (isnan(row_data[col]) )
02772                 {
02773                     for ( i = col -1 ; i <= col+1 ; i++ ) 
02774                     {
02775                         if ( i < 0 ) continue ;
02776                         if ( i >= ilx ) continue ;
02777                         corrected_row_data[i] = ZERO ; 
02778                     }
02779                     row_data[col] = 0. ;
02780                 }
02781         }
02782 
02783 
02784             /* ---------------------------------------------------------------
02785              * now we do the cubic spline interpolation to achieve the 
02786                fractional (see eclipse).
02787              */
02788             if ( -1 == sinfo_function1d_natural_spline(xnum,row_data, ilx, 
02789                                                        eval,corrected_row_data,
02790                                                        ilx ) )
02791         {
02792             sinfo_msg_error("error in spline interpolation\n") ;
02793         cpl_imagelist_delete(returnCube) ;
02794         return NULL ;
02795         }   
02796 
02797             new_sum = 0. ;
02798             for ( col = 0 ; col < ilx ; col++ )
02799             {
02800                 if (isnan(corrected_row_data[col])) continue ;
02801                 /* don't take the sinfo_edge points to calculate 
02802                    the scaling factor */
02803                 if ( col != 0 && col != ilx - 1 )
02804                 {
02805                     new_sum += corrected_row_data[col] ;
02806                 }
02807             }
02808             for ( col = 0 ; col < ilx ; col++ )
02809             {
02810                 /* ----------------------------------------------------------
02811                  * rescale the row data and fill the returned cube, 
02812                  * leave the left-most and right-most
02813                  * pixel column 
02814                  */
02815                 if ( col == 0 )
02816                 {
02817                     podata[col+row*ilx] = ZERO ;
02818                 }
02819                 else if ( col == ilx - 1 )
02820                 {
02821                     podata[col+row*ilx] = ZERO ;
02822                 }
02823                 else
02824                 {
02825                     if ( isnan(corrected_row_data[col]) ) 
02826                     {
02827                         podata[col+row*ilx] = ZERO ;
02828                     }
02829                     else
02830                     {
02831                         if (new_sum == 0.) new_sum = 1. ;
02832                      /* rescaling is commented out because it delivers 
02833                         wrong results
02834                         in case of appearance of blanks or bad pixels */
02835                   /*       corrected_row_data[col] *= sum / new_sum ; */
02836                         podata[col+row*ilx] = corrected_row_data[col] ;
02837                     }
02838                 }
02839             }
02840         }
02841     }       
02842 
02843     cpl_free(row_data) ;
02844     cpl_free(corrected_row_data) ;
02845     cpl_free(xnum) ;
02846     cpl_free(eval) ;
02847 
02848     return returnCube ;
02849 }
02850 
02868 float * 
02869 sinfo_new_calibrate_ns_test( cpl_image * ns_image,
02870                            int        n_slitlets,
02871                            int        halfWidth,
02872                            float      fwhm,
02873                            float      minDiff,
02874                            float      estimated_dist,
02875                            float      devtol,
02876                int        bottom,
02877                int        top )
02878 {
02879     int i, j, k, m, row, col, n, ni, na ;
02880     int position, counter, iters ;
02881     int xdim, ndat, its, numpar ;
02882     float sum, mean, maxval ;
02883     float tol, lab ;
02884     float * distances ;
02885     float * ret_distances ;
02886 
02887     float * xdat, * wdat ;
02888     int * mpar ;
02889 
02890     pixelvalue* row_buf=NULL ;
02891     float** distances_buf=NULL ;
02892     float* x_position=NULL ;
02893     int* found=NULL;
02894     int* found_clean=NULL ;
02895     int* found_cleanit=NULL ;
02896 
02897     Vector * line ;
02898     FitParams ** par ;
02899     int foundit, begin, end ;
02900     int zeroindicator ;
02901     int row_index ;
02902 
02903     int ilx=0;
02904     int ily=0;
02905     float* pidata=NULL;
02906 
02907     if ( ns_image == NULL )
02908     {
02909         sinfo_msg_error("sorry, no image given\n") ;
02910         return NULL ;
02911     }
02912     if ( n_slitlets < 1 )
02913     {
02914         sinfo_msg_error("wrong number of slitlets given\n") ;
02915         return NULL ;
02916     }
02917     if ( halfWidth < 0 || halfWidth >= estimated_dist )
02918     {
02919         sinfo_msg_error("wrong half width given\n") ;
02920         return NULL ;
02921     }
02922     if ( fwhm <= 0. )
02923     {
02924         sinfo_msg_error("wrong fwhm given\n") ;
02925         return NULL ;
02926     }
02927     if ( minDiff < 1. )
02928     {
02929         sinfo_msg_error("wrong minDiff given\n") ;
02930         return NULL ;
02931     }
02932 
02933     /* allocate memory for output array */
02934     if (NULL==(distances=(float *)cpl_calloc( n_slitlets , sizeof (float) ))) 
02935     {
02936         sinfo_msg_error("could not allocate memory\n") ;
02937         return NULL ;
02938     }
02939     /* allocate memory for output array */
02940     if (NULL == (ret_distances = (float *) cpl_calloc ( n_slitlets , 
02941                                            sizeof (float) ))) 
02942     {
02943         sinfo_msg_error("could not allocate memory\n") ;
02944         return NULL ;
02945     }
02946 
02947     ilx=cpl_image_get_size_x(ns_image);
02948     ily=cpl_image_get_size_y(ns_image);
02949     pidata=cpl_image_get_data_float(ns_image);
02950 
02951     row_buf=(pixelvalue*)cpl_calloc(ilx,sizeof(pixelvalue)) ;
02952     x_position=cpl_calloc(n_slitlets,sizeof(float)) ;
02953     found=cpl_calloc(3*n_slitlets,sizeof(int));
02954     found_clean=cpl_calloc(3*n_slitlets,sizeof(int)) ;
02955     found_cleanit=cpl_calloc(3*n_slitlets,sizeof(int)) ;
02956     distances_buf=sinfo_new_2Dfloatarray(ily,n_slitlets) ;
02957 
02958     /* go through the image rows */
02959     for ( row = 0 ; row < ily ; row++ )
02960     {
02961         zeroindicator = 0 ;
02962 
02963         /* initialize the distance buffer */
02964         for ( i = 0 ; i < n_slitlets ; i++ )
02965         {
02966             distances_buf[row][i] = ZERO ;
02967         }
02968 
02969         /* fill the row buffer array with image data */
02970         for ( col = 0 ; col < ilx ; col++ )
02971         {
02972             row_buf[col] = pidata[col + row*ilx] ;
02973         }
02974 
02975         /* determine the mean of the row data */
02976         sum = 0. ;
02977         n = 0 ;
02978         for ( i = 0 ; i < ilx ; i++ )
02979         {
02980             if ( isnan(row_buf[i]) )
02981             {
02982                 continue ;
02983             }
02984             sum += row_buf[i] ;
02985             n++ ;
02986         }
02987         mean = sum / (float)n ;
02988 
02989         /* store the positions of image values greater than the mean */
02990         n = 0 ;
02991         for ( i = 0 ; i < ilx ; i++ )
02992         {
02993             if (isnan(row_buf[i]))
02994             {
02995                 continue ;
02996             }
02997             if ( row_buf[i] > mean + ESTIMATE )
02998             {
02999                 found[n] = i ;
03000                 n++ ;
03001             }
03002         }
03003        
03004         if ( n < n_slitlets )
03005         {
03006             sinfo_msg_warning("t4 wrong number of intensity columns "
03007                               "found in row: %d, found number: %d", row, n) ;
03008             continue ;
03009         }
03010         else
03011         { 
03012             /* find the maximum value position around the found columns */
03013             na = 0 ;
03014             for ( i = 1 ; i < n ; i ++ )
03015             {
03016                 if ( found[i] - found[i-1] < halfWidth )
03017                 {
03018                     begin = found[i] - halfWidth ;
03019                     if ( begin < 0 )
03020                     {
03021                         begin = 0 ;
03022                     }
03023                     end = found[i] + halfWidth ;
03024                     if ( end >= ilx )
03025                     {
03026                         end = ilx - 1 ;
03027                     }
03028                     /* find the maximum value inside the box around 
03029                        the found positions*/
03030                     maxval = -FLT_MAX ;
03031                     foundit = 0 ;
03032                     for ( j = begin ; j <= end ; j++ )
03033                     {
03034                         /* do not consider boxes that contain bad pixels */
03035                         if (isnan(row_buf[j]))
03036                         {
03037                             continue ;
03038                         }
03039                         if (row_buf[j] >= maxval )
03040                         {
03041                             maxval = row_buf[j] ;
03042                             foundit = j ;
03043                         }
03044                     }
03045                     if (maxval == -FLT_MAX)
03046                     {
03047                         continue ;
03048                     }
03049                     for ( k = 0 ; k < na ; k++ )
03050                     {
03051                         if ( found_cleanit[k] >= begin && 
03052                              found_cleanit[k] < foundit )
03053                         {
03054                             na-- ;
03055                         }
03056                     }
03057                     for ( k = 0 ; k < n ; k++ )
03058                     {
03059                         if ( found[k] == foundit)
03060                         {
03061                             if ( found_cleanit[na-1] != found[k] )
03062                             {
03063                                 found_cleanit[na] = found[k] ;
03064                                 na++ ;
03065                             }
03066                         }
03067                     }
03068                 }
03069                 else
03070                 {
03071                     if ( i == 1 )
03072                     {
03073                         found_cleanit[na] = found[0] ;
03074                         na++ ;
03075                         found_cleanit[na] = found[1] ;
03076                         na++ ;
03077                     }
03078                     else
03079                     {
03080                         if ( found_cleanit[na-1] != found[i-1])
03081                         {
03082                             found_cleanit[na] = found[i-1] ;
03083                             na++ ;
03084                         }
03085                         if ( found_cleanit[na-1] != found[i])
03086                         {
03087                             found_cleanit[na] = found[i] ;
03088                             na++ ;
03089                         }
03090                     }
03091                 }
03092             }
03093 
03094             /* determine only one pixel position for each slitlet intensity */
03095             j = 1 ;
03096             for ( i = 1 ; i < na ; i++ )
03097             {
03098                 if ( (float)(found_cleanit[i] - found_cleanit[i-1]) < 
03099                              (estimated_dist - devtol) ||
03100                      (float)(found_cleanit[i] - found_cleanit[i-1]) > 
03101                              (estimated_dist + devtol) )
03102                 {
03103                     continue ;
03104                 }
03105                 else
03106                 {
03107                     found_clean[j-1] = found_cleanit[i-1] ;
03108                     found_clean[j]   = found_cleanit[i] ;
03109                     j++ ;
03110                 }
03111             }
03112         }
03113         if ( j > n_slitlets )
03114         {
03115             /* check the distance again */
03116             ni = 1 ;
03117             for ( i = 1 ; i < j ; i++ )
03118             {
03119                 if ( (float)(found_clean[i] - found_clean[i-1]) < 
03120                             (estimated_dist - devtol ) ||
03121                      (float)(found_clean[i] - found_clean[i-1]) > 
03122                             (estimated_dist + devtol ) )
03123                 { 
03124                     continue ;
03125                 }
03126                 else
03127                 {
03128                     found_clean[ni-1] = found_clean[i-1] ;
03129                     found_clean[ni]   = found_clean[i] ;
03130                     ni++ ;
03131                 }
03132             }
03133             if ( ni != n_slitlets )
03134             {
03135                 sinfo_msg_warning("t5 wrong number of intensity columns "
03136                                   "found in row: %d, found number: %d",
03137                                   row,ni) ;
03138                 continue ;
03139             }
03140             else 
03141             {
03142                 j = ni ;
03143             }
03144         }
03145         else if ( j < n_slitlets )
03146         {
03147             sinfo_msg_warning("t6 wrong number of intensity columns found "
03148                               "in row: %d , found number: %d\n", row, j) ;
03149             continue ;
03150         }
03151         counter = 0 ;
03152         /* go through the found intensity pixels in one row */
03153         for ( i = 0 ; i < j ; i++ )
03154         {
03155             /* allocate memory for the array where the line is fitted in */
03156             if ( NULL == (line = sinfo_new_vector (2*halfWidth + 1)) )
03157             {
03158                 sinfo_msg_error ("cannot allocate new Vector \n") ;
03159                 cpl_free(distances) ;
03160                 return NULL ;
03161             }
03162 
03163             /* allocate memory */
03164             xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
03165             wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
03166             mpar = (int *)   cpl_calloc( MAXPAR, sizeof (int) ) ;
03167             par = sinfo_new_fit_params(1) ;
03168 
03169             m = 0 ;
03170             for ( k = found_clean[i]-halfWidth ; 
03171                   k <= found_clean[i]+halfWidth ; k++ )
03172             {
03173                 if ( k < 0 )
03174                 {
03175                     k = 0. ;
03176                 }
03177                 else if ( k >= ilx )
03178                 {
03179                     k = ilx - 1 ;
03180                 }
03181                 else if ( isnan(row_buf[k]) )
03182                 {
03183                     zeroindicator = 1 ;
03184                     break ;
03185                 }
03186                 else
03187                 {
03188                     line -> data[m] = row_buf[k] ;
03189                     m++ ;
03190                 }
03191             }
03192             if ( zeroindicator == 1 )
03193             {
03194                 sinfo_new_destroy_vector(line) ;
03195                 cpl_free(xdat) ;
03196                 cpl_free(wdat) ;
03197                 cpl_free(mpar) ;
03198                 sinfo_new_destroy_fit_params(&par) ;
03199                 break ;
03200             }
03201 
03202             /*-----------------------------------------------------------------
03203              * go through the spectral sinfo_vector
03204              * determine the maximum pixel value in the spectral sinfo_vector
03205              */
03206             maxval = -FLT_MAX ;
03207             position = -INT32_MAX ;
03208             for ( k = 0 ; k < m ; k++ )
03209             {
03210                 xdat[k] = k ;
03211                 wdat[k] = 1.0 ;
03212                 if ( line -> data[k] >= maxval )
03213                 {
03214                     maxval = line -> data[k] ;
03215                     position = k ;
03216                 }
03217             }
03218 
03219             /* set initial values for the fitting routine */
03220             xdim     = XDIM ;
03221             ndat     = line -> n_elements ;
03222             numpar   = MAXPAR ;
03223             tol      = TOL ;
03224             lab      = LAB ;
03225             its      = ITS ;
03226             (*par) -> fit_par[1] = fwhm ;
03227             (*par) -> fit_par[2] = (float) position ;
03228             (*par) -> fit_par[3] = (float) (line -> data[0] + 
03229                                     line -> data[line->n_elements - 1]) / 2.0 ;
03230             (*par) -> fit_par[0]  = maxval - ((*par) -> fit_par[3]) ;
03231 
03232 
03233             /* exclude negative peaks and low signal cases */
03234             if ( (*par) -> fit_par[0] < minDiff )
03235             {
03236                 sinfo_msg_warning ("sorry, signal of line too low to fit "
03237                                    "in row: %d in slitlet %d\n", row, i) ;
03238                 sinfo_new_destroy_vector(line) ;
03239                 cpl_free(xdat) ;
03240                 cpl_free(wdat) ;
03241                 cpl_free(mpar) ;
03242                 sinfo_new_destroy_fit_params(&par) ;
03243                 continue ;
03244             }
03245 
03246             for ( k = 0 ; k < MAXPAR ; k++ )
03247             {
03248                 (*par) -> derv_par[k] = 0.0 ;
03249                 mpar[k] = 1 ;
03250             }
03251             /* finally, do the least square fit using a sinfo_gaussian */
03252             if ( 0 > ( iters = sinfo_new_lsqfit_c(xdat, &xdim, 
03253                                                   line -> data, wdat, 
03254                                                   &ndat, (*par) -> fit_par,
03255                                                   (*par) -> derv_par, mpar, 
03256                                                   &numpar, &tol, 
03257                                                   &its, &lab )) )
03258             {
03259           /*
03260                cpl_msg_debug ("sinfo_calibrate_ns_test:",
03261                               "sinfo_lsqfit_c: least squares fit failed,"
03262                               " error no.: %d in row: %d in slitlet %d\n", 
03263                               iters, row, i) ;
03264           */
03265                 sinfo_new_destroy_vector(line) ;
03266                 cpl_free(xdat) ;
03267                 cpl_free(wdat) ;
03268                 cpl_free(mpar) ;
03269                 sinfo_new_destroy_fit_params(&par) ;
03270                 continue ;
03271             }
03272 
03273             /* check for negative fit results */
03274             if ( (*par) -> fit_par[0] <= 0. || (*par) -> fit_par[1] <= 0. ||
03275                  (*par) -> fit_par[2] < 0. )
03276             {
03277                 sinfo_msg_warning ("negative parameters as fit result, not "
03278                                    "used! in row %d in slitlet %d", row, i) ;
03279                 sinfo_new_destroy_vector(line) ;
03280                 cpl_free(xdat) ;
03281                 cpl_free(wdat) ;
03282                 cpl_free(mpar) ;
03283                 sinfo_new_destroy_fit_params(&par) ;
03284                 continue ;
03285             }
03286 
03287             /* correct the fitted position for the given row of the line 
03288                in image coordinates */
03289             (*par) -> fit_par[2] =  (float) (found_clean[i] - halfWidth) + 
03290                                     (*par) -> fit_par[2] ;
03291             x_position[counter] = (*par) -> fit_par[2] ;
03292             counter ++ ;
03293 
03294             /* free memory */
03295             sinfo_new_destroy_fit_params(&par) ;
03296             sinfo_new_destroy_vector ( line ) ;
03297             cpl_free ( xdat ) ;
03298             cpl_free ( wdat ) ;
03299             cpl_free ( mpar ) ;
03300         }
03301         if (zeroindicator == 1)
03302         {
03303             sinfo_msg_debug ("bad pixel in fitting box in row: %d\n", row) ;
03304             continue ;
03305         }
03306 
03307         if ( counter != n_slitlets )
03308         {
03309             sinfo_msg_warning("wrong number of slitlets found "
03310                               "in row: %d", row) ;
03311             continue ;
03312         }
03313         /* store the distances between the sources and the slitlet centers */
03314         for ( i = 0 ; i < n_slitlets ; i++ )
03315         {
03316             distances_buf[row][i] = x_position[i] - (15.5 + 32.*(float)i) ;
03317         }
03318     }
03319 
03320     /* ----------------------------------------------------------------
03321      * go through the rows again and take the mean of the distances, 
03322      * throw away the runaways 
03323      */
03324     for ( i = 0 ; i < n_slitlets ; i++ )
03325     {
03326         n   = 0 ;
03327         sum = 0. ;
03328         for ( row = bottom ; row < top ; row++ )
03329         {
03330             if ( fabs( distances_buf[row][i] ) > devtol || 
03331                  isnan(distances_buf[row][i]) )
03332             {
03333              /*
03334           sinfo_msg("dist=%g devtol=%g isan=%d", 
03335             distances_buf[row][i],
03336             devtol,
03337             isnan(distances_buf[row][i]));
03338               */              
03339                 continue ;
03340             }
03341             sum += distances_buf[row][i] ;
03342             n++ ;
03343         }
03344         if ( n < 2 )
03345         {
03346             sinfo_msg_error("distances array could not be determined"
03347                             " completely!, deviations of distances from"
03348                             " devtol too big" ) ;
03349             cpl_free(distances) ;
03350             return NULL ;
03351         }
03352         else
03353         {
03354             distances[i] = sum / (float)n ;
03355         }
03356     }
03357 
03358     /* now sort the result according to the row sequence in the 
03359        reconstructed image*/
03360     for ( i = 0 ; i < n_slitlets ; i++ )
03361     {
03362         switch (i)
03363         {
03364             case 0:
03365                 row_index = 8 ;
03366                 break ;
03367             case 1:
03368                 row_index = 7 ;
03369                 break ;
03370             case 2:
03371                 row_index = 9 ;
03372                 break ;
03373             case 3:
03374                 row_index = 6 ;
03375                 break ;
03376             case 4:
03377                 row_index = 10 ;
03378                 break ;
03379             case 5:
03380                 row_index = 5 ;
03381                 break ;
03382             case 6:
03383                 row_index = 11 ;
03384                 break ;
03385             case 7:
03386                 row_index = 4 ;
03387                 break ;
03388             case 8:
03389                 row_index = 12 ;
03390                 break ;
03391             case 9:
03392                 row_index = 3 ;
03393                 break ;
03394             case 10:
03395                 row_index = 13 ;
03396                 break ;
03397             case 11:
03398                 row_index = 2 ;
03399                 break ;
03400             case 12:
03401                 row_index = 14 ;
03402                 break ;
03403             case 13:
03404                 row_index = 1 ;
03405                 break ;
03406             case 14:
03407                 row_index = 15 ;
03408                 break ;
03409             case 15:
03410                 row_index = 0 ;
03411                 break ;
03412             case 16:
03413                 row_index = 31 ;
03414                 break ;
03415             case 17:
03416                 row_index = 16 ;
03417                 break ;
03418             case 18:
03419                 row_index = 30 ;
03420                 break ;
03421             case 19:
03422                 row_index = 17 ;
03423                 break ;
03424             case 20:
03425                 row_index = 29 ;
03426                 break ;
03427             case 21:
03428                 row_index = 18 ;
03429                 break ;
03430             case 22:
03431                 row_index = 28 ;
03432                 break ;
03433             case 23:
03434                 row_index = 19 ;
03435                 break ;
03436             case 24:
03437                 row_index = 27 ;
03438                 break ;
03439             case 25:
03440                 row_index = 20 ;
03441                 break ;
03442             case 26:
03443                 row_index = 26 ;
03444                 break ;
03445             case 27:
03446                 row_index = 21 ;
03447                 break ;
03448             case 28:
03449                 row_index = 25 ;
03450                 break ;
03451             case 29:
03452                 row_index = 22 ;
03453                 break ;
03454             case 30:
03455                 row_index = 24 ;
03456                 break ;
03457             case 31:
03458                 row_index = 23 ;
03459                 break ;
03460             default:
03461                 sinfo_msg_error("wrong number of a slitlet\n") ;
03462                 cpl_free (distances) ;
03463                 return NULL ;
03464         }
03465         ret_distances[row_index] = distances[i] ;
03466     }
03467     cpl_free(distances) ;
03468 
03469     cpl_free(row_buf) ;
03470     cpl_free(x_position) ;
03471     cpl_free(found);
03472     cpl_free(found_clean) ;
03473     cpl_free(found_cleanit) ;
03474     sinfo_new_destroy_2Dfloatarray(&distances_buf,n_slitlets) ;
03475 
03476 
03477     return ret_distances ; 
03478 }
03503 cpl_image * 
03504 sinfo_new_make_true_resamp(cpl_image * calibImage, 
03505                            cpl_image * wavemap)
03506 {
03507     cpl_image * returnImage ;
03508     float edges[33] ;
03509     int imsize, kslit,i,j ;
03510     int slit_index ;
03511     int z, col, recol ;
03512     int wlx=0;
03513     int wly=0;
03514     int clx=0;
03515     int cly=0;
03516   
03517     float* pcdata=NULL;
03518     float* pwdata=NULL;
03519     float* prdata=NULL;
03520 
03521 
03522     wlx=cpl_image_get_size_x(wavemap);
03523     wly=cpl_image_get_size_y(wavemap);
03524     pwdata=cpl_image_get_data_float(wavemap);
03525 
03526     edges[0]=0;
03527     j=1;
03528     for(i=0;i<wlx-1;i++)
03529     {
03530         if((pwdata[i]-pwdata[i+1])>0.0025 || (pwdata[i]-pwdata[i+1])<-0.0025)
03531     {
03532         sinfo_msg_error("wavemap sinfo_edge %d", i+1);
03533         edges[j]=i+1;
03534         j++;
03535     }
03536     }
03537     edges[32]=2048;
03538 
03539     clx=cpl_image_get_size_x(calibImage);
03540     cly=cpl_image_get_size_y(calibImage);
03541     pcdata=cpl_image_get_data_float(calibImage);
03542  
03543     imsize = clx / N_SLITLETS ;
03544 
03545     /* allocate memory */  
03546     returnImage = cpl_image_new(clx,cly,CPL_TYPE_FLOAT);
03547     prdata=cpl_image_get_data_float(returnImage);
03548     for ( z = 0 ; z < cly ; z++ ) /* go through the z-axis */
03549     {
03550         for ( col = 0 ; col < clx ; col++ ) /* go through the image columns */
03551         prdata[col+z*clx]=ZERO;
03552     }
03553     
03554 
03555     /* now build the data cube out of the resampled image */
03556     for ( z = 0 ; z < cly ; z++ ) /* go through the z-axis */
03557     {
03558         kslit      = 0 ;
03559         slit_index = -1 ;
03560         recol      = -1 ;
03561         for ( col = 0 ; col < clx ; col++ ) /* go through the image columns */
03562         {
03563             /*if ( col % imsize == 0 )
03564             {*/
03565                 recol = 0 ;
03566                 /*kslit = col/imsize ;*/
03567         for(i=0;i<32;i++)
03568         {
03569           if(col>=sinfo_new_nint(edges[i]) && 
03570                      col<sinfo_new_nint(edges[i+1]))
03571             kslit=i;
03572         }
03573                 /* sort the slitlets in the right spiffi specific way */
03574                 switch (kslit)
03575                 {
03576                     case 0:
03577                         slit_index = 8 ;
03578                         break ;
03579                     case 1:
03580                         slit_index = 7 ;
03581                         break ;
03582                     case 2:
03583                         slit_index = 9 ;
03584                         break ;
03585                     case 3:
03586                         slit_index = 6 ;
03587                         break ;
03588                     case 4:
03589                         slit_index = 10 ;
03590                         break ;
03591                     case 5:
03592                         slit_index = 5 ;
03593                         break ;
03594                     case 6:
03595                         slit_index = 11 ;
03596                         break ;
03597                     case 7:
03598                         slit_index = 4 ;
03599                         break ;
03600                     case 8:
03601                         slit_index = 12 ;
03602                         break ;
03603                     case 9:
03604                         slit_index = 3 ;
03605                         break ;
03606                     case 10:
03607                         slit_index = 13 ;
03608                         break ;
03609                     case 11:
03610                         slit_index = 2 ;
03611                         break ;
03612                     case 12:
03613                         slit_index = 14 ;
03614                         break ;
03615                     case 13:
03616                         slit_index = 1 ;
03617                         break ;
03618                     case 14:
03619                         slit_index = 15 ;
03620                         break ;
03621                     case 15:
03622                         slit_index = 0 ;
03623                         break ;
03624                     case 16:
03625                         slit_index = 31 ;
03626                         break ;
03627                     case 17:
03628                         slit_index = 16 ;
03629                         break ;
03630                     case 18:
03631                         slit_index = 30 ;
03632                         break ;
03633                     case 19:
03634                         slit_index = 17 ;
03635                         break ;
03636                     case 20:
03637                         slit_index = 29 ;
03638                         break ;
03639                     case 21:
03640                         slit_index = 18 ;
03641                         break ;
03642                     case 22:
03643                         slit_index = 28 ;
03644                         break ;
03645                     case 23:
03646                         slit_index = 19 ;
03647                         break ;
03648                     case 24:
03649                         slit_index = 27 ;
03650                         break ;
03651                     case 25:
03652                         slit_index = 20 ;
03653                         break ;
03654                     case 26:
03655                         slit_index = 26 ;
03656                         break ;
03657                     case 27:
03658                         slit_index = 21 ;
03659                         break ;
03660                     case 28:
03661                         slit_index = 25 ;
03662                         break ;
03663                     case 29:
03664                         slit_index = 22 ;
03665                         break ;
03666                     case 30:
03667                         slit_index = 24 ;
03668                         break ;
03669                     case 31:
03670                         slit_index = 23 ;
03671                         break ;
03672                     default:
03673                         sinfo_msg_error("wrong slitlet index: couldn't be a "
03674                                         "spiffi image,  there must be 32 "
03675                                         "slitlets!") ;
03676                         break ;
03677                 }
03678 
03679             /*}*/
03680 
03681             /* fill each cube plane with one image row */
03682         if((col-sinfo_new_nint(edges[kslit]))>0 && 
03683                (col-sinfo_new_nint(edges[kslit]))<imsize-1 )
03684                 prdata[(col-sinfo_new_nint(edges[kslit]))+
03685                        slit_index*imsize+z*clx] = 
03686                       pcdata[col+z*clx] ;
03687             else
03688             prdata[(col-sinfo_new_nint(edges[kslit]))+
03689                         slit_index*imsize+z*clx] = ZERO;
03690             /*recol++ ;*/
03691 
03692         }
03693     }
03694     return returnImage ;
03695 }
03696 
03697 /*The old slitlet order*/
03698 /*switch (kslit)
03699                 {
03700                     case 0:
03701                         slit_index = 23 ;
03702                         break ;
03703                     case 1:
03704                         slit_index = 24 ;
03705                         break ;
03706                     case 2:
03707                         slit_index = 22 ;
03708                         break ;
03709                     case 3:
03710                         slit_index = 25 ;
03711                         break ;
03712                     case 4:
03713                         slit_index = 21 ;
03714                         break ;
03715                     case 5:
03716                         slit_index = 26 ;
03717                         break ;
03718                     case 6:
03719                         slit_index = 20 ;
03720                         break ;
03721                     case 7:
03722                         slit_index = 27 ;
03723                         break ;
03724                     case 8:
03725                         slit_index = 19 ;
03726                         break ;
03727                     case 9:
03728                         slit_index = 28 ;
03729                         break ;
03730                     case 10:
03731                         slit_index = 18 ;
03732                         break ;
03733                     case 11:
03734                         slit_index = 29 ;
03735                         break ;
03736                     case 12:
03737                         slit_index = 17 ;
03738                         break ;
03739                     case 13:
03740                         slit_index = 30 ;
03741                         break ;
03742                     case 14:
03743                         slit_index = 16 ;
03744                         break ;
03745                     case 15:
03746                         slit_index = 31 ;
03747                         break ;
03748                     case 16:
03749                         slit_index = 0 ;
03750                         break ;
03751                     case 17:
03752                         slit_index = 15 ;
03753                         break ;
03754                     case 18:
03755                         slit_index = 1 ;
03756                         break ;
03757                     case 19:
03758                         slit_index = 14 ;
03759                         break ;
03760                     case 20:
03761                         slit_index = 2 ;
03762                         break ;
03763                     case 21:
03764                         slit_index = 13 ;
03765                         break ;
03766                     case 22:
03767                         slit_index = 3 ;
03768                         break ;
03769                     case 23:
03770                         slit_index = 12 ;
03771                         break ;
03772                     case 24:
03773                         slit_index = 4 ;
03774                         break ;
03775                     case 25:
03776                         slit_index = 11 ;
03777                         break ;
03778                     case 26:
03779                         slit_index = 5 ;
03780                         break ;
03781                     case 27:
03782                         slit_index = 10 ;
03783                         break ;
03784                     case 28:
03785                         slit_index = 6 ;
03786                         break ;
03787                     case 29:
03788                         slit_index = 9 ;
03789                         break ;
03790                     case 30:
03791                         slit_index = 7 ;
03792                         break ;
03793                     case 31:
03794                         slit_index = 8 ;
03795                         break ;
03796                     default:
03797                         sinfo_msg_error("wrong slitlet index: couldn't "
03798                                         "be a spiffi image,  \
03799                                  there must be 32 slitlets!\n") ;
03800                         cpl_imagelist_delete(returnCube) ;
03801                         return NULL ;
03802                         break ;
03803                 }*/
03804 
03805 
03806 /*--------------------------------------------------------------------------*/

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