191 def run(self, mapperResults, exposure, **kwargs):
192 """Reduce a list of items produced by `ImageMapper`.
193
194 Either stitch the passed `mapperResults` list
195 together into a new Exposure (default) or pass it through
196 (if `self.config.reduceOperation` is 'none').
197
198 If `self.config.reduceOperation` is not 'none', then expect
199 that the `pipeBase.Struct`s in the `mapperResults` list
200 contain sub-exposures named 'subExposure', to be stitched back
201 into a single Exposure with the same dimensions, PSF, and mask
202 as the input `exposure`. Otherwise, the `mapperResults` list
203 is simply returned directly.
204
205 Parameters
206 ----------
207 mapperResults : `list`
208 list of `lsst.pipe.base.Struct` returned by `ImageMapper.run`.
209 exposure : `lsst.afw.image.Exposure`
210 the original exposure which is cloned to use as the
211 basis for the resulting exposure (if
212 ``self.config.mapper.reduceOperation`` is not 'None')
213 kwargs :
214 additional keyword arguments propagated from
215 `ImageMapReduceTask.run`.
216
217 Returns
218 -------
219 A `lsst.pipe.base.Struct` containing either an `lsst.afw.image.Exposure`
220 (named 'exposure') or a list (named 'result'),
221 depending on `config.reduceOperation`.
222
223 Notes
224 -----
225 1. This currently correctly handles overlapping sub-exposures.
226 For overlapping sub-exposures, use `config.reduceOperation='average'`.
227 2. This correctly handles varying PSFs, constructing the resulting
228 exposure's PSF via CoaddPsf (DM-9629).
229
230 Known issues
231
232 1. To be done: correct handling of masks (nearly there)
233 2. This logic currently makes *two* copies of the original exposure
234 (one here and one in `mapper.run()`). Possibly of concern
235 for large images on memory-constrained systems.
236 """
237
238 if self.config.reduceOperation == 'none':
239 return pipeBase.Struct(result=mapperResults)
240
241 if self.config.reduceOperation == 'coaddPsf':
242
243 coaddPsf = self._constructPsf(mapperResults, exposure)
244 return pipeBase.Struct(result=coaddPsf)
245
246 newExp = exposure.clone()
247 newMI = newExp.getMaskedImage()
248
249 reduceOp = self.config.reduceOperation
250 if reduceOp == 'copy':
251 weights = None
252 newMI.image.array[:, :] = np.nan
253 newMI.variance.array[:, :] = np.nan
254 else:
255 newMI.image.array[:, :] = 0.
256 newMI.variance.array[:, :] = 0.
257 if reduceOp == 'average':
258 weights = afwImage.ImageI(newMI.getBBox())
259
260 for item in mapperResults:
261 item = item.subExposure
262 if not (isinstance(item, afwImage.ExposureF) or isinstance(item, afwImage.ExposureI)
263 or isinstance(item, afwImage.ExposureU) or isinstance(item, afwImage.ExposureD)):
264 raise TypeError("""Expecting an Exposure type, got %s.
265 Consider using `reduceOperation="none".""" % str(type(item)))
266 subExp = newExp.Factory(newExp, item.getBBox())
267 subMI = subExp.getMaskedImage()
268 patchMI = item.getMaskedImage()
269 isValid = ~np.isnan(patchMI.image.array * patchMI.variance.array)
270
271 if reduceOp == 'copy':
272 subMI.image.array[isValid] = patchMI.image.array[isValid]
273 subMI.variance.array[isValid] = patchMI.variance.array[isValid]
274 subMI.mask.array[:, :] |= patchMI.mask.array
275
276 if reduceOp == 'sum' or reduceOp == 'average':
277 subMI.image.array[isValid] += patchMI.image.array[isValid]
278 subMI.variance.array[isValid] += patchMI.variance.array[isValid]
279 subMI.mask.array[:, :] |= patchMI.mask.array
280 if reduceOp == 'average':
281
282 wtsView = afwImage.ImageI(weights, item.getBBox())
283 wtsView.array[isValid] += 1
284
285
286 mask = newMI.mask
287 for m in self.config.badMaskPlanes:
288 mask.addMaskPlane(m)
289 bad = mask.getPlaneBitMask(self.config.badMaskPlanes)
290
291 isNan = np.where(np.isnan(newMI.image.array * newMI.variance.array))
292 if len(isNan[0]) > 0:
293
294 mask.array[isNan[0], isNan[1]] |= bad
295
296 if reduceOp == 'average':
297 wts = weights.array.astype(float)
298 self.log.info('AVERAGE: Maximum overlap: %f', np.nanmax(wts))
299 self.log.info('AVERAGE: Average overlap: %f', np.nanmean(wts))
300 self.log.info('AVERAGE: Minimum overlap: %f', np.nanmin(wts))
301 wtsZero = np.equal(wts, 0.)
302 wtsZeroInds = np.where(wtsZero)
303 wtsZeroSum = len(wtsZeroInds[0])
304 self.log.info('AVERAGE: Number of zero pixels: %f (%f%%)', wtsZeroSum,
305 wtsZeroSum * 100. / wtsZero.size)
306 notWtsZero = ~wtsZero
307 tmp = newMI.image.array
308 np.divide(tmp, wts, out=tmp, where=notWtsZero)
309 tmp = newMI.variance.array
310 np.divide(tmp, wts, out=tmp, where=notWtsZero)
311 if len(wtsZeroInds[0]) > 0:
312 newMI.image.array[wtsZeroInds] = np.nan
313 newMI.variance.array[wtsZeroInds] = np.nan
314
315
316 mask.array[wtsZeroInds] |= bad
317
318
319 if reduceOp == 'sum' or reduceOp == 'average':
320 psf = self._constructPsf(mapperResults, exposure)
321 newExp.setPsf(psf)
322
323 return pipeBase.Struct(exposure=newExp)
324