Interpolation flux conservation and IMEXAM measurement

Frank Valdes wrote on Jun 25, 1999

U = User, LD = Lindsey Davis

U>    I've been looking over the various alignment and interpolation
U> routines to allow stacking/coadding images.  I keep coming across
U> the issue of conservation of flux by using the Jacobian of the 
U> coordinate transform.  This makes sense if the scale is changed or
U> there are higher order terms in the transformation.  However, it is
U> not clear from any of the documents I've seen whether the 
U> interpolation algorithm is also included in this transformation that 
U> the Jacobian is taken of.  Dangling preposition or no, it seems to 
U> me that with nearly critically sampled data, interpolation at 
U> the pixel level is far more important to preservation of flux
U> than the terms you expect to find in a distortion correction.  

LD> The flux conservation option in the image resampling tasks is used to
LD> take care of the effect of scale changes in the transformation.
LD> Basically it is the ratio of the before and after area of the transformed
LD> pixel, and is a constant for simple linear transformations.

U>     So -- for example, with geotran -- do any of the interpolation
U> schemes truly conserve flux?  I had tried xregister with all 
U> the interpolation algorithms and found they were all fine for
U> the high S/N star images (checked quickly with imexamine's simple
U> aperture photometry), but were all systematically off when I
U> checked stars with fairly low S/N (execpt for "nearest neighbor", 
U> which isn't really an interpolation at all).  However, there are 
U> numerous ways I might have been misled, and I'd like to set the 
U> record straight since now I find I need a rotation term which 
U> won't work with "nearest neighbor".

LD> If the scale changes are large and in the sense that say only the 4th pixel
LD> in the input image is sampled then the interpolation algorithms are not
LD> flux conserving in the sense that not all the input image information is
LD> used, to compute the output image, even though in some sense the total flux
LD> in the image is conserved by the Jacobian. In this situation I suggest that
LD> people block average or block sum the image down close to the desired scale
LD> and then do the fine-tuning with interpolation.

LD> If the scale changes are in the other direction then all the information
LD> is used and the total flux should be conserved ... 

LD> In the case of simple shifts there is no scale change and the flux should
LD> be conserved. For example a fractional shift of 0.5,0.5 is like a block
LD> average. 

LD> In any case the total flux measured after a shift should not depend on the
LD> brightness of the measured object. If you measure through small apertures
LD> however you may see some differences as the interpolation changes the shape
LD> of the object, usually by broadening it and decreasing its amplitude.  In
LD> general the higher the order of the interpolant the better it will follow
LD> high frequency features, although the higher order interpolants will also
LD> ring in the presence of undersampled data. Nearest neighbor interpolation
LD> (in xregister and imshift)  is equivalent to shifting the image by an
LD> integer amount.

LD> I would first check the imexamine measurements. You need to be aware that
LD> by default imexamine tries to find the aperture which incloses "all" of the
LD> flux, and this might be different for different brightness objects, so that
LD> you may not be measuring through the same aperture for bright and faint
LD> cases.  To force a constant measurement radius, be sure to set the rimexam
LD> pset iterations parameter to 1. Also be sure to measure through a large
LD> enough aperture to see "all" the flux in the before and after images to
LD> avoid effects caused by the object shape change. One complication that can
LD> occur if the interpolant rings and the aperture edge hits a maxima or
LD> minima. Measuring through a large aperture should remove this effect.

U>     The alternative is to figure out how to apply this STSDAS
U> drizzle function, which I presume absolutely positively conserves
U> flux, but it seems an order of magnitude more to learn.  Any
U> comments?
U> 

LD> The next release of iraf will contain support for the sinc interpolant and
LD> drizzle resampling. Sinc is the best of the interpolants at preserving high
LD> frequency information and has had a lot of testing on mosaic data. Drizzle
LD> uses the basic ST resampling algorithm but has not had a lot of use here
LD> yet within the context of the magnify and geotran tasks. If you are
LD> interested in trying it or the sinc option out you can pull over the
LD> external addon package immatchx, which supports the new interpolants.

LD> The ST drizzle function should be truly flux conserving and well tested.
LD> The only complication with it is that it is tied to wfp2 image formats, and
LD> you have to fiddle with it a bit to use it on "generic" data. I have tested
LD> magnify and geotran against it for simple linear transformations and the
LD> agreement looked good, except for minor differences in how we define the
LD> pixel origin.

----

U>     You replied to a my query about flux conservation when
U> interplating image values (e.g. using geotran or xregister).   You
U> were curious about why I had seen systematic deviations in some
U> tests I had made, and mentioned a few pertinent facts about imexamine.
U>     I tried several new tests with that knowledge, and found much more
U> consistency.  In brief, in my earlier experiments I had been using
U> iterations=3, which caused much of the trouble.  However, in case
U> somebody else becomes worried about this, you might want to
U> pass my findings along.
U>     I found:

U> 1) With iterations=1 I get similar RMS measured magnitude deviations
U> with (a) shifting the image a fractional pixel shift with center=yes, 
U> and (b) shifting the imagecur positions fractional pixel shifts with 
U> center=no.  The magnitudes from the original image vs. the shifted 
U> images (or shifted imagecur positions) showed no systematic deviation, 
U> though the RMS of all the measurements of each star was larger for 
U> faint stars.  I expect that is from the slightly different samples 
U> of the sky noise, but I haven't looked in detail.  

U> 2) With iterations=1, comparing fluxes from the original image
U> and shifted images using a constant fractional pixel shift
U> and small rotation but different interpolation versions (linear,
U> poly3, poly5, spline3, sinc, drizzle, drizzle[0.5]), the same 
U> star will often (more often than random?) show the same sign of 
U> magnitude "error" for all the interpolated images compared with 
U> the original image.  However, eyeball averaging over all the stars 
U> shows no systematic deviation.  

U> 3) With iterations=3, a repeat of (2) gives more dramatic errors and
U> more significant systematic deviations for a given star -- i.e. for
U> many of the faint stars, the magnitudes from the shifted/interpolated 
U> images minus the original image's magnitude gave the same sign of 
U> "error", and it was large.  But again, some stars had negative 
U> magnitude "errors" and some positive, and some clustered around zero, 
U> so on average no systematic error is obvious to the eye.  

U>     I think (3) indicates that what I found in my earlier tests
U> was from small number statistics.  The large differences were because
U> I used iterations=3, and by luck the few "faint" stars (a whopping 2)
U> both had the same sign deviation compared with the original image.  
U> I had assumed that while I only had 2 of these faint stars, since I
U> got similar large deviations with all the interpolations (linear, poly3,
U> poly5, spline3) there was enough evidence to indicate a major problem. 
U> I decided only "nearest" was safe.  

U>     Anyway, this little exercise gave me confidence in the
U> interpolation mechanisms (esp. test 1 above), so I thank you for 
U> your comments.

---

Richard Hook wrote on Jun 29, 1999

Following the comments about flux conservation in the context of image
coaddition and distortion correction I thought a bit of clarification of
what Drizzle does and why would be useful. Andy Fruchter and I developed 
Drizzle mostly for the HST WFPC2 but it is actually quite general in most
respects. Please note that the Drizzle task and Lindsey's "drizzle interpolant"
option are distinct and aren't equivalent - although at heart they use the
same code for determining the overlap of a general quadrilateral with a pixel
grid. Drizzle is a "forward" method in which flux from each input is 
distributed onto an output image rather than an interpolation of an input.
"Forward" methods have the advantage that each input and output pixel can be 
weighted and the weights used to combine the flux in a statistically correct 
manner as well as propagating the weights appropriately.

To come to flux again - Drizzle does NOT conserve flux! Instead it reconstructs
surface intensity. For well-sampled data and linear transformations (ie,
constant Jacobian) this doesn't make any difference but if geometric distortion
is significant or the data is markedly undersampled then it does. This process 
is also intimately entangled with how flat-fielding is done. For HST/WFPC2 and
most other CCDs flat fielding is done in such a way that the observation of a
a patch of uniformly illuminated sky will, after application of the flat field,
result in a uniform output. However, if the Jacobian varies (eg, for WFPC2 the
pixels at the edge can be a few % smaller on the sky than near the centre) then
this way of doing things means that the photometry of point sources close to 
the edge will be in error. For WFPC2 correction factors to compensate for this
effect are available. When such a distorted, flat-fielded, image is drizzled
onto a uniform output grid then BOTH surface intensity AND point-source 
photometry are correct are correct - the flat background stays the same but
the area of small objects is corrected for the Jacobian effect and hence also
their total fluxes. 

I hope this makes sense. This is a rather confusing subject.

Regards, Richard

+--------------------------------------------------------------------+
| Richard Hook, ST-ECF                 E-mail:  rhook@eso.org        |
| European Southern Observatory,                                     |
| Karl-Schwarzschild-Str. 2,           Phone:   +49 89 32006-389     |
| D-85748 Garching, Germany            FAX:     +49 89 32006-480     |
+--------------------------------------------------------------------+

---

Last post on Jun 29, 1999