# An efficient 'drizzle' algorithm?

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


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])


This algorithm is easy to write in a for-loop fashion and hard to express as a series of vector operations, which makes is a perfect candidate for writing it in a lower-level language (like C) and just linking it to Python (for instance with Cython). Unless there is a drizzle lib for Python, this will be both fastest and easiest solution.

• Right, that's exactly the point - for loop is inefficient, so is there a C library that does this? numpy.bincount is pretty good, but maybe there's better? Dec 10, 2011 at 0:00
• But this is 5 lines of C plus few for interface... If you are looking for a pure Python, I would rather seek for some image manipulation-based solution.
– mbq
Dec 10, 2011 at 0:25