detlin.c

00001 /*******************************************************************************
00002 * E.S.O. - VLT project
00003 *
00004 *
00005 *
00006 * who       when      what
00007 * --------  --------  ----------------------------------------------
00008 * schreib  18/04/02  created
00009 */
00010 
00011 /************************************************************************
00012 *   NAME
00013 *        detlin.c -
00014 *        procedures to fit the pixel linearity of the detector
00015 *        to search for bad pixels
00016 *
00017 *   SYNOPSIS
00018 *   #include "detlin.h"
00019 * 1) OneCube * fitIntensityCourse( OneCube * flatStack,
00020 *                                  int       order,
00021 *                                  float     loReject,
00022 *                                  float     hiReject )
00023 *
00024 * 2) OneImage * searchBadPixels( OneCube *  coeffs,
00025 *                                double     threshSigmaFactor,
00026 *                                double     nonlinearThresh,
00027 *                                float      loReject,
00028 *                                float      hiReject )
00029 *
00030 * 3) OneImage * searchBadPixelsViaNoise( OneCube *  darks,
00031 *                                        float      threshSigmaFactor,
00032 *                                        float      loReject,
00033 *                                        float      hiReject )
00034 *
00035 * 4) int countBadPixels ( OneImage * bad )
00036 *
00037 * 5) OneImage * meanImageInSpec( OneImage * im, float fmedian )
00038 *
00039 * 6) OneImage * absDistImage( OneImage * im, float fmedian )
00040 *
00041 * 7) OneImage * localMedianImage( OneImage * im, 
00042 *                                 float fmedian, 
00043 *                                 float loReject,
00044 *                                 float hiReject,
00045 *                                 int half_box_size )
00046 *
00047 * 1) this routine searches for static bad pixel positions
00048 *    by searching the fitted first and second polynomial coefficients 
00049 *    of each pixel response (linear) for outliers.
00050 *    Pixel with high non-linear response are also declared as bad.
00051 *
00052 * 2) this routine fits polynomials along the z-axis of 
00053 *    a data cube and stores the resulting coefficients 
00054 *    in a data cube. The eclipse routine fit_1d_poly()
00055 *    is used for polynomial fitting
00056 *    The input is assumed to be a cube containing flatfield frames of
00057 *    different intensities (usually increasing or decreasing).
00058 *    for each pixel position on the detector, a curve is plotted
00059 *    of the pixel intensity in each plane against the clean mean
00060 *    intensity of the plane. Then a polynomial of user defined
00061 *    order is fitted to the curves.
00062 *
00063 * 3) this routine searches for static bad pixel positions
00064 *    This is done by building a cube of dark frames and examine 
00065 *    the noise variations in each pixel. If big deviations
00066 *    from a clean mean pixel noise occurr, the pixel is 
00067 *    declared as bad.
00068 *
00069 * 4) this routine counts the number of bad pixels
00070 *
00071 * 5)  mean filter, calculates the mean for an image
00072 *     by using the 4 closest pixels of every pixel in spectral direction
00073 *     (column).
00074 *     The values in the output image are determined according
00075 *     to the values of the input parameter.
00076 *     If fmedian = 0: always replace by mean
00077 *     if fmedian < 0: replace by mean if |pixel - mean| >
00078 *                                  -fmedian
00079 *     if fmedian > 0: replace by mean (fmedian as a factor of
00080 *                     the square root of the mean itself)
00081 *                     if |pixel - mean| >= fmedian * sqrt ( mean )
00082 *                     This can be used to consider photon noise.
00083 *                     This considers a dependence of the differences on the
00084 *                     pixel values themselves.
00085 *      Notice       :       it is assumed that most of the 4 nearest neighbor pixels
00086 *                           are not bad pixels!
00087 *                           blank pixels are not replaced!
00088 *
00089 * 6)  filter, calculates the absolute distances
00090 *     of the nearest neighbors for an image
00091 *     by using the 8 closest pixels of every pixel.
00092 *     The values in the output image are determined according
00093 *     to the values of the input parameter.
00094 *     If fmedian = 0: always replace by abs. distances
00095 *     if fmedian < 0: replace by abs. distances if |median_dist - dist| >
00096 *                     -fmedian
00097 *     if fmedian > 0: replace by abs. distances (fmedian as a factor of
00098 *                     the square root of the distance itself)
00099 *                     if |median_dist - dist| >= fmedian * sqrt ( dist )
00100 *                     This can be used to consider photon noise.
00101 *                     This considers a dependence of the differences on the
00102 *                     pixel values themselves.
00103 *     Notice     :    it is assumed that most of the 8 nearest neighbor pixels
00104 *                     are not bad pixels!
00105 *                     blank pixels are not replaced!
00106 *
00107 *  7) filter, calculates the local stdev in a moving box
00108 *     Then it calculates the difference of the pixel to the median
00109 *     of the nearest neighbors
00110 *     by using the 8 closest pixels of every pixel.
00111 *     The values in the output image are determined according
00112 *     to the values of the input parameter.
00113 *     If fmedian = 0: always replace by median
00114 *     if fmedian < 0: replace median if |median_dist - dist| >
00115 *                                        fmedian * stdev 
00116 *     if fmedian > 0: replace by median (fmedian as a factor of
00117 *                     the square root of the median itself)
00118 *     if |pixel - median| >= fmedian * sqrt ( median )
00119 *     This can be used to consider photon noise.
00120 *     This considers a dependence of the differences on the
00121 *     pixel values themselves.
00122 *     Notice       :       it is assumed that most of the 8 nearest neighbor pixels
00123 *                        are not bad pixels!
00124 *                        blank pixels are not replaced!
00125 *
00126 *   FILES
00127 *
00128 *   ENVIRONMENT
00129 *
00130 *   RETURN VALUES
00131 *
00132 *   CAUTIONS
00133 *
00134 *   EXAMPLES
00135 *
00136 *   SEE ALSO
00137 *
00138 *   BUGS
00139 *
00140 *------------------------------------------------------------------------
00141 */
00142 
00143 #define POSIX_SOURCE 1
00144 #include "vltPort.h"
00145 
00146 /*
00147  * System Headers
00148  */
00149 
00150 /*
00151  * Local Headers
00152  */
00153 
00154 #include "detlin.h"
00155 
00156 /*----------------------------------------------------------------------------
00157  *                          Function codes
00158  *--------------------------------------------------------------------------*/
00159 
00160 /*----------------------------------------------------------------------------
00161    Function     :       fitIntensityCourse()
00162    In           :       flatStack: flats with in/decreasing intensity
00163                                    stacked in a data cube
00164                         order:     order of the fit polynomial
00165    Out          :       data cube containing the fit result. 
00166                         the first polynomial coefficients in the first plane
00167                         and so on. 
00168    Job          :       this routine fits polynomials along the z-axis of 
00169                         a data cube and stores the resulting coefficients 
00170                         in a data cube. The eclipse routine fit_1d_poly()
00171                         is used for polynomial fitting
00172                         The input is assumed to be a cube containing flatfield frames of
00173                         different intensities (usually increasing or decreasing).
00174                         for each pixel position on the detector, a curve is plotted
00175                         of the pixel intensity in each plane against the clean mean
00176                         intensity of the plane. Then a polynomial of user defined
00177                         order is fitted to the curves.
00178  ---------------------------------------------------------------------------*/
00179 
00180 OneCube * fitIntensityCourse( OneCube * flatStack,
00181                               int       order,
00182                               float     loReject,
00183                               float     hiReject )
00184 {
00185     OneCube * retCube ;
00186     dpoint  * points ;
00187     int       i, z ;
00188     double * coeffs ;
00189     Stats  * stats[flatStack->np] ;
00190   
00191     if ( NullCube == flatStack )
00192     {
00193         cpl_msg_error("fitIntensityCourse:"," no input cube given!\n") ;
00194         return NullCube ;
00195     }
00196     if ( order <= 0 )
00197     {
00198         cpl_msg_error("fitIntensityCourse:","wrong order of polynomial given!\n") ;
00199         return NullCube ;
00200     }
00201     /* allocate memory for returned cube */
00202     if ( NullCube == ( retCube = newCube(flatStack->lx, flatStack->ly, order+1) ) ) 
00203     {
00204         cpl_msg_error("fitIntensityCourse:","could not allocate memory!\n") ;
00205         return NullCube ;
00206     }
00207 
00208     for ( z = 0 ; z < flatStack->np ; z++ )
00209     {
00210         stats[z] = imageStatsOnRectangle(flatStack->plane[z], loReject, hiReject, 
00211                                          0, 0, flatStack->lx-1, flatStack->ly-1) ;
00212         if ( stats[z] == NULL )
00213         {
00214             cpl_msg_error("fitIntensityCourse:","could not compute image statistics in plane: %d\n", z) ;
00215             destroy_cube(retCube) ;
00216             return NullCube ;
00217         }
00218     }
00219 
00220     /* go through the image plane and store the spectra in a double points data structure */
00221     for ( i = 0 ; i < (int)flatStack->plane[0]->nbpix ; i++ )
00222     {
00223         /* allocate dpoint object */
00224         if ( NULL == ( points = (dpoint*) cpl_calloc(flatStack->np, sizeof(dpoint)) ) )
00225         {
00226             cpl_msg_error("fitIntensityCourse:","could not allocate memory!\n") ;
00227             destroy_cube(retCube) ;
00228             return NullCube ;
00229         }
00230         for ( z = 0 ; z < flatStack->np ; z++ )
00231         {
00232             points[z].x = (double)stats[z]->cleanmean ;
00233             points[z].y = (double)flatStack->plane[z]->data[i] ;
00234         }
00235         if ( NULL == ( coeffs = fit_1d_poly(order, points, flatStack->np, NULL) ) )
00236         {
00237             cpl_msg_warning("could not fit spectrum of pixel:","%d\n", i) ;
00238             for ( z = 0 ; z < order+1 ; z++ )
00239             {
00240                 retCube->plane[z]->data[i] = ZERO ;
00241             }
00242         }
00243         else
00244         {
00245             for ( z = 0 ; z < order+1 ; z++ )
00246             {
00247                 retCube->plane[z]->data[i] = coeffs[z] ;
00248             }
00249         }
00250         cpl_free(points) ;
00251         free(coeffs) ;
00252     }
00253     for ( z = 0 ; z < flatStack->np ; z++ )
00254     {
00255         cpl_free (stats[z]) ;
00256     }
00257 
00258     return retCube ;
00259 }
00260 
00261 /*----------------------------------------------------------------------------
00262    Function     :       searchBadPixels()
00263    In           :       coeffs: fitted polynomial coefficients
00264                                 stored in a cube
00265                         threshSigmaFactor: factor of determined standard deviation
00266                                            of zero and slope coefficients to determine the 
00267                                            threshold for good or bad pixels
00268                         nonlinearThresh:   absolute threshold value of the found
00269                                            non-linear polynomial coefficients beyond
00270                                            which the pixels are declared as bad.
00271                         loReject, hiReject: percentage (0...100) of extreme pixel 
00272                                             values that is not considered for image
00273                                             statistics
00274    Out          :       Bad pixel mask image (1: good pixel, 0: bad pixel). 
00275    Job          :       this routine searches for static bad pixel positions
00276                         by searching the fitted first and second polynomial coefficients 
00277                         of each pixel response (linear) for outliers.
00278                         Pixel with high non-linear response are also declared as bad.
00279  ---------------------------------------------------------------------------*/
00280 
00281 OneImage * searchBadPixels( OneCube *  coeffs,
00282                             double     threshSigmaFactor,
00283                             double     nonlinearThresh,
00284                             float      loReject,
00285                             float      hiReject )
00286 {
00287     OneImage * retImage ;
00288     int i, z ;
00289     Stats * stats ;
00290 
00291     if ( NullCube == coeffs )
00292     {
00293         cpl_msg_error("searchBadPixels:","no input cube given!\n") ;
00294         return NullImage ;
00295     }
00296     if ( threshSigmaFactor <= 0. )
00297     {
00298         cpl_msg_error("searchBadPixels:","wrong sigma factor given, 0 or negativ!\n") ;
00299         return NullImage ;
00300     }
00301     if ( nonlinearThresh <= 0. )
00302     {
00303         cpl_msg_error("searchBadPixels:","wrong nonlinear threshold value given, 0 or negativ!\n") ;
00304         return NullImage ;
00305     }
00306     if ( coeffs->np <= 1 )
00307     {
00308         cpl_msg_error("searchBadPixels:","no cube given, only one plane!\n") ;
00309         return NullImage ;
00310     }
00311     /* allocate memory for return image */
00312     if ( NullImage == (retImage = new_image(coeffs->lx, coeffs->ly)) )
00313     {
00314         cpl_msg_error("searchBadPixels:","could not allocate memory!\n") ;
00315         return NullImage ;
00316     }
00317 
00318     /* first test the sensitivity deviations of each pixel */
00319     /* determine the clean mean and clean standard deviation in the whole image frame */
00320     stats = imageStatsOnRectangle(coeffs->plane[1], loReject, hiReject, 0, 0, coeffs->lx-1, coeffs->ly-1) ;
00321     if ( NULL == stats )
00322     {
00323         cpl_msg_error("searchBadPixels:","could not determine image statistics!\n") ;
00324         destroy_image(retImage) ;
00325         return NullImage ;
00326     }
00327     for ( i = 0 ; i < (int) coeffs->plane[1]->nbpix ; i++ )
00328     {
00329         if ( qfits_isnan(coeffs->plane[1]->data[i]) )
00330         {
00331             retImage->data[i] = 0. ;
00332         }
00333         else if ( stats->cleanmean - coeffs->plane[1]->data[i] > threshSigmaFactor*stats->cleanstdev )
00334         {
00335             retImage->data[i] = 0. ;
00336         }
00337         else
00338         {
00339            retImage->data[i] = 1. ;
00340         }
00341     }
00342     cpl_free(stats) ;
00343 
00344 
00345     /* -----------------------------------------------------
00346      * now test additionally the non-linearity if available. 
00347      * if a strong non-linearity occurs for pixels which are 
00348      * declared "good" so far (normal linear coefficients)
00349      * these pixels will be declared bad.    
00350      */
00351     if (coeffs->np > 1) 
00352     {
00353         for ( z = 2 ; z < coeffs->np ; z++ )
00354         {
00355             stats = imageStatsOnRectangle(coeffs->plane[z], loReject, hiReject, 0, 0, coeffs->lx-1, coeffs->ly-1) ;
00356             if ( NULL == stats )
00357             {
00358                 cpl_msg_error("searchBadPixels:","could not determine image statistics!\n") ;
00359                 destroy_image(retImage) ;
00360                 return NullImage ;
00361             }
00362             for ( i = 0 ; i < (int) coeffs->plane[z]->nbpix ; i++ )
00363             {
00364                 if ( retImage->data[i] == 1. && 
00365                      (fabs(coeffs->plane[z]->data[i] - stats->cleanmean) > 
00366                                      threshSigmaFactor*stats->cleanstdev ||
00367               fabs(coeffs->plane[z]->data[i]) > nonlinearThresh ) ) 
00368                 {
00369                     retImage->data[i] = 0. ;
00370                 }
00371             }
00372             cpl_free(stats) ;
00373         }
00374     }
00375     
00376     return retImage ;
00377 }
00378 
00379 /*----------------------------------------------------------------------------
00380    Function     :       searchBadPixelsViaNoise()
00381    In           :       darks: sequence of darks (NDIT = 1) 
00382                                stored in a cube, at least 10 to get good statistics
00383                         threshSigmaFactor: factor to determined standard deviation
00384                                            in each pixel to determine the threshold
00385                                            beyond which a pixel is declared as bad. 
00386                         loReject, hiReject: percentage (0...100) of extreme pixel 
00387                                             values that is not considered for image
00388                                             statistics
00389    Out          :       Bad pixel mask image (1: good pixel, 0: bad pixel). 
00390    Job          :       this routine searches for static bad pixel positions
00391                         This is done by building a cube of dark frames and examine 
00392                         the noise variations in each pixel. If big deviations
00393                         from a clean mean pixel noise occurr, the pixel is 
00394                         declared as bad.
00395  ---------------------------------------------------------------------------*/
00396 
00397 OneImage * searchBadPixelsViaNoise( OneCube *  darks,
00398                                     float      threshSigmaFactor,
00399                                     float      loReject,
00400                                     float      hiReject )
00401 {
00402     OneImage * bad ;
00403     int        z, n, i ;
00404     int        lx, ly ;
00405     int        row, col ;
00406     int        low_n, high_n ;
00407     float    * spectrum ;
00408     double     pix_sum ;
00409     double     sqr_sum ;
00410     Stats    * stats ;
00411     
00412     if ( NullCube == darks )
00413     {
00414         cpl_msg_error("searchBadPixelsViaNoise:","no input cube given!\n") ;
00415         return NullImage ;
00416     }
00417     if ( threshSigmaFactor <= 0. )
00418     {
00419         cpl_msg_error("searchBadPixelsViaNoise:","factor is smaller or equal zero!\n") ;
00420         return NullImage ;
00421     }
00422     if ( loReject < 0. || hiReject < 0. || (loReject + hiReject) >= 100.  )
00423     {
00424         cpl_msg_error("searchBadPixelsViaNoise:","wrong reject percentage values!\n") ;
00425         return NullImage ;
00426     }
00427     if ( darks->np < 1 )
00428     {
00429         cpl_msg_error("rr searchBadPixelsViaNoise:","not enough dark frames given for good statistics!\n") ;
00430         return NullImage ;
00431     }
00432     lx = darks->plane[0]->lx ;
00433     ly = darks->plane[0]->ly ;
00434     low_n  = (int)(loReject/100. *(float)darks->np) ;
00435     high_n = (int)(hiReject/100. *(float)darks->np) ;
00436     if (NullImage == (bad = new_image (lx, ly) ) )
00437     {
00438         cpl_msg_error("searchBadPixelsViaNoise:","could not allocate new memory!\n") ;
00439         return NullImage ;
00440     }
00441     if (NULL == (spectrum = (float*) cpl_calloc(darks->np, sizeof(float)) ) )
00442     {
00443         cpl_msg_error("searchBadPixelsViaNoise:","could not allocate new memory!\n") ;
00444         return NullImage ;
00445     }
00446     for ( col = 0 ; col < lx ; col++ )
00447     {
00448         for ( row = 0 ; row < ly ; row++ )
00449         {
00450             for ( z = 0 ; z < darks->np ; z++ )
00451             {
00452                 spectrum[z] = darks->plane[z]->data[col+lx*row] ;
00453             }
00454             pixel_qsort(spectrum, darks->np) ;
00455             n = 0  ;
00456             pix_sum = 0.; 
00457             sqr_sum = 0.; 
00458             for ( i = low_n ; i < darks->np - high_n ; i++ )
00459             {
00460                 pix_sum += (double)spectrum[i] ;
00461                 sqr_sum += ((double)spectrum[i]*(double)spectrum[i]) ;
00462                 n++ ;
00463             }
00464             /* compute the noise in each pixel */
00465             pix_sum /= (double)n ;
00466             sqr_sum /= (double)n ;
00467             bad->data[col+lx*row] = (float)sqrt(sqr_sum - pix_sum*pix_sum) ;
00468         }
00469     }
00470     cpl_free(spectrum) ;
00471     if ( NULL == (stats = imageStatsOnRectangle (bad, loReject, hiReject, 200, 200, 800, 800) ) )
00472     {
00473         cpl_msg_error("searchBadPixelsViaNoise:","could not get image statistics!\n") ;
00474         destroy_image (bad) ;
00475         return NullImage ;
00476     }
00477  
00478 
00479     /* now build the bad pixel mask */
00480     for ( col = 0 ; col < lx ; col++ )
00481     {
00482         for ( row = 0 ; row < ly ; row++ )
00483         {
00484             if (bad->data[col+lx*row] > stats->cleanmean+threshSigmaFactor*stats->cleanstdev ||
00485                 bad->data[col+lx*row] < stats->cleanmean-threshSigmaFactor*stats->cleanstdev) 
00486             {
00487                 bad->data[col+lx*row] = 0. ;
00488             }
00489             else
00490             {
00491                 bad->data[col+lx*row] = 1. ;
00492             }
00493         }
00494     }
00495     cpl_free (stats) ;
00496     return bad ;
00497 }
00498 
00499 /*----------------------------------------------------------------------------
00500    Function     :       countBadPixels()
00501    In           :       bad: bad pixel mask
00502    Out          :       number of bad pixels. 
00503    Job          :       this routine counts the number of bad pixels
00504  ---------------------------------------------------------------------------*/
00505 
00506 int countBadPixels ( OneImage * bad )
00507 {
00508     int i, n ;
00509    
00510     n = 0 ;
00511     for ( i = 0 ; i < (int) bad->nbpix ; i++ )
00512     {
00513         if ( bad->data[i] == 0 || qfits_isnan(bad->data[i]) )
00514         {
00515             n++ ;
00516         }
00517     }
00518     return n ;
00519 }
00520 
00521 
00522 /*----------------------------------------------------------------------------
00523    Function     :       meanImageInSpec()
00524    In           :       image, a threshold parameter
00525    Out          :       resulting image
00526    Job          :       mean filter, calculates the mean for an image
00527                         by using the 4 closest pixels of every pixel in spectral direction
00528                         (column).
00529                         The values in the output image are determined according
00530                         to the values of the input parameter.
00531                         If fmedian = 0: always replace by mean
00532                         if fmedian < 0: replace by mean if |pixel - mean| >
00533                                         -fmedian
00534                         if fmedian > 0: replace by mean (fmedian as a factor of
00535                                         the square root of the mean itself)
00536                                         if |pixel - mean| >= fmedian * sqrt ( mean )
00537                                         This can be used to consider photon noise.
00538                                         This considers a dependence of the differences on the
00539                                         pixel values themselves.
00540    Notice       :       it is assumed that most of the 4 nearest neighbor pixels
00541                         are not bad pixels!
00542                         blank pixels are not replaced!
00543  ---------------------------------------------------------------------------*/
00544 
00545 OneImage * meanImageInSpec( OneImage * im, float fmedian )
00546 {
00547     OneImage *   image       ;
00548     pixelvalue * value       ;
00549     pixelvalue   mean      ;
00550     int        * position    ;
00551     int          nposition   ;
00552     int          n, i, j ;
00553 
00554     if ( im == NullImage )
00555     {
00556         cpl_msg_error ("meanImageInSpec:","no image input\n") ;
00557         return NullImage ;
00558     }
00559 
00560     image = copy_image ( im ) ;
00561 
00562     /*----------------------------------------------------------------------
00563      * go through all pixels
00564      */
00565 
00566     for ( i = 0 ; i < (int) im -> nbpix ; i++ )
00567     {
00568         /* blank pixels are not replaced */
00569         if ( qfits_isnan(im -> data[i]) )
00570         {
00571             continue ;
00572         }
00573 
00574         /* initialize the buffer variables for the 2 nearest spectral neighbors */
00575         value = (pixelvalue * )cpl_calloc ( 4, sizeof ( pixelvalue * ) ) ;
00576         position = ( int * ) cpl_calloc ( 4, sizeof ( int * ) ) ;
00577 
00578         /*----------------------------------------------------------------------
00579          * determine the pixel position of the 8 nearest neighbors
00580          */
00581 
00582         position[0] = i + im -> lx     ; /* upper       */
00583         position[1] = i + 2*im -> lx   ; /* upper       */
00584         position[2] = i - im -> lx     ; /* lower       */
00585         position[3] = i - 2*im -> lx   ; /* lower       */
00586 
00587         /*----------------------------------------------------------------------
00588          * determine the positions of the image margins, top positions are changed
00589          * to low positions and vice versa. Right positions are changed to left
00590          * positions and vice versa.
00591          */
00592 
00593         if ( i >= 0 && i < im -> lx )    /* bottom line */
00594         {
00595             position[2] += 2 * im -> lx ;
00596             position[3] += 4 * im -> lx ;
00597         }
00598         else if ( i >= ((int) im -> nbpix - im -> lx ) && i < (int) im -> nbpix ) /* top line */
00599         {
00600             position[0] -= 2 * im -> lx ;
00601             position[1] -= 4 * im -> lx ;
00602         }
00603 
00604         /* ----------------------------------------------------------------------------
00605          * read the pixel values of the neighboring pixels,
00606          * blanks are not considered
00607          */
00608 
00609         nposition = 4 ;
00610         n = 0 ;
00611         for ( j = 0 ; j < nposition ; j ++ )
00612         {
00613             if ( !qfits_isnan(im -> data[position[j]]) )
00614             {
00615                 value[n] = im -> data[position[j]] ;
00616                 n ++ ;
00617             }
00618         }
00619         nposition = n ;
00620 
00621         if ( nposition < 1 )  /* all neighbors are blank */
00622         {
00623             image->data[i] = ZERO ;
00624             cpl_free(value) ;
00625             cpl_free(position) ;
00626             continue ;
00627         }
00628 
00629         /* determine the mean */
00630         mean = 0. ;
00631         for ( n = 0 ; n < nposition ; n++ )
00632         {
00633             mean += value[n] ;
00634         }
00635         mean /= (float) nposition ;
00636 
00637         /* -----------------------------------------------------------------
00638          * replace the pixel value by the median on conditions:
00639          * fmedian = 0:","always replace with mean.
00640          * fmedian < 0: interpret as absolute condition:
00641          *              if |pixel - mean| > -fmedian
00642          *              replace with mean.
00643          */
00644 
00645         if ( fmedian == 0 )
00646         {
00647             image -> data[i] = mean ;
00648         }
00649         else if ( fmedian < 0 &&
00650                   fabs ( mean - im -> data[i] ) >= -fmedian )
00651         {
00652             image -> data[i] = mean ;
00653         }
00654         else if ( fmedian > 0 &&
00655                   fabs ( mean - im -> data[i] ) >= fmedian * sqrt(fabs(mean)) )
00656         {
00657             image -> data[i] = mean ;
00658         }
00659         else
00660         {
00661             cpl_free (value) ;
00662             cpl_free (position) ;
00663             continue ;
00664         }
00665 
00666         cpl_free (value) ;
00667         cpl_free (position) ;
00668     }
00669     return image ;
00670 }
00671 
00672 
00673 /*----------------------------------------------------------------------------
00674    Function     :       absDistImage()
00675    In           :       image, a threshold parameter
00676    Out          :       resulting image
00677    Job          :       filter, calculates the absolute distances
00678                         of the nearest neighbors for an image
00679                         by using the 8 closest pixels of every pixel.
00680                         The values in the output image are determined according
00681                         to the values of the input parameter.
00682                         If fmedian = 0: always replace by abs. distances
00683                         if fmedian < 0: replace by abs. distances if |median_dist - dist| >
00684                                         -fmedian
00685                         if fmedian > 0: replace by abs. distances (fmedian as a factor of
00686                                         the square root of the distance itself)
00687                                         if |median_dist - dist| >= fmedian * sqrt ( dist )
00688                                         This can be used to consider photon noise.
00689                                         This considers a dependence of the differences on the
00690                                         pixel values themselves.
00691    Notice       :       it is assumed that most of the 8 nearest neighbor pixels
00692                         are not bad pixels!
00693                         blank pixels are not replaced!
00694  ---------------------------------------------------------------------------*/
00695 
00696 OneImage * absDistImage( OneImage * im, float fmedian )
00697 {
00698     OneImage *   image       ;
00699     pixelvalue * value       ;
00700     pixelvalue   dist      ;
00701     pixelvalue   median_dist      ;
00702     pixelvalue   pix_dist[im->nbpix] ;
00703     int        * position    ;
00704     int          nposition   ;
00705     int          n, m, i, j ;
00706     double       sum, sum2 ;
00707     double       stdev ;
00708 
00709     if ( im == NullImage )
00710      {
00711         cpl_msg_error ("absDistImage:","no image input\n") ;
00712         return NullImage ;
00713     }
00714 
00715     image = copy_image ( im ) ;
00716 
00717     /*----------------------------------------------------------------------
00718      * go through all pixels
00719      */
00720 
00721     sum = 0. ;
00722     sum2 = 0. ;
00723     m = 0 ;
00724     for ( i = 0 ;  i < (int) im -> nbpix ; i++ )
00725     {
00726         /* blank pixels are not replaced */
00727         if ( qfits_isnan(im -> data[i]) )
00728         {
00729             continue ;
00730         }
00731 
00732         /* initialize the buffer variables for the 8 nearest neighbors */
00733         value = (pixelvalue * )cpl_calloc ( 8, sizeof ( pixelvalue * ) ) ;
00734         position = ( int * ) cpl_calloc ( 8, sizeof ( int * ) ) ;
00735 
00736         /*----------------------------------------------------------------------
00737          * determine the pixel position of the 8 nearest neighbors
00738          */
00739 
00740         position[0] = i + im -> lx - 1 ; /* upper left  */
00741         position[1] = i + im -> lx     ; /* upper       */
00742         position[2] = i + im -> lx + 1 ; /* upper right */
00743         position[3] = i + 1            ; /* right       */
00744         position[4] = i - im -> lx + 1 ; /* lower right */
00745         position[5] = i - im -> lx     ; /* lower       */
00746         position[6] = i - im -> lx - 1 ; /* lower left  */
00747         position[7] = i - 1            ; /* left        */
00748 
00749         /*----------------------------------------------------------------------
00750          * determine the positions of the image margins, top positions are changed
00751          * to low positions and vice versa. Right positions are changed to left
00752          * positions and vice versa.
00753          */
00754 
00755         if ( i >= 0 && i < im -> lx )    /* bottom line */
00756         {
00757             position[4] += 2 * im -> lx ;
00758             position[5] += 2 * im -> lx ;
00759             position[6] += 2 * im -> lx ;
00760         }
00761         else if ( i >= ((int) im -> nbpix - im -> lx ) && i < (int) im -> nbpix ) /* top line */
00762         {
00763             position[0] -= 2 * im -> lx ;
00764             position[1] -= 2 * im -> lx ;
00765             position[2] -= 2 * im -> lx ;
00766         }
00767         else if ( i % im -> lx == 0 )    /* left side */
00768         {
00769             position[0] += 2 ;
00770             position[6] += 2 ;
00771             position[7] += 2 ;
00772         }
00773         else if ( i % im -> lx == im -> lx - 1 )    /* right side */
00774         {
00775             position[2] -= 2 ;
00776             position[3] -= 2 ;
00777             position[4] -= 2 ;
00778         }
00779 
00780         /* ----------------------------------------------------------------------------
00781          * read the pixel values of the neighboring pixels,
00782          * blanks are not considered
00783          */
00784 
00785         nposition = 8 ;
00786         n = 0 ;
00787         for ( j = 0 ; j < nposition ; j ++ )
00788         {
00789             if ( !qfits_isnan(im -> data[position[j]]) )
00790             {
00791                 value[n] = im -> data[position[j]] ;
00792                 n ++ ;
00793             }
00794         }
00795         nposition = n ;
00796 
00797         if ( nposition <= 1 )  /* almost all neighbors are blank */
00798         {
00799             image->data[i] = ZERO ;
00800             cpl_free(value) ;
00801             cpl_free(position) ;
00802             continue ;
00803         }
00804 
00805         /* determine the absolute distances */
00806         dist = 0. ;
00807         for ( n = 0 ; n < nposition ; n++ )
00808         {
00809             dist += (im->data[i] - value[n])*(im->data[i] - value[n]) ;    
00810         }
00811         dist = sqrt(dist)/(float) nposition ;
00812         pix_dist[m] = dist ;
00813         m++ ;
00814         sum += (double)dist ;
00815         sum2 += (double)dist * (double)dist ;
00816         cpl_free(value) ;
00817         cpl_free(position) ;
00818     }
00819     sum /= (double)m ;
00820     sum2 /= (double)m ;
00821     stdev = sqrt(sum2 - sum*sum) ;
00822 
00823     median_dist = median(pix_dist, m) ;
00824 
00825     for ( i = 0 ; i < (int) im -> nbpix ; i++ )
00826     {
00827         /* blank pixels are not replaced */
00828         if ( qfits_isnan(im -> data[i]) )
00829         {
00830             continue ;
00831         }
00832 
00833         /* initialize the buffer variables for the 8 nearest neighbors */
00834         value = (pixelvalue * )cpl_calloc ( 8, sizeof ( pixelvalue * ) ) ;
00835         position = ( int * ) cpl_calloc ( 8, sizeof ( int * ) ) ;
00836 
00837         /*----------------------------------------------------------------------
00838          * determine the pixel position of the 8 nearest neighbors
00839          */
00840 
00841         position[0] = i + im -> lx - 1 ; /* upper left  */
00842         position[1] = i + im -> lx     ; /* upper       */
00843         position[2] = i + im -> lx + 1 ; /* upper right */
00844         position[3] = i + 1            ; /* right       */
00845         position[4] = i - im -> lx + 1 ; /* lower right */
00846         position[5] = i - im -> lx     ; /* lower       */
00847         position[6] = i - im -> lx - 1 ; /* lower left  */
00848         position[7] = i - 1            ; /* left        */
00849 
00850         /*----------------------------------------------------------------------
00851          * determine the positions of the image margins, top positions are changed
00852          * to low positions and vice versa. Right positions are changed to left
00853          * positions and vice versa.
00854          */
00855 
00856         if ( i >= 0 && i < im -> lx )    /* bottom line */
00857         {
00858             position[4] += 2 * im -> lx ;
00859             position[5] += 2 * im -> lx ;
00860             position[6] += 2 * im -> lx ;
00861         }
00862         else if ( i >= ((int) im -> nbpix - im -> lx ) && i < (int) im -> nbpix ) /* top line */
00863         {
00864             position[0] -= 2 * im -> lx ;
00865             position[1] -= 2 * im -> lx ;
00866             position[2] -= 2 * im -> lx ;
00867         }
00868         else if ( i % im -> lx == 0 )    /* left side */
00869         {
00870             position[0] += 2 ;
00871             position[6] += 2 ;
00872             position[7] += 2 ;
00873         }
00874         else if ( i % im -> lx == im -> lx - 1 )    /* right side */
00875         {
00876             position[2] -= 2 ;
00877             position[3] -= 2 ;
00878             position[4] -= 2 ;
00879         }
00880 
00881         /* ----------------------------------------------------------------------------
00882          * read the pixel values of the neighboring pixels,
00883          * blanks are not considered
00884          */
00885 
00886         nposition = 8 ;
00887         n = 0 ;
00888         for ( j = 0 ; j < nposition ; j ++ )
00889         {
00890             if ( !qfits_isnan(im -> data[position[j]]) )
00891             {
00892                 value[n] = im -> data[position[j]] ;
00893                 n ++ ;
00894             }
00895         }
00896         nposition = n ;
00897 
00898         if ( nposition <= 1 )  /* almost all neighbors are blank */
00899         {
00900             image->data[i] = ZERO ;
00901             cpl_free(value) ;
00902             cpl_free(position) ;
00903             continue ;
00904         }
00905 
00906         /* determine the absolute distances */
00907         dist = 0. ;
00908         for ( n = 0 ; n < nposition ; n++ )
00909         {
00910             dist += (im->data[i] - value[n])*(im->data[i] - value[n]) ;    
00911         }
00912         dist = sqrt(dist)/(float) nposition ;
00913 
00914 
00915         /* -----------------------------------------------------------------
00916          * replace the pixel value by the median on conditions:
00917          * fmedian = 0: always replace with median.
00918          * fmedian < 0: interpret as absolute condition:
00919          *              if |pixel - median| > -fmedian
00920          *              replace with median.
00921          * fmedian > 0: replace by median (fmedian as a factor of
00922          *              the square root of the median itself)
00923          *              if |pixel - median| >= fmedian * sqrt ( median )
00924          *              considers a dependence on the pixel value.
00925          *              This can be used to consider photon noise.
00926          */
00927 
00928         if ( fmedian == 0 )
00929         {
00930             image -> data[i] = dist ;
00931         }
00932         else if ( fmedian < 0 &&
00933                   fabs ( median_dist - dist ) >= -fmedian*stdev )
00934         {
00935             image -> data[i] = dist ;
00936         }
00937         else if ( fmedian > 0 &&
00938                   fabs ( median_dist - dist ) >= fmedian*stdev * sqrt(fabs(dist)) )
00939         {
00940             image -> data[i] = dist ;
00941         }
00942         else
00943         {
00944             cpl_free (value) ;
00945             cpl_free (position) ;
00946             continue ;
00947         }
00948 
00949         cpl_free (value) ;
00950         cpl_free (position) ;
00951     }
00952     return image ;
00953 }
00954 
00955 /*----------------------------------------------------------------------------
00956    Function     :       localMedianImage()
00957    In           :       im: input image
00958                         fmedian:  a factor to the local standard deviation
00959                         loReject, hiReject: fraction of rejected values to determine
00960                                             a clean standard deviation
00961                         half_box_size: integer half size of the running box to determine
00962                                        the local clean standard deviation
00963    Out          :       resulting image
00964    Job          :       filter, calculates the local stdev in a moving box
00965                         Then it calculates the difference of the pixel to the median
00966                         of the nearest neighbors
00967                         by using the 8 closest pixels of every pixel.
00968                         The values in the output image are determined according
00969                         to the values of the input parameter.
00970                         If fmedian = 0: always replace by median
00971                         if fmedian < 0: replace median if |median_dist - dist| >
00972                                         fmedian * stdev 
00973                         if fmedian > 0: replace by median (fmedian as a factor of
00974                                         the square root of the median itself)
00975                                         if |pixel - median| >= fmedian * sqrt ( median )
00976                                         This can be used to consider photon noise.
00977                                         This considers a dependence of the differences on the
00978                                         pixel values themselves.
00979    Notice       :       it is assumed that most of the 8 nearest neighbor pixels
00980                         are not bad pixels!
00981                         blank pixels are not replaced!
00982  ---------------------------------------------------------------------------*/
00983 
00984 OneImage * localMedianImage( OneImage * im, 
00985                              float fmedian, 
00986                              float loReject,
00987                              float hiReject,
00988                              int half_box_size )
00989 {
00990     OneImage *   image       ;
00991     pixelvalue * value       ;
00992     pixelvalue   median      ;
00993     int        * position    ;
00994     int          nposition   ;
00995     int          n, i, j ;
00996     int          llx, lly, urx, ury ;
00997     Stats *      stats ;
00998 
00999     if ( im == NullImage )
01000     {
01001         cpl_msg_error ("localMedianImage:","no image input\n") ;
01002         return NullImage ;
01003     }
01004     if ( half_box_size < 0 )
01005     {
01006         cpl_msg_error ("localMedianImage:","negativ box_size given\n") ;
01007         return NullImage ;
01008     }
01009 
01010     image = copy_image ( im ) ;
01011 
01012     /*----------------------------------------------------------------------
01013      * go through all pixels
01014      */
01015 
01016     for ( i = 0 ; i < (int) im -> nbpix ; i++ )
01017     {
01018         /* blank pixels are not replaced */
01019         if ( qfits_isnan(im -> data[i]) )
01020         {
01021             continue ;
01022         }
01023 
01024         /* compute the image statistics in the box area */
01025         llx = i%im->lx - half_box_size ; 
01026         if ( llx < 0 ) llx = 0 ;
01027         lly = i%im->ly - half_box_size ;
01028         if ( lly < 0 ) lly = 0 ;
01029         urx = i%im->lx + half_box_size ; 
01030         if ( urx >= im->lx ) urx = im->lx - 1 ;
01031         ury = i%im->ly + half_box_size ;
01032         if ( ury >= im->ly ) ury = im->ly - 1 ;
01033 
01034         if ( NULL == (stats = imageStatsOnRectangle (im, loReject, hiReject, llx, lly, urx, ury)) ) 
01035         {
01036             cpl_msg_warning ("localMedianImage:","could not determine image statistics in pixel %d\n", i) ;
01037             continue ;
01038         }
01039 
01040         /* initialize the buffer variables for the 8 nearest neighbors */
01041         value = (pixelvalue * )cpl_calloc ( 8, sizeof ( pixelvalue * ) ) ;
01042         position = ( int * ) cpl_calloc ( 8, sizeof ( int * ) ) ;
01043 
01044         /*----------------------------------------------------------------------
01045          * determine the pixel position of the 8 nearest neighbors
01046          */
01047 
01048         position[0] = i + im -> lx - 1 ; /* upper left  */
01049         position[1] = i + im -> lx     ; /* upper       */
01050         position[2] = i + im -> lx + 1 ; /* upper right */
01051         position[3] = i + 1            ; /* right       */
01052         position[4] = i - im -> lx + 1 ; /* lower right */
01053         position[5] = i - im -> lx     ; /* lower       */
01054         position[6] = i - im -> lx - 1 ; /* lower left  */
01055         position[7] = i - 1            ; /* left        */
01056 
01057         /*----------------------------------------------------------------------
01058          * determine the positions of the image margins, top positions are changed
01059          * to low positions and vice versa. Right positions are changed to left
01060          * positions and vice versa.
01061          */
01062 
01063         if ( i >= 0 && i < im -> lx )    /* bottom line */
01064         {
01065             position[4] += 2 * im -> lx ;
01066             position[5] += 2 * im -> lx ;
01067             position[6] += 2 * im -> lx ;
01068         }
01069         else if ( i >= ((int) im -> nbpix - im -> lx ) && i < (int) im -> nbpix ) /* top line */
01070         {
01071             position[0] -= 2 * im -> lx ;
01072             position[1] -= 2 * im -> lx ;
01073             position[2] -= 2 * im -> lx ;
01074         }
01075         else if ( i % im -> lx == 0 )    /* left side */
01076         {
01077             position[0] += 2 ;
01078             position[6] += 2 ;
01079             position[7] += 2 ;
01080         }
01081         else if ( i % im -> lx == im -> lx - 1 )    /* right side */
01082         {
01083             position[2] -= 2 ;
01084             position[3] -= 2 ;
01085             position[4] -= 2 ;
01086         }
01087 
01088         /* ----------------------------------------------------------------------------
01089          * read the pixel values of the neighboring pixels,
01090          * blanks are not considered
01091          */
01092 
01093         nposition = 8 ;
01094         n = 0 ;
01095         for ( j = 0 ; j < nposition ; j ++ )
01096         {
01097             if ( !qfits_isnan(im -> data[position[j]]) )
01098             {
01099                 value[n] = im -> data[position[j]] ;
01100                 n ++ ;
01101             }
01102         }
01103         nposition = n ;
01104 
01105         if ( nposition <= 1 )  /* almost all neighbors are blank */
01106         {
01107             image->data[i] = ZERO ;
01108             cpl_free(value) ;
01109             cpl_free(position) ;
01110             cpl_free(stats) ;
01111             continue ;
01112         }
01113 
01114         /* sort the values and determine the median */
01115 
01116         pixel_qsort( value, nposition ) ;
01117         if ( nposition % 2 == 1 )
01118         {
01119             median = value [ nposition/2 ] ;
01120         }
01121         else
01122         {
01123             median = ( value [nposition/2 - 1] + value [nposition/2] ) / 2. ;
01124         }
01125 
01126         /* -----------------------------------------------------------------
01127          * replace the pixel value by the median on conditions:
01128          * fmedian = 0: always replace with median.
01129          * fmedian > 0: replace by median (fmedian as a factor of
01130          *              the square root of the median itself)
01131          *              if |pixel - median| >= fmedian * sqrt ( median )
01132          *              considers a dependence on the pixel value.
01133          *              This can be used to consider photon noise.
01134          */
01135 
01136         if ( fmedian == 0 )
01137         {
01138             image -> data[i] = median ;
01139         }
01140         else if ( fmedian < 0 &&
01141                   fabs ( median - im -> data[i] ) >= -fmedian * stats->cleanstdev)
01142         {
01143             image -> data[i] = median ;
01144         }
01145         else if ( fmedian > 0 &&
01146                   fabs ( median - im -> data[i] ) >= fmedian * sqrt(fabs(median)) )
01147         {
01148             image -> data[i] = median ;
01149         }
01150         else
01151         {
01152             cpl_free (value) ;
01153             cpl_free (position) ;
01154             cpl_free (stats) ;
01155             continue ;
01156         }
01157 
01158         cpl_free (value) ;
01159         cpl_free (position) ;
01160         cpl_free (stats) ;
01161     }
01162     return image ;
01163 }
01164 
01165 /*----------------------------------------------------------------------------
01166    Function     :       localMedianImage2()
01167    In           :       im: input image
01168                         fmedian:  a factor to the local standard deviation
01169                         loReject, hiReject: fraction of rejected values to determine
01170                                             a clean standard deviation
01171                         half_box_size: integer half size of the running box to determine
01172                                        the local clean standard deviation
01173    Out          :       resulting image
01174    Job          :       filter, calculates the local stdev in a moving box
01175                         Then it calculates the difference of the pixel to the 
01176                         clean mean of the moving box.
01177                         The values in the output image are determined according
01178                         to the values of the input parameter.
01179                         If fmedian = 0: always replace by clean mean
01180                         if fmedian < 0: replace clean mean if |cleanmean - dist| >
01181                                         fmedian * stdev 
01182                         if fmedian > 0: replace by clean mean (fmedian as a factor of
01183                                         the square root of the clean mean itself)
01184                                         if |pixel - cleanmean| >= fmedian * sqrt ( cleanmean )
01185                                         This can be used to consider photon noise.
01186                                         This considers a dependence of the differences on the
01187                                         pixel values themselves.
01188    Notice       :       it is assumed that most of the 8 nearest neighbor pixels
01189                         are not bad pixels!
01190                         blank pixels are not replaced!
01191  ---------------------------------------------------------------------------*/
01192 
01193 OneImage * localMedianImage2( OneImage * im, 
01194                              float fmedian, 
01195                              float loReject,
01196                              float hiReject,
01197                              int half_box_size )
01198 {
01199     OneImage *   image       ;
01200     pixelvalue   median      ;
01201     /*int          nposition   ;*/
01202     int          /*n,*/ i/*, j*/ ;
01203     int          llx, lly, urx, ury ;
01204     Stats *      stats ;
01205 
01206     if ( im == NullImage )
01207     {
01208         cpl_msg_error ("localMedianImage:","no image input\n") ;
01209         return NullImage ;
01210     }
01211     if ( half_box_size < 0 )
01212     {
01213         cpl_msg_error ("localMedianImage:","negativ box_size given\n") ;
01214         return NullImage ;
01215     }
01216 
01217     image = copy_image ( im ) ;
01218 
01219     /*----------------------------------------------------------------------
01220      * go through all pixels
01221      */
01222 
01223     for ( i = 0 ; i < (int) im -> nbpix ; i++ )
01224     {
01225         /* blank pixels are not replaced */
01226         if ( qfits_isnan(im -> data[i]) )
01227         {
01228             continue ;
01229         }
01230 
01231         /* compute the image statistics in the box area */
01232         llx = i%im->lx - half_box_size ; 
01233         if ( llx < 0 ) llx = 0 ;
01234         lly = i%im->ly - half_box_size ;
01235         if ( lly < 0 ) lly = 0 ;
01236         urx = i%im->lx + half_box_size ; 
01237         if ( urx >= im->lx ) urx = im->lx - 1 ;
01238         ury = i%im->ly + half_box_size ;
01239         if ( ury >= im->ly ) ury = im->ly - 1 ;
01240 
01241         if ( NULL == (stats = imageStatsOnRectangle (im, loReject, hiReject, llx, lly, urx, ury)) ) 
01242         {
01243             cpl_msg_warning ("localMedianImage:","could not determine image statistics in pixel %d\n", i) ;
01244             continue ;
01245         }
01246 
01247 
01248         /* -----------------------------------------------------------------
01249          * replace the pixel value by the median on conditions:
01250          * fmedian = 0: always replace with median.
01251          * fmedian > 0: replace by median (fmedian as a factor of
01252          *              the square root of the median itself)
01253          *              if |pixel - median| >= fmedian * sqrt ( median )
01254          *              considers a dependence on the pixel value.
01255          *              This can be used to consider photon noise.
01256          */
01257 
01258         if ( fmedian == 0 )
01259         {
01260             image -> data[i] = stats->cleanmean ;
01261         }
01262         else if ( fmedian < 0 &&
01263                   fabs ( stats->cleanmean - im -> data[i] ) >= -fmedian * stats->cleanstdev)
01264         {
01265             image -> data[i] = stats->cleanmean ;
01266         }
01267         else if ( fmedian > 0 &&
01268                   fabs ( stats->cleanmean - im -> data[i] ) >= fmedian * sqrt(fabs(median)) )
01269         {
01270             image -> data[i] = stats->cleanmean ;
01271         }
01272         else
01273         {
01274             cpl_free (stats) ;
01275             continue ;
01276         }
01277 
01278         cpl_free (stats) ;
01279     }
01280     return image ;
01281 }
01282 
01283 
01284 /*--------------------------------------------------------------------------*/

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