http://hmi.stanford.edu/doc/Tech_Notes/HMI-TN-2004-004-Interpolation/interpolate.html
On spatial interpolation of HMI filtergrams
This note describes the accuracy and speed of various methods for
spatial image interpolation when applied to a simulated HMI
filtergram. Such methods will be needed if it is deemed necessary to
co-register HMI filtergrams as part the level-1 observable calculation
to eliminate the effect of (differential) solar rotation.
It might also be necessary to interpolate the signal in time to
eliminate "acceleration effects" due to the solar p-modes signal on
the solar surface. That interpolation can probably be accomplished
with a constant FIR filter, optimized for signals with spectral
content similar to solar p-mode noise. This problem will not be
considered further here.
Methods tested
--------------
I implemented or optimized existing code for:
a) a convolutional cubic spline interpolator used in various MDI
processing (the function ccint2) due to Keys [1],
see Keys_cubic_interp.pdf.
b) a fourth order accurate version of Keys' method,
c) a true B-spline interpolator based on a fast filtering
algorithm developed by Unser et al [2-4]. The code currently
implements B-splines of order 2 through 5.
See Unser91.pdf for details.
A straightforward tensor product (separable) form of the filters is
used to extend the methods to 2-D.
Test procedure
--------------
To test the accuracy I generated a simulated "HMI filtergram" from a
512x512 La Palma image (taken at 4305 angstrom with 0.083" pixels)
provided by Ted Tarbell. I applied an MTF correction and pixel
integration corresponding to HMI's 14cm aperture at 6173 angstrom with
0.5" pixel size:
HMI_filtergram.png
Power spectra of these images after subtracting mean:
HMI_filtergram_pow.png
Filtergrams corresponding to the original images shifted in x and y by
delta=0,1,2,...,6 pixels were created. These were interpolated back to
delta=0 by the various methods above and compared to the true unshifted
filtergram.
Measured accuracy
-----------------
The relative rms interpolation errors (rms(interp-true)/rms(true)) as a
function of interpolation distance delta is plotted here:
interp_error_HMI.png
The maximal relative rms error corresponding to a shift of half a pixel
is tabulated here:
method max rms error (delta=0.5 pixels)
-----------------------------------------------------------------------
0 fcint(3) 2.9499e-03
1 fcint(4) 1.8778e-03
2 B-spline(2) 1.7922e-03
3 B-spline(3) 1.1648e-03
4 B-spline(4) 6.8938e-04
5 B-spline(5) 5.0010e-04
The input image was normalized to an average of 1, so the errors can be
compared to the rms photon shot noise of ~1/sqrt(1.5e5)= 2.6e-3. With
the existing MDI code the interpolation noise exceeds the shot noise
while the 5th order B-spline interpolator has a maximum rms error that
is about 20% of the shot noise.
Power spectra of the interpolation errors are plotted in the files below
(scroll down to see the power spectrum):
interp_err_pow_0.png
interp_err_pow_1.png
interp_err_pow_2.png
interp_err_pow_3.png
interp_err_pow_4.png
interp_err_pow_5.png
where the last digit before the filename extension refers to the row
number in the table above. All the schemes tested use piecewise
polynomial interpolants, so the interpolation error is predominantly
at high spatial frequency. The low error in the "corners" is simply
due to the MTF.
Timing results
--------------
Here are the timing results for interpolating a 4096x4096 image stored
in single precision floating point. The test was run on my workstation
(helios) which has a 3.06 GHz Pentium 4 CPU (helios). All programs were
compiled with the Intel C compiler version 7.0 with options "-O3 -ipo
-march=pentium4".
method time (sec.)
--------------------------------------
0 fcint(3) 1.38
1 fcint(4) 2.12
2 B-spline(2) 1.57+0.86
3 B-spline(3) 1.56+1.09
4 B-spline(4) 1.91+1.44
5 B-spline(5) 1.90+1.74
* ccint2(3) 2.34
Things to notice:
1. ccint2 and fcint(3) produce identical results, fcint(3) is simply
and hand tuned version of the ccint2 code from MDI's libM.
2. The time for the B-spline interpolator is reported in the form
X+Y, where X is the time taken to compute the B-spline coefficients
of the image and Y is the time to perform the actual interpolation
by convolution. The B-spline coefficients are independent of the
coordinate to which one interpolates and can therefore be computed
once and for all as each filtergram is read into the level-1
processing module. Only the Y cost is added for each additional
interpolation that is performed on the same filtergram.
3. The interpolation done here was a simple shift in x,y. For
Co-registering HMI filtergrams the transformation will obviously have
a more complicated form and some additional time must be added to
compute the pixel coordinates from orbit information, models of
differential solar rotation etc. However, if the deformation (or the
non-trivial part of it) is slowly varying it should be possible to
pre-compute (tabulat) pixel coordinates or even interpolation weights
or determine them by simple interpolation in pre-tabulated data.
[1] R. G. Keys,"Cubic Convolution Interpolation for Digital Image
Processing", IEEE Trans. Acoustics, Speech and Signal
Processing, Vol. 29 no. 6, December 1981, pp. 1153-1160.
[2] M. Unser,
"Splines: A Perfect Fit for Signal and Image Processing,"
IEEE Signal Processing Magazine, vol. 16, no. 6, pp. 22-38,
November 1999.
[3] M. Unser, A. Aldroubi and M. Eden,
"B-Spline Signal Processing: Part I--Theory,"
IEEE Transactions on Signal Processing, vol. 41, no. 2,
pp. 821-832, February 1993.
[4] M. Unser, A. Aldroubi and M. Eden,
"B-Spline Signal Processing: Part II--Efficient Design and
Applications," IEEE Transactions on Signal Processing, vol. 41,
no. 2, pp. 834-848, February 1993.