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