coltilt.c

00001 /*******************************************************************************
00002 * E.S.O. - VLT project
00003 *
00004 *
00005 *
00006 * who       when      what
00007 * --------  --------  ----------------------------------------------
00008 * schreib  19/12/00  created
00009 */
00010 
00011 /************************************************************************
00012 *   NAME
00013 *        coltilt.c -
00014 *        procedures to correct for tilted spectra 
00015 *
00016 *   SYNOPSIS
00017 *   1) float slopeOfSpectrum( OneImage * ns_image,
00018 *                             int        box_length,
00019 *                             float      fwhm,
00020 *                             float      minDiff )
00021 *   2) OneImage * shiftRows( OneImage * image,
00022 *                            float      slope )
00023 *   3) void parameterToAscii ( float * parameter, 
00024 *                              int n,
00025 *                              char * filename )
00026 *   4) float * asciiToParameter ( char * filename,
00027 *                                 int * n )
00028 *   5) double * curvatureOfSpectrum( OneImage * ns_image,
00029 *                            int        order,
00030 *                                    int        box_length,
00031 *                                    int        left_pos,
00032 *                                    int        right_pos,
00033 *                                    float      fwhm,
00034 *                                    float      minDiff )
00035 *
00036 *
00037 *   DESCRIPTION
00038 *   1)  determines the sub-pixel shifts of each row by using
00039 *       an image with at least one continuum spectrum of a pinhole
00040 *       this is done by searching the spectrum within the image
00041 *       then fitting the spectrum along the rows within a given box
00042 *       by a gaussian, so that the exact position is determined for
00043 *       each row. Afterwards, a straight line is fitted through the
00044 *       fitted positions. The slope of this linear fit is returned.
00045 *   2)  shifts the rows of a raw image by using the output of 
00046 *       slopeOfSpectrum and applying polynomial interpolation
00047 *   3)  stores parameters in an ASCII file 
00048 *   4)  writes parameters stored in an ASCII file in an float array
00049 *   5)  this routine determines the curvature of a spectrum by fitting
00050 *       a polynomial to a continuum spectrum. This is done by using
00051 *       an image with at least one continuum spectrum of a pinhole.
00052 *       this is done by searching the spectrum within the image
00053 *       then fitting the spectrum along the rows within a given box
00054 *       by a gaussian, so that the exact position is determined for
00055 *       each row. Afterwards, a polynomial is fitted through the
00056 *       found positions. The polynomial coefficients are returned.
00057 *
00058 *
00059 *   FILES
00060 *
00061 *   ENVIRONMENT
00062 *
00063 *   RETURN VALUES
00064 *
00065 *   CAUTIONS
00066 *
00067 *   EXAMPLES
00068 *
00069 *   SEE ALSO
00070 *
00071 *   BUGS
00072 *
00073 *------------------------------------------------------------------------
00074 */
00075 
00076 #define POSIX_SOURCE 1
00077 #include "vltPort.h"
00078 
00079 /*
00080  * System Headers
00081  */
00082 
00083 /*
00084  * Local Headers
00085  */
00086 
00087 #include "coltilt.h"
00088 
00089 /*----------------------------------------------------------------------------
00090  *                          Function codes
00091  *--------------------------------------------------------------------------*/
00092 
00093 /*----------------------------------------------------------------------------
00094    Function     :       slopeOfSpectrum()
00095    In           :       ns_image:   image with at least one continuum spectrum
00096                                     of a pinhole
00097                         box_length: width of the box in which the lines
00098                                     are fit by a gaussian
00099                         fwhm:       first guess of the full width at half maximum
00100                         minDiff:    minimum amplitude threshold below which 
00101                     the fit is not carried through
00102    Out          :       slope of a straight line fitted to the spectrum.
00103                         ZERO if something went wrong.
00104    Job          :       determines the sub-pixel shifts of each row by using
00105                         an image with at least one continuum spectrum of a pinhole
00106                         this is done by searching the spectrum within the image
00107                         then fitting the spectrum along the rows within a given box
00108                         by a gaussian, so that the exact position is determined for
00109                         each row. Afterwards, a straight line is fitted through the
00110                         fitted positions. The slope of this linear fit is returned.
00111  ---------------------------------------------------------------------------*/
00112 
00113 float slopeOfSpectrum( OneImage * ns_image,
00114                        int        box_length,
00115                        float      fwhm,
00116                        float      minDiff )
00117 {
00118     int i, k, row, col ;
00119     int counter, iters ;
00120     int xdim, ndat, its, numpar ;
00121     float maxval ;
00122     float tol, lab ;
00123     float col_value[ns_image->ly] ;
00124     float column_value[ns_image->ly] ;
00125     pixelvalue col_position[ns_image->ly] ;
00126     int column_position[ns_image->ly] ;
00127     int col_median ;
00128     float x_position[ns_image->ly] ;
00129     float * xdat, * wdat ;
00130     int * mpar ;
00131     Vector * line ;
00132     FitParams ** dec_par ;
00133     FitParams *  par ;
00134     float x[ns_image->lx], y[ns_image->lx] ; 
00135     float sig[ns_image->lx] ;
00136     int position ;
00137     int ndata, mwt ;
00138     float a, b, siga, sigb, chi2, q ;
00139     int bad_ind ;
00140 
00141 
00142     if ( ns_image == NullImage )
00143     {
00144         cpl_msg_error("slopeOfSpectrum:"," sorry, no image given\n") ;
00145         return FLAG ;
00146     }
00147     if ( box_length <= 1 || box_length >= (int) sqrt(ns_image -> lx) )
00148     {
00149         cpl_msg_error("slopeOfSpectrum:"," wrong box length given\n") ;
00150         return FLAG ;
00151     }
00152     if ( fwhm < 1. || fwhm > 10. )
00153     {
00154         cpl_msg_error("slopeOfSpectrum","wrong full width at half maximum given\n") ;
00155         return FLAG ;
00156     }
00157     if ( minDiff < 1. )
00158     {
00159         cpl_msg_error("slopeOfSpectrum","wrong amplitude threshold given\n") ;
00160         return FLAG ;
00161     }
00162 
00163     /* go through the image rows */
00164     for ( row = 0 ; row < ns_image -> ly ; row++ )
00165     {
00166         col_value[row] = -FLT_MAX ;
00167         col_position[row] = -1. ;
00168         /* find the maximum value in each row and store the found column */
00169         for ( col = 0 ; col < ns_image -> lx ; col++ )
00170         {
00171             if ( ns_image -> data[col+row*ns_image->lx] > col_value[row] )
00172             {
00173                 col_value[row]    = ns_image -> data[col+row*ns_image->lx] ;
00174                 col_position[row] = (pixelvalue)col ;
00175             }
00176         }
00177     }
00178 
00179     /* now determine the median of the found columns to be sure 
00180        to have the brightest spectrum */
00181     col_median = (int)median(col_position, ns_image->ly) ;
00182     cpl_msg_info ("slopeOfSpectrum","median column position of brightest spectrum %d\n", col_median) ;
00183 
00184     /* now find the peaks around col_median over the whole spectral range */
00185     for ( row = 0 ; row < ns_image -> ly ; row++ )
00186     {
00187         x_position[row] = 0. ;
00188         column_value[row] = -FLT_MAX ;
00189         column_position[row] = -1 ;
00190         for ( col = col_median - box_length ; col <= col_median + box_length ; col++ )
00191         {
00192             if ( ns_image -> data[col+row*ns_image->lx] > column_value[row] )
00193             {
00194                 column_value[row] = ns_image -> data[col+row*ns_image->lx] ; 
00195                 column_position[row] = col ;
00196             }
00197         }
00198 
00199         /* allocate memory for the array where the line is fitted in */
00200         if ( NULL == (line = newVector (2*box_length + 1)) )
00201         {
00202             cpl_msg_error ("slopeOfSpectrum","cannot allocate new Vector in row %d\n", row) ;
00203             return FLAG ;
00204         }
00205 
00206         /* allocate memory */
00207         xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00208         wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00209         mpar = (int *)   cpl_calloc( MAXPAR, sizeof (int) ) ;
00210         dec_par = newFitParams(1) ;
00211     par = dec_par[0];
00212 
00213         counter = 0 ;
00214         bad_ind = 0 ;
00215         /* store the values to fit in a Vector object */
00216         for ( col = column_position[row] - box_length ; col <= column_position[row] + box_length ; col++ )
00217         {
00218             if ( col < 0 || col >= ns_image -> lx )
00219             {
00220                 cpl_msg_error ("slopeOfSpectrum","wrong spectrum position or box_length given in row: %d\n", row) ;
00221                 cpl_free (xdat) ;
00222                 cpl_free (wdat) ;
00223                 cpl_free (mpar) ;
00224                 destroyFitParams(dec_par) ;
00225                 destroyVector( line ) ;
00226                 return FLAG ;
00227             }
00228             else if ( qfits_isnan(ns_image->data[col+row*ns_image->lx]) )
00229             {
00230                 bad_ind = 1 ;
00231             }
00232             else
00233             {
00234                 line -> data[counter] = ns_image -> data[col + row*ns_image->lx] ;
00235                 counter++ ;
00236             }
00237         }
00238  
00239         /* go to the next row if a bad pixel is inside the box */
00240         if ( bad_ind == 1 )
00241         {
00242             cpl_msg_warning ("slopeOfSpectrum","sorry, bad pixel inside fitting box in row: %d\n", row) ;
00243             cpl_free (xdat) ;
00244             cpl_free (wdat) ;
00245             cpl_free (mpar) ;
00246             destroyFitParams(dec_par) ;
00247             destroyVector( line ) ;
00248             continue ;
00249         }
00250 
00251         /*--------------------------------------------------------------------
00252          * go through the line vector
00253          * determine the maximum pixel value in the line vector
00254          */
00255         maxval = -FLT_MAX ;
00256         position = -INT32_MAX ;
00257         for ( i = 0 ; i < counter ; i++ )
00258         {
00259             xdat[i] = i ;
00260             wdat[i] = 1.0 ;
00261             if ( line -> data[i] >= maxval )
00262             {
00263                 maxval = line -> data[i] ;
00264                 position = i ;
00265             }
00266         }
00267 
00268         /* set initial values for the fitting routine */
00269         xdim     = XDIM ;
00270         ndat     = line -> n_elements ;
00271         numpar   = MAXPAR ;
00272         tol      = TOL ;
00273         lab      = LAB ;
00274         its      = ITS ;
00275         (*par).fit_par[1] = fwhm ;
00276         (*par).fit_par[2] = (float) position ;
00277         (*par).fit_par[3] = (float) (line -> data[0] + line -> data[line->n_elements - 1]) / 2.0 ;
00278 
00279         (*par).fit_par[0]  = maxval - ((*par).fit_par[3]) ;
00280         /* exclude negative peaks and low signal cases */
00281         if ( (*par).fit_par[0] < minDiff )
00282         {
00283             cpl_msg_warning ("slopeOfSpectrum","sorry, negative peak or signal of line too low to fit in row: %d\n", row) ;
00284             cpl_free (xdat) ;
00285             cpl_free (wdat) ;
00286             cpl_free (mpar) ;
00287             destroyFitParams(dec_par) ;
00288             destroyVector( line ) ;
00289             continue ;
00290         }
00291 
00292         for ( k = 0 ; k < MAXPAR ; k++ )
00293         {
00294             (*par).derv_par[k] = 0.0 ;
00295             mpar[k] = 1 ;
00296         }
00297 
00298         /* finally, do the least square fit using a gaussian */
00299         if ( 0 > ( iters = lsqfit_c( xdat, &xdim, line -> data, wdat, &ndat, (*par).fit_par,
00300                                      (*par).derv_par, mpar, &numpar, &tol, &its, &lab )) )
00301         {
00302             cpl_msg_warning ("slopeOfSpectrum","lsqfit_c: least squares fit failed in row: %d, error no.: %d\n", row, iters) ;
00303             cpl_free (xdat) ;
00304             cpl_free (wdat) ;
00305             cpl_free (mpar) ;
00306             destroyFitParams(dec_par) ;
00307             destroyVector( line ) ;
00308             continue ;
00309         }
00310 
00311         /* check for negative fit results */
00312         if ( (*par).fit_par[0] <= 0. || (*par).fit_par[1] <= 0. ||
00313              (*par).fit_par[2] <= 0. )
00314         {
00315             cpl_msg_warning ("slopeOfSpectrum","negative parameters as fit result, not used! in row: %d\n", row) ;
00316             cpl_free (xdat) ;
00317             cpl_free (wdat) ;
00318             cpl_free (mpar) ;
00319             destroyFitParams(dec_par) ;
00320             destroyVector( line ) ;
00321             continue ;
00322         }
00323 
00324         /* correct the fitted position for the given row of the line in image coordinates */
00325         x_position[row] =  (float) (column_position[row] - box_length) + (*par).fit_par[2] ;
00326 
00327         /* store the fit errors of the positions as weights for the later linear fit */
00328         sig[row] = (*par).derv_par[2] ;
00329 
00330         /* free memory */
00331         destroyFitParams(dec_par) ;
00332         destroyVector ( line ) ;
00333         cpl_free ( xdat ) ;
00334         cpl_free ( wdat ) ;
00335         cpl_free ( mpar ) ;
00336     }
00337 
00338     /* ------------------------------------------------------------------------------------
00339      * now that we have a sub-pixel resolved list of spectral maxima stored in x_position[row]
00340      * We can fit a flux weighted straight line to the positions to determine 
00341      * the spectral column shifts.
00342      */
00343     ndata = 0 ;
00344     for ( row = 0 ; row < ns_image -> lx ; row++ )
00345     {
00346         if ( x_position[row] == 0. || sig[row] == 0. )
00347         {
00348             continue ;
00349         }
00350         else 
00351         {
00352             y[ndata]   = x_position[row] ;
00353             x[ndata]   = (float)row ;
00354             sig[ndata] = sig[row] ;
00355             ndata++ ;
00356         }
00357     }
00358     if ( ndata < 10 )
00359     {
00360         cpl_msg_error("slopeOfSpectrum","not enough positions to do the linear fit\n") ;
00361         return FLAG ;
00362     }
00363 
00364     /* now do the fit and return the slope of the straight line */
00365     mwt = 0 ;
00366     myfit(x, y, ndata, sig, mwt, &a, &b, &siga, &sigb, &chi2, &q) ;
00367     return b ;
00368 }
00369 
00370 /*----------------------------------------------------------------------------
00371    Function     :       shiftRows()
00372    In           :       image:   raw image in which the rows should be shifted
00373                         slope:   slope of a fitted straight line along a spectrum
00374                                  output of slopeOfSpectrum
00375                         n_order: order of the interpolation polynom
00376    Out          :       resulting image with shifted rows
00377    Job          :       shifts the rows of a raw image by using the output of 
00378                         slopeOfSpectrum and applying polynomial interpolation
00379  ---------------------------------------------------------------------------*/
00380 
00381 OneImage * shiftRows( OneImage * image,
00382                       float      slope, 
00383                       int        n_order )
00384 {
00385     OneImage * returnImage ;
00386     float xshift ;
00387     int   intshift = 0 ;
00388     float sum, new_sum ;
00389     float xnum[n_order + 1] ;
00390     float row_data[image->lx] ;
00391     float corrected_row_data[image->lx] ;
00392     float eval/*, dy*/ ;
00393     float * imageptr ;
00394     int col, row ;
00395     int firstpos, n_points ;
00396     int i ;
00397     int flag;
00398 
00399     if ( image == NullImage )
00400     {
00401         cpl_msg_error("shiftRows","sorry, no image given\n") ;
00402         return NullImage ;
00403     }
00404 
00405     if ( slope == 0. )
00406     {
00407         cpl_msg_error("shiftRows","there is no need to shift the image rows!\n") ;
00408         return NullImage ;
00409     }
00410 
00411     if ( n_order <= 0 )
00412     {
00413         cpl_msg_error("shiftRows","wrong order of interpolation polynom given!\n") ;
00414         return NullImage ;
00415     }
00416 
00417     returnImage = copy_image( image ) ;
00418 
00419     n_points = n_order + 1 ;
00420     if ( n_points % 2 == 0 )
00421     {
00422         firstpos = (int)(n_points/2) - 1 ;
00423     }
00424     else
00425     {
00426         firstpos = (int)(n_points/2) ;
00427     }
00428 
00429     /* fill the xa[] array for the polint function */
00430     for ( i = 0 ; i < n_points ; i++ )
00431     {
00432         xnum[i] = i ;
00433     }
00434 
00435     /* go through the image rows */
00436     for ( row = 0 ; row < image -> ly ; row++ )
00437     {
00438         /* determine the shift for each row, the middle row is not shifted */
00439         xshift = slope * (float)( (image -> ly / 2) - row ) ;
00440 
00441         intshift = nint(xshift) ;
00442         xshift = xshift - (float)intshift ;
00443 
00444         for ( col = 0 ; col < image -> lx ; col++ )
00445         {
00446             corrected_row_data[col] = 0. ;
00447         }
00448         sum = 0. ; /* initialize flux for later rescaling */
00449         /* go through the image columns */
00450         for ( col = 0 ; col < image -> lx ; col++ )
00451         {
00452         /* consider integer pixel shifts */
00453         if ( intshift < 0 )
00454             {
00455         if ( col-intshift < image->lx )
00456         {
00457             row_data[col] = image->data[col-intshift+row*image->lx] ;
00458                 }
00459         else
00460         {
00461             row_data[col] = 0. ;
00462                 }
00463         }
00464         else if ( intshift > 0 )
00465         {
00466         if ( col - intshift >= 0 )
00467         {
00468             row_data[col] = image->data[col-intshift+row*image->lx] ;
00469                 }
00470         else
00471         {
00472             row_data[col] = 0. ;
00473         }
00474         }
00475         else
00476         {
00477                row_data[col] = image -> data[col+row*image->lx] ;
00478         }
00479 
00480             /* don't consider the edge pixels for flux calculation */
00481             if ( col != 0 && col != image->lx - 1 && !qfits_isnan(row_data[col]) )
00482             {
00483                 sum += row_data[col] ;
00484             }
00485             if (qfits_isnan(row_data[col])) 
00486             {
00487                 row_data[col] = 0. ;
00488                 for ( i = col - firstpos ; i < col - firstpos + n_points ; i++ )
00489                 {
00490                     if ( i < 0 ) continue ;
00491                     if ( i >= image->lx) continue ;
00492                     corrected_row_data[i] = ZERO ;
00493                 }
00494             }
00495         }
00496 
00497 
00498         /* --------------------------------------------------------------------
00499          * now we do the polynomial interpolation to achieve the fractional
00500          * shift that means call polint
00501          */
00502         new_sum = 0. ;
00503         for ( col = 0 ; col < image -> lx ; col++ )
00504         {
00505             /* ---------------------------------------------------------------
00506              * now determine the arrays of size n_points with which the
00507              * polynom is determined and determine the position eval
00508              * where the polynom is evaluated in polint of N.R..
00509              * Take care of the points near the row edges!
00510              */
00511             if ( qfits_isnan(corrected_row_data[col]) )
00512             {
00513                  continue ;
00514             }
00515             if ( col - firstpos < 0 )
00516             {
00517                 imageptr = &row_data[0] ;
00518                 eval     = (float)col - xshift ;
00519             }
00520             else if ( col - firstpos + n_points >= image -> lx )
00521             {
00522                 imageptr = &row_data[image->lx - n_points] ;
00523                 eval     = (float)(col + n_points - image->lx) - xshift ;
00524             }
00525             else
00526             {
00527                 imageptr = &row_data[col-firstpos] ;
00528                 eval     = (float)firstpos - xshift ;
00529             }
00530 
00531         flag=0;
00532         corrected_row_data[col]=nev_ille(xnum,imageptr,n_order,eval,&flag);
00533             /*polint( xnum - 1, imageptr, n_points, eval, &corrected_row_data[col], &dy ) ;*/
00534 
00535             /* don't take the edge points to calculate the scaling factor */
00536             if ( col != 0 && col != image->lx - 1 && !qfits_isnan(corrected_row_data[col]) )
00537             {
00538                 new_sum += corrected_row_data[col] ;
00539             }
00540         }
00541 
00542         if ( new_sum == 0. )
00543         {
00544             new_sum = 1. ;
00545         }
00546         for ( col = 0 ; col < image -> lx ; col++ )
00547         {
00548             if ( qfits_isnan(corrected_row_data[col]))
00549             {
00550                 returnImage->data[col+row*image->lx] = ZERO ;
00551             }
00552             else
00553             {
00554                 /* rescale the row data and fill the returned image */
00555               /* This gives sometimes inconsistent results if bad pixels are around */
00556               /* rescaling is commented out because it delivers wrong results
00557                  in case of appearance of blanks or bad pixels */
00558              /*   corrected_row_data[col] *= sum / new_sum ; */
00559                 returnImage->data[col+row*image->lx] = corrected_row_data[col] ;
00560             }
00561         }
00562     }      
00563     return returnImage ;
00564 }
00565 
00566 
00567 /*----------------------------------------------------------------------------
00568    Function     :       parameterToAscii()
00569    In           :       parameter: float parameter array to be stored in an ASCII file
00570                         number: number of parameters
00571                         filename: filename of ASCII file
00572    Out          :       nothing
00573    Job          :       stores parameters in an ASCII file 
00574  ---------------------------------------------------------------------------*/
00575 
00576 void parameterToAscii ( float * parameter, 
00577                         int n,
00578                         char * filename )
00579 {
00580     FILE * fp ;
00581     int     i ;
00582 
00583     if ( parameter == NULL || filename == NULL || n <= 0 )
00584     {
00585         cpl_msg_error ("parameterToAscii","input is missing or wrong!\n") ;
00586         return ;
00587     }
00588     
00589     if ( NULL == (fp = fopen ( filename, "w" ) ) )
00590     {
00591         cpl_msg_error("parameterToAscii","cannot open %s\n", filename) ;
00592         return ;
00593     }
00594 
00595     for ( i = 0 ; i < n ; i++ )
00596     {
00597         fprintf (fp, "%le\n", parameter[i] ) ;
00598     }
00599     fclose (fp ) ;
00600 }
00601 
00602 /*----------------------------------------------------------------------------
00603    Function     :       asciiToParameter()
00604    In           :       filename: filename of ASCII file
00605    Out          :       n: number of parameters in the parameter array
00606                         array of parameters
00607    Job          :       writes parameters stored in an ASCII file in an float array
00608  ---------------------------------------------------------------------------*/
00609 
00610 float * asciiToParameter ( char * filename,
00611                            int * n )
00612 {
00613     FILE * fp ;
00614     float * parameter ;
00615     int     i ;
00616 
00617     if ( filename == NULL || n == NULL )
00618     {
00619         cpl_msg_error ("asciiToParameter","input is missing or wrong!\n") ;
00620         return NULL ;
00621     }
00622     
00623     if ( NULL == (fp = fopen ( filename, "r" ) ) )
00624     {
00625         cpl_msg_error("asciiToParameter","cannot open %s\n", filename) ;
00626         return NULL ;
00627     }
00628 
00629     /* allocate memory */
00630 
00631     if ( NULL == ( parameter = (float*) cpl_calloc (ESTIMATE, sizeof(float)) ) )
00632     {
00633         cpl_msg_error ("asciiToParameter","cannot allocate memory\n") ;
00634         return NULL ;
00635     }
00636 
00637     i = 0 ;
00638     while ( fscanf(fp, "%g\n", &parameter[i]) != EOF )
00639     {
00640         i++ ;
00641     }
00642     *n = i ;
00643 
00644     fclose (fp ) ;
00645 
00646     return parameter ;
00647 }
00648 
00649 /*----------------------------------------------------------------------------
00650    Function     :       curvatureOfSpectrum()
00651    In           :       ns_image:   image with at least one continuum spectrum
00652                                     of a pinhole
00653             order:      order of the fit polynomial
00654                         box_length: width of the box in which the lines
00655                                     are fit by a gaussian
00656                         left_pos, right_pos: left and right positions between which
00657                          the spectrum should be located
00658                         fwhm:       first guess of the full width at half maximum
00659                         minDiff: minium amplitude threshold below which 
00660                  the fit is not carried through
00661    Out          :       resulting polynomial coefficients.
00662    Job          :       this routine determines the curvature of a spectrum by fitting
00663             a polynomial to a continuum spectrum. This is done by using
00664                         an image with at least one continuum spectrum of a pinhole.
00665                         this is done by searching the spectrum within the image
00666                         then fitting the spectrum along the rows within a given box
00667                         by a gaussian, so that the exact position is determined for
00668                         each row. Afterwards, a polynomial is fitted through the
00669                         found positions. The polynomial coefficients are returned.
00670  ---------------------------------------------------------------------------*/
00671 
00672 double * curvatureOfSpectrum( OneImage * ns_image,
00673                  int        order,
00674                          int        box_length,
00675                  int        left_pos,
00676              int        right_pos,
00677                          float      fwhm,
00678                          float      minDiff )
00679 {
00680     int i, k, row, col ;
00681     int counter, iters ;
00682     int xdim, ndat, its, numpar ;
00683     float maxval ;
00684     float tol, lab ;
00685     float col_value[ns_image->ly] ;
00686     float column_value[ns_image->ly] ;
00687     pixelvalue col_position[ns_image->ly] ;
00688     int column_position[ns_image->ly] ;
00689     int col_median ;
00690     float x_position[ns_image->ly] ;
00691     float * xdat, * wdat ;
00692     int * mpar ;
00693     Vector * line ;
00694     FitParams ** dec_par ;
00695     FitParams * par ;
00696     int position ;
00697     int ndata ;
00698     int bad_ind ;
00699     dpoint * list ;
00700     double * coeffs ;
00701     double offset ;
00702 
00703 
00704     if ( ns_image == NullImage )
00705     {
00706         cpl_msg_error("curvatureOfSpectrum","sorry, no image given\n") ;
00707         return NULL ;
00708     }
00709     if ( box_length <= 1 || box_length >= right_pos - left_pos )
00710     {
00711         cpl_msg_error("curvatureOfSpectrum","wrong box length given\n") ;
00712         return NULL ;
00713     }
00714     if ( fwhm < 1. || fwhm > 10. )
00715     {
00716         cpl_msg_error("curvatureOfSpectrum","wrong full width at half maximum given\n") ;
00717         return NULL ;
00718     }
00719     if ( left_pos < 0 || right_pos <= left_pos || right_pos > ns_image ->lx )
00720     {
00721         cpl_msg_error("curvatureOfSpectrum","wrong left and right positions\n") ;
00722         return NULL ;
00723     }
00724     if ( minDiff < 1. )
00725     {
00726         cpl_msg_error("curvatureOfSpectrum","wrong amplitude threshold given! \n") ;
00727         return NULL ;
00728     }
00729 
00730     /* go through the image rows */
00731     for ( row = 0 ; row < ns_image -> ly ; row++ )
00732     {
00733         col_value[row] = -FLT_MAX ;
00734         col_position[row] = -1. ;
00735         /* find the maximum value in each row and store the found column */
00736         for ( col = left_pos ; col < right_pos ; col++ )
00737         {
00738             if ( ns_image -> data[col+row*ns_image->lx] > col_value[row] )
00739             {
00740                 col_value[row]    = ns_image -> data[col+row*ns_image->lx] ;
00741                 col_position[row] = (pixelvalue)col ;
00742             }
00743         }
00744     }
00745 
00746     /* now determine the median of the found columns to be sure 
00747        to have the brightest spectrum */
00748     col_median = (int)median(col_position, right_pos - left_pos) ;
00749 
00750     /* now find the peaks around col_median over the whole spectral range */
00751     for ( row = 0 ; row < ns_image -> ly ; row++ )
00752     {
00753         x_position[row] = 0. ;
00754         column_value[row] = -FLT_MAX ;
00755         column_position[row] = -1 ;
00756         for ( col = col_median - box_length ; col <= col_median + box_length ; col++ )
00757         {
00758             if ( ns_image -> data[col+row*ns_image->lx] > column_value[row] )
00759             {
00760                 column_value[row] = ns_image -> data[col+row*ns_image->lx] ; 
00761                 column_position[row] = col ;
00762             }
00763         }
00764 
00765         /* allocate memory for the array where the line is fitted in */
00766         if ( NULL == (line = newVector (2*box_length + 1)) )
00767         {
00768             cpl_msg_error ("curvatureOfSpectrum","cannot allocate new Vector in row: %d\n", row) ;
00769             return NULL ;
00770         }
00771 
00772         /* allocate memory */
00773         xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00774         wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
00775         mpar = (int *)   cpl_calloc( MAXPAR, sizeof (int) ) ;
00776         dec_par = newFitParams(1) ;
00777         par = dec_par[0];
00778     
00779         counter = 0 ;
00780         bad_ind = 0 ;
00781         /* store the values to fit in a Vector object */
00782         for ( col = column_position[row] - box_length ; col <= column_position[row] + box_length ; col++ )
00783         {
00784             if ( col < 0 || col >= ns_image -> lx )
00785             {
00786                 cpl_msg_error ("curvatureOfSpectrum","wrong spectrum position or box_length given in row: %d\n", row) ;
00787                 cpl_free (xdat) ;
00788                 cpl_free (wdat) ;
00789                 cpl_free (mpar) ;
00790                 destroyFitParams(dec_par) ;
00791                 destroyVector( line ) ;
00792                 return NULL ;
00793             }
00794             else if ( qfits_isnan(ns_image->data[col+row*ns_image->lx]) )
00795             {
00796                 bad_ind = 1 ;
00797             }
00798             else
00799             {
00800                 line -> data[counter] = ns_image -> data[col + row*ns_image->lx] ;
00801                 counter++ ;
00802             }
00803         }
00804  
00805         /* go to the next row if a bad pixel is inside the box */
00806         if ( bad_ind == 1 )
00807         {
00808             cpl_msg_warning ("curvatureOfSpectrum","sorry, bad pixel inside fitting box in row: %d\n", row) ;
00809             cpl_free (xdat) ;
00810             cpl_free (wdat) ;
00811             cpl_free (mpar) ;
00812             destroyFitParams(dec_par) ;
00813             destroyVector( line ) ;
00814             continue ;
00815         }
00816 
00817         /*--------------------------------------------------------------------
00818          * go through the line vector
00819          * determine the maximum pixel value in the line vector
00820          */
00821         maxval = -FLT_MAX ;
00822         position = -INT32_MAX ;
00823         for ( i = 0 ; i < counter ; i++ )
00824         {
00825             xdat[i] = i ;
00826             wdat[i] = 1.0 ;
00827             if ( line -> data[i] >= maxval )
00828             {
00829                 maxval = line -> data[i] ;
00830                 position = i ;
00831             }
00832         }
00833 
00834         /* set initial values for the fitting routine */
00835         xdim     = XDIM ;
00836         ndat     = line -> n_elements ;
00837         numpar   = MAXPAR ;
00838         tol      = TOL ;
00839         lab      = LAB ;
00840         its      = ITS ;
00841         (*par).fit_par[1] = fwhm ;
00842         (*par).fit_par[2] = (float) position ;
00843         (*par).fit_par[3] = (float) (line -> data[0] + line -> data[line->n_elements - 1]) / 2.0 ;
00844 
00845         (*par).fit_par[0]  = maxval - ((*par).fit_par[3]) ;
00846         /* exclude negative peaks and low signal cases */
00847         if ( (*par).fit_par[0] < minDiff )
00848         {
00849             cpl_msg_warning ("curvatureOfSpectrum","sorry, negative peak or signal of line too low to fit in row: %d\n", row) ;
00850             cpl_free (xdat) ;
00851             cpl_free (wdat) ;
00852             cpl_free (mpar) ;
00853             destroyFitParams(dec_par) ;
00854             destroyVector( line ) ;
00855             continue ;
00856         }
00857 
00858         for ( k = 0 ; k < MAXPAR ; k++ )
00859         {
00860             (*par).derv_par[k] = 0.0 ;
00861             mpar[k] = 1 ;
00862         }
00863 
00864         /* finally, do the least square fit using a gaussian */
00865         if ( 0 > ( iters = lsqfit_c( xdat, &xdim, line -> data, wdat, &ndat, (*par).fit_par,
00866                                      (*par).derv_par, mpar, &numpar, &tol, &its, &lab )) )
00867         {
00868             cpl_msg_warning ("curvatureOfSpectrum","lsqfit_c: least squares fit failed in row: %d, error no.: %d\n", row, iters) ;
00869             cpl_free (xdat) ;
00870             cpl_free (wdat) ;
00871             cpl_free (mpar) ;
00872             destroyFitParams(dec_par) ;
00873             destroyVector( line ) ;
00874             continue ;
00875         }
00876 
00877         /* check for negative fit results */
00878         if ( (*par).fit_par[0] <= 0. || (*par).fit_par[1] <= 1. ||
00879              (*par).fit_par[2] <= 0. )
00880         {
00881             cpl_msg_warning ("curvatureOfSpectrum","negative parameters as fit result, not used! in row: %d\n", row) ;
00882             cpl_free (xdat) ;
00883             cpl_free (wdat) ;
00884             cpl_free (mpar) ;
00885             destroyFitParams(dec_par) ;
00886             destroyVector( line ) ;
00887             continue ;
00888         }
00889 
00890         /* correct the fitted position for the given row of the line in image coordinates */
00891         x_position[row] =  (float) (column_position[row] - box_length) + (*par).fit_par[2] ;
00892         printf("%d %f %f\n",row, (*par).fit_par[1], x_position[row] ) ;
00893 
00894         /* free memory */
00895         destroyFitParams(dec_par) ;
00896         destroyVector ( line ) ;
00897         cpl_free ( xdat ) ;
00898         cpl_free ( wdat ) ;
00899         cpl_free ( mpar ) ;
00900     }
00901     /* now allocate memory for the data to fit */
00902     if ( NULL == ( list = (dpoint*) cpl_calloc ( ns_image->ly, sizeof (dpoint)) ) ) 
00903     {
00904         cpl_msg_error("curvatureOfSpectrum","could not allocate memory!\n") ;
00905         return NULL ;
00906     }
00907 
00908     /* ------------------------------------------------------------------------------------
00909      * now that we have a sub-pixel resolved list of spectral maxima stored in x_position[row]
00910      * We can fit a flux weighted straight line to the positions to determine 
00911      * the spectral column shifts.
00912      */
00913     offset = (double) ns_image->ly/2. ;
00914     ndata = 0 ;
00915     for ( row = 0 ; row < ns_image -> ly ; row++ )
00916     {
00917         if ( x_position[row] == 0. )
00918         {
00919             continue ;
00920         }
00921         else 
00922         {
00923             list[ndata].y   = (double)x_position[row] ;
00924             list[ndata].x   = (double)row - offset ;
00925             ndata++ ;
00926         }
00927     }
00928     
00929 
00930     if ( NULL == (coeffs = fit_1d_poly(order, list, ndata, NULL)) )
00931     {
00932         cpl_msg_error("curvatureOfSpectrum","eclipse function fit_1d_poly() did not work!\n") ;
00933         return NULL ;
00934     }
00935     cpl_free ( list ) ;
00936      
00937     return coeffs ;
00938 }
00939 
00940 /*----------------------------------------------------------------------------
00941    Function     :       image_warp()
00942             unwarps SPIFFI spectra
00943  ---------------------------------------------------------------------------*/
00944  
00945 
00946 OneImage * image_warp( OneImage * image,
00947                       char      * kernel_type, 
00948                       char      * poly_table )
00949 {
00950         OneImage  * warped;
00951     char        ch,temp[150];
00952     char        poly_string[150];
00953     int     count=0;
00954     FILE    *   poly_in;
00955     /* Following are for polynomial transforms */
00956     poly2d  *   poly_u;         /* polynomial definition */
00957     poly2d  *   poly_v;         /* polynomial definition */
00958     
00959     if((poly_in=fopen(poly_table,"r"))==NULL) {
00960         cpl_msg_error("Cannot open file %s", poly_table);
00961         return NULL;
00962     }
00963     fscanf(poly_in,"%s %s %s %s",temp,temp,temp,temp);
00964     ch=fgetc(poly_in);
00965     while(ch!='p'){
00966         ch=fgetc(poly_in);
00967         poly_string[count]=ch;
00968         count++;
00969     }
00970     poly_string[count-2]='\0';
00971     /*fscanf(poly_in,"%s",poly_string);*/
00972     cpl_msg_info("%s",poly_string);
00973     
00974     fclose(poly_in);
00975     
00976     poly_u = poly2d_build_from_string(poly_string);
00977     
00978     /*poly_u = read_poly2d_from_table(poly_table) ;*/
00979         if (poly_u == NULL) {
00980             cpl_msg_error("image_warp","cannot read 2D poly from arc table") ;
00981             return NULL ;
00982         }
00983             poly_v=poly2d_build_from_string("0 1 1.0") ;
00984     
00985     warped = warp_image_generic(image, kernel_type, poly_u, poly_v);
00986     if (poly_u!=NULL) poly2d_free(poly_u);
00987     if (poly_v!=NULL) poly2d_free(poly_v);
00988     
00989     return warped;
00990 }
00991     
00992 #include <cpl_table.h> 
00993 OneImage * image_warp_fits( OneImage * image,
00994                       char      * kernel_type, 
00995                       char      * poly_table )
00996 {
00997         const char* _id="image_warp_fits";
00998         OneImage  * warped;
00999     char        poly_string[150];
01000     char        tmp_string[50];
01001         int             nraw=0;
01002     /* Following are for polynomial transforms */
01003     poly2d  *   poly_u;         /* polynomial definition */
01004     poly2d  *   poly_v;         /* polynomial definition */
01005     int*            degx;
01006     int*            degy;
01007     double*         coef;
01008         int             i=0;
01009         cpl_table* poly_tbl=NULL;
01010 
01011        if(is_fits_file(poly_table) != 1) {
01012          cpl_msg_error(_id,"Input file %s is not FITS",poly_table);
01013          return NULL;
01014        }
01015        poly_tbl = cpl_table_load(poly_table,1,0);
01016        nraw = cpl_table_get_nrow(poly_tbl);
01017        degx   = cpl_table_get_data_int(poly_tbl,"degx");
01018        degy   = cpl_table_get_data_int(poly_tbl,"degy");
01019        coef   = cpl_table_get_data_double(poly_tbl,"coeff");
01020 
01021     sprintf(poly_string,"%2d%2d%s%g ",degx[0],degy[0]," ",coef[0]);
01022         for(i=1;i<nraw;i++) {
01023       sprintf(tmp_string,"%2d%2d%s%g ",degx[i],degy[i]," ", coef[i]);
01024       strcat(poly_string,tmp_string);
01025         }
01026     /*fscanf(poly_in,"%s",poly_string);*/
01027     /* cpl_msg_info("%s",poly_string); */
01028     
01029     
01030     poly_u = poly2d_build_from_string(poly_string);
01031     
01032     /*poly_u = read_poly2d_from_table(poly_table) ;*/
01033         if (poly_u == NULL) {
01034             cpl_msg_error(_id,"cannot read 2D poly from arc table") ;
01035             return NULL ;
01036         }
01037         poly_v=poly2d_build_from_string("0 1 1.0") ;
01038         
01039     
01040     warped = warp_image_generic(image, kernel_type, poly_u, poly_v);
01041     if (poly_u!=NULL) poly2d_free(poly_u);
01042     if (poly_v!=NULL) poly2d_free(poly_v);
01043         cpl_table_delete(poly_tbl);
01044     return warped;
01045 }
01046     
01047  
01048 
01049 /*--------------------------------------------------------------------------*/
01050 
01051 /*--------------------------------------------------------------------------*/

Generated on Wed Oct 26 13:08:51 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001