cube_construct.c

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

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