/* * This file is part of the ESO SINFONI Pipeline * Copyright (C) 2004,2005 European Southern Observatory * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, 51 Franklin St, Fifth Floor, Boston, MA 02111-1307 USA */ /***************************************************************************** * E.S.O. - VLT project * * * * who when what * -------- -------- ---------------------------------------------- * schreib 14/11/00 created */ /**@{*/ /*---------------------------------------------------------------------------*/ /** * @defgroup sinfo_absolute_utils routines to determine the absolute positions of the slitlets out of an emission line frame */ /*---------------------------------------------------------------------------*/ /************************************************************************ * NAME * sinfo_absolute.c - routines to determine the absolute positions * of the slitlets out of an emission line frame *------------------------------------------------------------------------ */ #ifdef HAVE_CONFIG_H # include #endif #include "sinfo_vltPort.h" #include "sinfo_absolute.h" #include "sinfo_recipes.h" /*---------------------------------------------------------------------------- * Defines *--------------------------------------------------------------------------*/ static float sqrarg ; #define SQR(a) (sqrarg = (a) , sqrarg*sqrarg) #define XDIMA 1 /* dimension of the x values */ #define TOLA 0.001 /* fitting tolerance */ #define LABA 0.1 /* lambda parameter */ #define ITSA 200 /* maximum number of iterations */ #define LABFACA 10.0 /* lambda step factor */ #define LABMAXA 1.0e+10 /* maximum value for lambda */ #define LABMINA 1.0e-10 /* minimum value for lambda */ #define NPAR 4 /* number of fit parameters */ /*---------------------------------------------------------------------------- * Local variables *--------------------------------------------------------------------------*/ static double chi1 ; /* old reduced chi-squared */ static double chi2 ; /* new reduced chi-squared */ static double labda ; /* mixing parameter */ static double vec[NPAR] ; /* correction vector */ static double matrix1[NPAR][NPAR] ; /* original matrix */ static double matrix2[NPAR][NPAR] ; /* inverse of matrix1 */ static int nfree ; /* number of free parameters */ static int parptr[NPAR] ; /* parameter pointer */ static float slopewidth ; /* initial value for fit parameter 5: width of linear slope */ /*---------------------------------------------------------------------------- * Functions private to this module *--------------------------------------------------------------------------*/ static int sinfo_new_inv_mat_edge (void) ; static void sinfo_new_get_mat_edge ( float * xdat, int * xdim, float * ydat, float * wdat, int * ndat, float * fpar, float * epar/*, int * npar*/ ) ; static int sinfo_new_get_vec_edge ( float * xdat, int * xdim, float * ydat, float * wdat, int * ndat, float * fpar, float * epar, int * npar ) ; float sinfo_new_hat2 ( float * xdat, float * parlist/*, int * npar, int * ndat*/ ); float sinfo_new_hat1 ( float * xdat, float * parlist/*, int * npar, int * ndat*/ ); void sinfo_new_hat_deriv2(float * xdat, float * parlist, float * dervs/*, int * npar*/ ); void sinfo_new_hat_deriv1( float * xdat, float * parlist, float * dervs/*, int * npar*/ ); int sinfo_new_fit_slits1( cpl_image * lineImage, FitParams ** par, float ** sinfo_slit_pos, int box_length, float y_box ); int sinfo_new_fit_slits( cpl_image * lineImage, FitParams ** par, float ** sinfo_slit_pos, int box_length, float y_box, float slope_width ); int sinfo_new_fit_slits2( cpl_image * lineImage, FitParams ** par, float ** sinfo_slit_pos, int box_length, float y_box, float diff_tol ); /*---------------------------------------------------------------------------- * Function codes *--------------------------------------------------------------------------*/ /** @name sinfo_new_edge() @input: position array xdat, parameter list parlist, number of parameters in the list npar The parameters are: @param xdat data sampling values @param parlist @return function value of a linear slope function: Intensity edge function ^ | /------parlist(3) | / | / | / parlist(2) |-----------/ | ^ ^ |----------|----|---------------> X axis parlist(0)=pos1 pos2=parlist(1) @note: parlist(0) pos1 parlist(1) pos2 parlist(2) intensity left parlist(3) intensity right A linear slope function is a function with a constant intensity value for xdat values smaller than pos1, linear increasing between pos1 and pos2, constant intensity value for xdat values greater than pos2 (see illustration above) return 0 in case of invalid input. @memo This function calculates the value of a slope function with parameters parlist at the position xdat */ float sinfo_new_edge ( float * xdat, float * parlist/*, int * npar, int * ndat*/ ) { float return_value ; float slope1 ; cpl_ensure(xdat , CPL_ERROR_NULL_INPUT, 0.); cpl_ensure(parlist , CPL_ERROR_NULL_INPUT, 0.); /* compute the slopes */ slope1 = ( parlist[3] - parlist[2] ) / ( parlist[1] - parlist[0] ) ; /* now build the hat function out of the parameters */ if ( xdat[0] <= parlist[0] ) { return_value = parlist[2] ; } else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[1] ) { return_value = (xdat[0] - parlist[0]) * slope1 + parlist[2] ; } else if ( xdat[0] > parlist[1] ) { return_value = parlist[3] ; } else { return_value = 0. ; } return return_value ; } /** @name sinfo_new_hat2 @memo calculates the value of a hat function with parameters parlist at the position xdat @param xdat position array @param parlist parameter list @return function value of a linear hat function. Intensity edge function ^ | parlist(6) /------\ | / \ | / \ | / \------parlist(7) parlist(4) |-----------/ | ^ ^ ^ ^ |----------|----|------|---|------> X axis pos1 pos2 pos3 pos4 @note -The parameter list values are: -# parlist(0): pos1 -# parlist(1): pos2 -# parlist(2): pos3 -# parlist(3): pos4 -# parlist(4): background left -# parlist(5): background right -# parlist(6): intensity left -# parlist(7): intensity right @doc This function returns a function with a constant background value for xdat values smaller than pos1, linear increasing between pos1 and pos2, linear between pos2 and pos3, linear decreasing between pos3 and pos4, and constant background value for xdat values greater than pos4. (see figure above) TODO: NOT USED */ float sinfo_new_hat2 ( float * xdat, float * parlist/*, int * npar, int * ndat*/ ) { float return_value ; float slope1, slope2, slope3 ; /* compute the slopes */ slope1 = ( parlist[6] - parlist[4] ) / ( parlist[1] - parlist[0] ) ; slope2 = ( parlist[7] - parlist[5] ) / ( parlist[3] - parlist[2] ) ; slope3 = ( parlist[7] - parlist[6] ) / ( parlist[2] - parlist[1] ) ; /* now build the hat function out of the parameters */ if ( xdat[0] <= parlist[0] ) { return_value = parlist[4] ; } else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[1] ) { return_value = (xdat[0] - parlist[0]) * slope1 + parlist[4] ; } else if ( xdat[0] > parlist[1] && xdat[0] <= parlist[2] ) { return_value = (xdat[0] - parlist[1]) * slope3 + parlist[6] ; } else if ( xdat[0] > parlist[2] && xdat[0] <= parlist[3] ) { return_value = (parlist[3] - xdat[0]) * slope2 + parlist[5] ; } else if ( xdat[0] > parlist[3] ) { return_value = parlist[5] ; } else { return_value = 0. ; } return return_value ; } /** @name sinfo_new_hat1 @memo calculates the value of a hat function with parameters parlist at the position xdat @param xdat position array @param parlist parameter list @return function value of a linear hat function. @note -The parameter list values are: -# parlist(0): pos1 -# parlist(1): pos2 -# parlist(2): pos3 -# parlist(3): pos4 -# parlist(4): background left -# parlist(5): background right -# parlist(6): intensity left -# parlist(7): intensity right Intensity edge function ^ | parlist(6) /------\ | / \ | / \ | / \------parlist(7) parlist(4) |-----------/ | ^ ^ ^ ^ |----------|----|------|---|------> X axis pos1 pos2 pos3 pos4 @doc This function returns a function with a constant background value for xdat values smaller than pos1, linear increasing between pos1 and pos1+slope_width, constant value intensity between pos1+slope_width and pos2-slope_width, linear decreasing between pos2-slope_width and pos2, and constant background value for xdat values greater than pos2. (see figure above) TODO: NOT USED */ float sinfo_new_hat1 ( float * xdat, float * parlist/*, int * npar, int * ndat*/ ) { float return_value=0 ; float slope1, slope2 ; /* take only positive values for the fit parameters */ /* compute the slopes */ slope1 = (parlist[4] - parlist[2]) / slopewidth ; slope2 = (parlist[4] - parlist[3]) / slopewidth ; /* now build the hat function out of the parameters */ if ( xdat[0] <= parlist[0] ) { return_value = parlist[2] ; } else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[0] + slopewidth ) { return_value = (xdat[0] - parlist[0]) * slope1 + parlist[2] ; } else if ( xdat[0] > parlist[0] + slopewidth && xdat[0] <= parlist[1] - slopewidth ) { return_value = parlist[4] ; } else if ( xdat[0] > parlist[1] - slopewidth && xdat[0] <= parlist[1] ) { return_value = (parlist[1] - xdat[0]) * slope2 + parlist[3] ; } else if ( xdat[0] > parlist[1] ) { return_value = parlist[3] ; } return return_value ; } /** @name sinfo_new_edge_deriv() @brief calculates the partial derivatives for a slope function with parameters parlist at position xdat @param xdat position array @param parlist parameter list @param dervs parameter derivatives (accuracies) list @return nothing (void) Intensity edge function ^ | /------parlist(3) | / | / | / parlist(2) |-----------/ | ^ ^ |----------|----|---------------> X axis parlist(0)=pos1 pos2=parlist(1) @note -The parameter list values are: parlist -# parlist[0]: pos1 -# parlist[1]: pos2 -# parlist[2]: intensity left -# parlist[3]: intensity right -The derivative values of a hat function at position xdat: dervs -# dervs[0]: partial derivative by pos1 -# dervs[1]: partial derivative by pos2 -# dervs[2]: partial derivative by intensity left -# dervs[3]: partial derivative by intensity right */ cpl_error_code sinfo_new_edge_deriv( float * xdat, float * parlist, float * dervs/*, int * npar*/ ) { float deriv1_slope1 ; cpl_ensure_code(xdat , CPL_ERROR_NULL_INPUT); cpl_ensure_code(parlist , CPL_ERROR_NULL_INPUT); cpl_ensure_code(dervs , CPL_ERROR_NULL_INPUT); /* compute the slopes */ deriv1_slope1 =( parlist[3] - parlist[2] ) / SQR(parlist[1] - parlist[0]) ; /* now build the hat derivatives out of the parameters */ if ( xdat[0] <= parlist[0] ) { dervs[0] = 0. ; dervs[1] = 0. ; dervs[2] = 1. ; dervs[3] = 0. ; } else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[1] ) { dervs[0] = ( xdat[0] - parlist[1] ) * deriv1_slope1 ; dervs[1] = ( parlist[0] - xdat[0] ) * deriv1_slope1 ; dervs[2] = ( parlist[0] - xdat[0] ) / ( parlist[1] - parlist[0] ) + 1.; dervs[3] = ( xdat[0] - parlist[0] ) / ( parlist[1] - parlist[0] ) ; } else if ( xdat[0] > parlist[1] ) { dervs[0] = 0. ; dervs[1] = 0. ; dervs[2] = 0. ; dervs[3] = 1. ; } return cpl_error_get_code(); } /** @name sinfo_new_hat_deriv2 @memo calculates the partial derivatives for a hat function with parameters parlist at position xdat @param xdat position array @param parlist parameter list @param dervs parameter derivatives (accuracies) list @return nothing (void) @note -The parameter list values are: parlist -# parlist[0]: pos1 -# parlist[1]: pos2 -# parlist[2]: pos3 -# parlist[3]: pos4 -# parlist[4]: background left -# parlist[5]: background right -# parlist[6]: intensity left -# parlist[7]: intensity right -The derivative values of a hat function at position xdat: dervs -# dervs[0]: partial derivative by pos1 -# dervs[1]: partial derivative by pos2 -# dervs[2]: partial derivative by pos3 -# dervs[3]: partial derivative by pos4 -# dervs[4]: partial derivative by background left -# dervs[5]: partial derivative by background right -# dervs[6]: partial derivative by intensity left -# dervs[7]: partial derivative by intensity right Intensity edge function ^ | parlist(6) /------\ | / \ | / \ | / \------parlist(7) parlist(4) |-----------/ | ^ ^ ^ ^ |----------|----|------|---|------> X axis pos1 pos2 pos3 pos4 TODO: NOT USED */ void sinfo_new_hat_deriv2(float * xdat, float * parlist, float * dervs/*, int * npar*/ ) { float deriv1_slope1 ; float deriv1_slope2 ; float deriv1_slope3 ; /* compute the slopes */ deriv1_slope1 = ( parlist[6] - parlist[4] ) / SQR(parlist[1] - parlist[0]); deriv1_slope2 = ( parlist[7] - parlist[5] ) / SQR(parlist[3] - parlist[2]); deriv1_slope3 = ( parlist[7] - parlist[6] ) / SQR(parlist[2] - parlist[1]); /* now build the hat derivatives out of the parameters */ if ( xdat[0] <= parlist[0] ) { dervs[0] = 0. ; dervs[1] = 0. ; dervs[2] = 0. ; dervs[3] = 0. ; dervs[4] = 1. ; dervs[5] = 0. ; dervs[6] = 0. ; dervs[7] = 0. ; } else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[1] ) { dervs[0] = ( xdat[0] - parlist[1] ) * deriv1_slope1 ; dervs[1] = ( parlist[0] - xdat[0] ) * deriv1_slope1 ; dervs[2] = 0. ; dervs[3] = 0. ; dervs[4] = ( parlist[0] - xdat[0] ) / ( parlist[1] - parlist[0] ) + 1.; dervs[5] = 0. ; dervs[6] = ( xdat[0] - parlist[0] ) / ( parlist[1] - parlist[0] ) ; dervs[7] = 0. ; } else if ( xdat[0] > parlist[1] && xdat[0] <= parlist[2] ) { dervs[0] = 0. ; dervs[1] = (xdat[0] - parlist[2]) * deriv1_slope3 ; dervs[2] = (parlist[1] - xdat[0]) * deriv1_slope3 ; dervs[3] = 0. ; dervs[4] = 0. ; dervs[5] = 0. ; dervs[6] = (parlist[1] - xdat[0]) / (parlist[2] - parlist[1]) + 1. ; dervs[7] = (xdat[0] - parlist[1]) / (parlist[2] - parlist[1]) ; } else if ( xdat[0] > parlist[2] && xdat[0] <= parlist[3] ) { dervs[0] = 0. ; dervs[1] = 0. ; dervs[2] = ( parlist[3] - xdat[0] ) * deriv1_slope2 ; dervs[3] = ( xdat[0] - parlist[2] ) * deriv1_slope2 ; dervs[4] = 0. ; dervs[5] = ( xdat[0] - parlist[3] ) / ( parlist[3] - parlist[2] ) + 1.; dervs[6] = 0. ; dervs[7] = ( parlist[3] - xdat[0] ) / ( parlist[3] - parlist[2] ) ; } else if ( xdat[0] > parlist[3] ) { dervs[0] = 0. ; dervs[1] = 0. ; dervs[2] = 0. ; dervs[3] = 0. ; dervs[4] = 0. ; dervs[5] = 1. ; dervs[6] = 0. ; dervs[7] = 0. ; } } /** @name sinfo_new_hat_deriv1 @memo calculates the partial derivatives for a hat function with parameters parlist at position xdat @param xdat position array @param parlist parameter list @param dervs parameter derivatives (accuracies) list @return nothing (void) @note -The parameter list values are: parlist -# parlist[0]: pos1 -# parlist[1]: pos2 -# parlist[2]: background1 -# parlist[3]: background2 -# parlist[4]: intensity -The derivative values of a hat function at position xdat: dervs -# dervs[0]: partial derivative by pos1 -# dervs[1]: partial derivative by pos2 -# dervs[2]: partial derivative by background1 -# dervs[3]: partial derivative by background2 -# dervs[4]: partial derivative by intensity Intensity edge function ^ | /------parlist(3) | / | / | / parlist(2) |-----------/ | ^ ^ |----------|----|---------------> X axis parlist(0)=pos1 pos2=parlist(1) TODO: NOT USED */ void sinfo_new_hat_deriv1( float * xdat, float * parlist, float * dervs/*, int * npar*/ ) { /* now build the hat derivatives out of the parameters */ if ( xdat[0] <= parlist[0] ) { dervs[0] = 0. ; dervs[1] = 0. ; dervs[2] = 1. ; dervs[3] = 0. ; dervs[4] = 0. ; } else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[0] + slopewidth ) { dervs[0] = ( parlist[2] - parlist[4] ) / slopewidth ; dervs[1] = 0. ; dervs[2] = (( parlist[0] - xdat[0] ) / slopewidth ) + 1. ; dervs[3] = 0. ; dervs[4] = ( xdat[0] - parlist[0] ) / slopewidth ; } else if ( xdat[0] > parlist[0] + slopewidth && xdat[0] <= parlist[1] - slopewidth ) { dervs[0] = 0. ; dervs[1] = 0. ; dervs[2] = 0. ; dervs[3] = 0. ; dervs[4] = 1. ; } else if ( xdat[0] > parlist[1] - slopewidth && xdat[0] <= parlist[1] ) { dervs[0] = 0. ; dervs[1] = ( parlist[4] - parlist[3] ) / slopewidth ; dervs[2] = 0. ; dervs[3] = (( xdat[0] - parlist[1] ) / slopewidth ) + 1. ; dervs[4] = ( parlist[1] - xdat[0] ) / slopewidth ; } else if ( xdat[0] > parlist[1] ) { dervs[0] = 0. ; dervs[1] = 0. ; dervs[2] = 0. ; dervs[3] = 1. ; dervs[4] = 0. ; } } /** @name sinfo_new_inv_mat_edge @brief calculates the inverse of matrix2. @doc The algorithm used is the Gauss-Jordan algorithm described in Stoer, Numerische Mathematik, 1. Teil. @return integer (0 if it worked, -6 if determinant is zero) */ static int sinfo_new_inv_mat_edge (void) { double even ; double hv[NPAR] ; double mjk ; int evin ; int i, j, k ; int per[NPAR] ; /* set permutation array */ for ( i = 0 ; i < nfree ; i++ ) { per[i] = i ; } for ( j = 0 ; j < nfree ; j++ ) /* in j-th column */ { double rowmax ; int row; /* determine largest element of a row */ rowmax = fabs ( matrix2[j][j] ) ; row = j ; for ( i = j + 1 ; i < nfree ; i++ ) { if ( fabs ( matrix2[i][j] ) > rowmax ) { rowmax = fabs( matrix2[i][j] ) ; row = i ; } } /* determinant is zero! */ if ( matrix2[row][j] == 0.0 ) { return -6 ; } /* if the largest element is not on the diagonal, then permute rows */ if ( row > j ) { for ( k = 0 ; k < nfree ; k++ ) { even = matrix2[j][k] ; matrix2[j][k] = matrix2[row][k] ; matrix2[row][k] = even ; } /* keep track of permutation */ evin = per[j] ; per[j] = per[row] ; per[row] = evin ; } /* modify column */ even = 1.0 / matrix2[j][j] ; for ( i = 0 ; i < nfree ; i++ ) { matrix2[i][j] *= even ; } matrix2[j][j] = even ; for ( k = 0 ; k < j ; k++ ) { mjk = matrix2[j][k] ; for ( i = 0 ; i < j ; i++ ) { matrix2[i][k] -= matrix2[i][j] * mjk ; } for ( i = j + 1 ; i < nfree ; i++ ) { matrix2[i][k] -= matrix2[i][j] * mjk ; } matrix2[j][k] = -even * mjk ; } for ( k = j + 1 ; k < nfree ; k++ ) { mjk = matrix2[j][k] ; for ( i = 0 ; i < j ; i++ ) { matrix2[i][k] -= matrix2[i][j] * mjk ; } for ( i = j + 1 ; i < nfree ; i++ ) { matrix2[i][k] -= matrix2[i][j] * mjk ; } matrix2[j][k] = -even * mjk ; } } /* finally, repermute the columns */ for ( i = 0 ; i < nfree ; i++ ) { for ( k = 0 ; k < nfree ; k++ ) { hv[per[k]] = matrix2[i][k] ; } for ( k = 0 ; k < nfree ; k++ ) { matrix2[i][k] = hv[k] ; } } /* all is well */ return 0 ; } /** @name sinfo_new_get_mat_edge @memo builds the sinfo_matrix @param xdat position array @param xdim factor of the indexes of the position array @param ydat real data @param wdat weights @param ndat number of data points @param fpar function parameters @param epar partial derivatives of the function //@param npar number of function parameters @return nothing (void) */ static void sinfo_new_get_mat_edge ( float * xdat, int * xdim, float * ydat, float * wdat, int * ndat, float * fpar, float * epar/*, int * npar*/ ) { double wd ; double yd ; int i, j, n ; for ( j = 0 ; j < nfree ; j++ ) { vec[j] = 0.0 ; /* zero sinfo_vector */ for ( i = 0 ; i<= j ; i++ ) /* zero matrix only on and below diagonal */ { matrix1[j][i] = 0.0 ; } } chi2 = 0.0 ; /* reset reduced chi-squared */ /* loop through data points */ for ( n = 0 ; n < (*ndat) ; n++ ) { double wn = wdat[n] ; if ( wn > 0.0 ) /* legal weight ? */ { yd = ydat[n] - sinfo_new_edge( &xdat[(*xdim) * n], fpar/*, npar, ndat*/ ) ; sinfo_new_edge_deriv( &xdat[(*xdim) * n], fpar, epar/*, npar */) ; chi2 += yd * yd * wn ; /* add to chi-squared */ for ( j = 0 ; j < nfree ; j++ ) { wd = epar[parptr[j]] * wn ; /* weighted derivative */ vec[j] += yd * wd ; /* fill sinfo_vector */ for ( i = 0 ; i <= j ; i++ ) /* fill sinfo_matrix */ { matrix1[j][i] += epar[parptr[i]] * wd ; } } } } } /** @name sinfo_new_get_vec_edge @memo Calculates the correction vector. @param xdat position array @param xdim factor of the indexes of the position array @param ydat real data @param wdat weights @param ndat number of data points @param fpar function parameters @param epar partial derivatives of the function @param npar number of function parameters @return integer ( # 0 if it had worked, # -5 or -7 if diagonal element is wrong, or # -6, if determinant is zero ) @doc Calculates the correction vector. The matrix has been built by get_mat_edge(), we only have to re-scale it for the current value of lambda. The matrix is rescaled so that the diagonal gets the value 1 + lambda. Next we calculate the inverse of the matrix and then the correction vector. */ static int sinfo_new_get_vec_edge ( float * xdat, int * xdim, float * ydat, float * wdat, int * ndat, float * fpar, float * epar, int * npar ) { double dy ; double mii ; double mji ; double mjj ; int i, j, n, r ; /* loop to modify and scale the sinfo_matrix */ for ( j = 0 ; j < nfree ; j++ ) { mjj = matrix1[j][j] ; if ( mjj <= 0.0 ) /* diagonal element wrong */ { return -5 ; } mjj = sqrt( mjj ) ; for ( i = 0 ; i < j ; i++ ) { mji = matrix1[j][i] / mjj / sqrt( matrix1[i][i] ) ; matrix2[i][j] = matrix2[j][i] = mji ; } matrix2[j][j] = 1.0 + labda ; /* scaled value on diagonal */ } if ( (r = sinfo_new_inv_mat_edge()) ) /* sinfo_invert sinfo_matrix inlace */ { return r ; } for ( i = 0 ; i < (*npar) ; i ++ ) { epar[i] = fpar[i] ; } /* loop to calculate correction sinfo_vector */ for ( j = 0 ; j < nfree ; j++ ) { double dj = 0.0 ; mjj = matrix1[j][j] ; if ( mjj <= 0.0) /* not allowed */ { return -7 ; } mjj = sqrt ( mjj ) ; for ( i = 0 ; i < nfree ; i++ ) { mii = matrix1[i][i] ; if ( mii <= 0.0 ) { return -7 ; } mii = sqrt( mii ) ; dj += vec[i] * matrix2[j][i] / mjj / mii ; } epar[parptr[j]] += dj ; /* new parameters */ } chi1 = 0.0 ; /* reset reduced chi-squared */ /* loop through the data points */ for ( n = 0 ; n < (*ndat) ; n++ ) { double wn = wdat[n] ; /* get weight */ if ( wn > 0.0 ) /* legal weight */ { dy = ydat[n] - sinfo_new_edge( &xdat[(*xdim) * n], epar /*, npar, ndat*/) ; chi1 += wdat[n] * dy * dy ; } } return 0 ; } /** @name sinfo_new_lsqfit_edge @brief Least square fit of a function to a set of data points @param xdat position, coordinates of data points. xdat is 2 dimensional: XDAT ( XDIM, NDAT ) @param xdim dimension of fit @param ydat data points @param wdat weights for data points @param ndat number of data points @param fpar on input contains initial estimates of the parameters for non-linear fits, on output the fitted parameters. @param epar contains estimates of the errors in fitted parameters @param mpar logical mask telling which parameters are free (non-zero) and which parameters are fixed (0) @param npar number of function parameters ( free + fixed ) @param tol relative tolerance. sinfo_lsqfit_edge stops when successive iterations fail to produce a decrement in reduced chi-squared less than tol. If tol is less than the minimum tolerance possible, tol will be set to this value. This means that maximum accuracy can be obtained by setting tol = 0.0. @param its maximum number of iterations @param lab mixing parameter, lab determines the initial weight of steepest descent method relative to the Taylor method lab should be a small value (i.e. 0.01). lab can only be zero when the partial derivatives are independent of the parameters. In fact in this case lab should be exactly equal to zero. @returns returns number of iterations needed to achieve convergence according to tol. When this number is negative, the fitting was not continued because a fatal error occurred: # -1 too many free parameters, maximum is 32 # -2 no free parameters # -3 not enough degrees of freedom # -4 maximum number of iterations too small to obtain a solution which satisfies tol. # -5 diagonal of sinfo_matrix contains elements which are zero # -6 determinant of the coefficient sinfo_matrix is zero # -7 square root of a negative number @doc this is a routine for making a least-squares fit of a function to a set of data points. The method used is described in: Marquardt, J.Soc.Ind.Appl.Math. 11. 431 (1963). This method is a mixture of the steepest descent method and the Taylor method. */ int sinfo_new_lsqfit_edge ( float * xdat, int * xdim, float * ydat, float * wdat, int * ndat, float * fpar, float * epar, int * mpar, int * npar, float * tol , int * its , float * lab ) { int i, n, r ; int itc ; /* fate of fit */ int found ; /* fit converged: 1, not yet converged: 0 */ int nuse ; /* number of usable data points */ double tolerance ; /* accuracy */ itc = 0 ; /* fate of fit */ found = 0 ; /* reset */ nfree = 0 ; /* number of free parameters */ nuse = 0 ; /* number of legal data points */ if ( *tol < (FLT_EPSILON * 10.0 ) ) { tolerance = FLT_EPSILON * 10.0 ; /* default tolerance */ } else { tolerance = *tol ; /* tolerance */ } labda = fabs( *lab ) * LABFACA ; /* start value for mixing parameter */ for ( i = 0 ; i < (*npar) ; i++ ) { if ( mpar[i] ) { if ( nfree > NPAR ) /* too many free parameters */ { return -1 ; } parptr[nfree++] = i ; /* a free parameter */ } } if (nfree == 0) /* no free parameters */ { return -2 ; } for ( n = 0 ; n < (*ndat) ; n++ ) { if ( wdat[n] > 0.0 ) /* legal weight */ { nuse ++ ; } } if ( nfree >= nuse ) { return -3 ; /* no degrees of freedom */ } if ( labda == 0.0 ) /* linear fit */ { /* initialize fpar array */ for ( i = 0 ; i < nfree ; fpar[parptr[i++]] = 0.0 ) ; sinfo_new_get_mat_edge(xdat,xdim,ydat,wdat,ndat,fpar,epar/*, npar */) ; r = sinfo_new_get_vec_edge ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ; if ( r ) /* error */ { return r ; } for ( i = 0 ; i < (*npar) ; i++ ) { fpar[i] = epar[i] ; /* save new parameters */ epar[i] = 0.0 ; /* and set errors to zero */ } chi1 = sqrt( chi1 / (double) (nuse - nfree) ) ; for ( i = 0 ; i < nfree ; i++ ) { if ( (matrix1[i][i] <= 0.0 ) || (matrix2[i][i] <= 0.0) ) { return -7 ; } epar[parptr[i]] = chi1 * sqrt( matrix2[i][i] ) / sqrt( matrix1[i][i] ) ; } } else /* non-linear fit */ { /*---------------------------------------------------------------- * the non-linear fit uses the steepest descent method in combination * with the Taylor method. The mixing of these methods is controlled * by lambda. In the outer loop ( called the iteration loop ) we build * the sinfo_matrix and calculate the correction sinfo_vector. In the * inner loop * (called the interpolation loop) we check whether we have obtained a * better solution than the previous one. If so, we leave the inner loop * else we increase lambda ( give more weight to the steepest descent * method) calculate the correction sinfo_vector and check again. * After the inner loop * we do a final check on the goodness of the fit and if this satisfies * the tolerance we calculate the errors of the fitted parameters. */ while ( !found ) /* iteration loop */ { if ( itc++ == (*its) ) /* increase iteration counter */ { return -4 ; } sinfo_new_get_mat_edge( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar*/ ) ; /*------------------------------------------------------------- * here we decrease lambda since we may assume that each iteration * brings us closer to the answer. */ if ( labda > LABMINA ) { labda = labda / LABFACA ; /* decrease lambda */ } r = sinfo_new_get_vec_edge ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ; if ( (int)fpar[1] - (int)fpar[0] <= 0 && fpar[1] / fpar[0] > 0. ) { fpar[1] += 1. ; continue ; } if ( r ) /* error */ { return r ; } while ( chi1 >= chi2 ) /* interpolation loop */ { /*----------------------------------------------------------- * The next statement is based on experience, not on the * mathematics of the problem. It is assumed that we have * reached convergence when the pure steepest descent method * does not produce a better solution. */ if ( labda > LABMAXA ) /* assume solution found */ { break ; } labda = labda * LABFACA ; /* increase mixing parameter */ r = sinfo_new_get_vec_edge ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ; if ( (int)fpar[1] - (int)fpar[0] <= 0 && fpar[1] / fpar[0] > 0. ) { fpar[1] += 1. ; continue ; } if ( r ) /* error */ { return r ; } } if ( labda <= LABMAXA ) /* save old parameters */ { for ( i = 0 ; i < *npar ; i++ ) { fpar[i] = epar[i] ; } } if ( (fabs( chi2 - chi1 ) <= (tolerance * chi1)) || (labda > LABMAXA) ) { /*------------------------------------------------------------ * we have a satisfying solution, so now we need to calculate * the correct errors of the fitted parameters. This we do by * using the pure Taylor method because we are very close to * the real solution. */ labda = LABMINA ; /* for Taylor solution */ sinfo_new_get_mat_edge ( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar */) ; r = sinfo_new_get_vec_edge ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ; if ( r ) /* error */ { return r ; } for ( i = 0 ; i < (*npar) ; i++ ) { epar[i] = 0.0 ; /* set error to zero */ } chi2 = sqrt ( chi2 / (double) (nuse - nfree) ) ; for ( i = 0 ; i < nfree ; i++ ) { if ( (matrix1[i][i] <= 0.0) || (matrix2[i][i] <= 0.0) ) { return -7 ; } epar[parptr[i]] = chi2 * sqrt( matrix2[i][i] ) / sqrt( matrix1[i][i] ) ; } found = 1 ; /* we found a solution */ } } } return itc ; /* return number of iterations */ } /** @name sinfo_new_fit_slits1 @memo fits the beginning and end position of the slitlets by using non-linear least square fitting of a hat function @param lineImage emission line frame @param par fit parameter data structure of fitted lines @param sinfo_slit_pos allocated dummy array for the slitlet positions [32][2] @param box_length pixel length of the row box within the fit is done @param y_box small box in spectral direction within the slitlet may lie. @param sinfo_slit_pos (out) beginning and end position of the slitlets to sub-pixel accuracy @returns # 0 if it worked, # -1 if there was no line image given, # -2 if there were no line fit parameters given, # -3 if there was no dummy array for the slit positions allocated # -4 if the given box length is impossible # -5 if the given y box length is impossible # -6 if there were no emission lines found in the first image columns # -7 if not all slitlets could be found # -8 if the least squares fit failed TODO: NOT USED */ int sinfo_new_fit_slits1( cpl_image * lineImage, FitParams ** par, float ** sinfo_slit_pos, int box_length, float y_box ) { float* position=NULL ; int * sinfo_edge, * edgeclean ; int * dummyedge ; int * pos_row, * pos_rowclean ; Vector * box_buffer ; int col ; int i, j, k, n, ed ; /* int init1; int init2 ; */ int line ; int slit_length ; int agreed ; int bad_line ; int margin ; int xdim, ndat ; int numpar, its ; float tol, lab ; float fitpar[2*NPAR] ; float dervpar[NPAR] ; int ilx=0; //int ily=0; float* pidata=NULL; if ( NULL == lineImage ) { sinfo_msg_error("no line image given!" ) ; return -1 ; } ilx=cpl_image_get_size_x(lineImage); //ily=cpl_image_get_size_y(lineImage); pidata=cpl_image_get_data_float(lineImage); slit_length = (int) sqrt (ilx) ; if ( NULL == par ) { sinfo_msg_error("no line fit parameters given!") ; return -2 ; } if ( NULL == sinfo_slit_pos ) { sinfo_msg_error("no position array allocated!") ; return -3 ; } if ( box_length < slit_length || box_length >= 3*slit_length ) { sinfo_msg_error("wrong fitting box length given!" ) ; return -4 ; } if ( y_box <= 0. || y_box > 3. ) { sinfo_msg_error("wrong y box length given!" ) ; return -5 ; } /* allocate memory for the edges and the row positon of the slitlets */ sinfo_edge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ; dummyedge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ; edgeclean = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ; pos_row = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ; pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ; /* ------------------------------------------------------------------------ * go through the first image columns and the fit parameters and find * the line with the highest intensity */ agreed = -1 ; bad_line = -1 ; while( agreed == -1 ) { int found = -1 ; float max_intensity = -FLT_MAX ; for ( col = 0 ; col < box_length ; col++ ) { for ( i = 0 ; i < par[0]->n_params ; i++ ) { if ( par[i]->column == col && par[i]->line != bad_line ) { if ( par[i]->fit_par[0] > max_intensity ) { if ( par[i]->fit_par[1] > 0. ) { max_intensity = par[i]->fit_par[0] ; found = i ; } } } } } /* -------------------------------------------------------------------- * check if the found line is usable and if the neighbouring line * have intensity on near rows in neighbouring slitlets */ float row_pos ; int column ; line = par[found]->line ; column = par[found]->column ; row_pos = par[found]->fit_par[2] ; if ( found >= 0 && max_intensity > 0. ) { for ( i = 0 ; i < par[0]->n_params ; i++ ) { if ( par[i]->line == line-1 && par[i]->column == column + slit_length ) { if ( par[i]->fit_par[2] <= (row_pos + y_box) && par[i]->fit_par[2] >= (row_pos - y_box) ) { bad_line = line ; } } } if ( bad_line != line ) { agreed = 1 ; break ; } } else { sinfo_msg_error("no emission line found in the first image columns"); cpl_free( sinfo_edge ) ; cpl_free( pos_row ) ; cpl_free( edgeclean ) ; cpl_free( dummyedge ) ; cpl_free( pos_rowclean ) ; return -6 ; } } if ( agreed == -1 ) { sinfo_msg_error("no emission line found in the first image columns") ; return -6 ; } /* now find and store the raw sinfo_edge positions of the found slitlet */ n = 0 ; ed = 0 ; position=cpl_calloc(ilx,sizeof(float)); for ( col = 0 ; col < ilx ; col++ ) { for ( i = 0 ; i < par[0]->n_params ; i++ ) { if ( par[i]->column == col && par[i] -> line == line ) { if ( par[i]->fit_par[0] > 0. && par[i]->fit_par[1] > 0. ) { position[n] = par[i]->fit_par[2] ; if ( n > 0 && fabs(position[n] - position[n-1]) > y_box ) { sinfo_edge[ed] = col ; pos_row[ed] = sinfo_new_nint( position[n-1] ) ; ed++ ; if ( col >= ilx - slit_length - 5 ) { pos_row[ed] = sinfo_new_nint( position[n] ) ; } } n++ ; } } } } if ( ed < (slit_length - 1) ) { sinfo_msg_error("not enough slitlets found") ; cpl_free( position ); return -7 ; } /* now find the clean sinfo_edge and row positions of the slitlets */ for ( i = 1 ; i <= ed ; i ++ ) { if (dummyedge[i-1] != -1) { dummyedge[i-1] = sinfo_edge[i-1] ; } if ( (sinfo_edge[i] - sinfo_edge[i-1]) < slit_length - 3 || (sinfo_edge[i] - sinfo_edge[i-1]) > slit_length + 3 ) { dummyedge[i] = -1 ; } if ( (sinfo_edge[i+1] - sinfo_edge[i]) < slit_length - 3 || (sinfo_edge[i+1] - sinfo_edge[i]) > slit_length + 3 ) { dummyedge[i+1] = -1 ; } } k = 0 ; for ( i = 0 ; i < ed ; i++ ) { if ( dummyedge[i] != -1 && dummyedge[i] != 0 ) { edgeclean[k] = dummyedge[i] ; pos_rowclean[k] = pos_row[i] ; k++ ; if( edgeclean[k-1] > (ilx - slit_length - 6 ) ) { pos_rowclean[k] = pos_row[ed] ; } } } if ( k != slit_length - 1 ) { sinfo_msg_error("not enough clean slitlets found") ; return -7 ; } /* determine the margins of the fitting box outside the slitlets */ margin = ( box_length - slit_length ) / 2 ; /* now go through the slitlets and store the intensity in a buffer sinfo_vector */ for ( j = 0 ; j < k ; j++ ) { int m = 0 ; if ( j == 0 ) { box_buffer = sinfo_new_vector( edgeclean[0] + margin ) ; for( col = 0 ; col < edgeclean[0] + margin ; col++ ) { box_buffer->data[m] = pidata[col + ilx*pos_rowclean[0]] ; m++ ; } fitpar[0] = 3. ; fitpar[1] = 5. ; fitpar[2] = (float) edgeclean[0] - 1. ; fitpar[3] = (float) edgeclean[0] + 1. ; } else if ( j < k - 1 ) { box_buffer = sinfo_new_vector( edgeclean[j] - edgeclean[j-1] + 2*margin ) ; for ( col = edgeclean[j - 1] - margin ; col < edgeclean[j] + margin ; col++ ) { box_buffer->data[m] = pidata[col + ilx*pos_rowclean[j]] ; m++ ; } fitpar[0] = (float) margin - 1. ; fitpar[1] = (float) margin + 1. ; fitpar[2] = (float) (edgeclean[j] - edgeclean[j-1] + margin) - 1. ; fitpar[3] = (float) (edgeclean[j] - edgeclean[j-1] + margin) + 1. ; } /*else if ( j == k - 1 )*/ else { box_buffer = sinfo_new_vector( ilx - edgeclean[k-1] + margin ) ; for ( col = edgeclean[k - 1] - margin ; col < ilx ; col++ ) { box_buffer->data[m] = pidata[col + ilx*pos_rowclean[k]] ; m++ ; } fitpar[0] = (float) margin - 1. ; fitpar[1] = (float) margin + 1. ; fitpar[2] = (float) (ilx - edgeclean[k-1] + margin) - 3. ; fitpar[3] = (float) (ilx - edgeclean[k-1] + margin) - 1. ; } float* xdat = (float *) cpl_calloc(box_buffer -> n_elements,sizeof (float) ) ; float* wdat = (float *) cpl_calloc(box_buffer -> n_elements,sizeof (float) ) ; int* mpar = (int *) cpl_calloc(NPAR, sizeof (int) ) ; /* set initial values for the fitting routine */ float minval = FLT_MAX ; float maxval = -FLT_MAX ; for ( i = 0 ; i < box_buffer->n_elements ; i++ ) { xdat[i] = i ; wdat[i] = 1.0 ; if ( box_buffer -> data[i] < minval ) { minval = box_buffer -> data[i] ; } if ( box_buffer -> data[i] > maxval ) { maxval = box_buffer -> data[i] ; } } fitpar[4] = minval ; fitpar[5] = minval ; fitpar[6] = maxval ; fitpar[7] = maxval ; /* search for both positions of the half intensity of the hat within the buffer */ /* init1 = -1 ; init2 = -1 ; */ for ( i = 0 ; i < box_buffer->n_elements ; i++ ) { if ( box_buffer -> data[i] >= ( maxval - minval ) / 2. ) { //init1 = i ; break ; } } for ( i = box_buffer->n_elements - 1 ; i >= 0 ; i-- ) { if ( box_buffer -> data[i] >= ( maxval + minval ) / 2. ) { //init2 = i ; break ; } } /* determine the initial positions from the found values */ /* if ( init1 != -1 ) { fitpar[0] = init1 - 1. ; fitpar[1] = init1 + 1. ; } if ( init2 != -1 ) { fitpar[2] = init2 - 1. ; fitpar[3] = init2 + 1. ; } */ for ( i = 0 ; i < NPAR ; i++ ) { mpar[i] = 1 ; dervpar[i] = 0. ; } xdim = XDIMA ; ndat = box_buffer -> n_elements ; numpar = NPAR ; tol = TOLA ; lab = LABA ; its = ITSA ; int iters; /* finally, do the least squares fit over the buffer data */ if ( 0 > ( iters = sinfo_new_lsqfit_edge( xdat, &xdim, box_buffer -> data, wdat, &ndat, fitpar, dervpar, mpar, &numpar, &tol, &its, &lab )) ) { sinfo_msg_warning("least squares fit failed, error no.: %d", iters) ; return -8 ; } /* take care of the column position of the fit boxes to get the absolute positions */ if ( j == 0 ) { sinfo_slit_pos[0][0] = (fitpar[0] + fitpar[1]) / 2. ; sinfo_slit_pos[0][1] = (fitpar[2] + fitpar[3]) / 2. ; } else { sinfo_slit_pos[j][0] = (fitpar[0] + fitpar[1]) / 2. + (float)edgeclean[j-1] - (float)margin ; sinfo_slit_pos[j][1] = (fitpar[2] + fitpar[3]) / 2. + (float)edgeclean[j-1] - (float)margin ; } sinfo_slit_pos[k][0] = (fitpar[0] + fitpar[1]) / 2. + (float)edgeclean[k-1] - (float)margin ; sinfo_slit_pos[k][1] = (fitpar[2] + fitpar[3]) / 2. + (float)edgeclean[k-1] - (float)margin ; sinfo_new_destroy_vector ( box_buffer ) ; cpl_free( xdat ) ; cpl_free( wdat ) ; cpl_free( mpar ) ; } cpl_free( sinfo_edge ) ; cpl_free( pos_row ) ; cpl_free( edgeclean ) ; cpl_free( dummyedge ) ; cpl_free( pos_rowclean ) ; cpl_free( position ); return 0 ; } /** @name sinfo_new_fit_slits @memo fits the beginning and end position of the slitlets by using non-linear least square fitting of a hat function @param lineImage emission line frame @param par fit parameter data structure of fitted lines @param sinfo_slit_pos allocated dummy array for the slitlet positions [32][2] @param box_length pixel length of the row box within the fit is done @param y_box small box in spectral direction within the slitlet may lie. @param slope_width width of linear slope of the hat function, must be positive @param sinfo_slit_pos beginning and end position of the slitlets to sub-pixel accuracy @returns # 0 if it worked, # -1 if there was no line image given, # -2 if there were no line fit parameters given, # -3 if there was no dummy array for the slit positions # allocated # -4 if the given box length is impossible # -5 if the given y box length is impossible # -6 if the given width of the linear slope is wrong # -7 if there were no emission lines found in the first image columns # -8 if not all slitlets could be found # -9 if the least squares fit failed TODO: NOT USED */ int sinfo_new_fit_slits( cpl_image * lineImage, FitParams ** par, float ** sinfo_slit_pos, int box_length, float y_box, float slope_width ) { float* position=NULL ; int * sinfo_edge, * edgeclean ; int * dummyedge ; int * pos_row, * pos_rowclean ; Vector * box_buffer ; int col ; int i, j, k, n, ed ; int line ; int slit_length ; int agreed ; int bad_line ; int margin ; int xdim, ndat ; int numpar, its ; float tol, lab ; float fitpar[NPAR+1] ; float dervpar[NPAR] ; int ilx=0; //int ily=0; float* pidata=NULL; if ( NULL == lineImage ) { sinfo_msg_error("no line image given!" ) ; return -1 ; } ilx=cpl_image_get_size_x(lineImage); //ily=cpl_image_get_size_y(lineImage); pidata=cpl_image_get_data_float(lineImage); slit_length = (int) sqrt (ilx) ; if ( NULL == par ) { sinfo_msg_error("no line fit parameters given!" ) ; return -2 ; } if ( NULL == sinfo_slit_pos ) { sinfo_msg_error("no position array allocated!" ) ; return -3 ; } if ( box_length < slit_length || box_length >= 3*slit_length ) { sinfo_msg_error("wrong fitting box length given!" ) ; return -4 ; } if ( y_box <= 0. || y_box > 3. ) { sinfo_msg_error("wrong y box length given!" ) ; return -5 ; } if ( slope_width <= 0. ) { sinfo_msg_error("wrong width of linear slope given!") ; return -6 ; } /* initialise module global variable slopewidth */ slopewidth = slope_width ; /* allocate memory for the edges and the row position of the slitlets */ sinfo_edge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ; dummyedge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ; edgeclean = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ; pos_row = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ; pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ; /* ------------------------------------------------------------------------ * go through the first image columns and the fit parameters and find * the line with the highest intensity */ agreed = -1 ; bad_line = -1 ; while( agreed == -1 ) { int found = -1 ; float max_intensity = -FLT_MAX ; for ( col = 0 ; col < box_length ; col++ ) { for ( i = 0 ; i < par[0]->n_params ; i++ ) { if ( par[i]->column == col && par[i]->line != bad_line ) { if ( par[i]->fit_par[0] > max_intensity ) { if ( par[i]->fit_par[1] > 0. ) { max_intensity = par[i]->fit_par[0] ; found = i ; } } } } } /* -------------------------------------------------------------------- * check if the found line is usable and if the neighbouring line * have intensity on near rows in neighbouring slitlets */ float row_pos ; int column ; line = par[found]->line ; column = par[found]->column ; row_pos = par[found]->fit_par[2] ; if ( found >= 0 && max_intensity > 0. ) { for ( i = 0 ; i < par[0]->n_params ; i++ ) { if ( par[i]->line == line-1 && par[i]->column == column + slit_length ) { if ( par[i]->fit_par[2] <= (row_pos + y_box) && par[i]->fit_par[2] >= (row_pos - y_box) ) { bad_line = line ; } } } if ( bad_line != line ) { agreed = 1 ; break ; } } else { sinfo_msg_error("no emission line found in the first image columns"); cpl_free( sinfo_edge ) ; cpl_free( pos_row ) ; cpl_free( edgeclean ) ; cpl_free( dummyedge ) ; cpl_free( pos_rowclean ) ; return -7 ; } } if ( agreed == -1 ) { sinfo_msg_error("no emission line found in the first image columns") ; return -7 ; } /* now find and store the raw sinfo_edge positions of the found slitlet */ n = 0 ; ed = 0 ; position=cpl_calloc(ilx,sizeof(float)) ; for ( col = 0 ; col < ilx ; col++ ) { for ( i = 0 ; i < par[0]->n_params ; i++ ) { if ( par[i]->column == col && par[i] -> line == line ) { if ( par[i]->fit_par[0] > 0. && par[i]->fit_par[1] > 0. ) { position[n] = par[i]->fit_par[2] ; if ( n > 0 && fabs(position[n] - position[n-1]) > y_box ) { sinfo_edge[ed] = col ; pos_row[ed] = sinfo_new_nint( position[n-1] ) ; ed++ ; if ( col >= ilx - slit_length - 5 ) { pos_row[ed] = sinfo_new_nint( position[n] ) ; } } n++ ; } } } } if ( ed < (slit_length - 1) ) { sinfo_msg_error("not enough slitlets found") ; cpl_free( position ); return -8 ; } /* now find the clean sinfo_edge and row positions of the slitlets */ for ( i = 1 ; i <= ed ; i ++ ) { if (dummyedge[i-1] != -1) { dummyedge[i-1] = sinfo_edge[i-1] ; } if ( (sinfo_edge[i] - sinfo_edge[i-1]) < slit_length - 3 || (sinfo_edge[i] - sinfo_edge[i-1]) > slit_length + 3 ) { dummyedge[i] = -1 ; } if ( (sinfo_edge[i+1] - sinfo_edge[i]) < slit_length - 3 || (sinfo_edge[i+1] - sinfo_edge[i]) > slit_length + 3 ) { dummyedge[i+1] = -1 ; } } k = 0 ; for ( i = 0 ; i < ed ; i++ ) { if ( dummyedge[i] != -1 && dummyedge[i] != 0 ) { edgeclean[k] = dummyedge[i] ; pos_rowclean[k] = pos_row[i] ; k++ ; if( edgeclean[k-1] > (ilx - slit_length - 6 ) ) { pos_rowclean[k] = pos_row[ed] ; } } } if ( k != slit_length - 1 ) { sinfo_msg_error ("not enough clean slitlets found") ; return -7 ; } /* determine the margins of the fitting box outside the slitlets */ margin = ( box_length - slit_length ) / 2 ; /* now go through the slitlets and store the intensity in a buffer sinfo_vector */ for ( j = 0 ; j < k ; j++ ) { int m = 0 ; if ( j == 0 ) { box_buffer = sinfo_new_vector( edgeclean[0] + margin ) ; for( col = 0 ; col < edgeclean[0] + margin ; col++ ) { box_buffer->data[m] = pidata[col + ilx*pos_rowclean[0]] ; m++ ; } /* initial values for the fitted positions */ fitpar[0] = 1. ; fitpar[1] = (float)edgeclean[0] ; } else if ( j < k - 1 ) { box_buffer = sinfo_new_vector( edgeclean[j] - edgeclean[j-1] + 2*margin ) ; for ( col = edgeclean[j - 1] - margin ; col < edgeclean[j] + margin ; col++ ) { box_buffer->data[m] = pidata[col + ilx*pos_rowclean[j]] ; m++ ; } /* initial values for the fitted positions */ fitpar[0] = (float)margin ; fitpar[1] = (float)(edgeclean[j] - edgeclean[j-1] + margin ) ; } /*else if ( j == k - 1 )*/ else { box_buffer = sinfo_new_vector( ilx - edgeclean[k-1] + margin ) ; for ( col = edgeclean[k - 1] - margin ; col < ilx ; col++ ) { box_buffer->data[m] = pidata[col + ilx*pos_rowclean[k]] ; m++ ; } /* initial values for the fitted positions */ fitpar[0] = (float)margin ; fitpar[1] = (float)(m - 1) ; } float* xdat=(float *) cpl_calloc( box_buffer -> n_elements, sizeof (float) ) ; float* wdat=(float *) cpl_calloc( box_buffer -> n_elements, sizeof (float) ) ; int* mpar=(int *) cpl_calloc( NPAR, sizeof (int) ) ; /* set initial values for the fitting routine */ float minval, maxval ; minval = FLT_MAX ; maxval = -FLT_MAX ; for ( i = 0 ; i < box_buffer->n_elements ; i++ ) { xdat[i] = i ; wdat[i] = 1.0 ; if ( box_buffer -> data[i] < minval ) { minval = box_buffer -> data[i] ; } if ( box_buffer -> data[i] > maxval ) { maxval = box_buffer -> data[i] ; } } for ( i = 0 ; i < NPAR ; i++ ) { mpar[i] = 1 ; dervpar[i] = 0. ; } xdim = XDIMA ; ndat = box_buffer -> n_elements ; numpar = NPAR ; tol = TOLA ; lab = LABA ; its = ITSA ; fitpar[2] = minval ; fitpar[3] = minval ; fitpar[4] = maxval ; /* finally, do the least squares fit over the buffer data */ int iters; if ( 0 > ( iters = sinfo_new_lsqfit_edge( xdat, &xdim, box_buffer -> data, wdat, &ndat, fitpar, dervpar, mpar, &numpar, &tol, &its, &lab )) ) { sinfo_msg_warning("least squares fit failed, error no.: %d", iters) ; return -9 ; } /* take care of the column position of the fit boxes to get the absolute positions */ if ( j == 0 ) { sinfo_slit_pos[0][0] = fitpar[0] + slopewidth/2. ; sinfo_slit_pos[0][1] = fitpar[1] - slopewidth/2. ; } else { sinfo_slit_pos[j][0] = fitpar[0] + slopewidth/2. + (float)edgeclean[j-1] - (float)margin ; sinfo_slit_pos[j][1] = fitpar[1] - slopewidth/2. + (float)edgeclean[j-1] - (float)margin ; } sinfo_slit_pos[k][0] = fitpar[0] + slopewidth/2. + (float)edgeclean[k-1] - (float)margin ; sinfo_slit_pos[k][1] = fitpar[1] - slopewidth/2. + (float)edgeclean[k-1] - (float)margin ; sinfo_new_destroy_vector ( box_buffer ) ; cpl_free( xdat ) ; cpl_free( wdat ) ; cpl_free( mpar ) ; } cpl_free( sinfo_edge ) ; cpl_free( pos_row ) ; cpl_free( edgeclean ) ; cpl_free( dummyedge ) ; cpl_free( pos_rowclean ) ; cpl_free( position ); return 0 ; } /** @name sinfo_new_fit_slits2 @memo fits the beginning and end position of the slitlets by using non-linear least square fitting of a step function fits a step function to the slitlet edges exposed and indicated by the brightest emission lines. @param lineImage emission line frame @param par fit parameter data structure of fitted lines @param sinfo_slit_pos allocated dummy array for the slitlet positions [32][2] @param box_length pixel length of the row box within the fit is done @param y_box small box in spectral direction within the slitlet may lie. @param diff_tol maximum tolerable difference of the resulting fit position with respect to the expected position. If difference is greater the expected position is taken. @param sinfo_slit_pos beginning and end position of the slitlets to sub-pixel accuracy # 0 if it worked, # -1 if there was no line image given, # -2 if there were no line fit parameters given, # -3 if there was no dummy array for the slit positions # allocated # -4 if the given box length is impossible # -5 if the given y box length is impossible # -6 if the given difference tolerance is too small # -7 if there were no emission lines found in the first image columns # -8 if not all slitlets could be found fits the beginning and end position of the slitlets by using non-linear least square fitting of a step function fits a step function to the slitlet edges exposed and indicated by the brightest emission lines. To achieve this, the fit parameters are used to find the brightest emission line and to get its position for each column. The least squares fit is done by using a box bigger than the size of one slitlet and divides the box into two parts for both edges within the fit function is shifted. TODO: NOT USED */ int sinfo_new_fit_slits2( cpl_image * lineImage, FitParams ** par, float ** sinfo_slit_pos, int box_length, float y_box, float diff_tol ) { float* position=NULL ; int * sinfo_edge, * edgeclean ; int * dummyedge ; int * pos_row, * pos_rowclean ; Vector * box_buffer ; Vector * half_buffer ; int col ; int i, j, k, n, ed ; int init1 ; int line ; int nel, n_right, left_right ; int slit_length ; int agreed ; int bad_line ; int margin ; int iters, xdim, ndat ; int numpar, its ; int * mpar ; float * xdat, * wdat ; float tol, lab ; float fitpar[NPAR] ; float dervpar[NPAR] ; float minval, maxval ; float pos, last_pos ; int ilx=0; //int ily=0; float* pidata=NULL; if ( NULL == lineImage ) { sinfo_msg_error("no line image given!" ) ; return -1 ; } ilx=cpl_image_get_size_x(lineImage); //ily=cpl_image_get_size_y(lineImage); pidata=cpl_image_get_data_float(lineImage); slit_length = (int) sqrt (ilx) ; if ( NULL == par ) { sinfo_msg_error("no line fit parameters given!" ) ; return -2 ; } if ( NULL == sinfo_slit_pos ) { sinfo_msg_error("no position array allocated!" ) ; return -3 ; } if ( box_length < slit_length || box_length >= 3*slit_length ) { sinfo_msg_error("wrong fitting box length given!" ) ; return -4 ; } if ( y_box <= 0. || y_box > 3. ) { sinfo_msg_error("wrong y box length given!" ) ; return -5 ; } if ( diff_tol < 1. ) { sinfo_msg_error("diff_tol too small!" ) ; return -6 ; } /* allocate memory for the edges and the row position of the slitlets */ sinfo_edge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ; dummyedge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ; edgeclean = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ; pos_row = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ; pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ; /* ------------------------------------------------------------------------ * go through the first image columns and the fit parameters and find * the line with the highest intensity */ agreed = -1 ; bad_line = -1 ; while( agreed == -1 ) { int found = -1 ; float max_intensity = -FLT_MAX ; for ( col = 0 ; col < box_length ; col++ ) { for ( i = 0 ; i < par[0]->n_params ; i++ ) { if ( par[i]->column == col && par[i]->line != bad_line ) { if ( par[i]->fit_par[0] > max_intensity ) { if ( par[i]->fit_par[1] > 0. ) { max_intensity = par[i]->fit_par[0] ; found = i ; } } } } } /* -------------------------------------------------------------------- * check if the found line is usable and if the neighbouring line * have intensity on near rows in neighbouring slitlets */ int column ; line = par[found]->line ; column = par[found]->column ; float row_pos = par[found]->fit_par[2] ; if ( found >= 0 && max_intensity > 0. ) { for ( i = 0 ; i < par[0]->n_params ; i++ ) { if ( par[i]->line == line-1 && par[i]->column == column + slit_length ) { if ( par[i]->fit_par[2] <= (row_pos + y_box) && par[i]->fit_par[2] >= (row_pos - y_box) ) { bad_line = line ; } } } if ( bad_line != line ) { agreed = 1 ; break ; } } else { sinfo_msg_error("no emission line found in the first image columns"); cpl_free( sinfo_edge ) ; cpl_free( pos_row ) ; cpl_free( edgeclean ) ; cpl_free( dummyedge ) ; cpl_free( pos_rowclean ) ; return -7 ; } } if ( agreed == -1 ) { sinfo_msg_error("no emission line found in the first image columns") ; return -7 ; } /* now find and store the raw sinfo_edge positions of the found slitlet */ n = 0 ; ed = 0 ; position=cpl_calloc(ilx,sizeof(float)) ; for ( col = 0 ; col < ilx ; col++ ) { for ( i = 0 ; i < par[0]->n_params ; i++ ) { if ( par[i]->column == col && par[i]->line == line ) { if ( par[i]->fit_par[0] > 0. && par[i]->fit_par[1] > 0. ) { position[n] = par[i]->fit_par[2] ; if ( n > 0 && fabs(position[n] - position[n-1]) > y_box ) { sinfo_edge[ed] = col ; pos_row[ed] = sinfo_new_nint( position[n-1] ) ; ed++ ; if ( col >= ilx - slit_length - 5 ) { pos_row[ed] = sinfo_new_nint( position[n] ) ; } } n++ ; } } } } if ( ed < (slit_length - 1) ) { sinfo_msg_error("not enough slitlets found") ; cpl_free(position); return -8 ; } /* now find the clean sinfo_edge and row positions of the slitlets */ for ( i = 1 ; i <= ed ; i ++ ) { if (dummyedge[i-1] != -1) { dummyedge[i-1] = sinfo_edge[i-1] ; } if ( (sinfo_edge[i] - sinfo_edge[i-1]) < slit_length - 3 || (sinfo_edge[i] - sinfo_edge[i-1]) > slit_length + 3 ) { dummyedge[i] = -1 ; } if ( (sinfo_edge[i+1] - sinfo_edge[i]) < slit_length - 3 || (sinfo_edge[i+1] - sinfo_edge[i]) > slit_length + 3 ) { dummyedge[i+1] = -1 ; } } k = 0 ; for ( i = 0 ; i < ed ; i++ ) { if ( dummyedge[i] != -1 && dummyedge[i] != 0 ) { edgeclean[k] = dummyedge[i] ; pos_rowclean[k] = pos_row[i] ; k++ ; if( edgeclean[k-1] > (ilx - slit_length - 6 ) ) { pos_rowclean[k] = pos_row[ed] ; } } } if ( k != slit_length - 1 ) { sinfo_msg_error("not enough clean slitlets found") ; return -7 ; } /* determine the margins of the fitting box outside the slitlets */ margin = ( box_length - slit_length ) / 2 ; /* now go through the slitlets and store the intensity in a buffer sinfo_vector */ for ( j = 0 ; j <= k ; j++ ) { int m = 0 ; if ( j == 0 ) { box_buffer = sinfo_new_vector( edgeclean[0] + margin ) ; for( col = 0 ; col < edgeclean[0] + margin ; col++ ) { box_buffer->data[m] = pidata[col +ilx*pos_rowclean[0]] ; m++ ; } } else if ( j < k ) { box_buffer = sinfo_new_vector( edgeclean[j] - edgeclean[j-1] + 2*margin ) ; for ( col = edgeclean[j - 1] - margin ; col < edgeclean[j] + margin ; col++ ) { box_buffer->data[m] = pidata[col + ilx*pos_rowclean[j]] ; m++ ; } } else { box_buffer = sinfo_new_vector( ilx - edgeclean[k-1] + margin ) ; for ( col = edgeclean[k - 1] - margin ; col < ilx ; col++ ) { box_buffer->data[m] = pidata[col + ilx*pos_rowclean[k]] ; m++ ; } } for ( left_right = 0 ; left_right <= 1 ; left_right++ ) { nel = 0 ; if ( left_right == 0 ) { nel = box_buffer -> n_elements / 2 ; } else { if ( box_buffer -> n_elements % 2 == 0 ) { nel = box_buffer -> n_elements / 2 ; } else { nel = box_buffer -> n_elements / 2 + 1 ; } } /* now split the buffer in the midth in a left and right part for fitting */ half_buffer = sinfo_new_vector( nel ) ; if ( left_right == 0 ) { for ( i = 0 ; i < nel ; i++ ) { half_buffer -> data[i] = box_buffer -> data[i] ; } } else { n_right = 0 ; for ( i = box_buffer -> n_elements - 1 ; i >= box_buffer -> n_elements - nel ; i-- ) { half_buffer -> data[n_right] = box_buffer -> data[i] ; n_right++ ; } } xdat = (float *) cpl_calloc( nel, sizeof (float) ) ; wdat = (float *) cpl_calloc( nel, sizeof (float) ) ; mpar = (int *) cpl_calloc( NPAR, sizeof (int) ) ; /* set initial values for the fitting routine */ minval = FLT_MAX ; maxval = -FLT_MAX ; for ( i = 0 ; i < nel ; i++ ) { xdat[i] = i ; wdat[i] = 1.0 ; if ( half_buffer -> data[i] < minval ) { minval = half_buffer -> data[i] ; } if ( half_buffer -> data[i] > maxval ) { maxval = half_buffer -> data[i] ; } } fitpar[2] = minval ; fitpar[3] = maxval ; /* search for both positions of the half intensity of the hat within the buffer */ init1 = -1 ; for ( i = 0 ; i < nel ; i++ ) { if ( half_buffer -> data[i] >= ( maxval + minval ) / 2. ) { init1 = i ; break ; } } /* determine the initial positions from the found values */ if ( init1 != -1 ) { fitpar[0] = ((float)init1 - 1.) ; fitpar[1] = ((float)init1 + 1.) ; } for ( i = 0 ; i < NPAR ; i++ ) { mpar[i] = 1 ; dervpar[i] = 0. ; } xdim = XDIMA ; ndat = nel ; numpar = NPAR ; tol = TOLA ; lab = LABA ; its = ITSA ; /* finally, do the least squares fit over the buffer data */ if ( 0 > ( iters = sinfo_new_lsqfit_edge( xdat, &xdim, half_buffer -> data, wdat, &ndat, fitpar, dervpar, mpar, &numpar, &tol, &its, &lab )) ) { /* if the fit doesn't succeed the initial values are taken */ sinfo_msg_warning (" least squares fit failed, error" " no.: %d in slitlet: %d\n", iters, j) ; fitpar[0] = ((float)init1 - 1.) ; fitpar[1] = ((float)init1 + 1.) ; } pos = (fitpar[0] + fitpar[1]) / 2. ; /*---------------------------------------------------------------- * now discern the left and the right sinfo_edge fit of the * slitlets and associate the fit results with the absolute * positions in the image consider the difference of the fitted * slit position to the expected position and decide weather the * fit is taken or the expected value is taken. */ if ( left_right == 0 ) { /* take care of the column position of the fit boxes to get the absolute positions */ if ( j == 0 ) { if ( fabs(pos - ((float)edgeclean[0] - 1. - (float)slit_length)) < diff_tol ) { sinfo_slit_pos[0][0] = pos ; } else { sinfo_msg_warning("something wrong with fitted left" " position of slitlet 0") ; if ( (float) edgeclean[0] - 1. - (float)slit_length < 0. ) { sinfo_slit_pos[0][0] = 0. ; } else { sinfo_slit_pos[0][0] = (float)edgeclean[0] - 1. - (float)slit_length ; } } } else if ( j < k ) { if ( fabs( pos - (float)margin ) < diff_tol ) { sinfo_slit_pos[j][0] = pos + (float)edgeclean[j-1] - (float)margin ; } else { sinfo_msg_warning("something wrong with fitted left" " position of slitlet %d", j) ; sinfo_slit_pos[j][0] = (float)edgeclean[j-1] - 1. ; } } else { if ( fabs( pos - (float)margin ) < diff_tol ) { sinfo_slit_pos[k][0] = pos + (float)edgeclean[k-1] - (float)margin ; } else { sinfo_msg_warning("something wrong with fitted left" " position of slitlet %d", j) ; sinfo_slit_pos[k][0] = (float)edgeclean[k-1] - 1. ; } } } else { /* take care of the column position of the fit boxes to get the absolute positions */ if ( j == 0 ) { if ( fabs( (float)box_buffer->n_elements - pos - (float)edgeclean[0] ) < diff_tol ) { sinfo_slit_pos[0][1] = (float)(box_buffer->n_elements - 1) - pos ; } else { sinfo_msg_warning("something wrong with fitted " "right position of slitlet 0") ; sinfo_slit_pos[0][1] = (float)edgeclean[0] - 1. ; } } else if ( j < k ) { if ( fabs( (float)box_buffer->n_elements - pos + (float)edgeclean[j-1] - (float)margin - (float)edgeclean[j] ) < diff_tol ) { sinfo_slit_pos[j][1] = (float)(box_buffer->n_elements - 1) - pos + (float)edgeclean[j-1] - (float)margin ; } else { sinfo_msg_warning("something wrong with fitted right " "position of slitlet %d", j) ; sinfo_slit_pos[j][1] = (float)edgeclean[j] - 1. ; } } else { if ( edgeclean[k-1] + slit_length > ilx ) { last_pos = (float)(ilx - 1) ; } else { last_pos = (float)(edgeclean[k-1] - 1 + slit_length) ; } if ( fabs( (float)(box_buffer->n_elements - 1) - pos + (float)edgeclean[k-1] - (float)margin - last_pos ) < diff_tol ) { sinfo_slit_pos[k][1] = (float)(box_buffer->n_elements - 1) - pos + (float)edgeclean[k-1] - (float)margin ; } else { sinfo_msg_warning("something wrong with fitted right " "position of slitlet %d\n", j) ; sinfo_slit_pos[k][1] = last_pos ; } } } sinfo_new_destroy_vector ( half_buffer ) ; cpl_free( xdat ) ; cpl_free( wdat ) ; cpl_free( mpar ) ; } sinfo_new_destroy_vector ( box_buffer ) ; } cpl_free( sinfo_edge ) ; cpl_free( pos_row ) ; cpl_free( edgeclean ) ; cpl_free( dummyedge ) ; cpl_free( pos_rowclean ) ; cpl_free(position); return 0 ; } /** @name sinfo_new_fit_slits_edge() @param lineImage emission line frame @param par fit parameter data structure of fitted lines @param sinfo_slit_pos allocated dummy array for the slitlet positions [32][4] @param box_length pixel length of the row box within the fit is done @param y_box small box in spectral direction within the slitlet may lie. @param diff_tol maximum tolerable difference of the resulting fit position with respect to the expected position. If difference is greater the expected position is taken. @return sinfo_slit_pos beginning and end position of the slitlets to sub-pixel accuracy # 0 if it worked, # -1 if there was no line image given, # -2 if there were no line fit parameters given, # -3 if there was no dummy array for the slit positions allocated # -4 if the given box length is impossible # -5 if the given y box length is impossible # -6 if the given difference tolerance is too small # -7 if there were no emission lines found in the first image columns # -8 if not all slitlets could be found @doc fits the beginning and end position of the slitlets by using non-linear least square fitting of a hat function fits a step function to the slitlet edges exposed and indicated by the brightest emission lines. To achieve this, the fit parameters are used to find the brightest emission line and to get its position for each column. The least squares fit is done by using a box smaller than the size of two slitlets */ int sinfo_new_fit_slits_edge( cpl_image * lineImage, FitParams ** par, float ** sinfo_slit_pos, int box_length, float y_box, float diff_tol ) { float* position=NULL ; int * sinfo_edge, * edgeclean ; int * dummyedge ; int * pos_row, * pos_rowclean ; Vector * box_buffer ; Vector * half_buffer ; int row, col ; int i, j, k, n, ed ; int init1 ; int line ; int nel, n_right, left_right ; int slit_length ; int agreed ; int bad_line ; int margin ; int iters, xdim, ndat ; int numpar, its ; int * mpar ; float * xdat, * wdat ; float tol, lab ; float fitpar[NPAR] ; float dervpar[NPAR] ; float minval, maxval ; float pos, last_pos ; int ilx=0; //int ily=0; float* pidata=NULL; if ( NULL == lineImage ) { sinfo_msg_error(" no line image given!" ) ; return -1 ; } ilx=cpl_image_get_size_x(lineImage); //ily=cpl_image_get_size_y(lineImage); pidata=cpl_image_get_data_float(lineImage); slit_length = (int) sqrt (ilx) ; if ( NULL == par ) { sinfo_msg_error(" no line fit parameters given!" ) ; return -2 ; } if ( NULL == sinfo_slit_pos ) { sinfo_msg_error(" no position array allocated!" ) ; return -3 ; } if ( box_length < 4 || box_length >= 2*slit_length ) { sinfo_msg_error(" wrong fitting box length given!" ) ; sinfo_msg_error(" Must be 4 <= box_length < %d ",2*slit_length ) ; sinfo_msg_error(" You have chosen box_length = %d ",box_length); return -4 ; } if ( y_box <= 0. || y_box > 3. ) { sinfo_msg_error(" wrong y box length given!" ) ; sinfo_msg_error(" y_box=%f not in range (0,3]!",y_box); return -5 ; } if ( diff_tol < 1. ) { sinfo_msg_error(" diff_tol too small!" ) ; return -6 ; } /* allocate memory for the edges and the row position of the slitlets */ sinfo_edge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ; dummyedge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ; edgeclean = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ; pos_row = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ; pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ; /* ------------------------------------------------------------------------ * go through the first image columns and the fit parameters and find * the line with the highest intensity */ agreed = -1 ; bad_line = -1 ; while( agreed == -1 ) { int found = -1 ; float max_intensity = -FLT_MAX ; for ( col = 0 ; col < slit_length ; col++ ) { for ( i = 0 ; i < par[0]->n_params ; i++ ) { if ( par[i]->column == col && par[i] -> line != bad_line ) { if ( par[i]->fit_par[0] > max_intensity ) { if ( par[i]->fit_par[1] >= 1. && par[i]->fit_par[2] > 0. ) { max_intensity = par[i]->fit_par[0] ; found = i ; } } } } } /* -------------------------------------------------------------------- * check if the found line is usable and if the neighbouring line * have intensity on near rows in neighbouring slitlets */ line = par[found]->line ; int column = par[found]->column ; float row_pos = par[found]->fit_par[2] ; if ( found >= 0 && max_intensity > 0. ) { for ( i = 0 ; i < par[0]->n_params ; i++ ) { if ( par[i]->line == line-1 && par[i]->column == column + slit_length ) { if ( par[i]->fit_par[2] <= (row_pos + y_box) && par[i]->fit_par[2] >= (row_pos - y_box) ) { bad_line = line ; } } } if ( bad_line != line ) { agreed = 1 ; break ; } } else { sinfo_msg_error("no emission line found in " "the first image columns") ; cpl_free( sinfo_edge ) ; cpl_free( pos_row ) ; cpl_free( edgeclean ) ; cpl_free( dummyedge ) ; cpl_free( pos_rowclean ) ; return -7 ; } } if ( agreed == -1 ) { sinfo_msg_error(" no emission line found in the first image columns") ; cpl_free( sinfo_edge ) ; cpl_free( pos_row ) ; cpl_free( edgeclean ) ; cpl_free( dummyedge ) ; cpl_free( pos_rowclean ) ; return -7 ; } /* now find and store the raw sinfo_edge positions of the found slitlet */ n = 0 ; ed = 0 ; position=cpl_calloc(ilx,sizeof(float)) ; for ( col = 0 ; col < ilx ; col++ ) { for ( i = 0 ; i < par[0]->n_params ; i++ ) { if ( par[i]->column == col && par[i]->line == line ) { if ( par[i]->fit_par[0] > 0. && par[i]->fit_par[1] >= 1. && par[i]->fit_par[2] > 0. ) { position[n] = par[i]->fit_par[2] ; if ( n > 0 && fabs(position[n] - position[n-1]) > y_box ) { sinfo_edge[ed] = col ; pos_row[ed] = sinfo_new_nint( position[n-1] ) ; ed++ ; if ( col >= ilx - slit_length - 5 ) { pos_row[ed] = sinfo_new_nint( position[n] ) ; } } n++ ; } } } } if ( ed < (slit_length - 1) ) { sinfo_msg_error(" not enough slitlets found") ; cpl_free( sinfo_edge ) ; cpl_free( pos_row ) ; cpl_free( edgeclean ) ; cpl_free( dummyedge ) ; cpl_free( pos_rowclean ) ; cpl_free( position ); return -8 ; } /* now find the clean sinfo_edge and row positions of the slitlets */ for ( i = 1 ; i <= ed ; i ++ ) { if ( i == ed ) { if ( (sinfo_edge[i-1] - sinfo_edge[i-2]) < slit_length - 3 || (sinfo_edge[i-1] - sinfo_edge[i-2]) > slit_length + 3 ) { dummyedge[i-1] = -1 ; } } if (dummyedge[i-1] != -1) { dummyedge[i-1] = sinfo_edge[i-1] ; } else { continue ; } if ( i < ed ) { if ( (sinfo_edge[i] - sinfo_edge[i-1]) < slit_length - 3 || (sinfo_edge[i] - sinfo_edge[i-1]) > slit_length + 3 ) { dummyedge[i] = -1 ; } } if ( i + 1 < ed && dummyedge[i] != -1 ) { if ( (sinfo_edge[i+1] - sinfo_edge[i]) < slit_length - 3 || (sinfo_edge[i+1] - sinfo_edge[i]) > slit_length + 3 ) { dummyedge[i+1] = -1 ; } } } k = 0 ; for ( i = 0 ; i < ed ; i++ ) { if ( dummyedge[i] != -1 && dummyedge[i] != 0 ) { edgeclean[k] = dummyedge[i] ; pos_rowclean[k] = pos_row[i] ; k++ ; if( edgeclean[k-1] > (ilx - slit_length - 6 ) ) { pos_rowclean[k] = pos_row[ed] ; } } } if ( k != slit_length - 1 ) { sinfo_msg_error(" not enough clean slitlets found") ; cpl_free( sinfo_edge ) ; cpl_free( pos_row ) ; cpl_free( edgeclean ) ; cpl_free( dummyedge ) ; cpl_free( pos_rowclean ) ; return -8 ; } /* determine the margins of the fitting box outside the slitlets */ margin = box_length / 2 ; /* ------------------------------------------------------------------------ * now go through the slitlets, search along each column within a box with * half width y_box the maximum value and store these found values in a * buffer */ for ( j = 0 ; j <= k ; j++ ) { int m = 0 ; if ( j == 0 ) { box_buffer = sinfo_new_vector( edgeclean[0] + margin ) ; for( col = 0 ; col < edgeclean[0] + margin ; col++ ) { maxval = -FLT_MAX ; for ( row = pos_rowclean[0] - sinfo_new_nint(y_box) ; row <= pos_rowclean[0] + sinfo_new_nint(y_box) ; row++ ) { if ( maxval < pidata[col + ilx*row] ) { maxval = pidata[col + ilx*row] ; } } box_buffer->data[m] = maxval ; m++ ; } } else if ( j < k ) { box_buffer = sinfo_new_vector( edgeclean[j] - edgeclean[j-1] + 2*margin ) ; for ( col = edgeclean[j - 1] - margin ; col < edgeclean[j] + margin ; col++ ) { maxval = -FLT_MAX ; for ( row = pos_rowclean[j] - sinfo_new_nint(y_box) ; row <= pos_rowclean[j] + sinfo_new_nint(y_box) ; row++ ) { if ( maxval < pidata[col + ilx*row] ) { maxval = pidata[col + ilx*row] ; } } box_buffer->data[m] = maxval ; m++ ; } } else { box_buffer = sinfo_new_vector( ilx - edgeclean[k-1] + margin ) ; for ( col = edgeclean[k - 1] - margin ; col < ilx ; col++ ) { maxval = -FLT_MAX ; for ( row = pos_rowclean[k] - sinfo_new_nint(y_box) ; row <= pos_rowclean[k] + sinfo_new_nint(y_box) ; row++ ) { if ( maxval < pidata[col + ilx*row] ) { maxval = pidata[col + ilx*row] ; } } box_buffer->data[m] = maxval ; m++ ; } } for ( left_right = 0 ; left_right <= 1 ; left_right++ ) { nel = 0 ; if ( left_right == 0 ) { nel = box_buffer -> n_elements / 2 ; } else { if ( box_buffer -> n_elements % 2 == 0 ) { nel = box_buffer -> n_elements / 2 ; } else { nel = box_buffer -> n_elements / 2 + 1 ; } } /* now split the buffer in the midth in a left and right part for fitting */ half_buffer = sinfo_new_vector( nel ) ; if ( left_right == 0 ) { for ( i = 0 ; i < nel ; i++ ) { half_buffer -> data[i] = box_buffer -> data[i] ; } } else { n_right = 0 ; for ( i = box_buffer -> n_elements - 1 ; i >= box_buffer -> n_elements - nel ; i-- ) { half_buffer -> data[n_right] = box_buffer -> data[i] ; n_right++ ; } } xdat = (float *) cpl_calloc( nel, sizeof (float) ) ; wdat = (float *) cpl_calloc( nel, sizeof (float) ) ; mpar = (int *) cpl_calloc( NPAR, sizeof (int) ) ; /* set initial values for the fitting routine */ minval = FLT_MAX ; maxval = -FLT_MAX ; for ( i = 0 ; i < nel ; i++ ) { xdat[i] = i ; wdat[i] = 1.0 ; if ( half_buffer -> data[i] < minval ) { minval = half_buffer -> data[i] ; } if ( half_buffer -> data[i] > maxval ) { maxval = half_buffer -> data[i] ; } } fitpar[2] = minval ; fitpar[3] = maxval ; /* search for both positions of the half intensity of the hat within the buffer */ init1 = -1 ; for ( i = 0 ; i < nel ; i++ ) { if ( half_buffer -> data[i] >= ( maxval + minval ) / 2. ) { init1 = i ; break ; } } /* determine the initial positions from the found values */ if ( init1 != -1 ) { fitpar[0] = ((float)init1 - 1.) ; fitpar[1] = ((float)init1 + 1.) ; } for ( i = 0 ; i < NPAR ; i++ ) { mpar[i] = 1 ; dervpar[i] = 0. ; } xdim = XDIMA ; ndat = nel ; numpar = NPAR ; tol = TOLA ; lab = LABA ; its = ITSA ; /* finally, do the least squares fit over the buffer data */ if ( 0 > ( iters = sinfo_new_lsqfit_edge( xdat, &xdim, half_buffer -> data, wdat, &ndat, fitpar, dervpar, mpar, &numpar, &tol, &its, &lab )) ) { /* if the fit doesn't succeed the initial values are taken */ sinfo_msg_warning ("least squares fit failed, error " "no.: %d in slitlet: %d", iters, j) ; fitpar[0] = ((float)init1 - 1.) ; fitpar[1] = ((float)init1 + 1.) ; } pos = (fitpar[0] + fitpar[1]) / 2. ; /*----------------------------------------------------------------- * now discern the left and the right sinfo_edge fit of the * slitlets and associate the fit results with the absolute * positions in the image consider the difference of the fitted * slit position to the expected position and decide weather the * fit is taken or the expected value is taken. */ if ( left_right == 0 ) { /* take care of the column position of the fit boxes to get the absolute positions */ if ( j == 0 ) { if ( fabs(pos - ((float)edgeclean[0] - 1. - (float)slit_length)) < diff_tol ) { sinfo_slit_pos[0][0] = pos ; } else { sinfo_msg_warning("something wrong with fitted " "left position of slitlet 0") ; if ((float) edgeclean[0] - 1. - (float)slit_length < 0. ) { sinfo_slit_pos[0][0] = 0. ; } else { sinfo_slit_pos[0][0] = (float)edgeclean[0] - 1. - (float)slit_length ; } } } else if ( j < k ) { if ( fabs( pos - (float)margin ) < diff_tol ) { sinfo_slit_pos[j][0] = pos + (float)edgeclean[j-1] - (float)margin ; } else { sinfo_msg_warning("something wrong with fitted " "left position of slitlet %d", j) ; sinfo_slit_pos[j][0] = (float)edgeclean[j-1] - 1. ; } } else { if ( fabs( pos - (float)margin ) < diff_tol ) { sinfo_slit_pos[k][0] = pos + (float)edgeclean[k-1] - (float)margin ; } else { sinfo_msg_warning("something wrong with fitted left " "position of slitlet %d", j) ; sinfo_slit_pos[k][0] = (float)edgeclean[k-1] - 1. ; } } } else { /* take care of the column position of the fit boxes to get the absolute positions */ if ( j == 0 ) { if ( fabs( (float)box_buffer->n_elements - pos - (float)edgeclean[0] ) < diff_tol ) { sinfo_slit_pos[0][1] = (float)(box_buffer->n_elements - 1) - pos ; } else { sinfo_msg_warning("something wrong with fitted " "right position of slitlet 0") ; sinfo_slit_pos[0][1] = (float)edgeclean[0] - 1. ; } } else if ( j < k ) { if ( fabs( (float)box_buffer->n_elements - pos + (float)edgeclean[j-1] - (float)margin - (float)edgeclean[j] ) < diff_tol ) { sinfo_slit_pos[j][1] = (float)(box_buffer->n_elements - 1) - pos + (float)edgeclean[j-1] - (float)margin; } else { sinfo_msg_warning("something wrong with fitted " "right position of slitlet %d", j) ; sinfo_slit_pos[j][1] = (float)edgeclean[j] - 1. ; } } else { if ( edgeclean[k-1] + slit_length > ilx ) { last_pos = (float)(ilx - 1) ; } else { last_pos = (float)(edgeclean[k-1] - 1 + slit_length) ; } if ( fabs( (float)(box_buffer->n_elements - 1) - pos + (float)edgeclean[k-1] - (float)margin - last_pos ) < diff_tol ) { sinfo_slit_pos[k][1] = (float)(box_buffer->n_elements - 1) - pos + (float)edgeclean[k-1] - (float)margin ; } else { sinfo_msg_warning("something wrong with fitted " "right position of slitlet %d", j) ; sinfo_slit_pos[k][1] = last_pos ; } } } sinfo_new_destroy_vector ( half_buffer ) ; cpl_free( xdat ) ; cpl_free( wdat ) ; cpl_free( mpar ) ; } sinfo_new_destroy_vector ( box_buffer ) ; } cpl_free( sinfo_edge ) ; cpl_free( pos_row ) ; cpl_free( edgeclean ) ; cpl_free( dummyedge ) ; cpl_free( pos_rowclean ) ; cpl_free( position ); return 0 ; } /** @name sinfo_new_fit_slits_boltz_with_estimate() @param lineImage emission line frame @param sinfo_slit_pos estimation array for the slitlet positions [min32][2] @param box_length pixel length of the row box within the fit is done @param y_box small box in spectral direction within the slitlet may lie. @param low_pos @param high_pos pixel positions in spectral direction between which the line should be located. @return sinfo_slit_pos beginning and end position of the slitlets to sub-pixel accuracy 0 if it worked, -1 if it failed, @doc fits the beginning and end position of the slitlets by using non-linear least square fitting of a Boltzmann function fits a Boltzmann function to the slitlet edges exposed and indicated by the brightest emission lines. The slitlet is searched within user given positions. The least squares fit is done by using a box smaller than the size of two slitlets TODO: * 8) fits the beginning and end position of the slitlets * by using non-linear least square fitting of an edge function * fits a linear edge function to the slitlet edges exposed and indicated * by the brightest emission lines. The slitlet is searched within * user given positions. * The least squares fit is done by using a box smaller than * the size of two slitlets */ int sinfo_new_fit_slits_edge_with_estimate ( cpl_image * lineImage, float ** sinfo_slit_pos, int box_length, float y_box, float diff_tol, int low_pos, int high_pos ) { int* position=NULL ; Vector * box_buffer ; Vector * in_buffer ; int row, col ; int col_first, col_last ; int row_first, row_last ; int i, j, m; int init1 ; int left_right ; int n_buf, shift ; int slit_length ; int iters, xdim, ndat ; int numpar, its ; int * mpar ; float * xdat, * wdat ; float tol, lab ; float fitpar[NPAR] ; float dervpar[NPAR] ; float pos ; float new_pos ; int slitposition[SLITLENGTH] ; pixelvalue rowpos[SLITLENGTH] ; int ilx=0; int ily=0; float* pidata=NULL; /* slit_length = SLITLENGTH ; this is too high: 64 */ slit_length = N_SLITLETS ; /* this is better: 32 */ if ( NULL == lineImage ) { sinfo_msg_error(" no line image given!" ) ; return -1 ; } ilx=cpl_image_get_size_x(lineImage); ily=cpl_image_get_size_y(lineImage); pidata=cpl_image_get_data_float(lineImage); if ( NULL == sinfo_slit_pos ) { sinfo_msg_error(" no position array allocated!" ) ; return -1 ; } if ( box_length < 4 || box_length > 2*slit_length ) { sinfo_msg_error(" wrong fitting box length given!" ) ; sinfo_msg_error(" Must be 4 <= box_length < %d ",2*slit_length ) ; sinfo_msg_error(" You have chosen box_length = %d",box_length); return -1 ; } if ( y_box <= 0. || y_box > 6. ) { sinfo_msg_error("wrong y box length given!" ) ; sinfo_msg_error("You have chosen y_box=%f not in range (0,6]!",y_box); return -1 ; } if ( diff_tol <= 0. ) { sinfo_msg_error(" wrong diff_tol given!" ) ; return -1 ; } if ( low_pos >= high_pos || low_pos < 0 || high_pos <= 0 || high_pos > ily ) { sinfo_msg_error(" wrong user given search positions!" ) ; return -1 ; } /* now search for the maximum between the given positions for each col */ position=cpl_calloc(ilx,sizeof(int)) ; float maxval,minval; for ( col = 0 ; col < ilx ; col++ ) { int found_row = -1 ; maxval = -FLT_MAX ; for ( row = low_pos ; row <= high_pos ; row++ ) { if ( maxval < pidata[col+row*ilx] ) { maxval = pidata[col+row*ilx] ; found_row = row ; } } if ( maxval > -FLT_MAX && found_row > low_pos ) { position[col] = found_row ; } else { position[col] = 0 ; } } /* ------------------------------------------------------------------------ * now go through the slitlets, search along each column within a box with * half width y_box the maximum value and store these found values in a * buffer */ for ( j = 0 ; j < slit_length ; j++ ) { /* now go through the columns and determine the slitlet positions by * calculating the sinfo_median of the found positions */ int n = 0 ; for ( col = sinfo_new_nint(sinfo_slit_pos[j][0])+ 1 ; col < sinfo_new_nint(sinfo_slit_pos[j][1]) -1 ; col++ ) { rowpos[n] = (pixelvalue)position[col] ; n++ ; } slitposition[j] = (int)sinfo_new_median(rowpos, n) ; for ( left_right = 0 ; left_right <= 1 ; left_right++ ) { init1 = 0 ; col_first = sinfo_new_nint( sinfo_slit_pos[j][left_right] ) - box_length/2 ; col_last = sinfo_new_nint( sinfo_slit_pos[j][left_right] ) + box_length/2 ; if ( col_first < 0 ) { col_first = 0 ; init1 = 1 ; } if ( col_last > ilx ) { col_last = ilx ; init1 = 1 ; } if ( col_last - col_first <= 0 ) { sinfo_msg_error(" first and last column wrong!" ) ; return -1 ; } box_buffer = sinfo_new_vector( col_last - col_first ) ; m = 0 ; if ( left_right == 0 ) { for( col = col_first ; col < col_last ; col++ ) { row_first = slitposition[j] - sinfo_new_nint(y_box) ; row_last = slitposition[j] + sinfo_new_nint(y_box) ; if ( row_first < 0 ) { row_first = 0 ; } if ( row_last >= ily ) { row_last = ily - 1 ; } maxval = -FLT_MAX ; for ( row = row_first ; row <= row_last ; row++ ) { if ( maxval < pidata[col + ilx*row] ) { maxval = pidata[col + ilx*row] ; } } box_buffer->data[m] = maxval ; m++ ; } } else { for( col = col_last-1 ; col >= col_first ; col-- ) { row_first = slitposition[j] - sinfo_new_nint(y_box) ; row_last = slitposition[j] + sinfo_new_nint(y_box) ; if ( row_first < 0 ) { row_first = 0 ; } if ( row_last >= ily ) { row_last = ily - 1 ; } maxval = -FLT_MAX ; for ( row = row_first ; row <= row_last ; row++ ) { if ( maxval < pidata[col + ilx*row] ) { maxval = pidata[col + ilx*row] ; } } box_buffer->data[m] = maxval ; m++ ; } } xdat=(float *)cpl_calloc( box_buffer->n_elements, sizeof (float)); wdat=(float *)cpl_calloc( box_buffer->n_elements, sizeof (float)); mpar=(int *) cpl_calloc( NPAR, sizeof (int) ) ; /* set initial values for the fitting routine */ minval = FLT_MAX ; maxval = -FLT_MAX ; for ( i = 0 ; i < box_buffer->n_elements ; i++ ) { xdat[i] = i ; wdat[i] = 1.0 ; if ( box_buffer -> data[i] < minval ) { minval = box_buffer -> data[i] ; } if ( box_buffer -> data[i] > maxval ) { maxval = box_buffer -> data[i] ; } } fitpar[2] = minval ; fitpar[3] = maxval ; /*---------------------------------------------------------------- * if we have too few left background values (at the image edges) * the left margin of the buffer to fit is filled with the minimal * values in order to get a good fit */ if ( init1 == 1 ) { n_buf = box_buffer->n_elements + box_length/2 ; in_buffer = sinfo_new_vector( n_buf ) ; for ( i = 0 ; i < box_length/2 ; i++ ) { in_buffer -> data[i] = minval ; } shift = 0 ; for ( i = box_length/2 ; i < n_buf ; i++ ) { in_buffer -> data[i] = box_buffer -> data[shift] ; shift++ ; } sinfo_new_destroy_vector ( box_buffer ) ; box_buffer = sinfo_new_vector ( n_buf ) ; for ( i = 0 ; i < n_buf ; i++ ) { box_buffer -> data[i] = in_buffer -> data[i] ; } sinfo_new_destroy_vector ( in_buffer ) ; } /* determine the initial positions from the found values */ fitpar[0] = (float)box_buffer->n_elements/2. - 1. ; fitpar[1] = (float)box_buffer->n_elements/2. + 1. ; for ( i = 0 ; i < NPAR ; i++ ) { mpar[i] = 1 ; dervpar[i] = 0. ; } xdim = XDIMA ; ndat = box_buffer->n_elements ; numpar = NPAR ; tol = TOLA ; lab = LABA ; its = ITSA ; /* finally, do the least squares fit over the buffer data */ if ( 0 > ( iters = sinfo_new_lsqfit_edge( xdat, &xdim, box_buffer -> data, wdat, &ndat, fitpar, dervpar, mpar, &numpar, &tol, &its, &lab )) ) { sinfo_msg_warning (" least squares fit failed, error " "no.: %d in slitlet: %d\n", iters, j) ; sinfo_new_destroy_vector(box_buffer) ; cpl_free( xdat ) ; cpl_free( wdat ) ; cpl_free( mpar ) ; continue ; } if ( fitpar[1] <= fitpar[0] ) { sinfo_msg_warning ("fit failed due to negative slope of " "sinfo_new_edge function in slitlet: %d",j); sinfo_new_destroy_vector(box_buffer) ; cpl_free( xdat ) ; cpl_free( wdat ) ; cpl_free( mpar ) ; continue ; } pos = (fitpar[0] + fitpar[1])/2. ; if ( init1 == 1 ) { pos -= (float)box_length/2. ; } /*------------------------------------------------------------- * now compute the real slit positions using the guess positions * if the fit did not work the guess positions are taken * the same is done if the deviations are too big. */ if ( pos != 0. ) { if ( left_right == 0 ) { new_pos = (float)col_first + pos ; } else { new_pos = (float)col_last-1 - pos ; } if ( fabs(new_pos - sinfo_slit_pos[j][left_right]) < diff_tol ) { sinfo_slit_pos[j][left_right] = new_pos ; } else { sinfo_msg_warning (" deviation bigger than tolerance," " take the estimated slitlet positiona" " in slitlet: %d\n", j) ; } } cpl_free( xdat ) ; cpl_free( wdat ) ; cpl_free( mpar ) ; sinfo_new_destroy_vector ( box_buffer ) ; } } cpl_free(position); return 0 ; } /**@}*/