What efficient implementations of a 'drizzle' algorithm are available? The problem is, given a timestream of data in which each element is associated with a pixel in a map, how do you create that map? Each pixel may have many data points associated with it. Each data point may need to be weighted.
For example, in python/numpy, given a data array
d, a weight array
w, a map
m, a weight map
wm, and a mapping from d to m
xinds,yinds, you could do:
for jj,(xx,yy) in enumerate(xinds,yinds): m[xx,yy] += (d*w)[jj] wm[xx,yy] += w[jj] final_image = m/wm
w have the same length. Also,
yy are matrix
How can this be made more efficient? Are there tools in python or libraries in other languages to do this? Am I even calling the algorithm by its right name?
An efficient implementation in
IDL using the
histogram function is shown at David Fanning's website
After asking this question, I realized I had the answer...
numpy.bincount does exactly what I want in numpy. If the mapping
t = xinds + yinds*xsize where
xsize is the x-dimension of the map
# shapes of x,y indices need to be flat x,y = (a.ravel() for a in numpy.indices(m.shape)) dc = numpy.bincount(t,d*w) wc = numpy.bincount(t,w) m[x,y] = dc wm[x,y] = wc
It would still be useful to know of other implementations of this algorithm. Or perhaps ways to compute the
t mapping and different weighting schemes - I don't discuss that at all above, but I think in the context in which the term was coined (Hubble imaging) there are complexities involved in determining both variables.
EDIT2: Corollary - what if I want to median-drizzle? i.e., instead of averaging for the final map, median?
(this is not valid code, but vaguely pseudo-code... you can't have 2-dimensional lists in python, though nested lists are OK)
for jj,(xx,yy) in enumerate(xinds,yinds): m[xx,yy].append((d*w)[jj]) for jj,(xx,yy) in enumerate(xinds,yinds): m[xx,yy] = median(m[xx,yy])