[Numpy-discussion] nd_image.affine_transform edge effects

[Zachary Pincus]

Hello folks,

>> Hmm, this is worrisome. There really shouldn't be ringing on
>> continuous-tone images like Lena  -- right? (And at no step in an
>> image like that should gaussian filtering be necessary if you're
>> doing spline interpolation -- also right?)
>
> That's hard to say. Just because it's mainly a continuous-tone image
> doesn't necessarily mean it is well sampled everywhere. This depends
> both on the subject and the camera optics. Unlike the data I usually
> work with, I think everyday digital photographs (probably a photo
> scan in the case of Lena) do not generally have the detector sampling
> frequency matched to the optical resolution of the image. If that's
> true, the presence of aliasing in interpolated images depends on the
> structure of the subject and whether the scene has edges or high-
> frequency patterns in it.

Based on my reading of the two excellent Unser papers (both the one  
that ndimage's spline code is based on, and the one that Travis  
linked to), it seems like a major point of using splines for  
interpolation is *better* behavior in the case of non-band-limited  
data than the traditional 'anti-aliasing [e.g. lowpass] filter' +  
'resampling' + 'reconstruction [e.g. lowpass] filter' procedure.

That is, based on my (very limited) understanding (I understand the  
basics of the sampling theorem and can follow the Unser paper's  
update thereof, but not much more than that), in the spline case the  
spline fitting step *replaces* the anti-aliasing filter in the above  
procedure as the method for dealing with non-band-limited data. And  
the claim is that it in fact works better in many cases.

So it seems that something is wrong if the spline interpolation tools  
in ndimage only work properly in the band-limited case, or require  
lowpass filtering.


> Stefan's rotated Lena example is indeed a bit bizarre on zooming in!
> However, the artefacts are clearly localized to distinct edges, so I
> suspect this is indeed some kind of aliasing. Moreover, it looks like
> Lena has been decimated (reduced in size) prior to the rotation. That
> is definitely a good way to get artefacts, unless an anti-aliasing
> filter is applied before shrinking the image. My impression is that
> this image is probably somewhat undersampled (to understand exactly
> what that means, read up on the Sampling Theorem).

Agreed, it looks like aliasing. Nevertheless, any resampling  
procedure is supposed to deal with this internally, right? Either by  
lowpass filtering (traditional case), or by spline fitting (spline  
case as described by Unser and understood by me) -- it shouldn't be  
letting aliasing bubble through, correct?

>> The first was on Stefan's artificial data which had sharp edges, and
>> got very nasty ringing artifacts even with 3rd order splines. From
>> your recollection, is this expected behavior based on splines and the
>> nature of Stefan's image, or more likely to be a bug?
>
> Your question was aimed at Travis, so I don't want to discourage him
> from answering it :-), but looking at this in more detail, I do think
> the amplitude of the artefacts here is greater than I might expect due
> to ringing with a quadratic b-spline kernel, which I think has minima
> with amplitudes <10% of the central peak. There has to be SOME
> oscillation, but in Stefan's "rotate_artifacts" example it seems to be
> at the level of ~100%. Also, it is not present on one of the inner
> edges for some reason. So I do wonder if the algorithm in nd_image is
> making this worse than it needs to be.

Good observation! Here are cubic spline fits to a step and a delta  
function, which both have quite small amounts of ringing compared to  
what's observed -- at least on the color images. Maybe 10% ringing in  
each color channel adds up perceptually more than it does  
mathematically?

import numpy, Gnuplot, scipy.interpolate

# step
x = numpy.arange(-10, 10)
y = numpy.where(x > 0, 1, 0)
tck = scipy.interpolate.splrep(x, y, s=0, k=3)
nx = numpy.linspace(-10, 9, 200, True)
ny = scipy.interpolate.splev(nx, tck)
d = Gnuplot.Data(numpy.transpose([x, y]), with='points')
nd = Gnuplot.Data(numpy.transpose([nx,ny]), with='lines')
g = Gnuplot.Gnuplot()
g.plot(d, nd)
ny.min()
# -0.107490563074 -- largest ringing is ~10% of step size

# delta
x = numpy.arange(-10, 10)
y = numpy.where(x == 0, 1, 0)
tck = scipy.interpolate.splrep(x, y, s=0, k=3)
nx = numpy.linspace(-10, 9, 200, True)
ny = scipy.interpolate.splev(nx, tck)
d = Gnuplot.Data(numpy.transpose([x, y]), with='points')
nd = Gnuplot.Data(numpy.transpose([nx,ny]), with='lines')
g = Gnuplot.Gnuplot()
g.plot(d, nd)
ny.min()
# -0.136449610107 -- largest ringing is ~13% of impulse size

(sample plots are attached for reference).

Now let's use ndimage to rotate a step function by 45%, or zoom in on  
it:

import numpy, scipy.ndimage
o = numpy.ones((100, 50), dtype=float)
z = numpy.zeros((100, 50), dtype=float)
a = numpy.concatenate([o, z], axis=1)
b = scipy.ndimage.rotate(a, 45, order=3)
b.max()
# 1.08832325055
b.min()
# -0.0883232505527
c = scipy.ndimage.zoom(a, 5, order=3)
c.min()
# -0.107600148039
c.max()
# 1.10760014804

This (very simple) example looks exactly right. (Also, when I  
manually inspect the images in ImageJ -- which can load 32-bit  
floating point TIFFs, which is great -- they look correct). So what's  
going on with Lena and the other color/otherwise more complex images?


Zach


-------------- next part --------------
A non-text attachment was scrubbed...
Name: delta.pdf
Type: application/pdf
Size: 15861 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20070324/812e656a/attachment.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: step.pdf
Type: application/pdf
Size: 16996 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20070324/812e656a/attachment-0001.pdf>
-------------- next part --------------