sinfo_absolute.c

00001 /*
00002  * This file is part of the ESO SINFONI Pipeline
00003  * Copyright (C) 2004,2005 European Southern Observatory
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or
00008  * (at your option) any later version.
00009  *
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00018  */
00019 /*****************************************************************************
00020 * E.S.O. - VLT project
00021 *
00022 *
00023 *
00024 * who       when      what
00025 * --------  --------  ----------------------------------------------
00026 * schreib  14/11/00  created
00027 */
00028 /*---------------------------------------------------------------------------*/
00033 /*---------------------------------------------------------------------------*/
00037 /************************************************************************
00038 *   NAME
00039 *        sinfo_absolute.c - routines to determine the absolute positions 
00040 *        of the slitlets out of an emission line frame
00041 *
00042 *   SYNOPSIS
00043 *   #include "absolute.h"
00044 *
00045 *   1) float sinfo_new_edge( float * xdat, float * parlist, 
00046                              int * npar, int * ndat )
00047 *   2) void sinfo_new_edge_deriv( float * xdat, float * parlist, 
00048                                   float * dervs, int * npar )
00049 *   3) static int sinfo_new_inv_mat_edge (void)
00050 *   4) static new_void sinfo_get_mat_edge ( float * xdat,
00051 *                                 int   * xdim,
00052 *                                 float * ydat,
00053 *                                 float * wdat,
00054 *                                 int   * ndat,
00055 *                                 float * fpar,
00056 *                                 float * epar,
00057 *                                 int   * npar )
00058 *   5) static int sinfo_new_get_vec_edge ( float * xdat,
00059 *                                int   * xdim,
00060 *                                float * ydat,
00061 *                                float * wdat,
00062 *                                int   * ndat,
00063 *                                float * fpar,
00064 *                                float * epar,
00065 *                                int   * npar )
00066 *   6) int sinfo_new_lsqfit_edge ( float * xdat,
00067 *                        int   * xdim,
00068 *                        float * ydat,
00069 *                        float * wdat,
00070 *                        int   * ndat,
00071 *                        float * fpar,
00072 *                        float * epar,
00073 *                        int   * mpar,
00074 *                        int   * npar,
00075 *                        float * tol ,
00076 *                        int   * its ,
00077 *                        float * lab  )
00078 *   7) int sinfo_new_fit_slits_edge( cpl_image   * lineImage, 
00079 *                        FitParams ** par,
00080 *                        float     ** sinfo_slit_pos,
00081 *                        int          box_length,
00082 *                        float        y_box,
00083 *                        float        diff_tol )
00084 *   8) int sinfo_new_fit_slits_edge_with_estimate ( cpl_image   * lineImage,
00085 *                                     float     ** sinfo_slit_pos,
00086 *                                     int          box_length,
00087 *                                     float        y_box,
00088 *                                     float        diff_tol,
00089 *                                     int          low_pos,
00090 *                                     int          high_pos )
00091 *
00092 *   DESCRIPTION
00093 *   1) calculates the value of a slope function with parameters 
00094 *      parlist at the position xdat 
00095 *   2) calculates the partial derivatives for a slope function with
00096 *      parameters parlist at position xdat 
00097 *   3) calculates the inverse of matrix2. The algorithm used 
00098 *      is the Gauss-Jordan algorithm described in Stoer,
00099 *      Numerische Mathematik, 1. Teil.
00100 *   4) builds the sinfo_matrix 
00101 *   5) calculates the correction sinfo_vector. The sinfo_matrix has been
00102 *      built by get_mat_edge(), we only have to rescale it for the 
00103 *      current value of labda. The sinfo_matrix is rescaled so that
00104 *      the diagonal gets the value 1 + labda.
00105 *      Next we calculate the inverse of the sinfo_matrix and then
00106 *      the correction sinfo_vector.
00107 *   6) this is a routine for making a least-squares fit of a
00108 *      function to a set of data points. The method used is
00109 *      described in: Marquardt, J.Soc.Ind.Appl.Math. 11. 431 (1963).
00110 *      This method is a mixture of the steepest descent method 
00111 *      and the Taylor method.
00112 *   7) fits the beginning and end position of the slitlets
00113 *      by using non-linear least square fitting of a step function
00114 *      fits a step function to the slitlet edges exposed and indicated
00115 *      by the brightest emission lines. To achieve this, the fit
00116 *      parameters are used to find the brightest emission line
00117 *      and to get its position for each column.
00118 *      The least squares fit is done by using a box smaller than
00119 *      the size of two slitlets
00120 *   8) fits the beginning and end position of the slitlets
00121 *      by using non-linear least square fitting of an sinfo_edge  function
00122 *      fits a linear edge function to the slitlet edges exposed and indicated
00123 *      by the brightest emission lines. The slitlet is searched within
00124 *      user given positions.
00125 *      The least squares fit is done by using a box smaller than
00126 *      the size of two slitlets 
00127 *
00128 *   FILES
00129 *
00130 *   ENVIRONMENT
00131 *
00132 *   RETURN VALUES
00133 *
00134 *   CAUTIONS
00135 *
00136 *   EXAMPLES
00137 *
00138 *   SEE ALSO
00139 *
00140 *   BUGS
00141 *
00142 *------------------------------------------------------------------------
00143 */
00144 #ifdef HAVE_CONFIG_H
00145 #  include <config.h>
00146 #endif
00147 #include "sinfo_vltPort.h"
00148 /*
00149  * System Headers
00150  */
00151 
00152 /*
00153  * Local Headers
00154  */
00155 
00156 #include "sinfo_absolute.h"
00157 #include "sinfo_recipes.h"
00158 /*----------------------------------------------------------------------------
00159  *                              Defines
00160  *--------------------------------------------------------------------------*/
00161 static float  sqrarg ;
00162 #define SQR(a) (sqrarg = (a) , sqrarg*sqrarg)
00163 
00164 #define XDIMA         1         /* dimension of the x values */
00165 #define TOLA          0.001     /* fitting tolerance */
00166 #define LABA          0.1       /* labda parameter */
00167 #define ITSA          200       /* maximum number of iterations */
00168 #define LABFACA       10.0      /* labda step factor */
00169 #define LABMAXA       1.0e+10   /* maximum value for labda */
00170 #define LABMINA       1.0e-10   /* minimum value for labda */
00171 #define NPAR          4         /* number of fit parameters */
00172 
00173 /*----------------------------------------------------------------------------
00174  *                                 Local variables
00175  *--------------------------------------------------------------------------*/
00176 
00177 static double chi1 ;                    /* old reduced chi-squared */
00178 static double chi2 ;                    /* new reduced chi-squared */
00179 static double labda ;                   /* mixing parameter */
00180 static double vec[NPAR] ;               /* correction sinfo_vector */
00181 static double matrix1[NPAR][NPAR] ;     /* original sinfo_matrix */
00182 static double matrix2[NPAR][NPAR] ;     /* inverse of matrix1 */
00183 static int    nfree ;                   /* number of free parameters */
00184 static int    parptr[NPAR] ;            /* parameter pointer */
00185 static float  slopewidth ;              /* initial value for fit parameter 5: 
00186                                            width of linear slope */
00187 
00188 /*----------------------------------------------------------------------------
00189  *                  Functions private to this module
00190  *--------------------------------------------------------------------------*/
00191 static int sinfo_new_inv_mat_edge (void) ;
00192 
00193 static void sinfo_new_get_mat_edge ( float * xdat,
00194                            int   * xdim,
00195                            float * ydat,
00196                            float * wdat,
00197                            int   * ndat,
00198                            float * fpar,
00199                            float * epar/*,
00200                            int   * npar*/ ) ;
00201 
00202 static int sinfo_new_get_vec_edge ( float * xdat,
00203                           int   * xdim,
00204                           float * ydat,
00205                           float * wdat,
00206                           int   * ndat,
00207                           float * fpar,
00208                           float * epar,
00209                           int   * npar ) ;
00210 float 
00211 sinfo_new_hat2 ( float * xdat, float * parlist/*, int * npar, int * ndat*/ );
00212 
00213 float 
00214 sinfo_new_hat1 ( float * xdat, float * parlist/*, int * npar, int * ndat*/ );
00215 
00216 void 
00217 sinfo_new_hat_deriv2(float * xdat, float * parlist, 
00218                      float * dervs/*, int * npar*/ );
00219 
00220 void 
00221 sinfo_new_hat_deriv1( float * xdat, float * parlist, 
00222                       float * dervs/*, int * npar*/ );
00223 
00224 int 
00225 sinfo_new_fit_slits1( cpl_image   * lineImage, 
00226                FitParams ** par,
00227                float     ** sinfo_slit_pos,
00228                int          box_length,
00229               float        y_box );
00230 
00231 int 
00232 sinfo_new_fit_slits( cpl_image   * lineImage, 
00233               FitParams ** par,
00234               float     ** sinfo_slit_pos,
00235               int          box_length,
00236               float        y_box,
00237              float        slope_width );
00238 
00239 
00240 
00241 int 
00242 sinfo_new_fit_slits2( cpl_image   * lineImage, 
00243                FitParams ** par,
00244                float     ** sinfo_slit_pos,
00245                int          box_length,
00246                float        y_box,
00247               float        diff_tol );
00248 /*----------------------------------------------------------------------------
00249  *                          Function codes
00250  *--------------------------------------------------------------------------*/
00251 
00270 float 
00271 sinfo_new_edge ( float * xdat, float * parlist/*, int * npar, int * ndat*/ )
00272 {
00273     float return_value ;
00274     float slope1 ;
00275 
00276     /* compute the slopes */
00277     slope1 = ( parlist[3] - parlist[2] ) / ( parlist[1] - parlist[0] ) ;
00278 
00279     /* now build the hat function out of the parameters */
00280     if ( xdat[0] <= parlist[0] )
00281     {
00282         return_value = parlist[2] ;
00283     }
00284     else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[1] )
00285     {
00286         return_value = (xdat[0] - parlist[0]) * slope1 + parlist[2] ;
00287     }
00288     else if ( xdat[0] > parlist[1] )
00289     {
00290         return_value = parlist[3] ;
00291     }
00292     else
00293     {
00294         return_value = 0. ;
00295     }
00296     
00297     return return_value ;
00298 }
00299 
00329 float 
00330 sinfo_new_hat2 ( float * xdat, float * parlist/*, int * npar, int * ndat*/ )
00331 {
00332     float return_value ;
00333     float slope1, slope2, slope3 ;
00334 
00335     /* compute the slopes */
00336     slope1 = ( parlist[6] - parlist[4] ) / ( parlist[1] - parlist[0] ) ;
00337     slope2 = ( parlist[7] - parlist[5] ) / ( parlist[3] - parlist[2] ) ;
00338     slope3 = ( parlist[7] - parlist[6] ) / ( parlist[2] - parlist[1] ) ;
00339 
00340     /* now build the hat function out of the parameters */
00341     if ( xdat[0] <= parlist[0] )
00342     {
00343         return_value = parlist[4] ;
00344     }
00345     else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[1] )
00346     {
00347         return_value = (xdat[0] - parlist[0]) * slope1 + parlist[4] ;
00348     }
00349     else if ( xdat[0] > parlist[1] && xdat[0] <= parlist[2] )
00350     {
00351         return_value = (xdat[0] - parlist[1]) * slope3 + parlist[6] ;
00352     }
00353     else if ( xdat[0] > parlist[2] && xdat[0] <= parlist[3] )
00354     {
00355         return_value = (parlist[3] - xdat[0]) * slope2 + parlist[5] ;
00356     }
00357     else if ( xdat[0] > parlist[3] )
00358     {
00359         return_value = parlist[5] ;
00360     }
00361     else
00362     {
00363         return_value = 0. ;
00364     }
00365     
00366     return return_value ;
00367 }
00368 
00369 
00399 float 
00400 sinfo_new_hat1 ( float * xdat, float * parlist/*, int * npar, int * ndat*/ )
00401 {
00402     float return_value=0 ;
00403     float slope1, slope2 ;
00404 
00405    /* take only positive values for the fit parameters */
00406 
00407     /* compute the slopes */
00408     slope1 = (parlist[4] - parlist[2]) / slopewidth ;
00409     slope2 = (parlist[4] - parlist[3]) / slopewidth ;
00410 
00411     /* now build the hat function out of the parameters */
00412     if ( xdat[0] <= parlist[0] )
00413     {
00414         return_value = parlist[2] ;
00415     }
00416     else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[0] + slopewidth )
00417     {
00418         return_value = (xdat[0] - parlist[0]) * slope1 + parlist[2] ;
00419     }
00420     else if ( xdat[0] > parlist[0] + slopewidth && 
00421               xdat[0] <= parlist[1] - slopewidth )
00422     {
00423         return_value = parlist[4] ;
00424     }
00425     else if ( xdat[0] > parlist[1] - slopewidth && xdat[0] <= parlist[1] )
00426     {
00427         return_value = (parlist[1] - xdat[0]) * slope2 + parlist[3] ;
00428     }
00429     else if ( xdat[0] > parlist[1] )
00430     {
00431         return_value = parlist[3] ;
00432     }
00433     
00434     return return_value ;
00435 }
00436        
00437 
00464 void 
00465 sinfo_new_edge_deriv( float * xdat, float * parlist, 
00466                       float * dervs/*, int * npar*/ )
00467 {
00468     float deriv1_slope1 ;
00469 
00470     /* compute the slopes */
00471     deriv1_slope1 =( parlist[3] - parlist[2] ) / SQR(parlist[1] - parlist[0]) ;
00472 
00473     /* now build the hat derivatives out of the parameters */
00474     if ( xdat[0] <= parlist[0] )
00475     {
00476         dervs[0] = 0. ;
00477         dervs[1] = 0. ;
00478         dervs[2] = 1. ;
00479         dervs[3] = 0. ;
00480     }
00481     else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[1] )
00482     {
00483         dervs[0] = ( xdat[0] - parlist[1] ) * deriv1_slope1  ;
00484         dervs[1] = ( parlist[0] - xdat[0] ) * deriv1_slope1 ;
00485         dervs[2] = ( parlist[0] - xdat[0] ) / ( parlist[1] - parlist[0] ) + 1.;
00486         dervs[3] = ( xdat[0] - parlist[0] ) / ( parlist[1] - parlist[0] ) ;
00487     }
00488     else if ( xdat[0] > parlist[1] )
00489     {
00490         dervs[0] = 0. ;
00491         dervs[1] = 0. ;
00492         dervs[2] = 0. ;
00493         dervs[3] = 1. ;
00494     }
00495 }
00530 void 
00531 sinfo_new_hat_deriv2(float * xdat, float * parlist, 
00532                      float * dervs/*, int * npar*/ )
00533 {
00534     float deriv1_slope1 ;
00535     float deriv1_slope2 ;
00536     float deriv1_slope3 ;
00537 
00538     /* compute the slopes */
00539     deriv1_slope1 = ( parlist[6] - parlist[4] ) / SQR(parlist[1] - parlist[0]);
00540     deriv1_slope2 = ( parlist[7] - parlist[5] ) / SQR(parlist[3] - parlist[2]);
00541     deriv1_slope3 = ( parlist[7] - parlist[6] ) / SQR(parlist[2] - parlist[1]);
00542 
00543     /* now build the hat derivatives out of the parameters */
00544     if ( xdat[0] <= parlist[0] )
00545     {
00546         dervs[0] = 0. ;
00547         dervs[1] = 0. ;
00548         dervs[2] = 0. ;
00549         dervs[3] = 0. ;
00550         dervs[4] = 1. ;
00551         dervs[5] = 0. ;
00552         dervs[6] = 0. ;
00553         dervs[7] = 0. ;
00554     }
00555     else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[1] )
00556     {
00557         dervs[0] = ( xdat[0] - parlist[1] ) * deriv1_slope1  ;
00558         dervs[1] = ( parlist[0] - xdat[0] ) * deriv1_slope1 ;
00559         dervs[2] = 0. ;
00560         dervs[3] = 0. ;
00561         dervs[4] = ( parlist[0] - xdat[0] ) / ( parlist[1] - parlist[0] ) + 1.;
00562         dervs[5] = 0. ;
00563         dervs[6] = ( xdat[0] - parlist[0] ) / ( parlist[1] - parlist[0] ) ;
00564         dervs[7] = 0. ;
00565     }
00566     else if ( xdat[0] > parlist[1] && xdat[0] <= parlist[2] )
00567     {
00568         dervs[0] = 0. ;
00569         dervs[1] = (xdat[0] - parlist[2]) * deriv1_slope3 ;
00570         dervs[2] = (parlist[1] - xdat[0]) * deriv1_slope3 ;
00571         dervs[3] = 0. ;
00572         dervs[4] = 0. ;
00573         dervs[5] = 0. ;
00574         dervs[6] = (parlist[1] - xdat[0]) / (parlist[2] - parlist[1]) + 1. ;
00575         dervs[7] = (xdat[0] - parlist[1]) / (parlist[2] - parlist[1]) ;
00576     }
00577     else if ( xdat[0] > parlist[2] && xdat[0] <= parlist[3] )
00578     {
00579         dervs[0] = 0. ;
00580         dervs[1] = 0. ;
00581         dervs[2] = ( parlist[3] - xdat[0] ) * deriv1_slope2 ;
00582         dervs[3] = ( xdat[0] - parlist[2] ) * deriv1_slope2 ;
00583         dervs[4] = 0. ; 
00584         dervs[5] = ( xdat[0] - parlist[3] ) / ( parlist[3] - parlist[2] ) + 1.;
00585         dervs[6] = 0. ;
00586         dervs[7] = ( parlist[3] - xdat[0] ) / ( parlist[3] - parlist[2] ) ;
00587     }
00588     else if ( xdat[0] > parlist[3] )
00589     {
00590         dervs[0] = 0. ;
00591         dervs[1] = 0. ;
00592         dervs[2] = 0. ;
00593         dervs[3] = 0. ;
00594         dervs[4] = 0. ;
00595         dervs[5] = 1. ;
00596         dervs[6] = 0. ;
00597         dervs[7] = 0. ;
00598     }
00599 }
00600 
00629 void 
00630 sinfo_new_hat_deriv1( float * xdat, float * parlist, 
00631                       float * dervs/*, int * npar*/ )
00632 {
00633     /* now build the hat derivatives out of the parameters */
00634     if ( xdat[0] <= parlist[0] )
00635     {
00636         dervs[0] = 0. ;
00637         dervs[1] = 0. ;
00638         dervs[2] = 1. ;
00639         dervs[3] = 0. ;
00640         dervs[4] = 0. ;
00641     }
00642     else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[0] + slopewidth )
00643     {
00644         dervs[0] = ( parlist[2] - parlist[4] ) / slopewidth ;
00645         dervs[1] = 0. ;
00646         dervs[2] = (( parlist[0] - xdat[0] ) / slopewidth ) + 1. ;
00647         dervs[3] = 0. ;
00648         dervs[4] = ( xdat[0] - parlist[0] ) / slopewidth ;
00649     }
00650     else if ( xdat[0] > parlist[0] + slopewidth && xdat[0] <= 
00651               parlist[1] - slopewidth )
00652     {
00653         dervs[0] = 0. ;
00654         dervs[1] = 0. ;
00655         dervs[2] = 0. ;
00656         dervs[3] = 0. ;
00657         dervs[4] = 1. ;
00658     }
00659     else if ( xdat[0] > parlist[1] - slopewidth && xdat[0] <= parlist[1] )
00660     {
00661         dervs[0] = 0. ;
00662         dervs[1] = ( parlist[4] - parlist[3] ) / slopewidth ;
00663         dervs[2] = 0. ;
00664         dervs[3] = (( xdat[0] - parlist[1] ) / slopewidth ) + 1. ;
00665         dervs[4] = ( parlist[1] - xdat[0] ) / slopewidth ; 
00666     }
00667     else if ( xdat[0] > parlist[1] )
00668     {
00669         dervs[0] = 0. ;
00670         dervs[1] = 0. ;
00671         dervs[2] = 0. ;
00672         dervs[3] = 1. ;
00673         dervs[4] = 0. ;
00674     }
00675 }
00676    
00688 static int 
00689 sinfo_new_inv_mat_edge (void)
00690 {
00691     double even ;
00692     double hv[NPAR] ;
00693     double mjk ;
00694     double rowmax ;
00695     int evin ;
00696     int i, j, k, row ;
00697     int per[NPAR] ;
00698    
00699     /* set permutation array */
00700     for ( i = 0 ; i < nfree ; i++ )
00701     {
00702         per[i] = i ;
00703     }
00704     
00705     for ( j = 0 ; j < nfree ; j++ ) /* in j-th column */
00706     {
00707         /* determine largest element of a row */                           
00708         rowmax = fabs ( matrix2[j][j] ) ;
00709         row = j ;                         
00710 
00711         for ( i = j + 1 ; i < nfree ; i++ )
00712         {
00713             if ( fabs ( matrix2[i][j] ) > rowmax )
00714             {
00715                 rowmax = fabs( matrix2[i][j] ) ;
00716                 row = i ;
00717             }
00718         }
00719 
00720         /* determinant is zero! */
00721         if ( matrix2[row][j] == 0.0 )
00722         {
00723             return -6 ;
00724         }
00725         
00726         /* if the largest element is not on the diagonal, then permutate rows */
00727         if ( row > j )
00728         {
00729             for ( k = 0 ; k < nfree ; k++ )
00730             {
00731                 even = matrix2[j][k] ;
00732                 matrix2[j][k] = matrix2[row][k] ;
00733                 matrix2[row][k] = even ;
00734             }
00735             /* keep track of permutation */
00736             evin = per[j] ;
00737             per[j] = per[row] ;
00738             per[row] = evin ;
00739         }
00740         
00741         /* modify column */
00742         even = 1.0 / matrix2[j][j] ;
00743         for ( i = 0 ; i < nfree ; i++ )
00744         {
00745             matrix2[i][j] *= even ;
00746         }
00747         matrix2[j][j] = even ;
00748         
00749         for ( k = 0 ; k < j ; k++ )
00750         {
00751             mjk = matrix2[j][k] ;
00752             for ( i = 0 ; i < j ; i++ )
00753             {
00754                 matrix2[i][k] -= matrix2[i][j] * mjk ;
00755             }
00756             for ( i = j + 1 ; i < nfree ; i++ )
00757             {
00758                 matrix2[i][k] -= matrix2[i][j] * mjk ;
00759             }
00760             matrix2[j][k] = -even * mjk ;
00761         }
00762     
00763         for ( k = j + 1 ; k < nfree ; k++ )
00764         {
00765             mjk = matrix2[j][k] ;
00766             for ( i = 0 ; i < j ; i++ )
00767             {
00768                 matrix2[i][k]  -= matrix2[i][j] * mjk ;
00769             }
00770             for ( i = j + 1 ; i < nfree ; i++ )
00771             {
00772                 matrix2[i][k]  -= matrix2[i][j] * mjk ;
00773             }
00774             matrix2[j][k] = -even * mjk ;
00775         }
00776     }
00777     
00778     /* finally, repermute the columns */
00779     for ( i = 0 ; i < nfree ; i++ )
00780     {
00781         for ( k = 0 ; k < nfree ; k++ )
00782         {
00783             hv[per[k]] = matrix2[i][k] ;
00784         }
00785         for ( k = 0 ; k < nfree ; k++ )
00786         {
00787             matrix2[i][k] = hv[k] ;
00788         }
00789     }
00790         
00791     /* all is well */
00792     return 0 ;
00793 }
00794     
00813 static void 
00814 sinfo_new_get_mat_edge ( float * xdat,
00815                            int   * xdim,
00816                            float * ydat,
00817                            float * wdat,
00818                            int   * ndat,
00819                            float * fpar,
00820                            float * epar/*,
00821                            int   * npar*/ )
00822 {
00823     double wd ;
00824     double wn ;
00825     double yd ;
00826     int i, j, n ;
00827 
00828     for ( j = 0 ; j < nfree ; j++ )
00829     {
00830         vec[j] = 0.0 ; /* zero sinfo_vector */
00831         for ( i = 0 ; i<= j ; i++ )   /* zero matrix only on and 
00832                                          below diagonal */
00833         {
00834             matrix1[j][i] = 0.0 ;
00835         }
00836     }
00837     chi2 = 0.0 ;  /* reset reduced chi-squared */
00838     
00839     /* loop through data points */
00840     for ( n = 0 ; n < (*ndat) ; n++ )
00841     {
00842         wn = wdat[n] ;
00843         if ( wn > 0.0 )  /* legal weight ? */
00844         {
00845             yd = ydat[n] - sinfo_new_edge( &xdat[(*xdim) * n], 
00846                                            fpar/*, npar, ndat*/ ) ;
00847             sinfo_new_edge_deriv( &xdat[(*xdim) * n], fpar, epar/*, npar */) ;
00848             chi2 += yd * yd * wn ; /* add to chi-squared */
00849             for ( j = 0 ; j < nfree ; j++ )
00850             {
00851                 wd = epar[parptr[j]] * wn ;  /* weighted derivative */
00852                 vec[j] += yd * wd ;       /* fill sinfo_vector */
00853                 for ( i = 0 ; i <= j ; i++ ) /* fill sinfo_matrix */
00854                 {
00855                     matrix1[j][i] += epar[parptr[i]] * wd ;
00856                 }
00857             }
00858         }
00859     }                   
00860 }  
00861                 
00862             
00863 
00893 static int 
00894 sinfo_new_get_vec_edge ( float * xdat,
00895                           int   * xdim,
00896                           float * ydat,
00897                           float * wdat,
00898                           int   * ndat,
00899                           float * fpar,
00900                           float * epar,
00901                           int   * npar )
00902 {
00903     double dj ;
00904     double dy ;
00905     double mii ;
00906     double mji ;
00907     double mjj ;
00908     double wn ;
00909     int i, j, n, r ;
00910 
00911     /* loop to modify and scale the sinfo_matrix */
00912     for ( j = 0 ; j < nfree ; j++ )
00913     {
00914         mjj = matrix1[j][j] ;
00915         if ( mjj <= 0.0 )             /* diagonal element wrong */
00916         {
00917             return -5 ;
00918         }
00919         mjj = sqrt( mjj ) ;
00920         for ( i = 0 ; i < j ; i++ )
00921         {
00922             mji = matrix1[j][i] / mjj / sqrt( matrix1[i][i] ) ;
00923             matrix2[i][j] = matrix2[j][i] = mji ;
00924         }
00925         matrix2[j][j] = 1.0 + labda ;  /* scaled value on diagonal */
00926     }    
00927     
00928     if ( (r = sinfo_new_inv_mat_edge()) ) /* sinfo_invert sinfo_matrix inlace */
00929     {
00930         return r ;
00931     }
00932     
00933     for ( i = 0 ; i < (*npar) ; i ++ )
00934     {
00935         epar[i] = fpar[i] ;
00936     }
00937     
00938     /* loop to calculate correction sinfo_vector */
00939     for ( j = 0 ; j < nfree ; j++ )
00940     {
00941         dj = 0.0 ;
00942         mjj = matrix1[j][j] ;
00943         if ( mjj <= 0.0)               /* not allowed */
00944         {
00945             return -7 ;
00946         }
00947         mjj = sqrt ( mjj ) ;
00948         for ( i = 0 ; i < nfree ; i++ )
00949         {
00950             mii = matrix1[i][i] ;
00951             if ( mii <= 0.0 )
00952             {
00953                 return -7 ;
00954             }
00955             mii = sqrt( mii ) ;
00956             dj += vec[i] * matrix2[j][i] / mjj / mii ;
00957         }
00958         epar[parptr[j]] += dj ;       /* new parameters */
00959     }    
00960     chi1 = 0.0 ;                      /* reset reduced chi-squared */
00961  
00962     /* loop through the data points */
00963     for ( n = 0 ; n < (*ndat) ; n++ )
00964     {
00965         wn = wdat[n] ;               /* get weight */
00966         if ( wn > 0.0 )              /* legal weight */
00967         {
00968             dy = ydat[n] - sinfo_new_edge( &xdat[(*xdim) * n], epar
00969                                                 /*, npar, ndat*/) ;
00970             chi1 += wdat[n] * dy * dy ;
00971         }
00972     }
00973     return 0 ;
00974 }   
00975     
00976         
01030 int 
01031 sinfo_new_lsqfit_edge ( float * xdat,
01032                   int   * xdim,
01033                   float * ydat,
01034                   float * wdat,
01035                   int   * ndat,
01036                   float * fpar,
01037                   float * epar,
01038                   int   * mpar,
01039                   int   * npar,
01040                   float * tol ,
01041                   int   * its ,
01042                   float * lab  )
01043 {
01044     int i, n, r ;
01045     int itc ;                      /* fate of fit */
01046     int found ;                    /* fit converged: 1, not yet converged: 0 */
01047     int  nuse ;                    /* number of useable data points */
01048     double tolerance ;             /* accuracy */
01049 
01050     itc   = 0 ;                    /* fate of fit */
01051     found = 0 ;                    /* reset */
01052     nfree = 0 ;                    /* number of free parameters */
01053     nuse  = 0 ;                    /* number of legal data points */
01054 
01055     if ( *tol < (FLT_EPSILON * 10.0 ) )
01056     {
01057         tolerance = FLT_EPSILON * 10.0 ;  /* default tolerance */
01058     }
01059     else
01060     {
01061         tolerance = *tol ;                /* tolerance */
01062     }
01063     
01064     labda = fabs( *lab ) * LABFACA ;      /* start value for mixing parameter */
01065     for ( i = 0 ; i < (*npar) ; i++ )
01066     {
01067         if ( mpar[i] )
01068         {
01069             if ( nfree > NPAR )         /* too many free parameters */
01070             {
01071                 return -1 ;
01072             }
01073             parptr[nfree++] = i ;         /* a free parameter */
01074         }
01075     }
01076     
01077     if (nfree == 0)                       /* no free parameters */     
01078     {
01079         return -2 ;
01080     }
01081     
01082     for ( n = 0 ; n < (*ndat) ; n++ )
01083     {
01084         if ( wdat[n] > 0.0 )              /* legal weight */
01085         {
01086             nuse ++ ;
01087         }
01088     }
01089     
01090     if ( nfree >= nuse )
01091     {
01092         return -3 ;                       /* no degrees of freedom */
01093     }
01094     if ( labda == 0.0 )                   /* linear fit */
01095     {
01096         /* initialize fpar array */
01097         for ( i = 0 ; i < nfree ; fpar[parptr[i++]] = 0.0 ) ; 
01098         sinfo_new_get_mat_edge(xdat,xdim,ydat,wdat,ndat,fpar,epar/*, npar */) ;
01099         r =  sinfo_new_get_vec_edge ( xdat, xdim, ydat, wdat, ndat, 
01100                                       fpar, epar, npar ) ;
01101         if ( r )                         /* error */
01102         {
01103             return r ;
01104         }
01105         for ( i = 0 ; i < (*npar) ; i++ )
01106         {
01107             fpar[i] = epar[i] ;           /* save new parameters */
01108             epar[i] = 0.0 ;               /* and set errors to zero */
01109         }
01110         chi1 = sqrt( chi1 / (double) (nuse - nfree) ) ;
01111         for ( i = 0 ; i < nfree ; i++ )
01112         {
01113             if ( (matrix1[i][i] <= 0.0 ) || (matrix2[i][i] <= 0.0) )
01114             {
01115                 return -7 ;
01116             }
01117             epar[parptr[i]] = chi1 * sqrt( matrix2[i][i] ) / 
01118                                      sqrt( matrix1[i][i] ) ;
01119         }
01120     }
01121     else                                  /* non-linear fit */
01122     {
01123         /*----------------------------------------------------------------
01124          * the non-linear fit uses the steepest descent method in combination
01125          * with the Taylor method. The mixing of these methods is controlled
01126          * by labda. In the outer loop ( called the iteration loop ) we build
01127          * the sinfo_matrix and calculate the correction sinfo_vector. In the 
01128          * inner loop
01129          * (called the interpolation loop) we check whether we have obtained a
01130          * better solution than the previous one. If so, we leave the inner loop
01131          * else we increase labda ( give more weight to the steepest descent 
01132          * method) calculate the correction sinfo_vector and check again. 
01133          * After the inner loop
01134          * we do a final check on the goodness of the fit and if this satisfies
01135          * the tolerance we calculate the errors of the fitted parameters.
01136          */
01137         while ( !found )                  /* iteration loop */
01138         {      
01139             if ( itc++ == (*its) )        /* increase iteration counter */
01140             {
01141                 return -4 ;               
01142             }
01143             sinfo_new_get_mat_edge( xdat, xdim, ydat, wdat, ndat, 
01144                                     fpar, epar/*, npar*/ ) ;
01145             
01146             /*-------------------------------------------------------------
01147              * here we decrease labda since we may assume that each iteration
01148              * brings us closer to the answer.
01149              */
01150             if ( labda > LABMINA )
01151             {
01152                 labda = labda / LABFACA ;         /* decrease labda */
01153             }
01154             r = sinfo_new_get_vec_edge ( xdat, xdim, ydat, wdat, ndat, 
01155                                          fpar, epar, npar ) ;
01156             if ( (int)fpar[1] - (int)fpar[0] <= 0 && fpar[1] / fpar[0] > 0. )
01157             {
01158                 fpar[1] += 1. ;
01159                 continue ;
01160             } 
01161             if ( r )                      /* error */
01162             {
01163                 return r ;
01164             }
01165 
01166             while ( chi1 >= chi2 )        /* interpolation loop */
01167             {
01168                 /*-----------------------------------------------------------
01169                  * The next statement is based on experience, not on the 
01170                  * mathematics of the problem. It is assumed that we have 
01171                  * reached convergence when the pure steepest descent method 
01172                  * does not produce a better solution.
01173                  */
01174                 if ( labda > LABMAXA )    /* assume solution found */
01175                 {
01176                     break ;
01177                 }
01178                 labda = labda * LABFACA ;   /* increase mixing parameter */
01179                 r = sinfo_new_get_vec_edge ( xdat, xdim, ydat, wdat, 
01180                                              ndat, fpar, epar, npar ) ;
01181                 if ( (int)fpar[1] - (int)fpar[0] <= 0 && 
01182                      fpar[1] / fpar[0] > 0. )
01183                 {
01184                     fpar[1] += 1. ;
01185                     continue ;
01186                 } 
01187                 if ( r )                  /* error */
01188                 {
01189                     return r ;
01190                 }
01191             }
01192 
01193             if ( labda <= LABMAXA )        /* save old parameters */
01194             {
01195                 for ( i = 0 ; i < *npar ; i++ )
01196                 {
01197                     fpar[i] = epar[i] ;
01198                 }
01199             }
01200             if ( (fabs( chi2 - chi1 ) <= (tolerance * chi1)) || 
01201                   (labda > LABMAXA) )
01202             {
01203                 /*------------------------------------------------------------
01204                  * we have a satisfying solution, so now we need to calculate 
01205                  * the correct errors of the fitted parameters. This we do by 
01206                  * using the pure Taylor method because we are very close to 
01207                  * the real solution.
01208                  */
01209                 labda = LABMINA ;              /* for Taylor solution */
01210                 sinfo_new_get_mat_edge ( xdat, xdim, ydat, wdat, ndat, 
01211                                          fpar, epar/*, npar */) ;
01212                 r = sinfo_new_get_vec_edge ( xdat, xdim, ydat, wdat, 
01213                                              ndat, fpar, epar, npar ) ;
01214 
01215                 if ( r )                    /* error */
01216                 {
01217                     return r ;
01218                 }
01219                 for ( i = 0 ; i < (*npar) ; i++ )
01220                 {
01221                     epar[i] = 0.0 ;          /* set error to zero */
01222                 }
01223                 chi2 = sqrt ( chi2 / (double) (nuse - nfree) ) ;
01224 
01225                 for ( i = 0 ; i < nfree ; i++ )
01226                 {
01227                     if ( (matrix1[i][i] <= 0.0) || (matrix2[i][i] <= 0.0) )
01228                     {
01229                         return -7 ;
01230                     }
01231                     epar[parptr[i]] = chi2 * sqrt( matrix2[i][i] ) / 
01232                                              sqrt( matrix1[i][i] ) ;
01233                 }
01234                 found = 1 ;                  /* we found a solution */
01235             }
01236         }
01237     }
01238     return itc ;                             /* return number of iterations */
01239 }
01268 int 
01269 sinfo_new_fit_slits1( cpl_image   * lineImage, 
01270                FitParams ** par,
01271                float     ** sinfo_slit_pos,
01272                int          box_length,
01273                float        y_box )
01274 {
01275     float* position=NULL ;
01276     int   * sinfo_edge, * edgeclean ;
01277     int   * dummyedge ;
01278     int   * pos_row, * pos_rowclean ;
01279     Vector * box_buffer ;
01280     float max_intensity ;
01281     float row_pos ;
01282     int   col ;
01283     int   i, j, k, m, n, ed ;
01284     int   found, init1, init2 ;
01285     int   line ; 
01286     int   column ;
01287     int   slit_length ;
01288     int   agreed ;
01289     int   bad_line ;
01290     int   margin ;
01291     int   iters, xdim, ndat ;
01292     int   numpar, its ;
01293     int   * mpar ;
01294     float * xdat, * wdat ;
01295     float tol, lab ;
01296     float fitpar[NPAR] ;
01297     float dervpar[NPAR] ;
01298     float minval, maxval ;
01299     int ilx=0;
01300     int ily=0;
01301     float* pidata=NULL;
01302 
01303 
01304     if ( NULL == lineImage )
01305     {
01306         sinfo_msg_error("no line image given!" ) ;
01307         return -1 ;
01308     }
01309     ilx=cpl_image_get_size_x(lineImage);
01310     ily=cpl_image_get_size_y(lineImage);
01311     pidata=cpl_image_get_data_float(lineImage);
01312 
01313     slit_length = (int) sqrt (ilx) ;
01314     if ( NULL == par )
01315     {
01316         sinfo_msg_error("no line fit parameters given!") ;
01317         return -2 ;
01318     }
01319 
01320     if ( NULL == sinfo_slit_pos )
01321     {
01322         sinfo_msg_error("no position array allocated!") ;
01323         return -3 ;
01324     }
01325 
01326     if ( box_length <  slit_length ||
01327          box_length >= 3*slit_length )
01328     {
01329         sinfo_msg_error("wrong fitting box length given!" ) ;
01330         return -4 ;
01331     }
01332 
01333     if ( y_box <= 0.  || y_box > 3. )
01334     {
01335         sinfo_msg_error("wrong y box length given!" ) ;
01336         return -5 ;
01337     }
01338 
01339     /* allocate memory for the edges and the row positon of the slitlets */
01340     sinfo_edge         = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01341     dummyedge    = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01342     edgeclean    = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ;
01343     pos_row      = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01344     pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ;
01345 
01346     /* ------------------------------------------------------------------------
01347      * go through the first image columns and the fit parameters and find 
01348      *  the line with the highest intensity 
01349      */
01350     agreed = -1 ;
01351     bad_line = -1 ;
01352     while( agreed == -1 )
01353     {
01354         found = -1 ;
01355         max_intensity = -FLT_MAX ;
01356         for ( col = 0 ; col < box_length ; col++ )
01357         {
01358             for ( i = 0 ; i < par[0]->n_params ; i++ )
01359             {
01360                 if ( par[i]->column == col && par[i]->line != bad_line )
01361                 {
01362                     if ( par[i]->fit_par[0] > max_intensity )
01363                     {
01364                         if ( par[i]->fit_par[1] > 0. )
01365                         {
01366                             max_intensity = par[i]->fit_par[0] ;
01367                             found = i ;
01368                         }
01369                     }
01370                 }
01371             }  
01372         }
01373 
01374         /* --------------------------------------------------------------------
01375          * check if the found line is usable and if the neighbouring line 
01376          * have intensity on near rows in neighbouring slitlets 
01377          */
01378         line    = par[found]->line ;
01379         column  = par[found]->column ;
01380         row_pos = par[found]->fit_par[2] ;
01381         if ( found >= 0 && max_intensity > 0. )
01382         {
01383             for ( i = 0 ; i < par[0]->n_params ; i++ )
01384             {
01385                 if ( par[i]->line == line-1 && 
01386                      par[i]->column == column + slit_length )
01387                 {
01388                     if ( par[i]->fit_par[2] <= (row_pos + y_box) &&
01389                          par[i]->fit_par[2] >= (row_pos - y_box) )
01390                     {
01391                         bad_line = line ;
01392                     } 
01393                 }
01394             }
01395             if ( bad_line != line )
01396             {
01397                 agreed = 1 ;
01398                 break ;
01399             }
01400         }
01401         else 
01402         {
01403           sinfo_msg_error("no emission line found in the first image columns");
01404           return -6 ;
01405         }    
01406     }
01407 
01408  
01409     if ( agreed == -1 )
01410     {
01411         sinfo_msg_error("no emission line found in the first image columns") ;
01412         return -6 ;
01413     }    
01414  
01415     /* now find and store the raw sinfo_edge positions of the found slitlet */ 
01416     n  = 0 ;
01417     ed = 0 ;
01418   
01419 
01420     position=cpl_calloc(ilx,sizeof(float));
01421     for ( col = 0 ; col < ilx ; col++ )
01422     {
01423         for ( i = 0 ; i < par[0]->n_params ; i++ )
01424         {
01425             if ( par[i]->column == col && par[i] -> line == line )
01426             {
01427                 if ( par[i]->fit_par[0] > 0. && par[i]->fit_par[1] > 0. )
01428                 {
01429                     position[n] = par[i]->fit_par[2] ;
01430                     if ( n > 0 && fabs(position[n] - position[n-1]) > y_box )
01431                     {
01432                         sinfo_edge[ed] = col ; 
01433                         pos_row[ed] = sinfo_new_nint( position[n-1] ) ;
01434                         ed++ ;
01435                         if ( col >= ilx - slit_length - 5 ) 
01436                         {
01437                             pos_row[ed] =  sinfo_new_nint( position[n] ) ;
01438                         }
01439                     }
01440                     n++ ;
01441                 }
01442             }
01443         }
01444     }
01445     if ( ed < (slit_length - 1) )
01446     {
01447         sinfo_msg_error("not enough slitlets found") ;
01448         return -7 ;
01449     } 
01450 
01451     /* now find the clean sinfo_edge and row positions of the slitlets */
01452     for ( i = 1 ; i <= ed ; i ++ )
01453     {
01454         if (dummyedge[i-1] != -1)
01455         {
01456             dummyedge[i-1] = sinfo_edge[i-1] ;
01457         }
01458         if ( (sinfo_edge[i] - sinfo_edge[i-1]) < slit_length - 3 ||
01459              (sinfo_edge[i] - sinfo_edge[i-1]) > slit_length + 3 )
01460         {
01461             dummyedge[i]   = -1 ;
01462         }
01463         if ( (sinfo_edge[i+1] - sinfo_edge[i]) < slit_length - 3 ||
01464              (sinfo_edge[i+1] - sinfo_edge[i]) > slit_length + 3 )
01465         {
01466             dummyedge[i+1] = -1 ; 
01467         }
01468     }
01469     
01470     k = 0 ;
01471     for ( i = 0 ; i < ed ; i++ )
01472     {
01473         if ( dummyedge[i] != -1 && dummyedge[i] != 0 )
01474         {
01475             edgeclean[k] = dummyedge[i] ;
01476             pos_rowclean[k] = pos_row[i] ;
01477             k++ ;
01478             if( edgeclean[k-1] > (ilx - slit_length - 6 ) )
01479             {
01480                 pos_rowclean[k] = pos_row[ed] ;
01481             }
01482         }
01483     }
01484 
01485     if ( k != slit_length - 1 )
01486     {
01487         sinfo_msg_error("not enough clean slitlets found") ;
01488         return -7 ;
01489     } 
01490 
01491     /* determine the margins of the fitting box outside the slitlets */
01492     margin = ( box_length - slit_length ) / 2 ;
01493 
01494     /* now go through the slitlets and store the intensity in a 
01495        buffer sinfo_vector */
01496     for ( j = 0 ; j < k ; j++ )
01497     {
01498         m = 0 ;
01499         if ( j == 0 )
01500         {
01501             box_buffer = sinfo_new_vector( edgeclean[0] + margin ) ;
01502             for( col = 0 ; col < edgeclean[0] + margin ; col++ )
01503             {
01504                 box_buffer->data[m] = pidata[col + ilx*pos_rowclean[0]] ;
01505                 m++ ;
01506             }
01507             fitpar[0] = 3. ;
01508             fitpar[1] = 5. ;
01509             fitpar[2] = (float) edgeclean[0] - 1. ;
01510             fitpar[3] = (float) edgeclean[0] + 1. ;
01511           
01512         }
01513         else if ( j < k - 1 )
01514         {
01515             box_buffer = sinfo_new_vector( edgeclean[j] - 
01516                                            edgeclean[j-1] + 2*margin ) ;
01517             for ( col = edgeclean[j - 1] - margin ; 
01518                   col < edgeclean[j] + margin ; col++ )
01519             {
01520                 box_buffer->data[m] = pidata[col + ilx*pos_rowclean[j]] ;
01521                 m++ ;
01522             }
01523             fitpar[0] = (float) margin - 1. ;
01524             fitpar[1] = (float) margin + 1. ;
01525             fitpar[2] = (float) (edgeclean[j] - edgeclean[j-1] + margin) - 1. ;
01526             fitpar[3] = (float) (edgeclean[j] - edgeclean[j-1] + margin) + 1. ;
01527         }
01528         /*else if ( j == k - 1 )*/
01529         else
01530         {
01531             box_buffer = sinfo_new_vector( ilx - edgeclean[k-1] + margin ) ;
01532             for ( col = edgeclean[k - 1] - margin ; col < ilx ; col++ )
01533             {
01534                 box_buffer->data[m] = pidata[col + ilx*pos_rowclean[k]] ;
01535                 m++ ;
01536             }
01537             fitpar[0] = (float) margin - 1. ;
01538             fitpar[1] = (float) margin + 1. ;
01539             fitpar[2] = (float) (ilx - edgeclean[k-1] + margin) - 3. ;
01540             fitpar[3] = (float) (ilx - edgeclean[k-1] + margin) - 1. ;
01541         }
01542 
01543         xdat = (float *) cpl_calloc( box_buffer -> n_elements, sizeof (float) ) ;
01544         wdat = (float *) cpl_calloc( box_buffer -> n_elements, sizeof (float) ) ;
01545         mpar = (int *)   cpl_calloc( NPAR, sizeof (int) ) ;
01546 
01547         /* set initial values for the fitting routine */
01548         minval =  FLT_MAX ;
01549         maxval = -FLT_MAX ;
01550         for ( i = 0 ; i < box_buffer->n_elements ; i++ )
01551         {
01552             xdat[i] = i ;
01553             wdat[i] = 1.0 ;
01554             if ( box_buffer -> data[i] < minval )
01555             {
01556                 minval = box_buffer -> data[i] ;
01557             }
01558             if ( box_buffer -> data[i] > maxval )
01559             {
01560                 maxval = box_buffer -> data[i] ;
01561             }
01562         }
01563 
01564         fitpar[4] = minval ;
01565         fitpar[5] = minval ;
01566         fitpar[6] = maxval ; 
01567         fitpar[7] = maxval ;
01568 
01569         /* search for both positions of the half intensity of the 
01570            hat within the buffer */
01571         init1 = -1 ; 
01572         init2 = -1 ;
01573         for ( i = 0 ; i < box_buffer->n_elements ; i++ )
01574         {
01575             if ( box_buffer -> data[i] >= ( maxval - minval ) / 2. )
01576             {
01577                 init1 = i ;
01578                 break ;
01579             }
01580         }
01581 
01582         for ( i = box_buffer->n_elements - 1 ; i >= 0  ; i-- )
01583         {
01584             if ( box_buffer -> data[i] >= ( maxval + minval ) / 2. )
01585             {
01586                 init2 = i ;
01587                 break ;
01588             }
01589         }
01590   
01591         /* determine the initial positions from the found values */
01592       /*  if ( init1 != -1 )
01593         {
01594             fitpar[0] = init1 - 1. ;
01595             fitpar[1] = init1 + 1. ;
01596         }
01597         if ( init2 != -1 )
01598         {
01599             fitpar[2] = init2 - 1. ;
01600             fitpar[3] = init2 + 1. ;
01601         } */
01602 
01603         for ( i = 0 ; i < NPAR ; i++ )
01604         {
01605             mpar[i] = 1 ;
01606             dervpar[i] = 0. ;
01607         }
01608      
01609         xdim     = XDIMA ;
01610         ndat     = box_buffer -> n_elements ;
01611         numpar   = NPAR ;
01612         tol      = TOLA ;
01613         lab      = LABA ;
01614         its      = ITSA ;
01615         
01616         /* finally, do the least squares fit over the buffer data */
01617         if ( 0 > ( iters = sinfo_new_lsqfit_edge( xdat, &xdim, 
01618                                                   box_buffer -> data, 
01619                                                   wdat, &ndat, fitpar,
01620                                                   dervpar, mpar, 
01621                                                   &numpar, &tol, 
01622                                                   &its, &lab )) )
01623         {
01624             sinfo_msg_warning("least squares fit failed, error no.: %d", iters) ;
01625             return -8 ;
01626         }
01627 
01628 
01629         /* take care of the column position of the fit boxes to get 
01630            the absolute positions */
01631         if ( j == 0 )
01632         {
01633             sinfo_slit_pos[0][0] = (fitpar[0] + fitpar[1]) / 2. ;
01634             sinfo_slit_pos[0][1] = (fitpar[2] + fitpar[3]) / 2. ;
01635         }
01636         else
01637         {
01638             sinfo_slit_pos[j][0] = (fitpar[0] + fitpar[1]) / 2. + 
01639                                    (float)edgeclean[j-1] - (float)margin ;
01640             sinfo_slit_pos[j][1] = (fitpar[2] + fitpar[3]) / 2. + 
01641                                    (float)edgeclean[j-1] - (float)margin ;
01642         }
01643      
01644         sinfo_slit_pos[k][0] = (fitpar[0] + fitpar[1]) / 2.  + 
01645                                (float)edgeclean[k-1] - (float)margin ;
01646         sinfo_slit_pos[k][1] = (fitpar[2] + fitpar[3]) / 2. + 
01647                                (float)edgeclean[k-1] - (float)margin ;
01648 
01649         sinfo_new_destroy_vector ( box_buffer ) ;
01650         cpl_free( xdat ) ;
01651         cpl_free( wdat ) ;
01652         cpl_free( mpar ) ;
01653     }
01654         
01655     cpl_free( sinfo_edge ) ;
01656     cpl_free( pos_row ) ;
01657     cpl_free( edgeclean ) ;
01658     cpl_free( dummyedge ) ;
01659     cpl_free( pos_rowclean ) ;
01660     cpl_free( position );
01661 
01662     return 0 ;
01663 }
01664                               
01696 int 
01697 sinfo_new_fit_slits( cpl_image   * lineImage, 
01698               FitParams ** par,
01699               float     ** sinfo_slit_pos,
01700               int          box_length,
01701               float        y_box,
01702               float        slope_width )
01703 {
01704     float* position=NULL ;
01705 
01706     int   * sinfo_edge, * edgeclean ;
01707     int   * dummyedge ;
01708     int   * pos_row, * pos_rowclean ;
01709     Vector * box_buffer ;
01710     float max_intensity ;
01711     float row_pos ;
01712     int   col ;
01713     int   i, j, k, m, n, ed ;
01714     int   found ;
01715     int   line ; 
01716     int   column ;
01717     int   slit_length ;
01718     int   agreed ;
01719     int   bad_line ;
01720     int   margin ;
01721     int   iters, xdim, ndat ;
01722     int   numpar, its ;
01723     int   * mpar ;
01724     float * xdat, * wdat ;
01725     float tol, lab ;
01726     float fitpar[NPAR] ;
01727     float dervpar[NPAR] ;
01728     float minval, maxval ;
01729     int ilx=0;
01730     int ily=0;
01731     float* pidata=NULL;
01732 
01733 
01734     if ( NULL == lineImage )
01735     {
01736         sinfo_msg_error("no line image given!" ) ;
01737         return -1 ;
01738     }
01739     ilx=cpl_image_get_size_x(lineImage);
01740     ily=cpl_image_get_size_y(lineImage);
01741     pidata=cpl_image_get_data_float(lineImage);
01742 
01743     slit_length = (int) sqrt (ilx) ;
01744     if ( NULL == par )
01745     {
01746         sinfo_msg_error("no line fit parameters given!" ) ;
01747         return -2 ;
01748     }
01749 
01750     if ( NULL == sinfo_slit_pos )
01751     {
01752         sinfo_msg_error("no position array allocated!" ) ;
01753         return -3 ;
01754     }
01755 
01756     if ( box_length <  slit_length ||
01757          box_length >= 3*slit_length )
01758     {
01759         sinfo_msg_error("wrong fitting box length given!" ) ;
01760         return -4 ;
01761     }
01762 
01763     if ( y_box <= 0.  || y_box > 3. )
01764     {
01765         sinfo_msg_error("wrong y box length given!" ) ;
01766         return -5 ;
01767     }
01768 
01769     if ( slope_width <= 0. )
01770     {
01771         sinfo_msg_error("wrong width of linear slope given!") ;
01772         return -6 ;
01773     }
01774 
01775     /* initialize module global variable slopewidth */
01776     slopewidth = slope_width ;
01777 
01778     /* allocate memory for the edges and the row positon of the slitlets */
01779     sinfo_edge         = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01780     dummyedge    = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01781     edgeclean    = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ;
01782     pos_row      = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01783     pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ;
01784 
01785     /* ------------------------------------------------------------------------
01786      * go through the first image columns and the fit parameters and find 
01787      * the line with the highest intensity 
01788      */
01789     agreed = -1 ;
01790     bad_line = -1 ;
01791     while( agreed == -1 )
01792     {
01793         found = -1 ;
01794         max_intensity = -FLT_MAX ;
01795         for ( col = 0 ; col < box_length ; col++ )
01796         {
01797             for ( i = 0 ; i < par[0]->n_params ; i++ )
01798             {
01799                 if ( par[i]->column == col && par[i]->line != bad_line )
01800                 {
01801                     if ( par[i]->fit_par[0] > max_intensity )
01802                     {
01803                         if ( par[i]->fit_par[1] > 0. )
01804                         {
01805                             max_intensity = par[i]->fit_par[0] ;
01806                             found = i ;
01807                         }
01808                     }
01809                 }
01810             }  
01811         }
01812 
01813         /* --------------------------------------------------------------------
01814          * check if the found line is usable and if the neighbouring line 
01815          * have intensity on near rows in neighbouring slitlets 
01816          */
01817         line    = par[found]->line ;
01818         column  = par[found]->column ;
01819         row_pos = par[found]->fit_par[2] ;
01820         if ( found >= 0 && max_intensity > 0. )
01821         {
01822             for ( i = 0 ; i < par[0]->n_params ; i++ )
01823             {
01824                 if ( par[i]->line == line-1 && 
01825                      par[i]->column == column + slit_length )
01826                 {
01827                     if ( par[i]->fit_par[2] <= (row_pos + y_box) &&
01828                          par[i]->fit_par[2] >= (row_pos - y_box) )
01829                     {
01830                         bad_line = line ;
01831                     } 
01832                 }
01833             }
01834             if ( bad_line != line )
01835             {
01836                 agreed = 1 ;
01837                 break ;
01838             }
01839         }
01840         else 
01841         {
01842           sinfo_msg_error("no emission line found in the first image columns");
01843           return -7 ;
01844         }    
01845     }
01846 
01847  
01848     if ( agreed == -1 )
01849     {
01850         sinfo_msg_error("no emission line found in the first image columns") ;
01851         return -7 ;
01852     }    
01853  
01854     /* now find and store the raw sinfo_edge positions of the found slitlet */ 
01855     n  = 0 ;
01856     ed = 0 ;
01857     position=cpl_calloc(ilx,sizeof(float)) ;
01858 
01859     for ( col = 0 ; col < ilx ; col++ )
01860     {
01861         for ( i = 0 ; i < par[0]->n_params ; i++ )
01862         {
01863             if ( par[i]->column == col && par[i] -> line == line )
01864             {
01865                 if ( par[i]->fit_par[0] > 0. && par[i]->fit_par[1] > 0. )
01866                 {
01867                     position[n] = par[i]->fit_par[2] ;
01868                     if ( n > 0 && fabs(position[n] - position[n-1]) > y_box )
01869                     {
01870                         sinfo_edge[ed] = col ; 
01871                         pos_row[ed] = sinfo_new_nint( position[n-1] ) ;
01872                         ed++ ;
01873                         if ( col >= ilx - slit_length - 5 ) 
01874                         {
01875                             pos_row[ed] =  sinfo_new_nint( position[n] ) ;
01876                         }
01877                     }
01878                     n++ ;
01879                 }
01880             }
01881         }
01882     }
01883     if ( ed < (slit_length - 1) )
01884     {
01885         sinfo_msg_error("not enough slitlets found") ;
01886         return -8 ;
01887     } 
01888 
01889     /* now find the clean sinfo_edge and row positions of the slitlets */
01890     for ( i = 1 ; i <= ed ; i ++ )
01891     {
01892         if (dummyedge[i-1] != -1)
01893         {
01894             dummyedge[i-1] = sinfo_edge[i-1] ;
01895         }
01896         if ( (sinfo_edge[i] - sinfo_edge[i-1]) < slit_length - 3 ||
01897              (sinfo_edge[i] - sinfo_edge[i-1]) > slit_length + 3 )
01898         {
01899             dummyedge[i]   = -1 ;
01900         }
01901         if ( (sinfo_edge[i+1] - sinfo_edge[i]) < slit_length - 3 ||
01902              (sinfo_edge[i+1] - sinfo_edge[i]) > slit_length + 3 )
01903         {
01904             dummyedge[i+1] = -1 ; 
01905         }
01906     }
01907     
01908     k = 0 ;
01909     for ( i = 0 ; i < ed ; i++ )
01910     {
01911         if ( dummyedge[i] != -1 && dummyedge[i] != 0 )
01912         {
01913             edgeclean[k] = dummyedge[i] ;
01914             pos_rowclean[k] = pos_row[i] ;
01915             k++ ;
01916             if( edgeclean[k-1] > (ilx - slit_length - 6 ) )
01917             {
01918                 pos_rowclean[k] = pos_row[ed] ;
01919             }
01920         }
01921     }
01922 
01923     if ( k != slit_length - 1 )
01924     {
01925         sinfo_msg_error ("not enough clean slitlets found") ;
01926         return -7 ;
01927     } 
01928 
01929     /* determine the margins of the fitting box outside the slitlets */
01930     margin = ( box_length - slit_length ) / 2 ;
01931 
01932     /* now go through the slitlets and store the intensity in a buffer 
01933        sinfo_vector */
01934     for ( j = 0 ; j < k ; j++ )
01935     {
01936         m = 0 ;
01937         if ( j == 0 )
01938         {
01939             box_buffer = sinfo_new_vector( edgeclean[0] + margin ) ;
01940             for( col = 0 ; col < edgeclean[0] + margin ; col++ )
01941             {
01942                 box_buffer->data[m] = pidata[col + ilx*pos_rowclean[0]] ;
01943                 m++ ;
01944             }
01945             /* initial values for the fitted positions */
01946             fitpar[0] = 1. ;
01947             fitpar[1] = (float)edgeclean[0] ;
01948           
01949         }
01950         else if ( j < k - 1 )
01951         {
01952             box_buffer = sinfo_new_vector( edgeclean[j] - 
01953                                            edgeclean[j-1] + 2*margin ) ;
01954             for ( col = edgeclean[j - 1] - margin ; 
01955                   col < edgeclean[j] + margin ; col++ )
01956             {
01957                 box_buffer->data[m] = pidata[col + ilx*pos_rowclean[j]] ;
01958                 m++ ;
01959             }
01960             /* initial values for the fitted positions */
01961             fitpar[0] = (float)margin ;
01962             fitpar[1] = (float)(edgeclean[j] - edgeclean[j-1] + margin ) ;
01963         }
01964         /*else if ( j == k - 1 )*/
01965         else
01966         {
01967             box_buffer = sinfo_new_vector( ilx - edgeclean[k-1] + margin ) ;
01968             for ( col = edgeclean[k - 1] - margin ; col < ilx ; col++ )
01969             {
01970                 box_buffer->data[m] = pidata[col + ilx*pos_rowclean[k]] ;
01971                 m++ ;
01972             }
01973             /* initial values for the fitted positions */
01974             fitpar[0] = (float)margin ;
01975             fitpar[1] = (float)(m - 1) ;
01976         }
01977 
01978         xdat=(float *) cpl_calloc( box_buffer -> n_elements, sizeof (float) ) ;
01979         wdat=(float *) cpl_calloc( box_buffer -> n_elements, sizeof (float) ) ;
01980         mpar=(int *)   cpl_calloc( NPAR, sizeof (int) ) ;
01981 
01982         /* set initial values for the fitting routine */
01983         minval =  FLT_MAX ;
01984         maxval = -FLT_MAX ;
01985         for ( i = 0 ; i < box_buffer->n_elements ; i++ )
01986         {
01987             xdat[i] = i ;
01988             wdat[i] = 1.0 ;
01989             if ( box_buffer -> data[i] < minval )
01990             {
01991                 minval = box_buffer -> data[i] ;
01992             }
01993             if ( box_buffer -> data[i] > maxval )
01994             {
01995                 maxval = box_buffer -> data[i] ;
01996             }
01997         }
01998 
01999         for ( i = 0 ; i < NPAR ; i++ )
02000         {
02001             mpar[i] = 1 ;
02002             dervpar[i] = 0. ;
02003         }
02004      
02005         xdim     = XDIMA ;
02006         ndat     = box_buffer -> n_elements ;
02007         numpar   = NPAR ;
02008         tol      = TOLA ;
02009         lab      = LABA ;
02010         its      = ITSA ;
02011         
02012         fitpar[2] = minval ;
02013         fitpar[3] = minval ;
02014         fitpar[4] = maxval ;
02015 
02016         /* finally, do the least squares fit over the buffer data */
02017         if ( 0 > ( iters = sinfo_new_lsqfit_edge( xdat, &xdim, 
02018                                                   box_buffer -> data, 
02019                                                   wdat, &ndat, fitpar,
02020                                                   dervpar, mpar, &numpar, 
02021                                                   &tol, &its, &lab )) )
02022         {
02023             sinfo_msg_warning("least squares fit failed, error no.: %d", iters) ;
02024             return -9 ;
02025         }
02026 
02027 
02028         /* take care of the column position of the fit boxes 
02029            to get the absolute positions */
02030         if ( j == 0 )
02031         {
02032             sinfo_slit_pos[0][0] = fitpar[0] + slopewidth/2. ;
02033             sinfo_slit_pos[0][1] = fitpar[1] - slopewidth/2. ;
02034         }
02035         else
02036         {
02037             sinfo_slit_pos[j][0] = fitpar[0] + slopewidth/2. + 
02038                                    (float)edgeclean[j-1] - (float)margin ;
02039             sinfo_slit_pos[j][1] = fitpar[1] - slopewidth/2. + 
02040                                    (float)edgeclean[j-1] - (float)margin ;
02041         }
02042 
02043         sinfo_slit_pos[k][0] = fitpar[0] + slopewidth/2. + 
02044                                (float)edgeclean[k-1] - (float)margin ;
02045         sinfo_slit_pos[k][1] = fitpar[1] - slopewidth/2. + 
02046                                (float)edgeclean[k-1] - (float)margin ;
02047      
02048         sinfo_new_destroy_vector ( box_buffer ) ;
02049         cpl_free( xdat ) ;
02050         cpl_free( wdat ) ;
02051         cpl_free( mpar ) ;
02052     }
02053         
02054 
02055     cpl_free( sinfo_edge ) ;
02056     cpl_free( pos_row ) ;
02057     cpl_free( edgeclean ) ;
02058     cpl_free( dummyedge ) ;
02059     cpl_free( pos_rowclean ) ;
02060     cpl_free( position );
02061     return 0 ;
02062 }
02063                               
02108 int 
02109 sinfo_new_fit_slits2( cpl_image   * lineImage, 
02110                FitParams ** par,
02111                float     ** sinfo_slit_pos,
02112                int          box_length,
02113                float        y_box,
02114                float        diff_tol )
02115 {
02116     float* position=NULL ;
02117     int   * sinfo_edge, * edgeclean ;
02118     int   * dummyedge ;
02119     int   * pos_row, * pos_rowclean ;
02120     Vector * box_buffer ;
02121     Vector * half_buffer ;
02122     float max_intensity ;
02123     float row_pos ;
02124     int   col ;
02125     int   i, j, k, m, n, ed ;
02126     int   found, init1 ;
02127     int   line ; 
02128     int   nel, n_right, left_right ;
02129     int   column ;
02130     int   slit_length ;
02131     int   agreed ;
02132     int   bad_line ;
02133     int   margin ;
02134     int   iters, xdim, ndat ;
02135     int   numpar, its ;
02136     int   * mpar ;
02137     float * xdat, * wdat ;
02138     float tol, lab ;
02139     float fitpar[NPAR] ;
02140     float dervpar[NPAR] ;
02141     float minval, maxval ;
02142     float pos, last_pos ;
02143     int ilx=0;
02144     int ily=0;
02145     float* pidata=NULL;
02146 
02147 
02148     if ( NULL == lineImage )
02149     {
02150         sinfo_msg_error("no line image given!" ) ;
02151         return -1 ;
02152     }
02153     ilx=cpl_image_get_size_x(lineImage);
02154     ily=cpl_image_get_size_y(lineImage);
02155     pidata=cpl_image_get_data_float(lineImage);
02156 
02157     slit_length = (int) sqrt (ilx) ;
02158 
02159     if ( NULL == par )
02160     {
02161         sinfo_msg_error("no line fit parameters given!" ) ;
02162         return -2 ;
02163     }
02164 
02165     if ( NULL == sinfo_slit_pos )
02166     {
02167         sinfo_msg_error("no position array allocated!" ) ;
02168         return -3 ;
02169     }
02170 
02171     if ( box_length <  slit_length ||
02172          box_length >= 3*slit_length )
02173     {
02174         sinfo_msg_error("wrong fitting box length given!" ) ;
02175         return -4 ;
02176     }
02177 
02178     if ( y_box <= 0.  || y_box > 3. )
02179     {
02180         sinfo_msg_error("wrong y box length given!" ) ;
02181         return -5 ;
02182     }
02183 
02184     if ( diff_tol < 1. )
02185     {
02186         sinfo_msg_error("diff_tol too small!" ) ;
02187         return -6 ;
02188     }
02189 
02190     /* allocate memory for the edges and the row positon of the slitlets */
02191     sinfo_edge         = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
02192     dummyedge    = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
02193     edgeclean    = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ;
02194     pos_row      = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
02195     pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ;
02196 
02197     /* ------------------------------------------------------------------------
02198      * go through the first image columns and the fit parameters and find 
02199      * the line with the highest intensity 
02200      */
02201     agreed = -1 ;
02202     bad_line = -1 ;
02203     while( agreed == -1 )
02204     {
02205         found = -1 ;
02206         max_intensity = -FLT_MAX ;
02207         for ( col = 0 ; col < box_length ; col++ )
02208         {
02209             for ( i = 0 ; i < par[0]->n_params ; i++ )
02210             {
02211                 if ( par[i]->column == col && par[i]->line != bad_line )
02212                 {
02213                     if ( par[i]->fit_par[0] > max_intensity )
02214                     {
02215                         if ( par[i]->fit_par[1] > 0. )
02216                         {
02217                             max_intensity = par[i]->fit_par[0] ;
02218                             found = i ;
02219                         }
02220                     }
02221                 }
02222             }  
02223         }
02224 
02225         /* --------------------------------------------------------------------
02226          * check if the found line is usable and if the neighbouring line 
02227          * have intensity on near rows in neighbouring slitlets 
02228          */
02229         line    = par[found]->line ;
02230         column  = par[found]->column ;
02231         row_pos = par[found]->fit_par[2] ;
02232         if ( found >= 0 && max_intensity > 0. )
02233         {
02234             for ( i = 0 ; i < par[0]->n_params ; i++ )
02235             {
02236                 if ( par[i]->line == line-1 && 
02237                      par[i]->column == column + slit_length )
02238                 {
02239                     if ( par[i]->fit_par[2] <= (row_pos + y_box) &&
02240                          par[i]->fit_par[2] >= (row_pos - y_box) )
02241                     {
02242                         bad_line = line ;
02243                     } 
02244                 }
02245             }
02246             if ( bad_line != line )
02247             {
02248                 agreed = 1 ;
02249                 break ;
02250             }
02251         }
02252         else 
02253         {
02254           sinfo_msg_error("no emission line found in the first image columns");
02255           return -7 ;
02256         }    
02257     }
02258 
02259  
02260     if ( agreed == -1 )
02261     {
02262         sinfo_msg_error("no emission line found in the first image columns") ;
02263         return -7 ;
02264     }    
02265  
02266     /* now find and store the raw sinfo_edge positions of the found slitlet */ 
02267     n  = 0 ;
02268     ed = 0 ;
02269     position=cpl_calloc(ilx,sizeof(float)) ;
02270 
02271     for ( col = 0 ; col < ilx ; col++ )
02272     {
02273         for ( i = 0 ; i < par[0]->n_params ; i++ )
02274         {
02275             if ( par[i]->column == col && par[i]->line == line )
02276             {
02277                 if ( par[i]->fit_par[0] > 0. && par[i]->fit_par[1] > 0. )
02278                 {
02279                     position[n] = par[i]->fit_par[2] ;
02280                     if ( n > 0 && fabs(position[n] - position[n-1]) > y_box )
02281                     {
02282                         sinfo_edge[ed] = col ; 
02283                         pos_row[ed] = sinfo_new_nint( position[n-1] ) ;
02284                         ed++ ;
02285                         if ( col >= ilx - slit_length - 5 ) 
02286                         {
02287                             pos_row[ed] =  sinfo_new_nint( position[n] ) ;
02288                         }
02289                     }
02290                     n++ ;
02291                 }
02292             }
02293         }
02294     }
02295     if ( ed < (slit_length - 1) )
02296     {
02297         sinfo_msg_error("not enough slitlets found") ;
02298         return -8 ;
02299     } 
02300 
02301     /* now find the clean sinfo_edge and row positions of the slitlets */
02302     for ( i = 1 ; i <= ed ; i ++ )
02303     {
02304         if (dummyedge[i-1] != -1)
02305         {
02306             dummyedge[i-1] = sinfo_edge[i-1] ;
02307         }
02308         if ( (sinfo_edge[i] - sinfo_edge[i-1]) < slit_length - 3 ||
02309              (sinfo_edge[i] - sinfo_edge[i-1]) > slit_length + 3 )
02310         {
02311             dummyedge[i]   = -1 ;
02312         }
02313         if ( (sinfo_edge[i+1] - sinfo_edge[i]) < slit_length - 3 ||
02314              (sinfo_edge[i+1] - sinfo_edge[i]) > slit_length + 3 )
02315         {
02316             dummyedge[i+1] = -1 ; 
02317         }
02318     }
02319     
02320     k = 0 ;
02321     for ( i = 0 ; i < ed ; i++ )
02322     {
02323         if ( dummyedge[i] != -1 && dummyedge[i] != 0 )
02324         {
02325             edgeclean[k] = dummyedge[i] ;
02326             pos_rowclean[k] = pos_row[i] ;
02327             k++ ;
02328             if( edgeclean[k-1] > (ilx - slit_length - 6 ) )
02329             {
02330                 pos_rowclean[k] = pos_row[ed] ;
02331             }
02332         }
02333     }
02334 
02335     if ( k != slit_length - 1 )
02336     {
02337         sinfo_msg_error("not enough clean slitlets found") ;
02338         return -7 ;
02339     } 
02340 
02341     /* determine the margins of the fitting box outside the slitlets */
02342     margin = ( box_length - slit_length ) / 2 ;
02343 
02344     /* now go through the slitlets and store the intensity in a 
02345        buffer sinfo_vector */
02346     for ( j = 0 ; j <= k ; j++ )
02347     {
02348         m = 0 ;
02349         if ( j == 0 )
02350         {
02351             box_buffer = sinfo_new_vector( edgeclean[0] + margin ) ;
02352             for( col = 0 ; col < edgeclean[0] + margin ; col++ )
02353             {
02354                 box_buffer->data[m] = pidata[col +ilx*pos_rowclean[0]] ;
02355                 m++ ;
02356             }
02357         }
02358         else if ( j < k )
02359         {
02360             box_buffer = sinfo_new_vector( edgeclean[j] - 
02361                                            edgeclean[j-1] + 2*margin ) ;
02362             for ( col = edgeclean[j - 1] - margin ; 
02363                   col < edgeclean[j] + margin ; col++ )
02364             {
02365                 box_buffer->data[m] = pidata[col + ilx*pos_rowclean[j]] ;
02366                 m++ ;
02367             }
02368         }
02369         else 
02370         {
02371             box_buffer = sinfo_new_vector( ilx - edgeclean[k-1] + margin ) ;
02372             for ( col = edgeclean[k - 1] - margin ; col < ilx ; col++ )
02373             {
02374                 box_buffer->data[m] = pidata[col + ilx*pos_rowclean[k]] ;
02375                 m++ ;
02376             }
02377         }
02378 
02379         for ( left_right = 0 ; left_right <= 1 ; left_right++ )
02380         { 
02381             nel = 0 ;
02382             if ( left_right == 0 )
02383             {
02384                 nel = box_buffer -> n_elements / 2 ;
02385             }
02386             else
02387             {
02388                 if ( box_buffer -> n_elements % 2 == 0 )
02389                 {
02390                     nel = box_buffer -> n_elements / 2 ;
02391                 }
02392                 else
02393                 {
02394                     nel = box_buffer -> n_elements / 2 + 1 ;
02395                 }
02396             }
02397 
02398             /* now split the buffer in the midth in a left and right 
02399                part for fitting */
02400             half_buffer = sinfo_new_vector( nel ) ;
02401             if ( left_right == 0 )
02402             {
02403                 for ( i = 0 ; i < nel ; i++ )
02404                 {
02405                     half_buffer -> data[i] = box_buffer -> data[i] ;
02406                 }
02407             }
02408             else
02409             {
02410                 n_right = 0 ;
02411                 for ( i = box_buffer -> n_elements - 1 ; 
02412                       i >= box_buffer -> n_elements - nel ; i-- )
02413                 {
02414                     half_buffer -> data[n_right] = box_buffer -> data[i] ;
02415                     n_right++ ;
02416                 }
02417             }
02418 
02419             xdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
02420             wdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
02421             mpar = (int *)   cpl_calloc( NPAR, sizeof (int) ) ;
02422 
02423             /* set initial values for the fitting routine */
02424             minval =  FLT_MAX ;
02425             maxval = -FLT_MAX ;
02426             for ( i = 0 ; i < nel ; i++ )
02427             {
02428                 xdat[i] = i ;
02429                 wdat[i] = 1.0 ;
02430                 if ( half_buffer -> data[i] < minval )
02431                 {
02432                     minval = half_buffer -> data[i] ;
02433                 }
02434                 if ( half_buffer -> data[i] > maxval )
02435                 {
02436                     maxval = half_buffer -> data[i] ;
02437                 }
02438             }
02439 
02440             fitpar[2] = minval ;
02441             fitpar[3] = maxval ; 
02442 
02443             /* search for both positions of the half intensity of 
02444                the hat within the buffer */
02445             init1 = -1 ; 
02446             for ( i = 0 ; i < nel ; i++ )
02447             {
02448                 if ( half_buffer -> data[i] >= ( maxval + minval ) / 2. )
02449                 {
02450                     init1 = i ;
02451                     break ;
02452                 }
02453             }
02454 
02455             /* determine the initial positions from the found values */
02456             if ( init1 != -1 )
02457             {
02458                 fitpar[0] = ((float)init1 - 1.) ;
02459                 fitpar[1] = ((float)init1 + 1.) ;
02460             }
02461 
02462             for ( i = 0 ; i < NPAR ; i++ )
02463             {
02464                 mpar[i] = 1 ;
02465                 dervpar[i] = 0. ;
02466             }
02467       
02468             xdim     = XDIMA ;
02469             ndat     = nel ;
02470             numpar   = NPAR ;
02471             tol      = TOLA ;
02472             lab      = LABA ;
02473             its      = ITSA ;
02474          
02475             /* finally, do the least squares fit over the buffer data */
02476             if ( 0 > ( iters = sinfo_new_lsqfit_edge( xdat, &xdim, 
02477                                                       half_buffer -> data, 
02478                                                       wdat, &ndat, fitpar,
02479                                                       dervpar, mpar, &numpar, 
02480                                                       &tol, &its, &lab )) )
02481             { 
02482                 /* if the fit doesn't succeed the initial values are taken */
02483                 sinfo_msg_warning (" least squares fit failed, error"
02484                                    " no.: %d in slitlet: %d\n", iters, j) ;
02485                 fitpar[0] = ((float)init1 - 1.) ;
02486                 fitpar[1] = ((float)init1 + 1.) ;
02487             }
02488 
02489             pos = (fitpar[0] + fitpar[1]) / 2. ;
02490 
02491             /*---------------------------------------------------------------- 
02492              * now discern the left and the right sinfo_edge fit of the 
02493              * slitlets and associate the fit results with the absolute 
02494              * positions in the image consider the difference of the fitted 
02495              * slit position to the expected position and decide wether the 
02496              * fit is taken or the expected value is taken.
02497              */
02498             if ( left_right == 0 )
02499             {
02500                 /* take care of the column position of the fit boxes 
02501                    to get the absolute positions */
02502                 if ( j == 0 )
02503                 {
02504                     if ( fabs(pos - ((float)edgeclean[0] - 1. - 
02505                          (float)slit_length)) < diff_tol )
02506                     {
02507                         sinfo_slit_pos[0][0] = pos ;
02508                     }
02509                     else
02510                     {
02511                         sinfo_msg_warning("something wrong with fitted left"
02512                                           " position of slitlet 0") ;
02513                         if ( (float) edgeclean[0] - 1. - 
02514                              (float)slit_length < 0. )
02515                         {
02516                             sinfo_slit_pos[0][0] = 0. ;
02517                         }
02518                         else
02519                         {
02520                             sinfo_slit_pos[0][0] = (float)edgeclean[0] - 1. - 
02521                             (float)slit_length ;
02522                         }
02523                     }
02524                 }
02525                 else if ( j < k )
02526                 {
02527                     if ( fabs( pos - (float)margin ) < diff_tol )
02528                     {
02529                         sinfo_slit_pos[j][0] = pos + (float)edgeclean[j-1] - 
02530                                                (float)margin ;
02531                     }
02532                     else
02533                     {
02534                         sinfo_msg_warning("something wrong with fitted left"
02535                                           " position of slitlet %d", j) ;
02536                         sinfo_slit_pos[j][0] = (float)edgeclean[j-1] - 1. ;
02537                     }
02538                 }
02539                 else
02540                 {
02541                     if ( fabs( pos - (float)margin ) < diff_tol )
02542                     {
02543                         sinfo_slit_pos[k][0] = pos + (float)edgeclean[k-1] - 
02544                                                (float)margin ;
02545                     }
02546                     else
02547                     {
02548                         sinfo_msg_warning("something wrong with fitted left"
02549                                           " position of slitlet %d", j) ;
02550                         sinfo_slit_pos[k][0] = (float)edgeclean[k-1] - 1. ;
02551                     }
02552                 }
02553             }
02554             else
02555             {
02556                 /* take care of the column position of the fit boxes 
02557                    to get the absolute positions */
02558                 if ( j == 0 )
02559                 {
02560                     if ( fabs( (float)box_buffer->n_elements - pos - 
02561                                (float)edgeclean[0] ) < diff_tol )
02562                     {
02563                         sinfo_slit_pos[0][1] = (float)(box_buffer->n_elements - 
02564                                                 1) - pos ;
02565                     }
02566                     else
02567                     {
02568                         sinfo_msg_warning("something wrong with fitted "
02569                                           "right position of slitlet 0") ;
02570                         sinfo_slit_pos[0][1] = (float)edgeclean[0] - 1. ;
02571                     }
02572                 }
02573                 else if ( j < k )
02574                 {
02575                     if ( fabs( (float)box_buffer->n_elements - pos
02576                              + (float)edgeclean[j-1] - (float)margin - 
02577                                (float)edgeclean[j] ) < diff_tol )
02578                     {
02579                         sinfo_slit_pos[j][1] = (float)(box_buffer->n_elements -
02580                                         1) - pos
02581                                        + (float)edgeclean[j-1] - (float)margin ;
02582                     }
02583                     else
02584                     {
02585                         sinfo_msg_warning("something wrong with fitted right "
02586                                           "position of slitlet %d", j) ;
02587                         sinfo_slit_pos[j][1] = (float)edgeclean[j] - 1. ;
02588                     }
02589                 }
02590                 else
02591                 {
02592                     if ( edgeclean[k-1] + slit_length > ilx )
02593                     {
02594                         last_pos = (float)(ilx - 1) ;
02595                     }
02596                     else
02597                     {
02598                         last_pos = (float)(edgeclean[k-1] - 1 + slit_length) ;
02599                     }
02600                     if ( fabs( (float)(box_buffer->n_elements - 1) - pos
02601                              + (float)edgeclean[k-1] - (float)margin - 
02602                                last_pos ) < diff_tol )
02603                     {
02604                         sinfo_slit_pos[k][1] = (float)(box_buffer->n_elements -
02605                                                1) - pos
02606                                        + (float)edgeclean[k-1] - (float)margin ;
02607                     }
02608                     else
02609                     {
02610                         sinfo_msg_warning("something wrong with fitted right "
02611                                           "position of slitlet %d\n", j) ;
02612                         sinfo_slit_pos[k][1] = last_pos ;
02613                     }
02614                 }
02615             }
02616 
02617             sinfo_new_destroy_vector ( half_buffer ) ;
02618             cpl_free( xdat ) ;
02619             cpl_free( wdat ) ;
02620             cpl_free( mpar ) ;
02621         }
02622         sinfo_new_destroy_vector ( box_buffer ) ;
02623     }
02624         
02625     cpl_free( sinfo_edge ) ;
02626     cpl_free( pos_row ) ;
02627     cpl_free( edgeclean ) ;
02628     cpl_free( dummyedge ) ;
02629     cpl_free( pos_rowclean ) ;
02630     cpl_free(position);
02631     return 0 ;
02632 }
02633                               
02668 int 
02669 sinfo_new_fit_slits_edge( cpl_image   * lineImage, 
02670                   FitParams ** par,
02671                   float     ** sinfo_slit_pos,
02672                   int          box_length,
02673                   float        y_box,
02674                   float        diff_tol )
02675 {
02676     float* position=NULL ;
02677     int   * sinfo_edge, * edgeclean ;
02678     int   * dummyedge ;
02679     int   * pos_row, * pos_rowclean ;
02680     Vector * box_buffer ;
02681     Vector * half_buffer ;
02682     float max_intensity ;
02683     float row_pos ;
02684     int   row, col ;
02685     int   i, j, k, m, n, ed ;
02686     int   found, init1 ;
02687     int   line ; 
02688     int   nel, n_right, left_right ;
02689     int   column ;
02690     int   slit_length ;
02691     int   agreed ;
02692     int   bad_line ;
02693     int   margin ;
02694     int   iters, xdim, ndat ;
02695     int   numpar, its ;
02696     int   * mpar ;
02697     float * xdat, * wdat ;
02698     float tol, lab ;
02699     float fitpar[NPAR] ;
02700     float dervpar[NPAR] ;
02701     float minval, maxval ;
02702     float pos, last_pos ;
02703     int ilx=0;
02704     int ily=0;
02705     float* pidata=NULL;
02706 
02707 
02708     if ( NULL == lineImage )
02709     {
02710         sinfo_msg_error(" no line image given!" ) ;
02711         return -1 ;
02712     }
02713     ilx=cpl_image_get_size_x(lineImage);
02714     ily=cpl_image_get_size_y(lineImage);
02715     pidata=cpl_image_get_data_float(lineImage);
02716 
02717     slit_length = (int) sqrt (ilx) ;
02718 
02719     if ( NULL == par )
02720     {
02721         sinfo_msg_error(" no line fit parameters given!" ) ;
02722         return -2 ;
02723     }
02724 
02725     if ( NULL == sinfo_slit_pos )
02726     {
02727         sinfo_msg_error(" no position array allocated!" ) ;
02728         return -3 ;
02729     }
02730 
02731     if ( box_length <  4 ||
02732          box_length >= 2*slit_length )
02733     {
02734         sinfo_msg_error(" wrong fitting box length given!" ) ;
02735         sinfo_msg_error(" Must be 4 <= box_length < %d ",2*slit_length ) ;
02736         sinfo_msg_error(" You have chosen box_length = %d ",box_length);
02737         return -4 ;
02738     }
02739 
02740     if ( y_box <= 0.  || y_box > 3. )
02741     {
02742         sinfo_msg_error(" wrong y box length given!" ) ;
02743         sinfo_msg_error(" y_box=%f not in range (0,3]!",y_box);
02744         return -5 ;
02745     }
02746 
02747     if ( diff_tol < 1. )
02748     {
02749         sinfo_msg_error(" diff_tol too small!" ) ;
02750         return -6 ;
02751     }
02752 
02753     /* allocate memory for the edges and the row positon of the slitlets */
02754     sinfo_edge         = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
02755     dummyedge    = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
02756     edgeclean    = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ;
02757     pos_row      = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
02758     pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ;
02759 
02760     /* ------------------------------------------------------------------------
02761      * go through the first image columns and the fit parameters and find 
02762      * the line with the highest intensity 
02763      */
02764     agreed = -1 ;
02765     bad_line = -1 ;
02766     while( agreed == -1 )
02767     {
02768         found = -1 ;
02769         max_intensity = -FLT_MAX ;
02770         for ( col = 0 ; col < slit_length ; col++ )
02771         {
02772             for ( i = 0 ; i < par[0]->n_params ; i++ )
02773             {
02774                 if ( par[i]->column == col && par[i] -> line != bad_line )
02775                 {
02776                     if ( par[i]->fit_par[0] > max_intensity )
02777                     {
02778                         if ( par[i]->fit_par[1] >= 1. && 
02779                              par[i]->fit_par[2] > 0. )
02780                         {
02781                             max_intensity = par[i]->fit_par[0] ;
02782                             found = i ;
02783                         }
02784                     }
02785                 }
02786             }  
02787         }
02788 
02789         /* --------------------------------------------------------------------
02790          * check if the found line is usable and if the neighbouring line 
02791          * have intensity on near rows in neighbouring slitlets 
02792          */
02793         line    = par[found]->line ;
02794         column  = par[found]->column ;
02795         row_pos = par[found]->fit_par[2] ;
02796         if ( found >= 0 && max_intensity > 0. )
02797         {
02798             for ( i = 0 ; i < par[0]->n_params ; i++ )
02799             {
02800                 if ( par[i]->line == line-1 && 
02801                      par[i]->column == column + slit_length )
02802                 {
02803                     if ( par[i]->fit_par[2] <= (row_pos + y_box) &&
02804                          par[i]->fit_par[2] >= (row_pos - y_box) )
02805                     {
02806                         bad_line = line ;
02807                     } 
02808                 }
02809             }
02810             if ( bad_line != line )
02811             {
02812                 agreed = 1 ;
02813                 break ;
02814             }
02815         }
02816         else 
02817         {
02818             sinfo_msg_error("no emission line found in "
02819                             "the first image columns") ;
02820             cpl_free( sinfo_edge ) ;
02821             cpl_free( pos_row ) ;
02822             cpl_free( edgeclean ) ;
02823             cpl_free( dummyedge ) ;
02824             cpl_free( pos_rowclean ) ;
02825             return -7 ;
02826         }    
02827     }
02828 
02829  
02830     if ( agreed == -1 )
02831     {
02832         sinfo_msg_error(" no emission line found in the first image columns") ;
02833         cpl_free( sinfo_edge ) ;
02834         cpl_free( pos_row ) ;
02835         cpl_free( edgeclean ) ;
02836         cpl_free( dummyedge ) ;
02837         cpl_free( pos_rowclean ) ;
02838         return -7 ;
02839     }    
02840  
02841     /* now find and store the raw sinfo_edge positions of the found slitlet */ 
02842     n  = 0 ;
02843     ed = 0 ;
02844     position=cpl_calloc(ilx,sizeof(float)) ;
02845 
02846     for ( col = 0 ; col < ilx ; col++ )
02847     {
02848         for ( i = 0 ; i < par[0]->n_params ; i++ )
02849         {
02850             if ( par[i]->column == col && par[i]->line == line )
02851             {
02852                 if ( par[i]->fit_par[0] > 0. && 
02853                      par[i]->fit_par[1] >= 1.  && 
02854                      par[i]->fit_par[2] > 0. )
02855                 {
02856                     position[n] = par[i]->fit_par[2] ;
02857                     if ( n > 0 && fabs(position[n] - position[n-1]) > y_box )
02858                     {
02859                         sinfo_edge[ed] = col ; 
02860                         pos_row[ed] = sinfo_new_nint( position[n-1] ) ;
02861                         ed++ ;
02862                         if ( col >= ilx - slit_length - 5 ) 
02863                         {
02864                             pos_row[ed] =  sinfo_new_nint( position[n] ) ;
02865                         }
02866                     }
02867                     n++ ;
02868                 }
02869             }
02870         }
02871     }
02872     if ( ed < (slit_length - 1) )
02873     {
02874         sinfo_msg_error(" not enough slitlets found") ;
02875         cpl_free( sinfo_edge ) ;
02876         cpl_free( pos_row ) ;
02877         cpl_free( edgeclean ) ;
02878         cpl_free( dummyedge ) ;
02879         cpl_free( pos_rowclean ) ;
02880         return -8 ;
02881     } 
02882 
02883     /* now find the clean sinfo_edge and row positions of the slitlets */
02884     for ( i = 1 ; i <= ed ; i ++ )
02885     {
02886         if ( i == ed )
02887         {
02888             if ( (sinfo_edge[i-1] - sinfo_edge[i-2]) < slit_length - 3 ||
02889                  (sinfo_edge[i-1] - sinfo_edge[i-2]) > slit_length + 3 )
02890             {
02891                 dummyedge[i-1]   = -1 ;
02892             }
02893             
02894         }
02895         if (dummyedge[i-1] != -1)
02896         {
02897             dummyedge[i-1] = sinfo_edge[i-1] ;
02898         }
02899         else
02900         {
02901             continue ;
02902         }
02903         if ( i < ed )
02904         {
02905             if ( (sinfo_edge[i] - sinfo_edge[i-1]) < slit_length - 3 ||
02906                  (sinfo_edge[i] - sinfo_edge[i-1]) > slit_length + 3 )
02907             {
02908                 dummyedge[i]   = -1 ;
02909             }
02910         }
02911         if ( i + 1 < ed && dummyedge[i] != -1 )
02912         {
02913             if ( (sinfo_edge[i+1] - sinfo_edge[i]) < slit_length - 3 ||
02914                  (sinfo_edge[i+1] - sinfo_edge[i]) > slit_length + 3 )
02915             {
02916                 dummyedge[i+1] = -1 ; 
02917             }
02918         }
02919     }
02920     
02921     k = 0 ;
02922     for ( i = 0 ; i < ed ; i++ )
02923     {
02924         if ( dummyedge[i] != -1 && dummyedge[i] != 0 )
02925         {
02926             edgeclean[k] = dummyedge[i] ;
02927             pos_rowclean[k] = pos_row[i] ;
02928             k++ ;
02929             if( edgeclean[k-1] > (ilx - slit_length - 6 ) )
02930             {
02931                 pos_rowclean[k] = pos_row[ed] ;
02932             }
02933         }
02934     }
02935 
02936     if ( k != slit_length - 1 )
02937     {
02938         sinfo_msg_error(" not enough clean slitlets found") ;
02939         cpl_free( sinfo_edge ) ;
02940         cpl_free( pos_row ) ;
02941         cpl_free( edgeclean ) ;
02942         cpl_free( dummyedge ) ;
02943         cpl_free( pos_rowclean ) ;
02944         return -8 ;
02945     } 
02946 
02947     /* determine the margins of the fitting box outside the slitlets */
02948     margin = box_length / 2 ;
02949 
02950     /* ------------------------------------------------------------------------
02951      * now go through the slitlets, search along each column within a box with 
02952      * half width y_box the maximum value and store these found values in a 
02953      * buffer
02954      */
02955     for ( j = 0 ; j <= k ; j++ )
02956     {
02957         m = 0 ;
02958         if ( j == 0 )
02959         {
02960             box_buffer = sinfo_new_vector( edgeclean[0] + margin ) ;
02961             for( col = 0 ; col < edgeclean[0] + margin ; col++ )
02962             {
02963                 maxval = -FLT_MAX ;
02964                 for ( row = pos_rowclean[0] - sinfo_new_nint(y_box) ; 
02965                       row <= pos_rowclean[0] + sinfo_new_nint(y_box) ; row++ )
02966                 {
02967                     if ( maxval < pidata[col + ilx*row] )
02968                     {
02969                         maxval = pidata[col + ilx*row] ;
02970                     }
02971                 }
02972                 box_buffer->data[m] = maxval ;
02973                 m++ ;
02974             }
02975         }
02976         else if ( j < k )
02977         {
02978             box_buffer = sinfo_new_vector( edgeclean[j] - 
02979                                            edgeclean[j-1] + 2*margin ) ;
02980             for ( col = edgeclean[j - 1] - margin ; 
02981                   col < edgeclean[j] + margin ; col++ )
02982             {
02983                 maxval = -FLT_MAX ;
02984                 for ( row = pos_rowclean[j] - sinfo_new_nint(y_box) ; 
02985                       row <= pos_rowclean[j] + sinfo_new_nint(y_box) ; row++ )
02986                 {
02987                     if ( maxval < pidata[col + ilx*row] )
02988                     {
02989                         maxval = pidata[col + ilx*row] ;
02990                     }
02991                 }
02992                 box_buffer->data[m] = maxval ;
02993                 m++ ;
02994             }
02995         }
02996         else 
02997         {
02998             box_buffer = sinfo_new_vector( ilx - edgeclean[k-1] + margin ) ;
02999             for ( col = edgeclean[k - 1] - margin ; col < ilx ; col++ )
03000             {
03001                 maxval = -FLT_MAX ;
03002                 for ( row = pos_rowclean[k] - sinfo_new_nint(y_box) ; 
03003                       row <= pos_rowclean[k] + sinfo_new_nint(y_box) ; row++ )
03004                 {
03005                     if ( maxval < pidata[col + ilx*row] )
03006                     {
03007                         maxval = pidata[col + ilx*row] ;
03008                     }
03009                 }
03010                 box_buffer->data[m] = maxval ;
03011                 m++ ;
03012             }
03013         }
03014 
03015         for ( left_right = 0 ; left_right <= 1 ; left_right++ )
03016         { 
03017             nel = 0 ;
03018             if ( left_right == 0 )
03019             {
03020                 nel = box_buffer -> n_elements / 2 ;
03021             }
03022             else
03023             {
03024                 if ( box_buffer -> n_elements % 2 == 0 )
03025                 {
03026                     nel = box_buffer -> n_elements / 2 ;
03027                 }
03028                 else
03029                 {
03030                     nel = box_buffer -> n_elements / 2 + 1 ;
03031                 }
03032             }
03033 
03034             /* now split the buffer in the midth in a left and right 
03035                part for fitting */
03036             half_buffer = sinfo_new_vector( nel ) ;
03037             if ( left_right == 0 )
03038             {
03039                 for ( i = 0 ; i < nel ; i++ )
03040                 {
03041                     half_buffer -> data[i] = box_buffer -> data[i] ;
03042                 }
03043             }
03044             else
03045             {
03046                 n_right = 0 ;
03047                 for ( i = box_buffer -> n_elements - 1 ; 
03048                       i >= box_buffer -> n_elements - nel ; i-- )
03049                 {
03050                     half_buffer -> data[n_right] = box_buffer -> data[i] ;
03051                     n_right++ ;
03052                 }
03053             }
03054 
03055             xdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
03056             wdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
03057             mpar = (int *)   cpl_calloc( NPAR, sizeof (int) ) ;
03058 
03059             /* set initial values for the fitting routine */
03060             minval =  FLT_MAX ;
03061             maxval = -FLT_MAX ;
03062             for ( i = 0 ; i < nel ; i++ )
03063             {
03064                 xdat[i] = i ;
03065                 wdat[i] = 1.0 ;
03066                 if ( half_buffer -> data[i] < minval )
03067                 {
03068                     minval = half_buffer -> data[i] ;
03069                 }
03070                 if ( half_buffer -> data[i] > maxval )
03071                 {
03072                     maxval = half_buffer -> data[i] ;
03073                 }
03074             }
03075 
03076             fitpar[2] = minval ;
03077             fitpar[3] = maxval ; 
03078 
03079             /* search for both positions of the half intensity of 
03080                the hat within the buffer */
03081             init1 = -1 ; 
03082             for ( i = 0 ; i < nel ; i++ )
03083             {
03084                 if ( half_buffer -> data[i] >= ( maxval + minval ) / 2. )
03085                 {
03086                     init1 = i ;
03087                     break ;
03088                 }
03089             }
03090 
03091             /* determine the initial positions from the found values */
03092             if ( init1 != -1 )
03093             {
03094                 fitpar[0] = ((float)init1 - 1.) ;
03095                 fitpar[1] = ((float)init1 + 1.) ;
03096             }
03097 
03098             for ( i = 0 ; i < NPAR ; i++ )
03099             {
03100                 mpar[i] = 1 ;
03101                 dervpar[i] = 0. ;
03102             }
03103       
03104             xdim     = XDIMA ;
03105             ndat     = nel ;
03106             numpar   = NPAR ;
03107             tol      = TOLA ;
03108             lab      = LABA ;
03109             its      = ITSA ;
03110          
03111             /* finally, do the least squares fit over the buffer data */
03112             if ( 0 > ( iters = sinfo_new_lsqfit_edge( xdat, &xdim, 
03113                                                       half_buffer -> data, 
03114                                                       wdat, &ndat, fitpar,
03115                                                       dervpar, mpar, &numpar, 
03116                                                       &tol, &its, &lab )) )
03117             { 
03118                 /* if the fit doesn't succeed the initial values are taken */
03119                 sinfo_msg_warning ("least squares fit failed, error "
03120                                    "no.: %d in slitlet: %d", iters, j) ;
03121                 fitpar[0] = ((float)init1 - 1.) ;
03122                 fitpar[1] = ((float)init1 + 1.) ;
03123             }
03124 
03125             pos = (fitpar[0] + fitpar[1]) / 2. ;
03126 
03127             /*----------------------------------------------------------------- 
03128              * now discern the left and the right sinfo_edge fit of the 
03129              * slitlets and associate the fit results with the absolute 
03130              * positions in the image consider the difference of the fitted 
03131              * slit position to the expected position and decide wether the 
03132              * fit is taken or the expected value is taken.
03133              */
03134             if ( left_right == 0 )
03135             {
03136                 /* take care of the column position of the fit boxes to get 
03137                    the absolute positions */
03138                 if ( j == 0 )
03139                 {
03140                     if ( fabs(pos - ((float)edgeclean[0] - 1. - 
03141                          (float)slit_length)) < diff_tol )
03142                     {
03143                         sinfo_slit_pos[0][0] = pos ;
03144                     }
03145                     else
03146                     {
03147                         sinfo_msg_warning("something wrong with fitted "
03148                                           "left position of slitlet 0") ;
03149                         if ((float) edgeclean[0] - 1. - (float)slit_length < 0. )
03150                         {
03151                             sinfo_slit_pos[0][0] = 0. ;
03152                         }
03153                         else
03154                         {
03155                             sinfo_slit_pos[0][0] = (float)edgeclean[0] - 1. - 
03156                                                    (float)slit_length ;
03157                         }
03158                     }
03159                 }
03160                 else if ( j < k )
03161                 {
03162                     if ( fabs( pos - (float)margin ) < diff_tol )
03163                     {
03164                         sinfo_slit_pos[j][0] = pos + (float)edgeclean[j-1] - 
03165                                                      (float)margin ;
03166                     }
03167                     else
03168                     {
03169                         sinfo_msg_warning("something wrong with fitted "
03170                                           "left position of slitlet %d", j) ;
03171                         sinfo_slit_pos[j][0] = (float)edgeclean[j-1] - 1. ;
03172                     }
03173                 }
03174                 else
03175                 {
03176                     if ( fabs( pos - (float)margin ) < diff_tol )
03177                     {
03178                         sinfo_slit_pos[k][0] = pos + (float)edgeclean[k-1] - 
03179                                                      (float)margin ;
03180                     }
03181                     else
03182                     {
03183                         sinfo_msg_warning("something wrong with fitted left "
03184                                           "position of slitlet %d", j) ;
03185                         sinfo_slit_pos[k][0] = (float)edgeclean[k-1] - 1. ;
03186                     }
03187                 }
03188             }
03189             else
03190             {
03191                 /* take care of the column position of the fit boxes to 
03192                    get the absolute positions */
03193                 if ( j == 0 )
03194                 {
03195                     if ( fabs( (float)box_buffer->n_elements - pos - 
03196                                (float)edgeclean[0] ) < diff_tol )
03197                     {
03198                         sinfo_slit_pos[0][1] = (float)(box_buffer->n_elements - 
03199                                                1) - pos ;
03200                     }
03201                     else
03202                     {
03203                         sinfo_msg_warning("something wrong with fitted "
03204                                           "right position of slitlet 0") ;
03205                         sinfo_slit_pos[0][1] = (float)edgeclean[0] - 1. ;
03206                     }
03207                 }
03208                 else if ( j < k )
03209                 {
03210                     if ( fabs( (float)box_buffer->n_elements - pos
03211                              + (float)edgeclean[j-1] - (float)margin - 
03212                                (float)edgeclean[j] ) < diff_tol )
03213                     {
03214                         sinfo_slit_pos[j][1] = (float)(box_buffer->n_elements -
03215                                                1) - pos
03216                                        + (float)edgeclean[j-1] - (float)margin;
03217                     }
03218                     else
03219                     {
03220                         sinfo_msg_warning("something wrong with fitted "
03221                                           "right position of slitlet %d", j) ;
03222                         sinfo_slit_pos[j][1] = (float)edgeclean[j] - 1. ;
03223                     }
03224                 }
03225                 else
03226                 {
03227                     if ( edgeclean[k-1] + slit_length > ilx )
03228                     {
03229                         last_pos = (float)(ilx - 1) ;
03230                     }
03231                     else
03232                     {
03233                         last_pos = (float)(edgeclean[k-1] - 1 + slit_length) ;
03234                     }
03235                     if ( fabs( (float)(box_buffer->n_elements - 1) - pos
03236                              + (float)edgeclean[k-1] - (float)margin - 
03237                                last_pos ) < diff_tol )
03238                     {
03239                         sinfo_slit_pos[k][1] = (float)(box_buffer->n_elements -
03240                                                1) - pos
03241                                        + (float)edgeclean[k-1] - (float)margin ;
03242                     }
03243                     else
03244                     {
03245                         sinfo_msg_warning("something wrong with fitted "
03246                                           "right position of slitlet %d", j) ;
03247                         sinfo_slit_pos[k][1] = last_pos ;
03248                     }
03249                 }
03250             }
03251 
03252             sinfo_new_destroy_vector ( half_buffer ) ;
03253             cpl_free( xdat ) ;
03254             cpl_free( wdat ) ;
03255             cpl_free( mpar ) ;
03256         }
03257         sinfo_new_destroy_vector ( box_buffer ) ;
03258     }
03259         
03260     cpl_free( sinfo_edge ) ;
03261     cpl_free( pos_row ) ;
03262     cpl_free( edgeclean ) ;
03263     cpl_free( dummyedge ) ;
03264     cpl_free( pos_rowclean ) ;
03265     cpl_free( position );
03266     return 0 ;
03267 }
03268 
03291 int 
03292 sinfo_new_fit_slits_edge_with_estimate ( cpl_image   * lineImage,
03293                                float     ** sinfo_slit_pos,
03294                                int          box_length,
03295                                float        y_box,
03296                                float        diff_tol,
03297                                int          low_pos,
03298                                int          high_pos )
03299 {
03300     int*   position=NULL ;
03301     Vector * box_buffer ;
03302     Vector * in_buffer ;
03303     int   found_row ;
03304     int   row, col ;
03305     int   col_first, col_last ;
03306     int   row_first, row_last ;
03307     int   i, j, m, n ;
03308     int   init1 ;
03309     int   left_right ;
03310     int   n_buf, shift ;
03311     int   slit_length ;
03312     int   iters, xdim, ndat ;
03313     int   numpar, its ;
03314     int   * mpar ;
03315     float * xdat, * wdat ;
03316     float tol, lab ;
03317     float fitpar[NPAR] ;
03318     float dervpar[NPAR] ;
03319     float minval, maxval ;
03320     float pos ;
03321     float new_pos ;
03322     int   slitposition[SLITLENGTH] ;
03323     pixelvalue rowpos[SLITLENGTH] ;
03324     int ilx=0;
03325     int ily=0;
03326     float* pidata=NULL;
03327 
03328     slit_length = SLITLENGTH ; /* this is too high: 64 */
03329     slit_length = N_SLITLETS ; /* this is better: 32 */
03330 
03331 
03332     if ( NULL == lineImage )
03333     {
03334         sinfo_msg_error(" no line image given!" ) ;
03335         return -1 ;
03336     }
03337     ilx=cpl_image_get_size_x(lineImage);
03338     ily=cpl_image_get_size_y(lineImage);
03339     pidata=cpl_image_get_data_float(lineImage);
03340 
03341     if ( NULL == sinfo_slit_pos )
03342     {
03343         sinfo_msg_error(" no position array allocated!" ) ;
03344         return -1 ;
03345     }
03346 
03347     if ( box_length < 4 ||
03348          box_length > 2*slit_length )
03349     {
03350         sinfo_msg_error(" wrong fitting box length given!" ) ;
03351         sinfo_msg_error(" Must be 4 <= box_length < %d ",2*slit_length ) ;
03352         sinfo_msg_error(" You have chosen box_length = %d",box_length);
03353 
03354 
03355         return -1 ;
03356     }
03357 
03358     if ( y_box <= 0. || y_box > 6. )
03359     {
03360         sinfo_msg_error("wrong y box length given!" ) ;
03361         sinfo_msg_error("You have chosen y_box=%f not in range (0,6]!",y_box);
03362         return -1 ;
03363     }
03364     if ( diff_tol <= 0. )
03365     {
03366         sinfo_msg_error(" wrong diff_tol given!" ) ;
03367         return -1 ;
03368     }
03369 
03370         if ( low_pos >= high_pos || low_pos < 0 || 
03371              high_pos <= 0 || high_pos > ily )
03372     {
03373         sinfo_msg_error(" wrong user given search positions!" ) ;
03374         return -1 ;
03375     }
03376 
03377     /* now search for the maximum between the given positions for each col */
03378     position=cpl_calloc(ilx,sizeof(int)) ;
03379 
03380     for ( col = 0 ; col < ilx ; col++ )
03381     {
03382         found_row = -1 ;
03383         maxval = -FLT_MAX ;
03384         for ( row = low_pos ; row <= high_pos ; row++ )
03385         {
03386             if ( maxval < pidata[col+row*ilx] )
03387             {
03388                 maxval = pidata[col+row*ilx] ;
03389                 found_row = row ;
03390             }
03391         }
03392         if ( maxval > -FLT_MAX && found_row > low_pos )
03393         {
03394             position[col] = found_row ;
03395         }
03396         else
03397         {
03398             position[col] = 0 ;
03399         }
03400     }
03401 
03402     /* ------------------------------------------------------------------------
03403      * now go through the slitlets, search along each column within a box with
03404      * half width y_box the maximum value and store these found values in a 
03405      * buffer
03406      */
03407     for ( j = 0 ; j < slit_length ; j++ )
03408     {
03409         /* now go through the columns and determine the slitlet positions by
03410          * calculating the sinfo_median of the found positions
03411          */
03412         n = 0 ;
03413         for ( col = sinfo_new_nint(sinfo_slit_pos[j][0])+ 1 ; 
03414               col < sinfo_new_nint(sinfo_slit_pos[j][1]) -1 ; col++ )
03415         {
03416             rowpos[n] = (pixelvalue)position[col] ;
03417             n++ ;
03418         }
03419         slitposition[j] = (int)sinfo_new_median(rowpos, n) ;
03420         for ( left_right = 0 ; left_right <= 1 ; left_right++ )
03421 
03422         {
03423             init1 = 0 ;
03424             col_first = sinfo_new_nint( sinfo_slit_pos[j][left_right] ) - 
03425                                         box_length/2 ;
03426             col_last  = sinfo_new_nint( sinfo_slit_pos[j][left_right] ) + 
03427                         box_length/2 ;
03428             if ( col_first < 0 )
03429             {
03430                 col_first = 0 ;
03431                 init1 = 1 ;
03432             }
03433             if ( col_last > ilx )
03434             {
03435                 col_last = ilx ;
03436                 init1 = 1 ;
03437             }
03438             if ( col_last - col_first <= 0 )
03439             {
03440                 sinfo_msg_error(" first and last column wrong!" ) ;
03441                 return -1 ;
03442             }
03443             box_buffer = sinfo_new_vector( col_last - col_first ) ;
03444             m = 0 ;
03445             if ( left_right == 0 )
03446             {
03447                 for( col = col_first ; col < col_last ; col++ )
03448                 {
03449                     row_first = slitposition[j] - sinfo_new_nint(y_box) ;
03450                     row_last  = slitposition[j] + sinfo_new_nint(y_box) ;
03451                     if ( row_first < 0 )
03452                     {
03453                         row_first = 0 ;
03454                     }
03455                     if ( row_last >= ily  )
03456                     {
03457                         row_last = ily - 1 ;
03458                     }
03459                     maxval = -FLT_MAX ;
03460                     for ( row = row_first ; row <= row_last ; row++ )
03461                     {
03462                         if ( maxval < pidata[col + ilx*row] )
03463                         {
03464                             maxval = pidata[col + ilx*row] ;
03465                         }
03466                     }
03467                     box_buffer->data[m] = maxval ;
03468                     m++ ;
03469                 }
03470             }
03471             else
03472             {
03473                 for( col = col_last-1 ; col >= col_first ; col-- )
03474                 {
03475                     row_first = slitposition[j] - sinfo_new_nint(y_box) ;
03476                     row_last  = slitposition[j] + sinfo_new_nint(y_box) ;
03477                     if ( row_first < 0 )
03478                     {
03479                         row_first = 0 ;
03480                     }
03481                     if ( row_last >= ily  )
03482                     {
03483                         row_last = ily - 1 ;
03484                     }
03485                     maxval = -FLT_MAX ;
03486                     for ( row = row_first ; row <= row_last ; row++ )
03487                     {
03488                         if ( maxval < pidata[col + ilx*row] )
03489                         {
03490                             maxval = pidata[col + ilx*row] ;
03491                         }
03492                     }
03493                     box_buffer->data[m] = maxval ;
03494                     m++ ;
03495                 }
03496             }
03497 
03498             xdat=(float *)cpl_calloc( box_buffer->n_elements, sizeof (float));
03499             wdat=(float *)cpl_calloc( box_buffer->n_elements, sizeof (float));
03500             mpar=(int *)  cpl_calloc( NPAR, sizeof (int) ) ;
03501 
03502             /* set initial values for the fitting routine */
03503             minval =  FLT_MAX ;
03504             maxval = -FLT_MAX ;
03505             for ( i = 0 ; i < box_buffer->n_elements ; i++ )
03506             {
03507                 xdat[i] = i ;
03508                 wdat[i] = 1.0 ;
03509                 if ( box_buffer -> data[i] < minval )
03510                 {
03511                     minval = box_buffer -> data[i] ;
03512                 }
03513                 if ( box_buffer -> data[i] > maxval )
03514                 {
03515                     maxval = box_buffer -> data[i] ;
03516                 }
03517             }
03518             fitpar[2] = minval ;
03519             fitpar[3] = maxval ;
03520             /*----------------------------------------------------------------
03521              * if we have too few left background values (at the image edges)
03522              * the left margin of the buffer to fit is filled with the minimal
03523              * values in order to get a good fit
03524              */
03525 
03526             if ( init1 == 1 )
03527             {
03528                 n_buf = box_buffer->n_elements + box_length/2 ;
03529                 in_buffer = sinfo_new_vector( n_buf ) ;
03530                 for ( i = 0 ; i < box_length/2 ; i++ )
03531                 {
03532                     in_buffer -> data[i] = minval ;
03533                 }
03534                 shift = 0 ;
03535                 for ( i = box_length/2 ; i < n_buf ; i++ )
03536                 {
03537                     in_buffer -> data[i] = box_buffer -> data[shift] ;
03538                     shift++ ;
03539                 }
03540                 sinfo_new_destroy_vector ( box_buffer ) ;
03541                 box_buffer = sinfo_new_vector ( n_buf ) ;
03542                 for ( i = 0 ; i < n_buf ; i++ )
03543                 {
03544                     box_buffer -> data[i] = in_buffer -> data[i] ;
03545                 }
03546                 sinfo_new_destroy_vector ( in_buffer ) ;
03547             }
03548             /* determine the initial positions from the found values */
03549             fitpar[0] = (float)box_buffer->n_elements/2. - 1. ;
03550             fitpar[1] = (float)box_buffer->n_elements/2. + 1. ;
03551 
03552             for ( i = 0 ; i < NPAR ; i++ )
03553             {
03554                 mpar[i] = 1 ;
03555                 dervpar[i] = 0. ;
03556             }
03557 
03558             xdim     = XDIMA ;
03559             ndat     = box_buffer->n_elements ;
03560             numpar   = NPAR ;
03561             tol      = TOLA ;
03562             lab      = LABA ;
03563             its      = ITSA ;
03564 
03565             /* finally, do the least squares fit over the buffer data */
03566             if ( 0 > ( iters = sinfo_new_lsqfit_edge( xdat, &xdim, 
03567                                                       box_buffer -> data, 
03568                                                       wdat, &ndat, fitpar,
03569                                                       dervpar, mpar, &numpar, 
03570                                                       &tol, &its, &lab )) )
03571             {
03572                 sinfo_msg_warning (" least squares fit failed, error "
03573                                    "no.: %d in slitlet: %d\n", iters, j) ;
03574                 sinfo_new_destroy_vector(box_buffer) ;
03575                 cpl_free( xdat ) ;
03576                 cpl_free( wdat ) ;
03577                 cpl_free( mpar ) ;
03578                 continue ;
03579             }
03580             if ( fitpar[1] <= fitpar[0] )
03581             {
03582                 sinfo_msg_warning ("fit failed due to negative slope of "
03583                                    "sinfo_new_edge function in slitlet: %d",j);
03584                 sinfo_new_destroy_vector(box_buffer) ;
03585                 cpl_free( xdat ) ;
03586                 cpl_free( wdat ) ;
03587                 cpl_free( mpar ) ;
03588                 continue ;
03589             }
03590 
03591             pos = (fitpar[0] + fitpar[1])/2. ;
03592             if ( init1 == 1 )
03593             {
03594                 pos -= (float)box_length/2. ;
03595             }
03596 
03597             /*-------------------------------------------------------------
03598              * now compute the real slit positions using the guess positions
03599              * if the fit did not work the guess positions are taken
03600              * the same is done if the deviations are too big.
03601              */
03602             if ( pos != 0. )
03603             {
03604                 if ( left_right == 0 )
03605                 {
03606                     new_pos = (float)col_first + pos ;
03607                 }
03608                 else
03609                 {
03610                     new_pos = (float)col_last-1 - pos ;
03611                 }
03612                 if ( fabs(new_pos - sinfo_slit_pos[j][left_right]) < diff_tol )
03613                 {
03614                     sinfo_slit_pos[j][left_right] = new_pos ;
03615                 }
03616                 else
03617                 {
03618                     sinfo_msg_warning (" deviation bigger than tolerance,"
03619                                        " take the estimated slitlet positiona"
03620                                        " in slitlet: %d\n", j) ;
03621                 }
03622             }
03623 
03624             cpl_free( xdat ) ;
03625             cpl_free( wdat ) ;
03626             cpl_free( mpar ) ;
03627             sinfo_new_destroy_vector ( box_buffer ) ;
03628         }
03629     }
03630     cpl_free(position);
03631     return 0 ;
03632 }
03633 
03634 /*--------------------------------------------------------------------------*/

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