LSSTDESC / descwl_coadd

Tools for coadding using the LSST DM software stack
GNU General Public License v3.0
2 stars 2 forks source link

Warp and coadd "on-the-fly" #44

Closed esheldon closed 3 years ago

esheldon commented 3 years ago

I learned from @erykoff that we can do warping and adding of images to the coadd on the fly.

This doesn't matter for the sim code where each image is realized when it is passed in. But the stack passes around objects that have not yet staged their data, and this can save huge amounts of memory.

I don't remember how to do it though.

erykoff commented 3 years ago

Right now what the official Stack does is warp ccds into the coadd frame and persist them (and it does one warp per visit rather than one per ccd, but this is an implementation detail). The coaddition is performed by by the afwMath.statisticsStack method which can perform all sorts of statistics (mean, weighted mean, median, sigma-clipping, etc). Of course in order to do median and sigma-clipped means at every position you need to know the value for every input to the coadd stack. To avoid the memory blowing up too much the stack will read in sub-regions of each input warp and do the statistics on the tower in chunks.

However, if we don't want medians or sigma clipping and simply want a mean or weighted mean, then we don't need to (a) use afwMath.statisticsStack or (b) have everything in memory. You simply need to take each ccd, warp it, and add the warped values (or weight*value) to the pixel accumulator and add 1 (or weight) to the denominator accumulator (or whatever it should be called). Then you can discard the warp and go onto the next input. In the end you divide sum(weight*value)/sum(weight) and you're done.

Originally when we discussed this I thought it would be easy because all of the inputs to the cell-based coadds have no edges, but in fact it's easy because we want a simple weighted mean, and of course can be used to save memory in any (weighted) mean stacking. (And in fact is essentially what I do in https://github.com/LSSTDESC/supreme/ ).

esheldon commented 3 years ago

thanks Eli. Where should I look there?

erykoff commented 3 years ago

Oh, in suprême, it's not directly applicable, but the main loop is https://github.com/LSSTDESC/supreme/blob/master/supreme/patch_mapper.py#L245-L246 and it just accumulates map pixels from each input image. The point being is that to compute mean, weighted mean, min, max, and sum you don't need to have the full stack in memory.