Coverage for src/ufig/plugins/background_subtract.py: 85%

1# Copyright (c) 2016 ETH Zurich, Institute of Astronomy, Tomasz Kacprzak 

2# <tomasz.kacprzak@phys.ethz.ch> 

3""" 

4Created on Jan 5, 2016 

5@author: Tomasz Kacprzak 

6 

7""" 

8 

9import numpy as np 

10from ivy.plugin.base_plugin import BasePlugin 

11 

12__all__ = ["sigma_clip"] 

13 

14 

15class Plugin(BasePlugin): 

16 """ 

17 Subtracting background from an image, following recipe in SWarp 

18 https://www.astromatic.net/pubsvn/software/swarp/trunk/doc/swarp.pdf 

19 

20 :param bgsubract_n_iter: number if iterations for sigma clipping, higher 

21 bgsubract_n_iter means more precision and less speed 

22 :param bgsubract_n_downsample: use every n_downsample pixel, higher n_downsample 

23 means faster and less accurate background estimation 

24 

25 :return: Image with background subtracted 

26 """ 

27 

28 def __call__(self): 

29 if hasattr(self.ctx, "image_mask"): 29 ↛ 30line 29 didn't jump to line 30 because the condition on line 29 was never true

30 img_flat = self.ctx.image[self.ctx.image_mask].ravel() 

31 else: 

32 img_flat = self.ctx.image.ravel() 

33 data = img_flat[:: self.ctx.parameters.bgsubract_n_downsample] 

34 data_clipped = sigma_clip( 

35 data, sig=3, iters=self.ctx.parameters.bgsubract_n_iter, maout="inplace" 

36 ) 

37 background_level = 2.5 * np.ma.median(data_clipped) - 1.5 * np.ma.mean( 

38 data_clipped 

39 ) 

40 self.ctx.image -= background_level 

41 

42 if hasattr(self.ctx, "image_mask"): 42 ↛ 43line 42 didn't jump to line 43 because the condition on line 42 was never true

43 self.ctx.image[~self.ctx.image_mask] = 0 

44 

45 def __str__(self): 

46 return "subtract background" 

47 

48 

49# modified sigma clip function 

50def sigma_clip(data, sig=3, iters=1, cenfunc=np.median, varfunc=np.var, maout=False): 

51 """ 

52 Perform sigma-clipping on the provided data. 

53 

54 This performs the sigma clipping algorithm - i.e. the data will be iterated 

55 over, each time rejecting points that are more than a specified number of 

56 standard deviations discrepant. 

57 

58 .. note:: 

59 `scipy.stats.sigmaclip` provides a subset of the functionality in this 

60 function. 

61 

62 :param data: array-like The data to be sigma-clipped (any shape). 

63 :param sig: float The number of standard deviations (*not* variances) to use as the 

64 clipping limit. 

65 :param iters: int or NoneThe number of iterations to perform clipping for, or None 

66 to clip until convergence is achieved (i.e. continue until the last 

67 iteration clips nothing). 

68 :param cenfunc: callable The technique to compute the center for the clipping. Must 

69 be a callable that takes in a 1D data array and outputs the central 

70 value. Defaults to the median. 

71 :param varfunc: callable The technique to compute the variance about the center. 

72 Must be a callable that takes in a 1D data array and outputs the 

73 width estimator that will be interpreted as a variance. Defaults to 

74 the variance. 

75 :param maout: bool or 'copy' If True, a masked array will be returned. If the 

76 special string 'inplace', the masked array will contain the same array 

77 as `data`, otherwise the array data will be copied. 

78 

79 :return: filtereddata : `numpy.ndarray` or `numpy.masked.MaskedArray` 

80 If `maout` is True, this is a masked array with a shape matching the 

81 input that is masked where the algorithm has rejected those values. 

82 Otherwise, a 1D array of values including only those that are not 

83 clipped. 

84 :return: mask : boolean array 

85 Only present if `maout` is False. A boolean array with a shape matching 

86 the input `data` that is False for rejected values and True for all 

87 others. 

88 

89 Examples 

90 

91 This will generate random variates from a Gaussian distribution and return 

92 only the points that are within 2 *sample* standard deviation from the 

93 median:: 

94 

95 >>> from astropy.stats import sigma_clip 

96 >>> from numpy.random import randn 

97 >>> randvar = randn(10000) 

98 >>> data,mask = sigma_clip(randvar, 2, 1) 

99 

100 This will clipping on a similar distribution, but for 3 sigma relative to 

101 the sample *mean*, will clip until converged, and produces a 

102 `numpy.masked.MaskedArray`:: 

103 

104 >>> from astropy.stats import sigma_clip 

105 >>> from numpy.random import randn 

106 >>> from numpy import mean 

107 >>> randvar = randn(10000) 

108 >>> maskedarr = sigma_clip(randvar, 3, None, mean, maout=True) 

109 

110 """ 

111 

112 data = np.array(data, copy=False) 

113 oldshape = data.shape 

114 data = data.ravel() 

115 data_centered = np.zeros_like(data) 

116 

117 mask = np.ones(data.size, bool) 

118 

119 for _ in range(iters): 

120 std_nsig = np.sqrt(varfunc(data[mask]) * sig**2) 

121 mean = cenfunc(data[mask]) 

122 data_centered[:] = data[:] 

123 data_centered -= mean 

124 mask = data_centered <= std_nsig 

125 

126 if maout: 126 ↛ 129line 126 didn't jump to line 129 because the condition on line 126 was always true

127 return np.ma.MaskedArray(data, ~mask, copy=maout != "inplace") 

128 else: 

129 return data[mask], mask.reshape(oldshape)