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
where xx,yy,d, and w have the same length. Also, xx and yy are matrix x,y locations.
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
EDIT:
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
(a function implementing this)
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])
