[Numpy-discussion] nd_image.affine_transform edge effects
Zachary Pincus
zpincus@stanford....
Sat Mar 24 14:56:05 CDT 2007
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
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
More information about the Numpy-discussion
mailing list