# Presentation: 016074

## A fast flux-conserving resampling algorithm

**Procedure**
- Figure 1 illustrates the typical resampling problem.  A grid of resampling pixels (“rixels”) is superposed over some original grid of image pixels.  For a flux-conserving resampling, we want to integrate each rixel over the area it covers.   If we assume that the flux within each pixel is constant, the resampling problem reduces to finding the set of pixels overlapped by each rixel and then for a given rixel-pixel pair calculating the overlap area. The total flux in the rixel is then just the sum of contributions from the overlapped pixels.
- Calculating the set of candidate overlapping pixels for a given rixel is illustrated in Figure 2.  For resampling of comparably sized grids, about 1-10 candidate pixels will be found for each rixel.
- Each candidate pixel is treated as a clipping box for the rixel.  A Sutherland-Hodgeman clipping (1974) is applied as illustrated in Figures 3 and 4. The rixel is originally a trapezoid described as a list of  vertices.  As it is clipped by each of the bounds of the pixel, vertices can be added or deleted from the vertex list.
- One key simplifying assumption is that the rixel is convex.  Clipping a convex area with a convex clipping region (the pixel) must return another convex area (or nothing).  If the rixel were non-convex, it could be split into several distinct polygons.  The rixel is clipped in turn by all four boundaries of the pixel.

**Limits and issues**
*Are rixel boundaries really straight?*
- One of the most commonly use projections is the  gnomonic  or tangent projection.  This is the natural projection for data taken in a flat focal plane.  In this projection straight lines correspond to great circles on the sky, so that a straight line in one tangent projection is straight in any other.  In this common case where both the rixels and pixels are defined using a tangent projection there is no approximation.
- When other projections, or specific telescope distortions are involved, these typically have curvatures with scales of a radian.  The actual curvature of the rixel boundaries means that the algorithm assigns areas of pixels to the incorrect rixels.  The fraction of the area of the rixel that is incorrectly assigned will normally be something like the square of the size of the pixel in radians.  For 1” pixels this is <10-10.  This error does not affect the conservation of flux, it incorrectly assigns a tiny fraction of the flux to the wrong rixels.
*How do we handle intensive quantities like temperature?*
- Intensive quantities can be computed in the same fashion, but after summing the contributions from all pixels, the rixel should be divided by its area to give a weighted average over the rixel.
*Convexity of the rixel trapezoid seems to be important.  Can rixels be concave?*
- Yes, but  generally only near a singularity in one of the projections.  Looking at the rixel to the right, the sides A and B are perpendicular in the projection in which the rixels form a rectangular grid.  For them to  concave there must be a local distortion of angles (not just a rotation) of more than 90 degrees between the rixel and pixel frames.  The basic algorithm described here can be adapted to concave rixels, but the current implementation does not do so.
*What if the flux is not uniformly distributed in the pixel?*
- If the integral of the flux within an arbitrary triangle is easily computed, then this method could be adapted to handle non-uniform flux.  However if the function being resampled is known analytically, then other methods of interpolation are likely to be faster and more effective.
*Can the algorithm handling a ‘drizzle’ resampling?  Smoothing?*
- This is easily accommodated so long as the drizzle region for the pixel is rectangular with the same orientation as the original pixel.  The clip region for each pixel is set to the drizzle region, rather than the pixel boundaries.  A kind of box-car averaging might be implemented by setting the clipping region  larger than the pixel.  This would naturally allow non-integral values for the box-car size.
*Where can one get a copy of the algorithm?*
- A simple Java library for the resampling is available at
- http://skyview.gsfc.nasa.gov/polysamp
- This includes only the resampling itself.  Other libraries are available for pixel reprojection.

- Thomas A. McGlynn
*NASA/GSFC*

**Figure 1**.  A blue grid of resampling pixels (‘rixels’) overlays the black grid of the original image.

**Figure 2**.  For a given single rixel we calculate the bounding-box.
- The shaded pixels, which overlap the bounding box, are candidates for contributions to the rixel.

- If there is any overlap, the result of clipping by all four edges of the pixel is just the vertex list of the overlap of the rixel and pixel. It is easy to decompose the overlap into a set of triangles  (Figure 5).   Standard geometric formulae are used to compute the area of the triangles and hence the contribution of this pixel to the rixel flux.
- The algorithm iterates over the over all candidate pixels for a given rixel, and over each rixel in the resampled image.  Since the algorithm is local, can easily be parallelized to run on multiple processors.

**Figure 3**.  Given an input vertex list, the table indicates the 0-2 vertices to be output for each input depending upon whether it and the previous vertex satisfy the constraint.  Crossover means to output the point at which the current edge crosses over the constraint boundary.

- The image above of the region around 3c273 has been flipped, rescaled and rotated on the left using clipping resampling.  The integral over the image is conserved to one part in 1013.  A bilinear interpolation conserved flux to about 0.4%.

**Examples**
- A set of lines and points were resampled using nearest neighbor resampling and the clipping algorithm after a rescaling and rotation.  Features are introduced in the NN resampling.  Bilinear interpolation produces images similar in appearance to the clipping resampling, but flux is typically conserved only to about 1%.  Clipping conserves flux limited only by roundoff errors.

**References**
- Sutherland, I.E., and Hodgeman, G.W., 1974, **CACM** *17*, 32.
- Liang, Y.-D.,  and Barsky, B.A., 1984, **CACM** *27*, 151.
- Fruchter, A.S.,  and Hook, R.N., 2002, **PASP** *114*, 144.

**Figure 4.**  The rixel is clipped on each side in turn.  A more efficient implementation might notice the third clipping is not needed.

**Figure 5**. Calculating the overlap area is done by decomposing the overlap into triangles.

- The overlap is a convex region and the vertices are listed in order.   We can be split the overlap into triangles using vertices 1,2,3; 1,3,4; 1,4,5 (and 1,5,6, ... if there were more vertices).

**Efficiency**
- While the computational cost of this algorithm is larger than for the much simpler nearest neighbor or bilinear interpolation, the CPU cost is not vastly greater than the costs for projecting the rixels into the pixel frame.  Since reprojection is required for all  interpolation methods, clipping resampling ends up being 5-10 times slower than the other methods overall.  No transcendental functions are required for the resampling, though the calculation of the area of each triangle does require a square root.
- Resampling of of a 500x500 pixel  3C 273 image into an 850x850 grid (on the right) required about 3.0 seconds on a 1 GHz Linux machine in a pure Java implementation.  This compares with about 0.5-0.7 seconds for nearest neighbor and bilinear interpolations.   For comparable sizes of the rixels and pixels, this simple implementation resamples at about 100,000 rixels per second per Gigahertz of processor.  Thus even for megapixel regions, this algorithm can give reasonable interactive response.  The code is easily parallelized, but less able to take advantage of vector processors.
- There are several approaches that may substantially speed up the algorithm.  Direct coding in C may address language inefficiencies in Java.  Simple modifications to the code may make modest improvements (e.g., see figure 4).  More sophisticated clipping techniques, e.g., the Liang-Barsky (1984) may require fewer comparisons.  While I am unaware of this algorithm being used elsewhere in astronomy it seems not unlikely that similar or better techniques have been developed in the image processing arena.

- Original

- Nearest neighbor

- Clipping

- 1

- 2

- 3

- 4

- 5