/* * This file is part of the ESO SINFONI Pipeline * Copyright (C) 2004,2005 European Southern Observatory * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, 51 Franklin St, Fifth Floor, Boston, MA 02111-1307 USA */ /******************************************************************************* * E.S.O. - VLT project * * * * who when what * -------- -------- ---------------------------------------------- * amodigli 18/04/02 created */ #ifdef HAVE_CONFIG_H # include #endif #define POSIX_SOURCE 1 #include "sinfo_vltPort.h" /* * System Headers */ /* * Local Headers */ #include "sinfo_detlin.h" #include "sinfo_recipes.h" #include "sinfo_fit_curve.h" /**@{*/ /** * @defgroup sinfo_detlin Detector Linearity Determination Functions * * TBD */ /*---------------------------------------------------------------------------- * Function codes *--------------------------------------------------------------------------*/ /** @brief this routine fits polynomials along the z-axis of a data cube and stores the resulting coefficients in a data cube. @param flatStack: flats with in/decreasing intensity stacked in a data cube @param order: order of the fit polynomial @return data cube containing the fit result. the first polynomial coefficients in the first plane and so on. @doc this routine fits polynomials along the z-axis of a data cube and stores the resulting coefficients in a data cube. The eclipse routine sinfo_fit_1d_poly() is used for polynomial fitting The input is assumed to be a cube containing flatfield frames of different intensities (usually increasing or decreasing). for each pixel position on the detector, a curve is plotted of the pixel intensity in each plane against the clean mean intensity of the plane. Then a polynomial of user defined order is fitted to the curves. */ cpl_imagelist * sinfo_new_fit_intensity_course(cpl_imagelist * flatStack, int order, float loReject, float hiReject ) { cpl_imagelist * ret_iml ; int i, z ; double * coeffs ; Stats ** stats=NULL ; int sx; int sy; int sz; float* psrcdata; float* presdata; cpl_image* img_tmp=NULL; sx=cpl_image_get_size_x(cpl_imagelist_get(flatStack,0)); sy=cpl_image_get_size_y(cpl_imagelist_get(flatStack,0)); sz=cpl_imagelist_get_size(flatStack); stats=(Stats**) cpl_calloc(sz,sizeof(Stats*)) ; if ( NULL == flatStack ) { sinfo_msg_error("no input cube given!") ; cpl_free(stats); return NULL ; } if ( order <= 0 ) { sinfo_msg_error("wrong order of polynomial given!") ; cpl_free(stats); return NULL ; } /* allocate memory for returned cube */ ret_iml = cpl_imagelist_new(); for ( z = 0 ; z < order+1 ; z++ ) { img_tmp=cpl_image_new(sx,sy,CPL_TYPE_FLOAT); cpl_imagelist_set(ret_iml,img_tmp,z); } for ( z = 0 ; z < sz ; z++ ) { stats[z]= sinfo_new_image_stats_on_rectangle(cpl_imagelist_get(flatStack,z), loReject, hiReject, 0, 0, sx-1, sy-1) ; if ( stats[z] == NULL ) { sinfo_msg_error("could not compute image statistics " "in plane: %d", z) ; cpl_imagelist_delete(ret_iml) ; return NULL ; } } /* go through the image plane and store the spectra in a double points data structure */ for ( i = 0 ; i < sx*sy ; i++ ) { /* allocate dpoint object */ dpoint * points ; if ( NULL == ( points = (dpoint*) cpl_calloc(sz, sizeof(dpoint)) ) ) { sinfo_msg_error("could not allocate memory!\n") ; cpl_imagelist_delete(ret_iml) ; return NULL ; } for ( z = 0 ; z < sz ; z++ ) { if(NULL==(img_tmp = cpl_imagelist_get(flatStack,z))) { sinfo_msg_error("could not get image!"); cpl_imagelist_delete(ret_iml) ; cpl_free(points) ; return NULL; } else { psrcdata=cpl_image_get_data_float(img_tmp); points[z].x = (double)stats[z]->cleanmean ; points[z].y = (double)psrcdata[i] ; } } if ( NULL == ( coeffs = sinfo_fit_1d_poly(order, points, sz, NULL) ) ) { sinfo_msg_warning("could not fit spectrum of pixel: %d\n", i) ; for ( z = 0 ; z < order+1 ; z++ ) { presdata=cpl_image_get_data_float(cpl_imagelist_get(ret_iml,z)); presdata[i] = ZERO ; } } else { for ( z = 0 ; z < order+1 ; z++ ) { if(NULL==(img_tmp = cpl_imagelist_get(ret_iml,z))) { sinfo_msg_error("could not get image!"); cpl_imagelist_delete(ret_iml) ; return NULL; } else { presdata=cpl_image_get_data_float(img_tmp); presdata[i] = coeffs[z] ; } } } cpl_free(points) ; cpl_free(coeffs) ; } for ( z = 0 ; z < sz ; z++ ) { cpl_free (stats[z]) ; } cpl_free(stats); return ret_iml ; } /** @brief search bad pixels @param coeffs: fitted polynomial coefficients stored in a cube @param threshSigmaFactor: factor of determined standard deviation of zero and slope coefficients to determine the threshold for good or bad pixels @param nonlinearThresh: absolute threshold value of the found non-linear polynomial coefficients beyond which the pixels are declared as bad. @param loReject, @param hiReject: percentage (0...100) of extreme pixel values that is not considered for image statistics @return Bad pixel mask image (1: good pixel, 0: bad pixel). Job this routine searches for static bad pixel positions by searching the fitted first and second polynomial coefficients of each pixel response (linear) for outliers. Pixel with high non-linear response are also declared as bad. */ cpl_image * sinfo_new_search_bad_pixels( cpl_imagelist * coeffs, double threshSigmaFactor, double nonlinearThresh, float loReject, float hiReject ) { int i ; Stats * stats ; int sx=0; int sy=0; int sz=0; cpl_image * img_res ; cpl_image* img_src=NULL; float* psrcdata=NULL; float* presdata=NULL; if ( NULL == coeffs ) { sinfo_msg_error("no input cube given!\n") ; return NULL ; } if ( threshSigmaFactor <= 0. ) { sinfo_msg_error("wrong sigma factor given, 0 or negativ!\n") ; return NULL ; } if ( nonlinearThresh <= 0. ) { sinfo_msg_error("wrong nonlinear threshold value given, " "0 or negative!") ; return NULL ; } sz=cpl_imagelist_get_size(coeffs); if ( sz <= 1 ) { sinfo_msg_error("no cube given, only one plane!\n") ; return NULL ; } /* Note that we refer to image #1! */ img_src=cpl_imagelist_get(coeffs,1); sx=cpl_image_get_size_x(img_src); sy=cpl_image_get_size_y(img_src); /* allocate memory for return image */ if ( NULL == (img_res = cpl_image_new(sx, sy,CPL_TYPE_FLOAT)) ) { sinfo_msg_error("could not allocate memory!\n") ; return NULL ; } /* first test the sensitivity deviations of each pixel */ /* determine the clean mean and clean standard deviation in the whole image frame */ stats = sinfo_new_image_stats_on_rectangle(img_src, loReject, hiReject, 0, 0, sx-1, sy-1) ; if ( NULL == stats ) { sinfo_msg_error("could not determine image statistics!\n") ; cpl_image_delete(img_res) ; return NULL ; } psrcdata=cpl_image_get_data_float(img_src); presdata=cpl_image_get_data_float(img_res); for ( i = 0 ; i < (int) sx*sy ; i++ ) { if ( isnan(psrcdata[i]) ) { presdata[i] = 0. ; } else if ( stats->cleanmean - psrcdata[i] > threshSigmaFactor*stats->cleanstdev ) { presdata[i] = 0. ; } else { presdata[i] = 1. ; } } cpl_free(stats) ; /* ----------------------------------------------------- * now test additionally the non-linearity if available. * if a strong non-linearity occurs for pixels which are * declared "good" so far (normal linear coefficients) * these pixels will be declared bad. */ if (sz > 1) { int z; for ( z = 2 ; z < sz ; z++ ) { img_src=cpl_imagelist_get(coeffs,z); sx=cpl_image_get_size_x(img_src); sy=cpl_image_get_size_y(img_src); psrcdata=cpl_image_get_data_float(img_src); stats = sinfo_new_image_stats_on_rectangle(img_src, loReject, hiReject, 0, 0, sx-1, sy-1) ; if ( NULL == stats ) { sinfo_msg_error("could not determine image statistics!\n") ; cpl_image_delete(img_res) ; return NULL ; } presdata=cpl_image_get_data_float(img_res); for ( i = 0 ; i < (int) sx*sy ; i++ ) { if ( presdata[i] == 1. && (fabs(psrcdata[i] - stats->cleanmean) > threshSigmaFactor*stats->cleanstdev || fabs(psrcdata[i]) > nonlinearThresh ) ) { presdata[i] = 0. ; } } cpl_free(stats) ; } } return img_res ; } /** @brief search bad pixels @param darks: sequence of darks (NDIT = 1) stored in a cube, at least 10 to get good statistics @param threshSigmaFactor: factor to determined standard deviation in each pixel to determine the threshold beyond which a pixel is declared as bad. @param loReject, @param hiReject: percentage (0...100) of extreme pixel values that is \ not considered for image statistics @return Bad pixel mask image (1: good pixel, 0: bad pixel). Job this routine searches for static bad pixel positions This is done by building a cube of sinfo_dark frames and examine the noise variations in each pixel. If big deviations from a clean mean pixel noise occurr, the pixel is declared as bad. */ cpl_image * sinfo_new_search_bad_pixels_via_noise(cpl_imagelist * darks, float threshSigmaFactor, float loReject, float hiReject ) { cpl_image * bp_map ; int z, n, i ; int lx, ly ; int row, col ; int low_n, high_n ; float * spectrum ; double pix_sum ; double sqr_sum ; Stats * stats ; cpl_image* img_src=NULL; float* psrcdata=NULL; float* pbpdata=NULL; int lz=0; if ( NULL == darks ) { sinfo_msg_error("no input cube given!\n") ; return NULL ; } if ( threshSigmaFactor <= 0. ) { sinfo_msg_error("factor is smaller or equal zero!\n") ; return NULL ; } if ( loReject < 0. || hiReject < 0. || (loReject + hiReject) >= 100. ) { sinfo_msg_error("wrong reject percentage values!\n") ; return NULL ; } lz=cpl_imagelist_get_size(darks); if ( lz < 1 ) { sinfo_msg_error("not enough dark frames given for good statistics!") ; return NULL ; } img_src=cpl_imagelist_get(darks,0); lx = cpl_image_get_size_x(img_src) ; ly = cpl_image_get_size_y(img_src) ; low_n = (int)(loReject/100. *(float)lz) ; high_n = (int)(hiReject/100. *(float)lz) ; if (NULL == (bp_map = cpl_image_new (lx, ly,CPL_TYPE_FLOAT) ) ) { sinfo_msg_error("could not allocate new memory!\n") ; return NULL ; } pbpdata=cpl_image_get_data(bp_map); if (NULL == (spectrum = (float*) cpl_calloc(lz, sizeof(float)) ) ) { sinfo_msg_error("could not allocate new memory!\n") ; return NULL ; } for ( row = 0 ; row < ly ; row++ ) { for ( col = 0 ; col < lx ; col++ ) { for ( z = 0 ; z < lz ; z++ ) { img_src=cpl_imagelist_get(darks,z); psrcdata=cpl_image_get_data(img_src); spectrum[z] = psrcdata[col+lx*row] ; } sinfo_pixel_qsort(spectrum, lz) ; n = 0 ; pix_sum = 0.; sqr_sum = 0.; for ( i = low_n ; i < lz - high_n ; i++ ) { pix_sum += (double)spectrum[i] ; sqr_sum += ((double)spectrum[i]*(double)spectrum[i]) ; n++ ; } /* compute the noise in each pixel */ pix_sum /= (double)n ; sqr_sum /= (double)n ; pbpdata[col+lx*row] = (float)sqrt(sqr_sum - pix_sum*pix_sum) ; } } cpl_free(spectrum) ; if ( NULL == (stats = sinfo_new_image_stats_on_rectangle (bp_map, loReject, hiReject, 200, 200, 800, 800) ) ) { sinfo_msg_error("could not get image statistics!\n") ; cpl_image_delete (bp_map) ; return NULL ; } /* now build the bad pixel mask */ for ( row = 0 ; row < ly ; row++ ) { for ( col = 0 ; col < lx ; col++ ) { if (pbpdata[col+lx*row] > stats->cleanmean+threshSigmaFactor*stats->cleanstdev || pbpdata[col+lx*row] < stats->cleanmean-threshSigmaFactor*stats->cleanstdev) { pbpdata[col+lx*row] = 0. ; } else { pbpdata[col+lx*row] = 1. ; } } } cpl_free (stats) ; return bp_map ; } /** @brief this routine counts the number of bad pixels @param bad: bad pixel mask @return number of bad pixels. */ int sinfo_new_count_bad_pixels (cpl_image * bad ) { int i, n ; int sx=cpl_image_get_size_x(bad); int sy=cpl_image_get_size_y(bad); float* pbpdata=cpl_image_get_data(bad); n = 0 ; for ( i = 0 ; i < (int) sx*sy ; i++ ) { if ( pbpdata[i] == 0 || isnan(pbpdata[i]) ) { n++ ; } } return n ; } /** @brief filter, calculates the absolute distances of the nearest neighbors for an image by using the 8 closest pixels of every pixel. @param image, a threshold parameter @return resulting image @doc filter, calculates the absolute distances of the nearest neighbors for an image by using the 8 closest pixels of every pixel. The values in the output image are determined according to the values of the input parameter. If fmedian = 0: always replace by abs. distances if fmedian < 0: replace by abs. distances if |median_dist - dist| > -fmedian if fmedian > 0: replace by abs. distances (fmedian as a factor of the square root of the distance itself) if |median_dist - dist| >= fmedian * sqrt ( dist ) This can be used to consider photon noise. This considers a dependence of the differences on the pixel values themselves. @Note it is assumed that most of the 8 nearest neighbor pixels are not bad pixels! blank pixels are not replaced! */ cpl_image * sinfo_new_abs_dist_image(cpl_image * im, float fmedian ) { cpl_image * image ; pixelvalue * value ; pixelvalue dist ; pixelvalue median_dist ; pixelvalue* pix_dist=NULL ; int * position ; int nposition ; int n, m, i, j ; double sum, sum2 ; double stdev ; float* pdata=NULL; int lx=0; int ly=0; if ( im == NULL ) { sinfo_msg_error ("no image input\n") ; return NULL ; } image = cpl_image_duplicate ( im ) ; /*---------------------------------------------------------------------- * go through all pixels */ sum = 0. ; sum2 = 0. ; m = 0 ; pdata = cpl_image_get_data(im); lx=cpl_image_get_size_x(im); ly=cpl_image_get_size_y(im); pix_dist=(pixelvalue*)cpl_calloc(lx*ly,sizeof(pixelvalue)) ; for ( i = 0 ; i < (int) lx*ly ; i++ ) { /* blank pixels are not replaced */ if ( isnan(pdata[i]) ) { continue ; } /* initialize the buffer variables for the 8 nearest neighbors */ value = (pixelvalue * )cpl_calloc ( 8, sizeof ( pixelvalue * ) ) ; position = ( int * ) cpl_calloc ( 8, sizeof ( int * ) ) ; /*-------------------------------------------------------------------- * determine the pixel position of the 8 nearest neighbors */ position[0] = i + lx - 1 ; /* upper left */ position[1] = i + lx ; /* upper */ position[2] = i + lx + 1 ; /* upper right */ position[3] = i + 1 ; /* right */ position[4] = i - lx + 1 ; /* lower right */ position[5] = i - lx ; /* lower */ position[6] = i - lx - 1 ; /* lower left */ position[7] = i - 1 ; /* left */ /*-------------------------------------------------------------------- * determine the positions of the image margins, top positions are * changed to low positions and vice versa. Right positions are * changed to left positions and vice versa. */ if ( i >= 0 && i < lx ) /* bottom line */ { position[4] += 2 * lx ; position[5] += 2 * lx ; position[6] += 2 * lx ; } else if ( i >= ((int) lx*ly - lx ) && i < (int) lx*ly ) /* top line */ { position[0] -= 2 * lx ; position[1] -= 2 * lx ; position[2] -= 2 * lx ; } else if ( i % lx == 0 ) /* left side */ { position[0] += 2 ; position[6] += 2 ; position[7] += 2 ; } else if ( i % lx == lx - 1 ) /* right side */ { position[2] -= 2 ; position[3] -= 2 ; position[4] -= 2 ; } /* ------------------------------------------------------------------- * read the pixel values of the neighboring pixels, * blanks are not considered */ nposition = 8 ; n = 0 ; for ( j = 0 ; j < nposition ; j ++ ) { if ( !isnan(pdata[position[j]]) ) { value[n] = pdata[position[j]] ; n ++ ; } } nposition = n ; if ( nposition <= 1 ) /* almost all neighbors are blank */ { pdata[i] = ZERO ; cpl_free(value) ; cpl_free(position) ; continue ; } /* determine the absolute distances */ dist = 0. ; for ( n = 0 ; n < nposition ; n++ ) { dist += (pdata[i] - value[n])*(pdata[i] - value[n]) ; } dist = sqrt(dist)/(float) nposition ; pix_dist[m] = dist ; m++ ; sum += (double)dist ; sum2 += (double)dist * (double)dist ; cpl_free(value) ; cpl_free(position) ; } sum /= (double)m ; sum2 /= (double)m ; stdev = sqrt(sum2 - sum*sum) ; median_dist = sinfo_new_median(pix_dist, m) ; for ( i = 0 ; i < (int) lx*ly ; i++ ) { /* blank pixels are not replaced */ if ( isnan(pdata[i]) ) { continue ; } /* initialize the buffer variables for the 8 nearest neighbors */ value = (pixelvalue * )cpl_calloc ( 8, sizeof ( pixelvalue * ) ) ; position = ( int * ) cpl_calloc ( 8, sizeof ( int * ) ) ; /*------------------------------------------------------------------- * determine the pixel position of the 8 nearest neighbors */ position[0] = i + lx - 1 ; /* upper left */ position[1] = i + lx ; /* upper */ position[2] = i + lx + 1 ; /* upper right */ position[3] = i + 1 ; /* right */ position[4] = i - lx + 1 ; /* lower right */ position[5] = i - lx ; /* lower */ position[6] = i - lx - 1 ; /* lower left */ position[7] = i - 1 ; /* left */ /*------------------------------------------------------------- * determine the positions of the image margins, top positions are * changed to low positions and vice versa. Right positions are * changed to left positions and vice versa. */ if ( i >= 0 && i < lx ) /* bottom line */ { position[4] += 2 * lx ; position[5] += 2 * lx ; position[6] += 2 * lx ; } else if ( i >= ((int) lx*ly - lx ) && i < (int) lx*ly ) /* top line */ { position[0] -= 2 * lx ; position[1] -= 2 * lx ; position[2] -= 2 * lx ; } else if ( i % lx == 0 ) /* left side */ { position[0] += 2 ; position[6] += 2 ; position[7] += 2 ; } else if ( i % lx == lx - 1 ) /* right side */ { position[2] -= 2 ; position[3] -= 2 ; position[4] -= 2 ; } /* ------------------------------------------------------------------- * read the pixel values of the neighboring pixels, * blanks are not considered */ nposition = 8 ; n = 0 ; for ( j = 0 ; j < nposition ; j ++ ) { if ( !isnan(pdata[position[j]]) ) { value[n] = pdata[position[j]] ; n ++ ; } } nposition = n ; if ( nposition <= 1 ) /* almost all neighbors are blank */ { pdata[i] = ZERO ; cpl_free(value) ; cpl_free(position) ; continue ; } /* determine the absolute distances */ dist = 0. ; for ( n = 0 ; n < nposition ; n++ ) { dist += (pdata[i] - value[n])*(pdata[i] - value[n]) ; } dist = sqrt(dist)/(float) nposition ; /* ----------------------------------------------------------------- * replace the pixel value by the sinfo_median on conditions: * fmedian = 0: always replace with sinfo_median. * fmedian < 0: interpret as absolute condition: * if |pixel - sinfo_median| > -fmedian * replace with sinfo_median. * fmedian > 0: replace by sinfo_median (fmedian as a factor of * the square root of the sinfo_median itself) * if |pixel - median| >= fmedian * sqrt ( median ) * considers a dependence on the pixel value. * This can be used to consider photon noise. */ if ( fmedian == 0 ) { pdata[i] = dist ; } else if ( fmedian < 0 && fabs ( median_dist - dist ) >= -fmedian*stdev ) { pdata[i] = dist ; } else if ( fmedian > 0 && fabs ( median_dist - dist ) >= fmedian*stdev * sqrt(fabs(dist)) ) { pdata[i] = dist ; } else { cpl_free (value) ; cpl_free (position) ; continue ; } cpl_free (value) ; cpl_free (position) ; } cpl_free(pix_dist); return image ; } /*---------------------------------------------------------------------------- Function: sinfo_new_local_median_image() In : im: input image fmedian: a factor to the local standard deviation loReject, hiReject: fraction of rejected values to determine a clean standard deviation half_box_size: integer half size of the running box to determine the local clean standard deviation Out : resulting image Job : filter, calculates the local stdev in a moving box Then it calculates the difference of the pixel to the median of the nearest neighbors by using the 8 closest pixels of every pixel. The values in the output image are determined according to the values of the input parameter. If fmedian = 0: always replace by median if fmedian < 0: replace median if |median_dist - dist| > fmedian * stdev if fmedian > 0: replace by median (fmedian as a factor of the square root of the median itself) if |pixel - median| >= fmedian*sqrt(median) This can be used to consider photon noise. This considers a dependence of the differences on the pixel values themselves. Notice : it is assumed that most of the 8 nearest neighbor pixels are not bad pixels! blank pixels are not replaced! ---------------------------------------------------------------------------*/ cpl_image * sinfo_new_local_median_image( cpl_image * im, float fmedian, float loReject, float hiReject, int half_box_size ) { cpl_image * image ; pixelvalue * value ; pixelvalue median ; int * position ; int nposition ; int n, i, j ; int llx, lly, urx, ury ; Stats * stats ; int lx=0; int ly=0; float* pidata=NULL; float* podata=NULL; if ( im == NULL ) { sinfo_msg_error ("no image input") ; return NULL ; } if ( half_box_size < 0 ) { sinfo_msg_error ("negativ box_size given") ; return NULL ; } image = cpl_image_duplicate ( im ) ; lx=cpl_image_get_size_x(im); ly=cpl_image_get_size_y(im); pidata=cpl_image_get_data(im); podata=cpl_image_get_data(image); /*---------------------------------------------------------------------- * go through all pixels */ for ( i = 0 ; i < (int) lx*ly ; i++ ) { /* blank pixels are not replaced */ if ( isnan(pidata[i]) ) { continue ; } /* compute the image statistics in the box area */ llx = i%lx - half_box_size ; if ( llx < 0 ) llx = 0 ; lly = i%ly - half_box_size ; if ( lly < 0 ) lly = 0 ; urx = i%lx + half_box_size ; if ( urx >= lx ) urx = lx - 1 ; ury = i%ly + half_box_size ; if ( ury >= ly ) ury = ly - 1 ; if ( NULL == (stats = sinfo_new_image_stats_on_rectangle (im, loReject, hiReject, llx, lly, urx, ury)) ) { sinfo_msg_warning("could not determine image statistics "); sinfo_msg_warning("in pixel %d", i) ; continue ; } /* initialize the buffer variables for the 8 nearest neighbors */ value = (pixelvalue * )cpl_calloc ( 8, sizeof ( pixelvalue * ) ) ; position = ( int * ) cpl_calloc ( 8, sizeof ( int * ) ) ; /*---------------------------------------------------------------------- * determine the pixel position of the 8 nearest neighbors */ position[0] = i + lx - 1 ; /* upper left */ position[1] = i + lx ; /* upper */ position[2] = i + lx + 1 ; /* upper right */ position[3] = i + 1 ; /* right */ position[4] = i - lx + 1 ; /* lower right */ position[5] = i - lx ; /* lower */ position[6] = i - lx - 1 ; /* lower left */ position[7] = i - 1 ; /* left */ /*--------------------------------------------------------------------- * determine the positions of the image margins, top positions are * changed to low positions and vice versa. Right positions are * changed to left positions and vice versa. */ if ( i >= 0 && i < lx ) /* bottom line */ { position[4] += 2 * lx ; position[5] += 2 * lx ; position[6] += 2 * lx ; } else if ( i >= ((int) lx*ly - lx ) && i < (int) lx*ly ) /* top line */ { position[0] -= 2 * lx ; position[1] -= 2 * lx ; position[2] -= 2 * lx ; } else if ( i % lx == 0 ) /* left side */ { position[0] += 2 ; position[6] += 2 ; position[7] += 2 ; } else if ( i % lx == lx - 1 ) /* right side */ { position[2] -= 2 ; position[3] -= 2 ; position[4] -= 2 ; } /* ------------------------------------------------------------------ * read the pixel values of the neighboring pixels, * blanks are not considered */ nposition = 8 ; n = 0 ; for ( j = 0 ; j < nposition ; j ++ ) { if ( !isnan(pidata[position[j]]) ) { value[n] = pidata[position[j]] ; n ++ ; } } nposition = n ; if ( nposition <= 1 ) /* almost all neighbors are blank */ { podata[i] = ZERO ; cpl_free(value) ; cpl_free(position) ; cpl_free(stats) ; continue ; } /* sort the values and determine the median */ sinfo_pixel_qsort( value, nposition ) ; if ( nposition % 2 == 1 ) { median = value [ nposition/2 ] ; } else { median = ( value [nposition/2 - 1] + value [nposition/2] ) / 2. ; } /* ----------------------------------------------------------------- * replace the pixel value by the median on conditions: * fmedian = 0: always replace with median. * fmedian > 0: replace by median (fmedian as a factor of * the square root of the median itself) * if |pixel - median| >= fmedian * sqrt ( median ) * considers a dependence on the pixel value. * This can be used to consider photon noise. */ if ( fmedian == 0 ) { podata[i] = median ; } else if ( fmedian < 0 && fabs ( median - pidata[i] ) >= -fmedian * stats->cleanstdev) { podata[i] = median ; } else if ( fmedian > 0 && fabs ( median - pidata[i] ) >= fmedian * sqrt(fabs(median)) ) { podata[i] = median ; } else { cpl_free (value) ; cpl_free (position) ; cpl_free (stats) ; continue ; } cpl_free (value) ; cpl_free (position) ; cpl_free (stats) ; } return image ; } /*---------------------------------------------------------------------------- Function: sinfo_new_mean_image_in_spec() In : image, a threshold parameter Out : resulting image Job : mean filter, calculates the mean for an image by using the 4 closest pixels of every pixel in spectral direction (column). The values in the output image are determined according to the values of the input parameter. If fmedian = 0: always replace by mean if fmedian < 0: replace by mean if |pixel - mean| > -fmedian if fmedian > 0: replace by mean (fmedian as a factor of the square root of the mean itself) if |pixel - mean| >= fmedian * sqrt ( mean ) This can be used to consider photon noise. This considers a dependence of the differences on the pixel values themselves. Notice : it is assumed that most of the 4 nearest neighbor pixels are not bad pixels! blank pixels are not replaced! ---------------------------------------------------------------------------*/ cpl_image * sinfo_new_mean_image_in_spec( cpl_image * im, float fmedian ) { cpl_image * image ; pixelvalue * value ; pixelvalue mean ; int * position ; int nposition ; int n, i, j ; int lx=0; int ly=0; float* pidata=NULL; float* podata=NULL; if ( im == NULL ) { sinfo_msg_error ("no image input") ; return NULL ; } image = cpl_image_duplicate ( im ) ; lx=cpl_image_get_size_x(im); ly=cpl_image_get_size_y(im); pidata=cpl_image_get_data(im); podata=cpl_image_get_data(image); /*---------------------------------------------------------------------- * go through all pixels */ for ( i = 0 ; i < (int) lx*ly ; i++ ) { /* blank pixels are not replaced */ if ( isnan(pidata[i]) ) { continue ; } /* initialize the buffer variables for the 2 nearest spectral neighbors */ value = (pixelvalue * )cpl_calloc ( 4, sizeof ( pixelvalue * ) ) ; position = ( int * ) cpl_calloc ( 4, sizeof ( int * ) ) ; /*-------------------------------------------------------------------- * determine the pixel position of the 8 nearest neighbors */ position[0] = i + lx ; /* upper */ position[1] = i + 2*lx ; /* upper */ position[2] = i - lx ; /* lower */ position[3] = i - 2*lx ; /* lower */ /*------------------------------------------------------------------- * determine the positions of the image margins, top positions are * changed to low positions and vice versa. Right positions are changed * to left positions and vice versa. */ if ( i >= 0 && i < lx ) /* bottom line */ { position[2] += 2 * lx ; position[3] += 4 * lx ; } else if ( i >= ((int) lx*ly - lx ) && i < (int) lx*ly ) /* top line */ { position[0] -= 2 * lx ; position[1] -= 4 * lx ; } /* ------------------------------------------------------------------- * read the pixel values of the neighboring pixels, * blanks are not considered */ nposition = 4 ; n = 0 ; for ( j = 0 ; j < nposition ; j ++ ) { if ( !isnan(pidata[position[j]]) ) { value[n] = pidata[position[j]] ; n ++ ; } } nposition = n ; if ( nposition < 1 ) /* all neighbors are blank */ { podata[i] = ZERO ; cpl_free(value) ; cpl_free(position) ; continue ; } /* determine the mean */ mean = 0. ; for ( n = 0 ; n < nposition ; n++ ) { mean += value[n] ; } mean /= (float) nposition ; /* ----------------------------------------------------------------- * replace the pixel value by the median on conditions: * fmedian = 0:","always replace with mean. * fmedian < 0: interpret as absolute condition: * if |pixel - mean| > -fmedian * replace with mean. */ if ( fmedian == 0 ) { podata[i] = mean ; } else if ( fmedian < 0 && fabs ( mean - pidata[i] ) >= -fmedian ) { podata[i] = mean ; } else if ( fmedian > 0 && fabs ( mean - pidata[i] ) >= fmedian * sqrt(fabs(mean)) ) { podata[i] = mean ; } else { cpl_free (value) ; cpl_free (position) ; continue ; } cpl_free (value) ; cpl_free (position) ; } return image ; } /**@}*/