LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
accumulator_mean_stack.py
Go to the documentation of this file.
1# This file is part of meas_algorithms.
2#
3# LSST Data Management System
4# This product includes software developed by the
5# LSST Project (http://www.lsst.org/).
6# See COPYRIGHT file at the top of the source tree.
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <https://www.lsstcorp.org/LegalNotices/>.
21#
22import warnings
23
24import numpy as np
25
26
27__all__ = ['AccumulatorMeanStack']
28
29
31 """Stack masked images.
32
33 Parameters
34 ----------
35 shape : `tuple`
36 Shape of the input and output images.
37 bit_mask_value : `int`
38 Bit mask to flag for "bad" inputs that should not be stacked.
39 mask_threshold_dict : `dict` [`int`: `float`], optional
40 Dictionary of mapping from bit number to threshold for flagging.
41 Only bad bits (in bit_mask_value) which mask fractional weight
42 greater than this threshold will be flagged in the output image.
43 mask_map : `list` [`tuple`], optional
44 Mapping from input image bits to aggregated coadd bits.
45 no_good_pixels_mask : `int`, optional
46 Bit mask to set when there are no good pixels in the stack.
47 If not set then will set coadd masked image 'NO_DATA' bit.
48 calc_error_from_input_variance : `bool`, optional
49 Calculate the error from the input variance?
50 compute_n_image : `bool`, optional
51 Calculate the n_image map as well as stack?
52 """
53 def __init__(self, shape,
54 bit_mask_value, mask_threshold_dict={},
55 mask_map=[], no_good_pixels_mask=None,
56 calc_error_from_input_variance=True,
57 compute_n_image=False):
58 self.shape = shape
59 self.bit_mask_value = bit_mask_value
60 self.mask_map = mask_map
61 self.no_good_pixels_mask = no_good_pixels_mask
62 self.calc_error_from_input_variance = calc_error_from_input_variance
63 self.compute_n_image = compute_n_image
64
65 # Only track threshold bits that are in the bad bit_mask_value.
67 for bit in mask_threshold_dict:
68 if (self.bit_mask_value & 2**bit) > 0:
69 self.mask_threshold_dict[bit] = mask_threshold_dict[bit]
70
71 # sum_weight holds the sum of weights for each pixel.
72 self.sum_weight = np.zeros(shape, dtype=np.float64)
73 # sum_wdata holds the sum of weight*data for each pixel.
74 self.sum_wdata = np.zeros(shape, dtype=np.float64)
75
76 if calc_error_from_input_variance:
77 # sum_w2var holds the sum of weight**2 * variance for each pixel.
78 self.sum_w2var = np.zeros(shape, dtype=np.float64)
79 else:
80 # sum_weight2 holds the sum of weight**2 for each pixel.
81 self.sum_weight2 = np.zeros(shape, dtype=np.float64)
82 # sum_wdata2 holds the sum of weight * data**2 for each pixel.
83 self.sum_wdata2 = np.zeros(shape, dtype=np.float64)
84
85 self.or_mask = np.zeros(shape, dtype=np.int64)
87 for bit in self.mask_threshold_dict:
88 self.rejected_weights_by_bit[bit] = np.zeros(shape, dtype=np.float64)
89
90 self.masked_pixels_mask = np.zeros(shape, dtype=np.int64)
91
92 if self.compute_n_image:
93 self.n_image = np.zeros(shape, dtype=np.int32)
94
95 def add_masked_image(self, masked_image, weight=1.0):
96 """Add a masked image to the stack.
97
98 Parameters
99 ----------
100 masked_image : `lsst.afw.image.MaskedImage`
101 Masked image to add to the stack.
102 weight : `float` or `np.ndarray`, optional
103 Weight to apply for weighted mean. If an array,
104 must be same size and shape as input masked_image.
105 """
106 good_pixels = np.where(((masked_image.mask.array & self.bit_mask_value) == 0)
107 & np.isfinite(masked_image.mask.array))
108
109 self.sum_weight[good_pixels] += weight
110 self.sum_wdata[good_pixels] += weight*masked_image.image.array[good_pixels]
111
112 if self.compute_n_image:
113 self.n_image[good_pixels] += 1
114
116 self.sum_w2var[good_pixels] += (weight**2.)*masked_image.variance.array[good_pixels]
117 else:
118 self.sum_weight2[good_pixels] += weight**2.
119 self.sum_wdata2[good_pixels] += weight*(masked_image.image.array[good_pixels]**2.)
120
121 # Mask bits are propagated for good pixels
122 self.or_mask[good_pixels] |= masked_image.mask.array[good_pixels]
123
124 # Bad pixels are only tracked if they cross a threshold
125 for bit in self.mask_threshold_dict:
126 bad_pixels = ((masked_image.mask.array & 2**bit) > 0)
127 self.rejected_weights_by_bit[bit][bad_pixels] += weight
128 self.masked_pixels_mask[bad_pixels] |= 2**bit
129
130 def fill_stacked_masked_image(self, stacked_masked_image):
131 """Fill the stacked mask image after accumulation.
132
133 Parameters
134 ----------
135 stacked_masked_image : `lsst.afw.image.MaskedImage`
136 Total masked image.
137 """
138 with warnings.catch_warnings():
139 # Let the NaNs through and flag bad pixels below
140 warnings.simplefilter("ignore")
141
142 # The image plane is sum(weight*data)/sum(weight)
143 stacked_masked_image.image.array[:, :] = self.sum_wdata/self.sum_weight
144
146 mean_var = self.sum_w2var/(self.sum_weight**2.)
147 else:
148 # Compute the biased estimator
149 variance = self.sum_wdata2/self.sum_weight - stacked_masked_image.image.array[:, :]**2.
150 # De-bias
151 variance *= (self.sum_weight**2.)/(self.sum_weight**2. - self.sum_weight2)
152
153 # Compute the mean variance
154 mean_var = variance*self.sum_weight2/(self.sum_weight**2.)
155
156 stacked_masked_image.variance.array[:, :] = mean_var
157
158 # Propagate bits when they cross the threshold
159 for bit in self.mask_threshold_dict:
160 hypothetical_total_weight = self.sum_weight + self.rejected_weights_by_bit[bit]
161 self.rejected_weights_by_bit[bit] /= hypothetical_total_weight
162 propagate = np.where(self.rejected_weights_by_bit[bit] > self.mask_threshold_dict[bit])
163 self.or_mask[propagate] |= 2**bit
164
165 # Map mask planes to new bits for pixels that had at least one
166 # bad input rejected and are in the mask_map.
167 for mask_tuple in self.mask_map:
168 self.or_mask[(self.masked_pixels_mask & mask_tuple[0]) > 0] |= mask_tuple[1]
169
170 stacked_masked_image.mask.array[:, :] = self.or_mask
171
172 if self.no_good_pixels_mask is None:
173 mask_dict = stacked_masked_image.mask.getMaskPlaneDict()
174 no_good_pixels_mask = 2**(mask_dict['NO_DATA'])
175 else:
176 no_good_pixels_mask = self.no_good_pixels_mask
177
178 bad_pixels = (self.sum_weight <= 0.0)
179 stacked_masked_image.mask.array[bad_pixels] |= no_good_pixels_mask
180
181 def add_image(self, image, weight=1.0):
182 """Add an image to the stack.
183
184 No bit-filtering is performed when adding an image.
185
186 Parameters
187 ----------
188 image : `lsst.afw.image.Image`
189 Image to add to the stack.
190 weight : `float` or `np.ndarray`, optional
191 Weight to apply for weighted mean. If an array,
192 must be same size and shape as input image.
193 """
194 self.sum_weight[:, :] += weight
195 self.sum_wdata[:, :] += weight*image.array[:]
196
197 if self.compute_n_image:
198 self.n_image[:, :] += 1
199
200 def fill_stacked_image(self, stacked_image):
201 """Fill the image after accumulation.
202
203 Parameters
204 ----------
205 stacked_image : `lsst.afw.image.Image`
206 Total image.
207 """
208 with warnings.catch_warnings():
209 # Let the NaNs through, this should only happen
210 # if we're stacking with no inputs.
211 warnings.simplefilter("ignore")
212
213 # The image plane is sum(weight*data)/sum(weight)
214 stacked_image.array[:, :] = self.sum_wdata/self.sum_weight
215
216 @staticmethod
218 """Convert stats control to threshold dict.
219
220 Parameters
221 ----------
222 stats_ctrl : `lsst.afw.math.StatisticsControl`
223
224 Returns
225 -------
226 threshold_dict : `dict`
227 Dict mapping from bit to propagation threshold.
228 """
229 threshold_dict = {}
230 for bit in range(64):
231 threshold_dict[bit] = stats_ctrl.getMaskPropagationThreshold(bit)
232
233 return threshold_dict
__init__(self, shape, bit_mask_value, mask_threshold_dict={}, mask_map=[], no_good_pixels_mask=None, calc_error_from_input_variance=True, compute_n_image=False)