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 /*--------------------------------------------------------------------------*/
1.2.13.1 written by Dimitri van Heesch,
© 1997-2001