00001 /******************************************************************************
00002 * E.S.O. - VLT project
00003 *
00004 *
00005 *
00006 * who when what
00007 * -------- -------- ----------------------------------------------
00008 * schreib 13/07/00 created
00009 */
00010
00011 /************************************************************************
00012 * NAME
00013 * wave_calibration.c -
00014 * routines needed for wavelength calibration
00015 *
00016 * SYNOPSIS
00017 *
00018 * 1) FitParams ** newFitParams( int n_params )
00019 *
00020 * 2) void destroyFitParams ( FitParams ** params )
00021 *
00022 * 3) void dumpFitParamsToAscii ( FitParams ** params, char * filename )
00023 *
00024 * 4) void dumpAsciiToFitParams ( FitParams ** params, char * filename )
00025 *
00026 * 5) int findLines( OneImage * lineImage,
00027 * float * wave_position,
00028 * float * wave_intensity,
00029 * int n_lines,
00030 * int ** row_clean,
00031 * float ** wavelength_clean,
00032 * float beginWave,
00033 * float dispersion,
00034 * float mindiff,
00035 * int halfWidth,
00036 * int * n_found_lines,
00037 * float sigma,
00038 * int * sum_lines )
00039 *
00040 * 6) int readList( char * listname, float * lineCenter, float * lineIntensity )
00041 *
00042 *
00043 * 7) int linefit ( OneImage * mergedImage,
00044 * FitParams * par,
00045 * float fwhm,
00046 * int lineInd,
00047 * int column,
00048 * int halfWidth,
00049 * int lineRow,
00050 * float min_amplitude )
00051 *
00052 * 8) int fitLines ( OneImage * line_image,
00053 * FitParams ** allParams,
00054 * float fwhm,
00055 * int * n_lines,
00056 * int ** row,
00057 * float ** wavelength,
00058 * int width,
00059 * float min_amplitude )
00060 *
00061 * 9) float polyfit( FitParams ** par,
00062 * int column,
00063 * int n_lines,
00064 * int n_rows,
00065 * float dispersion,
00066 * float max_residual,
00067 * float * acoefs,
00068 * float * dacoefs,
00069 * int * n_reject,
00070 * int n_fitcoefs )
00071 *
00072 * 10) float coefsCrossFit ( int n_columns,
00073 * float * acoefs,
00074 * float * dacoefs,
00075 * float * bcoefs,
00076 * int n_fitcoefs,
00077 * float sigma_factor )
00078 *
00079 *
00080 * 11) OneImage * waveMap( OneImage * lineImage,
00081 * float ** bcoefs,
00082 * int n_a_fitcoefs,
00083 * int n_b_fitcoefs,
00084 * float * wavelength,
00085 * float * intensity,
00086 * int n_lines,
00087 * int magFactor,
00088 * int * bad_column_mask,
00089 * int n_bad_columns )
00090 *
00091 * 12) int wavelengthCalibration( OneImage * image,
00092 * FitParams ** par ,
00093 * float ** bcoefs,
00094 * float * wave,
00095 * int n_lines,
00096 * int ** row_clean,
00097 * float ** wavelength_clean,
00098 * int * n_found_lines,
00099 * float dispersion,
00100 * int halfWidth,
00101 * float minAmplitude,
00102 * float max_residual,
00103 * float fwhm,
00104 * int n_a_fitcoefs,
00105 * int n_b_fitcoefs,
00106 * float sigmaFactor )
00107 *
00108 * 13) OneImage * convolveImageByGauss( OneImage * lineImage,
00109 * int hw )
00110 *
00111 * 14) OneImage * definedResampling( OneImage * image,
00112 * OneImage * calimage,
00113 * int n_params,
00114 * int n_rows,
00115 * double * dispersion,
00116 * float * minval,
00117 * float * maxval,
00118 * double * centralLambda,
00119 * int * centralpix )
00120 *
00121 * DESCRIPTION
00122 *
00123 * 1) allocates memory for a new vector of
00124 * FitParams data structures
00125 * 2) frees memory of a vector of FitParams data structures
00126 * 3) dumps the fit parameters to an ASCII file
00127 * 4) dumps ASCII information to an allocated FitParams data structure
00128 * 5) determines the pixel shift between the line list
00129 * and the real image by using the beginning wavelength
00130 * on the detector and the dispersion estimate.
00131 * 6) reads the line data of the calibration lamps
00132 * 7) fits a gaussian to a 1-dimensional slice of an image,
00133 * this routine uses the routine lsqfit_c as a non-linear
00134 * least square fit method (Levenberg-Marquardt).
00135 * 8) calculates and stores the fit parameters of the neon
00136 * emission lines of a neon frame by using the linefit
00137 * routine.
00138 * 9) fits a second order polynom
00139 * lambda[i] = a1 + a2*pos[i] + a3*pos[i]^2
00140 * to determine the connection between the listed wave-
00141 * length values and the gauss-fitted positions for each
00142 * image column using the singular value decomposition
00143 * method.
00144 * 10) Fits the each single parameter of the three fit parameters
00145 * acoefs from polyfit through the image columns
00146 * 11) this routine determines a wavelength calibration map
00147 * frame associating a wavelength value to each pixel
00148 * by using the fit coefficients determined before.
00149 * 12) this routine takes an image from a calibration
00150 * emission lamp and delivers the fit coefficients of
00151 * a polynomial fit across the columns
00152 * of the coefficients of the polynomial line position
00153 * fits as output. Furthermore it delivers an array of the fit parameters
00154 * as output. This routine expects Nyquist sampled spectra
00155 * (either an interleaved image or an image convolved with an
00156 * appropriate function in spectral direction)
00157 * 13) convolves an emission line image with a gaussian
00158 * with user given integer half width by using the eclipse
00159 * routine function1d_filter_lowpass().
00160 * 14) Given a source image and a corresponding wavelength
00161 * calibration file this routine produces an image
00162 * in which elements in a given row are associated
00163 * with a single wavelength. It thus corrects for
00164 * the wavelength shifts between adjacent elements
00165 * in the rows of the input image. The output image
00166 * is larger in the wavelength domain than the input
00167 * image with pixels in each column corresponding to
00168 * undefined (blank, ZERO) values. The distribution
00169 * of these undefined values varies from column to
00170 * column. The input image is resampled at discrete
00171 * wavelength intervals using a polynomial interpolation
00172 * routine.
00173 * The wavelength intervals (dispersion) and the
00174 * central wavelength are defined and stable for each
00175 * used grating. Thus, each row has a defined wavelength
00176 * for each grating. Only the number of rows can be
00177 * changed by the user.
00178 *
00179 * FILES
00180 *
00181 * ENVIRONMENT
00182 *
00183 * RETURN VALUES
00184 *
00185 * CAUTIONS
00186 *
00187 * EXAMPLES
00188 *
00189 * SEE ALSO
00190 *
00191 * BUGS
00192 *
00193 *------------------------------------------------------------------------
00194 */
00195
00196 #include "vltPort.h"
00197
00198 /*
00199 * System Headers
00200 */
00201
00202 /*
00203 * Local Headers
00204 */
00205
00206 #include "wave_calibration.h"
00207 #include "solve_poly_root.h"
00208
00209
00210 /*----------------------------------------------------------------------------
00211 Function : newFitParams()
00212 In : number of spectral lines that will be fitted
00213 Out : allocated array of FitParams data structures
00214 Job : allocates memory for a new array of
00215 FitParams data structures
00216 ---------------------------------------------------------------------------*/
00217
00218 FitParams ** newFitParams( int n_params )
00219 {
00220 FitParams ** new_params =NULL;
00221 FitParams * temp_params =NULL;
00222 float * temp_fit_mem =NULL;
00223 float * temp_derv_mem=NULL;
00224 int i ;
00225
00226
00227 if ( n_params <= 0 )
00228 {
00229 cpl_msg_error ("newFitParams:"," wrong number of lines to fit\n") ;
00230 return NULL ;
00231 }
00232
00233 if ( NULL == (new_params = (FitParams **) cpl_calloc ( n_params , sizeof (FitParams*) ) ) )
00234 {
00235 cpl_msg_error ("newFitParams:"," could not allocate memory\n") ;
00236 return NULL ;
00237 }
00238 if ( NULL == (temp_params = cpl_calloc ( n_params , sizeof (FitParams) ) ) )
00239 {
00240 cpl_msg_error ("newFitParams:"," could not allocate memory\n") ;
00241 return NULL ;
00242 }
00243 if ( NULL == (temp_fit_mem = (float *) cpl_calloc( n_params*MAXPAR, sizeof (float) ) ) )
00244 {
00245 cpl_msg_error ("newFitParams:"," could not allocate memory\n") ;
00246 return NULL ;
00247 }
00248
00249
00250 if ( NULL == (temp_derv_mem = (float *) cpl_calloc( n_params*MAXPAR, sizeof (float) ) ) )
00251 {
00252 cpl_msg_error ("newFitParams:"," could not allocate memory\n") ;
00253 return NULL ;
00254 }
00255
00256 for ( i = 0 ; i < n_params ; i ++ )
00257 {
00258 new_params[i] = temp_params+i;
00259 new_params[i] -> fit_par = temp_fit_mem+i*MAXPAR;
00260 new_params[i] -> derv_par = temp_derv_mem+i*MAXPAR;
00261 new_params[i] -> column = 0 ;
00262 new_params[i] -> line = 0 ;
00263 new_params[i] -> wavelength = 0. ;
00264 new_params[i] -> n_params = n_params ;
00265 }
00266
00267 return new_params ;
00268 }
00269
00270 /*----------------------------------------------------------------------------
00271 Function : destroyFitParams()
00272 In : fit params to destroy
00273 Out : nothing
00274 Job : frees memory of an array of FitParams data structures
00275 ---------------------------------------------------------------------------*/
00276
00277 void destroyFitParams ( FitParams ** params )
00278 {
00279
00280 if ( params == NULL )
00281 {
00282 return ;
00283 }
00284
00285 cpl_free ( params[0] -> fit_par ) ;
00286 cpl_free ( params[0] -> derv_par ) ;
00287 cpl_free ( params[0] ) ;
00288 cpl_free ( params ) ;
00289 }
00290
00291 /*----------------------------------------------------------------------------
00292 Function : dumpFitParamsToAscii()
00293 In : params: fit params to dump
00294 filename
00295 Out : filled ASCII file
00296 Job : dumps the fit parameters to an ASCII file
00297 ---------------------------------------------------------------------------*/
00298
00299 void dumpFitParamsToAscii ( FitParams ** params, char * filename )
00300 {
00301 FILE * fp ;
00302 int i ;
00303
00304 if ( NULL == params )
00305 {
00306 cpl_msg_error ("dumpFitParamsToAscii:"," no fit parameters available!\n") ;
00307 return ;
00308 }
00309
00310 if ( NULL == filename )
00311 {
00312 cpl_msg_error ("dumpFitParamsToAscii:"," no filename available!\n") ;
00313 return ;
00314 }
00315
00316 if ( NULL == (fp = fopen ( filename, "w" ) ) )
00317 {
00318 cpl_msg_error("dumpFitParamsToAscii:"," cannot open %s\n", filename) ;
00319 return ;
00320 }
00321
00322 for ( i = 0 ; i < params[0] -> n_params ; i++ )
00323 {
00324 fprintf(fp, "%d %d %d %f %f %f %f %f %f %f %f %f\n",
00325 params[i]->n_params, params[i]->column, params[i]->line,
00326 params[i]->wavelength, params[i]->fit_par[0], params[i]->fit_par[1],
00327 params[i]->fit_par[2], params[i]->fit_par[3],
00328 params[i]->derv_par[0], params[i]->derv_par[1],
00329 params[i]->derv_par[2], params[i]->derv_par[3] ) ;
00330 }
00331 fclose(fp) ;
00332 }
00333
00334 /*----------------------------------------------------------------------------
00335 Function : dumpAsciiToFitParams()
00336 In : allocated dummy for the fit params
00337 filename
00338 Out : params: filled FitParams object
00339 Job : dumps ASCII information to an allocated FitParams data structure
00340 ---------------------------------------------------------------------------*/
00341
00342 void dumpAsciiToFitParams ( FitParams ** params, char * filename )
00343 {
00344 FILE * fp ;
00345 int i ;
00346
00347 if ( NULL == params )
00348 {
00349 cpl_msg_error ("dumpAsciiToFitParams:"," no fit parameters available!\n") ;
00350 return ;
00351 }
00352
00353 if ( NULL == filename )
00354 {
00355 cpl_msg_error ("dumpAsciiToFitParams:"," no filename available!\n") ;
00356 return ;
00357 }
00358
00359 if ( NULL == (fp = fopen ( filename, "r" ) ) )
00360 {
00361 cpl_msg_error("dumpAsciiToFitParams:"," cannot open %s\n", filename) ;
00362 return ;
00363 }
00364
00365 for ( i = 0 ; i < params[0]->n_params ; i++ )
00366 {
00367 fscanf(fp, "%d %d %d %f %f %f %f %f %f %f %f %f\n",
00368 ¶ms[i]->n_params, ¶ms[i]->column, ¶ms[i]->line,
00369 ¶ms[i]->wavelength, ¶ms[i]->fit_par[0], ¶ms[i]->fit_par[1],
00370 ¶ms[i]->fit_par[2], ¶ms[i]->fit_par[3],
00371 ¶ms[i]->derv_par[0], ¶ms[i]->derv_par[1],
00372 ¶ms[i]->derv_par[2], ¶ms[i]->derv_par[3] ) ;
00373 }
00374 fclose(fp) ;
00375 }
00376
00377 /*----------------------------------------------------------------------------
00378 Function : findLines()
00379 In : lineImage: merged emission line image,
00380 wave_position: wavelength list in Angstroems
00381 wave_intensity: corresponding intensity list
00382 n_lines: number of lines in list
00383 row_clean: resulting list of the row indices but without the
00384 lines that are too close to each other for the fit
00385 wavelength_clean: corrected wavelength list corresponding to
00386 the row_clean array
00387 beginWave: beginning wavelength on detector in microns
00388 dispersion: dispersion of the grating on the detector
00389 (microns per pixel, attention: merged image).
00390 mindiff: minimal difference of mean and median column
00391 intensity to do the correlation.
00392 This is done to avoid correlations in columns
00393 without emission line intensity.
00394 halfWidth: half width of the box where the line must sit,
00395 n_found_lines: number of found and correlated
00396 emission lines in a column.
00397 sigma: sigma of Gaussian that is convolved with the artificial spectrum
00398 Out : 0 if all went o.k.
00399 row: resulting list of the row indices of the line positions
00400 row_clean: resulting list of the row indices but without the
00401 lines that are too close to each other for the fit
00402 wavelength: wavelength from the list corresponding to the
00403 found row positions
00404 wavelength_clean: corrected wavelength list corresponding to
00405 the row_clean array
00406 n_found_lines: number of found and correlated
00407 emission lines in a column.
00408 sum_lines: total sum of found and correlated emission lines
00409 -1 if something has gone wrong
00410 Job : determines the pixel shift between the line list
00411 and the real image by using the beginning wavelength
00412 on the detector and the dispersion estimate.
00413 ---------------------------------------------------------------------------*/
00414
00415 int findLines( OneImage * lineImage,
00416 float * wave_position,
00417 float * wave_intensity,
00418 int n_lines,
00419 int ** row_clean,
00420 float ** wavelength_clean,
00421 float beginWave,
00422 float dispersion1,
00423 float dispersion2,
00424 float mindiff,
00425 int halfWidth,
00426 int * n_found_lines,
00427 float sigma,
00428 int * sum_lines )
00429 {
00430 int ** row ;
00431 float ** wavelength ;
00432 float buf1, buf2 ;
00433 float meanval ;
00434 float colmedian ;
00435 float * column, * tempcol ;
00436 float * lines ;
00437 float * conv_lines ;
00438 float * wave_buffer ;
00439 float * wave_mem;
00440 int * dummy_row ;
00441 int i, j, k, m ;
00442 int position ;
00443 int gmax, gmin ;
00444 int col ;
00445 int * row_mem;
00446 float sum ;
00447 float angst ;
00448
00449 if ( NullImage == lineImage )
00450 {
00451 cpl_msg_error ("findLines:"," no image given\n") ;
00452 return -1 ;
00453 }
00454
00455 if ( n_lines <= 0 || NULL == wave_position )
00456 {
00457 cpl_msg_error("findLines:"," no line list given\n") ;
00458 return -1 ;
00459 }
00460 if ( NULL == wave_intensity )
00461 {
00462 cpl_msg_error("findLines:"," no intensity list given\n") ;
00463 return -1 ;
00464 }
00465
00466 if ( dispersion1 == 0. )
00467 {
00468 cpl_msg_error("findLines:"," wrong dispersion given\n") ;
00469 return -1 ;
00470 }
00471
00472 if ( row_clean == NULL )
00473 {
00474 cpl_msg_error("findLines:"," row_clean array is not allocated\n") ;
00475 return -1 ;
00476 }
00477
00478 if ( wavelength_clean == NULL )
00479 {
00480 cpl_msg_error("findLines:"," wavelength_clean array is not allocated\n") ;
00481 return -1 ;
00482 }
00483
00484 if ( beginWave <= 0. )
00485 {
00486 cpl_msg_error ("findLines:"," impossible beginWave given\n") ;
00487 return -1 ;
00488 }
00489 if ( mindiff < -100. )
00490 {
00491 cpl_msg_error ("findLines:"," wrong mindiff value\n") ;
00492 return -1 ;
00493 }
00494
00495 if ( halfWidth <= 0 )
00496 {
00497 cpl_msg_error("findLines:"," wrong half width of fit box given\n") ;
00498 return -1 ;
00499 }
00500
00501 if ( n_found_lines == NULL )
00502 {
00503 cpl_msg_error("findLines:"," n_found_lines not allocated\n") ;
00504 return -1 ;
00505 }
00506
00507 if ( sigma <= 0. || sigma >= lineImage -> ly / 2)
00508 {
00509 cpl_msg_error("findLines:"," wrong sigma given\n") ;
00510 return -1 ;
00511 }
00512
00513 /* allocate memory */
00514 row = (int**) cpl_calloc( lineImage->lx, sizeof(int*)) ;
00515 wavelength = (float**) cpl_calloc( lineImage->lx, sizeof(float*)) ;
00516 row_mem = cpl_calloc( n_lines*lineImage->lx, sizeof(int) ) ;
00517 wave_mem = cpl_calloc( n_lines*lineImage->lx, sizeof(float) ) ;
00518 for ( i = 0 ; i <lineImage -> lx ; i++ )
00519 {
00520 row[i] = row_mem + i*n_lines;
00521 wavelength[i] = wave_mem + i*n_lines;
00522 }
00523
00524 /* find if the wavelength is given in microns, nanometers or Angstroem */
00525 if ( wave_position[0] > 10000. )
00526 {
00527 /* Angstroem */
00528 angst = 10000. ;
00529 }
00530 else if ( wave_position[0] > 1000. && wave_position[0] < 10000. )
00531 {
00532 /* nanometers */
00533 angst = 1000. ;
00534 }
00535 else
00536 {
00537 /* microns */
00538 angst = 1. ;
00539 }
00540
00541 /*-------------------------------------------------------------------------------
00542 * compute the mean and median intensity value in the given column
00543 * and determine if there is enough intensity in the column to do the correlation
00544 */
00545 tempcol = (float*) cpl_calloc(lineImage->ly, sizeof(float)) ;
00546 *sum_lines = 0 ;
00547 buf1 = 0. ;
00548 buf2 = 0. ;
00549 /* allocate memory */
00550 column = (float*) cpl_calloc(lineImage -> ly, sizeof (float)) ;
00551 lines = (float*) cpl_calloc(lineImage -> ly, sizeof (float)) ;
00552 conv_lines = (float*) cpl_calloc(lineImage -> ly, sizeof (float)) ;
00553 wave_buffer = (float*) cpl_calloc(lineImage -> ly, sizeof (float)) ;
00554 dummy_row = (int*) cpl_calloc(n_lines, sizeof(int)) ;
00555
00556 for ( col = 0 ; col < lineImage->lx ; col++ )
00557 {
00558 n_found_lines[col] = 0 ;
00559 sum = 0. ;
00560 for ( i = 0 ; i < lineImage -> ly ; i++ )
00561 {
00562 if (qfits_isnan(lineImage->data[col + i*lineImage->lx]) )
00563 {
00564 tempcol[i] = 0. ;
00565 continue ;
00566 }
00567
00568 sum = sum + lineImage -> data[col + i*lineImage -> lx] ;
00569 tempcol[i] = lineImage -> data[col + i*lineImage -> lx];
00570 }
00571 meanval = sum / lineImage -> ly ;
00572 /* lets assume the median is the background */
00573 colmedian = median ( tempcol, lineImage->ly);
00574
00575 if ( meanval - colmedian < mindiff )
00576 {
00577 cpl_msg_warning("findLines:"," sorry, not enough intensity (mean: %f, diff: %f) in image column: %d to correlate emission lines\n", meanval, meanval - colmedian,col) ;
00578 continue ;
00579 }
00580
00581 for ( i = 0 ; i < lineImage -> ly ; i++ )
00582 {
00583 conv_lines[i]=0;
00584 wave_buffer[i]=0;
00585 }
00586 for ( i = 0 ; i < n_lines ; i++ )
00587 {
00588 dummy_row[i] = 0;
00589 }
00590
00591 /* go through the column with index col */
00592 for ( i = 0 ; i < lineImage -> ly ; i++ )
00593 {
00594 if ( qfits_isnan(lineImage->data[col+i*lineImage->lx]) )
00595 {
00596 column[i] = 0. ;
00597 }
00598 else
00599 {
00600 column[i] = lineImage -> data[col + i*lineImage -> lx] ;
00601 }
00602
00603 /* determine the line position on the pixels */
00604 /*lines[i] = (dispersion * (float) i + beginWave) * angst ;*/
00605 lines[i] = (dispersion1 * (float) (i-lineImage->ly/2) + dispersion2 * (float) (i-lineImage->ly/2) * (float) (i-lineImage->ly/2) + beginWave) * angst ;
00606
00607 /* ---------------------------------------------------------------
00608 * find the nearest line position for each wavelength in the list
00609 * and set this position to the given line intensity as weight
00610 */
00611 for ( j = 0 ; j < n_lines ; j ++ )
00612 {
00613 buf1 = 0. ;
00614 buf2 = 0. ;
00615 if ( (wave_position[j] >= (lines[i] - fabs(dispersion1)/2.*angst)) &&
00616 (wave_position[j] <= (lines[i] + fabs(dispersion1)/2.*angst)) )
00617 {
00618 buf1 = wave_intensity[j] ; /* set the given line intensity as weight */
00619 buf2 = wave_position[j] ;
00620 break ;
00621 }
00622 }
00623 lines[i] = buf1 ;
00624 wave_buffer[i] = buf2 ; /* get the wavelength associated with the corresponding
00625 found emission line */
00626
00627 /* convolve the artificial spectrum by a Gaussian with given sigma value */
00628 if ( lines[i] != 0. )
00629 {
00630 /* consider only +- 2 sigma */
00631 gmin = nint((float) i - 2. * sigma) ;
00632 gmax = nint((float) i + 2. * sigma) ;
00633
00634 /* be aware of image margins */
00635 if ( gmin < 0 )
00636 {
00637 gmin = 0 ;
00638 }
00639 if ( gmax >= lineImage -> ly )
00640 {
00641 gmax = lineImage -> ly - 1 ;
00642 }
00643 for ( j = gmin ; j <= gmax ; j++ )
00644 {
00645 conv_lines[j] +=
00646 lines[i] * exp( (double)( -0.5*((j - i)*(j - i)))/(double) (sigma*sigma) ) ;
00647 }
00648 }
00649 }
00650
00651 /* do the cross correlation */
00652 position = INT32_MAX ;
00653 position = correlation(column, conv_lines, lineImage -> ly ) ;
00654 if ( abs (position) > lineImage -> ly / 4 )
00655 {
00656 cpl_msg_warning("findLines:"," sorry, shift of artificial data relative to image (%d) seems to be too high in column: %d\n", position, col) ;
00657 continue ;
00658 }
00659
00660 j = 0 ;
00661 for ( i = 0 ; i < lineImage -> ly ; i ++ )
00662 {
00663 if ( lines[i] != 0.0 )
00664 {
00665 if ( (i - position) >= 0 && (i - position) < lineImage -> ly )
00666 {
00667 row[col][j] = i - position ;
00668 /* get the wavelength corresponding to found line row index */
00669 wavelength[col][j] = wave_buffer[i] / angst ;
00670 j++ ;
00671 }
00672 }
00673 }
00674
00675
00676 /* -------------------------------------------------------------------------------------------
00677 * determine the row_clean array, that means, take only the row values if the distance between
00678 * adjacent lines is large enough for the fit
00679 */
00680 for ( k = 1 ; k <= j ; k ++ )
00681 {
00682 if (dummy_row[k-1] != -1)
00683 {
00684 dummy_row[k-1] = row[col][k-1] ;
00685 }
00686 if ( (row[col][k] - row[col][k-1]) <= 2*halfWidth )
00687 {
00688 dummy_row[k-1] = -1 ;
00689 if (k<n_lines)
00690 dummy_row[k] = -1 ;
00691 }
00692 if ( (row[col][k+1] - row[col][k]) <= 2*halfWidth)
00693 {
00694 if (k<n_lines)
00695 dummy_row[k] = -1 ;
00696 if (k+1<n_lines)
00697 dummy_row[k+1] = -1 ;
00698 }
00699 }
00700
00701 m = 0 ;
00702 for ( k = 0 ; k < j ; k ++ )
00703 {
00704 if ( dummy_row[k] != -1 && dummy_row[k] != 0 )
00705 {
00706 row_clean[col][m] = dummy_row[k] ;
00707 wavelength_clean[col][m] = wavelength[col][k] ;
00708 m ++ ;
00709 }
00710 }
00711
00712 n_found_lines[col] = m ;
00713
00714 *sum_lines += n_found_lines[col] ;
00715 }
00716 cpl_free (column) ;
00717 cpl_free (lines) ;
00718 cpl_free (conv_lines) ;
00719 cpl_free (dummy_row) ;
00720 cpl_free (wave_buffer) ;
00721 cpl_free (row_mem) ;
00722 cpl_free (wave_mem) ;
00723 cpl_free (tempcol) ;
00724 cpl_free (row) ;
00725 cpl_free (wavelength) ;
00726
00727 return 0 ;
00728 }
00729
00730 /*----------------------------------------------------------------------------
00731 Function : readList()
00732 In : name of the list file, arrays to store the wavelength
00733 and the intensities of the emission lines
00734 Out : number of lines in file
00735 Job : reads the line data of the calibration lamps
00736 ---------------------------------------------------------------------------*/
00737
00738 int readList( char * listname, float * lineCenter, float * lineIntensity )
00739 {
00740 FILE * fp ;
00741 int i, n_lines ;
00742
00743 if ( NULL == lineCenter )
00744 {
00745 cpl_msg_error("readList:"," lineCenter array is not allocated\n") ;
00746 return -1 ;
00747 }
00748
00749 if ( NULL == lineIntensity )
00750 {
00751 cpl_msg_error("readList:"," lineIntensity array is not allocated\n") ;
00752 return -1 ;
00753 }
00754
00755 if ( NULL == (fp = fopen ( listname, "r" ) ) )
00756 {
00757 cpl_msg_error("readList:"," cannot open %s\n", listname) ;
00758 return -1 ;
00759 }
00760
00761 i = 0 ;
00762 while ( fscanf( fp, "%f %f", &lineCenter[i], &lineIntensity[i] ) != EOF )
00763 {
00764 i ++ ;
00765 }
00766 n_lines = i ;
00767 fclose(fp) ;
00768
00769 return n_lines ;
00770 }
00771
00772
00773 /*----------------------------------------------------------------------------
00774 Function : linefit()
00775 In : mergedImage: image of a calibration emission lamp,
00776 par: dummys for the resulting fitting parameters,
00777 fwhm: guess value for full width of half maximum of gaussian
00778 lineInd: index of the emission line,
00779 column: present index of the image column,
00780 halfWidth: half width of the box where the line must sit,
00781 lineRow: row index where the line is expected,
00782 min_amplitude: minimum line amplitude with respect to the background
00783 to do the fit
00784 Out : the fitting parameter data structure containing the
00785 resulting parameters.
00786 integers:
00787 number of iterations if all was ok,
00788 -8 if no input image was given,
00789 -9 if no dummy for the fit parameters is given,
00790 -10 if the wrong column index was given,
00791 -11 if the wrong box width was given,
00792 -12 if the wrong row index was given,
00793 -13 if a wrong minimum amplitude factor was given
00794 -14 if the spectral vector data structure memory
00795 could not be allocated
00796 -15 wrong row index or box width was given,
00797 -16 signal too low to fit,
00798 -17 least squares fit failed
00799 Job : fits a gaussian to a 1-dimensional slice of an image,
00800 this routine uses the routine lsqfit_c as a non-linear
00801 least square fit method (Levenberg-Marquardt).
00802 ---------------------------------------------------------------------------*/
00803
00804 int linefit ( OneImage * mergedImage,
00805 FitParams * par,
00806 float fwhm,
00807 int lineInd,
00808 int column,
00809 int halfWidth,
00810 int lineRow,
00811 float min_amplitude,
00812 Vector * line,
00813 int * mpar,
00814 float * xdat,
00815 float * wdat )
00816 {
00817 int i, j ;
00818 int iters, xdim, ndat ;
00819 int numpar, its ;
00820 int position ;
00821 float maxval, tol, lab ;
00822
00823 if ( mergedImage == NullImage )
00824 {
00825 cpl_msg_error ("linefit:"," no image given as input\n") ;
00826 return -8 ;
00827 }
00828 if ( par == NULL )
00829 {
00830 cpl_msg_error("linefit:"," fit parameters not given\n") ;
00831 return -9 ;
00832 }
00833 if ( column < 0 || column > mergedImage -> lx )
00834 {
00835 cpl_msg_error ("linefit:"," wrong column number\n") ;
00836 return -10 ;
00837 }
00838 if ( halfWidth < 0 || halfWidth > mergedImage -> ly )
00839 {
00840 cpl_msg_error ("linefit:"," wrong width given\n") ;
00841 return -11 ;
00842 }
00843 if ( lineRow < 0 || lineRow > mergedImage -> ly )
00844 {
00845 cpl_msg_error ("linefit:"," wrong number of row of the line given\n") ;
00846 return -12 ;
00847 }
00848 if ( min_amplitude < 1. )
00849 {
00850 cpl_msg_error ("linefit:"," wrong minimum amplitude\n") ;
00851 return -13 ;
00852 }
00853
00854 /* initialise the Vector */
00855 for ( i = 0 ; i < line -> n_elements ; i++)
00856 {
00857 line->data[i] = 0;
00858 }
00859
00860 par -> column = column ;
00861 par -> line = lineInd ;
00862
00863 /* determine the values of the spectral vector given as input */
00864 /* go through the chosen column */
00865
00866 j = 0 ;
00867 for ( i = lineRow-halfWidth ; i <= lineRow+halfWidth ; i++ )
00868 {
00869 if ( i < 0 || i >= mergedImage -> ly )
00870 {
00871 cpl_msg_error ("linefit:"," wrong line position or width given\n") ;
00872 return -15 ;
00873 }
00874 else
00875 {
00876 line -> data[j] = mergedImage -> data[column + i*mergedImage -> lx] ;
00877 j ++ ;
00878 }
00879 }
00880
00881 /*--------------------------------------------------------------------
00882 * go through the spectral vector
00883 * determine the maximum pixel value in the spectral vector
00884 */
00885 maxval = -FLT_MAX ;
00886 position = -INT32_MAX ;
00887 for ( i = 0 ; i < line -> n_elements ; i++ )
00888 {
00889 xdat[i] = i ;
00890 wdat[i] = 1.0 ;
00891 if ( line -> data[i] >= maxval )
00892 {
00893 maxval = line -> data[i] ;
00894 position = i ;
00895 }
00896 }
00897
00898 /* set initial values for the fitting routine */
00899 xdim = XDIM ;
00900 ndat = line -> n_elements ;
00901 numpar = MAXPAR ;
00902 tol = TOL ;
00903 lab = LAB ;
00904 its = ITS ;
00905 par -> fit_par[1] = fwhm ;
00906 par -> fit_par[2] = (float) position ;
00907 par -> fit_par[3] = (float) (line -> data[0] + line -> data[line->n_elements - 1]) / 2.0 ;
00908 par -> fit_par[0] = maxval - (par -> fit_par[3]) ;
00909
00910 /* exclude low signal cases */
00911 if ( par->fit_par[0] < min_amplitude )
00912 {
00913 cpl_msg_debug ("linefit:"," sorry, amplitude of line too low to fit: %f\n",par->fit_par[0] ) ;
00914 return -16 ;
00915 }
00916
00917 for ( i = 0 ; i < MAXPAR ; i++ )
00918 {
00919 par -> derv_par[i] = 0.0 ;
00920 mpar[i] = 1 ;
00921 }
00922
00923 /* finally, do the least square fit using a gaussian */
00924 if ( 0 > ( iters = lsqfit_c( xdat, &xdim, line -> data, wdat, &ndat, par -> fit_par,
00925 par -> derv_par, mpar, &numpar, &tol, &its, &lab )) )
00926 {
00927 cpl_msg_debug ("linefit:"," lsqfit_c: least squares fit failed, error no.: %d\n", iters) ;
00928 return -17 ;
00929 }
00930
00931 /* correct the fitted position for the given row of the line in image coordinates */
00932 par -> fit_par[2] = (float) (lineRow - halfWidth) + par -> fit_par[2] ;
00933
00934 /* all was o.k. */
00935 return iters ;
00936 }
00937
00938 /*----------------------------------------------------------------------------
00939 Function : fitLines
00940 In : line_image: merged image of a calibration lamp ,
00941 allParams: allocated vector of FitParams data structures,
00942 fwhm: guess value for full width of half maximum of gaussian
00943 n_lines: number of neon lines that will be fitted in one column ,
00944 row: list of the rows of the fitted lines
00945 wavelength: list of wavelength corresponding to the found line rows
00946 width: half width of a box around the found rows within the line is fitted
00947 min_amplitude: minimum line amplitude with respect to the background
00948 to do the fit
00949 Out : filled FitParams data structure vector,
00950 number of successfully fitted lines,
00951 errors: negative integers resulting from the linefit
00952 routine and:
00953 -18: no image given,
00954 -19: number of emission lines or number of slitlets is wrong,
00955 -20: vector of the slitlet boundaries or of the line rows or
00956 of the half width are empty.
00957 -21: no wavelength array given.
00958 Job : calculates and stores the fit parameters of the neon
00959 emission lines of a neon frame by using the linefit
00960 routine.
00961 ---------------------------------------------------------------------------*/
00962
00963 int fitLines ( OneImage * line_image,
00964 FitParams ** allParams,
00965 float fwhm,
00966 int * n_lines,
00967 int ** row,
00968 float ** wavelength,
00969 int width,
00970 float min_amplitude )
00971 {
00972 int i, k, l ;
00973 int result ;
00974 Vector * line;
00975 int * mpar;
00976 float * xdat, * wdat;
00977
00978 if ( line_image == NULL )
00979 {
00980 cpl_msg_error ("fitLines:"," no image given\n") ;
00981 return -18 ;
00982 }
00983 if ( n_lines == NULL )
00984 {
00985 cpl_msg_error ("fitLines:"," no counter of emission lines\n") ;
00986 return -19 ;
00987 }
00988 if ( row == NULL || width <= 0 )
00989 {
00990 cpl_msg_error ("fitLines:"," row or width vectors are empty\n") ;
00991 return -20 ;
00992 }
00993 if ( wavelength == NULL )
00994 {
00995 cpl_msg_error ("fitLines:"," no wavelength array given\n") ;
00996 return -21 ;
00997 }
00998
00999 k = 0 ;
01000
01001 /* allocate memory for the spectral vector */
01002 line = newVector (2*width + 1) ;
01003 /* allocate memory */
01004 xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
01005 wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
01006 mpar = (int *) cpl_calloc( MAXPAR, sizeof (int) ) ;
01007
01008 /* go through the columns */
01009 for ( i = 0 ; i < line_image -> lx ; i++ )
01010 {
01011 if ( n_lines[i] == 0 )
01012 {
01013 continue ;
01014 }
01015 /* go through the emission lines in a column */
01016 for ( l = 0 ; l < n_lines[i] ; l++ )
01017 {
01018 if ( row[i][l] <= 0 )
01019 {
01020 continue ;
01021 }
01022
01023 /* --------------------------------------------------------------
01024 * fit the single lines using linefit and store the parameters in
01025 * an array of the FitParams data structure allParams[].
01026 */
01027 if ( (result = linefit ( line_image, allParams[k], fwhm, l, i,
01028 width, row[i][l], min_amplitude,line,mpar,xdat,wdat ) ) < 0 )
01029 {
01030 cpl_msg_debug ("fitLines:"," linefit failed, error no.: %d, column: %d, row: %d, line: %d\n", result, i, row[i][l], l) ;
01031 continue ;
01032 }
01033 if ( (allParams[k] -> fit_par[0] <= 0.) || (allParams[k] -> fit_par[1] <= 0.)
01034 || (allParams[k] -> fit_par[2] <= 0.) )
01035 {
01036 cpl_msg_warning ("fitLines:"," negative fit parameters in column: %d, line: %d\n", i, l) ;
01037 continue ;
01038 }
01039 allParams[k] -> wavelength = wavelength[i][l] ;
01040 k++ ;
01041 }
01042 }
01043
01044 /* free memory */
01045 destroyVector(line);
01046 cpl_free(xdat);
01047 cpl_free(wdat);
01048 cpl_free(mpar);
01049
01050 /* all is o.k. */
01051 return k ;
01052 }
01053
01054 /*----------------------------------------------------------------------------
01055 Function : polyfit()
01056 In : par: filled array of fit parameter structure
01057 column: image column index
01058 n_lines: number of found lines in column
01059 n_rows: number of image rows
01060 dispersion: microns per pixel
01061 max_residual: maximum residual value, beyond that value
01062 the fit is rejected.
01063 acoefs: array of the 3 coefficients of the fitted
01064 parabola
01065 dacoefs: array of standard deviations of the 3 coefficients
01066 n_reject: rejected number of fits due to high residuals
01067 n_fitcoefs: number of polynomial coefficients to fit
01068 Out : chisq, the three fit coefficients acoefs[i] and their
01069 standard deviations dacoefs[i], the rejected number of fits
01070 due to too high residuals: n_reject
01071 FLT_MAX in case of error
01072 Job : fits a polynomial
01073 lambda[i] = a1 + a2*pos[i] + a3*pos[i]^2 +...
01074 to determine the connection between the listed wave-
01075 length values and the gauss-fitted positions for each
01076 image column using the singular value decomposition
01077 method.
01078 ---------------------------------------------------------------------------*/
01079
01080 float polyfit( FitParams ** par,
01081 int column,
01082 int n_lines,
01083 int n_rows,
01084 float dispersion,
01085 float max_residual,
01086 float * acoefs,
01087 float * dacoefs,
01088 int * n_reject,
01089 int n_fitcoefs )
01090 {
01091 float ** ucoefs, ** vcoefs, ** covar ;
01092 float *mem;
01093 float wcoefs[n_fitcoefs] ;
01094 float * lambda, * posit ;
01095 float * weight, * resid ;
01096 float * newlam, * newpos, * newwet ;
01097 float chisq, result ;
01098 float offset ;
01099 int num, found ;
01100 int i, j, k, n ;
01101
01102 /* reset the fit coefficients and their errors */
01103 for ( i = 0 ; i < n_fitcoefs ; i++ )
01104 {
01105 acoefs[i] = 0. ;
01106 dacoefs[i] = 0. ;
01107 }
01108 if ( NULL == par )
01109 {
01110 cpl_msg_error("polyfit:"," no fit params given\n");
01111 return FLT_MAX ;
01112 }
01113
01114 if ( 0 >= n_lines )
01115 {
01116 /*
01117 cpl_msg_warning ("polyfit:"," sorry, number of lines is wrong") ;
01118 */
01119 return FLT_MAX ;
01120 }
01121 if ( 0 >= n_rows )
01122 {
01123 cpl_msg_error ("polyfit:"," sorry, number of rows is wrong") ;
01124 return FLT_MAX ;
01125 }
01126 if ( dispersion == 0. )
01127 {
01128 cpl_msg_error ("polyfit:"," sorry, wrong dispersion given") ;
01129 return FLT_MAX ;
01130 }
01131
01132 offset = (float)(n_rows - 1)/2. ;
01133
01134 /* allocate memory */
01135
01136 mem = (float*) cpl_calloc( n_lines*7, sizeof (float) ) ;
01137 lambda = mem;
01138 posit = mem + n_lines;
01139 weight = mem + n_lines*2;
01140 resid = mem + n_lines*3;
01141 newlam = mem + n_lines*4;
01142 newpos = mem + n_lines*5;
01143 newwet = mem + n_lines*6;
01144
01145 /*lambda = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01146 posit = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01147 weight = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01148 resid = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01149 newlam = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01150 newpos = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01151 newwet = (float*) cpl_calloc( n_lines, sizeof (float) ) ;*/
01152
01153 /* allocate coefficient matrices*/
01154 ucoefs = matrix ( 1, n_lines, 1, n_fitcoefs ) ;
01155 vcoefs = matrix ( 1, n_lines, 1, n_fitcoefs ) ;
01156 covar = matrix ( 1, n_fitcoefs, 1, n_fitcoefs ) ;
01157
01158 /* go through all fit parameters */
01159 n = 0 ;
01160 for ( i = 0 ; i < (par[0] -> n_params) ; i ++ )
01161 {
01162 found = -1 ;
01163 /* find the given column and go through the lines in that column */
01164 for ( j = 0 ; j < n_lines ; j ++ )
01165 {
01166 if ( (par[i] -> column == column) && (par[i] -> line == j) )
01167 {
01168 found = i ;
01169 }
01170 else
01171 {
01172 continue ;
01173 }
01174
01175 /* store only fit params with reasonable values */
01176 if ( par[found] -> derv_par[2] != 0. && par[found] -> fit_par[2] > 0. &&
01177 par[found] -> wavelength > 0. && par[found] -> fit_par[1] > 0. &&
01178 par[found] -> fit_par[0] > 0. )
01179 {
01180 /* ---------------------------------------------------------------------------
01181 * store the found position, error of the position as weight and the associated
01182 * wavelength values of the fitted lines
01183 */
01184 posit[n] = par[found] -> fit_par[2] ;
01185 weight[n] = par[found] -> derv_par[2] ;
01186 lambda[n] = par[found] -> wavelength ;
01187 n ++ ;
01188 }
01189 else
01190 {
01191 continue ;
01192 }
01193 }
01194
01195 }
01196
01197 num = n ;
01198 if ( num < n_fitcoefs )
01199 {
01200 cpl_msg_warning("polyfit:"," not enough lines found in column %d to determine the three coefficients.\n", column) ;
01201 for ( i = 0 ; i < n_fitcoefs ; i++ )
01202 {
01203 acoefs[i] = ZERO ;
01204 dacoefs[i] = ZERO ;
01205 }
01206 free_matrix ( ucoefs, 1/*, n_lines*/, 1/*, n_fitcoefs*/ ) ;
01207 free_matrix ( vcoefs, 1/*, n_lines*/, 1/*, n_fitcoefs*/ ) ;
01208 free_matrix ( covar, 1/*, n_fitcoefs*/, 1/*, n_fitcoefs*/ ) ;
01209 /*cpl_free (lambda) ;
01210 cpl_free (posit) ;
01211 cpl_free (weight) ;
01212 cpl_free (resid) ;
01213 cpl_free (newlam) ;
01214 cpl_free (newpos) ;
01215 cpl_free (newwet) ;*/
01216 cpl_free (mem);
01217 return FLT_MAX ;
01218 }
01219
01220 /*--------------------------------------------------------------------------------
01221 * scale the pixel position values to smaller than 1 and transform the weights to
01222 * wavelength units
01223 */
01224
01225 for ( i = 0 ; i < num ; i ++ )
01226 {
01227 posit[i] = (posit[i] - offset)/offset ;
01228 weight[i] *= fabs(dispersion) ;
01229 }
01230
01231 /* do the fit using the singular value decomposition method */
01232 svd_fitting( posit - 1, lambda - 1, weight - 1, num, acoefs-1, n_fitcoefs,
01233 ucoefs, vcoefs, wcoefs-1, covar, &chisq, fpol ) ;
01234
01235 /* scale the linear and the quadratic coefficient */
01236 for ( i = 1 ; i < n_fitcoefs ; i++ )
01237 {
01238 acoefs[i] /= pow(offset, i) ;
01239 }
01240
01241 /* now that we have determined the fit coefficients, find the residuals */
01242 *n_reject = 0 ;
01243
01244 j = 0 ;
01245 for ( i = 0 ; i < num ; i++ )
01246 {
01247 result = 0. ;
01248 for ( k = 0 ; k < n_fitcoefs ; k++ )
01249 {
01250 result += acoefs[k] * pow(posit[i], k) ;
01251 }
01252
01253 resid[i] = lambda[i] - result ;
01254
01255 if ( fabs( resid[i] ) > max_residual)
01256 {
01257 (*n_reject) ++ ;
01258 }
01259 else
01260 {
01261 newlam[j] = lambda[i] ;
01262 newpos[j] = posit[i] ;
01263 newwet[j] = weight[i] ;
01264 j++ ;
01265 }
01266 }
01267
01268 num = j ;
01269 if ( num >= n_fitcoefs )
01270 {
01271 svd_fitting( newpos - 1, newlam - 1, newwet - 1, num, acoefs-1, n_fitcoefs, ucoefs,
01272 vcoefs, wcoefs-1, covar, &chisq, fpol ) ;
01273
01274 /* scale the resulting coefficients */
01275 for ( i = 0 ; i < n_fitcoefs ; i++ )
01276 {
01277 acoefs[i] /= pow(offset, i) ;
01278 dacoefs[i] = sqrt( (double) covar[i+1][i+1] ) / pow(offset, i) ;
01279 }
01280 }
01281 else
01282 {
01283 cpl_msg_warning ("polyfit:"," too many lines rejected (number: %d) due to high residuals, fit coefficients are set zero, in column: %d\n", *n_reject, column) ;
01284 for ( i = 0 ; i < n_fitcoefs ; i++ )
01285 {
01286 acoefs[i] = ZERO ;
01287 dacoefs[i] = ZERO ;
01288 }
01289 }
01290
01291 free_matrix ( ucoefs, 1/*, n_lines*/, 1/*, n_fitcoefs*/ ) ;
01292 free_matrix ( vcoefs, 1/*, n_lines*/, 1/*, n_fitcoefs*/ ) ;
01293 free_matrix ( covar, 1/*, n_fitcoefs*/, 1/*, n_fitcoefs*/ ) ;
01294 /*cpl_free (lambda) ;
01295 cpl_free (posit) ;
01296 cpl_free (weight) ;
01297 cpl_free (resid) ;
01298 cpl_free (newlam) ;
01299 cpl_free (newpos) ;
01300 cpl_free (newwet) ;*/
01301 cpl_free (mem);
01302
01303 return chisq ;
01304 }
01305
01306 /*----------------------------------------------------------------------------
01307 Function : coefsCrossFit()
01308 In : n_columns: number of image columns
01309 acoefs: coeffs fitted in polyfit
01310 note: this is a vector of coefficients with
01311 the same index for all columns
01312 dacoefs: fit errors of the corresponding acoefs
01313 bcoefs: the fitted coefs
01314 n_fitcoefs: number of fit coefficients
01315 sigma_factor: factor of sigma beyond which the
01316 column coefficients are discarded for the fit
01317 Out : chisq, the found fit coefficients
01318 Job : Fits each single polnomial coefficient acoefs resulting from polyfit
01319 across the image columns
01320 ---------------------------------------------------------------------------*/
01321
01322 float coefsCrossFit ( int n_columns,
01323 float * acoefs,
01324 float * dacoefs,
01325 float * bcoefs,
01326 int n_fitcoefs,
01327 float sigma_factor )
01328 {
01329 float col_index, sub_col_index[n_columns] ;
01330 float sub_acoefs[n_columns], sub_dacoefs[n_columns] ;
01331 float wcoefs[n_fitcoefs] ;
01332 float ** ucoefs, **vcoefs, **covar ;
01333 float chisq ;
01334 float * acoefsclean ;
01335 double sum, sumq, mean ;
01336 double sigma ;
01337 double cliphi, cliplo ;
01338 float offset ;
01339 int i, n, num, ndata ;
01340 int nc ;
01341
01342 if ( n_columns < 1 )
01343 {
01344 cpl_msg_error("coefsCrossFit:"," wrong number of image columns given\n") ;
01345 return FLT_MAX ;
01346 }
01347 if ( acoefs == NULL || dacoefs == NULL )
01348 {
01349 cpl_msg_error("coefsCrossFit:"," coeffs or errors of coefficients are not given\n") ;
01350 return FLT_MAX ;
01351 }
01352 if ( bcoefs == NULL )
01353 {
01354 cpl_msg_error("coefsCrossFit:"," coeffs are not allocated\n") ;
01355 return FLT_MAX ;
01356 }
01357
01358 if ( n_fitcoefs < 1 )
01359 {
01360 cpl_msg_error("coefsCrossFit:"," wrong number of fit coefficients\n") ;
01361 return FLT_MAX ;
01362 }
01363 if ( sigma_factor <= 0. )
01364 {
01365 cpl_msg_error("coefsCrossFit:"," impossible sigma_factor given!\n") ;
01366 return FLT_MAX ;
01367 }
01368
01369 offset = (float)(n_columns - 1) / 2. ;
01370
01371 /* ----------------------------------------------------------
01372 * determine the clean mean and sigma value of the coefficients,
01373 * that means reject 10 % of the extreme low and high values
01374 */
01375 nc = 0 ;
01376 for ( i = 0 ; i < n_columns ; i++ )
01377 {
01378 if ( qfits_isnan(acoefs[i]) || acoefs[i] == 0. || dacoefs[i] == 0. )
01379 {
01380 continue ;
01381 }
01382 else
01383 {
01384 nc++ ;
01385 }
01386 }
01387 acoefsclean = (float*) cpl_calloc(nc , sizeof(float)) ;
01388 nc = 0 ;
01389 for ( i = 0 ; i < n_columns ; i++ )
01390 {
01391 if ( qfits_isnan(acoefs[i]) || acoefs[i] == 0. || dacoefs[i] == 0. )
01392 {
01393 continue ;
01394 }
01395 else
01396 {
01397 acoefsclean[nc] = acoefs[i] ;
01398 nc++ ;
01399 }
01400 }
01401 pixel_qsort(acoefsclean, nc) ;
01402 sum = 0. ;
01403 sumq = 0. ;
01404 mean = 0. ;
01405 sigma = 0. ;
01406 n = 0 ;
01407 for ( i = (int)((float)nc*LOW_REJECT) ; i < (int)((float)nc*HIGH_REJECT) ; i++ )
01408 {
01409 sum += (double)acoefsclean[i] ;
01410 sumq += ((double)acoefsclean[i] * (double)acoefsclean[i]) ;
01411 n ++ ;
01412 }
01413 mean = sum/(double)n ;
01414 sigma = sqrt( sumq/(double)n - (mean * mean) ) ;
01415 cliphi = mean + sigma * (double)sigma_factor ;
01416 cliplo = mean - sigma * (double)sigma_factor ;
01417
01418 /* fit only the reasonnable values */
01419 num = 0 ;
01420 for ( i = 0 ; i < n_columns ; i++ )
01421 {
01422 /* associate the column indices to the corresponding array */
01423 col_index = (float) i ;
01424
01425 /* take only the reasonnable coefficients */
01426 if ( !qfits_isnan(acoefs[i]) && (acoefs[i] <= cliphi) && (acoefs[i] >= cliplo) &&
01427 (dacoefs[i] != 0. ) && (acoefs[i] != 0.) )
01428 {
01429 sub_acoefs[num] = acoefs[i] ;
01430 sub_dacoefs[num] = dacoefs[i] ;
01431 sub_col_index[num] = col_index ;
01432 num ++ ;
01433 }
01434 }
01435 ndata = num ;
01436
01437 if ( ndata < n_fitcoefs )
01438 {
01439 cpl_msg_error("coefsCrossFit:"," not enough data found to determine the fit coefficients.\n") ;
01440
01441 return FLT_MAX ;
01442 }
01443
01444 /* allocate coefficient matrices */
01445 ucoefs = matrix(1, ndata, 1, n_fitcoefs) ;
01446 vcoefs = matrix(1, ndata, 1, n_fitcoefs) ;
01447 covar = matrix ( 1, n_fitcoefs, 1, n_fitcoefs ) ;
01448
01449 /* scale the x-values for the fit */
01450 for ( i = 0 ; i < ndata ; i++ )
01451 {
01452 sub_col_index[i] = (sub_col_index[i] - offset) / offset ;
01453 }
01454
01455 /* finally, do the singular value decomposition fit */
01456 svd_fitting ( sub_col_index-1, sub_acoefs-1, sub_dacoefs-1, ndata, bcoefs-1,
01457 n_fitcoefs, ucoefs, vcoefs, wcoefs-1, covar, &chisq, fpol ) ;
01458
01459 /* scale the found coefficients */
01460 for ( i = 0 ; i < n_fitcoefs ; i ++ )
01461 {
01462 bcoefs[i] /= pow(offset, i) ;
01463 }
01464
01465 /* free memory */
01466 cpl_free (acoefsclean) ;
01467 free_matrix( ucoefs, 1/*, ndata*/, 1/*, n_fitcoefs */) ;
01468 free_matrix( vcoefs, 1/*, ndata*/, 1/*, n_fitcoefs */) ;
01469 free_matrix ( covar, 1/*, n_fitcoefs*/, 1/*, n_fitcoefs*/ ) ;
01470
01471 return chisq ;
01472 }
01473
01474
01475 /*----------------------------------------------------------------------------
01476 Function : waveMap()
01477 In : lineImage: image from a calibration emission lamp,
01478 bcoefs: transformed fit coefficients
01479 n_a_fitcoefs: number of fit coefficients for the single
01480 column fits lambda-position
01481 n_b_fitcoefs: number of fit coefficients for the fits of
01482 the single a coefficients across the columns
01483 wavelength: wavelength list from lamp file
01484 intensity: corresponding line intensity from line list
01485 n_lines: number of lines in the list
01486 magFactor: magnifying factor of the image for FFT
01487 Out : wavelength calibration map image.
01488 Job : this routine determines a wavelength calibration map
01489 frame associating a wavelength value to each pixel
01490 by using the fit coefficients determined before.
01491 ---------------------------------------------------------------------------*/
01492
01493 OneImage * waveMap( OneImage * lineImage,
01494 float ** bcoefs,
01495 int n_a_fitcoefs,
01496 int n_b_fitcoefs,
01497 float * wavelength,
01498 float * intensity,
01499 int n_lines,
01500 int magFactor)
01501 {
01502 OneImage * retImage ;
01503 float cenpos, cenpix ;
01504 float centreval, centrepix, wavelag ;
01505 float pixvalue ;
01506 double a[n_a_fitcoefs], wave[n_lines] ;
01507 float a_initial ;
01508 int i, j, k, l/*, m*/, line, col, row, found, sign ;
01509 int var, maxlag, cmin, cmax, offset ;
01510 float emline[2*magFactor*lineImage -> ly], spec[2*magFactor*lineImage -> ly] ;
01511 double * result ;
01512 float col_off ;
01513 float angst ;
01514 double xcorr_max ;
01515 int delta ;
01516 double z[2*(n_a_fitcoefs - 1)] ;
01517 gsl_poly_complex_workspace * w ;
01518
01519 if ( NullImage == lineImage )
01520 {
01521 cpl_msg_error("waveMap:","no image given\n") ;
01522 return NullImage ;
01523 }
01524
01525 if ( NULL == wavelength || n_lines <= 0 )
01526 {
01527 cpl_msg_error("waveMap:","no wavelength list given\n") ;
01528 return NullImage ;
01529 }
01530
01531 if ( NULL == intensity )
01532 {
01533 cpl_msg_error("waveMap:","no intensity list given\n") ;
01534 return NullImage ;
01535 }
01536
01537 if ( NULL == bcoefs )
01538 {
01539 cpl_msg_error("waveMap:","no coefficients given\n") ;
01540 return NullImage ;
01541 }
01542
01543 if ( magFactor <= 1 )
01544 {
01545 cpl_msg_error("waveMap:","wrong magnifying factor given\n") ;
01546 return NullImage ;
01547 }
01548
01549 /* allocate memory */
01550 if ( NULL == ( retImage = new_image( lineImage -> lx, lineImage -> ly ) ))
01551 {
01552 cpl_msg_error("waveMap:","cannot allocate a new image\n");
01553 return NullImage ;
01554 }
01555
01556 var = (magFactor - 1)*(magFactor - 1) ;
01557 offset = lineImage -> ly * (magFactor/4 + 1) ;
01558
01559 /* find out if Angstroem or microns are used */
01560 if ( wavelength[0] > 10000. )
01561 {
01562 /* Angstroem */
01563 angst = 10000. ;
01564 }
01565 else if ( wavelength[0] > 1000. && wavelength[0] < 10000. )
01566 {
01567 /* nanometers */
01568 angst = 1000. ;
01569 }
01570 else
01571 {
01572 /* microns */
01573 angst = 1. ;
01574 }
01575
01576 /* go through the image columns */
01577 for ( col = 0 ; col < lineImage -> lx ; col++ )
01578 {
01579 /* initialize the emline array for each column */
01580 for ( i = 0 ; i < 2*magFactor*lineImage -> ly ; i++ )
01581 {
01582 emline[i] = 0. ;
01583 }
01584 col_off = (float)col - (float)(lineImage->lx-1)/2. ;
01585
01586 /* determine the coefficients by using the given bcoefs */
01587 for ( i = 0 ; i < n_a_fitcoefs ; i++ )
01588 {
01589 /* initialize coefficients and solution */
01590 a[i] = 0. ;
01591 if (i < n_a_fitcoefs-1)
01592 {
01593 z[2*i] = 0. ;
01594 z[2*i+1] = 0. ;
01595 }
01596 for ( j = 0 ; j < n_b_fitcoefs ; j++ )
01597 {
01598 a[i] += bcoefs[i][j] * pow(col_off, j) ;
01599 }
01600 }
01601 a_initial = a[0] ;
01602
01603 /* go through the lines and generate an artificial spectrum */
01604 for ( line = 0 ; line < n_lines ; line++ )
01605 {
01606 /* go from Angstroem to micron */
01607 wave[line] = wavelength[line]/angst ;
01608
01609 /* ---------------------------------------------------------------------
01610 * solve the polynomial for the exact offset of the line that means
01611 * find the root of the polynomial of order n_fitcoefs - 1
01612 */
01613 a[0] = a_initial - wave[line] ;
01614
01615 if (NULL == (w = gsl_poly_complex_workspace_alloc(n_a_fitcoefs)))
01616 {
01617 cpl_msg_error("waveMap:","could not allocate complex workspace!") ;
01618 destroy_image(retImage) ;
01619 return NullImage ;
01620 }
01621 if (-1 == gsl_poly_complex_solve(a, n_a_fitcoefs, w, z))
01622 {
01623 cpl_msg_error("waveMap:","gsl_poly_complex_solve did not work!") ;
01624 destroy_image(retImage) ;
01625 return NullImage ;
01626 }
01627 gsl_poly_complex_workspace_free(w) ;
01628
01629
01630 j = 0 ;
01631 found = -1 ;
01632 for ( i = 0 ; i < n_a_fitcoefs - 1 ; i++ )
01633 {
01634 /* test for appropriate solution */
01635 if( z[2*i] > (-1.)*(float) lineImage->ly/2. &&
01636 z[2*i] < (float)lineImage -> ly/2. && z[2*i+1] == 0. )
01637 {
01638 found = 2*i ;
01639 j ++ ;
01640 }
01641 else
01642 {
01643 continue ;
01644 }
01645 }
01646 if ( j == 0 )
01647 {
01648 cpl_msg_warning("waveMap:","no offset solution found for line %d in column %d\n", line, col) ;
01649 continue ;
01650 }
01651 else if ( j == 1 )
01652 {
01653 cenpos = z[found] + (float) lineImage->ly /2. ;
01654 }
01655 else
01656 {
01657 cpl_msg_warning("waveMap:","two or more offset solutions found for line %d in column %d\n", line, col) ;
01658 continue ;
01659 }
01660
01661 /*---------------------------------------------------------------------------------
01662 * magnify image by the given factor add an additional offset
01663 */
01664 cenpix = cenpos * (float) magFactor + (float) offset ;
01665
01666 /* determine max and min pixel limits over which line should be convolved */
01667 cmin = (nint(cenpix) - (var-1)) > 0 ? nint(cenpix) - (var-1) : 0 ;
01668 cmax = (nint(cenpix) + (var-1)) < 2*magFactor*lineImage -> ly ?
01669 nint(cenpix) + (var-1) : 2*magFactor*lineImage -> ly ;
01670
01671 /* convolve neon lines with Gaussian function */
01672 for ( j = cmin ; j < cmax ; j++ )
01673 {
01674 emline[j] += intensity[line] * exp((double)(-0.5*(j-cenpix)*(j-cenpix))/(double)var) ;
01675 }
01676 }
01677
01678 /*-----------------------------------------------------------------------------------------
01679 * for each column, map the image data points onto an magFactor times bigger element grid for FFT
01680 * in the cross correlation, first initialize the two helping arrays for each new column.
01681 */
01682 for ( k = 0 ; k < 2*magFactor*lineImage -> ly ; k++ )
01683 {
01684 spec[k] = 0. ;
01685 }
01686
01687 /* now take the image data points of the column and put them into the spec array */
01688 for ( row = 0 ; row < lineImage -> ly ; row++ ) /* go through the column */
01689 {
01690 /* insert 8 values for each image row (magnification) and add same offset as for emline array */
01691 for ( l = 0 ; l < magFactor ; l++ )
01692 {
01693 /* set bad pixels or negative values to zero */
01694 if (!qfits_isnan(lineImage -> data[col + row*lineImage -> lx]) &&
01695 (lineImage -> data[col + row*lineImage -> lx] > 0.))
01696 {
01697 spec[offset + l + (row * magFactor)] =
01698 lineImage -> data[col + row*lineImage -> lx] ;
01699 }
01700 else
01701 {
01702 spec[offset + l + (row * magFactor)] = 0. ;
01703 }
01704 }
01705 }
01706
01707 /* now call the cross correlation routine */
01708 if (NULL == (result = xcorrel(spec, 2*magFactor*lineImage->ly, emline, 2*magFactor*lineImage -> ly,
01709 magFactor*lineImage->ly, &delta, &maxlag, &xcorr_max)) )
01710 {
01711 cpl_msg_warning ("cross-correl","cross correlation did not work, col: %d is set ZERO\n", col) ;
01712 for ( row = 0 ; row < lineImage -> ly ; row++ )
01713 {
01714 retImage -> data[col + row*lineImage -> lx] = ZERO ;
01715 }
01716 continue ;
01717 }
01718
01719 if ( xcorr_max <= 0. )
01720 {
01721 cpl_msg_warning ("cross-correl","cross correlation sum is negative, col: %d is set ZERO\n", col) ;
01722 for ( row = 0 ; row < lineImage -> ly ; row++ )
01723 {
01724 retImage -> data[col + row*lineImage -> lx] = ZERO ;
01725 }
01726 cpl_free(result) ;
01727 continue ;
01728 }
01729
01730 wavelag = (float) -delta / (float) magFactor ;
01731 if ( fabs(wavelag) > (float)lineImage->ly/20. )
01732 {
01733 cpl_msg_warning ("cross-correl","wave lag too big, col: %d is set ZERO\n", col) ;
01734 for ( row = 0 ; row < lineImage -> ly ; row++ )
01735 {
01736 retImage -> data[col + row*lineImage -> lx] = ZERO ;
01737 }
01738 cpl_free(result) ;
01739 continue ;
01740 }
01741
01742 /*-----------------------------------------------------------------------------
01743 * determine new zero order coefficient centreval, of which the formula is
01744 * determined by setting equal a polynomial shifted by wavelag with the same
01745 * higher order coefficients and set the new zero order coefficient to
01746 * get both sides of the equation approximately equal.
01747 */
01748 centreval = a_initial ;
01749 for ( i = 1 ; i < n_a_fitcoefs ; i++ )
01750 {
01751 if ( i%2 == 0 )
01752 {
01753 sign = -1 ;
01754 }
01755 else
01756 {
01757 sign = 1 ;
01758 }
01759 centreval += (float)sign * a[i]*pow(wavelag, i) ;
01760 }
01761
01762 /* prepare to write out wavelength as pixel values */
01763 for ( row = 0 ; row < lineImage -> ly ; row++ )
01764 {
01765 centrepix = (float)row - ((float)lineImage->ly - 1.)/2. ;
01766 pixvalue = 0. ;
01767 for ( i = 1 ; i < n_a_fitcoefs ; i++ )
01768 {
01769 pixvalue += a[i]*pow(centrepix, i) ;
01770 }
01771 retImage -> data[col + row*lineImage -> lx] =
01772 centreval + pixvalue ;
01773 }
01774 cpl_free(result) ;
01775 }
01776 return retImage ;
01777 }
01778
01779 /*----------------------------------------------------------------------------
01780 Function : wavelengthCalibration()
01781 In : image: merged image from a calibration emission lamp,
01782 wave: wavelength array read from the wavelength list
01783 n_lines: number of lines in the wavelength list
01784 row_clean: resulting list of the row indices but without the
01785 lines that are too close to each other for the fit
01786 output of findLines()
01787 wavelength_clean: corrected wavelength list corresponding to
01788 the row_clean array
01789 output of findLines()
01790 n_found_lines: output of findLines(): total number of found emission lines
01791 dispersion: dispersion of spectrum: micron per pixel
01792 halfWidth: half width of the box where the line must sit
01793 minAmplitude: minimum amplitude of the Gaussian to do the fit
01794 max_residual: maximum residual value, beyond that value
01795 the polynomial lambda-position fit is rejected.
01796 fwhm: first guess for the full width of half maximum
01797 of the gaussian line fit
01798 n_a_fitcoefs: number of fit coefficients for the single
01799 column fits: lambda-position
01800 n_b_fitcoefs: number of fit coefficients for the fits of
01801 the single a coefficients across the columns
01802 sigmaFactor: factor of the standard deviation of the determined
01803 polynomial coefficients of the columns beyond
01804 which these coefficients are not used to carry out
01805 the polynomial fit across the columns.
01806 pixel_tolerance: maximum tolerated difference between estimated
01807 and fitted line positions.
01808 Out : 0 if all went o.k., -1 if something went wrong.
01809 bcoefs: array of cooefficients of the polynomial fit across the columns.
01810 par: array of the resulting FitParams data structure
01811 Job : this routine takes an image from a calibration
01812 emission lamp and delivers the fit coefficients of
01813 a polynomial fit across the columns
01814 of the coefficients of the polynomial line position
01815 fits as output. Furthermore it delivers an array of the fit parameters
01816 as output. This routine expects Nyquist sampled spectra
01817 (either an interleaved image or an image convolved with an
01818 appropriate function in spectral direction)
01819 ---------------------------------------------------------------------------*/
01820
01821 int wavelengthCalibration( OneImage * image,
01822 FitParams ** par ,
01823 float ** bcoefs,
01824 float * wave,
01825 int n_lines,
01826 int ** row_clean,
01827 float ** wavelength_clean,
01828 int * n_found_lines,
01829 float dispersion,
01830 int halfWidth,
01831 float minAmplitude,
01832 float max_residual,
01833 float fwhm,
01834 int n_a_fitcoefs,
01835 int n_b_fitcoefs,
01836 float sigmaFactor,
01837 float pixel_tolerance )
01838
01839 {
01840 int i, j, k ;
01841 int n_fit ;
01842 int n_reject ;
01843 float * acoefs ;
01844 float * dacoefs ;
01845 float ** abuf ;
01846 float ** dabuf ;
01847 float chisq_poly, chisq_cross ;
01848 int zeroind ;
01849 /*float * mem ;*/
01850
01851
01852 if ( NullImage == image )
01853 {
01854 cpl_msg_error("wavelengthCalibration:","no image given\n") ;
01855 return -1 ;
01856 }
01857 if ( par == NULL )
01858 {
01859 cpl_msg_error("wavelengthCalibration:","no fit parameter data structure given\n") ;
01860 return -1 ;
01861 }
01862 if ( wave == NULL )
01863 {
01864 cpl_msg_error("wavelengthCalibration:","no wavelength list given\n") ;
01865 return -1 ;
01866 }
01867 if ( n_lines <= 0 )
01868 {
01869 cpl_msg_error("wavelengthCalibration:","impossible number of lines in line list given\n") ;
01870 return -1 ;
01871 }
01872 if ( row_clean == NULL )
01873 {
01874 cpl_msg_error("wavelengthCalibration:","no row_clean array given\n") ;
01875 return -1 ;
01876 }
01877 if ( wavelength_clean == NULL )
01878 {
01879 cpl_msg_error("wavelengthCalibration:","no wavelength_clean array given\n") ;
01880 return -1 ;
01881 }
01882
01883 if ( dispersion == 0. )
01884 {
01885 cpl_msg_error("wavelengthCalibration:","impossible dispersion given\n") ;
01886 return -1 ;
01887 }
01888
01889 if ( halfWidth <= 0 || halfWidth > image->ly/2 )
01890 {
01891 cpl_msg_error("wavelengthCalibration:","impossible half width of the fitting box given\n") ;
01892 return -1 ;
01893 }
01894 if ( minAmplitude < 1. )
01895 {
01896 cpl_msg_error("wavelengthCalibration:","impossible minimal amplitude\n") ;
01897 return -1 ;
01898 }
01899
01900 if ( max_residual <= 0. || max_residual > 1. )
01901 {
01902 cpl_msg_error("wavelengthCalibration:","impossible max_residual given\n") ;
01903 return -1 ;
01904 }
01905
01906 if ( fwhm <= 0. || fwhm > 10. )
01907 {
01908 cpl_msg_error("wavelengthCalibration:","impossible fwhm given\n") ;
01909
01910 return -1 ;
01911 }
01912
01913 if ( n_a_fitcoefs <= 0 || n_a_fitcoefs > 9 )
01914 {
01915 cpl_msg_error("wavelengthCalibration:","unrealistic n_a_fitcoefs given\n") ;
01916 return -1 ;
01917 }
01918
01919 if ( n_b_fitcoefs <= 0 || n_b_fitcoefs > 9 )
01920 {
01921 cpl_msg_error("wavelengthCalibration:"," unrealistic n_b_fitcoefs given\n") ;
01922 return -1 ;
01923 }
01924 if ( sigmaFactor <= 0. )
01925 {
01926 cpl_msg_error("wavelengthCalibration:"," impossible sigmaFactor given\n") ;
01927 return -1 ;
01928 }
01929
01930 /* initialize the variables */
01931 n_reject = 0 ;
01932 n_fit = 0 ;
01933
01934 /* fit each found line by using a gaussian function and determine the exact position */
01935 if ( 0 > (n_fit = fitLines( image , par, fwhm, n_found_lines, row_clean, wavelength_clean,
01936 halfWidth, minAmplitude )) )
01937 {
01938 cpl_msg_error("wavelengthCalibration:"," cannot fit the lines, error code of fitLines: %d\n", n_fit) ;
01939 return -1 ;
01940 }
01941
01942 /* first check for faked lines like bad pixels */
01943 if ( -1 == checkForFakeLines (par, dispersion, wavelength_clean, row_clean, n_found_lines,
01944 image->lx, pixel_tolerance) )
01945 {
01946 cpl_msg_error("waveCal:"," cannot fit the lines, error code of fitLines: %d\n", n_fit) ;
01947 return -1 ;
01948 }
01949
01950 /* allocate memory */
01951 if ( NULL == (acoefs = (float*) cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
01952 NULL == (dacoefs = (float*) cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
01953 NULL == (abuf = (float**) cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
01954 NULL == (dabuf = (float**) cpl_calloc (n_a_fitcoefs, sizeof(float*))) )
01955 {
01956 cpl_msg_error("wavelengthCalibration:"," cannot allocate memory\n") ;
01957 return -1 ;
01958 }
01959
01960 for ( i = 0 ; i < n_a_fitcoefs ; i++ )
01961 {
01962 if ( NULL == (abuf[i] = (float*) cpl_calloc(image -> lx, sizeof(float))) ||
01963 NULL == (dabuf[i] = (float*) cpl_calloc(image -> lx, sizeof(float))) )
01964 {
01965 cpl_msg_error("wavelengthCalibration:"," cannot allocate memory\n") ;
01966 cpl_free(abuf) ;
01967 cpl_free(dabuf) ;
01968 return -1 ;
01969 }
01970 }
01971
01972 /* fit wavelengths to the corresponding found positions for each column */
01973 k = 0 ;
01974
01975 for ( i = 0 ; i < image -> lx ; i++ )
01976 {
01977 zeroind = 0 ;
01978 if ( FLT_MAX == (chisq_poly = polyfit( par, i, n_found_lines[i], image -> ly, dispersion,
01979 max_residual, acoefs, dacoefs, &n_reject, n_a_fitcoefs)) )
01980 {
01981 /*
01982 cpl_msg_warning ("wavelengthCalibration:"," error in polyfitt in column: %d\n", i) ;
01983 */
01984 for ( j = 0 ; j < n_a_fitcoefs ; j++ )
01985 {
01986 acoefs[j] = ZERO ;
01987 dacoefs[j] = ZERO ;
01988 }
01989 }
01990
01991 for ( j = 0 ; j < n_a_fitcoefs ; j++ )
01992 {
01993 if ( acoefs[0] <= 0. || acoefs[1] ==0. ||
01994 dacoefs[j] == 0. || qfits_isnan(acoefs[j]) )
01995 {
01996 zeroind = 1 ;
01997
01998 }
01999 }
02000 for ( j = 0 ; j < n_a_fitcoefs ; j++ )
02001 {
02002 if ( zeroind == 0 )
02003 {
02004 abuf[j][i] = acoefs[j] ;
02005 dabuf[j][i] = dacoefs[j] ;
02006 }
02007 else
02008 {
02009 abuf[j][i] = ZERO ;
02010 dabuf[j][i] = ZERO ;
02011 }
02012 }
02013 }
02014
02015 /* fit each acoefs across the columns to smooth the result */
02016 for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02017 {
02018 if ( FLT_MAX == (chisq_cross = coefsCrossFit( image -> lx, abuf[i], dabuf[i],
02019 bcoefs[i], n_b_fitcoefs, sigmaFactor )) )
02020 {
02021 cpl_msg_error ("wavelengthCalibration:"," cannot carry out the fitting of coefficients across the columns, for the coefficient with index: %d\n", i) ;
02022 for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02023 {
02024 cpl_free (abuf[i]) ;
02025 cpl_free (dabuf[i]) ;
02026 }
02027 cpl_free ( acoefs ) ;
02028 cpl_free ( dacoefs ) ;
02029 cpl_free ( abuf ) ;
02030 cpl_free ( dabuf ) ;
02031 return -1 ;
02032 }
02033 }
02034
02035 /* free all allocated memory */
02036 for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02037 {
02038 cpl_free (abuf[i]) ;
02039 cpl_free (dabuf[i]) ;
02040 }
02041 cpl_free ( acoefs ) ;
02042 cpl_free ( dacoefs ) ;
02043 cpl_free ( abuf ) ;
02044 cpl_free ( dabuf ) ;
02045
02046 return 0 ;
02047 }
02048
02049
02050 /*----------------------------------------------------------------------------
02051 Function : convolveImageByGauss()
02052 In : lineImage: emission line image
02053 hw: kernel half width of the gaussian
02054 response function
02055 Out : emission line image convolved with a gaussian
02056 Job : convolves an emission line image with a gaussian
02057 with user given integer half width by using the eclipse
02058 routine function1d_filter_lowpass().
02059 ---------------------------------------------------------------------------*/
02060
02061 OneImage * convolveImageByGauss( OneImage * lineImage,
02062 int hw )
02063 {
02064 OneImage * returnImage ;
02065 float column_buffer[lineImage->ly] ;
02066 float * filter ;
02067 int col, row ;
02068
02069 if ( lineImage == NullImage )
02070 {
02071 cpl_msg_error("convolveImageByGauss:"," no input image given!\n") ;
02072 return NullImage ;
02073 }
02074
02075 if ( hw < 1 )
02076 {
02077 cpl_msg_error("convolveImageByGauss:"," wrong half width given!\n") ;
02078 return NullImage ;
02079 }
02080
02081 /* allocate memory for returned image */
02082 if ( NULL == ( returnImage = new_image( lineImage -> lx, lineImage -> ly ) ))
02083 {
02084 cpl_msg_error("convolveImageByGauss:"," cannot allocate a new image\n");
02085 return NullImage ;
02086 }
02087
02088 /* go through the image columns and save them in a buffer */
02089 for ( col = 0 ; col < lineImage -> lx ; col++ )
02090 {
02091 for ( row = 0 ; row < lineImage -> ly ; row++ )
02092 {
02093 column_buffer[row] = lineImage -> data[col + row*lineImage->lx] ;
02094 }
02095
02096 /*---------------------------------------------------------------------
02097 * now low pass filter the columns by the gaussian and fill the return
02098 * image.
02099 */
02100 filter = function1d_filter_lowpass( column_buffer,
02101 lineImage -> ly,
02102 LOW_PASS_GAUSSIAN,
02103 hw ) ;
02104 for ( row = 0 ; row < lineImage -> ly ; row++ )
02105 {
02106 returnImage -> data[col + row*lineImage->lx] = filter[row] ;
02107 }
02108 function1d_del(filter) ;
02109 }
02110
02111 return returnImage ;
02112 }
02113
02114 /*----------------------------------------------------------------------------
02115 Function : definedResampling()
02116 In : image: source image to be calibrated
02117 calimage: wavelength map image
02118 n_params: number of fit parameters for the polynomial interpolation
02119 standard should be 3
02120 that means order of polynom + 1
02121 n_rows: desired number of rows for the final image,
02122 this will be the final number of spectral pixels
02123 in the final data cube.
02124 dispersion: dummy for the resulting dispersion
02125 minval: dummy for minimal wavelength value,
02126 maxval: dummy for maximal wavelength value
02127 centralLambda: dummy for the final central wavelength
02128 Out : wavelength calibrated source image,
02129 dispersion: resulting spectral dispersion (microns/pixel)
02130 is chosen as the minimum dispersion found in
02131 the wavelength map - 2% of this value
02132 minval: minimal wavelength value,
02133 maxval: maximal wavelength value
02134 centralLambda: final central wavelength value
02135 centralpix: row of central wavelength
02136 (in image coordinates!)
02137 Job : Given a source image and a corresponding wavelength
02138 calibration file this routine produces an image
02139 in which elements in a given row are associated
02140 with a single wavelength. It thus corrects for
02141 the wavelength shifts between adjacent elements
02142 in the rows of the input image. The output image
02143 is larger in the wavelength domain than the input
02144 image with pixels in each column corresponding to
02145 undefined (blank, ZERO) values. The distribution
02146 of these undefined values varies from column to
02147 column. The input image is resampled at discrete
02148 wavelength intervals using a polynomial interpolation
02149 routine.
02150 The wavelength intervals (dispersion) and the
02151 central wavelength are defined and stable for each
02152 used grating. Thus, each row has a defined wavelength
02153 for each grating. Only the number of rows can be
02154 changed by the user.
02155 ---------------------------------------------------------------------------*/
02156
02157 OneImage * definedResampling( OneImage * image,
02158 OneImage * calimage,
02159 int n_params,
02160 int* n_rows,
02161 double * dispersion,
02162 float * minval,
02163 float * maxval,
02164 double * centralLambda,
02165 int * centralpix )
02166 {
02167 OneImage * retImage ;
02168 OneImage * tempCalImage ;
02169 OneImage * tempImage ;
02170 float lambda ;
02171 float dif, lambda_renorm ;
02172 float imagecol[image->ly] ;
02173 float retimagecol[2560] ; /* retimagecol[n_rows] ; */
02174 float calcol[calimage->ly] ;
02175 float x_renorm[n_params] ;
02176 float * imageptr ;
02177 float sum, new_sum ;
02178 float disp, mindisp ;
02179 int calcolpos[2560];
02180 int i/*, j*/, col, row, testrow ;
02181 int half_width, firstpos ;
02182 int dispInd ;
02183 int n ;
02184 int flag;
02185 float temprow;
02186 float minLambda = 0. ;
02187 /*dpoint list[n_params] ;*/
02188 /*double * polycoeffs ;*/
02189 double poly ;
02190 /*float error;*/
02191 int zeroind ;
02192
02193
02194 if ( NULL == image )
02195 {
02196 cpl_msg_error("definedResampling:"," source image not given\n") ;
02197 return NullImage ;
02198 }
02199
02200 if ( NULL == calimage )
02201 {
02202 cpl_msg_error("definedResampling:"," wavelength map image not given\n") ;
02203 return NullImage ;
02204 }
02205
02206 if ( image -> lx != calimage -> lx ||
02207 image -> ly != calimage -> ly )
02208 {
02209 cpl_msg_error("definedResampling:"," source image and wavelength map image are not compatible in size\n") ;
02210 return NullImage ;
02211 }
02212
02213 if ( n_params < 1 )
02214 {
02215 cpl_msg_error ("definedResampling:"," wrong number of fit parameters given\n") ;
02216 return NullImage ;
02217 }
02218
02219 if ( n_params > 4 )
02220 {
02221 cpl_msg_warning("definedResampling:"," attention: very high number of fit parameters given, not tested !!!\n") ;
02222 }
02223
02224 /*if ( n_rows <= calimage -> ly)
02225 {
02226 cpl_msg_error ("definedResampling:"," number of rows of resampled image will be smaller than in wavelength calibration map, information would get lost!\n") ;
02227 return NullImage ;
02228 }*/
02229
02230 dispInd = 0 ;
02231 /* first determine the dispersion direction */
02232 for ( col = 0 ; col < calimage->lx ; col++ )
02233 {
02234 if ( qfits_isnan(calimage -> data[col]) || calimage -> data[col] <= 0. )
02235 {
02236 continue ;
02237 }
02238 if ((calimage->data[col] - calimage->data[col+(calimage->lx)*(calimage->ly-1)]) > 0. )
02239 {
02240 dispInd-- ;
02241 }
02242 else if ((calimage->data[col] - calimage->data[col+(calimage->lx)*(calimage->ly-1)]) < 0. )
02243 {
02244 dispInd++ ;
02245 }
02246 else
02247 {
02248 continue ;
02249 }
02250 }
02251 if ( dispInd == 0 )
02252 {
02253 cpl_msg_error("definedResampling:"," zero dispersion?\n");
02254 return NullImage ;
02255 }
02256
02257 /* mirror the wavelength map and the raw image if the dispersion is negative */
02258 if ( dispInd < 0 )
02259 {
02260 /* allocate a temp image */
02261 if ( NULL == ( tempCalImage = new_image( calimage->lx, calimage->ly ) ))
02262 {
02263 cpl_msg_error("definedResampling:"," cannot allocate a new image\n");
02264 return NullImage ;
02265 }
02266 if ( NULL == ( tempImage = new_image( image->lx, image->ly ) ))
02267 {
02268 cpl_msg_error("definedResampling:"," cannot allocate a new image\n");
02269 destroy_image(tempCalImage) ;
02270 return NullImage ;
02271 }
02272 for ( col = 0 ; col < calimage->lx ; col++ )
02273 {
02274 n = calimage->ly - 1 ;
02275 for ( row = 0 ; row < calimage->ly ; row++ )
02276 {
02277 tempCalImage->data[col+row*calimage->lx] = calimage->data[col+n*calimage->lx] ;
02278 tempImage->data[col+row*calimage->lx] = image->data[col+n*calimage->lx] ;
02279 n-- ;
02280 }
02281 }
02282 for ( i = 0 ; i < (int) image->nbpix ; i++ )
02283 {
02284 image->data[i] = tempImage->data[i] ;
02285 calimage->data[i] = tempCalImage->data[i] ;
02286 }
02287 destroy_image(tempCalImage) ;
02288 destroy_image(tempImage) ;
02289 }
02290
02291 /* determine the max and min pixel value in the first and the last row */
02292 *maxval = -FLT_MAX ;
02293 *minval = FLT_MAX ;
02294 mindisp = FLT_MAX ;
02295 for ( col = 0 ; col < calimage->lx ; col++ )
02296 {
02297 if ( qfits_isnan(calimage -> data[col]) || calimage -> data[col] <= 0. )
02298 {
02299 continue ;
02300 }
02301 disp = (calimage->data[col+(calimage->lx)*((calimage->ly)-1)]
02302 - calimage->data[col]) / (float)calimage->ly ;
02303 if ( mindisp > disp )
02304 {
02305 mindisp = disp ;
02306 }
02307 if ( *minval >= calimage -> data[col] )
02308 {
02309 *minval = calimage -> data[col] ;
02310 }
02311 if ( *maxval <= calimage -> data[col + (calimage->lx)*((calimage->ly)-1)] )
02312 {
02313 *maxval = calimage -> data[col + (calimage->lx)*((calimage->ly)-1)] ;
02314 }
02315 }
02316 /* find the used grating and set the dispersion to the defined value */
02317 if (*minval > 1.9 )
02318 {
02319 if ( calimage -> ly > 1024 && calimage -> ly < 3000)
02320 {
02321 *dispersion = DISPERSION_K_DITH ;
02322 *centralLambda = CENTRALLAMBDA_K ;
02323 }
02324 else if ( calimage -> ly < 2000)
02325 {
02326 *dispersion = DISPERSION_K ;
02327 *centralLambda = CENTRALLAMBDA_K ;
02328 }
02329 else
02330 {
02331 *dispersion = DISPERSION_K_DITH/2 ;
02332 *centralLambda = CENTRALLAMBDA_K ;
02333 }
02334 }
02335 else if (*minval < 1.2 )
02336 {
02337 if ( calimage -> ly > 1024 )
02338 {
02339 *dispersion = DISPERSION_J_DITH ;
02340 *centralLambda = CENTRALLAMBDA_J ;
02341 }
02342 else
02343 {
02344 *dispersion = DISPERSION_J ;
02345 *centralLambda = CENTRALLAMBDA_J ;
02346 }
02347 }
02348 else
02349 {
02350 if ( *maxval > 2.3 )
02351 {
02352 if ( calimage -> ly > 1024 )
02353 {
02354 *dispersion = DISPERSION_HK_DITH ;
02355 *centralLambda = CENTRALLAMBDA_HK ;
02356 }
02357 else
02358 {
02359 *dispersion = DISPERSION_HK ;
02360 *centralLambda = CENTRALLAMBDA_HK ;
02361 }
02362 }
02363 else
02364 {
02365 if ( calimage -> ly > 1024 )
02366 {
02367 *dispersion = DISPERSION_H_DITH ;
02368 *centralLambda = CENTRALLAMBDA_H ;
02369 }
02370 else
02371 {
02372 *dispersion = DISPERSION_H ;
02373 *centralLambda = CENTRALLAMBDA_H ;
02374 }
02375 }
02376 }
02377 /*if ( *minval + (float)n_rows * *dispersion < *maxval )
02378 {
02379 cpl_msg_error("definedResampling:"," given number of rows too small!\n");
02380 return NullImage ;
02381 }*/
02382 if ( (*maxval - *minval) / *dispersion < (float)calimage->ly )
02383 {
02384 cpl_msg_error("definedResampling:"," must be something wrong with the wavelength map!\n");
02385 return NullImage ;
02386 }
02387
02388 /* determine the central pixel and the lambda in the first image row */
02389 *n_rows = floor(floor(0.5+(*maxval - *minval) / *dispersion)/2+0.5)*2;
02390 *centralpix = *n_rows / 2 ;
02391 minLambda = *centralLambda - *dispersion * (float)*centralpix ;
02392 /*if ( (minLambda + *dispersion * n_rows) < *maxval )
02393 {
02394 cpl_msg_error("definedResampling:"," not enough rows defined \n");
02395 return NullImage ;
02396 }
02397 if ( minLambda > *minval )
02398 {
02399 cpl_msg_error("definedResampling:"," not enough rows defined \n");
02400 return NullImage ;
02401 }*/
02402
02403 /* allocate memory */
02404 if ( NULL == ( retImage = new_image( image -> lx, *n_rows ) ))
02405 {
02406 cpl_msg_error("definedResampling:"," cannot allocate a new image\n");
02407 return NullImage ;
02408 }
02409
02410 /* now go through the columns */
02411 for ( col = 0 ; col < retImage -> lx ; col++ )
02412 {
02413
02414 /*------------------------------------------------------------------
02415 * copy the columns of the source image and the wavemap image into
02416 * buffer arrays to speed things up
02417 */
02418 sum = 0. ;
02419 for ( row = 0 ; row < image -> ly ; row++ )
02420 {
02421 imagecol[row] = image -> data[col + row*image->lx] ;
02422 if (!qfits_isnan(imagecol[row]))
02423 {
02424 sum += imagecol[row] ;
02425 }
02426 calcol[row] = calimage -> data[col + row*calimage->lx] ;
02427 }
02428 for ( row = 0 ; row < retImage -> ly ; row++ )
02429 {
02430 retimagecol[row] = 0. ;
02431 calcolpos[row] = -1;
02432 }
02433
02434 for ( row=0 ; row < calimage->ly ; row++)
02435 {
02436 temprow = (calcol[row]- minLambda)/ *dispersion;
02437 if (temprow >= 0 && temprow < retImage->ly)
02438 calcolpos[(int) temprow] = row;
02439 }
02440
02441 zeroind = 0 ;
02442
02443
02444
02445 for ( row = 0 ; row < retImage -> ly ; row++ )
02446 {
02447 lambda = minLambda + *dispersion * (float) row ;
02448
02449 /*---------------------------------------------------------------
02450 * lambda must lie between the two available wavelength extremes
02451 * otherwise the image pixels are set to ZERO
02452 */
02453 if ( row < calimage->ly )
02454 {
02455 if ( qfits_isnan(calcol[row]) )
02456 {
02457 zeroind = 1 ;
02458 }
02459 }
02460
02461 if ( (lambda < calcol[0]) || (lambda > calcol[(calimage->ly)-1]) || zeroind == 1 )
02462 {
02463 retimagecol[row] = ZERO ;
02464 continue ;
02465 }
02466
02467 /*testrow = 0 ;
02468 while ( lambda > calcol[testrow] )
02469 {
02470 testrow++ ;
02471 }*/
02472 if (calcolpos[row]==-1) {
02473 if(row>= (*n_rows-1)) calcolpos[row] = calcolpos[row-1];
02474 if(row< (*n_rows-1)) calcolpos[row] = calcolpos[row+1];
02475 }
02476 if (lambda-calcol[calcolpos[row]-1]==0.) calcolpos[row]=calcolpos[row]-1;
02477 testrow = calcolpos[row];
02478
02479 /*--------------------------------------------------------------------
02480 * at this point calcol[testrow-1] < lambda <= calcol[testrow]
02481 * now determine the box position in which the polint fit is carried through.
02482 * the half width of the box is half the number of fit parameters.
02483 * Now we determine the start position of the fitting box and treat
02484 * the special case of being near the edge.
02485 */
02486
02487 if ( n_params % 2 == 0 )
02488 {
02489 half_width = (int)(n_params/2) - 1 ;
02490 }
02491 else
02492 {
02493 half_width = (int)(n_params/2) ;
02494 }
02495
02496
02497 if ( qfits_isnan(imagecol[testrow]) )
02498 {
02499 for ( i = row-half_width ; i < row-half_width+n_params ; i++ )
02500 {
02501 if (i < 0) continue ;
02502 if ( i >= retImage->ly ) continue ;
02503 retimagecol[i] = ZERO ;
02504 }
02505 imagecol[testrow] = 0. ;
02506 }
02507
02508 }
02509
02510
02511
02512 /* now loop over the rows and establish the lambda for each row */
02513 new_sum = 0. ;
02514 for ( row = 0 ; row < retImage -> ly ; row++ )
02515 {
02516 if ( qfits_isnan(retimagecol[row]) )
02517 {
02518 continue ;
02519 }
02520 lambda = minLambda + *dispersion * (float) row ;
02521
02522 /*---------------------------------------------------------------
02523 * lambda must lie between the two available wavelength extremes
02524 * otherwise the image pixels are set to ZERO
02525 */
02526 if ( (lambda < calcol[0]) || (lambda > calcol[(calimage->ly)-1]) )
02527 {
02528 retimagecol[row] = ZERO ;
02529 continue ;
02530 }
02531 /*testrow = 0 ;
02532 while ( lambda > calcol[testrow] )
02533 {
02534 testrow++ ;
02535 }*/
02536 if (calcolpos[row]==-1) {
02537 if(row >= (*n_rows-1)) calcolpos[row] = calcolpos[row-1];
02538 if(row < (*n_rows-1)) calcolpos[row] = calcolpos[row+1];
02539 }
02540 testrow = calcolpos[row];
02541
02542 /*--------------------------------------------------------------------
02543 * at this point calcol[testrow-1] < lambda <= calcol[testrow]
02544 * now determine the box position in which the polynomial interpolation is carried through.
02545 * the half width of the box is half the number of fit parameters.
02546 * Now we determine the start position of the fitting box and treat
02547 * the special case of being near the edge.
02548 */
02549
02550 if ( n_params % 2 == 0 )
02551 {
02552 half_width = (int)(n_params/2) - 1 ;
02553 }
02554 else
02555 {
02556 half_width = (int)(n_params/2) ;
02557 }
02558
02559 firstpos = testrow - half_width ;
02560 if ( firstpos < 0 )
02561 {
02562 firstpos = 0 ;
02563 }
02564 else if ( firstpos > ((calimage->ly)-n_params) )
02565 {
02566 firstpos = calimage->ly - n_params ;
02567 }
02568
02569 if ( qfits_isnan(imagecol[firstpos]) )
02570 {
02571 retimagecol[row] = ZERO ;
02572 continue ;
02573 }
02574
02575
02576 /* we must rescale the x-values (smaller than 1) for the fitting routine */
02577 dif = calcol[firstpos+n_params-1] - calcol[firstpos] ;
02578 for ( i = 0 ; i < n_params ; i++ )
02579 {
02580 x_renorm[i] = (calcol[firstpos + i] - calcol[firstpos]) / dif ;
02581 }
02582
02583
02584
02585 lambda_renorm = ( lambda - calcol[firstpos] ) / dif ;
02586
02587 imageptr = &imagecol[firstpos] ;
02588
02589 flag = 0;
02590 poly=nev_ille(x_renorm, imageptr, n_params-1, lambda_renorm, &flag);
02591
02592 new_sum += poly ;
02593 retimagecol[row] = poly ;
02594 }
02595
02596
02597 /* now renorm the total flux */
02598 for ( row = 0 ; row < retImage -> ly ; row++ )
02599 {
02600 if ( new_sum == 0. ) new_sum = 1. ;
02601 if ( qfits_isnan(retimagecol[row]) )
02602 {
02603 retImage->data[col+row*retImage->lx] = ZERO ;
02604 }
02605 else
02606 {
02607 /* rescaling is commented out because it delivers wrong results
02608 in case of appearance of blanks or bad pixels */
02609 retImage->data[col+row*retImage->lx] = retimagecol[row] /* * sum/new_sum*/ ;
02610 }
02611 }
02612
02613 }
02614
02615 return retImage ;
02616 }
02617
02618
02619 /*___oOo___*/
1.2.13.1 written by Dimitri van Heesch,
© 1997-2001