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
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218 #ifdef HAVE_CONFIG_H
00219 # include <config.h>
00220 #endif
00221 #include "sinfo_vltPort.h"
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231 #include "sinfo_function_1d.h"
00232 #include "sinfo_wave_calibration.h"
00233 #include "sinfo_solve_poly_root.h"
00234 #include "sinfo_recipes.h"
00235 #include "sinfo_svd.h"
00236
00237
00253 FitParams ** sinfo_new_fit_params( int n_params )
00254 {
00255 FitParams ** new_params =NULL;
00256 FitParams * temp_params =NULL;
00257 float * temp_fit_mem =NULL;
00258 float * temp_derv_mem=NULL;
00259 int i ;
00260
00261
00262 if ( n_params <= 0 )
00263 {
00264 sinfo_msg_error (" wrong number of lines to fit\n") ;
00265 return NULL ;
00266 }
00267
00268 if (NULL==(new_params=(FitParams **) cpl_calloc ( n_params ,
00269 sizeof (FitParams*) ) ) )
00270 {
00271 sinfo_msg_error (" could not allocate memory\n") ;
00272 return NULL ;
00273 }
00274 if ( NULL == (temp_params = cpl_calloc ( n_params , sizeof (FitParams) ) ) )
00275 {
00276 sinfo_msg_error (" could not allocate memory\n") ;
00277 return NULL ;
00278 }
00279 if ( NULL == (temp_fit_mem = (float *) cpl_calloc( n_params*MAXPAR,
00280 sizeof (float) ) ) )
00281 {
00282 sinfo_msg_error (" could not allocate memory\n") ;
00283 return NULL ;
00284 }
00285
00286
00287 if ( NULL == (temp_derv_mem =
00288 (float *) cpl_calloc( n_params*MAXPAR, sizeof (float) ) ) )
00289 {
00290 sinfo_msg_error (" could not allocate memory\n") ;
00291 return NULL ;
00292 }
00293
00294 for ( i = 0 ; i < n_params ; i ++ )
00295 {
00296 new_params[i] = temp_params+i;
00297 new_params[i] -> fit_par = temp_fit_mem+i*MAXPAR;
00298 new_params[i] -> derv_par = temp_derv_mem+i*MAXPAR;
00299 new_params[i] -> column = 0 ;
00300 new_params[i] -> line = 0 ;
00301 new_params[i] -> wavelength = 0. ;
00302 new_params[i] -> n_params = n_params ;
00303 }
00304
00305 return new_params ;
00306 }
00307
00315 void sinfo_new_destroy_fit_params ( FitParams *** params )
00316 {
00317
00318 if ( *params == NULL )
00319 {
00320 return ;
00321 }
00322
00323 cpl_free ( (*params)[0] -> fit_par ) ;
00324 (*params)[0] -> fit_par=NULL;
00325 cpl_free ( (*params)[0] -> derv_par ) ;
00326 (*params)[0] -> derv_par=NULL;
00327 cpl_free ( (*params)[0] ) ;
00328 (*params)[0]=NULL;
00329 cpl_free ( (*params) ) ;
00330 (*params)=NULL;
00331 }
00332
00341 void sinfo_new_dump_fit_params_to_ascii ( FitParams ** params, char * filename )
00342 {
00343 FILE * fp ;
00344 int i ;
00345
00346 if ( NULL == params )
00347 {
00348 sinfo_msg_error (" no fit parameters available!\n") ;
00349 return ;
00350 }
00351
00352 if ( NULL == filename )
00353 {
00354 sinfo_msg_error (" no filename available!\n") ;
00355 return ;
00356 }
00357
00358 if ( NULL == (fp = fopen ( filename, "w" ) ) )
00359 {
00360 sinfo_msg_error(" cannot open %s\n", filename) ;
00361 return ;
00362 }
00363
00364 for ( i = 0 ; i < params[0] -> n_params ; i++ )
00365 {
00366 fprintf(fp, "%d %d %d %f %f %f %f %f %f %f %f %f\n",
00367 params[i]->n_params,
00368 params[i]->column,
00369 params[i]->line,
00370 params[i]->wavelength,
00371 params[i]->fit_par[0],
00372 params[i]->fit_par[1],
00373 params[i]->fit_par[2],
00374 params[i]->fit_par[3],
00375 params[i]->derv_par[0],
00376 params[i]->derv_par[1],
00377 params[i]->derv_par[2],
00378 params[i]->derv_par[3] ) ;
00379 }
00380 fclose(fp) ;
00381 }
00382
00391 void
00392 sinfo_new_dump_ascii_to_fit_params ( FitParams ** params, char * filename )
00393 {
00394 FILE * fp ;
00395 int i ;
00396
00397 if ( NULL == params )
00398 {
00399 sinfo_msg_error (" no fit parameters available!\n") ;
00400 return ;
00401 }
00402
00403 if ( NULL == filename )
00404 {
00405 sinfo_msg_error (" no filename available!\n") ;
00406 return ;
00407 }
00408
00409 if ( NULL == (fp = fopen ( filename, "r" ) ) )
00410 {
00411 sinfo_msg_error(" cannot open %s\n", filename) ;
00412 return ;
00413 }
00414
00415 for ( i = 0 ; i < params[0]->n_params ; i++ )
00416 {
00417 fscanf(fp, "%d %d %d %f %f %f %f %f %f %f %f %f\n",
00418 ¶ms[i]->n_params,
00419 ¶ms[i]->column,
00420 ¶ms[i]->line,
00421 ¶ms[i]->wavelength,
00422 ¶ms[i]->fit_par[0],
00423 ¶ms[i]->fit_par[1],
00424 ¶ms[i]->fit_par[2],
00425 ¶ms[i]->fit_par[3],
00426 ¶ms[i]->derv_par[0],
00427 ¶ms[i]->derv_par[1],
00428 ¶ms[i]->derv_par[2],
00429 ¶ms[i]->derv_par[3] ) ;
00430 }
00431 fclose(fp) ;
00432 }
00433
00472 int sinfo_new_find_lines(cpl_image * lineImage,
00473 float * wave_position,
00474 float * wave_intensity,
00475 int n_lines,
00476 int ** row_clean,
00477 float ** wavelength_clean,
00478 float beginWave,
00479 float dispersion1,
00480 float dispersion2,
00481 float mindiff,
00482 int halfWidth,
00483 int * n_found_lines,
00484 float sigma,
00485 int * sum_lines )
00486 {
00487 int ** row ;
00488 float ** wavelength ;
00489 float buf1, buf2 ;
00490 float meanval ;
00491 float colmedian ;
00492 float * column, * tempcol ;
00493 float * lines ;
00494 float * conv_lines ;
00495 float * wave_buffer ;
00496 float * wave_mem;
00497 int * dummy_row ;
00498 int i, j, k, m ;
00499 int position ;
00500 int gmax, gmin ;
00501 int col ;
00502 int * row_mem;
00503 float sum ;
00504 float angst ;
00505
00506 int lx=0;
00507 int ly=0;
00508 float* pdata=NULL;
00509
00510 if ( NULL == lineImage )
00511 {
00512 sinfo_msg_error (" no image given\n") ;
00513 return -1 ;
00514 }
00515
00516 lx=cpl_image_get_size_x(lineImage);
00517 ly=cpl_image_get_size_y(lineImage);
00518 pdata=cpl_image_get_data_float(lineImage);
00519
00520 if ( n_lines <= 0 || NULL == wave_position )
00521 {
00522 sinfo_msg_error(" no line list given\n") ;
00523 return -1 ;
00524 }
00525 if ( NULL == wave_intensity )
00526 {
00527 sinfo_msg_error(" no intensity list given\n") ;
00528 return -1 ;
00529 }
00530
00531 if ( dispersion1 == 0. )
00532 {
00533 sinfo_msg_error(" wrong dispersion given\n") ;
00534 return -1 ;
00535 }
00536
00537 if ( row_clean == NULL )
00538 {
00539 sinfo_msg_error(" row_clean array is not allocated\n") ;
00540 return -1 ;
00541 }
00542
00543 if ( wavelength_clean == NULL )
00544 {
00545 sinfo_msg_error(" wavelength_clean array is not allocated\n") ;
00546 return -1 ;
00547 }
00548
00549 if ( beginWave <= 0. )
00550 {
00551 sinfo_msg_error (" impossible beginWave given\n") ;
00552 return -1 ;
00553 }
00554 if ( mindiff < -100. )
00555 {
00556 sinfo_msg_error (" wrong mindiff value\n") ;
00557 return -1 ;
00558 }
00559
00560 if ( halfWidth <= 0 )
00561 {
00562 sinfo_msg_error(" wrong half width of fit box given\n") ;
00563 return -1 ;
00564 }
00565
00566 if ( n_found_lines == NULL )
00567 {
00568 sinfo_msg_error(" n_found_lines not allocated\n") ;
00569 return -1 ;
00570 }
00571
00572 if ( sigma <= 0. || sigma >= ly / 2)
00573 {
00574 sinfo_msg_error(" wrong sigma given\n") ;
00575 return -1 ;
00576 }
00577
00578
00579 row = (int**) cpl_calloc( lx, sizeof(int*)) ;
00580 wavelength = (float**) cpl_calloc( lx, sizeof(float*)) ;
00581 row_mem = cpl_calloc( n_lines*lx, sizeof(int) ) ;
00582 wave_mem = cpl_calloc( n_lines*lx, sizeof(float) ) ;
00583 for ( i = 0 ; i <lx ; i++ )
00584 {
00585 row[i] = row_mem + i*n_lines;
00586 wavelength[i] = wave_mem + i*n_lines;
00587 }
00588
00589
00590 if ( wave_position[0] > 10000. )
00591 {
00592
00593 angst = 10000. ;
00594 }
00595 else if ( wave_position[0] > 1000. && wave_position[0] < 10000. )
00596 {
00597
00598 angst = 1000. ;
00599 }
00600 else
00601 {
00602
00603 angst = 1. ;
00604 }
00605
00606
00607
00608
00609
00610
00611 tempcol = (float*) cpl_calloc(ly, sizeof(float)) ;
00612 *sum_lines = 0 ;
00613 buf1 = 0. ;
00614 buf2 = 0. ;
00615
00616 column = (float*) cpl_calloc(ly, sizeof (float)) ;
00617 lines = (float*) cpl_calloc(ly, sizeof (float)) ;
00618 conv_lines = (float*) cpl_calloc(ly, sizeof (float)) ;
00619 wave_buffer = (float*) cpl_calloc(ly, sizeof (float)) ;
00620 dummy_row = (int*) cpl_calloc(n_lines, sizeof(int)) ;
00621
00622 for ( col = 0 ; col < lx ; col++ )
00623 {
00624 n_found_lines[col] = 0 ;
00625 sum = 0. ;
00626 for ( i = 0 ; i < ly ; i++ )
00627 {
00628 if (isnan(pdata[col + i*lx]) )
00629 {
00630 tempcol[i] = 0. ;
00631 continue ;
00632 }
00633
00634 sum = sum + pdata[col + i*lx] ;
00635 tempcol[i] = pdata[col + i*lx];
00636 }
00637 meanval = sum / ly ;
00638
00639 colmedian = sinfo_new_median ( tempcol, ly);
00640
00641 if ( meanval - colmedian < mindiff )
00642 {
00643 sinfo_msg_warning(" sorry, not enough intensity "
00644 "(mean: %f, diff: %f) in image column: "
00645 "%d to correlate emission lines\n",
00646 meanval, meanval - colmedian,col) ;
00647 continue ;
00648 }
00649
00650 for ( i = 0 ; i < ly ; i++ )
00651 {
00652 conv_lines[i]=0;
00653 wave_buffer[i]=0;
00654 }
00655 for ( i = 0 ; i < n_lines ; i++ )
00656 {
00657 dummy_row[i] = 0;
00658 }
00659
00660
00661 for ( i = 0 ; i < ly ; i++ )
00662 {
00663 if ( isnan(pdata[col+i*lx]) )
00664 {
00665 column[i] = 0. ;
00666 }
00667 else
00668 {
00669 column[i] = pdata[col + i*lx] ;
00670 }
00671
00672
00673
00674 lines[i] = (dispersion1 * (float) (i-ly/2) +
00675 dispersion2 * (float) (i-ly/2) *
00676 (float) (i-ly/2) +
00677 beginWave) * angst ;
00678
00679
00680
00681
00682
00683 for ( j = 0 ; j < n_lines ; j ++ )
00684 {
00685 buf1 = 0. ;
00686 buf2 = 0. ;
00687 if ( (wave_position[j] >= (lines[i] -
00688 fabs(dispersion1)/2.*angst)) &&
00689 (wave_position[j] <= (lines[i] +
00690 fabs(dispersion1)/2.*angst)) )
00691 {
00692 buf1 = wave_intensity[j] ;
00693
00694 buf2 = wave_position[j] ;
00695 break ;
00696 }
00697 }
00698 lines[i] = buf1 ;
00699 wave_buffer[i] = buf2 ;
00700
00701
00702
00703
00704
00705 if ( lines[i] != 0. )
00706 {
00707
00708 gmin = sinfo_new_nint((float) i - 2. * sigma) ;
00709 gmax = sinfo_new_nint((float) i + 2. * sigma) ;
00710
00711
00712 if ( gmin < 0 )
00713 {
00714 gmin = 0 ;
00715 }
00716 if ( gmax >= ly )
00717 {
00718 gmax = ly - 1 ;
00719 }
00720 for ( j = gmin ; j <= gmax ; j++ )
00721 {
00722 conv_lines[j] +=
00723 lines[i] * exp( (double)( -0.5*((j - i)*(j - i)))/
00724 (double) (sigma*sigma) ) ;
00725 }
00726 }
00727 }
00728
00729
00730 position = INT32_MAX ;
00731 position = sinfo_new_correlation(column, conv_lines, ly ) ;
00732 if ( abs (position) > ly / 4 )
00733 {
00734 sinfo_msg_warning(" sorry, shift of artificial data relative to"
00735 " image (%d) seems to be too high in column: %d",
00736 position, col) ;
00737 continue ;
00738 }
00739
00740 j = 0 ;
00741 for ( i = 0 ; i < ly ; i ++ )
00742 {
00743 if ( lines[i] != 0.0 )
00744 {
00745 if ( (i - position) >= 0 && (i - position) < ly )
00746 {
00747 row[col][j] = i - position ;
00748
00749
00750 wavelength[col][j] = wave_buffer[i] / angst ;
00751 j++ ;
00752 }
00753 }
00754 }
00755
00756
00757
00758
00759
00760
00761
00762 for ( k = 1 ; k <= j && k<(lx-1); k ++ )
00763 {
00764 if (dummy_row[k-1] != -1)
00765 {
00766 dummy_row[k-1] = row[col][k-1] ;
00767 }
00768 if ( (row[col][k] - row[col][k-1]) <= 2*halfWidth )
00769 {
00770 dummy_row[k-1] = -1 ;
00771 if (k<n_lines)
00772 dummy_row[k] = -1 ;
00773 }
00774
00775 if ( (row[col][k+1] - row[col][k]) <= 2*halfWidth)
00776 {
00777 if (k<n_lines)
00778 dummy_row[k] = -1 ;
00779 if (k+1<n_lines)
00780 dummy_row[k+1] = -1 ;
00781 }
00782 }
00783
00784 m = 0 ;
00785 for ( k = 0 ; k < j ; k ++ )
00786 {
00787 if ( dummy_row[k] != -1 && dummy_row[k] != 0 )
00788 {
00789 row_clean[col][m] = dummy_row[k] ;
00790 wavelength_clean[col][m] = wavelength[col][k] ;
00791 m ++ ;
00792 }
00793 }
00794
00795 n_found_lines[col] = m ;
00796
00797 *sum_lines += n_found_lines[col] ;
00798 }
00799 cpl_free (column) ;
00800 cpl_free (lines) ;
00801 cpl_free (conv_lines) ;
00802 cpl_free (dummy_row) ;
00803 cpl_free (wave_buffer) ;
00804 cpl_free (row_mem) ;
00805 cpl_free (wave_mem) ;
00806 cpl_free (tempcol) ;
00807 cpl_free (row) ;
00808 cpl_free (wavelength) ;
00809
00810 return 0 ;
00811 }
00812
00822 int
00823 sinfo_new_read_list( char * listname,
00824 float * lineCenter,
00825 float * lineIntensity )
00826 {
00827 FILE * fp ;
00828 int i, n_lines ;
00829
00830 if ( NULL == lineCenter )
00831 {
00832 sinfo_msg_error(" lineCenter array is not allocated\n") ;
00833 return -1 ;
00834 }
00835
00836 if ( NULL == lineIntensity )
00837 {
00838 sinfo_msg_error(" lineIntensity array is not allocated\n") ;
00839 return -1 ;
00840 }
00841
00842 if ( NULL == (fp = fopen ( listname, "r" ) ) )
00843 {
00844 sinfo_msg_error(" cannot open %s\n", listname) ;
00845 return -1 ;
00846 }
00847
00848 i = 0 ;
00849 while ( fscanf( fp, "%f %f", &lineCenter[i], &lineIntensity[i] ) != EOF )
00850 {
00851 i ++ ;
00852 }
00853 n_lines = i ;
00854 fclose(fp) ;
00855
00856 return n_lines ;
00857 }
00858
00859
00891 int sinfo_new_line_fit ( cpl_image * mergedImage,
00892 FitParams * par,
00893 float fwhm,
00894 int lineInd,
00895 int column,
00896 int halfWidth,
00897 int lineRow,
00898 float min_amplitude,
00899 Vector * line,
00900 int * mpar,
00901 float * xdat,
00902 float * wdat )
00903 {
00904 int i, j ;
00905 int iters, xdim, ndat ;
00906 int numpar, its ;
00907 int position ;
00908 float maxval, tol, lab ;
00909 int lx=0;
00910 int ly=0;
00911 float* pdata=NULL;
00912
00913 if ( mergedImage == NULL )
00914 {
00915 sinfo_msg_error (" no image given as input\n") ;
00916 return -8 ;
00917 }
00918 lx=cpl_image_get_size_x(mergedImage);
00919 ly=cpl_image_get_size_y(mergedImage);
00920 pdata=cpl_image_get_data_float(mergedImage);
00921
00922
00923 if ( par == NULL )
00924 {
00925 sinfo_msg_error(" fit parameters not given\n") ;
00926 return -9 ;
00927 }
00928 if ( column < 0 || column > lx )
00929 {
00930 sinfo_msg_error (" wrong column number\n") ;
00931 return -10 ;
00932 }
00933 if ( halfWidth < 0 || halfWidth > ly )
00934 {
00935 sinfo_msg_error (" wrong width given\n") ;
00936 return -11 ;
00937 }
00938 if ( lineRow < 0 || lineRow > ly )
00939 {
00940 sinfo_msg_error (" wrong number of row of the line given\n") ;
00941 return -12 ;
00942 }
00943 if ( min_amplitude < 1. )
00944 {
00945 sinfo_msg_error (" wrong minimum amplitude\n") ;
00946 return -13 ;
00947 }
00948
00949
00950 for ( i = 0 ; i < line -> n_elements ; i++)
00951 {
00952 line->data[i] = 0;
00953 }
00954
00955 par -> column = column ;
00956 par -> line = lineInd ;
00957
00958
00959
00960
00961 j = 0 ;
00962 for ( i = lineRow-halfWidth ; i <= lineRow+halfWidth ; i++ )
00963 {
00964 if ( i < 0 || i >= ly )
00965 {
00966 sinfo_msg_error (" wrong line position or width given\n") ;
00967 return -15 ;
00968 }
00969 else
00970 {
00971 line -> data[j] = pdata[column + i*lx] ;
00972 j ++ ;
00973 }
00974 }
00975
00976
00977
00978
00979
00980 maxval = -FLT_MAX ;
00981 position = -INT32_MAX ;
00982 for ( i = 0 ; i < line -> n_elements ; i++ )
00983 {
00984 xdat[i] = i ;
00985 wdat[i] = 1.0 ;
00986 if ( line -> data[i] >= maxval )
00987 {
00988 maxval = line -> data[i] ;
00989 position = i ;
00990 }
00991 }
00992
00993
00994 xdim = XDIM ;
00995 ndat = line -> n_elements ;
00996 numpar = MAXPAR ;
00997 tol = TOL ;
00998 lab = LAB ;
00999 its = ITS ;
01000 par -> fit_par[1] = fwhm ;
01001 par -> fit_par[2] = (float) position ;
01002 par -> fit_par[3] = (float) (line -> data[0] +
01003 line -> data[line->n_elements - 1]) / 2.0 ;
01004 par -> fit_par[0] = maxval - (par -> fit_par[3]) ;
01005
01006
01007 if ( par->fit_par[0] < min_amplitude )
01008 {
01009 cpl_msg_debug ("sinfo_linefit:",
01010 " sorry, amplitude of line too low to fit: %f",
01011 par->fit_par[0] ) ;
01012 return -16 ;
01013 }
01014
01015 for ( i = 0 ; i < MAXPAR ; i++ )
01016 {
01017 par -> derv_par[i] = 0.0 ;
01018 mpar[i] = 1 ;
01019 }
01020
01021
01022 if ( 0 > ( iters = sinfo_new_lsqfit_c( xdat, &xdim,
01023 line -> data, wdat,
01024 &ndat, par -> fit_par,
01025 par -> derv_par, mpar,
01026 &numpar, &tol, &its, &lab )) )
01027 {
01028 cpl_msg_debug ("sinfo_linefit:",
01029 " sinfo_new_lsqfit_c: least squares fit failed,"
01030 " error no.: %d\n", iters) ;
01031 return -17 ;
01032 }
01033
01034
01035
01036 par -> fit_par[2] = (float) (lineRow - halfWidth) + par -> fit_par[2] ;
01037
01038
01039 return iters ;
01040 }
01041
01066 int sinfo_new_fit_lines ( cpl_image * line_image,
01067 FitParams ** allParams,
01068 float fwhm,
01069 int * n_lines,
01070 int ** row,
01071 float ** wavelength,
01072 int width,
01073 float min_amplitude )
01074 {
01075 int i, k, l ;
01076 int result ;
01077 Vector * line;
01078 int * mpar;
01079 float * xdat, * wdat;
01080 int lx=0;
01081 int ly=0;
01082 float* pdata=NULL;
01083
01084 if ( line_image == NULL )
01085 {
01086 sinfo_msg_error (" no image given\n") ;
01087 return -18 ;
01088 }
01089 lx=cpl_image_get_size_x(line_image);
01090 ly=cpl_image_get_size_y(line_image);
01091 pdata=cpl_image_get_data_float(line_image);
01092
01093 if ( n_lines == NULL )
01094 {
01095 sinfo_msg_error (" no counter of emission lines\n") ;
01096 return -19 ;
01097 }
01098 if ( row == NULL || width <= 0 )
01099 {
01100 sinfo_msg_error (" row or width vectors are empty\n") ;
01101 return -20 ;
01102 }
01103 if ( wavelength == NULL )
01104 {
01105 sinfo_msg_error (" no wavelength array given\n") ;
01106 return -21 ;
01107 }
01108
01109 k = 0 ;
01110
01111
01112 line = sinfo_new_vector (2*width + 1) ;
01113
01114 xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
01115 wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
01116 mpar = (int *) cpl_calloc( MAXPAR, sizeof (int) ) ;
01117
01118
01119 for ( i = 0 ; i < lx ; i++ )
01120 {
01121 if ( n_lines[i] == 0 )
01122 {
01123 continue ;
01124 }
01125
01126 for ( l = 0 ; l < n_lines[i] ; l++ )
01127 {
01128 if ( row[i][l] <= 0 )
01129 {
01130 continue ;
01131 }
01132
01133
01134
01135
01136
01137
01138 if ( (result = sinfo_new_line_fit ( line_image,
01139 allParams[k], fwhm, l, i,
01140 width, row[i][l],
01141 min_amplitude,line,mpar,
01142 xdat,wdat ) ) < 0 )
01143 {
01144 cpl_msg_debug ("sinfo_fitLines:",
01145 " sinfo_linefit failed, error no.: %d,"
01146 " column: %d, row: %d, line: %d\n",
01147 result, i, row[i][l], l) ;
01148 continue ;
01149 }
01150 if ( (allParams[k] -> fit_par[0] <= 0.) ||
01151 (allParams[k] -> fit_par[1] <= 0.)
01152 || (allParams[k] -> fit_par[2] <= 0.) )
01153 {
01154 sinfo_msg_warning (" negative fit parameters in column: %d,"
01155 " line: %d\n", i, l) ;
01156 continue ;
01157 }
01158 allParams[k] -> wavelength = wavelength[i][l] ;
01159 k++ ;
01160 }
01161 }
01162
01163
01164 sinfo_new_destroy_vector(line);
01165 cpl_free(xdat);
01166 cpl_free(wdat);
01167 cpl_free(mpar);
01168
01169
01170 return k ;
01171 }
01172
01200 float sinfo_new_polyfit( FitParams ** par,
01201 int column,
01202 int n_lines,
01203 int n_rows,
01204 float dispersion,
01205 float max_residual,
01206 float * acoefs,
01207 float * dacoefs,
01208 int * n_reject,
01209 int n_fitcoefs )
01210 {
01211 float ** ucoefs, ** vcoefs, ** covar ;
01212 float *mem;
01213 float * lambda, * posit ;
01214 float * weight, * resid ;
01215 float * newlam, * newpos, * newwet ;
01216 float * wcoefs=NULL ;
01217 float chisq, result ;
01218 float offset ;
01219 int num, found ;
01220 int i, j, k, n ;
01221
01222
01223 for ( i = 0 ; i < n_fitcoefs ; i++ )
01224 {
01225 acoefs[i] = 0. ;
01226 dacoefs[i] = 0. ;
01227 }
01228 if ( NULL == par )
01229 {
01230 sinfo_msg_error(" no fit params given\n");
01231 return FLT_MAX ;
01232 }
01233
01234 if ( 0 >= n_lines )
01235 {
01236
01237
01238
01239 return FLT_MAX ;
01240 }
01241 if ( 0 >= n_rows )
01242 {
01243 sinfo_msg_error (" sorry, number of rows is wrong") ;
01244 return FLT_MAX ;
01245 }
01246 if ( dispersion == 0. )
01247 {
01248 sinfo_msg_error (" sorry, wrong dispersion given") ;
01249 return FLT_MAX ;
01250 }
01251
01252 offset = (float)(n_rows - 1)/2. ;
01253
01254
01255
01256 mem = (float*) cpl_calloc( n_lines*7, sizeof (float) ) ;
01257 lambda = mem;
01258 posit = mem + n_lines;
01259 weight = mem + n_lines*2;
01260 resid = mem + n_lines*3;
01261 newlam = mem + n_lines*4;
01262 newpos = mem + n_lines*5;
01263 newwet = mem + n_lines*6;
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274 ucoefs = sinfo_matrix ( 1, n_lines, 1, n_fitcoefs ) ;
01275 vcoefs = sinfo_matrix ( 1, n_lines, 1, n_fitcoefs ) ;
01276 covar = sinfo_matrix ( 1, n_fitcoefs, 1, n_fitcoefs ) ;
01277 wcoefs=cpl_calloc(n_fitcoefs,sizeof(float)) ;
01278
01279
01280 n = 0 ;
01281 for ( i = 0 ; i < (par[0] -> n_params) ; i ++ )
01282 {
01283 found = -1 ;
01284
01285 for ( j = 0 ; j < n_lines ; j ++ )
01286 {
01287 if ( (par[i] -> column == column) && (par[i] -> line == j) )
01288 {
01289 found = i ;
01290 }
01291 else
01292 {
01293 continue ;
01294 }
01295
01296
01297 if ( par[found] -> derv_par[2] != 0. &&
01298 par[found] -> fit_par[2] > 0. &&
01299 par[found] -> wavelength > 0. &&
01300 par[found] -> fit_par[1] > 0. &&
01301 par[found] -> fit_par[0] > 0. )
01302 {
01303
01304
01305
01306
01307
01308 posit[n] = par[found] -> fit_par[2] ;
01309 weight[n] = par[found] -> derv_par[2] ;
01310 lambda[n] = par[found] -> wavelength ;
01311 n ++ ;
01312 }
01313 else
01314 {
01315 continue ;
01316 }
01317 }
01318
01319 }
01320
01321 num = n ;
01322 if ( num < n_fitcoefs )
01323 {
01324 sinfo_msg_warning("not enough lines found in column %d to "
01325 "determine the three coefficients.\n", column) ;
01326 for ( i = 0 ; i < n_fitcoefs ; i++ )
01327 {
01328 acoefs[i] = ZERO ;
01329 dacoefs[i] = ZERO ;
01330 }
01331 sinfo_free_matrix ( ucoefs, 1, 1 ) ;
01332 sinfo_free_matrix ( vcoefs, 1, 1 ) ;
01333 sinfo_free_matrix ( covar, 1, 1 ) ;
01334
01335
01336
01337
01338
01339
01340
01341 cpl_free (mem);
01342 return FLT_MAX ;
01343 }
01344
01345
01346
01347
01348
01349
01350 for ( i = 0 ; i < num ; i ++ )
01351 {
01352 posit[i] = (posit[i] - offset)/offset ;
01353 weight[i] *= fabs(dispersion) ;
01354 }
01355
01356
01357 sinfo_svd_fitting( posit - 1, lambda - 1,
01358 weight - 1, num, acoefs-1, n_fitcoefs,
01359 ucoefs, vcoefs, wcoefs-1, covar, &chisq, sinfo_fpol ) ;
01360
01361
01362 for ( i = 1 ; i < n_fitcoefs ; i++ )
01363 {
01364 acoefs[i] /= pow(offset, i) ;
01365 }
01366
01367
01368 *n_reject = 0 ;
01369
01370 j = 0 ;
01371 for ( i = 0 ; i < num ; i++ )
01372 {
01373 result = 0. ;
01374 for ( k = 0 ; k < n_fitcoefs ; k++ )
01375 {
01376 result += acoefs[k] * pow(posit[i], k) ;
01377 }
01378
01379 resid[i] = lambda[i] - result ;
01380
01381 if ( fabs( resid[i] ) > max_residual)
01382 {
01383 (*n_reject) ++ ;
01384 }
01385 else
01386 {
01387 newlam[j] = lambda[i] ;
01388 newpos[j] = posit[i] ;
01389 newwet[j] = weight[i] ;
01390 j++ ;
01391 }
01392 }
01393
01394 num = j ;
01395 if ( num >= n_fitcoefs )
01396 {
01397 sinfo_svd_fitting( newpos - 1, newlam - 1,
01398 newwet - 1, num, acoefs-1, n_fitcoefs, ucoefs,
01399 vcoefs, wcoefs-1, covar, &chisq, sinfo_fpol ) ;
01400
01401
01402 for ( i = 0 ; i < n_fitcoefs ; i++ )
01403 {
01404 acoefs[i] /= pow(offset, i) ;
01405 dacoefs[i] = sqrt( (double) covar[i+1][i+1] ) / pow(offset, i) ;
01406 }
01407 }
01408 else
01409 {
01410 sinfo_msg_warning (" too many lines rejected (number: %d) "
01411 "due to high residuals, fit coefficients are set "
01412 "zero, in column: %d\n", *n_reject, column) ;
01413 for ( i = 0 ; i < n_fitcoefs ; i++ )
01414 {
01415 acoefs[i] = ZERO ;
01416 dacoefs[i] = ZERO ;
01417 }
01418 }
01419
01420 sinfo_free_matrix ( ucoefs, 1, 1 ) ;
01421 sinfo_free_matrix ( vcoefs, 1, 1 ) ;
01422 sinfo_free_matrix ( covar, 1, 1 ) ;
01423
01424
01425
01426
01427
01428
01429
01430 cpl_free (mem);
01431 cpl_free(wcoefs) ;
01432
01433 return chisq ;
01434 }
01435
01452 float sinfo_new_coefs_cross_fit ( int n_columns,
01453 float * acoefs,
01454 float * dacoefs,
01455 float * bcoefs,
01456 int n_fitcoefs,
01457 float sigma_factor )
01458 {
01459 float col_index;
01460 float* sub_col_index=NULL ;
01461 float* sub_acoefs=NULL ;
01462 float* sub_dacoefs=NULL ;
01463 float* wcoefs=NULL ;
01464 float ** ucoefs, **vcoefs, **covar ;
01465 float chisq ;
01466 float * acoefsclean ;
01467 double sum, sumq, mean ;
01468 double sigma ;
01469 double cliphi, cliplo ;
01470 float offset ;
01471 int i, n, num, ndata ;
01472 int nc ;
01473
01474
01475 if ( n_columns < 1 )
01476 {
01477 sinfo_msg_error(" wrong number of image columns given\n") ;
01478 return FLT_MAX ;
01479 }
01480 if ( acoefs == NULL || dacoefs == NULL )
01481 {
01482 sinfo_msg_error(" coeffs or errors of coefficients are not given\n") ;
01483 return FLT_MAX ;
01484 }
01485 if ( bcoefs == NULL )
01486 {
01487 sinfo_msg_error(" coeffs are not allocated\n") ;
01488 return FLT_MAX ;
01489 }
01490
01491 if ( n_fitcoefs < 1 )
01492 {
01493 sinfo_msg_error(" wrong number of fit coefficients\n") ;
01494 return FLT_MAX ;
01495 }
01496 if ( sigma_factor <= 0. )
01497 {
01498 sinfo_msg_error(" impossible sigma_factor given!\n") ;
01499 return FLT_MAX ;
01500 }
01501
01502 offset = (float)(n_columns - 1) / 2. ;
01503
01504
01505
01506
01507
01508
01509 wcoefs=cpl_calloc(n_fitcoefs,sizeof(float)) ;
01510
01511 nc = 0 ;
01512 for ( i = 0 ; i < n_columns ; i++ )
01513 {
01514 if ( isnan(acoefs[i]) || acoefs[i] == 0. || dacoefs[i] == 0. )
01515 {
01516 continue ;
01517 }
01518 else
01519 {
01520 nc++ ;
01521 }
01522 }
01523 acoefsclean = (float*) cpl_calloc(nc , sizeof(float)) ;
01524 nc = 0 ;
01525 for ( i = 0 ; i < n_columns ; i++ )
01526 {
01527 if ( isnan(acoefs[i]) || acoefs[i] == 0. || dacoefs[i] == 0. )
01528 {
01529 continue ;
01530 }
01531 else
01532 {
01533 acoefsclean[nc] = acoefs[i] ;
01534 nc++ ;
01535 }
01536 }
01537 sinfo_pixel_qsort(acoefsclean, nc) ;
01538 sum = 0. ;
01539 sumq = 0. ;
01540 mean = 0. ;
01541 sigma = 0. ;
01542 n = 0 ;
01543 for ( i = (int)((float)nc*LOW_REJECT) ;
01544 i < (int)((float)nc*HIGH_REJECT) ; i++ )
01545 {
01546 sum += (double)acoefsclean[i] ;
01547 sumq += ((double)acoefsclean[i] * (double)acoefsclean[i]) ;
01548 n ++ ;
01549 }
01550 mean = sum/(double)n ;
01551 sigma = sqrt( sumq/(double)n - (mean * mean) ) ;
01552 cliphi = mean + sigma * (double)sigma_factor ;
01553 cliplo = mean - sigma * (double)sigma_factor ;
01554
01555 sub_col_index=cpl_calloc(n_columns,sizeof(float)) ;
01556 sub_acoefs=cpl_calloc(n_columns,sizeof(float));
01557 sub_dacoefs=cpl_calloc(n_columns,sizeof(float)) ;
01558
01559
01560 num = 0 ;
01561 for ( i = 0 ; i < n_columns ; i++ )
01562 {
01563
01564 col_index = (float) i ;
01565
01566
01567 if ( !isnan(acoefs[i]) &&
01568 (acoefs[i] <= cliphi) && (acoefs[i] >= cliplo) &&
01569 (dacoefs[i] != 0. ) && (acoefs[i] != 0.) )
01570 {
01571 sub_acoefs[num] = acoefs[i] ;
01572 sub_dacoefs[num] = dacoefs[i] ;
01573 sub_col_index[num] = col_index ;
01574 num ++ ;
01575 }
01576 }
01577 ndata = num ;
01578
01579 if ( ndata < n_fitcoefs )
01580 {
01581 sinfo_msg_error("not enough data found to determine "
01582 "the fit coefficients.\n") ;
01583
01584 return FLT_MAX ;
01585 }
01586
01587
01588 ucoefs = sinfo_matrix(1, ndata, 1, n_fitcoefs) ;
01589 vcoefs = sinfo_matrix(1, ndata, 1, n_fitcoefs) ;
01590 covar = sinfo_matrix ( 1, n_fitcoefs, 1, n_fitcoefs ) ;
01591
01592
01593 for ( i = 0 ; i < ndata ; i++ )
01594 {
01595 sub_col_index[i] = (sub_col_index[i] - offset) / offset ;
01596 }
01597
01598
01599 sinfo_svd_fitting ( sub_col_index-1, sub_acoefs-1,
01600 sub_dacoefs-1, ndata, bcoefs-1,
01601 n_fitcoefs, ucoefs, vcoefs,
01602 wcoefs-1, covar, &chisq, sinfo_fpol ) ;
01603
01604
01605 for ( i = 0 ; i < n_fitcoefs ; i ++ )
01606 {
01607 bcoefs[i] /= pow(offset, i) ;
01608 }
01609
01610
01611 cpl_free (acoefsclean) ;
01612 sinfo_free_matrix( ucoefs, 1, 1) ;
01613 sinfo_free_matrix( vcoefs, 1, 1) ;
01614 sinfo_free_matrix ( covar, 1, 1 ) ;
01615
01616 cpl_free(sub_col_index) ;
01617 cpl_free(sub_acoefs) ;
01618 cpl_free(sub_dacoefs) ;
01619 cpl_free(wcoefs) ;
01620
01621 return chisq ;
01622 }
01623
01624
01644 cpl_image * sinfo_new_wave_map( cpl_image * lineImage,
01645 float ** bcoefs,
01646 int n_a_fitcoefs,
01647 int n_b_fitcoefs,
01648 float * wavelength,
01649 float * intensity,
01650 int n_lines,
01651 int magFactor)
01652 {
01653 cpl_image * retImage ;
01654 float cenpos, cenpix ;
01655 float centreval, centrepix, wavelag ;
01656 float pixvalue ;
01657 float a_initial ;
01658 int i, j, k, l, line, col, row, found, sign ;
01659 int var, maxlag, cmin, cmax, offset ;
01660 double * result ;
01661 float col_off ;
01662 float angst ;
01663 double xcorr_max ;
01664 int delta ;
01665
01666 double* z=NULL ;
01667 double* a=NULL ;
01668 double* wave=NULL ;
01669 float* emline=NULL ;
01670 float* spec=NULL ;
01671 int ilx=0;
01672 int ily=0;
01673 int olx=0;
01674 int oly=0;
01675 float* pidata=NULL;
01676 float* podata=NULL;
01677
01678
01679 gsl_poly_complex_workspace * w ;
01680
01681 if ( NULL == lineImage )
01682 {
01683 sinfo_msg_error("no image given\n") ;
01684 return NULL ;
01685 }
01686 ilx=cpl_image_get_size_x(lineImage);
01687 ily=cpl_image_get_size_y(lineImage);
01688 pidata=cpl_image_get_data_float(lineImage);
01689
01690 if ( NULL == wavelength || n_lines <= 0 )
01691 {
01692 sinfo_msg_error("no wavelength list given\n") ;
01693 return NULL ;
01694 }
01695
01696 if ( NULL == intensity )
01697 {
01698 sinfo_msg_error("no intensity list given\n") ;
01699 return NULL ;
01700 }
01701
01702 if ( NULL == bcoefs )
01703 {
01704 sinfo_msg_error("no coefficients given\n") ;
01705 return NULL ;
01706 }
01707
01708 if ( magFactor <= 1 )
01709 {
01710 sinfo_msg_error("wrong magnifying factor given\n") ;
01711 return NULL ;
01712 }
01713
01714
01715 if ( NULL == ( retImage = cpl_image_new( ilx, ily,CPL_TYPE_FLOAT ) ))
01716 {
01717 sinfo_msg_error("cannot allocate a new image\n");
01718 return NULL ;
01719 }
01720 olx=cpl_image_get_size_x(retImage);
01721 oly=cpl_image_get_size_y(retImage);
01722 podata=cpl_image_get_data_float(retImage);
01723
01724
01725 var = (magFactor - 1)*(magFactor - 1) ;
01726 offset = ily * (magFactor/4 + 1) ;
01727
01728
01729 if ( wavelength[0] > 10000. )
01730 {
01731
01732 angst = 10000. ;
01733 }
01734 else if ( wavelength[0] > 1000. && wavelength[0] < 10000. )
01735 {
01736
01737 angst = 1000. ;
01738 }
01739 else
01740 {
01741
01742 angst = 1. ;
01743 }
01744
01745 z=cpl_calloc(2*(n_a_fitcoefs - 1),sizeof(double)) ;
01746 a=cpl_calloc(n_a_fitcoefs,sizeof(double));
01747 wave=cpl_calloc(n_lines,sizeof(double)) ;
01748 emline=cpl_calloc(2*magFactor*ily,sizeof(float));
01749 spec=cpl_calloc(2*magFactor*ily,sizeof(float)) ;
01750
01751
01752 for ( col = 0 ; col < ilx ; col++ )
01753 {
01754
01755 for ( i = 0 ; i < 2*magFactor*ily ; i++ )
01756 {
01757 emline[i] = 0. ;
01758 }
01759 col_off = (float)col - (float)(ilx-1)/2. ;
01760
01761
01762 for ( i = 0 ; i < n_a_fitcoefs ; i++ )
01763 {
01764
01765 a[i] = 0. ;
01766 if (i < n_a_fitcoefs-1)
01767 {
01768 z[2*i] = 0. ;
01769 z[2*i+1] = 0. ;
01770 }
01771 for ( j = 0 ; j < n_b_fitcoefs ; j++ )
01772 {
01773 a[i] += bcoefs[i][j] * pow(col_off, j) ;
01774 }
01775 }
01776 a_initial = a[0] ;
01777
01778
01779 for ( line = 0 ; line < n_lines ; line++ )
01780 {
01781
01782 wave[line] = wavelength[line]/angst ;
01783
01784
01785
01786
01787
01788 a[0] = a_initial - wave[line] ;
01789
01790 if(NULL==(w=sinfo_gsl_poly_complex_workspace_alloc(n_a_fitcoefs)))
01791 {
01792 sinfo_msg_error("could not allocate complex workspace!") ;
01793 cpl_image_delete(retImage) ;
01794 return NULL ;
01795 }
01796 if (-1 == sinfo_gsl_poly_complex_solve(a, n_a_fitcoefs, w, z))
01797 {
01798 sinfo_msg_error("sinfo_gsl_poly_complex_solve did not work!") ;
01799 cpl_image_delete(retImage) ;
01800 return NULL ;
01801 }
01802 sinfo_gsl_poly_complex_workspace_free(w) ;
01803
01804
01805 j = 0 ;
01806 found = -1 ;
01807 for ( i = 0 ; i < n_a_fitcoefs - 1 ; i++ )
01808 {
01809
01810 if( z[2*i] > (-1.)*(float) ily/2. &&
01811 z[2*i] < (float)ily/2. && z[2*i+1] == 0. )
01812 {
01813 found = 2*i ;
01814 j ++ ;
01815 }
01816 else
01817 {
01818 continue ;
01819 }
01820 }
01821 if ( j == 0 )
01822 {
01823 sinfo_msg_warning("no offset solution found "
01824 "for line %d in column %d\n", line, col) ;
01825 continue ;
01826 }
01827 else if ( j == 1 )
01828 {
01829 cenpos = z[found] + (float) ily /2. ;
01830 }
01831 else
01832 {
01833 sinfo_msg_warning("two or more offset solutions found "
01834 "for line %d in column %d\n", line, col) ;
01835 continue ;
01836 }
01837
01838
01839
01840
01841 cenpix = cenpos * (float) magFactor + (float) offset ;
01842
01843
01844
01845 cmin = (sinfo_new_nint(cenpix) - (var-1)) > 0 ?
01846 sinfo_new_nint(cenpix) - (var-1) : 0 ;
01847 cmax = (sinfo_new_nint(cenpix) + (var-1)) < 2*magFactor * ily ?
01848 sinfo_new_nint(cenpix) + (var-1) : 2*magFactor * ily ;
01849
01850
01851 for ( j = cmin ; j < cmax ; j++ )
01852 {
01853 emline[j] += intensity[line] *
01854 exp((double)(-0.5*(j-cenpix)*(j-cenpix))/(double)var) ;
01855 }
01856 }
01857
01858
01859
01860
01861
01862
01863 for ( k = 0 ; k < 2*magFactor * ily ; k++ )
01864 {
01865 spec[k] = 0. ;
01866 }
01867
01868
01869
01870 for ( row = 0 ; row < ily ; row++ )
01871 {
01872
01873
01874 for ( l = 0 ; l < magFactor ; l++ )
01875 {
01876
01877 if (!isnan(pidata[col + row * ilx]) &&
01878 (pidata[col + row * ilx] > 0.))
01879 {
01880 spec[offset + l + (row * magFactor)] =
01881 pidata[col + row * ilx] ;
01882 }
01883 else
01884 {
01885 spec[offset + l + (row * magFactor)] = 0. ;
01886 }
01887 }
01888 }
01889
01890
01891 if (NULL == (result = sinfo_new_xcorrel(spec, 2*magFactor * ily,
01892 emline, 2*magFactor * ily,
01893 magFactor * ily, &delta,
01894 &maxlag, &xcorr_max)) )
01895 {
01896 sinfo_msg_warning ("cross sinfo_new_correlation did not work,"
01897 " col: %d is set ZERO\n", col) ;
01898 for ( row = 0 ; row < ily ; row++ )
01899 {
01900 podata[col + row * ilx] = ZERO ;
01901 }
01902 continue ;
01903 }
01904
01905 if ( xcorr_max <= 0. )
01906 {
01907 sinfo_msg_warning ("cross sinfo_new_correlation sum is negative,"
01908 " col: %d is set ZERO\n", col) ;
01909 for ( row = 0 ; row < ily ; row++ )
01910 {
01911 podata[col + row * ilx] = ZERO ;
01912 }
01913 cpl_free(result) ;
01914 continue ;
01915 }
01916
01917 wavelag = (float) -delta / (float) magFactor ;
01918 if ( fabs(wavelag) > (float)ily/20. )
01919 {
01920 sinfo_msg_warning("wave lag too big, col: %d is set ZERO\n", col) ;
01921 for ( row = 0 ; row < ily ; row++ )
01922 {
01923 podata[col + row * ilx] = ZERO ;
01924 }
01925 cpl_free(result) ;
01926 continue ;
01927 }
01928
01929
01930
01931
01932
01933
01934
01935
01936 centreval = a_initial ;
01937 for ( i = 1 ; i < n_a_fitcoefs ; i++ )
01938 {
01939 if ( i%2 == 0 )
01940 {
01941 sign = -1 ;
01942 }
01943 else
01944 {
01945 sign = 1 ;
01946 }
01947 centreval += (float)sign * a[i]*pow(wavelag, i) ;
01948 }
01949
01950
01951 for ( row = 0 ; row < ily ; row++ )
01952 {
01953 centrepix = (float)row - ((float)ily - 1.)/2. ;
01954 pixvalue = 0. ;
01955 for ( i = 1 ; i < n_a_fitcoefs ; i++ )
01956 {
01957 pixvalue += a[i]*pow(centrepix, i) ;
01958 }
01959 podata[col + row * ilx] = centreval + pixvalue ;
01960 }
01961 cpl_free(result) ;
01962 }
01963
01964
01965
01966 cpl_free(z) ;
01967 cpl_free(a) ;
01968 cpl_free(wave) ;
01969 cpl_free(emline) ;
01970 cpl_free(spec) ;
01971
01972 return retImage ;
01973 }
01974
02019 int sinfo_new_wavelength_calibration( cpl_image * image,
02020 FitParams ** par ,
02021 float ** bcoefs,
02022 float * wave,
02023 int n_lines,
02024 int ** row_clean,
02025 float ** wavelength_clean,
02026 int * n_found_lines,
02027 float dispersion,
02028 int halfWidth,
02029 float minAmplitude,
02030 float max_residual,
02031 float fwhm,
02032 int n_a_fitcoefs,
02033 int n_b_fitcoefs,
02034 float sigmaFactor,
02035 float pixel_tolerance )
02036
02037 {
02038 int i, j, k ;
02039 int n_fit ;
02040 int n_reject ;
02041 float * acoefs ;
02042 float * dacoefs ;
02043 float ** abuf ;
02044 float ** dabuf ;
02045 float chisq_poly, chisq_cross ;
02046 int zeroind ;
02047
02048 int lx=0;
02049 int ly=0;
02050 float* pdata=NULL;
02051
02052 if ( NULL == image )
02053 {
02054 sinfo_msg_error("no image given\n") ;
02055 return -1 ;
02056 }
02057 lx=cpl_image_get_size_x(image);
02058 ly=cpl_image_get_size_y(image);
02059 pdata=cpl_image_get_data_float(image);
02060
02061 if ( par == NULL )
02062 {
02063 sinfo_msg_error("no fit parameter data structure given\n") ;
02064 return -1 ;
02065 }
02066 if ( wave == NULL )
02067 {
02068 sinfo_msg_error("no wavelength list given\n") ;
02069 return -1 ;
02070 }
02071 if ( n_lines <= 0 )
02072 {
02073 sinfo_msg_error("impossible number of lines in line list given\n") ;
02074 return -1 ;
02075 }
02076 if ( row_clean == NULL )
02077 {
02078 sinfo_msg_error("no row_clean array given\n") ;
02079 return -1 ;
02080 }
02081 if ( wavelength_clean == NULL )
02082 {
02083 sinfo_msg_error("no wavelength_clean array given\n") ;
02084 return -1 ;
02085 }
02086
02087 if ( dispersion == 0. )
02088 {
02089 sinfo_msg_error("impossible dispersion given\n") ;
02090 return -1 ;
02091 }
02092
02093 if ( halfWidth <= 0 || halfWidth > ly/2 )
02094 {
02095 sinfo_msg_error("impossible half width of the fitting box given\n") ;
02096 return -1 ;
02097 }
02098 if ( minAmplitude < 1. )
02099 {
02100 sinfo_msg_error("impossible minimal amplitude\n") ;
02101 return -1 ;
02102 }
02103
02104 if ( max_residual <= 0. || max_residual > 1. )
02105 {
02106 sinfo_msg_error("impossible max_residual given\n") ;
02107 return -1 ;
02108 }
02109
02110 if ( fwhm <= 0. || fwhm > 10. )
02111 {
02112 sinfo_msg_error("impossible fwhm given\n") ;
02113
02114 return -1 ;
02115 }
02116
02117 if ( n_a_fitcoefs <= 0 || n_a_fitcoefs > 9 )
02118 {
02119 sinfo_msg_error("unrealistic n_a_fitcoefs given\n") ;
02120 return -1 ;
02121 }
02122
02123 if ( n_b_fitcoefs <= 0 || n_b_fitcoefs > 9 )
02124 {
02125 sinfo_msg_error(" unrealistic n_b_fitcoefs given\n") ;
02126 return -1 ;
02127 }
02128 if ( sigmaFactor <= 0. )
02129 {
02130 sinfo_msg_error(" impossible sigmaFactor given\n") ;
02131 return -1 ;
02132 }
02133
02134
02135 n_reject = 0 ;
02136 n_fit = 0 ;
02137
02138
02139
02140 if ( 0 > (n_fit = sinfo_new_fit_lines( image , par, fwhm,
02141 n_found_lines, row_clean,
02142 wavelength_clean,
02143 halfWidth, minAmplitude )) )
02144 {
02145 sinfo_msg_error(" cannot fit the lines, "
02146 "error code of sinfo_fitLines: %d\n", n_fit) ;
02147 return -1 ;
02148 }
02149
02150
02151 if ( -1 == sinfo_new_check_for_fake_lines (par, dispersion,
02152 wavelength_clean, row_clean,
02153 n_found_lines,
02154 lx, pixel_tolerance) )
02155 {
02156 sinfo_msg_error("cannot fit the lines, "
02157 "error code of sinfo_fitLines: %d", n_fit) ;
02158 return -1 ;
02159 }
02160
02161
02162 if (NULL == (acoefs = (float*) cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
02163 NULL == (dacoefs = (float*)cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
02164 NULL == (abuf = (float**) cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
02165 NULL == (dabuf = (float**) cpl_calloc (n_a_fitcoefs, sizeof(float*))) )
02166 {
02167 sinfo_msg_error(" cannot allocate memory\n") ;
02168 return -1 ;
02169 }
02170
02171 for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02172 {
02173 if ( NULL == (abuf[i] = (float*) cpl_calloc(lx, sizeof(float))) ||
02174 NULL == (dabuf[i] = (float*) cpl_calloc(lx, sizeof(float))) )
02175 {
02176 sinfo_msg_error(" cannot allocate memory\n") ;
02177 cpl_free(abuf) ;
02178 cpl_free(dabuf) ;
02179 return -1 ;
02180 }
02181 }
02182
02183
02184 k = 0 ;
02185
02186 for ( i = 0 ; i < lx ; i++ )
02187 {
02188 zeroind = 0 ;
02189 if ( FLT_MAX == (chisq_poly = sinfo_new_polyfit( par, i,
02190 n_found_lines[i],
02191 ly, dispersion,
02192 max_residual, acoefs,
02193 dacoefs, &n_reject,
02194 n_a_fitcoefs)) )
02195 {
02196
02197
02198
02199 for ( j = 0 ; j < n_a_fitcoefs ; j++ )
02200 {
02201 acoefs[j] = ZERO ;
02202 dacoefs[j] = ZERO ;
02203 }
02204 }
02205
02206 for ( j = 0 ; j < n_a_fitcoefs ; j++ )
02207 {
02208 if ( acoefs[0] <= 0. || acoefs[1] ==0. ||
02209 dacoefs[j] == 0. || isnan(acoefs[j]) )
02210 {
02211 zeroind = 1 ;
02212
02213 }
02214 }
02215 for ( j = 0 ; j < n_a_fitcoefs ; j++ )
02216 {
02217 if ( zeroind == 0 )
02218 {
02219 abuf[j][i] = acoefs[j] ;
02220 dabuf[j][i] = dacoefs[j] ;
02221 }
02222 else
02223 {
02224 abuf[j][i] = ZERO ;
02225 dabuf[j][i] = ZERO ;
02226 }
02227 }
02228 }
02229
02230
02231 for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02232 {
02233 if ( FLT_MAX == (chisq_cross = sinfo_new_coefs_cross_fit(lx,
02234 abuf[i],
02235 dabuf[i],
02236 bcoefs[i],
02237 n_b_fitcoefs,
02238 sigmaFactor)))
02239 {
02240 sinfo_msg_error (" cannot carry out the fitting of coefficients"
02241 " across the columns, for the coefficient with"
02242 " index: %d\n", i) ;
02243 for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02244 {
02245 cpl_free (abuf[i]) ;
02246 cpl_free (dabuf[i]) ;
02247 }
02248 cpl_free ( acoefs ) ;
02249 cpl_free ( dacoefs ) ;
02250 cpl_free ( abuf ) ;
02251 cpl_free ( dabuf ) ;
02252 return -1 ;
02253 }
02254 }
02255
02256
02257 for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02258 {
02259 cpl_free (abuf[i]) ;
02260 cpl_free (dabuf[i]) ;
02261 }
02262 cpl_free ( acoefs ) ;
02263 cpl_free ( dacoefs ) ;
02264 cpl_free ( abuf ) ;
02265 cpl_free ( dabuf ) ;
02266
02267 return 0 ;
02268 }
02269
02270
02282 cpl_image * sinfo_new_convolve_image_by_gauss( cpl_image * lineImage,
02283 int hw )
02284 {
02285 cpl_image * returnImage ;
02286 float* column_buffer=NULL ;
02287 float * filter ;
02288 int col, row ;
02289 int ilx=0;
02290 int ily=0;
02291 int olx=0;
02292 int oly=0;
02293 float* pidata=NULL;
02294 float* podata=NULL;
02295
02296 if ( lineImage == NULL )
02297 {
02298 sinfo_msg_error(" no input image given!\n") ;
02299 return NULL ;
02300 }
02301 ilx=cpl_image_get_size_x(lineImage);
02302 ily=cpl_image_get_size_y(lineImage);
02303 pidata=cpl_image_get_data_float(lineImage);
02304
02305 if ( hw < 1 )
02306 {
02307 sinfo_msg_error(" wrong half width given!\n") ;
02308 return NULL ;
02309 }
02310
02311
02312 if ( NULL == ( returnImage = cpl_image_new(ilx,ily,CPL_TYPE_FLOAT ) ))
02313 {
02314 sinfo_msg_error(" cannot allocate a new image\n");
02315 return NULL ;
02316 }
02317 olx=cpl_image_get_size_x(returnImage);
02318 oly=cpl_image_get_size_y(returnImage);
02319 podata=cpl_image_get_data_float(returnImage);
02320
02321 column_buffer=cpl_calloc(ily,sizeof(float)) ;
02322
02323
02324 for ( col = 0 ; col < ilx ; col++ )
02325 {
02326 for ( row = 0 ; row < ily ; row++ )
02327 {
02328 column_buffer[row] = pidata[col + row*ilx] ;
02329 }
02330
02331
02332
02333
02334
02335 filter = sinfo_function1d_filter_lowpass( column_buffer,
02336 ily,
02337 LOW_PASS_GAUSSIAN,
02338 hw ) ;
02339 for ( row = 0 ; row < ily ; row++ )
02340 {
02341 podata[col + row*ilx] = filter[row] ;
02342 }
02343 sinfo_function1d_del(filter) ;
02344 }
02345
02346 cpl_free(column_buffer);
02347 return returnImage ;
02348 }
02349
02389 cpl_image * sinfo_new_defined_resampling( cpl_image * image,
02390 cpl_image * calimage,
02391 int n_params,
02392 int* n_rows,
02393 double * dispersion,
02394 float * minval,
02395 float * maxval,
02396 double * centralLambda,
02397 int * centralpix )
02398 {
02399 cpl_image * retImage ;
02400 cpl_image * tempCalImage ;
02401 cpl_image * tempImage ;
02402 float lambda ;
02403 float dif, lambda_renorm ;
02404 float retimagecol[2560] ;
02405
02406 float* imagecol=NULL ;
02407 float* calcol=NULL ;
02408 float* x_renorm=NULL ;
02409
02410 float * imageptr ;
02411 float sum, new_sum ;
02412 float disp, mindisp ;
02413 int calcolpos[2560];
02414 int i, col, row, testrow ;
02415 int half_width, firstpos ;
02416 int dispInd ;
02417 int n ;
02418 int flag;
02419 float temprow;
02420 float minLambda = 0. ;
02421
02422
02423 double poly ;
02424
02425 int zeroind ;
02426 int ilx=0;
02427 int ily=0;
02428 int clx=0;
02429 int cly=0;
02430 int olx=0;
02431 int oly=0;
02432
02433 float* podata=NULL;
02434 float* pidata=NULL;
02435 float* pcdata=NULL;
02436 float* ptidata=NULL;
02437 float* ptcdata=NULL;
02438
02439
02440 if ( NULL == image )
02441 {
02442 sinfo_msg_error(" source image not given\n") ;
02443 return NULL ;
02444 }
02445 ilx=cpl_image_get_size_x(image);
02446 ily=cpl_image_get_size_y(image);
02447 pidata=cpl_image_get_data_float(image);
02448
02449
02450 if ( NULL == calimage )
02451 {
02452 sinfo_msg_error(" wavelength map image not given\n") ;
02453 return NULL ;
02454 }
02455 clx=cpl_image_get_size_x(calimage);
02456 cly=cpl_image_get_size_y(calimage);
02457 pcdata=cpl_image_get_data_float(calimage);
02458
02459 if ( ilx != clx ||
02460 ily != cly )
02461 {
02462 sinfo_msg_error("source image and wavelength map image "
02463 "are not compatible in size\n") ;
02464 return NULL ;
02465 }
02466
02467 if ( n_params < 1 )
02468 {
02469 sinfo_msg_error (" wrong number of fit parameters given\n") ;
02470 return NULL ;
02471 }
02472
02473 if ( n_params > 4 )
02474 {
02475 sinfo_msg_warning(" attention: very high number of fit "
02476 "parameters given, not tested !!!\n") ;
02477 }
02478
02479
02480 imagecol=cpl_calloc(ily,sizeof(float)) ;
02481 calcol=cpl_calloc(cly,sizeof(float)) ;
02482 x_renorm=cpl_calloc(n_params,sizeof(float)) ;
02483
02484
02485
02486
02487
02488
02489
02490
02491
02492
02493 dispInd = 0 ;
02494
02495 for ( col = 0 ; col < clx ; col++ )
02496 {
02497 if ( isnan(pcdata[col]) || pcdata[col] <= 0. )
02498 {
02499 continue ;
02500 }
02501 if ((pcdata[col] - pcdata[col+(clx)*(cly-1)]) > 0. )
02502 {
02503 dispInd-- ;
02504 }
02505 else if ((pcdata[col] - pcdata[col+(clx)*(cly-1)]) < 0. )
02506 {
02507 dispInd++ ;
02508 }
02509 else
02510 {
02511 continue ;
02512 }
02513 }
02514
02515 if ( dispInd == 0 )
02516 {
02517 sinfo_msg_error(" zero dispersion?\n");
02518 return NULL ;
02519 }
02520
02521
02522
02523 if ( dispInd < 0 )
02524 {
02525
02526 if ( NULL == ( tempCalImage = cpl_image_new(clx,cly,CPL_TYPE_FLOAT)))
02527 {
02528 sinfo_msg_error(" cannot allocate a new image\n");
02529 return NULL ;
02530 }
02531 ptcdata=cpl_image_get_data_float(tempCalImage);
02532 if ( NULL == ( tempImage = cpl_image_new( ilx, ily,CPL_TYPE_FLOAT)))
02533 {
02534 sinfo_msg_error(" cannot allocate a new image\n");
02535 cpl_image_delete(tempCalImage) ;
02536 return NULL ;
02537 }
02538 ptidata=cpl_image_get_data_float(tempImage);
02539
02540 for ( col = 0 ; col < clx ; col++ )
02541 {
02542 n = cly - 1 ;
02543 for ( row = 0 ; row < cly ; row++ )
02544 {
02545 ptcdata[col+row*clx] = pcdata[col+n*clx] ;
02546 ptidata[col+row*clx] = pidata[col+n*clx] ;
02547 n-- ;
02548 }
02549 }
02550 for ( i = 0 ; i < (int) ilx*ily ; i++ )
02551 {
02552 pidata[i] = ptidata[i] ;
02553 pcdata[i] = ptcdata[i] ;
02554 }
02555 cpl_image_delete(tempCalImage) ;
02556 cpl_image_delete(tempImage) ;
02557 }
02558
02559
02560 *maxval = -FLT_MAX ;
02561 *minval = FLT_MAX ;
02562 mindisp = FLT_MAX ;
02563 for ( col = 0 ; col < clx ; col++ )
02564 {
02565 if ( isnan(pcdata[col]) || pcdata[col] <= 0. )
02566 {
02567 continue ;
02568 }
02569 disp = (pcdata[col+(clx)*((cly)-1)]
02570 - pcdata[col]) / (float)cly ;
02571 if ( mindisp > disp )
02572 {
02573 mindisp = disp ;
02574 }
02575 if ( *minval >= pcdata[col] )
02576 {
02577 *minval = pcdata[col] ;
02578 }
02579 if ( *maxval <= pcdata[col + (clx)*((cly)-1)] )
02580 {
02581 *maxval = pcdata[col + (clx)*((cly)-1)] ;
02582 }
02583 }
02584
02585 if (*minval > 1.9 )
02586 {
02587 if ( cly > 1024 && cly < 3000)
02588 {
02589 *dispersion = DISPERSION_K_DITH ;
02590 *centralLambda = CENTRALLAMBDA_K ;
02591 }
02592 else if ( cly < 2000)
02593 {
02594 *dispersion = DISPERSION_K ;
02595 *centralLambda = CENTRALLAMBDA_K ;
02596 }
02597 else
02598 {
02599 *dispersion = DISPERSION_K_DITH/2 ;
02600 *centralLambda = CENTRALLAMBDA_K ;
02601 }
02602 }
02603 else if (*minval < 1.2 )
02604 {
02605 if ( cly > 1024 )
02606 {
02607 *dispersion = DISPERSION_J_DITH ;
02608 *centralLambda = CENTRALLAMBDA_J ;
02609 }
02610 else
02611 {
02612 *dispersion = DISPERSION_J ;
02613 *centralLambda = CENTRALLAMBDA_J ;
02614 }
02615 }
02616 else
02617 {
02618 if ( *maxval > 2.3 )
02619 {
02620 if ( cly > 1024 )
02621 {
02622 *dispersion = DISPERSION_HK_DITH ;
02623 *centralLambda = CENTRALLAMBDA_HK ;
02624 }
02625 else
02626 {
02627 *dispersion = DISPERSION_HK ;
02628 *centralLambda = CENTRALLAMBDA_HK ;
02629 }
02630 }
02631 else
02632 {
02633 if ( cly > 1024 )
02634 {
02635 *dispersion = DISPERSION_H_DITH ;
02636 *centralLambda = CENTRALLAMBDA_H ;
02637 }
02638 else
02639 {
02640 *dispersion = DISPERSION_H ;
02641 *centralLambda = CENTRALLAMBDA_H ;
02642 }
02643 }
02644 }
02645
02646
02647
02648
02649
02650 if ( (*maxval - *minval) / *dispersion < (float)cly )
02651 {
02652 sinfo_msg_error(" must be something wrong with the wavelength map!\n");
02653 return NULL ;
02654 }
02655
02656
02657 *n_rows = floor(floor(0.5+(*maxval - *minval) / *dispersion)/2+0.5)*2;
02658 *centralpix = *n_rows / 2 ;
02659 minLambda = *centralLambda - *dispersion * (float)*centralpix ;
02660
02661
02662
02663
02664
02665
02666
02667
02668
02669
02670
02671
02672 if ( NULL == ( retImage = cpl_image_new( ilx, *n_rows,CPL_TYPE_FLOAT ) ))
02673 {
02674 sinfo_msg_error(" cannot allocate a new image\n");
02675 return NULL ;
02676 }
02677 podata=cpl_image_get_data_float(retImage);
02678 olx=cpl_image_get_size_x(retImage);
02679 oly=cpl_image_get_size_y(retImage);
02680
02681
02682 for ( col = 0 ; col < olx ; col++ )
02683 {
02684
02685
02686
02687
02688
02689 sum = 0. ;
02690 for ( row = 0 ; row < ily ; row++ )
02691 {
02692 imagecol[row] = pidata[col + row*ilx] ;
02693 if (!isnan(imagecol[row]))
02694 {
02695 sum += imagecol[row] ;
02696 }
02697 calcol[row] = pcdata[col + row*clx] ;
02698 }
02699 for ( row = 0 ; row < oly ; row++ )
02700 {
02701 retimagecol[row] = 0. ;
02702 calcolpos[row] = -1;
02703 }
02704
02705 for ( row=0 ; row < cly ; row++)
02706 {
02707 temprow = (calcol[row]- minLambda)/ *dispersion;
02708 if (temprow >= 0 && temprow < oly)
02709 calcolpos[(int) temprow] = row;
02710 }
02711
02712 zeroind = 0 ;
02713
02714
02715
02716 for ( row = 0 ; row < oly ; row++ )
02717 {
02718 lambda = minLambda + *dispersion * (float) row ;
02719
02720
02721
02722
02723
02724 if ( row < cly )
02725 {
02726 if ( isnan(calcol[row]) )
02727 {
02728 zeroind = 1 ;
02729 }
02730 }
02731
02732 if ( (lambda < calcol[0]) ||
02733 (lambda > calcol[(cly)-1]) || zeroind == 1 )
02734 {
02735 retimagecol[row] = ZERO ;
02736 continue ;
02737 }
02738
02739
02740
02741
02742
02743
02744 if (calcolpos[row]==-1) {
02745 if(row>= (*n_rows-1)) calcolpos[row] = calcolpos[row-1];
02746 if(row< (*n_rows-1)) calcolpos[row] = calcolpos[row+1];
02747 }
02748 if (lambda-calcol[calcolpos[row]-1]==0.) {
02749 calcolpos[row]=calcolpos[row]-1;
02750 }
02751 testrow = calcolpos[row];
02752
02753
02754
02755
02756
02757
02758
02759
02760
02761
02762 if ( n_params % 2 == 0 )
02763 {
02764 half_width = (int)(n_params/2) - 1 ;
02765 }
02766 else
02767 {
02768 half_width = (int)(n_params/2) ;
02769 }
02770
02771
02772 if ( isnan(imagecol[testrow]) )
02773 {
02774 for ( i = row-half_width ; i < row-half_width+n_params ; i++ )
02775 {
02776 if (i < 0) continue ;
02777 if ( i >= oly ) continue ;
02778 retimagecol[i] = ZERO ;
02779 }
02780 imagecol[testrow] = 0. ;
02781 }
02782
02783 }
02784
02785
02786
02787
02788 new_sum = 0. ;
02789 for ( row = 0 ; row < oly ; row++ )
02790 {
02791 if ( isnan(retimagecol[row]) )
02792 {
02793 continue ;
02794 }
02795 lambda = minLambda + *dispersion * (float) row ;
02796
02797
02798
02799
02800
02801 if ( (lambda < calcol[0]) || (lambda > calcol[(cly)-1]) )
02802 {
02803 retimagecol[row] = ZERO ;
02804 continue ;
02805 }
02806
02807
02808
02809
02810
02811 if (calcolpos[row]==-1) {
02812 if(row >= (*n_rows-1)) calcolpos[row] = calcolpos[row-1];
02813 if(row < (*n_rows-1)) calcolpos[row] = calcolpos[row+1];
02814 }
02815 testrow = calcolpos[row];
02816
02817
02818
02819
02820
02821
02822
02823
02824
02825
02826 if ( n_params % 2 == 0 )
02827 {
02828 half_width = (int)(n_params/2) - 1 ;
02829 }
02830 else
02831 {
02832 half_width = (int)(n_params/2) ;
02833 }
02834
02835 firstpos = testrow - half_width ;
02836 if ( firstpos < 0 )
02837 {
02838 firstpos = 0 ;
02839 }
02840 else if ( firstpos > ((cly)-n_params) )
02841 {
02842 firstpos = cly - n_params ;
02843 }
02844
02845 if ( isnan(imagecol[firstpos]) )
02846 {
02847 retimagecol[row] = ZERO ;
02848 continue ;
02849 }
02850
02851
02852
02853
02854 dif = calcol[firstpos+n_params-1] - calcol[firstpos] ;
02855 for ( i = 0 ; i < n_params ; i++ )
02856 {
02857 x_renorm[i] = (calcol[firstpos + i] - calcol[firstpos]) / dif ;
02858 }
02859
02860
02861
02862 lambda_renorm = ( lambda - calcol[firstpos] ) / dif ;
02863
02864 imageptr = &imagecol[firstpos] ;
02865
02866 flag = 0;
02867 poly=sinfo_new_nev_ille(x_renorm, imageptr,
02868 n_params-1, lambda_renorm, &flag);
02869
02870 new_sum += poly ;
02871 retimagecol[row] = poly ;
02872 }
02873
02874
02875
02876 for ( row = 0 ; row < oly ; row++ )
02877 {
02878 if ( new_sum == 0. ) new_sum = 1. ;
02879 if ( isnan(retimagecol[row]) )
02880 {
02881 podata[col+row*olx] = ZERO ;
02882 }
02883 else
02884 {
02885
02886
02887 podata[col+row*olx] = retimagecol[row] ;
02888 }
02889 }
02890
02891 }
02892
02893 cpl_free(imagecol) ;
02894 cpl_free(calcol) ;
02895 cpl_free(x_renorm) ;
02896
02897 return retImage ;
02898 }
02899
02901