00001 /*******************************************************************************
00002 * E.S.O. - VLT project
00003 *
00004 *
00005 *
00006 * who when what
00007 * -------- -------- ----------------------------------------------
00008 * schreib 27/02/01 created
00009 */
00010
00011 /************************************************************************
00012 * NAME
00013 * boltzmann.c -
00014 * routines to determine the absolute positions of the slitlets out
00015 * of an emission line frame
00016 *
00017 * SYNOPSIS
00018 * #include "absolute.h"
00019 *
00020 * 1) float boltz ( float * xdat, float * parlist )
00021 * 2) void boltz_deriv( float * xdat, float * parlist, float * dervs )
00022 * 3) static int inv_mat (void)
00023 * 4) static void get_mat ( float * xdat,
00024 * int * xdim,
00025 * float * ydat,
00026 * float * wdat,
00027 * int * ndat,
00028 * float * fpar,
00029 * float * epar,
00030 * int * npar )
00031 * 5) static int get_vec ( float * xdat,
00032 * int * xdim,
00033 * float * ydat,
00034 * float * wdat,
00035 * int * ndat,
00036 * float * fpar,
00037 * float * epar,
00038 * int * npar )
00039 * 6) int lsqfit ( float * xdat,
00040 * int * xdim,
00041 * float * ydat,
00042 * float * wdat,
00043 * int * ndat,
00044 * float * fpar,
00045 * float * epar,
00046 * int * mpar,
00047 * int * npar,
00048 * float * tol ,
00049 * int * its ,
00050 * float * lab )
00051 * 7) int fitSlitsBoltz( OneImage * lineImage,
00052 * FitParams ** par,
00053 * float ** slit_pos,
00054 * int box_length,
00055 * float y_box,
00056 * float diff_tol )
00057 * 8) int fitSlitsBoltzSingleLine ( OneImage * lineImage,
00058 * float ** slit_pos,
00059 * int box_length,
00060 * float y_box,
00061 * int low_pos,
00062 * int high_pos )
00063 * 9) int fitSlitsBoltzWithEstimate ( OneImage * lineImage,
00064 * float ** slit_pos,
00065 * int box_length,
00066 * float y_box,
00067 * float diff_tol,
00068 * int low_pos,
00069 * int high_pos )
00070 *
00071 * DESCRIPTION
00072 * 1) calculates the value of a Boltzmann function with parameters
00073 * parlist at the position xdat
00074 * 2) calculates the partial derivatives for a Boltzmann function with
00075 * parameters parlist at position xdat
00076 * 3) calculates the inverse of matrix2. The algorithm used
00077 * is the Gauss-Jordan algorithm described in Stoer,
00078 * Numerische Mathematik, 1. Teil.
00079 * 4) builds the matrix
00080 * 5) calculates the correction vector. The matrix has been
00081 * built by get_mat(), we only have to rescale it for the
00082 * current value of labda. The matrix is rescaled so that
00083 * the diagonal gets the value 1 + labda.
00084 * Next we calculate the inverse of the matrix and then
00085 * the correction vector.
00086 * 6) this is a routine for making a least-squares fit of a
00087 * function to a set of data points. The method used is
00088 * described in: Marquardt, J.Soc.Ind.Appl.Math. 11. 431 (1963).
00089 * This method is a mixture of the steepest descent method
00090 * and the Taylor method.
00091 * 7) fits the beginning and end position of the slitlets
00092 * by using non-linear least square fitting of Boltzmann function
00093 * fits a Boltzmann function to the slitlet edges exposed and indicated
00094 * by the brightest emission lines. To achieve this, the fit
00095 * parameters are used to find the brightest emission line
00096 * and to get its position for each column.
00097 * The least squares fit is done by using a box smaller than
00098 * the size of two slitlets
00099 * 8) fits the beginning and end position of the slitlets
00100 * by using non-linear least square fitting of a Boltzmann function
00101 * fits a Boltzmann function to the slitlet edges exposed and indicated
00102 * by the brightest emission lines. The slitlet is searched within
00103 * user given positions.
00104 * The least squares fit is done by using a box smaller than
00105 * the size of two slitlets
00106 * 9) fits the beginning and end position of the slitlets
00107 * by using non-linear least square fitting of a Boltzmann function
00108 * fits a Boltzmann function to the slitlet edges exposed and indicated
00109 * by the brightest emission lines. The slitlet is searched within
00110 * user given positions.
00111 * The least squares fit is done by using a box smaller than
00112 * the size of two slitlets
00113 *
00114 * FILES
00115 *
00116 * ENVIRONMENT
00117 *
00118 * RETURN VALUES
00119 *
00120 * CAUTIONS
00121 *
00122 * EXAMPLES
00123 *
00124 * SEE ALSO
00125 *
00126 * BUGS
00127 *
00128 *------------------------------------------------------------------------
00129 */
00130
00131 #include "vltPort.h"
00132
00133 /*
00134 * System Headers
00135 */
00136
00137 /*
00138 * Local Headers
00139 */
00140
00141 #include "absolute.h"
00142 #include "recipes.h"
00143
00144 /*----------------------------------------------------------------------------
00145 * Defines
00146 *--------------------------------------------------------------------------*/
00147
00148 #define XDIMA 1 /* dimension of the x values */
00149 #define TOLA 0.001 /* fitting tolerance */
00150 #define LABA 0.1 /* labda parameter */
00151 #define ITSA 200 /* maximum number of iterations */
00152 #define LABFACA 10.0 /* labda step factor */
00153 #define LABMAXA 1.0e+10 /* maximum value for labda */
00154 #define LABMINA 1.0e-10 /* minimum value for labda */
00155 #define NPAR 4 /* number of fit parameters */
00156
00157 /*----------------------------------------------------------------------------
00158 * Local variables
00159 *--------------------------------------------------------------------------*/
00160
00161 static double chi1 ; /* old reduced chi-squared */
00162 static double chi2 ; /* new reduced chi-squared */
00163 static double labda ; /* mixing parameter */
00164 static double vec[NPAR] ; /* correction vector */
00165 static double matrix1[NPAR][NPAR] ; /* original matrix */
00166 static double matrix2[NPAR][NPAR] ; /* inverse of matrix1 */
00167 static int nfree ; /* number of free parameters */
00168 static int parptr[NPAR] ; /* parameter pointer */
00169
00170 /*----------------------------------------------------------------------------
00171 * Functions private to this module
00172 *--------------------------------------------------------------------------*/
00173
00174 static int inv_mat (void) ;
00175
00176 static void get_mat ( float * xdat,
00177 int * xdim,
00178 float * ydat,
00179 float * wdat,
00180 int * ndat,
00181 float * fpar,
00182 float * epar/*,
00183 int * npar*/ ) ;
00184
00185 static int get_vec ( float * xdat,
00186 int * xdim,
00187 float * ydat,
00188 float * wdat,
00189 int * ndat,
00190 float * fpar,
00191 float * epar,
00192 int * npar ) ;
00193
00194 /*----------------------------------------------------------------------------
00195 * Function codes
00196 *--------------------------------------------------------------------------*/
00197
00198 /*----------------------------------------------------------------------------
00199 Function : boltz()
00200 In : position array xdat, parameter list parlist
00201 The parameters are:
00202 parlist(0): background1
00203 parlist(1): background2
00204 parlist(2): central position
00205 parlist(3): width
00206 Out : function value of a Boltzmann function
00207 that is
00208 y = (parlist(0) - parlist(1)) / (1+exp((x-parlist(2))/parlist(3))) + parlist(1)
00209 Job : calculates the value of a Boltzmann function with parameters
00210 parlist at the position xdat
00211 ---------------------------------------------------------------------------*/
00212
00213 float boltz ( float * xdat, float * parlist )
00214 {
00215 float return_value ;
00216
00217 /* now build the boltzman function out of the parameters */
00218 return_value =
00219 (parlist[0] - parlist[1]) / (1 + exp(( xdat[0] - parlist[2] ) / parlist[3])) + parlist[1] ;
00220
00221 return return_value ;
00222 }
00223
00224 /*----------------------------------------------------------------------------
00225 Function : boltz_deriv()
00226 In : position array xdat, parameter list parlist
00227 The parameters are:
00228 parlist(0): background1
00229 parlist(1): background2
00230 parlist(2): central position
00231 parlist(3): width
00232 derivative value of a Boltzmann function at position xdat: dervs
00233 dervs[0]: partial derivative by background1
00234 dervs[1]: partial derivative by background2
00235 dervs[2]: partial derivative by central position
00236 dervs[3]: partial derivative by the width
00237 Out : nothing, void
00238 Job : calculates the partial derivatives for a Boltzmann function with
00239 parameters parlist at position xdat
00240 ----------------------------------------------------------------------------*/
00241
00242 void boltz_deriv( float * xdat, float * parlist, float * dervs )
00243 {
00244 float subst ;
00245
00246 subst = (xdat[0] - parlist[2]) / parlist[3] ;
00247
00248 dervs[0] = 1. / ( 1. + exp(subst) ) ;
00249
00250 dervs[1] = -1. / ( 1. + exp(subst) ) + 1. ;
00251
00252 dervs[2] = ( (parlist[0] - parlist[1]) / parlist[3] * exp(subst) ) /
00253 ( (1. + exp(subst)) * (1. + exp(subst)) ) ;
00254
00255 dervs[3] = ( (parlist[0] - parlist[1]) * (xdat[0] - parlist[2]) /
00256 (parlist[3]*parlist[3]) * exp(subst) ) /
00257 ( (1. + exp(subst)) * (1. + exp(subst)) ) ;
00258 }
00259
00260 /*----------------------------------------------------------------------------
00261 Function : inv_mat()
00262 In : void
00263 Out : integer (0 if it worked, -6 if determinant is zero)
00264 Job : calculates the inverse of matrix2. The algorithm used
00265 is the Gauss-Jordan algorithm described in Stoer,
00266 Numerische Mathematik, 1. Teil.
00267 ---------------------------------------------------------------------------*/
00268
00269 static int inv_mat (void)
00270 {
00271 double even ;
00272 double hv[NPAR] ;
00273 double mjk ;
00274 double rowmax ;
00275 int evin ;
00276 int i, j, k, row ;
00277 int per[NPAR] ;
00278
00279 /* set permutation array */
00280 for ( i = 0 ; i < nfree ; i++ )
00281 {
00282 per[i] = i ;
00283 }
00284
00285 for ( j = 0 ; j < nfree ; j++ ) /* in j-th column */
00286 {
00287 /* determine largest element of a row */
00288 rowmax = fabs ( matrix2[j][j] ) ;
00289 row = j ;
00290
00291 for ( i = j + 1 ; i < nfree ; i++ )
00292 {
00293 if ( fabs ( matrix2[i][j] ) > rowmax )
00294 {
00295 rowmax = fabs( matrix2[i][j] ) ;
00296 row = i ;
00297 }
00298 }
00299
00300 /* determinant is zero! */
00301 if ( matrix2[row][j] == 0.0 )
00302 {
00303 return -6 ;
00304 }
00305
00306 /* if the largest element is not on the diagonal, then permutate rows */
00307 if ( row > j )
00308 {
00309 for ( k = 0 ; k < nfree ; k++ )
00310 {
00311 even = matrix2[j][k] ;
00312 matrix2[j][k] = matrix2[row][k] ;
00313 matrix2[row][k] = even ;
00314 }
00315 /* keep track of permutation */
00316 evin = per[j] ;
00317 per[j] = per[row] ;
00318 per[row] = evin ;
00319 }
00320
00321 /* modify column */
00322 even = 1.0 / matrix2[j][j] ;
00323 for ( i = 0 ; i < nfree ; i++ )
00324 {
00325 matrix2[i][j] *= even ;
00326 }
00327 matrix2[j][j] = even ;
00328
00329 for ( k = 0 ; k < j ; k++ )
00330 {
00331 mjk = matrix2[j][k] ;
00332 for ( i = 0 ; i < j ; i++ )
00333 {
00334 matrix2[i][k] -= matrix2[i][j] * mjk ;
00335 }
00336 for ( i = j + 1 ; i < nfree ; i++ )
00337 {
00338 matrix2[i][k] -= matrix2[i][j] * mjk ;
00339 }
00340 matrix2[j][k] = -even * mjk ;
00341 }
00342
00343 for ( k = j + 1 ; k < nfree ; k++ )
00344 {
00345 mjk = matrix2[j][k] ;
00346 for ( i = 0 ; i < j ; i++ )
00347 {
00348 matrix2[i][k] -= matrix2[i][j] * mjk ;
00349 }
00350 for ( i = j + 1 ; i < nfree ; i++ )
00351 {
00352 matrix2[i][k] -= matrix2[i][j] * mjk ;
00353 }
00354 matrix2[j][k] = -even * mjk ;
00355 }
00356 }
00357
00358 /* finally, repermute the columns */
00359 for ( i = 0 ; i < nfree ; i++ )
00360 {
00361 for ( k = 0 ; k < nfree ; k++ )
00362 {
00363 hv[per[k]] = matrix2[i][k] ;
00364 }
00365 for ( k = 0 ; k < nfree ; k++ )
00366 {
00367 matrix2[i][k] = hv[k] ;
00368 }
00369 }
00370
00371 /* all is well */
00372 return 0 ;
00373 }
00374
00375 /*----------------------------------------------------------------------------
00376 Function : get_mat()
00377 In : xdat: position
00378 xdim: factor of the indices of the position array
00379 ydat: real data
00380 wdat: weights
00381 ndat: number of data points
00382 fpar: function parameters
00383 epar: partial derivatives of the function
00384 npar: number of function parameters
00385 Out : void
00386 Job : builds the matrix
00387 ---------------------------------------------------------------------------*/
00388
00389 static void get_mat ( float * xdat,
00390 int * xdim,
00391 float * ydat,
00392 float * wdat,
00393 int * ndat,
00394 float * fpar,
00395 float * epar/*,
00396 int * npar*/ )
00397 {
00398 double wd ;
00399 double wn ;
00400 double yd ;
00401 int i, j, n ;
00402
00403 for ( j = 0 ; j < nfree ; j++ )
00404 {
00405 vec[j] = 0.0 ; /* zero vector */
00406 for ( i = 0 ; i<= j ; i++ ) /* zero matrix only on and below diagonal */
00407 {
00408 matrix1[j][i] = 0.0 ;
00409 }
00410 }
00411 chi2 = 0.0 ; /* reset reduced chi-squared */
00412
00413 /* loop through data points */
00414 for ( n = 0 ; n < (*ndat) ; n++ )
00415 {
00416 wn = wdat[n] ;
00417 if ( wn > 0.0 ) /* legal weight ? */
00418 {
00419 yd = ydat[n] - boltz( &xdat[(*xdim) * n], fpar ) ;
00420 boltz_deriv( &xdat[(*xdim) * n], fpar, epar ) ;
00421 chi2 += yd * yd * wn ; /* add to chi-squared */
00422 for ( j = 0 ; j < nfree ; j++ )
00423 {
00424 wd = epar[parptr[j]] * wn ; /* weighted derivative */
00425 vec[j] += yd * wd ; /* fill vector */
00426 for ( i = 0 ; i <= j ; i++ ) /* fill matrix */
00427 {
00428 matrix1[j][i] += epar[parptr[i]] * wd ;
00429 }
00430 }
00431 }
00432 }
00433 }
00434
00435
00436 /*----------------------------------------------------------------------------
00437 Function : get_vec()
00438 In : xdat: position
00439 xdim: factor of the indices of the position array
00440 ydat: real data
00441 wdat: weights
00442 ndat: number of data points
00443 fpar: function parameters
00444 epar: partial derivatives of the function
00445 npar: number of function parameters
00446 Out : integer (0 if it had worked, -5 or -7
00447 if diagonal element is wrong, or -6,
00448 if determinant is zero )
00449 Job : calculates the correction vector. The matrix has been
00450 built by get_mat(), we only have to rescale it for the
00451 current value of labda. The matrix is rescaled so that
00452 the diagonal gets the value 1 + labda.
00453 Next we calculate the inverse of the matrix and then
00454 the correction vector.
00455 ---------------------------------------------------------------------------*/
00456
00457 static int get_vec ( float * xdat,
00458 int * xdim,
00459 float * ydat,
00460 float * wdat,
00461 int * ndat,
00462 float * fpar,
00463 float * epar,
00464 int * npar )
00465 {
00466 double dj ;
00467 double dy ;
00468 double mii ;
00469 double mji ;
00470 double mjj ;
00471 double wn ;
00472 int i, j, n, r ;
00473
00474 /* loop to modify and scale the matrix */
00475 for ( j = 0 ; j < nfree ; j++ )
00476 {
00477 mjj = matrix1[j][j] ;
00478 if ( mjj <= 0.0 ) /* diagonal element wrong */
00479 {
00480 return -5 ;
00481 }
00482 mjj = sqrt( mjj ) ;
00483 for ( i = 0 ; i < j ; i++ )
00484 {
00485 mji = matrix1[j][i] / mjj / sqrt( matrix1[i][i] ) ;
00486 matrix2[i][j] = matrix2[j][i] = mji ;
00487 }
00488 matrix2[j][j] = 1.0 + labda ; /* scaled value on diagonal */
00489 }
00490
00491 if ( (r = inv_mat()) ) /* invert matrix inlace */
00492 {
00493 return r ;
00494 }
00495
00496 for ( i = 0 ; i < (*npar) ; i ++ )
00497 {
00498 epar[i] = fpar[i] ;
00499 }
00500
00501 /* loop to calculate correction vector */
00502 for ( j = 0 ; j < nfree ; j++ )
00503 {
00504 dj = 0.0 ;
00505 mjj = matrix1[j][j] ;
00506 if ( mjj <= 0.0) /* not allowed */
00507 {
00508 return -7 ;
00509 }
00510 mjj = sqrt ( mjj ) ;
00511 for ( i = 0 ; i < nfree ; i++ )
00512 {
00513 mii = matrix1[i][i] ;
00514 if ( mii <= 0.0 )
00515 {
00516 return -7 ;
00517 }
00518 mii = sqrt( mii ) ;
00519 dj += vec[i] * matrix2[j][i] / mjj / mii ;
00520 }
00521 epar[parptr[j]] += dj ; /* new parameters */
00522 }
00523 chi1 = 0.0 ; /* reset reduced chi-squared */
00524
00525 /* loop through the data points */
00526 for ( n = 0 ; n < (*ndat) ; n++ )
00527 {
00528 wn = wdat[n] ; /* get weight */
00529 if ( wn > 0.0 ) /* legal weight */
00530 {
00531 dy = ydat[n] - boltz( &xdat[(*xdim) * n], epar) ;
00532 chi1 += wdat[n] * dy * dy ;
00533 }
00534 }
00535 return 0 ;
00536 }
00537
00538
00539
00540 /*----------------------------------------------------------------------------
00541 Function : lsqfit()
00542 In : xdat: position, coordinates of data points.
00543 xdat is 2 dimensional: XDAT ( XDIM, NDAT )
00544 xdim: dimension of fit
00545 ydat: data points
00546 wdat: weights for data points
00547 ndat: number of data points
00548 fpar: on input contains initial estimates of the
00549 parameters for non-linear fits, on output the
00550 fitted parameters.
00551 epar: contains estimates of the errors in fitted parameters
00552 mpar: logical mask telling which parameters are free (non-zero)
00553 and which parameters are fixed (0)
00554 npar: number of function parameters ( free + fixed )
00555 tol: relative tolerance. lsqfit stops when successive iterations
00556 fail to produce a decrement in reduced chi-squared less
00557 than tol. If tol is less than the minimum tolerance
00558 possible, tol will be set to this value. This means
00559 that maximum accuracy can be obtained by setting
00560 tol = 0.0.
00561 its: maximum number of iterations
00562 lab: mixing parameter, lab determines the initial weight
00563 of steepest descent method relative to the Taylor method
00564 lab should be a small value (i.e. 0.01). lab can only
00565 be zero when the partial derivatives are independent
00566 of the parameters. In fact in this case lab should be
00567 exactly equal to zero.
00568 Out : returns number of iterations needed to achieve convergence
00569 according to tol. When this number is negative, the fitting
00570 was not continued because a fatal error occurred:
00571 -1 too many free parameters, maximum is 32
00572 -2 no free parameters
00573 -3 not enough degrees of freedom
00574 -4 maximum number of iterations too small to obtain
00575 a solution which satisfies tol.
00576 -5 diagonal of matrix contains elements which are zero
00577 -6 determinant of the coefficient matrix is zero
00578 -7 square root of a negative number
00579 Job : this is a routine for making a least-squares fit of a
00580 function to a set of data points. The method used is
00581 described in: Marquardt, J.Soc.Ind.Appl.Math. 11. 431 (1963).
00582 This method is a mixture of the steepest descent method
00583 and the Taylor method.
00584 ---------------------------------------------------------------------------*/
00585
00586 int lsqfit ( float * xdat,
00587 int * xdim,
00588 float * ydat,
00589 float * wdat,
00590 int * ndat,
00591 float * fpar,
00592 float * epar,
00593 int * mpar,
00594 int * npar,
00595 float * tol ,
00596 int * its ,
00597 float * lab )
00598 {
00599 int i, n, r ;
00600 int itc ; /* fate of fit */
00601 int found ; /* fit converged: 1, not yet converged: 0 */
00602 int nuse ; /* number of useable data points */
00603 double tolerance ; /* accuracy */
00604
00605 itc = 0 ; /* fate of fit */
00606 found = 0 ; /* reset */
00607 nfree = 0 ; /* number of free parameters */
00608 nuse = 0 ; /* number of legal data points */
00609
00610 if ( *tol < (FLT_EPSILON * 10.0 ) )
00611 {
00612 tolerance = FLT_EPSILON * 10.0 ; /* default tolerance */
00613 }
00614 else
00615 {
00616 tolerance = *tol ; /* tolerance */
00617 }
00618
00619 labda = fabs( *lab ) * LABFACA ; /* start value for mixing parameter */
00620 for ( i = 0 ; i < (*npar) ; i++ )
00621 {
00622 if ( mpar[i] )
00623 {
00624 if ( nfree > NPAR ) /* too many free parameters */
00625 {
00626 return -1 ;
00627 }
00628 parptr[nfree++] = i ; /* a free parameter */
00629 }
00630 }
00631
00632 if (nfree == 0) /* no free parameters */
00633 {
00634 return -2 ;
00635 }
00636
00637 for ( n = 0 ; n < (*ndat) ; n++ )
00638 {
00639 if ( wdat[n] > 0.0 ) /* legal weight */
00640 {
00641 nuse ++ ;
00642 }
00643 }
00644
00645 if ( nfree >= nuse )
00646 {
00647 return -3 ; /* no degrees of freedom */
00648 }
00649 if ( labda == 0.0 ) /* linear fit */
00650 {
00651 for ( i = 0 ; i < nfree ; fpar[parptr[i++]] = 0.0 ) ; /* initialize fpar array */
00652 get_mat ( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar*/ ) ;
00653 r = get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00654 if ( r ) /* error */
00655 {
00656 return r ;
00657 }
00658 for ( i = 0 ; i < (*npar) ; i++ )
00659 {
00660 fpar[i] = epar[i] ; /* save new parameters */
00661 epar[i] = 0.0 ; /* and set errors to zero */
00662 }
00663 chi1 = sqrt( chi1 / (double) (nuse - nfree) ) ;
00664 for ( i = 0 ; i < nfree ; i++ )
00665 {
00666 if ( (matrix1[i][i] <= 0.0 ) || (matrix2[i][i] <= 0.0) )
00667 {
00668 return -7 ;
00669 }
00670 epar[parptr[i]] = chi1 * sqrt( matrix2[i][i] ) / sqrt( matrix1[i][i] ) ;
00671 }
00672 }
00673 else /* non-linear fit */
00674 {
00675 /*----------------------------------------------------------------
00676 * the non-linear fit uses the steepest descent method in combination
00677 * with the Taylor method. The mixing of these methods is controlled
00678 * by labda. In the outer loop ( called the iteration loop ) we build
00679 * the matrix and calculate the correction vector. In the inner loop
00680 * (called the interpolation loop) we check whether we have obtained a
00681 * better solution than the previous one. If so, we leave the inner loop
00682 * else we increase labda ( give more weight to the steepest descent method)
00683 * calculate the correction vector and check again. After the inner loop
00684 * we do a final check on the goodness of the fit and if this satisfies
00685 * the tolerance we calculate the errors of the fitted parameters.
00686 */
00687 while ( !found ) /* iteration loop */
00688 {
00689 if ( itc++ == (*its) ) /* increase iteration counter */
00690 {
00691 return -4 ;
00692 }
00693 get_mat( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar */) ;
00694
00695 /*-------------------------------------------------------------
00696 * here we decrease labda since we may assume that each iteration
00697 * brings us closer to the answer.
00698 */
00699 if ( labda > LABMINA )
00700 {
00701 labda = labda / LABFACA ; /* decrease labda */
00702 }
00703 r = get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00704
00705 if ( r ) /* error */
00706 {
00707 return r ;
00708 }
00709
00710 while ( chi1 >= chi2 ) /* interpolation loop */
00711 {
00712 /*-----------------------------------------------------------
00713 * The next statement is based on experience, not on the mathematics
00714 * of the problem. It is assumed that we have reached convergence
00715 * when the pure steepest descent method does not produce a better
00716 * solution.
00717 */
00718 if ( labda > LABMAXA ) /* assume solution found */
00719 {
00720 break ;
00721 }
00722 labda = labda * LABFACA ; /* increase mixing parameter */
00723 r = get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00724
00725 if ( r ) /* error */
00726 {
00727 return r ;
00728 }
00729 }
00730
00731 if ( labda <= LABMAXA ) /* save old parameters */
00732 {
00733 for ( i = 0 ; i < *npar ; i++ )
00734 {
00735 fpar[i] = epar[i] ;
00736 }
00737 }
00738 if ( (fabs( chi2 - chi1 ) <= (tolerance * chi1)) || (labda > LABMAXA) )
00739 {
00740 /*---------------------------------------------------------------------
00741 * we have a satisfying solution, so now we need to calculate the correct
00742 * errors of the fitted parameters. This we do by using the pure Taylor
00743 * method because we are very close to the real solution.
00744 */
00745 labda = LABMINA ; /* for Taylor solution */
00746 get_mat ( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar */) ;
00747 r = get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00748
00749 if ( r ) /* error */
00750 {
00751 return r ;
00752 }
00753 for ( i = 0 ; i < (*npar) ; i++ )
00754 {
00755 epar[i] = 0.0 ; /* set error to zero */
00756 }
00757 chi2 = sqrt ( chi2 / (double) (nuse - nfree) ) ;
00758
00759 for ( i = 0 ; i < nfree ; i++ )
00760 {
00761 if ( (matrix1[i][i] <= 0.0) || (matrix2[i][i] <= 0.0) )
00762 {
00763 return -7 ;
00764 }
00765 epar[parptr[i]] = chi2 * sqrt( matrix2[i][i] ) / sqrt( matrix1[i][i] ) ;
00766 }
00767 found = 1 ; /* we found a solution */
00768 }
00769 }
00770 }
00771 return itc ; /* return number of iterations */
00772 }
00773
00774 /*----------------------------------------------------------------------------
00775 Function : fitSlitsBoltz()
00776 In : lineImage: emission line frame
00777 par: fit parameter data structure of fitted lines
00778 slit_pos: allocated dummy array for the slitlet positions [32][2]
00779 box_length: pixel length of the row box within the fit is done
00780 y_box: small box in spectral direction within the slitlet
00781 may lie.
00782 diff_tol: maximum tolerable difference of the resulting fit position
00783 with respect to the expected position. If difference is
00784 greater the expected position is taken.
00785 Out : slit_pos: beginning and end position of the slitlets to
00786 sub-pixel accuracy
00787 0 if it worked,
00788 -1 if there was no line image given,
00789 -2 if there were no line fit parameters given,
00790 -3 if there was no dummy array for the slit positions
00791 allocated
00792 -4 if the given box length is impossible
00793 -5 if the given y box length is impossible
00794 -6 if the given difference tolerance is too small
00795 -7 if there were no emission lines found in the first image columns
00796 -8 if not all slitlets could be found
00797 Job : fits the beginning and end position of the slitlets
00798 by using non-linear least square fitting of a Boltzmann function
00799 fits a Boltzmann function to the slitlet edges exposed and indicated
00800 by the brightest emission lines. To achieve this, the fit
00801 parameters are used to find the brightest emission line
00802 and to get its position for each column.
00803 The least squares fit is done by using a box smaller than
00804 the size of two slitlets
00805 ---------------------------------------------------------------------------*/
00806
00807 int fitSlitsBoltz ( OneImage * lineImage,
00808 FitParams ** par,
00809 float ** slit_pos,
00810 int box_length,
00811 float y_box,
00812 float diff_tol )
00813 {
00814 float position[lineImage -> lx] ;
00815 int * edge, * edgeclean ;
00816 int * dummyedge ;
00817 int * pos_row, * pos_rowclean ;
00818 Vector * box_buffer ;
00819 Vector * half_buffer ;
00820 Vector * in_buffer ;
00821 float max_intensity ;
00822 float row_pos ;
00823 int row, col ;
00824 int i, j, k, m, n, ed ;
00825 int found, init1 ;
00826 int line ;
00827 int nel, n_right, left_right ;
00828 int n_buf, edge_ind, shift ;
00829 int column ;
00830 int slit_length ;
00831 int agreed ;
00832 int bad_line ;
00833 int margin ;
00834 int iters, xdim, ndat ;
00835 int numpar, its ;
00836 int * mpar ;
00837 float * xdat, * wdat ;
00838 float tol, lab ;
00839 float fitpar[NPAR] ;
00840 float dervpar[NPAR] ;
00841 float minval, maxval ;
00842 float min ;
00843 float pos, last_pos ;
00844 int old_col=0;
00845 int old_pos=0;
00846
00847
00848 slit_length = SLITLENGTH ;
00849 if ( NullImage == lineImage )
00850 {
00851 cpl_msg_error( "fitSlitsBoltz:"," no line image given!\n" ) ;
00852 return -1 ;
00853 }
00854
00855 if ( NULL == par )
00856 {
00857 cpl_msg_error( "fitSlitsBoltz:"," no line fit parameters given!\n" ) ;
00858 return -2 ;
00859 }
00860
00861 if ( NULL == slit_pos )
00862 {
00863 cpl_msg_error( "fitSlitsBoltz:"," no position array allocated!\n" ) ;
00864 return -3 ;
00865 }
00866
00867 if ( box_length < 4 ||
00868 box_length > 2*slit_length )
00869 {
00870 cpl_msg_error( "fitSlitsBoltz:"," wrong fitting box length given!\n" ) ;
00871 return -4 ;
00872 }
00873
00874 if ( y_box <= 0. || y_box > 6. )
00875 {
00876 cpl_msg_error( "fitSlitsBoltz:"," wrong y box length given!\n" ) ;
00877 return -5 ;
00878 }
00879
00880 if ( diff_tol < 1. )
00881 {
00882 cpl_msg_error( "fitSlitsBoltz:"," diff_tol too small!\n" ) ;
00883 return -6 ;
00884 }
00885
00886 /* allocate memory for the edges and the row positon of the slitlets */
00887 edge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
00888 dummyedge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
00889 edgeclean = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ;
00890 pos_row = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
00891 pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ;
00892
00893
00894 /* --------------------------------------------------------------------------
00895 * go through the first image columns and the fit parameters and find the line
00896 * with the highest intensity
00897 */
00898 agreed = -1 ;
00899 bad_line = -1 ;
00900 while( agreed == -1 )
00901 {
00902 found = -1 ;
00903 max_intensity = -FLT_MAX ;
00904 for ( col = 0 ; col < slit_length ; col++ )
00905 {
00906 for ( i = 0 ; i < par[0]->n_params ; i++ )
00907 {
00908 if ( par[i]->column == col && par[i]->line != bad_line )
00909 {
00910 if ( par[i]->fit_par[0] > max_intensity )
00911 {
00912 if ( par[i]->fit_par[1] >= 1. && par[i]->fit_par[2] > 0. )
00913 {
00914 max_intensity = par[i]->fit_par[0] ;
00915 found = i ;
00916 }
00917 }
00918 }
00919 }
00920 }
00921
00922 /* --------------------------------------------------------------------
00923 * check if the found line is usable and if the neighbouring line
00924 * have intensity on near rows in neighbouring slitlets
00925 */
00926 line = par[found]->line ;
00927 column = par[found]->column ;
00928 row_pos = par[found]->fit_par[2] ;
00929 if ( found >= 0 && max_intensity > 0. )
00930 {
00931 for ( i = 0 ; i < par[0]->n_params ; i++ )
00932 {
00933 if ( par[i]->line == line-1 &&
00934 par[i]->column == column + slit_length )
00935 {
00936 if ( par[i]->fit_par[2] <= (row_pos + y_box) &&
00937 par[i]->fit_par[2] >= (row_pos - y_box) )
00938 {
00939 bad_line = line ;
00940 }
00941 }
00942 }
00943 if ( bad_line != line )
00944 {
00945 agreed = 1 ;
00946 break ;
00947 }
00948 }
00949 else
00950 {
00951 cpl_msg_error("fitSlitsBoltz:"," no emission line found in the first image columns\n") ;
00952 cpl_free( edge ) ;
00953 cpl_free( pos_row ) ;
00954 cpl_free( edgeclean ) ;
00955 cpl_free( dummyedge ) ;
00956 cpl_free( pos_rowclean ) ;
00957 return -7 ;
00958 }
00959 }
00960
00961
00962
00963 if ( agreed == -1 )
00964 {
00965 cpl_msg_error("fitSlitsBoltz:"," no emission line found in the first image columns\n") ;
00966 cpl_free( edge ) ;
00967 cpl_free( pos_row ) ;
00968 cpl_free( edgeclean ) ;
00969 cpl_free( dummyedge ) ;
00970 cpl_free( pos_rowclean ) ;
00971 return -7 ;
00972 }
00973
00974
00975 /* now find and store the raw edge positions of the found slitlet */
00976 n = 0 ;
00977 ed = 0 ;
00978 /* was for ( col = 0 ; col < lineImage -> lx - slit_length/2 ; col++ ) */
00979 for ( col = slit_length/2 ; col < lineImage -> lx - slit_length/2 ; col++ )
00980 {
00981 for ( i = 0 ; i < par[0]->n_params ; i++ )
00982 {
00983 /*
00984 printf("p1=%f p2=%f p3=%f\n",
00985 par[i]->fit_par[0],par[i]->fit_par[1],par[i]->fit_par[2]);
00986 */
00987 if ( par[i]->column == col && par[i]->line == line )
00988 {
00989 if ( par[i]->fit_par[0] > 0. &&
00990 par[i]->fit_par[1] >= 1. &&
00991 par[i]->fit_par[2] > 0. )
00992 {
00993 position[n] = par[i]->fit_par[2] ;
00994 old_pos=position[n];
00995 if ( n > 0 &&
00996 fabs(position[n] - position[n-1]) > y_box &&
00997 (col-old_col) > (slit_length-SLIT_POS_ERR) )
00998 {
00999
01000 old_col=col;
01001 edge[ed] = col ;
01002 pos_row[ed] = nint( position[n-1] ) ;
01003 /* printf("edge[%d]=%d , pos_row=%d\n",ed,edge[ed],pos_row[ed]); */
01004 ed++ ;
01005 if ( col >= lineImage -> lx - slit_length - SLIT_POS_ERR ) {
01006 pos_row[ed] = nint( position[n] ) ;
01007 }
01008 } else if ( ((col-old_col) > (slit_length+SLIT_POS_ERR)) && (col>120) ) {
01009 old_col=col;
01010 edge[ed] = col ;
01011 pos_row[ed] = nint( position[n-1] ) ;
01012 cpl_msg_warning("fitSlitsBoltz:","added1 slitlet edge[%d]=%d, pos_row=%d",
01013 ed,edge[ed],pos_row[ed]);
01014 ed++ ;
01015 if ( col >= lineImage -> lx - slit_length - SLIT_POS_ERR ) {
01016 pos_row[ed] = nint( position[n] ) ;
01017 }
01018 }
01019 n++ ;
01020 }
01021 } else if ( ((col-old_col) > (slit_length+SLIT_POS_ERR)) &&
01022 (col>120) ) {
01023 /*
01024 printf("check col=%d col-old_col=%d check=%d\n",
01025 col,(col-old_col),(slit_length+SLIT_POS_ERR));
01026 */
01027 position[n] = old_pos ;
01028
01029 old_col+=slit_length;
01030 edge[ed] = old_col; ;
01031 pos_row[ed] = nint( position[n-1] ) ;
01032
01033
01034 cpl_msg_warning("fitSlitsBoltz:","added2 slitlet edge[%d]=%d, pos_row=%d",
01035 ed,edge[ed],pos_row[ed]);
01036 ed++ ;
01037 if ( old_col >= lineImage -> lx - slit_length - SLIT_POS_ERR ) {
01038 pos_row[ed] = old_pos ;
01039 }
01040 n++;
01041 }
01042 }
01043 }
01044
01045
01046 if ( ed < (N_SLITLETS - 1) )
01047 {
01048 cpl_msg_error("fitSlitsBoltz:"," not enough slitlets, found: %d\n", ed) ;
01049 cpl_free( edge ) ;
01050 cpl_free( pos_row ) ;
01051 cpl_free( edgeclean ) ;
01052 cpl_free( dummyedge ) ;
01053 cpl_free( pos_rowclean ) ;
01054 return -8 ;
01055 }
01056
01057 /* now find the clean edge and row positions of the slitlets */
01058 /* printf("ed=%d\n",ed); */
01059 for ( i = 1 ; i <= ed ; i ++ )
01060 {
01061 if ( i == ed )
01062 {
01063 if ( (edge[i-1] - edge[i-2]) < slit_length - SLIT_LEN_ERR ||
01064 (edge[i-1] - edge[i-2]) > slit_length + SLIT_LEN_ERR )
01065 {
01066 /* printf("e(i-1)=%d e(i-2)=%d i=%d\n",edge[i-1], edge[i-2],i); */
01067 dummyedge[i-1] = -1 ;
01068 }
01069 }
01070 if (dummyedge[i-1] != -1)
01071 {
01072 dummyedge[i-1] = edge[i-1] ;
01073 }
01074 else
01075 {
01076 continue ;
01077 }
01078 if ( i < ed )
01079 {
01080 if ( (edge[i] - edge[i-1]) < slit_length - SLIT_LEN_ERR ||
01081 (edge[i] - edge[i-1]) > slit_length + SLIT_LEN_ERR )
01082 {
01083 /* printf("e(i)=%d e(i-1)=%d i=%d\n",edge[i], edge[i-1],i); */
01084 dummyedge[i] = -1 ;
01085 }
01086 }
01087 if ( i+1 < ed && dummyedge[i] != -1 )
01088 {
01089 if ( (edge[i+1] - edge[i]) < slit_length - SLIT_LEN_ERR ||
01090 (edge[i+1] - edge[i]) > slit_length + SLIT_LEN_ERR )
01091 {
01092 /* printf("e(i+1)=%d e(i)=%d i=%d\n",edge[i+1], edge[i],i); */
01093 dummyedge[i+1] = -1 ;
01094 }
01095 }
01096 }
01097
01098 k = 0 ;
01099 for ( i = 0 ; i < ed ; i++ )
01100 {
01101 if ( dummyedge[i] != -1 && dummyedge[i] != 0 )
01102 {
01103 edgeclean[k] = dummyedge[i] ;
01104 pos_rowclean[k] = pos_row[i] ;
01105 k++ ;
01106 if( edgeclean[k-1] > (lineImage -> lx - slit_length -2*SLIT_LEN_ERR ) )
01107 {
01108 pos_rowclean[k] = pos_row[ed] ;
01109 }
01110 }
01111 }
01112 /*
01113 for ( i = 0 ; i < k ; i++ )
01114 {
01115 cpl_msg_warning("fitSlitsBoltz:","%d %d", edgeclean[i], pos_rowclean[i]);
01116 }
01117 */
01118 if ( k != N_SLITLETS - 1 )
01119 {
01120 cpl_msg_error("fitSlitsBoltz:"," wrong number of clean slitlets found: %d\n", k+1) ;
01121 cpl_free( edge ) ;
01122 cpl_free( pos_row ) ;
01123 cpl_free( edgeclean ) ;
01124 cpl_free( dummyedge ) ;
01125 cpl_free( pos_rowclean ) ;
01126 return -7 ;
01127 }
01128
01129 /* determine the margins of the fitting box outside the slitlets */
01130 margin = box_length / 2 ;
01131
01132 /* --------------------------------------------------------------------------
01133 * now go through the slitlets, search along each column within a box with
01134 * half width y_box the maximum value and store these found values in a buffer
01135 */
01136 if(
01137 ( (pos_rowclean[0]-nint(y_box)) < 0 ) ||
01138 ( (pos_rowclean[0]+nint(y_box)) >lineImage->ly )
01139 ) {
01140
01141 cpl_msg_error("fitSlitsBoltz:",
01142 "pos_rowclean[0] <0 something wrong!\n") ;
01143 cpl_free( edge ) ;
01144 cpl_free( pos_row ) ;
01145 cpl_free( edgeclean ) ;
01146 cpl_free( dummyedge ) ;
01147 cpl_free( pos_rowclean ) ;
01148 return -7 ;
01149
01150 }
01151
01152 for ( j = 0 ; j <= k ; j++ )
01153 {
01154 m = 0 ;
01155 if ( j == 0 )
01156 {
01157 box_buffer = newVector( edgeclean[0] + margin ) ;
01158 for( col = 0 ; col < edgeclean[0] + margin ; col++ )
01159 {
01160 maxval = -FLT_MAX ;
01161 for ( row = pos_rowclean[0] - nint(y_box) ; row <= pos_rowclean[0] + nint(y_box) ; row++ )
01162 {
01163 if ( maxval < lineImage->data[col + lineImage->lx*row] )
01164 {
01165 maxval = lineImage -> data[col + lineImage->lx*row] ;
01166 }
01167 }
01168 box_buffer->data[m] = maxval ;
01169 m++ ;
01170 }
01171 }
01172 else if ( j < k )
01173 {
01174 box_buffer = newVector( edgeclean[j] - edgeclean[j-1] + 2*margin ) ;
01175 for ( col = edgeclean[j - 1] - margin ; col < edgeclean[j] + margin ; col++ )
01176 {
01177 maxval = -FLT_MAX ;
01178 for ( row = pos_rowclean[j] - nint(y_box) ; row <= pos_rowclean[j] + nint(y_box) ; row++ )
01179 {
01180 if ( maxval < lineImage->data[col + lineImage->lx*row] )
01181 {
01182 maxval = lineImage -> data[col + lineImage->lx*row] ;
01183 }
01184 }
01185 box_buffer->data[m] = maxval ;
01186 m++ ;
01187 }
01188 }
01189 else
01190 {
01191 box_buffer = newVector( lineImage -> lx - edgeclean[k-1] + margin ) ;
01192 for ( col = edgeclean[k - 1] - margin ; col < lineImage -> lx ; col++ )
01193 {
01194 maxval = -FLT_MAX ;
01195 for ( row = pos_rowclean[k-2] - nint(y_box) ; row <= pos_rowclean[k-2] + nint(y_box) ; row++ )
01196 {
01197 if ( maxval < lineImage->data[col + lineImage->lx*row] )
01198 {
01199 maxval = lineImage -> data[col + lineImage->lx*row] ;
01200 }
01201 }
01202 if(maxval>0) box_buffer->data[m] = maxval ;
01203 else box_buffer->data[m] = 0;
01204 m++ ;
01205 }
01206 }
01207
01208 /* determine the minimum value in the box to get background1 value for the edge slitlets */
01209 min = FLT_MAX ;
01210 for ( i = 0 ; i < box_buffer->n_elements ; i++ )
01211 {
01212 if ( box_buffer -> data[i] < min )
01213 {
01214 min = box_buffer -> data[i] ;
01215 }
01216 }
01217
01218 for ( left_right = 0 ; left_right <= 1 ; left_right++ )
01219 {
01220 nel = 0 ;
01221 if ( left_right == 0 )
01222 {
01223 nel = box_buffer -> n_elements / 2 ;
01224 }
01225 else
01226 {
01227 if ( box_buffer -> n_elements % 2 == 0 )
01228 {
01229 nel = box_buffer -> n_elements / 2 ;
01230 }
01231 else
01232 {
01233 nel = box_buffer -> n_elements / 2 + 1 ;
01234 }
01235 }
01236
01237 /* now split the buffer in the midth in a left and right part for fitting */
01238 half_buffer = newVector( nel ) ;
01239 if ( left_right == 0 )
01240 {
01241 for ( i = 0 ; i < nel ; i++ )
01242 {
01243 half_buffer -> data[i] = box_buffer -> data[i] ;
01244 }
01245 }
01246 else
01247 {
01248 n_right = 0 ;
01249 for ( i = box_buffer -> n_elements - 1 ; i >= box_buffer -> n_elements - nel ; i-- )
01250 {
01251 half_buffer -> data[n_right] = box_buffer -> data[i] ;
01252 n_right++ ;
01253 }
01254 }
01255
01256 xdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
01257 wdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
01258 mpar = (int *) cpl_calloc( NPAR, sizeof (int) ) ;
01259
01260 /* set initial values for the fitting routine */
01261 minval = FLT_MAX ;
01262 maxval = -FLT_MAX ;
01263 for ( i = 0 ; i < nel ; i++ )
01264 {
01265 xdat[i] = i ;
01266 wdat[i] = 1.0 ;
01267 if ( half_buffer -> data[i] < minval )
01268 {
01269 minval = half_buffer -> data[i] ;
01270 }
01271 if ( half_buffer -> data[i] > maxval )
01272 {
01273 maxval = half_buffer -> data[i] ;
01274 }
01275 }
01276
01277 fitpar[0] = minval ;
01278 fitpar[1] = maxval ;
01279
01280 /* search for both positions of the half intensity of the hat within the buffer */
01281 init1 = -1 ;
01282 for ( i = 0 ; i < nel ; i++ )
01283 {
01284 if ( half_buffer -> data[i] >= ( maxval + minval ) / 2. )
01285 {
01286 init1 = i ;
01287 break ;
01288 }
01289 }
01290
01291 /*------------------------------------------------------------------
01292 * if we have too few left background values (at the image edges)
01293 * the left margin of the buffer to fit is filled with the minimal
01294 * values in order to get a good fit
01295 */
01296
01297 edge_ind = 0 ;
01298 if ( init1 < 3 )
01299 {
01300 n_buf = half_buffer->n_elements + margin ;
01301 in_buffer = newVector( n_buf ) ;
01302 for ( i = 0 ; i < margin ; i++ )
01303 {
01304 in_buffer -> data[i] = min ;
01305 }
01306 shift = 0 ;
01307 for ( i = margin ; i < n_buf ; i++ )
01308 {
01309 in_buffer -> data[i] = half_buffer -> data[shift] ;
01310 shift++ ;
01311 }
01312 destroyVector ( half_buffer ) ;
01313 half_buffer = newVector ( n_buf ) ;
01314 for ( i = 0 ; i < n_buf ; i++ )
01315 {
01316 half_buffer -> data[i] = in_buffer -> data[i] ;
01317 }
01318 edge_ind = 1 ;
01319 init1 += margin ;
01320 destroyVector ( in_buffer ) ;
01321 }
01322
01323 /* determine the initial positions from the found values */
01324 if ( init1 != -1 )
01325 {
01326 fitpar[2] = (float)init1 ;
01327 }
01328 fitpar[3] = 1. ;
01329
01330 for ( i = 0 ; i < NPAR ; i++ )
01331 {
01332 mpar[i] = 1 ;
01333 dervpar[i] = 0. ;
01334 }
01335
01336 xdim = XDIMA ;
01337 ndat = nel ;
01338 numpar = NPAR ;
01339 tol = TOLA ;
01340 lab = LABA ;
01341 its = ITSA ;
01342
01343 /* finally, do the least squares fit over the buffer data */
01344 if ( 0 > ( iters = lsqfit( xdat, &xdim, half_buffer -> data, wdat, &ndat, fitpar,
01345 dervpar, mpar, &numpar, &tol, &its, &lab )) )
01346 {
01347 /* if the fit doesn't succeed the initial values are taken */
01348 cpl_msg_warning ("lsqfit:"," least squares fit failed, error no.: %d in slitlet: %d\n", iters, j) ;
01349 fitpar[2] = (float)init1 ;
01350 }
01351
01352 pos = fitpar[2] ;
01353 if ( edge_ind == 1 )
01354 {
01355 pos -= (float)margin ;
01356 }
01357
01358 /*------------------------------------------------------------------------------------------
01359 * now discern the left and the right edge fit of the slitlets and associate the fit results
01360 * with the absolute positions in the image
01361 * consider the difference of the fitted slit position to the expected position
01362 * and decide wether the fit is taken or the expected value is taken.
01363 */
01364 if ( left_right == 0 )
01365 {
01366 /* take care of the column position of the fit boxes to get the absolute positions */
01367 if ( j == 0 )
01368 {
01369 if ( fabs(pos - ((float)edgeclean[0] - 1. - (float)slit_length)) < diff_tol )
01370 {
01371 slit_pos[0][0] = pos ;
01372 }
01373 else
01374 {
01375 cpl_msg_warning("lsqfit:","something wrong with fitted left position of slitlet 0\n") ;
01376 if ( (float) edgeclean[0] - 1. - (float)slit_length < 0. )
01377 {
01378 slit_pos[0][0] = 0. ;
01379 }
01380 else
01381 {
01382 slit_pos[0][0] = (float)edgeclean[0] - 1. - (float)slit_length ;
01383 }
01384 }
01385 }
01386 else if ( j < k )
01387 {
01388 if ( fabs( pos - (float)margin ) < diff_tol )
01389 {
01390 slit_pos[j][0] = pos + (float)edgeclean[j-1] - (float)margin ;
01391 }
01392 else
01393 {
01394 cpl_msg_warning("lsqfit:","something wrong with fitted left position of slitlet %d\n", j) ;
01395 slit_pos[j][0] = (float)edgeclean[j-1] - 1. ;
01396 }
01397 }
01398 else
01399 {
01400 if ( fabs( pos - (float)margin ) < diff_tol )
01401 {
01402 slit_pos[k][0] = pos + (float)edgeclean[k-1] - (float)margin ;
01403 }
01404 else
01405 {
01406 cpl_msg_warning("lsqfit:","something wrong with fitted left position of slitlet %d\n", j) ;
01407 slit_pos[k][0] = (float)edgeclean[k-1] - 1. ;
01408 }
01409 }
01410 }
01411 else
01412 {
01413 /* take care of the column position of the fit boxes to get the absolute positions */
01414 if ( j == 0 )
01415 {
01416 if ( fabs( (float)box_buffer->n_elements - pos - (float)edgeclean[0] ) < diff_tol )
01417 {
01418 slit_pos[0][1] = (float)(box_buffer->n_elements - 1) - pos ;
01419 }
01420 else
01421 {
01422 cpl_msg_warning("lsqfit:","something wrong with fitted right position of slitlet 0\n") ;
01423 slit_pos[0][1] = (float)edgeclean[0] - 1. ;
01424 }
01425 }
01426 else if ( j < k )
01427 {
01428 if ( fabs( (float)box_buffer->n_elements - pos
01429 + (float)edgeclean[j-1] - (float)margin - (float)edgeclean[j] ) < diff_tol )
01430 {
01431 slit_pos[j][1] = (float)(box_buffer->n_elements - 1) - pos
01432 + (float)edgeclean[j-1] - (float)margin ;
01433 }
01434 else
01435 {
01436 cpl_msg_warning("lsqfit:","something wrong with fitted right position of slitlet %d\n", j) ;
01437 slit_pos[j][1] = (float)edgeclean[j] - 1. ;
01438 }
01439 }
01440 else
01441 {
01442 if ( edgeclean[k-1] + slit_length > lineImage -> lx )
01443 {
01444 last_pos = (float)(lineImage -> lx - 1) ;
01445 }
01446 else
01447 {
01448 last_pos = (float)(edgeclean[k-1] - 1 + slit_length) ;
01449 }
01450 if ( fabs( (float)(box_buffer->n_elements - 1) - pos
01451 + (float)edgeclean[k-1] - (float)margin - last_pos ) < diff_tol )
01452 {
01453 slit_pos[k][1] = (float)(box_buffer->n_elements - 1) - pos
01454 + (float)edgeclean[k-1] - (float)margin ;
01455 }
01456 else
01457 {
01458 cpl_msg_warning("lsqfit:","something wrong with fitted right position of slitlet %d\n", j) ;
01459 slit_pos[k][1] = last_pos ;
01460 }
01461 }
01462 }
01463
01464 destroyVector ( half_buffer ) ;
01465 cpl_free( xdat ) ;
01466 cpl_free( wdat ) ;
01467 cpl_free( mpar ) ;
01468 }
01469 destroyVector ( box_buffer ) ;
01470 }
01471
01472
01473 cpl_free( edge ) ;
01474 cpl_free( pos_row ) ;
01475 cpl_free( edgeclean ) ;
01476 cpl_free( dummyedge ) ;
01477 cpl_free( pos_rowclean ) ;
01478 return 0 ;
01479 }
01480
01481 /*----------------------------------------------------------------------------
01482 Function : fitSlitsBoltzSingleLine()
01483 In : lineImage: emission line frame
01484 slit_pos: allocated dummy array for the slitlet positions [min32][2]
01485 box_length: pixel length of the row box within the fit is done
01486 y_box: small box in spectral direction within the slitlet
01487 may lie.
01488 low_pos, high_pos: pixel positions in spectral direction between which the line
01489 should be located.
01490 Out : slit_pos: beginning and end position of the slitlets to
01491 sub-pixel accuracy
01492 0 if it worked,
01493 -1 if it failed,
01494 Job : fits the beginning and end position of the slitlets
01495 by using non-linear least square fitting of a Boltzmann function
01496 fits a Boltzmann function to the slitlet edges exposed and indicated
01497 by the brightest emission lines. The slitlet is searched within
01498 user given positions.
01499 The least squares fit is done by using a box smaller than
01500 the size of two slitlets
01501 ---------------------------------------------------------------------------*/
01502
01503 int fitSlitsBoltzSingleLine ( OneImage * lineImage,
01504 float ** slit_pos,
01505 int box_length,
01506 float y_box,
01507 int low_pos,
01508 int high_pos )
01509 {
01510 int position[lineImage->lx] ;
01511 int * edge, * edgeclean ;
01512 int * dummyedge ;
01513 int * pos_row, * pos_rowclean ;
01514 Vector * box_buffer ;
01515 Vector * half_buffer ;
01516 Vector * in_buffer ;
01517 int found_row ;
01518 int row, col ;
01519 int i, j, k, m, ed ;
01520 int init1 ;
01521 int nel, n_right, left_right ;
01522 int n_buf, edge_ind, shift ;
01523 int slit_length ;
01524 int margin ;
01525 int iters, xdim, ndat ;
01526 int numpar, its ;
01527 int * mpar ;
01528 float * xdat, * wdat ;
01529 float tol, lab ;
01530 float fitpar[NPAR] ;
01531 float dervpar[NPAR] ;
01532 float minval, maxval ;
01533 float min ;
01534 float pos, last_pos ;
01535
01536 slit_length = SLITLENGTH ;
01537
01538 if ( NullImage == lineImage )
01539 {
01540 cpl_msg_error( "fitSlitsBoltzSingleLine:"," no line image given!\n" ) ;
01541 return -1 ;
01542 }
01543
01544 if ( NULL == slit_pos )
01545 {
01546 cpl_msg_error( "fitSlitsBoltzSingleLine:"," no position array allocated!\n" ) ;
01547 return -1 ;
01548 }
01549
01550 if ( box_length < 4 ||
01551 box_length > 2*slit_length )
01552 {
01553 cpl_msg_error( "fitSlitsBoltzSingleLine:"," wrong fitting box length given!\n" ) ;
01554 return -1 ;
01555 }
01556
01557 if ( y_box <= 0. || y_box > 6. )
01558 {
01559 cpl_msg_error( "fitSlitsBoltzSingleLine:"," wrong y box length given!\n" ) ;
01560 return -1 ;
01561 }
01562
01563 if ( low_pos >= high_pos || low_pos < 0 || high_pos <= 0 || high_pos >= lineImage->lx )
01564 {
01565 cpl_msg_error( "fitSlitsBoltzSingleLine:"," wrong user given search positions!\n" ) ;
01566 return -1 ;
01567 }
01568
01569 /* allocate memory for the edges and the row position of the slitlets */
01570 edge = (int*) cpl_calloc( lineImage->lx/2, sizeof(int) ) ;
01571 dummyedge = (int*) cpl_calloc( lineImage->lx/2, sizeof(int) ) ;
01572 edgeclean = (int*) cpl_calloc( lineImage->lx/2, sizeof(int) ) ;
01573 pos_row = (int*) cpl_calloc( lineImage->lx/2, sizeof(int) ) ;
01574 pos_rowclean = (int*) cpl_calloc( lineImage->lx/2, sizeof(int) ) ;
01575
01576 /* now search for the maximum between the given positions for each col */
01577 for ( col = 0 ; col < lineImage -> lx ; col++ )
01578 {
01579 found_row = -1 ;
01580 maxval = -FLT_MAX ;
01581 for ( row = low_pos ; row <= high_pos ; row++ )
01582 {
01583 if ( maxval < lineImage->data[col+row*lineImage->lx] )
01584 {
01585 maxval = lineImage->data[col+row*lineImage->lx] ;
01586 found_row = row ;
01587 }
01588 }
01589 if ( maxval > -FLT_MAX && found_row > low_pos )
01590 {
01591 position[col] = found_row ;
01592 }
01593 else
01594 {
01595 position[col] = 0 ;
01596 }
01597 }
01598
01599 /* now find and store the raw edge positions of the found slitlet */
01600 ed = 0 ;
01601 for ( col = 0 ; col < (lineImage->lx) - 1 ; col++ )
01602 {
01603 if ( position[col] > 0 && position[col+1] > 0 &&
01604 abs(position[col+1] - position[col]) > 10 )
01605 {
01606 edge[ed] = col ;
01607 pos_row[ed] = position[col] ;
01608 ed++ ;
01609 }
01610
01611 }
01612 if (ed <= 1)
01613 {
01614 cpl_msg_error( "fitSlitsBoltzSingleLine:"," no slitlets found!\n" ) ;
01615 cpl_free( edge ) ;
01616 cpl_free( pos_row ) ;
01617 cpl_free( edgeclean ) ;
01618 cpl_free( dummyedge ) ;
01619 cpl_free( pos_rowclean ) ;
01620 return -1 ;
01621 }
01622
01623 /* now find the clean edge and row positions of the slitlets */
01624 for ( i = 1 ; i <= ed ; i ++ )
01625 {
01626 if ( i == ed )
01627 {
01628 if ( (edge[i-1] - edge[i-2]) < slit_length - SLIT_LEN_ERR ||
01629 (edge[i-1] - edge[i-2]) > slit_length + SLIT_LEN_ERR )
01630 {
01631 dummyedge[i-1] = -1 ;
01632 }
01633 }
01634 if (dummyedge[i-1] != -1)
01635 {
01636 dummyedge[i-1] = edge[i-1] ;
01637 }
01638 else
01639 {
01640 continue ;
01641 }
01642 if ( i < ed )
01643 {
01644 if ( (edge[i] - edge[i-1]) < slit_length - SLIT_LEN_ERR ||
01645 (edge[i] - edge[i-1]) > slit_length + SLIT_LEN_ERR )
01646 {
01647 dummyedge[i] = -1 ;
01648 }
01649 }
01650 if ( i+1 < ed && dummyedge[i] != -1 )
01651 {
01652 if ( (edge[i+1] - edge[i]) < slit_length - SLIT_LEN_ERR ||
01653 (edge[i+1] - edge[i]) > slit_length + SLIT_LEN_ERR )
01654 {
01655 dummyedge[i+1] = -1 ;
01656 }
01657 }
01658 }
01659
01660 k = 0 ;
01661 for ( i = 0 ; i < ed ; i++ )
01662 {
01663 if ( dummyedge[i] != -1 && dummyedge[i] != 0 )
01664 {
01665 edgeclean[k] = dummyedge[i] ;
01666 pos_rowclean[k] = pos_row[i] ;
01667 k++ ;
01668 if( edgeclean[k-1] > (lineImage -> lx - slit_length - 2*SLIT_LEN_ERR ) )
01669 {
01670 pos_rowclean[k] = pos_row[ed] ;
01671 }
01672 }
01673 }
01674
01675 /* determine the margins of the fitting box outside the slitlets */
01676 margin = box_length / 2 ;
01677
01678 /* --------------------------------------------------------------------------
01679 * now go through the slitlets, search along each column within a box with
01680 * half width y_box the maximum value and store these found values in a buffer
01681 */
01682 for ( j = 0 ; j <= k ; j++ )
01683 {
01684 m = 0 ;
01685 if ( j == 0 )
01686 {
01687 box_buffer = newVector( edgeclean[0] + margin ) ;
01688 for( col = 0 ; col < edgeclean[0] + margin ; col++ )
01689 {
01690 maxval = -FLT_MAX ;
01691 for ( row = pos_rowclean[0] - nint(y_box) ; row <= pos_rowclean[0] + nint(y_box) ; row++ )
01692 {
01693 if ( maxval < lineImage->data[col + lineImage->lx*row] )
01694 {
01695 maxval = lineImage -> data[col + lineImage->lx*row] ;
01696 }
01697 }
01698 box_buffer->data[m] = maxval ;
01699 m++ ;
01700 }
01701 }
01702 else if ( j < k )
01703 {
01704 box_buffer = newVector( edgeclean[j] - edgeclean[j-1] + 2*margin ) ;
01705 for ( col = edgeclean[j - 1] - margin ; col < edgeclean[j] + margin ; col++ )
01706 {
01707 maxval = -FLT_MAX ;
01708 for ( row = pos_rowclean[j] - nint(y_box) ; row <= pos_rowclean[j] + nint(y_box) ; row++ )
01709 {
01710 if ( maxval < lineImage->data[col + lineImage->lx*row] )
01711 {
01712 maxval = lineImage -> data[col + lineImage->lx*row] ;
01713 }
01714 }
01715 box_buffer->data[m] = maxval ;
01716 m++ ;
01717 }
01718 }
01719 else
01720 {
01721 box_buffer = newVector( lineImage -> lx - edgeclean[k-1] + margin ) ;
01722 for ( col = edgeclean[k - 1] - margin ; col < lineImage -> lx ; col++ )
01723 {
01724 if ( col < 0 )
01725 {
01726 col = 0 ;
01727 }
01728
01729 maxval = -FLT_MAX ;
01730 for ( row = pos_rowclean[k] - nint(y_box) ; row <= pos_rowclean[k] + nint(y_box) ; row++ )
01731 {
01732 if ( row < 0 )
01733 {
01734 continue ;
01735 }
01736 if ( maxval < lineImage->data[col + row * lineImage->lx] )
01737 {
01738 maxval = lineImage -> data[col + row * lineImage->lx] ;
01739 }
01740 }
01741 box_buffer->data[m] = maxval ;
01742 m++ ;
01743 }
01744 }
01745
01746 /* determine the minimum value in the box to get background1 value for the edge slitlets */
01747 min = FLT_MAX ;
01748 for ( i = 0 ; i < box_buffer->n_elements ; i++ )
01749 {
01750 if ( box_buffer -> data[i] < min )
01751 {
01752 min = box_buffer -> data[i] ;
01753 }
01754 }
01755
01756 for ( left_right = 0 ; left_right <= 1 ; left_right++ )
01757 {
01758 nel = 0 ;
01759 if ( left_right == 0 )
01760 {
01761 nel = box_buffer -> n_elements / 2 ;
01762 }
01763 else
01764 {
01765 if ( box_buffer -> n_elements % 2 == 0 )
01766 {
01767 nel = box_buffer -> n_elements / 2 ;
01768 }
01769 else
01770 {
01771 nel = box_buffer -> n_elements / 2 + 1 ;
01772 }
01773 }
01774
01775 /* now split the buffer in the midth in a left and right part for fitting */
01776 half_buffer = newVector( nel ) ;
01777 if ( left_right == 0 )
01778 {
01779 for ( i = 0 ; i < nel ; i++ )
01780 {
01781 half_buffer -> data[i] = box_buffer -> data[i] ;
01782 }
01783 }
01784 else
01785 {
01786 n_right = 0 ;
01787 for ( i = box_buffer -> n_elements - 1 ; i >= box_buffer -> n_elements - nel ; i-- )
01788 {
01789 half_buffer -> data[n_right] = box_buffer -> data[i] ;
01790 n_right++ ;
01791 }
01792 }
01793
01794 xdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
01795 wdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
01796 mpar = (int *) cpl_calloc( NPAR, sizeof (int) ) ;
01797
01798 /* set initial values for the fitting routine */
01799 minval = FLT_MAX ;
01800 maxval = -FLT_MAX ;
01801 for ( i = 0 ; i < nel ; i++ )
01802 {
01803 xdat[i] = i ;
01804 wdat[i] = 1.0 ;
01805 if ( half_buffer -> data[i] < minval )
01806 {
01807 minval = half_buffer -> data[i] ;
01808 }
01809 if ( half_buffer -> data[i] > maxval )
01810 {
01811 maxval = half_buffer -> data[i] ;
01812 }
01813 }
01814 fitpar[0] = minval ;
01815 fitpar[1] = maxval ;
01816
01817 /* search for both positions of the half intensity of the hat within the buffer */
01818 init1 = -1 ;
01819 for ( i = 0 ; i < nel ; i++ )
01820 {
01821 if ( half_buffer -> data[i] >= ( maxval + minval ) / 2. )
01822 {
01823 init1 = i ;
01824 break ;
01825 }
01826 }
01827
01828 /*------------------------------------------------------------------
01829 * if we have too few left background values (at the image edges)
01830 * the left margin of the buffer to fit is filled with the minimal
01831 * values in order to get a good fit
01832 */
01833
01834 edge_ind = 0 ;
01835 if ( init1 < 3 )
01836 {
01837 n_buf = half_buffer->n_elements + margin ;
01838 in_buffer = newVector( n_buf ) ;
01839 for ( i = 0 ; i < margin ; i++ )
01840 {
01841 in_buffer -> data[i] = min ;
01842 }
01843 shift = 0 ;
01844 for ( i = margin ; i < n_buf ; i++ )
01845 {
01846 in_buffer -> data[i] = half_buffer -> data[shift] ;
01847 shift++ ;
01848 }
01849 destroyVector ( half_buffer ) ;
01850 half_buffer = newVector ( n_buf ) ;
01851 for ( i = 0 ; i < n_buf ; i++ )
01852 {
01853 half_buffer -> data[i] = in_buffer -> data[i] ;
01854 }
01855 edge_ind = 1 ;
01856 init1 += margin ;
01857 destroyVector ( in_buffer ) ;
01858 }
01859
01860 /* determine the initial positions from the found values */
01861 if ( init1 != -1 )
01862 {
01863 fitpar[2] = (float)init1 ;
01864 }
01865 fitpar[3] = 1. ;
01866
01867 for ( i = 0 ; i < NPAR ; i++ )
01868 {
01869 mpar[i] = 1 ;
01870 dervpar[i] = 0. ;
01871 }
01872
01873 xdim = XDIMA ;
01874 ndat = nel ;
01875 numpar = NPAR ;
01876 tol = TOLA ;
01877 lab = LABA ;
01878 its = ITSA ;
01879
01880 /* finally, do the least squares fit over the buffer data */
01881 if ( 0 > ( iters = lsqfit( xdat, &xdim, half_buffer -> data, wdat, &ndat, fitpar,
01882 dervpar, mpar, &numpar, &tol, &its, &lab )) )
01883 {
01884 cpl_msg_warning ("lsqfit:"," least squares fit failed, error no.: %d in slitlet: %d\n", iters, j) ;
01885 fitpar[2] = 0. ;
01886 }
01887 if ( fitpar[3] <=0. )
01888 {
01889 cpl_msg_warning ("fitSlitsBoltzSingleLine:"," fit failed due to negative width of boltzmann function in slitlet: %d\n", j) ;
01890 fitpar[2] = 0. ;
01891 }
01892
01893 pos = fitpar[2] ;
01894 if ( edge_ind == 1 )
01895 {
01896 pos -= (float)margin ;
01897 }
01898
01899 /*------------------------------------------------------------------------------------------
01900 * now discern the left and the right edge fit of the slitlets and associate the fit results
01901 * with the absolute positions in the image
01902 * consider the difference of the fitted slit position to the expected position
01903 * and decide wether the fit is taken or the expected value is taken.
01904 */
01905 if ( left_right == 0 )
01906 {
01907 /* take care of the column position of the fit boxes to get the absolute positions */
01908 if ( j == 0 )
01909 {
01910 slit_pos[0][0] = pos ;
01911 if ( slit_pos[0][0] - (int) slit_pos[0][0] == 0.)
01912 {
01913 slit_pos[0][0] = 0. ;
01914 }
01915 }
01916 else if ( j < k )
01917 {
01918 slit_pos[j][0] = pos + (float)edgeclean[j-1] - (float)margin ;
01919 if ( slit_pos[j][0] - (int) slit_pos[j][0] == 0.)
01920 {
01921 slit_pos[j][0] = 0. ;
01922 }
01923 }
01924 else
01925 {
01926 slit_pos[k][0] = pos + (float)edgeclean[k-1] - (float)margin ;
01927 if ( slit_pos[k][0] - (int) slit_pos[k][0] == 0.)
01928 {
01929 slit_pos[k][0] = 0. ;
01930 }
01931 }
01932 }
01933 else
01934 {
01935 /* take care of the column position of the fit boxes to get the absolute positions */
01936 if ( j == 0 )
01937 {
01938 slit_pos[0][1] = (float)(box_buffer->n_elements - 1) - pos ;
01939 if ( slit_pos[0][1] - (int) slit_pos[0][1] == 0.)
01940 {
01941 slit_pos[0][1] = 0. ;
01942 }
01943 }
01944 else if ( j < k )
01945 {
01946 slit_pos[j][1] = (float)(box_buffer->n_elements - 1) - pos
01947 + (float)edgeclean[j-1] - (float)margin ;
01948 if ( slit_pos[j][1] - (int) slit_pos[j][1] == 0.)
01949 {
01950 slit_pos[j][1] = 0. ;
01951 }
01952 }
01953 else
01954 {
01955 last_pos = (float)(edgeclean[k-1] - 1 + slit_length) ;
01956 slit_pos[k][1] = (float)(box_buffer->n_elements - 1) - pos
01957 + (float)edgeclean[k-1] - (float)margin ;
01958 if ( slit_pos[k][1] - (int) slit_pos[k][1] == 0.)
01959 {
01960 slit_pos[k][1] = 0. ;
01961 }
01962 }
01963 }
01964
01965 destroyVector ( half_buffer ) ;
01966 cpl_free( xdat ) ;
01967 cpl_free( wdat ) ;
01968 cpl_free( mpar ) ;
01969 }
01970 destroyVector ( box_buffer ) ;
01971 }
01972
01973 cpl_free( edge ) ;
01974 cpl_free( pos_row ) ;
01975 cpl_free( edgeclean ) ;
01976 cpl_free( dummyedge ) ;
01977 cpl_free( pos_rowclean ) ;
01978 return 0 ;
01979 }
01980
01981 /*----------------------------------------------------------------------------
01982 Function : fitSlitsBoltzWithEstimate()
01983 In : lineImage: emission line frame
01984 slit_pos: estimation array for the slitlet positions [min32][2]
01985 box_length: pixel length of the row box within the fit is done
01986 y_box: small box in spectral direction within the slitlet
01987 may lie.
01988 low_pos, high_pos: pixel positions in spectral direction between which the line
01989 should be located.
01990 Out : slit_pos: beginning and end position of the slitlets to
01991 sub-pixel accuracy
01992 0 if it worked,
01993 -1 if it failed,
01994 Job : fits the beginning and end position of the slitlets
01995 by using non-linear least square fitting of a Boltzmann function
01996 fits a Boltzmann function to the slitlet edges exposed and indicated
01997 by the brightest emission lines. The slitlet is searched within
01998 user given positions.
01999 The least squares fit is done by using a box smaller than
02000 the size of two slitlets
02001 ---------------------------------------------------------------------------*/
02002
02003 int fitSlitsBoltzWithEstimate ( OneImage * lineImage,
02004 float ** slit_pos,
02005 int box_length,
02006 float y_box,
02007 float diff_tol,
02008 int low_pos,
02009 int high_pos )
02010 {
02011 int position[lineImage->lx] ;
02012 Vector * box_buffer ;
02013 Vector * in_buffer ;
02014 int found_row ;
02015 int row, col ;
02016 int col_first, col_last ;
02017 int row_first, row_last ;
02018 int i, j, m, n ;
02019 int init1 ;
02020 int left_right ;
02021 int n_buf, shift ;
02022 int slit_length ;
02023 int iters, xdim, ndat ;
02024 int numpar, its ;
02025 int * mpar ;
02026 float * xdat, * wdat ;
02027 float tol, lab ;
02028 float fitpar[NPAR] ;
02029 float dervpar[NPAR] ;
02030 float minval, maxval ;
02031 float pos ;
02032 float new_pos ;
02033 int slitposition[SLITLENGTH] ;
02034 pixelvalue rowpos[SLITLENGTH] ;
02035
02036 slit_length = SLITLENGTH ; /* this setting is too much 64 */
02037 slit_length = N_SLITLETS ; /* this setting is better: 32 */
02038
02039 if ( NullImage == lineImage )
02040 {
02041 cpl_msg_error( "fitSlitsBoltzWithEstimate:"," no line image given!\n" ) ;
02042 return -1 ;
02043 }
02044
02045 if ( NULL == slit_pos )
02046 {
02047 cpl_msg_error( "fitSlitsBoltzWithEstimate:"," no position array allocated!\n" ) ;
02048 return -1 ;
02049 }
02050
02051 if ( box_length < 4 ||
02052 box_length > 2*slit_length )
02053 {
02054 cpl_msg_error( "fitSlitsBoltzWithEstimate:"," wrong fitting box length given!\n" ) ;
02055 return -1 ;
02056 }
02057
02058 if ( y_box <= 0. || y_box > 6. )
02059 {
02060 cpl_msg_error( "fitSlitsBoltzWithEstimate:"," wrong y box length given!\n" ) ;
02061 return -1 ;
02062 }
02063 if ( diff_tol <= 0. )
02064 {
02065 cpl_msg_error( "fitSlitsBoltzWithEstimate:"," wrong diff_tol given!\n" ) ;
02066 return -1 ;
02067 }
02068
02069 if ( low_pos >= high_pos || low_pos < 0 || high_pos <= 0 || high_pos > lineImage->ly )
02070 {
02071 cpl_msg_error( "fitSlitsBoltzWithEstimate:"," wrong user given search positions!\n" ) ;
02072 return -1 ;
02073 }
02074
02075 /* now search for the maximum between the given positions for each col */
02076 for ( col = 0 ; col < lineImage -> lx ; col++ )
02077 {
02078 found_row = -1 ;
02079 maxval = -FLT_MAX ;
02080 for ( row = low_pos ; row <= high_pos ; row++ )
02081 {
02082 if ( maxval < lineImage->data[col+row*lineImage->lx] )
02083 {
02084 maxval = lineImage->data[col+row*lineImage->lx] ;
02085 found_row = row ;
02086 }
02087 }
02088 if ( maxval > -FLT_MAX && found_row > low_pos )
02089 {
02090 position[col] = found_row ;
02091 }
02092 else
02093 {
02094 position[col] = 0 ;
02095 }
02096 }
02097
02098 /* --------------------------------------------------------------------------
02099 * now go through the slitlets, search along each column within a box with
02100 * half width y_box the maximum value and store these found values in a buffer
02101 */
02102 for ( j = 0 ; j < slit_length ; j++ )
02103 {
02104 /* now go through the columns and determine the slitlet positions by
02105 * calculating the median of the found positions
02106 */
02107 n = 0 ;
02108
02109 for ( col = nint(slit_pos[j][0])+ 1 ; col < nint(slit_pos[j][1]) -1 ; col++ )
02110 {
02111 rowpos[n] = (pixelvalue)position[col] ;
02112 n++ ;
02113 }
02114
02115 slitposition[j] = (int)median(rowpos, n) ;
02116 for ( left_right = 0 ; left_right <= 1 ; left_right++ )
02117 {
02118 init1 = 0 ;
02119 col_first = nint( slit_pos[j][left_right] ) - box_length/2 ;
02120 col_last = nint( slit_pos[j][left_right] ) + box_length/2 ;
02121 if ( col_first < 0 )
02122 {
02123 col_first = 0 ;
02124 init1 = 1 ;
02125 }
02126 if ( col_last > lineImage->lx )
02127 {
02128 col_last = lineImage->lx ;
02129 init1 = 1 ;
02130 }
02131 if ( col_last - col_first <= 0 )
02132 {
02133 cpl_msg_error( "fitSlitsBoltzWithEstimate:"," first and last column wrong!\n" ) ;
02134 return -1 ;
02135 }
02136 box_buffer = newVector( col_last - col_first ) ;
02137 m = 0 ;
02138
02139
02140 if ( left_right == 0 )
02141 {
02142 for( col = col_first ; col < col_last ; col++ )
02143 {
02144 row_first = slitposition[j] - nint(y_box) ;
02145 row_last = slitposition[j] + nint(y_box) ;
02146 if ( row_first < 0 )
02147 {
02148 row_first = 0 ;
02149 }
02150 if ( row_last >= lineImage->ly )
02151 {
02152 row_last = lineImage->ly - 1 ;
02153 }
02154 maxval = -FLT_MAX ;
02155 for ( row = row_first ; row <= row_last ; row++ )
02156 {
02157 if ( maxval < lineImage->data[col + lineImage->lx*row] )
02158 {
02159 maxval = lineImage -> data[col + lineImage->lx*row] ;
02160 }
02161 }
02162 box_buffer->data[m] = maxval ;
02163 m++ ;
02164 }
02165
02166 }
02167 else
02168 {
02169
02170 for( col = col_last-1 ; col >= col_first ; col-- )
02171 {
02172 row_first = slitposition[j] - nint(y_box) ;
02173 row_last = slitposition[j] + nint(y_box) ;
02174 if ( row_first < 0 )
02175 {
02176 row_first = 0 ;
02177 }
02178 if ( row_last >= lineImage->ly )
02179 {
02180 row_last = lineImage->ly - 1 ;
02181 }
02182 maxval = -FLT_MAX ;
02183 for ( row = row_first ; row <= row_last ; row++ )
02184 {
02185 if ( maxval < lineImage->data[col + lineImage->lx*row] )
02186 {
02187 maxval = lineImage -> data[col + lineImage->lx*row] ;
02188 }
02189 }
02190 box_buffer->data[m] = maxval ;
02191 m++ ;
02192 }
02193
02194 }
02195
02196 xdat = (float *) cpl_calloc( box_buffer->n_elements, sizeof (float) ) ;
02197 wdat = (float *) cpl_calloc( box_buffer->n_elements, sizeof (float) ) ;
02198 mpar = (int *) cpl_calloc( NPAR, sizeof (int) ) ;
02199
02200 /* set initial values for the fitting routine */
02201 minval = FLT_MAX ;
02202 maxval = -FLT_MAX ;
02203
02204 for ( i = 0 ; i < box_buffer->n_elements ; i++ )
02205 {
02206 xdat[i] = i ;
02207 wdat[i] = 1.0 ;
02208 if ( box_buffer -> data[i] < minval )
02209 {
02210 minval = box_buffer -> data[i] ;
02211 }
02212 if ( box_buffer -> data[i] > maxval )
02213 {
02214 maxval = box_buffer -> data[i] ;
02215 }
02216 }
02217 fitpar[0] = minval ;
02218 fitpar[1] = maxval ;
02219
02220 /*------------------------------------------------------------------
02221 * if we have too few left background values (at the image edges)
02222 * the left margin of the buffer to fit is filled with the minimal
02223 * values in order to get a good fit
02224 */
02225
02226
02227 if ( init1 == 1 )
02228 {
02229 n_buf = box_buffer->n_elements + box_length/2 ;
02230 in_buffer = newVector( n_buf ) ;
02231 for ( i = 0 ; i < box_length/2 ; i++ )
02232 {
02233 in_buffer -> data[i] = minval ;
02234 }
02235 shift = 0 ;
02236 for ( i = box_length/2 ; i < n_buf ; i++ )
02237 {
02238 in_buffer -> data[i] = box_buffer -> data[shift] ;
02239 shift++ ;
02240 }
02241 destroyVector ( box_buffer ) ;
02242 box_buffer = newVector ( n_buf ) ;
02243 for ( i = 0 ; i < n_buf ; i++ )
02244 {
02245 box_buffer -> data[i] = in_buffer -> data[i] ;
02246 }
02247 destroyVector ( in_buffer ) ;
02248 }
02249
02250 /* determine the initial positions from the found values */
02251 fitpar[2] = (float)box_buffer->n_elements/2. ;
02252 fitpar[3] = 1. ;
02253
02254 for ( i = 0 ; i < NPAR ; i++ )
02255 {
02256 mpar[i] = 1 ;
02257 dervpar[i] = 0. ;
02258 }
02259
02260 xdim = XDIMA ;
02261 ndat = box_buffer->n_elements ;
02262 numpar = NPAR ;
02263 tol = TOLA ;
02264 lab = LABA ;
02265 its = ITSA ;
02266
02267 /* finally, do the least squares fit over the buffer data */
02268 if ( 0 > ( iters = lsqfit( xdat, &xdim, box_buffer -> data, wdat, &ndat, fitpar,
02269 dervpar, mpar, &numpar, &tol, &its, &lab )) )
02270 {
02271 cpl_msg_warning ("lsqfit:"," least squares fit failed, error no.: %d in slitlet: %d\n", iters, j) ;
02272 destroyVector(box_buffer) ;
02273 cpl_free( xdat ) ;
02274 cpl_free( wdat ) ;
02275 cpl_free( mpar ) ;
02276 continue ;
02277 }
02278
02279 if ( fitpar[3] <=0. )
02280 {
02281 cpl_msg_warning ("fitSlitsBoltzWithEstimate:"," fit failed due to negative width of boltzmann function in slitlet: %d\n", j) ;
02282 destroyVector(box_buffer) ;
02283 cpl_free( xdat ) ;
02284 cpl_free( wdat ) ;
02285 cpl_free( mpar ) ;
02286 continue ;
02287 }
02288 pos = fitpar[2] ;
02289 if ( init1 == 1 )
02290 {
02291 pos -= (float)box_length/2. ;
02292 }
02293
02294 /*-------------------------------------------------------------
02295 * now compute the real slit positions using the guess positions
02296 * if the fit did not work the guess positions are taken
02297 * the same is done if the deviations are too big.
02298 */
02299 if ( pos != 0. )
02300 {
02301 if ( left_right == 0 )
02302 {
02303 new_pos = (float)col_first + pos ;
02304 }
02305 else
02306 {
02307 new_pos = (float)col_last-1 - pos ;
02308 }
02309 if ( fabs(new_pos - slit_pos[j][left_right]) < diff_tol )
02310 {
02311 slit_pos[j][left_right] = new_pos ;
02312 }
02313 else
02314 {
02315 cpl_msg_warning ("fitSlitsBoltzWithEstimate:"," deviation bigger than tolerance, take the estimated slitlet position in slitlet: %d\n", j) ;
02316 }
02317 }
02318
02319 cpl_free( xdat ) ;
02320 cpl_free( wdat ) ;
02321 cpl_free( mpar ) ;
02322 destroyVector ( box_buffer ) ;
02323
02324 }
02325 }
02326
02327 return 0 ;
02328 }
02329
02330 /*--------------------------------------------------------------------------*/
1.2.13.1 written by Dimitri van Heesch,
© 1997-2001