590 likes | 748 Views
Pragmatic approaches to image resampling. Tom McGlynn NASA/GSFC. Outline. Why do we resample? Isn’t resampling a solved problem? Some approaches to resampling. A new algorithm for exact-area resampling using clipping. How should astronomers choose resampling algorithms?.
E N D
Pragmatic approaches to image resampling Tom McGlynn NASA/GSFC
Outline • Why do we resample? • Isn’t resampling a solved problem? • Some approaches to resampling. • A new algorithm for exact-area resampling using clipping. • How should astronomers choose resampling algorithms? IPAM Workshop Jan 27, 2004
Why do we need to resample? • Display – transform image into ‘standard’ form • Undo warps and distortions • Transform to standard frame • Resizing • Magnification and minification • Rotation • Image comparison – transform one image to match another • Similar operations as in display, but final frame may not be ‘standard’ • Mosaicking • Building sky region and all sky images • Image arithmetic • Dither additions, image differencing, speckle analysis • Can have most extreme requirements on accuracy of reconstruction but often involves very similar images. Different problems levy different requirements for accuracy and robustness. Many resampling problems involve substantial changes in geometry. IPAM Workshop Jan 27, 2004
Resampling examples SkyView transforms the EGRET all sky map in Galactic coordinates to Equatorial coordinates. IPAM Workshop Jan 27, 2004
WMAP data must be transformed from HEALpix formats to human friendly ones. IPAM Workshop Jan 27, 2004
Transforming SOHO/EIT images to Cartesian coordinates. And that’s before handling differential rotation… On Resampling of Solar Images Solar Physics, 2003, C.E. DeForest IPAM Workshop Jan 27, 2004
SkyView dynamically mosaics surveys data for display or image comparison IPAM Workshop Jan 27, 2004
SkyMorph image addition. Three co-added NEAT images taken at 20 minute intervals. (1999 discovery images for asteroid 9460) DSS image (1950) Adding NEAT images for asteroid discovery using SkyMorph. Other survey datasets can also be searched for pre-discovery IPAM Workshop Jan 27, 2004
Classical image resampling • Image reconstruction The continuous image is regenerated by some method, often by interpolating between the grid point by using an interpolation kernel. • Sampling The new grid is constructed by sampling the reconstructed function appropriately. • Filters and transformations Intermediate filters and transformations may be applied at various points in the process. IPAM Workshop Jan 27, 2004
Ideal Reconstruction • A band-limited function can be exactly reconstructed by convolution with the appropriate kernel. IPAM Workshop Jan 27, 2004
Shannon says…By hypothesis f… …has a Fourier transform, F bandlimited to 1/2T. -1/(2T) 1/(2T) Pixelate: Multiply by Comb function) with separation of T… … means convolve the Fourier transform with the Fourier Transform of Comb, but this is just another Comb function with separation of 1/T which just relplicates the Fourier Transform. 5T So mutiplying by a boxcar gets back original F … … but that’s equivalent to convolving in space with Fourier transform of boxcar, sinc(). So we just need to use sinc interpolation to get back original image.
Solved problem? • So sinc function is the optimal filter. Just reconstruct the image using sinc interpolation and resample as needed! IPAM Workshop Jan 27, 2004
But… real images are finite • Chips typically have 1-4K pixels with larger detectors using arrays of chips. • Edge effects need to be considered … • Copying image, image reversal, constant value, 0 beyond edge • How do we handle missing pixels? • Typical CD song has ~100K samples before we begin hearing the music. IPAM Workshop Jan 27, 2004
…Is the image well sampled? • Is sampling rate adequate? How do we tell? • What can cause high frequencies? • Noise and photon statistics • Features • Processing • Filtering and compression – flat can mean high frequencies. • The sky • 1 Gigapixel Nyquist-sampled image for 0.001” resolution covers only 16”. • Do users prefilter images? IPAM Workshop Jan 27, 2004
…The sky isn’t a plane • Proper basis for sphere is spherical harmonics. • CMB studies use HEALPix pixels to make this easier. • Are HTM grids going to be used? • Derivation of sinc function kernel assumes rectangular grid. Other functions are appropriate for other samplings. • What are the effects of ignoring curvature? IPAM Workshop Jan 27, 2004
…Point versus Area resampling • Pixels have finite extent. • In 2-d case the relationship between the original and resampling pixels can be complex. • 1-D and 2-D resampling share ‘calculus’ but 2-D geometry is more complex. • If input and output pixels are similar we can presume the function being reconstructed is the flux convolved with the pixel mask. • How do we handle resampling pixels with complex shapes? IPAM Workshop Jan 27, 2004
What makes resampling hard? • Calculus • Undersampling • Noise, features, missing pixels • Constraints on the output range (avoidance of negative values, integer valued functions) • Geometry • Differing projections and coordinate systems and orientations. • Scale changes • Non rectangular pixels. • Variations over the image IPAM Workshop Jan 27, 2004
Alternatives • Image regeneration • Is the image a derived product? • Pixels added? • Models? • Retake the image • Adapt requirements to minimize projection issues • Often not feasible or desirable. IPAM Workshop Jan 27, 2004
Some reconstruction kernels The background image is reconstructed by convolution with the given kernel. The sinc function is the ‘optimal’ reconstruction kernel for a well sampled image. IPAM Workshop Jan 27, 2004
Nearest Neighbor Universally panned… • Bad geometric distortions • Suspect photometry • Poor for magnification but even the least elegant techniques may have their place. • Fast • Histogram of values can be preserved. • Sometimes exactly invertible. • Integers stay integers • Classifications (does vegetation+water=city) • Photon counting algorithms • Retention of discontinuities • Man-made boundaries • Constellation maps IPAM Workshop Jan 27, 2004
Bilinear Interpolation • Fast • Better at preserving astrometry and photometry But • Smooths and blurs the image IPAM Workshop Jan 27, 2004
‘Higher order’ interpolation • Best job reproducing well sampled images • Limits blurring But… • Sometimes slower • Can be sensitive to ‘features’ or undersampling. - negative values when resampling. • Harder to handle missing pixels • Gazillions of choices • Interpolating or approximating • Polynomial approximations • Splines • Local support polynomials • Truncated sinc • Windowed sinc: Lanczos, Hamming, … • Image support for calculation of kernel • Image support for resampling calculations • Very high order methods (e.g., exact sinc) can be very slow IPAM Workshop Jan 27, 2004
More choices: Interpolation versus Redistribution • Interpolation • Output flux at a pixel is computed as a sum of weighted input pixels using reconstruction kernel. • Tends to assure differential properties of image (continuity, derivatives) • Redistribution (Drizzle, exact area sampling) • Input pixel’s flux is distributed over output pixels using redistribution kernel. • Can ensure integral properties of image (flux or flat field). • Easy to handle missing pixels or other discontinuities • Global integrals can be conserved to the limits of arithmetic precision IPAM Workshop Jan 27, 2004
Exact Area Resampling Calculate the resampling pixels as the weighted averages of the pixels they cover (or weighted sums for extensive data) weighting IPAM Workshop Jan 27, 2004
Approaches • Delaunay Triangulation • Using suitable collection of points one finds a set of triangles where each triangle belongs to just one input and output pixel. • Hideously slow… • Girard’s theorem (Montage: Berriman and Good) • Work on the celestial sphere and use Girard’s theorem to calculate area – a `practical’ example of using parallel transport of vectors. • Very, very slow. • Clipping • How much of the resampling pixels can you view through the window of the input pixels? • Comparable speed to other high-order methods. • Equivalent to algorithm used within Drizzle? IPAM Workshop Jan 27, 2004
Exact area resampling using clipping. For each resampling pixel… First find all the pixels that may overlap the resampling pixel, by looking at the range of the resampling pixel corners. Now for each candidate original pixel… IPAM Workshop Jan 27, 2004
Clip the resampling pixel on each edge. Clip the resampling pixel by each of four clipping boundaries. (Sutherland-Hodgman algorithm but the extensive clipping literature suggests more efficient approaches, e.g., Liang-Barsky). 1 2 3 4 IPAM Workshop Jan 27, 2004
Clipping by a infinite line v3 c2 For i = 1 to n if vi-1 is outside if vi is inside emit crossing point (e.g., c2) and vi else if vi is outside emit crossing point (e.g., c1) else emit vi v4 v2 c1 v1 Inside Outside Input (v1,v2,v3,v4) Output (v1,c1,c2,v3,v4) IPAM Workshop Jan 27, 2004
Triangulating the overlap D C The area of the overlap polygon ABCDE is easily computed as the sum of the triangles ABC, ACD, ADE E A B IPAM Workshop Jan 27, 2004
Normalization • Intensive images • Track total flux and overlapped area for each pixel and use ratio for pixel value. • Preserves flat fielding • Extensive images • Just add flux. • Preserves total flux IPAM Workshop Jan 27, 2004
Features • Convexity of pixels allows simplification of algorithm (since convex regions clip convex regions to convex regions). • Clipping on rectangular grid is especially easy but either grid could be triangular or hexagonal or even discontiguous. • Don’t have to make clipping window the same as pixel boundaries. • If clipping window is smaller we get Drizzle-like algorithm • If clipping window is larger than pixels we have a variable box-car along with resampling. • Easy to accommodate ‘bad’ pixels or regions. • Symmetry between resampling and original grid • Can resample in the `convenient’ direction if transformation is easier one way or the other. Not the most ‘accurate’ algorithm, but it can be extremely robust. “Best worst-case resampler” IPAM Workshop Jan 27, 2004
Exact area sampling Drizzle redistribution Boxcar smoothing with resampling Step pyramid kernel with different fraction of flux in each box Possible clipping ‘kernels’ Adjusting the size of the rectangle changes how the algorithm considers the flux in the pixel to be distributed. IPAM Workshop Jan 27, 2004
How can astronomers decide which resampling algorithm to use? • Accuracy • Does resampling affect astronomical measurements? • Detection, Astrometry, Photometry, Morphology, Resolution • Cost (CPU, memory) • Complexity • Coding and comprehension costs • Robustness • Point versus area resampling • How well does it work when in hard resampling situations? IPAM Workshop Jan 27, 2004
Need for a roadmap Goals: Mosaicking, Resizing, Displaying, Subtracting, Dithering, Comparing, Undistorting, Projecting Issues: Undersampling Pixel distortions Spherical geometry Edges Prefiltering Bad pixel regions Compression … Users Methods: Nearest neighbor Linear interpolation Splines Polynomial kernels Gaussian Sinc Hamming Lanczos Exact area … Traits: Accuracy Speed Comprehensibility Availability Robustness … IPAM Workshop Jan 27, 2004
Analysis in other regimes • Quantitative Comparison of Sinc-Approximating Kernels for Medical Image Interpolation • Erik H. W. Meijering, Wiro J. Niessen, Josien P. W. Pluim, Max A. Viergever • Round-robin resampling – a sequence of resamplings leading back to the original image. • Good for determining the ‘best’, but not for assessing the cost • Truncated sinc resamplers are the worst. • Variety of higher order resamplers do pretty well. • But the images don’t look much like typical astronomical data) IPAM Workshop Jan 27, 2004
… • Image Reconstruction by Convolution with • Symmetrical Piecewise nth-Order • Polynomial Kernels • Erik H. W. Meijering, Karel J. Zuiderveld, Max A. Viergever • IEEE Transactions on Image Processing, vol. 8, no. 2, February 1999, pp. 192, 201. • Going beyond third order kernels doesn’t seem to buy one anything. IPAM Workshop Jan 27, 2004
Comparing image resamplers via a model of the human vision system Richard Harvey, Stephen King, Richard Aldridgeand J. Andrew Banghamh A big driver in commercial applications is rendering of text.
Even astronomical references often compare samplers in qualitative terms. Features in Lanczos resampling. NN LI Lcz2 Lcz3 Lcz4 SWarp v2:0 User's guide E. Bertin Moire patterns in linearly interpolated resampled images IPAM Workshop Jan 27, 2004
Measuring accuracy quantitatively • Single resampling • Astronomicalish images • How does noise affect resampling? • Two resampling scenarios: • Small Pixel: Rotation, final pixel size/original pixel size=1.1 • Big Pixel: Rotation, final pixel size/original pixel size=2.5 • Use Sextractor to estimate parameters of original and resampled objects. • Model image with 100 gaussian objects and variable noise • Use DSS image for reality check IPAM Workshop Jan 27, 2004
Test Images 0: 0 (Number: Noise) 5: 0.006 10: 0.03 15: 0.2 20: 1.0
North Pole image DSS image IPAM Workshop Jan 27, 2004
Samplers Tested • Interpolation algorithms • NN – Nearest Neighbor • LI – Linear Interpolation • LCZn – A Lanczos n-lobe interpolator (n=3,4) • SPn – An n’th order spline (n=3,4) • Redistribution algorithms • CL – Clipping exact area • CL0.5 – Clipping using a window of half the size of the pixel (similar to Drizzle) • MN –Montage exact area algorithm, resampling done on celestial sphere. IPAM Workshop Jan 27, 2004
Small Pixel Detection Did we detect all of the objects detected in the unresampled data? Boxes give the results for the real image looking at the first six samplers listed. The dashed line gives the measurement in the original unresampled image. IPAM Workshop Jan 27, 2004
Small Pixel Astrometry What was the average offset of the resampled image from the measurement in the original image? IPAM Workshop Jan 27, 2004
Small Pixel Photometry How much did the flux change when we resampled? IPAM Workshop Jan 27, 2004
Small Pixel Morphology The image modeled circular gaussians. What is the average axis ratio (a-b)/a measured in the resampled data? Note that graphs compare with model and not measured values. IPAM Workshop Jan 27, 2004
Small Pixel Blurring Is the resampled image blurred by the resampling? How much larger is the it? IPAM Workshop Jan 27, 2004
Big Pixel Detections IPAM Workshop Jan 27, 2004
Big Pixel Astrometry IPAM Workshop Jan 27, 2004
Big Pixel Photometry IPAM Workshop Jan 27, 2004