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