00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00014
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
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 #include "vltPort.h"
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134 #include "absolute.h"
00135 #include "recipes.h"
00136
00137
00138
00139
00140 static float sqrarg ;
00141 #define SQR(a) (sqrarg = (a) , sqrarg*sqrarg)
00142
00143 #define XDIMA 1
00144 #define TOLA 0.001
00145 #define LABA 0.1
00146 #define ITSA 200
00147 #define LABFACA 10.0
00148 #define LABMAXA 1.0e+10
00149 #define LABMINA 1.0e-10
00150 #define NPAR 4
00151
00152
00153
00154
00155
00156 static double chi1 ;
00157 static double chi2 ;
00158 static double labda ;
00159 static double vec[NPAR] ;
00160 static double matrix1[NPAR][NPAR] ;
00161 static double matrix2[NPAR][NPAR] ;
00162 static int nfree ;
00163 static int parptr[NPAR] ;
00164 static float slopewidth ;
00165
00166
00167
00168
00169 static int inv_mat_edge (void) ;
00170
00171 static void get_mat_edge ( float * xdat,
00172 int * xdim,
00173 float * ydat,
00174 float * wdat,
00175 int * ndat,
00176 float * fpar,
00177 float * epar
00178 ) ;
00179
00180 static int get_vec_edge ( float * xdat,
00181 int * xdim,
00182 float * ydat,
00183 float * wdat,
00184 int * ndat,
00185 float * fpar,
00186 float * epar,
00187 int * npar ) ;
00188
00189
00190
00191
00192
00193
00215
00216 float edge ( float * xdat, float * parlist )
00217 {
00218 float return_value ;
00219 float slope1 ;
00220
00221
00222 slope1 = ( parlist[3] - parlist[2] ) / ( parlist[1] - parlist[0] ) ;
00223
00224
00225 if ( xdat[0] <= parlist[0] )
00226 {
00227 return_value = parlist[2] ;
00228 }
00229 else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[1] )
00230 {
00231 return_value = (xdat[0] - parlist[0]) * slope1 + parlist[2] ;
00232 }
00233 else if ( xdat[0] > parlist[1] )
00234 {
00235 return_value = parlist[3] ;
00236 }
00237 else
00238 {
00239 return_value = 0. ;
00240 }
00241
00242 return return_value ;
00243 }
00244
00245
00272
00273
00274
00275 float hat2 ( float * xdat, float * parlist )
00276 {
00277 float return_value ;
00278 float slope1, slope2, slope3 ;
00279
00280
00281 slope1 = ( parlist[6] - parlist[4] ) / ( parlist[1] - parlist[0] ) ;
00282 slope2 = ( parlist[7] - parlist[5] ) / ( parlist[3] - parlist[2] ) ;
00283 slope3 = ( parlist[7] - parlist[6] ) / ( parlist[2] - parlist[1] ) ;
00284
00285
00286 if ( xdat[0] <= parlist[0] )
00287 {
00288 return_value = parlist[4] ;
00289 }
00290 else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[1] )
00291 {
00292 return_value = (xdat[0] - parlist[0]) * slope1 + parlist[4] ;
00293 }
00294 else if ( xdat[0] > parlist[1] && xdat[0] <= parlist[2] )
00295 {
00296 return_value = (xdat[0] - parlist[1]) * slope3 + parlist[6] ;
00297 }
00298 else if ( xdat[0] > parlist[2] && xdat[0] <= parlist[3] )
00299 {
00300 return_value = (parlist[3] - xdat[0]) * slope2 + parlist[5] ;
00301 }
00302 else if ( xdat[0] > parlist[3] )
00303 {
00304 return_value = parlist[5] ;
00305 }
00306 else
00307 {
00308 return_value = 0. ;
00309 }
00310
00311 return return_value ;
00312 }
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336 float hat1 ( float * xdat, float * parlist )
00337 {
00338 float return_value=0 ;
00339 float slope1, slope2 ;
00340
00341
00342
00343
00344 slope1 = (parlist[4] - parlist[2]) / slopewidth ;
00345 slope2 = (parlist[4] - parlist[3]) / slopewidth ;
00346
00347
00348 if ( xdat[0] <= parlist[0] )
00349 {
00350 return_value = parlist[2] ;
00351 }
00352 else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[0] + slopewidth )
00353 {
00354 return_value = (xdat[0] - parlist[0]) * slope1 + parlist[2] ;
00355 }
00356 else if ( xdat[0] > parlist[0] + slopewidth && xdat[0] <= parlist[1] - slopewidth )
00357 {
00358 return_value = parlist[4] ;
00359 }
00360 else if ( xdat[0] > parlist[1] - slopewidth && xdat[0] <= parlist[1] )
00361 {
00362 return_value = (parlist[1] - xdat[0]) * slope2 + parlist[3] ;
00363 }
00364 else if ( xdat[0] > parlist[1] )
00365 {
00366 return_value = parlist[3] ;
00367 }
00368
00369 return return_value ;
00370 }
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391 void edge_deriv( float * xdat, float * parlist, float * dervs )
00392 {
00393 float deriv1_slope1 ;
00394
00395
00396 deriv1_slope1 = ( parlist[3] - parlist[2] ) / SQR(parlist[1] - parlist[0]) ;
00397
00398
00399 if ( xdat[0] <= parlist[0] )
00400 {
00401 dervs[0] = 0. ;
00402 dervs[1] = 0. ;
00403 dervs[2] = 1. ;
00404 dervs[3] = 0. ;
00405 }
00406 else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[1] )
00407 {
00408 dervs[0] = ( xdat[0] - parlist[1] ) * deriv1_slope1 ;
00409 dervs[1] = ( parlist[0] - xdat[0] ) * deriv1_slope1 ;
00410 dervs[2] = ( parlist[0] - xdat[0] ) / ( parlist[1] - parlist[0] ) + 1. ;
00411 dervs[3] = ( xdat[0] - parlist[0] ) / ( parlist[1] - parlist[0] ) ;
00412 }
00413 else if ( xdat[0] > parlist[1] )
00414 {
00415 dervs[0] = 0. ;
00416 dervs[1] = 0. ;
00417 dervs[2] = 0. ;
00418 dervs[3] = 1. ;
00419 }
00420 }
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449 void hat_deriv2( float * xdat, float * parlist, float * dervs )
00450 {
00451 float deriv1_slope1 ;
00452 float deriv1_slope2 ;
00453 float deriv1_slope3 ;
00454
00455
00456 deriv1_slope1 = ( parlist[6] - parlist[4] ) / SQR(parlist[1] - parlist[0]) ;
00457 deriv1_slope2 = ( parlist[7] - parlist[5] ) / SQR(parlist[3] - parlist[2]) ;
00458 deriv1_slope3 = ( parlist[7] - parlist[6] ) / SQR(parlist[2] - parlist[1]) ;
00459
00460
00461 if ( xdat[0] <= parlist[0] )
00462 {
00463 dervs[0] = 0. ;
00464 dervs[1] = 0. ;
00465 dervs[2] = 0. ;
00466 dervs[3] = 0. ;
00467 dervs[4] = 1. ;
00468 dervs[5] = 0. ;
00469 dervs[6] = 0. ;
00470 dervs[7] = 0. ;
00471 }
00472 else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[1] )
00473 {
00474 dervs[0] = ( xdat[0] - parlist[1] ) * deriv1_slope1 ;
00475 dervs[1] = ( parlist[0] - xdat[0] ) * deriv1_slope1 ;
00476 dervs[2] = 0. ;
00477 dervs[3] = 0. ;
00478 dervs[4] = ( parlist[0] - xdat[0] ) / ( parlist[1] - parlist[0] ) + 1. ;
00479 dervs[5] = 0. ;
00480 dervs[6] = ( xdat[0] - parlist[0] ) / ( parlist[1] - parlist[0] ) ;
00481 dervs[7] = 0. ;
00482 }
00483 else if ( xdat[0] > parlist[1] && xdat[0] <= parlist[2] )
00484 {
00485 dervs[0] = 0. ;
00486 dervs[1] = (xdat[0] - parlist[2]) * deriv1_slope3 ;
00487 dervs[2] = (parlist[1] - xdat[0]) * deriv1_slope3 ;
00488 dervs[3] = 0. ;
00489 dervs[4] = 0. ;
00490 dervs[5] = 0. ;
00491 dervs[6] = (parlist[1] - xdat[0]) / (parlist[2] - parlist[1]) + 1. ;
00492 dervs[7] = (xdat[0] - parlist[1]) / (parlist[2] - parlist[1]) ;
00493 }
00494 else if ( xdat[0] > parlist[2] && xdat[0] <= parlist[3] )
00495 {
00496 dervs[0] = 0. ;
00497 dervs[1] = 0. ;
00498 dervs[2] = ( parlist[3] - xdat[0] ) * deriv1_slope2 ;
00499 dervs[3] = ( xdat[0] - parlist[2] ) * deriv1_slope2 ;
00500 dervs[4] = 0. ;
00501 dervs[5] = ( xdat[0] - parlist[3] ) / ( parlist[3] - parlist[2] ) + 1. ;
00502 dervs[6] = 0. ;
00503 dervs[7] = ( parlist[3] - xdat[0] ) / ( parlist[3] - parlist[2] ) ;
00504 }
00505 else if ( xdat[0] > parlist[3] )
00506 {
00507 dervs[0] = 0. ;
00508 dervs[1] = 0. ;
00509 dervs[2] = 0. ;
00510 dervs[3] = 0. ;
00511 dervs[4] = 0. ;
00512 dervs[5] = 1. ;
00513 dervs[6] = 0. ;
00514 dervs[7] = 0. ;
00515 }
00516 }
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540 void hat_deriv1( float * xdat, float * parlist, float * dervs )
00541 {
00542
00543 if ( xdat[0] <= parlist[0] )
00544 {
00545 dervs[0] = 0. ;
00546 dervs[1] = 0. ;
00547 dervs[2] = 1. ;
00548 dervs[3] = 0. ;
00549 dervs[4] = 0. ;
00550 }
00551 else if ( xdat[0] > parlist[0] && xdat[0] <= parlist[0] + slopewidth )
00552 {
00553 dervs[0] = ( parlist[2] - parlist[4] ) / slopewidth ;
00554 dervs[1] = 0. ;
00555 dervs[2] = (( parlist[0] - xdat[0] ) / slopewidth ) + 1. ;
00556 dervs[3] = 0. ;
00557 dervs[4] = ( xdat[0] - parlist[0] ) / slopewidth ;
00558 }
00559 else if ( xdat[0] > parlist[0] + slopewidth && xdat[0] <= parlist[1] - slopewidth )
00560 {
00561 dervs[0] = 0. ;
00562 dervs[1] = 0. ;
00563 dervs[2] = 0. ;
00564 dervs[3] = 0. ;
00565 dervs[4] = 1. ;
00566 }
00567 else if ( xdat[0] > parlist[1] - slopewidth && xdat[0] <= parlist[1] )
00568 {
00569 dervs[0] = 0. ;
00570 dervs[1] = ( parlist[4] - parlist[3] ) / slopewidth ;
00571 dervs[2] = 0. ;
00572 dervs[3] = (( xdat[0] - parlist[1] ) / slopewidth ) + 1. ;
00573 dervs[4] = ( parlist[1] - xdat[0] ) / slopewidth ;
00574 }
00575 else if ( xdat[0] > parlist[1] )
00576 {
00577 dervs[0] = 0. ;
00578 dervs[1] = 0. ;
00579 dervs[2] = 0. ;
00580 dervs[3] = 1. ;
00581 dervs[4] = 0. ;
00582 }
00583 }
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594 static int inv_mat_edge (void)
00595 {
00596 double even ;
00597 double hv[NPAR] ;
00598 double mjk ;
00599 double rowmax ;
00600 int evin ;
00601 int i, j, k, row ;
00602 int per[NPAR] ;
00603
00604
00605 for ( i = 0 ; i < nfree ; i++ )
00606 {
00607 per[i] = i ;
00608 }
00609
00610 for ( j = 0 ; j < nfree ; j++ )
00611 {
00612
00613 rowmax = fabs ( matrix2[j][j] ) ;
00614 row = j ;
00615
00616 for ( i = j + 1 ; i < nfree ; i++ )
00617 {
00618 if ( fabs ( matrix2[i][j] ) > rowmax )
00619 {
00620 rowmax = fabs( matrix2[i][j] ) ;
00621 row = i ;
00622 }
00623 }
00624
00625
00626 if ( matrix2[row][j] == 0.0 )
00627 {
00628 return -6 ;
00629 }
00630
00631
00632 if ( row > j )
00633 {
00634 for ( k = 0 ; k < nfree ; k++ )
00635 {
00636 even = matrix2[j][k] ;
00637 matrix2[j][k] = matrix2[row][k] ;
00638 matrix2[row][k] = even ;
00639 }
00640
00641 evin = per[j] ;
00642 per[j] = per[row] ;
00643 per[row] = evin ;
00644 }
00645
00646
00647 even = 1.0 / matrix2[j][j] ;
00648 for ( i = 0 ; i < nfree ; i++ )
00649 {
00650 matrix2[i][j] *= even ;
00651 }
00652 matrix2[j][j] = even ;
00653
00654 for ( k = 0 ; k < j ; k++ )
00655 {
00656 mjk = matrix2[j][k] ;
00657 for ( i = 0 ; i < j ; i++ )
00658 {
00659 matrix2[i][k] -= matrix2[i][j] * mjk ;
00660 }
00661 for ( i = j + 1 ; i < nfree ; i++ )
00662 {
00663 matrix2[i][k] -= matrix2[i][j] * mjk ;
00664 }
00665 matrix2[j][k] = -even * mjk ;
00666 }
00667
00668 for ( k = j + 1 ; k < nfree ; k++ )
00669 {
00670 mjk = matrix2[j][k] ;
00671 for ( i = 0 ; i < j ; i++ )
00672 {
00673 matrix2[i][k] -= matrix2[i][j] * mjk ;
00674 }
00675 for ( i = j + 1 ; i < nfree ; i++ )
00676 {
00677 matrix2[i][k] -= matrix2[i][j] * mjk ;
00678 }
00679 matrix2[j][k] = -even * mjk ;
00680 }
00681 }
00682
00683
00684 for ( i = 0 ; i < nfree ; i++ )
00685 {
00686 for ( k = 0 ; k < nfree ; k++ )
00687 {
00688 hv[per[k]] = matrix2[i][k] ;
00689 }
00690 for ( k = 0 ; k < nfree ; k++ )
00691 {
00692 matrix2[i][k] = hv[k] ;
00693 }
00694 }
00695
00696
00697 return 0 ;
00698 }
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714 static void get_mat_edge ( float * xdat,
00715 int * xdim,
00716 float * ydat,
00717 float * wdat,
00718 int * ndat,
00719 float * fpar,
00720 float * epar
00721 )
00722 {
00723 double wd ;
00724 double wn ;
00725 double yd ;
00726 int i, j, n ;
00727
00728 for ( j = 0 ; j < nfree ; j++ )
00729 {
00730 vec[j] = 0.0 ;
00731 for ( i = 0 ; i<= j ; i++ )
00732 {
00733 matrix1[j][i] = 0.0 ;
00734 }
00735 }
00736 chi2 = 0.0 ;
00737
00738
00739 for ( n = 0 ; n < (*ndat) ; n++ )
00740 {
00741 wn = wdat[n] ;
00742 if ( wn > 0.0 )
00743 {
00744 yd = ydat[n] - edge( &xdat[(*xdim) * n], fpar ) ;
00745 edge_deriv( &xdat[(*xdim) * n], fpar, epar) ;
00746 chi2 += yd * yd * wn ;
00747 for ( j = 0 ; j < nfree ; j++ )
00748 {
00749 wd = epar[parptr[j]] * wn ;
00750 vec[j] += yd * wd ;
00751 for ( i = 0 ; i <= j ; i++ )
00752 {
00753 matrix1[j][i] += epar[parptr[i]] * wd ;
00754 }
00755 }
00756 }
00757 }
00758 }
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782 static int get_vec_edge ( float * xdat,
00783 int * xdim,
00784 float * ydat,
00785 float * wdat,
00786 int * ndat,
00787 float * fpar,
00788 float * epar,
00789 int * npar )
00790 {
00791 double dj ;
00792 double dy ;
00793 double mii ;
00794 double mji ;
00795 double mjj ;
00796 double wn ;
00797 int i, j, n, r ;
00798
00799
00800 for ( j = 0 ; j < nfree ; j++ )
00801 {
00802 mjj = matrix1[j][j] ;
00803 if ( mjj <= 0.0 )
00804 {
00805 return -5 ;
00806 }
00807 mjj = sqrt( mjj ) ;
00808 for ( i = 0 ; i < j ; i++ )
00809 {
00810 mji = matrix1[j][i] / mjj / sqrt( matrix1[i][i] ) ;
00811 matrix2[i][j] = matrix2[j][i] = mji ;
00812 }
00813 matrix2[j][j] = 1.0 + labda ;
00814 }
00815
00816 if ( (r = inv_mat_edge()) )
00817 {
00818 return r ;
00819 }
00820
00821 for ( i = 0 ; i < (*npar) ; i ++ )
00822 {
00823 epar[i] = fpar[i] ;
00824 }
00825
00826
00827 for ( j = 0 ; j < nfree ; j++ )
00828 {
00829 dj = 0.0 ;
00830 mjj = matrix1[j][j] ;
00831 if ( mjj <= 0.0)
00832 {
00833 return -7 ;
00834 }
00835 mjj = sqrt ( mjj ) ;
00836 for ( i = 0 ; i < nfree ; i++ )
00837 {
00838 mii = matrix1[i][i] ;
00839 if ( mii <= 0.0 )
00840 {
00841 return -7 ;
00842 }
00843 mii = sqrt( mii ) ;
00844 dj += vec[i] * matrix2[j][i] / mjj / mii ;
00845 }
00846 epar[parptr[j]] += dj ;
00847 }
00848 chi1 = 0.0 ;
00849
00850
00851 for ( n = 0 ; n < (*ndat) ; n++ )
00852 {
00853 wn = wdat[n] ;
00854 if ( wn > 0.0 )
00855 {
00856 dy = ydat[n] - edge( &xdat[(*xdim) * n], epar) ;
00857 chi1 += wdat[n] * dy * dy ;
00858 }
00859 }
00860 return 0 ;
00861 }
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911 int lsqfit_edge ( float * xdat,
00912 int * xdim,
00913 float * ydat,
00914 float * wdat,
00915 int * ndat,
00916 float * fpar,
00917 float * epar,
00918 int * mpar,
00919 int * npar,
00920 float * tol ,
00921 int * its ,
00922 float * lab )
00923 {
00924 int i, n, r ;
00925 int itc ;
00926 int found ;
00927 int nuse ;
00928 double tolerance ;
00929
00930 itc = 0 ;
00931 found = 0 ;
00932 nfree = 0 ;
00933 nuse = 0 ;
00934
00935 if ( *tol < (FLT_EPSILON * 10.0 ) )
00936 {
00937 tolerance = FLT_EPSILON * 10.0 ;
00938 }
00939 else
00940 {
00941 tolerance = *tol ;
00942 }
00943
00944 labda = fabs( *lab ) * LABFACA ;
00945 for ( i = 0 ; i < (*npar) ; i++ )
00946 {
00947 if ( mpar[i] )
00948 {
00949 if ( nfree > NPAR )
00950 {
00951 return -1 ;
00952 }
00953 parptr[nfree++] = i ;
00954 }
00955 }
00956
00957 if (nfree == 0)
00958 {
00959 return -2 ;
00960 }
00961
00962 for ( n = 0 ; n < (*ndat) ; n++ )
00963 {
00964 if ( wdat[n] > 0.0 )
00965 {
00966 nuse ++ ;
00967 }
00968 }
00969
00970 if ( nfree >= nuse )
00971 {
00972 return -3 ;
00973 }
00974 if ( labda == 0.0 )
00975 {
00976 for ( i = 0 ; i < nfree ; fpar[parptr[i++]] = 0.0 ) ;
00977 get_mat_edge ( xdat, xdim, ydat, wdat, ndat, fpar, epar) ;
00978 r = get_vec_edge ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00979 if ( r )
00980 {
00981 return r ;
00982 }
00983 for ( i = 0 ; i < (*npar) ; i++ )
00984 {
00985 fpar[i] = epar[i] ;
00986 epar[i] = 0.0 ;
00987 }
00988 chi1 = sqrt( chi1 / (double) (nuse - nfree) ) ;
00989 for ( i = 0 ; i < nfree ; i++ )
00990 {
00991 if ( (matrix1[i][i] <= 0.0 ) || (matrix2[i][i] <= 0.0) )
00992 {
00993 return -7 ;
00994 }
00995 epar[parptr[i]] = chi1 * sqrt( matrix2[i][i] ) / sqrt( matrix1[i][i] ) ;
00996 }
00997 }
00998 else
00999 {
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012 while ( !found )
01013 {
01014 if ( itc++ == (*its) )
01015 {
01016 return -4 ;
01017 }
01018 get_mat_edge( xdat, xdim, ydat, wdat, ndat, fpar, epar ) ;
01019
01020
01021
01022
01023
01024 if ( labda > LABMINA )
01025 {
01026 labda = labda / LABFACA ;
01027 }
01028 r = get_vec_edge ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
01029 if ( (int)fpar[1] - (int)fpar[0] <= 0 && fpar[1] / fpar[0] > 0. )
01030 {
01031 fpar[1] += 1. ;
01032 continue ;
01033 }
01034 if ( r )
01035 {
01036 return r ;
01037 }
01038
01039 while ( chi1 >= chi2 )
01040 {
01041
01042
01043
01044
01045
01046
01047 if ( labda > LABMAXA )
01048 {
01049 break ;
01050 }
01051 labda = labda * LABFACA ;
01052 r = get_vec_edge ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
01053 if ( (int)fpar[1] - (int)fpar[0] <= 0 && fpar[1] / fpar[0] > 0. )
01054 {
01055 fpar[1] += 1. ;
01056 continue ;
01057 }
01058 if ( r )
01059 {
01060 return r ;
01061 }
01062 }
01063
01064 if ( labda <= LABMAXA )
01065 {
01066 for ( i = 0 ; i < *npar ; i++ )
01067 {
01068 fpar[i] = epar[i] ;
01069 }
01070 }
01071 if ( (fabs( chi2 - chi1 ) <= (tolerance * chi1)) || (labda > LABMAXA) )
01072 {
01073
01074
01075
01076
01077
01078 labda = LABMINA ;
01079 get_mat_edge ( xdat, xdim, ydat, wdat, ndat, fpar, epar) ;
01080 r = get_vec_edge ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
01081
01082 if ( r )
01083 {
01084 return r ;
01085 }
01086 for ( i = 0 ; i < (*npar) ; i++ )
01087 {
01088 epar[i] = 0.0 ;
01089 }
01090 chi2 = sqrt ( chi2 / (double) (nuse - nfree) ) ;
01091
01092 for ( i = 0 ; i < nfree ; i++ )
01093 {
01094 if ( (matrix1[i][i] <= 0.0) || (matrix2[i][i] <= 0.0) )
01095 {
01096 return -7 ;
01097 }
01098 epar[parptr[i]] = chi2 * sqrt( matrix2[i][i] ) / sqrt( matrix1[i][i] ) ;
01099 }
01100 found = 1 ;
01101 }
01102 }
01103 }
01104 return itc ;
01105 }
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131 int fitSlits1( OneImage * lineImage,
01132 FitParams ** par,
01133 float ** slit_pos,
01134 int box_length,
01135 float y_box )
01136 {
01137 float position[lineImage -> lx] ;
01138 int * edge, * edgeclean ;
01139 int * dummyedge ;
01140 int * pos_row, * pos_rowclean ;
01141 Vector * box_buffer ;
01142 float max_intensity ;
01143 float row_pos ;
01144 int col ;
01145 int i, j, k, m, n, ed ;
01146 int found, init1, init2 ;
01147 int line ;
01148 int column ;
01149 int slit_length ;
01150 int agreed ;
01151 int bad_line ;
01152 int margin ;
01153 int iters, xdim, ndat ;
01154 int numpar, its ;
01155 int * mpar ;
01156 float * xdat, * wdat ;
01157 float tol, lab ;
01158 float fitpar[NPAR] ;
01159 float dervpar[NPAR] ;
01160 float minval, maxval ;
01161
01162 slit_length = (int) sqrt (lineImage->lx) ;
01163
01164 if ( NullImage == lineImage )
01165 {
01166 cpl_msg_error( "fitSlits1","no line image given!\n" ) ;
01167 return -1 ;
01168 }
01169
01170 if ( NULL == par )
01171 {
01172 cpl_msg_error( "fitSlits1","no line fit parameters given!\n" ) ;
01173 return -2 ;
01174 }
01175
01176 if ( NULL == slit_pos )
01177 {
01178 cpl_msg_error( "fitSlits1","no position array allocated!\n" ) ;
01179 return -3 ;
01180 }
01181
01182 if ( box_length < slit_length ||
01183 box_length >= 3*slit_length )
01184 {
01185 cpl_msg_error( "fitSlits1","wrong fitting box length given!\n" ) ;
01186 return -4 ;
01187 }
01188
01189 if ( y_box <= 0. || y_box > 3. )
01190 {
01191 cpl_msg_error( "fitSlits1","wrong y box length given!\n" ) ;
01192 return -5 ;
01193 }
01194
01195
01196 edge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01197 dummyedge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01198 edgeclean = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ;
01199 pos_row = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01200 pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ;
01201
01202
01203
01204
01205
01206 agreed = -1 ;
01207 bad_line = -1 ;
01208 while( agreed == -1 )
01209 {
01210 found = -1 ;
01211 max_intensity = -FLT_MAX ;
01212 for ( col = 0 ; col < box_length ; col++ )
01213 {
01214 for ( i = 0 ; i < par[0]->n_params ; i++ )
01215 {
01216 if ( par[i]->column == col && par[i]->line != bad_line )
01217 {
01218 if ( par[i]->fit_par[0] > max_intensity )
01219 {
01220 if ( par[i]->fit_par[1] > 0. )
01221 {
01222 max_intensity = par[i]->fit_par[0] ;
01223 found = i ;
01224 }
01225 }
01226 }
01227 }
01228 }
01229
01230
01231
01232
01233
01234 line = par[found]->line ;
01235 column = par[found]->column ;
01236 row_pos = par[found]->fit_par[2] ;
01237 if ( found >= 0 && max_intensity > 0. )
01238 {
01239 for ( i = 0 ; i < par[0]->n_params ; i++ )
01240 {
01241 if ( par[i]->line == line-1 &&
01242 par[i]->column == column + slit_length )
01243 {
01244 if ( par[i]->fit_par[2] <= (row_pos + y_box) &&
01245 par[i]->fit_par[2] >= (row_pos - y_box) )
01246 {
01247 bad_line = line ;
01248 }
01249 }
01250 }
01251 if ( bad_line != line )
01252 {
01253 agreed = 1 ;
01254 break ;
01255 }
01256 }
01257 else
01258 {
01259 cpl_msg_error("fitSlits1","no emission line found in the first image columns\n") ;
01260 return -6 ;
01261 }
01262 }
01263
01264
01265 if ( agreed == -1 )
01266 {
01267 cpl_msg_error("fitSlits1","no emission line found in the first image columns\n") ;
01268 return -6 ;
01269 }
01270
01271
01272 n = 0 ;
01273 ed = 0 ;
01274 for ( col = 0 ; col < lineImage -> lx ; col++ )
01275 {
01276 for ( i = 0 ; i < par[0]->n_params ; i++ )
01277 {
01278 if ( par[i]->column == col && par[i] -> line == line )
01279 {
01280 if ( par[i]->fit_par[0] > 0. && par[i]->fit_par[1] > 0. )
01281 {
01282 position[n] = par[i]->fit_par[2] ;
01283 if ( n > 0 && fabs(position[n] - position[n-1]) > y_box )
01284 {
01285 edge[ed] = col ;
01286 pos_row[ed] = nint( position[n-1] ) ;
01287 ed++ ;
01288 if ( col >= lineImage -> lx - slit_length - 5 )
01289 {
01290 pos_row[ed] = nint( position[n] ) ;
01291 }
01292 }
01293 n++ ;
01294 }
01295 }
01296 }
01297 }
01298 if ( ed < (slit_length - 1) )
01299 {
01300 cpl_msg_error("fitSlits1","not enough slitlets found\n") ;
01301 return -7 ;
01302 }
01303
01304
01305 for ( i = 1 ; i <= ed ; i ++ )
01306 {
01307 if (dummyedge[i-1] != -1)
01308 {
01309 dummyedge[i-1] = edge[i-1] ;
01310 }
01311 if ( (edge[i] - edge[i-1]) < slit_length - 3 ||
01312 (edge[i] - edge[i-1]) > slit_length + 3 )
01313 {
01314 dummyedge[i] = -1 ;
01315 }
01316 if ( (edge[i+1] - edge[i]) < slit_length - 3 ||
01317 (edge[i+1] - edge[i]) > slit_length + 3 )
01318 {
01319 dummyedge[i+1] = -1 ;
01320 }
01321 }
01322
01323 k = 0 ;
01324 for ( i = 0 ; i < ed ; i++ )
01325 {
01326 if ( dummyedge[i] != -1 && dummyedge[i] != 0 )
01327 {
01328 edgeclean[k] = dummyedge[i] ;
01329 pos_rowclean[k] = pos_row[i] ;
01330 k++ ;
01331 if( edgeclean[k-1] > (lineImage -> lx - slit_length - 6 ) )
01332 {
01333 pos_rowclean[k] = pos_row[ed] ;
01334 }
01335 }
01336 }
01337
01338 if ( k != slit_length - 1 )
01339 {
01340 cpl_msg_error("fitSlits1","not enough clean slitlets found\n") ;
01341 return -7 ;
01342 }
01343
01344
01345 margin = ( box_length - slit_length ) / 2 ;
01346
01347
01348 for ( j = 0 ; j < k ; j++ )
01349 {
01350 m = 0 ;
01351 if ( j == 0 )
01352 {
01353 box_buffer = newVector( edgeclean[0] + margin ) ;
01354 for( col = 0 ; col < edgeclean[0] + margin ; col++ )
01355 {
01356 box_buffer->data[m] = lineImage -> data[col + lineImage->lx*pos_rowclean[0]] ;
01357 m++ ;
01358 }
01359 fitpar[0] = 3. ;
01360 fitpar[1] = 5. ;
01361 fitpar[2] = (float) edgeclean[0] - 1. ;
01362 fitpar[3] = (float) edgeclean[0] + 1. ;
01363
01364 }
01365 else if ( j < k - 1 )
01366 {
01367 box_buffer = newVector( edgeclean[j] - edgeclean[j-1] + 2*margin ) ;
01368 for ( col = edgeclean[j - 1] - margin ; col < edgeclean[j] + margin ; col++ )
01369 {
01370 box_buffer->data[m] = lineImage -> data[col + lineImage->lx*pos_rowclean[j]] ;
01371 m++ ;
01372 }
01373 fitpar[0] = (float) margin - 1. ;
01374 fitpar[1] = (float) margin + 1. ;
01375 fitpar[2] = (float) (edgeclean[j] - edgeclean[j-1] + margin) - 1. ;
01376 fitpar[3] = (float) (edgeclean[j] - edgeclean[j-1] + margin) + 1. ;
01377 }
01378
01379 else
01380 {
01381 box_buffer = newVector( lineImage -> lx - edgeclean[k-1] + margin ) ;
01382 for ( col = edgeclean[k - 1] - margin ; col < lineImage -> lx ; col++ )
01383 {
01384 box_buffer->data[m] = lineImage -> data[col + lineImage->lx*pos_rowclean[k]] ;
01385 m++ ;
01386 }
01387 fitpar[0] = (float) margin - 1. ;
01388 fitpar[1] = (float) margin + 1. ;
01389 fitpar[2] = (float) (lineImage -> lx - edgeclean[k-1] + margin) - 3. ;
01390 fitpar[3] = (float) (lineImage -> lx - edgeclean[k-1] + margin) - 1. ;
01391 }
01392
01393 xdat = (float *) cpl_calloc( box_buffer -> n_elements, sizeof (float) ) ;
01394 wdat = (float *) cpl_calloc( box_buffer -> n_elements, sizeof (float) ) ;
01395 mpar = (int *) cpl_calloc( NPAR, sizeof (int) ) ;
01396
01397
01398 minval = FLT_MAX ;
01399 maxval = -FLT_MAX ;
01400 for ( i = 0 ; i < box_buffer->n_elements ; i++ )
01401 {
01402 xdat[i] = i ;
01403 wdat[i] = 1.0 ;
01404 if ( box_buffer -> data[i] < minval )
01405 {
01406 minval = box_buffer -> data[i] ;
01407 }
01408 if ( box_buffer -> data[i] > maxval )
01409 {
01410 maxval = box_buffer -> data[i] ;
01411 }
01412 }
01413
01414 fitpar[4] = minval ;
01415 fitpar[5] = minval ;
01416 fitpar[6] = maxval ;
01417 fitpar[7] = maxval ;
01418
01419
01420 init1 = -1 ;
01421 init2 = -1 ;
01422 for ( i = 0 ; i < box_buffer->n_elements ; i++ )
01423 {
01424 if ( box_buffer -> data[i] >= ( maxval - minval ) / 2. )
01425 {
01426 init1 = i ;
01427 break ;
01428 }
01429 }
01430
01431 for ( i = box_buffer->n_elements - 1 ; i >= 0 ; i-- )
01432 {
01433 if ( box_buffer -> data[i] >= ( maxval + minval ) / 2. )
01434 {
01435 init2 = i ;
01436 break ;
01437 }
01438 }
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452 for ( i = 0 ; i < NPAR ; i++ )
01453 {
01454 mpar[i] = 1 ;
01455 dervpar[i] = 0. ;
01456 }
01457
01458 xdim = XDIMA ;
01459 ndat = box_buffer -> n_elements ;
01460 numpar = NPAR ;
01461 tol = TOLA ;
01462 lab = LABA ;
01463 its = ITSA ;
01464
01465
01466 if ( 0 > ( iters = lsqfit_edge( xdat, &xdim, box_buffer -> data, wdat, &ndat, fitpar,
01467 dervpar, mpar, &numpar, &tol, &its, &lab )) )
01468 {
01469 cpl_msg_error ("lsqfit:"," least squares fit failed, error no.: %d\n", iters) ;
01470 return -8 ;
01471 }
01472
01473
01474
01475 if ( j == 0 )
01476 {
01477 slit_pos[0][0] = (fitpar[0] + fitpar[1]) / 2. ;
01478 slit_pos[0][1] = (fitpar[2] + fitpar[3]) / 2. ;
01479 }
01480 else
01481 {
01482 slit_pos[j][0] = (fitpar[0] + fitpar[1]) / 2. + (float)edgeclean[j-1] - (float)margin ;
01483 slit_pos[j][1] = (fitpar[2] + fitpar[3]) / 2. + (float)edgeclean[j-1] - (float)margin ;
01484 }
01485
01486 slit_pos[k][0] = (fitpar[0] + fitpar[1]) / 2. + (float)edgeclean[k-1] - (float)margin ;
01487 slit_pos[k][1] = (fitpar[2] + fitpar[3]) / 2. + (float)edgeclean[k-1] - (float)margin ;
01488
01489 destroyVector ( box_buffer ) ;
01490 cpl_free( xdat ) ;
01491 cpl_free( wdat ) ;
01492 cpl_free( mpar ) ;
01493 }
01494
01495 cpl_free( edge ) ;
01496 cpl_free( pos_row ) ;
01497 cpl_free( edgeclean ) ;
01498 cpl_free( dummyedge ) ;
01499 cpl_free( pos_rowclean ) ;
01500 return 0 ;
01501 }
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529 int fitSlits( OneImage * lineImage,
01530 FitParams ** par,
01531 float ** slit_pos,
01532 int box_length,
01533 float y_box,
01534 float slope_width )
01535 {
01536 float position[lineImage -> lx] ;
01537 int * edge, * edgeclean ;
01538 int * dummyedge ;
01539 int * pos_row, * pos_rowclean ;
01540 Vector * box_buffer ;
01541 float max_intensity ;
01542 float row_pos ;
01543 int col ;
01544 int i, j, k, m, n, ed ;
01545 int found ;
01546 int line ;
01547 int column ;
01548 int slit_length ;
01549 int agreed ;
01550 int bad_line ;
01551 int margin ;
01552 int iters, xdim, ndat ;
01553 int numpar, its ;
01554 int * mpar ;
01555 float * xdat, * wdat ;
01556 float tol, lab ;
01557 float fitpar[NPAR] ;
01558 float dervpar[NPAR] ;
01559 float minval, maxval ;
01560
01561 slit_length = (int) sqrt (lineImage->lx) ;
01562
01563 if ( NullImage == lineImage )
01564 {
01565 cpl_msg_error( "fitSlits","no line image given!\n" ) ;
01566 return -1 ;
01567 }
01568
01569 if ( NULL == par )
01570 {
01571 cpl_msg_error( "fitSlits","no line fit parameters given!\n" ) ;
01572 return -2 ;
01573 }
01574
01575 if ( NULL == slit_pos )
01576 {
01577 cpl_msg_error( "fitSlits","no position array allocated!\n" ) ;
01578 return -3 ;
01579 }
01580
01581 if ( box_length < slit_length ||
01582 box_length >= 3*slit_length )
01583 {
01584 cpl_msg_error( "fitSlits","wrong fitting box length given!\n" ) ;
01585 return -4 ;
01586 }
01587
01588 if ( y_box <= 0. || y_box > 3. )
01589 {
01590 cpl_msg_error( "fitSlits","wrong y box length given!\n" ) ;
01591 return -5 ;
01592 }
01593
01594 if ( slope_width <= 0. )
01595 {
01596 cpl_msg_error("fitSlits","wrong width of linear slope given!\n") ;
01597 return -6 ;
01598 }
01599
01600
01601 slopewidth = slope_width ;
01602
01603
01604 edge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01605 dummyedge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01606 edgeclean = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ;
01607 pos_row = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01608 pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ;
01609
01610
01611
01612
01613
01614 agreed = -1 ;
01615 bad_line = -1 ;
01616 while( agreed == -1 )
01617 {
01618 found = -1 ;
01619 max_intensity = -FLT_MAX ;
01620 for ( col = 0 ; col < box_length ; col++ )
01621 {
01622 for ( i = 0 ; i < par[0]->n_params ; i++ )
01623 {
01624 if ( par[i]->column == col && par[i]->line != bad_line )
01625 {
01626 if ( par[i]->fit_par[0] > max_intensity )
01627 {
01628 if ( par[i]->fit_par[1] > 0. )
01629 {
01630 max_intensity = par[i]->fit_par[0] ;
01631 found = i ;
01632 }
01633 }
01634 }
01635 }
01636 }
01637
01638
01639
01640
01641
01642 line = par[found]->line ;
01643 column = par[found]->column ;
01644 row_pos = par[found]->fit_par[2] ;
01645 if ( found >= 0 && max_intensity > 0. )
01646 {
01647 for ( i = 0 ; i < par[0]->n_params ; i++ )
01648 {
01649 if ( par[i]->line == line-1 &&
01650 par[i]->column == column + slit_length )
01651 {
01652 if ( par[i]->fit_par[2] <= (row_pos + y_box) &&
01653 par[i]->fit_par[2] >= (row_pos - y_box) )
01654 {
01655 bad_line = line ;
01656 }
01657 }
01658 }
01659 if ( bad_line != line )
01660 {
01661 agreed = 1 ;
01662 break ;
01663 }
01664 }
01665 else
01666 {
01667 cpl_msg_error("fitSlits","no emission line found in the first image columns\n") ;
01668 return -7 ;
01669 }
01670 }
01671
01672
01673 if ( agreed == -1 )
01674 {
01675 cpl_msg_error("fitSlits","no emission line found in the first image columns\n") ;
01676 return -7 ;
01677 }
01678
01679
01680 n = 0 ;
01681 ed = 0 ;
01682 for ( col = 0 ; col < lineImage -> lx ; col++ )
01683 {
01684 for ( i = 0 ; i < par[0]->n_params ; i++ )
01685 {
01686 if ( par[i]->column == col && par[i] -> line == line )
01687 {
01688 if ( par[i]->fit_par[0] > 0. && par[i]->fit_par[1] > 0. )
01689 {
01690 position[n] = par[i]->fit_par[2] ;
01691 if ( n > 0 && fabs(position[n] - position[n-1]) > y_box )
01692 {
01693 edge[ed] = col ;
01694 pos_row[ed] = nint( position[n-1] ) ;
01695 ed++ ;
01696 if ( col >= lineImage -> lx - slit_length - 5 )
01697 {
01698 pos_row[ed] = nint( position[n] ) ;
01699 }
01700 }
01701 n++ ;
01702 }
01703 }
01704 }
01705 }
01706 if ( ed < (slit_length - 1) )
01707 {
01708 cpl_msg_error("fitSlits","not enough slitlets found\n") ;
01709 return -8 ;
01710 }
01711
01712
01713 for ( i = 1 ; i <= ed ; i ++ )
01714 {
01715 if (dummyedge[i-1] != -1)
01716 {
01717 dummyedge[i-1] = edge[i-1] ;
01718 }
01719 if ( (edge[i] - edge[i-1]) < slit_length - 3 ||
01720 (edge[i] - edge[i-1]) > slit_length + 3 )
01721 {
01722 dummyedge[i] = -1 ;
01723 }
01724 if ( (edge[i+1] - edge[i]) < slit_length - 3 ||
01725 (edge[i+1] - edge[i]) > slit_length + 3 )
01726 {
01727 dummyedge[i+1] = -1 ;
01728 }
01729 }
01730
01731 k = 0 ;
01732 for ( i = 0 ; i < ed ; i++ )
01733 {
01734 if ( dummyedge[i] != -1 && dummyedge[i] != 0 )
01735 {
01736 edgeclean[k] = dummyedge[i] ;
01737 pos_rowclean[k] = pos_row[i] ;
01738 k++ ;
01739 if( edgeclean[k-1] > (lineImage -> lx - slit_length - 6 ) )
01740 {
01741 pos_rowclean[k] = pos_row[ed] ;
01742 }
01743 }
01744 }
01745
01746 if ( k != slit_length - 1 )
01747 {
01748 cpl_msg_error ("fitSlits","not enough clean slitlets found\n") ;
01749 return -7 ;
01750 }
01751
01752
01753 margin = ( box_length - slit_length ) / 2 ;
01754
01755
01756 for ( j = 0 ; j < k ; j++ )
01757 {
01758 m = 0 ;
01759 if ( j == 0 )
01760 {
01761 box_buffer = newVector( edgeclean[0] + margin ) ;
01762 for( col = 0 ; col < edgeclean[0] + margin ; col++ )
01763 {
01764 box_buffer->data[m] = lineImage -> data[col + lineImage->lx*pos_rowclean[0]] ;
01765 m++ ;
01766 }
01767
01768 fitpar[0] = 1. ;
01769 fitpar[1] = (float)edgeclean[0] ;
01770
01771 }
01772 else if ( j < k - 1 )
01773 {
01774 box_buffer = newVector( edgeclean[j] - edgeclean[j-1] + 2*margin ) ;
01775 for ( col = edgeclean[j - 1] - margin ; col < edgeclean[j] + margin ; col++ )
01776 {
01777 box_buffer->data[m] = lineImage -> data[col + lineImage->lx*pos_rowclean[j]] ;
01778 m++ ;
01779 }
01780
01781 fitpar[0] = (float)margin ;
01782 fitpar[1] = (float)(edgeclean[j] - edgeclean[j-1] + margin ) ;
01783 }
01784
01785 else
01786 {
01787 box_buffer = newVector( lineImage -> lx - edgeclean[k-1] + margin ) ;
01788 for ( col = edgeclean[k - 1] - margin ; col < lineImage -> lx ; col++ )
01789 {
01790 box_buffer->data[m] = lineImage -> data[col + lineImage->lx*pos_rowclean[k]] ;
01791 m++ ;
01792 }
01793
01794 fitpar[0] = (float)margin ;
01795 fitpar[1] = (float)(m - 1) ;
01796 }
01797
01798 xdat = (float *) cpl_calloc( box_buffer -> n_elements, sizeof (float) ) ;
01799 wdat = (float *) cpl_calloc( box_buffer -> n_elements, sizeof (float) ) ;
01800 mpar = (int *) cpl_calloc( NPAR, sizeof (int) ) ;
01801
01802
01803 minval = FLT_MAX ;
01804 maxval = -FLT_MAX ;
01805 for ( i = 0 ; i < box_buffer->n_elements ; i++ )
01806 {
01807 xdat[i] = i ;
01808 wdat[i] = 1.0 ;
01809 if ( box_buffer -> data[i] < minval )
01810 {
01811 minval = box_buffer -> data[i] ;
01812 }
01813 if ( box_buffer -> data[i] > maxval )
01814 {
01815 maxval = box_buffer -> data[i] ;
01816 }
01817 }
01818
01819 for ( i = 0 ; i < NPAR ; i++ )
01820 {
01821 mpar[i] = 1 ;
01822 dervpar[i] = 0. ;
01823 }
01824
01825 xdim = XDIMA ;
01826 ndat = box_buffer -> n_elements ;
01827 numpar = NPAR ;
01828 tol = TOLA ;
01829 lab = LABA ;
01830 its = ITSA ;
01831
01832 fitpar[2] = minval ;
01833 fitpar[3] = minval ;
01834 fitpar[4] = maxval ;
01835
01836
01837 if ( 0 > ( iters = lsqfit_edge( xdat, &xdim, box_buffer -> data, wdat, &ndat, fitpar,
01838 dervpar, mpar, &numpar, &tol, &its, &lab )) )
01839 {
01840 cpl_msg_error ("lsqfit:"," least squares fit failed, error no.: %d\n", iters) ;
01841 return -9 ;
01842 }
01843
01844
01845
01846 if ( j == 0 )
01847 {
01848 slit_pos[0][0] = fitpar[0] + slopewidth/2. ;
01849 slit_pos[0][1] = fitpar[1] - slopewidth/2. ;
01850 }
01851 else
01852 {
01853 slit_pos[j][0] = fitpar[0] + slopewidth/2. + (float)edgeclean[j-1] - (float)margin ;
01854 slit_pos[j][1] = fitpar[1] - slopewidth/2. + (float)edgeclean[j-1] - (float)margin ;
01855 }
01856
01857 slit_pos[k][0] = fitpar[0] + slopewidth/2. + (float)edgeclean[k-1] - (float)margin ;
01858 slit_pos[k][1] = fitpar[1] - slopewidth/2. + (float)edgeclean[k-1] - (float)margin ;
01859
01860 destroyVector ( box_buffer ) ;
01861 cpl_free( xdat ) ;
01862 cpl_free( wdat ) ;
01863 cpl_free( mpar ) ;
01864 }
01865
01866 cpl_free( edge ) ;
01867 cpl_free( pos_row ) ;
01868 cpl_free( edgeclean ) ;
01869 cpl_free( dummyedge ) ;
01870 cpl_free( pos_rowclean ) ;
01871 return 0 ;
01872 }
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908 int fitSlits2( OneImage * lineImage,
01909 FitParams ** par,
01910 float ** slit_pos,
01911 int box_length,
01912 float y_box,
01913 float diff_tol )
01914 {
01915 float position[lineImage -> lx] ;
01916 int * edge, * edgeclean ;
01917 int * dummyedge ;
01918 int * pos_row, * pos_rowclean ;
01919 Vector * box_buffer ;
01920 Vector * half_buffer ;
01921 float max_intensity ;
01922 float row_pos ;
01923 int col ;
01924 int i, j, k, m, n, ed ;
01925 int found, init1 ;
01926 int line ;
01927 int nel, n_right, left_right ;
01928 int column ;
01929 int slit_length ;
01930 int agreed ;
01931 int bad_line ;
01932 int margin ;
01933 int iters, xdim, ndat ;
01934 int numpar, its ;
01935 int * mpar ;
01936 float * xdat, * wdat ;
01937 float tol, lab ;
01938 float fitpar[NPAR] ;
01939 float dervpar[NPAR] ;
01940 float minval, maxval ;
01941 float pos, last_pos ;
01942
01943 slit_length = (int) sqrt (lineImage->lx) ;
01944
01945 if ( NullImage == lineImage )
01946 {
01947 cpl_msg_error( "fitSlits2","no line image given!\n" ) ;
01948 return -1 ;
01949 }
01950
01951 if ( NULL == par )
01952 {
01953 cpl_msg_error( "fitSlits2","no line fit parameters given!\n" ) ;
01954 return -2 ;
01955 }
01956
01957 if ( NULL == slit_pos )
01958 {
01959 cpl_msg_error( "fitSlits2","no position array allocated!\n" ) ;
01960 return -3 ;
01961 }
01962
01963 if ( box_length < slit_length ||
01964 box_length >= 3*slit_length )
01965 {
01966 cpl_msg_error( "fitSlits2","wrong fitting box length given!\n" ) ;
01967 return -4 ;
01968 }
01969
01970 if ( y_box <= 0. || y_box > 3. )
01971 {
01972 cpl_msg_error( "fitSlits2","wrong y box length given!\n" ) ;
01973 return -5 ;
01974 }
01975
01976 if ( diff_tol < 1. )
01977 {
01978 cpl_msg_error( "fitSlits2","diff_tol too small!\n" ) ;
01979 return -6 ;
01980 }
01981
01982
01983 edge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01984 dummyedge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01985 edgeclean = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ;
01986 pos_row = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
01987 pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ;
01988
01989
01990
01991
01992
01993 agreed = -1 ;
01994 bad_line = -1 ;
01995 while( agreed == -1 )
01996 {
01997 found = -1 ;
01998 max_intensity = -FLT_MAX ;
01999 for ( col = 0 ; col < box_length ; col++ )
02000 {
02001 for ( i = 0 ; i < par[0]->n_params ; i++ )
02002 {
02003 if ( par[i]->column == col && par[i]->line != bad_line )
02004 {
02005 if ( par[i]->fit_par[0] > max_intensity )
02006 {
02007 if ( par[i]->fit_par[1] > 0. )
02008 {
02009 max_intensity = par[i]->fit_par[0] ;
02010 found = i ;
02011 }
02012 }
02013 }
02014 }
02015 }
02016
02017
02018
02019
02020
02021 line = par[found]->line ;
02022 column = par[found]->column ;
02023 row_pos = par[found]->fit_par[2] ;
02024 if ( found >= 0 && max_intensity > 0. )
02025 {
02026 for ( i = 0 ; i < par[0]->n_params ; i++ )
02027 {
02028 if ( par[i]->line == line-1 &&
02029 par[i]->column == column + slit_length )
02030 {
02031 if ( par[i]->fit_par[2] <= (row_pos + y_box) &&
02032 par[i]->fit_par[2] >= (row_pos - y_box) )
02033 {
02034 bad_line = line ;
02035 }
02036 }
02037 }
02038 if ( bad_line != line )
02039 {
02040 agreed = 1 ;
02041 break ;
02042 }
02043 }
02044 else
02045 {
02046 cpl_msg_error("fitSlits2","no emission line found in the first image columns\n") ;
02047 return -7 ;
02048 }
02049 }
02050
02051
02052 if ( agreed == -1 )
02053 {
02054 cpl_msg_error("fitSlits2","no emission line found in the first image columns\n") ;
02055 return -7 ;
02056 }
02057
02058
02059 n = 0 ;
02060 ed = 0 ;
02061 for ( col = 0 ; col < lineImage -> lx ; col++ )
02062 {
02063 for ( i = 0 ; i < par[0]->n_params ; i++ )
02064 {
02065 if ( par[i]->column == col && par[i]->line == line )
02066 {
02067 if ( par[i]->fit_par[0] > 0. && par[i]->fit_par[1] > 0. )
02068 {
02069 position[n] = par[i]->fit_par[2] ;
02070 if ( n > 0 && fabs(position[n] - position[n-1]) > y_box )
02071 {
02072 edge[ed] = col ;
02073 pos_row[ed] = nint( position[n-1] ) ;
02074 ed++ ;
02075 if ( col >= lineImage -> lx - slit_length - 5 )
02076 {
02077 pos_row[ed] = nint( position[n] ) ;
02078 }
02079 }
02080 n++ ;
02081 }
02082 }
02083 }
02084 }
02085 if ( ed < (slit_length - 1) )
02086 {
02087 cpl_msg_error("fitSlits2","not enough slitlets found\n") ;
02088 return -8 ;
02089 }
02090
02091
02092 for ( i = 1 ; i <= ed ; i ++ )
02093 {
02094 if (dummyedge[i-1] != -1)
02095 {
02096 dummyedge[i-1] = edge[i-1] ;
02097 }
02098 if ( (edge[i] - edge[i-1]) < slit_length - 3 ||
02099 (edge[i] - edge[i-1]) > slit_length + 3 )
02100 {
02101 dummyedge[i] = -1 ;
02102 }
02103 if ( (edge[i+1] - edge[i]) < slit_length - 3 ||
02104 (edge[i+1] - edge[i]) > slit_length + 3 )
02105 {
02106 dummyedge[i+1] = -1 ;
02107 }
02108 }
02109
02110 k = 0 ;
02111 for ( i = 0 ; i < ed ; i++ )
02112 {
02113 if ( dummyedge[i] != -1 && dummyedge[i] != 0 )
02114 {
02115 edgeclean[k] = dummyedge[i] ;
02116 pos_rowclean[k] = pos_row[i] ;
02117 k++ ;
02118 if( edgeclean[k-1] > (lineImage -> lx - slit_length - 6 ) )
02119 {
02120 pos_rowclean[k] = pos_row[ed] ;
02121 }
02122 }
02123 }
02124
02125 if ( k != slit_length - 1 )
02126 {
02127 cpl_msg_error("fitSlits2","not enough clean slitlets found\n") ;
02128 return -7 ;
02129 }
02130
02131
02132 margin = ( box_length - slit_length ) / 2 ;
02133
02134
02135 for ( j = 0 ; j <= k ; j++ )
02136 {
02137 m = 0 ;
02138 if ( j == 0 )
02139 {
02140 box_buffer = newVector( edgeclean[0] + margin ) ;
02141 for( col = 0 ; col < edgeclean[0] + margin ; col++ )
02142 {
02143 box_buffer->data[m] = lineImage -> data[col + lineImage->lx*pos_rowclean[0]] ;
02144 m++ ;
02145 }
02146 }
02147 else if ( j < k )
02148 {
02149 box_buffer = newVector( edgeclean[j] - edgeclean[j-1] + 2*margin ) ;
02150 for ( col = edgeclean[j - 1] - margin ; col < edgeclean[j] + margin ; col++ )
02151 {
02152 box_buffer->data[m] = lineImage -> data[col + lineImage->lx*pos_rowclean[j]] ;
02153 m++ ;
02154 }
02155 }
02156 else
02157 {
02158 box_buffer = newVector( lineImage -> lx - edgeclean[k-1] + margin ) ;
02159 for ( col = edgeclean[k - 1] - margin ; col < lineImage -> lx ; col++ )
02160 {
02161 box_buffer->data[m] = lineImage -> data[col + lineImage->lx*pos_rowclean[k]] ;
02162 m++ ;
02163 }
02164 }
02165
02166 for ( left_right = 0 ; left_right <= 1 ; left_right++ )
02167 {
02168 nel = 0 ;
02169 if ( left_right == 0 )
02170 {
02171 nel = box_buffer -> n_elements / 2 ;
02172 }
02173 else
02174 {
02175 if ( box_buffer -> n_elements % 2 == 0 )
02176 {
02177 nel = box_buffer -> n_elements / 2 ;
02178 }
02179 else
02180 {
02181 nel = box_buffer -> n_elements / 2 + 1 ;
02182 }
02183 }
02184
02185
02186 half_buffer = newVector( nel ) ;
02187 if ( left_right == 0 )
02188 {
02189 for ( i = 0 ; i < nel ; i++ )
02190 {
02191 half_buffer -> data[i] = box_buffer -> data[i] ;
02192 }
02193 }
02194 else
02195 {
02196 n_right = 0 ;
02197 for ( i = box_buffer -> n_elements - 1 ; i >= box_buffer -> n_elements - nel ; i-- )
02198 {
02199 half_buffer -> data[n_right] = box_buffer -> data[i] ;
02200 n_right++ ;
02201 }
02202 }
02203
02204 xdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
02205 wdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
02206 mpar = (int *) cpl_calloc( NPAR, sizeof (int) ) ;
02207
02208
02209 minval = FLT_MAX ;
02210 maxval = -FLT_MAX ;
02211 for ( i = 0 ; i < nel ; i++ )
02212 {
02213 xdat[i] = i ;
02214 wdat[i] = 1.0 ;
02215 if ( half_buffer -> data[i] < minval )
02216 {
02217 minval = half_buffer -> data[i] ;
02218 }
02219 if ( half_buffer -> data[i] > maxval )
02220 {
02221 maxval = half_buffer -> data[i] ;
02222 }
02223 }
02224
02225 fitpar[2] = minval ;
02226 fitpar[3] = maxval ;
02227
02228
02229 init1 = -1 ;
02230 for ( i = 0 ; i < nel ; i++ )
02231 {
02232 if ( half_buffer -> data[i] >= ( maxval + minval ) / 2. )
02233 {
02234 init1 = i ;
02235 break ;
02236 }
02237 }
02238
02239
02240 if ( init1 != -1 )
02241 {
02242 fitpar[0] = ((float)init1 - 1.) ;
02243 fitpar[1] = ((float)init1 + 1.) ;
02244 }
02245
02246 for ( i = 0 ; i < NPAR ; i++ )
02247 {
02248 mpar[i] = 1 ;
02249 dervpar[i] = 0. ;
02250 }
02251
02252 xdim = XDIMA ;
02253 ndat = nel ;
02254 numpar = NPAR ;
02255 tol = TOLA ;
02256 lab = LABA ;
02257 its = ITSA ;
02258
02259
02260 if ( 0 > ( iters = lsqfit_edge( xdat, &xdim, half_buffer -> data, wdat, &ndat, fitpar,
02261 dervpar, mpar, &numpar, &tol, &its, &lab )) )
02262 {
02263
02264 cpl_msg_warning ("lsqfit:"," least squares fit failed, error no.: %d in slitlet: %d\n", iters, j) ;
02265 fitpar[0] = ((float)init1 - 1.) ;
02266 fitpar[1] = ((float)init1 + 1.) ;
02267 }
02268
02269 pos = (fitpar[0] + fitpar[1]) / 2. ;
02270
02271
02272
02273
02274
02275
02276
02277 if ( left_right == 0 )
02278 {
02279
02280 if ( j == 0 )
02281 {
02282 if ( fabs(pos - ((float)edgeclean[0] - 1. - (float)slit_length)) < diff_tol )
02283 {
02284 slit_pos[0][0] = pos ;
02285 }
02286 else
02287 {
02288 cpl_msg_warning("lsqfit","something wrong with fitted left position of slitlet 0\n") ;
02289 if ( (float) edgeclean[0] - 1. - (float)slit_length < 0. )
02290 {
02291 slit_pos[0][0] = 0. ;
02292 }
02293 else
02294 {
02295 slit_pos[0][0] = (float)edgeclean[0] - 1. - (float)slit_length ;
02296 }
02297 }
02298 }
02299 else if ( j < k )
02300 {
02301 if ( fabs( pos - (float)margin ) < diff_tol )
02302 {
02303 slit_pos[j][0] = pos + (float)edgeclean[j-1] - (float)margin ;
02304 }
02305 else
02306 {
02307 cpl_msg_warning("lsqfit","something wrong with fitted left position of slitlet %d\n", j) ;
02308 slit_pos[j][0] = (float)edgeclean[j-1] - 1. ;
02309 }
02310 }
02311 else
02312 {
02313 if ( fabs( pos - (float)margin ) < diff_tol )
02314 {
02315 slit_pos[k][0] = pos + (float)edgeclean[k-1] - (float)margin ;
02316 }
02317 else
02318 {
02319 cpl_msg_warning("lsqfit","something wrong with fitted left position of slitlet %d\n", j) ;
02320 slit_pos[k][0] = (float)edgeclean[k-1] - 1. ;
02321 }
02322 }
02323 }
02324 else
02325 {
02326
02327 if ( j == 0 )
02328 {
02329 if ( fabs( (float)box_buffer->n_elements - pos - (float)edgeclean[0] ) < diff_tol )
02330 {
02331 slit_pos[0][1] = (float)(box_buffer->n_elements - 1) - pos ;
02332 }
02333 else
02334 {
02335 cpl_msg_warning("lsqfit","something wrong with fitted right position of slitlet 0\n") ;
02336 slit_pos[0][1] = (float)edgeclean[0] - 1. ;
02337 }
02338 }
02339 else if ( j < k )
02340 {
02341 if ( fabs( (float)box_buffer->n_elements - pos
02342 + (float)edgeclean[j-1] - (float)margin - (float)edgeclean[j] ) < diff_tol )
02343 {
02344 slit_pos[j][1] = (float)(box_buffer->n_elements - 1) - pos
02345 + (float)edgeclean[j-1] - (float)margin ;
02346 }
02347 else
02348 {
02349 cpl_msg_warning("lsqfit","something wrong with fitted right position of slitlet %d\n", j) ;
02350 slit_pos[j][1] = (float)edgeclean[j] - 1. ;
02351 }
02352 }
02353 else
02354 {
02355 if ( edgeclean[k-1] + slit_length > lineImage -> lx )
02356 {
02357 last_pos = (float)(lineImage -> lx - 1) ;
02358 }
02359 else
02360 {
02361 last_pos = (float)(edgeclean[k-1] - 1 + slit_length) ;
02362 }
02363 if ( fabs( (float)(box_buffer->n_elements - 1) - pos
02364 + (float)edgeclean[k-1] - (float)margin - last_pos ) < diff_tol )
02365 {
02366 slit_pos[k][1] = (float)(box_buffer->n_elements - 1) - pos
02367 + (float)edgeclean[k-1] - (float)margin ;
02368 }
02369 else
02370 {
02371 cpl_msg_warning("lsqfit","something wrong with fitted right position of slitlet %d\n", j) ;
02372 slit_pos[k][1] = last_pos ;
02373 }
02374 }
02375 }
02376
02377 destroyVector ( half_buffer ) ;
02378 cpl_free( xdat ) ;
02379 cpl_free( wdat ) ;
02380 cpl_free( mpar ) ;
02381 }
02382 destroyVector ( box_buffer ) ;
02383 }
02384
02385 cpl_free( edge ) ;
02386 cpl_free( pos_row ) ;
02387 cpl_free( edgeclean ) ;
02388 cpl_free( dummyedge ) ;
02389 cpl_free( pos_rowclean ) ;
02390 return 0 ;
02391 }
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409
02410
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02423
02424
02425
02426
02427 int fitSlitsEdge( OneImage * lineImage,
02428 FitParams ** par,
02429 float ** slit_pos,
02430 int box_length,
02431 float y_box,
02432 float diff_tol )
02433 {
02434 float position[lineImage -> lx] ;
02435 int * edge, * edgeclean ;
02436 int * dummyedge ;
02437 int * pos_row, * pos_rowclean ;
02438 Vector * box_buffer ;
02439 Vector * half_buffer ;
02440 float max_intensity ;
02441 float row_pos ;
02442 int row, col ;
02443 int i, j, k, m, n, ed ;
02444 int found, init1 ;
02445 int line ;
02446 int nel, n_right, left_right ;
02447 int column ;
02448 int slit_length ;
02449 int agreed ;
02450 int bad_line ;
02451 int margin ;
02452 int iters, xdim, ndat ;
02453 int numpar, its ;
02454 int * mpar ;
02455 float * xdat, * wdat ;
02456 float tol, lab ;
02457 float fitpar[NPAR] ;
02458 float dervpar[NPAR] ;
02459 float minval, maxval ;
02460 float pos, last_pos ;
02461
02462 slit_length = (int) sqrt (lineImage->lx) ;
02463
02464 if ( NullImage == lineImage )
02465 {
02466 cpl_msg_error( "fitSlitsEdge:"," no line image given!\n" ) ;
02467 return -1 ;
02468 }
02469
02470 if ( NULL == par )
02471 {
02472 cpl_msg_error( "fitSlitsEdge:"," no line fit parameters given!\n" ) ;
02473 return -2 ;
02474 }
02475
02476 if ( NULL == slit_pos )
02477 {
02478 cpl_msg_error( "fitSlitsEdge:"," no position array allocated!\n" ) ;
02479 return -3 ;
02480 }
02481
02482 if ( box_length < 4 ||
02483 box_length >= 2*slit_length )
02484 {
02485 cpl_msg_error( "fitSlitsEdge:"," wrong fitting box length given!\n" ) ;
02486 cpl_msg_error( "fitSlitsEdge:"," Must be 4 <= box_length < %d \n",
02487 2*slit_length ) ;
02488 cpl_msg_error( "fitSlitsEdge:"," You have chosen box_length = %d \n",
02489 box_length);
02490
02491 return -4 ;
02492 }
02493
02494 if ( y_box <= 0. || y_box > 3. )
02495 {
02496 cpl_msg_error( "fitSlitsEdge:"," wrong y box length given!\n" ) ;
02497 cpl_msg_error( "fitSlitsEdge:"," y_box=%f not in range (0,3]!\n",y_box);
02498 return -5 ;
02499 }
02500
02501 if ( diff_tol < 1. )
02502 {
02503 cpl_msg_error( "fitSlitsEdge:"," diff_tol too small!\n" ) ;
02504 return -6 ;
02505 }
02506
02507
02508 edge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
02509 dummyedge = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
02510 edgeclean = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ;
02511 pos_row = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
02512 pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ;
02513
02514
02515
02516
02517
02518 agreed = -1 ;
02519 bad_line = -1 ;
02520 while( agreed == -1 )
02521 {
02522 found = -1 ;
02523 max_intensity = -FLT_MAX ;
02524 for ( col = 0 ; col < slit_length ; col++ )
02525 {
02526 for ( i = 0 ; i < par[0]->n_params ; i++ )
02527 {
02528 if ( par[i]->column == col && par[i] -> line != bad_line )
02529 {
02530 if ( par[i]->fit_par[0] > max_intensity )
02531 {
02532 if ( par[i]->fit_par[1] >= 1. && par[i]->fit_par[2] > 0. )
02533 {
02534 max_intensity = par[i]->fit_par[0] ;
02535 found = i ;
02536 }
02537 }
02538 }
02539 }
02540 }
02541
02542
02543
02544
02545
02546 line = par[found]->line ;
02547 column = par[found]->column ;
02548 row_pos = par[found]->fit_par[2] ;
02549 if ( found >= 0 && max_intensity > 0. )
02550 {
02551 for ( i = 0 ; i < par[0]->n_params ; i++ )
02552 {
02553 if ( par[i]->line == line-1 &&
02554 par[i]->column == column + slit_length )
02555 {
02556 if ( par[i]->fit_par[2] <= (row_pos + y_box) &&
02557 par[i]->fit_par[2] >= (row_pos - y_box) )
02558 {
02559 bad_line = line ;
02560 }
02561 }
02562 }
02563 if ( bad_line != line )
02564 {
02565 agreed = 1 ;
02566 break ;
02567 }
02568 }
02569 else
02570 {
02571 cpl_msg_error("fitSlitsEdge:"," no emission line found in the first image columns\n") ;
02572 cpl_free( edge ) ;
02573 cpl_free( pos_row ) ;
02574 cpl_free( edgeclean ) ;
02575 cpl_free( dummyedge ) ;
02576 cpl_free( pos_rowclean ) ;
02577 return -7 ;
02578 }
02579 }
02580
02581
02582 if ( agreed == -1 )
02583 {
02584 cpl_msg_error("fitSlitsEdge:"," no emission line found in the first image columns\n") ;
02585 cpl_free( edge ) ;
02586 cpl_free( pos_row ) ;
02587 cpl_free( edgeclean ) ;
02588 cpl_free( dummyedge ) ;
02589 cpl_free( pos_rowclean ) ;
02590 return -7 ;
02591 }
02592
02593
02594 n = 0 ;
02595 ed = 0 ;
02596 for ( col = 0 ; col < lineImage -> lx ; col++ )
02597 {
02598 for ( i = 0 ; i < par[0]->n_params ; i++ )
02599 {
02600 if ( par[i]->column == col && par[i]->line == line )
02601 {
02602 if ( par[i]->fit_par[0] > 0. && par[i]->fit_par[1] >= 1. && par[i]->fit_par[2] > 0. )
02603 {
02604 position[n] = par[i]->fit_par[2] ;
02605 if ( n > 0 && fabs(position[n] - position[n-1]) > y_box )
02606 {
02607 edge[ed] = col ;
02608 pos_row[ed] = nint( position[n-1] ) ;
02609 ed++ ;
02610 if ( col >= lineImage -> lx - slit_length - 5 )
02611 {
02612 pos_row[ed] = nint( position[n] ) ;
02613 }
02614 }
02615 n++ ;
02616 }
02617 }
02618 }
02619 }
02620 if ( ed < (slit_length - 1) )
02621 {
02622 cpl_msg_error("fitSlitsEdge:"," not enough slitlets found\n") ;
02623 cpl_free( edge ) ;
02624 cpl_free( pos_row ) ;
02625 cpl_free( edgeclean ) ;
02626 cpl_free( dummyedge ) ;
02627 cpl_free( pos_rowclean ) ;
02628 return -8 ;
02629 }
02630
02631
02632 for ( i = 1 ; i <= ed ; i ++ )
02633 {
02634 if ( i == ed )
02635 {
02636 if ( (edge[i-1] - edge[i-2]) < slit_length - 3 ||
02637 (edge[i-1] - edge[i-2]) > slit_length + 3 )
02638 {
02639 dummyedge[i-1] = -1 ;
02640 }
02641
02642 }
02643 if (dummyedge[i-1] != -1)
02644 {
02645 dummyedge[i-1] = edge[i-1] ;
02646 }
02647 else
02648 {
02649 continue ;
02650 }
02651 if ( i < ed )
02652 {
02653 if ( (edge[i] - edge[i-1]) < slit_length - 3 ||
02654 (edge[i] - edge[i-1]) > slit_length + 3 )
02655 {
02656 dummyedge[i] = -1 ;
02657 }
02658 }
02659 if ( i + 1 < ed && dummyedge[i] != -1 )
02660 {
02661 if ( (edge[i+1] - edge[i]) < slit_length - 3 ||
02662 (edge[i+1] - edge[i]) > slit_length + 3 )
02663 {
02664 dummyedge[i+1] = -1 ;
02665 }
02666 }
02667 }
02668
02669 k = 0 ;
02670 for ( i = 0 ; i < ed ; i++ )
02671 {
02672 if ( dummyedge[i] != -1 && dummyedge[i] != 0 )
02673 {
02674 edgeclean[k] = dummyedge[i] ;
02675 pos_rowclean[k] = pos_row[i] ;
02676 k++ ;
02677 if( edgeclean[k-1] > (lineImage -> lx - slit_length - 6 ) )
02678 {
02679 pos_rowclean[k] = pos_row[ed] ;
02680 }
02681 }
02682 }
02683
02684 if ( k != slit_length - 1 )
02685 {
02686 cpl_msg_error("fitSlitsEdge:"," not enough clean slitlets found\n") ;
02687 cpl_free( edge ) ;
02688 cpl_free( pos_row ) ;
02689 cpl_free( edgeclean ) ;
02690 cpl_free( dummyedge ) ;
02691 cpl_free( pos_rowclean ) ;
02692 return -8 ;
02693 }
02694
02695
02696 margin = box_length / 2 ;
02697
02698
02699
02700
02701
02702 for ( j = 0 ; j <= k ; j++ )
02703 {
02704 m = 0 ;
02705 if ( j == 0 )
02706 {
02707 box_buffer = newVector( edgeclean[0] + margin ) ;
02708 for( col = 0 ; col < edgeclean[0] + margin ; col++ )
02709 {
02710 maxval = -FLT_MAX ;
02711 for ( row = pos_rowclean[0] - nint(y_box) ; row <= pos_rowclean[0] + nint(y_box) ; row++ )
02712 {
02713 if ( maxval < lineImage->data[col + lineImage->lx*row] )
02714 {
02715 maxval = lineImage -> data[col + lineImage->lx*row] ;
02716 }
02717 }
02718 box_buffer->data[m] = maxval ;
02719 m++ ;
02720 }
02721 }
02722 else if ( j < k )
02723 {
02724 box_buffer = newVector( edgeclean[j] - edgeclean[j-1] + 2*margin ) ;
02725 for ( col = edgeclean[j - 1] - margin ; col < edgeclean[j] + margin ; col++ )
02726 {
02727 maxval = -FLT_MAX ;
02728 for ( row = pos_rowclean[j] - nint(y_box) ; row <= pos_rowclean[j] + nint(y_box) ; row++ )
02729 {
02730 if ( maxval < lineImage->data[col + lineImage->lx*row] )
02731 {
02732 maxval = lineImage -> data[col + lineImage->lx*row] ;
02733 }
02734 }
02735 box_buffer->data[m] = maxval ;
02736 m++ ;
02737 }
02738 }
02739 else
02740 {
02741 box_buffer = newVector( lineImage -> lx - edgeclean[k-1] + margin ) ;
02742 for ( col = edgeclean[k - 1] - margin ; col < lineImage -> lx ; col++ )
02743 {
02744 maxval = -FLT_MAX ;
02745 for ( row = pos_rowclean[k] - nint(y_box) ; row <= pos_rowclean[k] + nint(y_box) ; row++ )
02746 {
02747 if ( maxval < lineImage->data[col + lineImage->lx*row] )
02748 {
02749 maxval = lineImage -> data[col + lineImage->lx*row] ;
02750 }
02751 }
02752 box_buffer->data[m] = maxval ;
02753 m++ ;
02754 }
02755 }
02756
02757 for ( left_right = 0 ; left_right <= 1 ; left_right++ )
02758 {
02759 nel = 0 ;
02760 if ( left_right == 0 )
02761 {
02762 nel = box_buffer -> n_elements / 2 ;
02763 }
02764 else
02765 {
02766 if ( box_buffer -> n_elements % 2 == 0 )
02767 {
02768 nel = box_buffer -> n_elements / 2 ;
02769 }
02770 else
02771 {
02772 nel = box_buffer -> n_elements / 2 + 1 ;
02773 }
02774 }
02775
02776
02777 half_buffer = newVector( nel ) ;
02778 if ( left_right == 0 )
02779 {
02780 for ( i = 0 ; i < nel ; i++ )
02781 {
02782 half_buffer -> data[i] = box_buffer -> data[i] ;
02783 }
02784 }
02785 else
02786 {
02787 n_right = 0 ;
02788 for ( i = box_buffer -> n_elements - 1 ; i >= box_buffer -> n_elements - nel ; i-- )
02789 {
02790 half_buffer -> data[n_right] = box_buffer -> data[i] ;
02791 n_right++ ;
02792 }
02793 }
02794
02795 xdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
02796 wdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
02797 mpar = (int *) cpl_calloc( NPAR, sizeof (int) ) ;
02798
02799
02800 minval = FLT_MAX ;
02801 maxval = -FLT_MAX ;
02802 for ( i = 0 ; i < nel ; i++ )
02803 {
02804 xdat[i] = i ;
02805 wdat[i] = 1.0 ;
02806 if ( half_buffer -> data[i] < minval )
02807 {
02808 minval = half_buffer -> data[i] ;
02809 }
02810 if ( half_buffer -> data[i] > maxval )
02811 {
02812 maxval = half_buffer -> data[i] ;
02813 }
02814 }
02815
02816 fitpar[2] = minval ;
02817 fitpar[3] = maxval ;
02818
02819
02820 init1 = -1 ;
02821 for ( i = 0 ; i < nel ; i++ )
02822 {
02823 if ( half_buffer -> data[i] >= ( maxval + minval ) / 2. )
02824 {
02825 init1 = i ;
02826 break ;
02827 }
02828 }
02829
02830
02831 if ( init1 != -1 )
02832 {
02833 fitpar[0] = ((float)init1 - 1.) ;
02834 fitpar[1] = ((float)init1 + 1.) ;
02835 }
02836
02837 for ( i = 0 ; i < NPAR ; i++ )
02838 {
02839 mpar[i] = 1 ;
02840 dervpar[i] = 0. ;
02841 }
02842
02843 xdim = XDIMA ;
02844 ndat = nel ;
02845 numpar = NPAR ;
02846 tol = TOLA ;
02847 lab = LABA ;
02848 its = ITSA ;
02849
02850
02851 if ( 0 > ( iters = lsqfit_edge( xdat, &xdim, half_buffer -> data, wdat, &ndat, fitpar,
02852 dervpar, mpar, &numpar, &tol, &its, &lab )) )
02853 {
02854
02855 cpl_msg_warning ("lsqfit:"," least squares fit failed, error no.: %d in slitlet: %d\n", iters, j) ;
02856 fitpar[0] = ((float)init1 - 1.) ;
02857 fitpar[1] = ((float)init1 + 1.) ;
02858 }
02859
02860 pos = (fitpar[0] + fitpar[1]) / 2. ;
02861
02862
02863
02864
02865
02866
02867
02868 if ( left_right == 0 )
02869 {
02870
02871 if ( j == 0 )
02872 {
02873 if ( fabs(pos - ((float)edgeclean[0] - 1. - (float)slit_length)) < diff_tol )
02874 {
02875 slit_pos[0][0] = pos ;
02876 }
02877 else
02878 {
02879 cpl_msg_warning("lsqfit:","something wrong with fitted left position of slitlet 0\n") ;
02880 if ( (float) edgeclean[0] - 1. - (float)slit_length < 0. )
02881 {
02882 slit_pos[0][0] = 0. ;
02883 }
02884 else
02885 {
02886 slit_pos[0][0] = (float)edgeclean[0] - 1. - (float)slit_length ;
02887 }
02888 }
02889 }
02890 else if ( j < k )
02891 {
02892 if ( fabs( pos - (float)margin ) < diff_tol )
02893 {
02894 slit_pos[j][0] = pos + (float)edgeclean[j-1] - (float)margin ;
02895 }
02896 else
02897 {
02898 cpl_msg_warning("lsqfit:","something wrong with fitted left position of slitlet %d\n", j) ;
02899 slit_pos[j][0] = (float)edgeclean[j-1] - 1. ;
02900 }
02901 }
02902 else
02903 {
02904 if ( fabs( pos - (float)margin ) < diff_tol )
02905 {
02906 slit_pos[k][0] = pos + (float)edgeclean[k-1] - (float)margin ;
02907 }
02908 else
02909 {
02910 cpl_msg_warning("lsqfit:","something wrong with fitted left position of slitlet %d\n", j) ;
02911 slit_pos[k][0] = (float)edgeclean[k-1] - 1. ;
02912 }
02913 }
02914 }
02915 else
02916 {
02917
02918 if ( j == 0 )
02919 {
02920 if ( fabs( (float)box_buffer->n_elements - pos - (float)edgeclean[0] ) < diff_tol )
02921 {
02922 slit_pos[0][1] = (float)(box_buffer->n_elements - 1) - pos ;
02923 }
02924 else
02925 {
02926 cpl_msg_warning("lsqfit:","something wrong with fitted right position of slitlet 0\n") ;
02927 slit_pos[0][1] = (float)edgeclean[0] - 1. ;
02928 }
02929 }
02930 else if ( j < k )
02931 {
02932 if ( fabs( (float)box_buffer->n_elements - pos
02933 + (float)edgeclean[j-1] - (float)margin - (float)edgeclean[j] ) < diff_tol )
02934 {
02935 slit_pos[j][1] = (float)(box_buffer->n_elements - 1) - pos
02936 + (float)edgeclean[j-1] - (float)margin ;
02937 }
02938 else
02939 {
02940 cpl_msg_warning("lsqfit:","something wrong with fitted right position of slitlet %d\n", j) ;
02941 slit_pos[j][1] = (float)edgeclean[j] - 1. ;
02942 }
02943 }
02944 else
02945 {
02946 if ( edgeclean[k-1] + slit_length > lineImage -> lx )
02947 {
02948 last_pos = (float)(lineImage -> lx - 1) ;
02949 }
02950 else
02951 {
02952 last_pos = (float)(edgeclean[k-1] - 1 + slit_length) ;
02953 }
02954 if ( fabs( (float)(box_buffer->n_elements - 1) - pos
02955 + (float)edgeclean[k-1] - (float)margin - last_pos ) < diff_tol )
02956 {
02957 slit_pos[k][1] = (float)(box_buffer->n_elements - 1) - pos
02958 + (float)edgeclean[k-1] - (float)margin ;
02959 }
02960 else
02961 {
02962 cpl_msg_warning("lsqfit:","something wrong with fitted right position of slitlet %d\n", j) ;
02963 slit_pos[k][1] = last_pos ;
02964 }
02965 }
02966 }
02967
02968 destroyVector ( half_buffer ) ;
02969 cpl_free( xdat ) ;
02970 cpl_free( wdat ) ;
02971 cpl_free( mpar ) ;
02972 }
02973 destroyVector ( box_buffer ) ;
02974 }
02975
02976 cpl_free( edge ) ;
02977 cpl_free( pos_row ) ;
02978 cpl_free( edgeclean ) ;
02979 cpl_free( dummyedge ) ;
02980 cpl_free( pos_rowclean ) ;
02981 return 0 ;
02982 }
02983
02984
02985
02986
02987
02988
02989
02990
02991
02992
02993
02994
02995
02996
02997
02998
02999
03000
03001
03002
03003
03004
03005
03006 int fitSlitsEdgeWithEstimate ( OneImage * lineImage,
03007 float ** slit_pos,
03008 int box_length,
03009 float y_box,
03010 float diff_tol,
03011 int low_pos,
03012 int high_pos )
03013 {
03014 int position[lineImage->lx] ;
03015 Vector * box_buffer ;
03016 Vector * in_buffer ;
03017 int found_row ;
03018 int row, col ;
03019 int col_first, col_last ;
03020 int row_first, row_last ;
03021 int i, j, m, n ;
03022 int init1 ;
03023 int left_right ;
03024 int n_buf, shift ;
03025 int slit_length ;
03026 int iters, xdim, ndat ;
03027 int numpar, its ;
03028 int * mpar ;
03029 float * xdat, * wdat ;
03030 float tol, lab ;
03031 float fitpar[NPAR] ;
03032 float dervpar[NPAR] ;
03033 float minval, maxval ;
03034 float pos ;
03035 float new_pos ;
03036 int slitposition[SLITLENGTH] ;
03037 pixelvalue rowpos[SLITLENGTH] ;
03038
03039 slit_length = SLITLENGTH ;
03040 slit_length = N_SLITLETS ;
03041
03042
03043 if ( NullImage == lineImage )
03044 {
03045 cpl_msg_error( "fitSlitsEdgeWithEstimate:"," no line image given!\n" ) ;
03046 return -1 ;
03047 }
03048
03049 if ( NULL == slit_pos )
03050 {
03051 cpl_msg_error( "fitSlitsEdgeWithEstimate:"," no position array allocated!\n" ) ;
03052 return -1 ;
03053 }
03054
03055 if ( box_length < 4 ||
03056 box_length > 2*slit_length )
03057 {
03058 cpl_msg_error( "fitSlitsEdgeWithEstimate:"," wrong fitting box length given!\n" ) ;
03059 cpl_msg_error( "fitSlitsEdgeWithEstimate:"," Must be 4 <= box_length < %d \n",
03060 2*slit_length ) ;
03061 cpl_msg_error( "fitSlitsEdgeWithEstimate:"," You have chosen box_length = %d \n",
03062 box_length);
03063
03064
03065 return -1 ;
03066 }
03067
03068 if ( y_box <= 0. || y_box > 6. )
03069 {
03070 cpl_msg_error( "fitSlitsEdgeWithEstimate:"," wrong y box length given!\n" ) ;
03071 cpl_msg_error( "fitSlitsEdgeWithEstimate:"," You have chosen y_box=%f not in range (0,6]!\n",y_box ) ;
03072 return -1 ;
03073 }
03074 if ( diff_tol <= 0. )
03075 {
03076 cpl_msg_error( "fitSlitsEdgeWithEstimate:"," wrong diff_tol given!\n" ) ;
03077 return -1 ;
03078 }
03079
03080 if ( low_pos >= high_pos || low_pos < 0 || high_pos <= 0 || high_pos > lineImage->ly )
03081 {
03082 cpl_msg_error( "fitSlitsEdgeWithEstimate:"," wrong user given search positions!\n" ) ;
03083 return -1 ;
03084 }
03085
03086
03087 for ( col = 0 ; col < lineImage -> lx ; col++ )
03088 {
03089 found_row = -1 ;
03090 maxval = -FLT_MAX ;
03091 for ( row = low_pos ; row <= high_pos ; row++ )
03092 {
03093 if ( maxval < lineImage->data[col+row*lineImage->lx] )
03094 {
03095 maxval = lineImage->data[col+row*lineImage->lx] ;
03096 found_row = row ;
03097 }
03098 }
03099 if ( maxval > -FLT_MAX && found_row > low_pos )
03100 {
03101 position[col] = found_row ;
03102 }
03103 else
03104 {
03105 position[col] = 0 ;
03106 }
03107 }
03108
03109
03110
03111
03112
03113 for ( j = 0 ; j < slit_length ; j++ )
03114 {
03115
03116
03117
03118 n = 0 ;
03119 for ( col = nint(slit_pos[j][0])+ 1 ; col < nint(slit_pos[j][1]) -1 ; col++ )
03120 {
03121 rowpos[n] = (pixelvalue)position[col] ;
03122 n++ ;
03123 }
03124 slitposition[j] = (int)median(rowpos, n) ;
03125 for ( left_right = 0 ; left_right <= 1 ; left_right++ )
03126
03127 {
03128 init1 = 0 ;
03129 col_first = nint( slit_pos[j][left_right] ) - box_length/2 ;
03130 col_last = nint( slit_pos[j][left_right] ) + box_length/2 ;
03131 if ( col_first < 0 )
03132 {
03133 col_first = 0 ;
03134 init1 = 1 ;
03135 }
03136 if ( col_last > lineImage->lx )
03137 {
03138 col_last = lineImage->lx ;
03139 init1 = 1 ;
03140 }
03141 if ( col_last - col_first <= 0 )
03142 {
03143 cpl_msg_error( "fitSlitsEdgeWithEstimate:"," first and last column wrong!\n" ) ;
03144 return -1 ;
03145 }
03146 box_buffer = newVector( col_last - col_first ) ;
03147 m = 0 ;
03148 if ( left_right == 0 )
03149 {
03150 for( col = col_first ; col < col_last ; col++ )
03151 {
03152 row_first = slitposition[j] - nint(y_box) ;
03153 row_last = slitposition[j] + nint(y_box) ;
03154 if ( row_first < 0 )
03155 {
03156 row_first = 0 ;
03157 }
03158 if ( row_last >= lineImage->ly )
03159 {
03160 row_last = lineImage->ly - 1 ;
03161 }
03162 maxval = -FLT_MAX ;
03163 for ( row = row_first ; row <= row_last ; row++ )
03164 {
03165 if ( maxval < lineImage->data[col + lineImage->lx*row] )
03166 {
03167 maxval = lineImage -> data[col + lineImage->lx*row] ;
03168 }
03169 }
03170 box_buffer->data[m] = maxval ;
03171 m++ ;
03172 }
03173 }
03174 else
03175 {
03176 for( col = col_last-1 ; col >= col_first ; col-- )
03177 {
03178 row_first = slitposition[j] - nint(y_box) ;
03179 row_last = slitposition[j] + nint(y_box) ;
03180 if ( row_first < 0 )
03181 {
03182 row_first = 0 ;
03183 }
03184 if ( row_last >= lineImage->ly )
03185 {
03186 row_last = lineImage->ly - 1 ;
03187 }
03188 maxval = -FLT_MAX ;
03189 for ( row = row_first ; row <= row_last ; row++ )
03190 {
03191 if ( maxval < lineImage->data[col + lineImage->lx*row] )
03192 {
03193 maxval = lineImage -> data[col + lineImage->lx*row] ;
03194 }
03195 }
03196 box_buffer->data[m] = maxval ;
03197 m++ ;
03198 }
03199 }
03200
03201 xdat = (float *) cpl_calloc( box_buffer->n_elements, sizeof (float) ) ;
03202 wdat = (float *) cpl_calloc( box_buffer->n_elements, sizeof (float) ) ;
03203 mpar = (int *) cpl_calloc( NPAR, sizeof (int) ) ;
03204
03205
03206 minval = FLT_MAX ;
03207 maxval = -FLT_MAX ;
03208 for ( i = 0 ; i < box_buffer->n_elements ; i++ )
03209 {
03210 xdat[i] = i ;
03211 wdat[i] = 1.0 ;
03212 if ( box_buffer -> data[i] < minval )
03213 {
03214 minval = box_buffer -> data[i] ;
03215 }
03216 if ( box_buffer -> data[i] > maxval )
03217 {
03218 maxval = box_buffer -> data[i] ;
03219 }
03220 }
03221 fitpar[2] = minval ;
03222 fitpar[3] = maxval ;
03223
03224
03225
03226
03227
03228
03229 if ( init1 == 1 )
03230 {
03231 n_buf = box_buffer->n_elements + box_length/2 ;
03232 in_buffer = newVector( n_buf ) ;
03233 for ( i = 0 ; i < box_length/2 ; i++ )
03234 {
03235 in_buffer -> data[i] = minval ;
03236 }
03237 shift = 0 ;
03238 for ( i = box_length/2 ; i < n_buf ; i++ )
03239 {
03240 in_buffer -> data[i] = box_buffer -> data[shift] ;
03241 shift++ ;
03242 }
03243 destroyVector ( box_buffer ) ;
03244 box_buffer = newVector ( n_buf ) ;
03245 for ( i = 0 ; i < n_buf ; i++ )
03246 {
03247 box_buffer -> data[i] = in_buffer -> data[i] ;
03248 }
03249 destroyVector ( in_buffer ) ;
03250 }
03251
03252 fitpar[0] = (float)box_buffer->n_elements/2. - 1. ;
03253 fitpar[1] = (float)box_buffer->n_elements/2. + 1. ;
03254
03255 for ( i = 0 ; i < NPAR ; i++ )
03256 {
03257 mpar[i] = 1 ;
03258 dervpar[i] = 0. ;
03259 }
03260
03261 xdim = XDIMA ;
03262 ndat = box_buffer->n_elements ;
03263 numpar = NPAR ;
03264 tol = TOLA ;
03265 lab = LABA ;
03266 its = ITSA ;
03267
03268
03269 if ( 0 > ( iters = lsqfit_edge( xdat, &xdim, box_buffer -> data, wdat, &ndat, fitpar,
03270 dervpar, mpar, &numpar, &tol, &its, &lab )) )
03271 {
03272 cpl_msg_warning ("lsqfit_edge:"," least squares fit failed, error no.: %d in slitlet: %d\n", iters, j) ;
03273 destroyVector(box_buffer) ;
03274 cpl_free( xdat ) ;
03275 cpl_free( wdat ) ;
03276 cpl_free( mpar ) ;
03277 continue ;
03278 }
03279 if ( fitpar[1] <= fitpar[0] )
03280 {
03281 cpl_msg_warning ("fitSlitsEdgeWithEstimate:"," fit failed due to negative slope of edge function in slitlet: %d\n", j) ;
03282 destroyVector(box_buffer) ;
03283 cpl_free( xdat ) ;
03284 cpl_free( wdat ) ;
03285 cpl_free( mpar ) ;
03286 continue ;
03287 }
03288
03289 pos = (fitpar[0] + fitpar[1])/2. ;
03290 if ( init1 == 1 )
03291 {
03292 pos -= (float)box_length/2. ;
03293 }
03294
03295
03296
03297
03298
03299
03300 if ( pos != 0. )
03301 {
03302 if ( left_right == 0 )
03303 {
03304 new_pos = (float)col_first + pos ;
03305 }
03306 else
03307 {
03308 new_pos = (float)col_last-1 - pos ;
03309 }
03310 if ( fabs(new_pos - slit_pos[j][left_right]) < diff_tol )
03311 {
03312 slit_pos[j][left_right] = new_pos ;
03313 }
03314 else
03315 {
03316 cpl_msg_warning ("fitSlitsEdgeWithEstimate:"," deviation bigger than tolerance, take the estimated slitlet position in slitlet: %d\n", j) ;
03317 }
03318 }
03319
03320 cpl_free( xdat ) ;
03321 cpl_free( wdat ) ;
03322 cpl_free( mpar ) ;
03323 destroyVector ( box_buffer ) ;
03324 }
03325 }
03326
03327 return 0 ;
03328 }
03329
03330