Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-04-05 Thread Stefan van der Walt
Hi James

On Wed, Apr 04, 2007 at 08:29:50PM -0400, James Turner wrote:
   It looks like the last output value is produced by reflecting the
   input and then interpolating, but presumably then the first value
   should be 3.9, for consistency, not 3.1? Does that make sense?
 
 Aargh. I think I see what's happening now. The input is supposed to
 be interpolated and then reflected like this:
 
   [4  3  2  1]  -  [3.1  3.1  2.1  1.1  1.1]
 
 The problem is that there is still one value too many being
 interpolated here before the reflection takes place. Do the
 sections beginning at lines 168  178 need changing in a similar way
 to their counterparts at lines 129  139? I started looking into
 this, but I don't understand the code well enough to be sure I'm
 making the right changes...

Thanks for spotting that.  When I fix those lines, I see:

[[ 3.901   3.099   2.099   1.1002  1.8998  2.901 ]
 [ 3.901   3.099   2.099   1.1002  1.8998  2.901 ]]

I'll submit to SVN later today.  Note that I also enabled 'mirror'
mode, which works almost the same way as reflect:

Reflect: 1 2 3 4 - 1 2 3 4 4 3 2 1
Mirror:  1 2 3 4 - 1 2 3 4 3 2 1

Cheers
Stéfan
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-04-04 Thread James Turner
  OK, that was a one-line patch.  Please test to see if there are any
  subtle conditions on the border that I may have missed. I know of one
  already, but I'd be glad if you can find any others :)

Thanks, Stefan! That looks much better.

Today I finally had time to figure out the basics of SVN, make a patch and
apply your changes to my numarray version of nd_image (I'll switch to numpy
as soon as STScI does a full numpy-based release...). Your integer clipping
changes wouldn't compile under numarray unmodified, but my data are floating
point anyway, so I just applied and tested the array indexing changes.

It looks like there may still be some edge effects due to the mirroring
properties of the spline algorithm for higher orders, but the gross problem
of extrapolating 1 pixel into the mirrored data has gone away :-). I think
that's good enough for nd_image to be useful for me, but if I can find time
later it would be good to look into improving the behaviour further.

For my real data, mode=constant now seems to work well, but I also tested
some simple examples (like in this thread) using reflect and wrap. I'm
not 100% sure from the numarray manual what their correct behaviour is
supposed to be, but I noticed some things that seem anomalous. For example:

-

import numarray as N
import numarray.nd_image as ndi

I = N.zeros((2,4),N.Float32)
I[:,:] = N.arange(4.0, 0.0, -1.0)

trmatrix = N.array([[1,0],[0,1]])
troffset1 = (0.0, -0.1)

I_off1 = ndi.affine_transform(I, trmatrix, troffset1, order=1,
   mode='reflect', output_shape=(2,6))

print I
print I_off1

-

produces

[[ 4.  3.  2.  1.]
  [ 4.  3.  2.  1.]]
[[ 3.099   3.099   2.099   1.1002  1.8998
1.8998]
  [ 3.099   3.099   2.099   1.1002  1.8998
1.8998]]

It looks like the last output value is produced by reflecting the
input and then interpolating, but presumably then the first value
should be 3.9, for consistency, not 3.1? Does that make sense?

Cheers,

James.
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-28 Thread Stefan van der Walt
Hi,

I notice now that we've been having this discussion on the wrong list
-- oops!  We're nearly done, though.

On Mon, Mar 26, 2007 at 04:16:51AM -0400, James Turner wrote:
 For what it's worth, I'd agree with both of you that the numeric
 overflow should be documented if not fixed. It sounds like Stefan has
 figured out a solution for it though. If you make sense of the code in
 ni_interpolation.c, Stefan, I'd be very interested in how to make it
 calculate one less value at the edges :-).

I fixed the overflow issue in SVN.  I'd appreciate it if you could
test and let me know whether everything works as expected.  The output
values are now simply clipped to the bounds of the numeric type,
i.e. UInt8 to [0,255] etc.

As for the values at the edges, I'm still working on it.

Since we have a fundamental problem with the spline filter approach at
the borders, we may also want to look at subdivision interpolation
schemes, like the one described in this article:

http://citeseer.ist.psu.edu/devilliers03dubucdeslauriers.html

Regards
Stéfan
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-28 Thread Stefan van der Walt
On Wed, Mar 28, 2007 at 05:14:59PM +0200, Stefan van der Walt wrote:
 As for the values at the edges, I'm still working on it.

OK, that was a one-line patch.  Please test to see if there are any
subtle conditions on the border that I may have missed.  I know of one
already, but I'd be glad if you can find any others :)

Cheers
Stéfan
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-26 Thread Zachary Pincus
Thanks for the information and the paper link, James. I certainly  
appreciate the perspective, and now see why the anti-aliasing and  
reconstruction filtering might best be left to clients of a  
resampling procedure.

Hopefully at least some of the kinks in the spline interpolation (to  
date: obligate mirror boundary conditions, extra edge values, integer  
overflow) can be smoothed out.
I can't guarantee that I'll have the time or expertise to deal with  
ni_interpolation.c, but I'll try to give it a shot some time here.

Zach


On Mar 26, 2007, at 1:16 AM, James Turner wrote:

 PS... (sorry for all the posts, for anyone who isn't interested...)

 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?

 In the general case, I don't think it is appropriate for the  
 resampling
 procedure to use low-pass filtering internally to avoid artefacts,
 except perhaps when downsampling. It probably makes sense for computer
 graphics work, but there are cases where the input data are band
 limited to begin with and any degradation in resolution is  
 unacceptable.
 Where needed, I think low-pass filtering should either be the
 responsibility of the main program or an option. It's not even  
 possible
 for the resampling procedure to prevent artefacts in every case, since
 the aliasing in a badly undersampled image cannot be removed post
 factum (this is for undersampled photos rather than artificial  
 graphics,
 which I think are fundamentally different because everything is  
 defined
 on the grid, although I haven't sat down and proved it  
 mathematically).
 I'm also not sure how the procedure could decide on the level of
 smoothing needed for a given dataset without external information.

 Of course intermediate-order splines will probably keep everyone  
 happy,
 being reasonably robust against ringing effects without causing much
 smoothing or interpolation error :-).

 By the way, I think you and Stefan might be interested in a medical
 imaging paper by Lehmann et al. (1999), which gives a very nice  
 overview
 of the properties of different interpolation kernels:

 http://ieeexplore.ieee.org/Xplore/login.jsp?url=/ 
 iel5/42/17698/00816070.pdf?arnumber=816070

 For what it's worth, I'd agree with both of you that the numeric
 overflow should be documented if not fixed. It sounds like Stefan has
 figured out a solution for it though. If you make sense of the code in
 ni_interpolation.c, Stefan, I'd be very interested in how to make it
 calculate one less value at the edges :-).

 Cheers,

 James.

 ___
 Numpy-discussion mailing list
 Numpy-discussion@scipy.org
 http://projects.scipy.org/mailman/listinfo/numpy-discussion

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-24 Thread 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

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 = 

Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-24 Thread Stefan van der Walt
On Sat, Mar 24, 2007 at 01:41:21AM -0400, James Turner wrote:
 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.

Agreed, but the aliasing effects isn't not the problem here, as it
should be visible in the input image as well.  I'd expect a
third-order spline interpolation to be more smooth than a first-order
interpolant, but in the resulting images this isn't the case.

See

http://mentat.za.net/results/lena_small.png
http://mentat.za.net/results/img_rot_30_1.png (1st order spline)
http://mentat.za.net/results/img_rot_30_3.png (3rd order spline)

 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).

The artefacts arn't visible in the source image (url above).  The
image definately is a scaled down version of the original Lena -- very
interesting, btw, see

http://www.cs.cmu.edu/~chuck/lennapg/lenna.shtml

 investigating it further. One experiment might be to blur the original
 Lena with a Gaussian whose sigma is 1 pixel of the shrunken image
 before actually shrinking her, then do the rotation.

A rotation should take place without significant shifts in colour.
This almost looks like a value overflow problem.

 So I do wonder if the algorithm in nd_image is making this worse
 than it needs to be.

That is my suspicion, too.

 compare the results? I just tried doing a similar rotation in PyRAF on
 a monochrome image with a bicubic spline, and see considerably smaller
 artefacts (just a compact overshoot of probably a few % at the
 edge).

Could you apply the PyRAF rotation on the Lena given above and post
the result?

I always thought we could simply revert to using bilinear and bicubic
polygon interpolation (instead of spline interpolation), but now I
read on wikipedia:


In the mathematical subfield of numerical analysis, spline
interpolation is a form of interpolation where the interpolant is a
special type of piecewise polynomial called a spline. Spline
interpolation is preferred over polynomial interpolation because the
interpolation error can be made small even when using low degree
polynomials for the spline. Thus, spline interpolation avoids the
problem of Runge's phenomenon which occurs when using high degree
polynomials.


http://en.wikipedia.org/wiki/Spline_interpolation

also take a look at

http://en.wikipedia.org/wiki/Runge%27s_phenomenon

So much for side-stepping Runge's phenomenon :)

Cheers
Stéfan
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-24 Thread Stefan van der Walt
On Sat, Mar 24, 2007 at 03:25:38PM -0700, Zachary Pincus wrote:
 If Lena is converted to floating-point before the rotation is  
 applied, and then the intensity range is clipped to [0,255] and  
 converted back to uint8 before saving, everything looks fine.

Thanks, Zachary!  I can confirm that.

 So, is this a bug? Well, I still think so. Given that small ringing  
 is going to happen on all but the very smoothest images, and given  
 that ndimage is going to be used on non-floating-point types, it  
 would be good if there were some explicit internal clipping to the  
 data type's range. Otherwise, the ndimage resampling tools are unfit  
 for use on non-floating-point data that resides near the edges of the  
 range of the data type.

I agree.

 Though I'm not quite sure how one would structure the calculations so  
 that it would be possible to tell when over/underflow happened... it  
 might not be possible. In which case, either the tools should use  
 floating-point math at some of the steps internally (as few as  
 possible) before clipping and converting to the required data type,  
 or explicit warnings should be added to the documentation.

I think the spline interpolation already uses floating point math?
Looks like we are seeing a type conversion without range checking:

In [47]: x = N.array([1.2,200,255,255.6,270])

In [48]: x.astype(N.uint8)
Out[48]: array([  1, 200, 255, 255,  14], dtype=uint8)

I'll have to go and take a look at the code, but this shouldn't be
hard to fix -- clipping is fairly fast (don't we have a fast clipping
method by David C now?), even for extremely large images.

Regards
Stéfan
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-23 Thread Travis Oliphant
Stefan van der Walt wrote:

On Thu, Mar 22, 2007 at 04:33:53PM -0700, Travis Oliphant wrote:
  

I would rather opt for changing the spline fitting algorithm than for
padding with zeros.
 

  

 From what I understand, the splines used in ndimage have the implicit 
mirror-symmetric boundary condition which also allows them to be 
computed rapidly.  There may be ways to adapt other boundary conditions 
and maintain rapid evaluation, but it is not trivial as far as I know.  
Standard spline-fitting allows multiple boundary conditions because 
matrix inversion is used.  I think the spline-fitting done in ndimage 
relies on equal-spacing and mirror-symmetry to allow simple IIR filters 
to be used to compute the spline coefficients very rapidly.



Thanks, Travis.  I wasn't aware of these restrictions.

Would it be possible to call fitpack to do the spline fitting?  I
noticed that it doesn't exhibit the same mirror-property:
  

It doesn't.   It uses standard spline-fitting techniques (which are also 
slower).   So, yes you could call fitpack but for a large 
multi-dimensional array it's going to take longer. 

-Travis

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-23 Thread Travis Oliphant
James Turner wrote:

By the way, ringing at sharp edges is an intrinsic feature of higher-
order spline interpolation, right? I believe this kind of interpolant
is really intended for smooth (band-limited) data. I'm not sure why
the pre-filtering makes a difference though; I don't yet understand
well enough what the pre-filter actually does.
  

Yes, ringing at edges is an intrinsic feature of higher-order spline 
interpolation.   Eventually, the spline interpolant becomes a 
band-limited sinc-interpolator which will assume that edges are really 
points sampled from a sinc.  So, if you re-sample at different points 
you get the ringing effect.  But, you shouldn't see a lot of ringing 
below order 7.

The pre-filter obtains the spline-interpolation coefficients.  The 
spline assumes the continuous function represented by the samples at x_i is

f(x,y) = sum(c_ij beta^o(x-x_i) beta^o(y-y_i))

The pre-filter is computing the coefficients c_ij.  You then evaluate 
at any point you like using the continuous function implied by the 
spline.   The function beta^o is a spline function and depends on the 
order of the spline.

I'm not sure what people normally do in computer graphics, since I'm
working more with natural band-limited images, but in the case of
Stefan's rotate_artifacts example, wouldn't it be appropriate to
use the 1st- or 0th-order spline instead of the 2nd order? If your
real-life data are smooth enough, however, then in theory the
ringing with higher orders should go away.
  

Yes, that is true.  But, you really shouldn't see much ringing with a 
3rd order spline, though.

I studied these splines for a while, but my memory can fail me.  You can 
look at the papers of M. Unser for much more information.  Here is a 
link to a good review paper.

http://bigwww.epfl.ch/publications/unser0001.pdf


-Travis

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-23 Thread Zachary Pincus
Hello all,

On Mar 23, 2007, at 3:04 AM, Stefan van der Walt wrote:

 On Thu, Mar 22, 2007 at 11:20:37PM -0700, Zachary Pincus wrote:
 The actual transform operators then use these coefficients to
 (properly) compute pixel values at different locations. I just
 assumed that the pre-filtering was broken (even on natural images
 with smooth variations) because images pre-filtered with the spline
 filter didn't look like traditional pre-filtered images ought.

 IMO, ticket 213 should be closed as PEBCAK (I'll do that forthwith);
 further I'm sorry to have caused Peter to be further bothered about
 this non-issue.

 The ticket is more concerned with the usability of ndimage.
 Currently, the principle of least surprise is violated.  If you use
 the default arguments of ndimage.rotate (except for axes, which should
 still be fixed to be (1,0) by default) and rotate Lena by 30 degrees,
 you get something with green splotches on (attached to the ticket).

 IMO, we should either change the default parameters or update the
 documentation before closing the ticket.

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?)

I assumed the problem was a non-issue after some testing with a few  
natural images didn't introduce any artifacts like those, but now  
with the Lena example I'm not sure. This is frustrating, to be sure,  
mostly because of my limited knowledge on the subject -- just enough  
to be dangerous.

Zach



___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-22 Thread Stefan van der Walt
On Thu, Mar 22, 2007 at 02:41:52PM -0400, Anne Archibald wrote:
 On 22/03/07, James Turner [EMAIL PROTECTED] wrote:
 
  So, its not really a bug, its a undesired feature...
 
 It is curable, though painful - you can pad the image out, given an
 estimate of the size of the window. Yes, this sucks.

I would rather opt for changing the spline fitting algorithm than for
padding with zeros.

Stéfan
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-22 Thread Travis Oliphant
Stefan van der Walt wrote:

On Thu, Mar 22, 2007 at 02:41:52PM -0400, Anne Archibald wrote:
  

On 22/03/07, James Turner [EMAIL PROTECTED] wrote:



So, its not really a bug, its a undesired feature...
  

It is curable, though painful - you can pad the image out, given an
estimate of the size of the window. Yes, this sucks.



I would rather opt for changing the spline fitting algorithm than for
padding with zeros.
  

 From what I understand, the splines used in ndimage have the implicit 
mirror-symmetric boundary condition which also allows them to be 
computed rapidly.  There may be ways to adapt other boundary conditions 
and maintain rapid evaluation, but it is not trivial as far as I know.  
Standard spline-fitting allows multiple boundary conditions because 
matrix inversion is used.  I think the spline-fitting done in ndimage 
relies on equal-spacing and mirror-symmetry to allow simple IIR filters 
to be used to compute the spline coefficients very rapidly.

This assumption would have to be re-visited for any changes to be made.

-Travis

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-22 Thread Stefan van der Walt
On Thu, Mar 22, 2007 at 04:33:53PM -0700, Travis Oliphant wrote:
 I would rather opt for changing the spline fitting algorithm than for
 padding with zeros.
   
 
  From what I understand, the splines used in ndimage have the implicit 
 mirror-symmetric boundary condition which also allows them to be 
 computed rapidly.  There may be ways to adapt other boundary conditions 
 and maintain rapid evaluation, but it is not trivial as far as I know.  
 Standard spline-fitting allows multiple boundary conditions because 
 matrix inversion is used.  I think the spline-fitting done in ndimage 
 relies on equal-spacing and mirror-symmetry to allow simple IIR filters 
 to be used to compute the spline coefficients very rapidly.

Thanks, Travis.  I wasn't aware of these restrictions.

Would it be possible to call fitpack to do the spline fitting?  I
noticed that it doesn't exhibit the same mirror-property:

In [24]: z = scipy.interpolate.splrep([0,1,2,3,4],[0,4,3,2,1])

In [25]: scipy.interpolate.splev([0,1,2,3,4,5],z)
Out[25]: 
array([ -1.32724622e-16,   4.e+00,   3.e+00,
 2.e+00,   1.e+00,  -1.2500e+00])

Regards
Stéfan
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-22 Thread Zachary Pincus
Hello all,

 By the way, ringing at sharp edges is an intrinsic feature of higher-
 order spline interpolation, right? I believe this kind of interpolant
 is really intended for smooth (band-limited) data. I'm not sure why
 the pre-filtering makes a difference though; I don't yet understand
 well enough what the pre-filter actually does.

 I'm not sure what people normally do in computer graphics, since I'm
 working more with natural band-limited images, but in the case of
 Stefan's rotate_artifacts example, wouldn't it be appropriate to
 use the 1st- or 0th-order spline instead of the 2nd order? If your
 real-life data are smooth enough, however, then in theory the
 ringing with higher orders should go away.

James is indeed correct. This whole thing is my mistake based on a  
mis-interpretation of terminology. I've more carefully re-read the  
paper on which the ndimage spline code is based ('Splines: A Perfect  
Fit for Signal/Image Processing' by Michael Unser; it's on citeseer).  
I now (hopefully) understand that the spline pre-filter, while  
explicitly analogous to a traditional anti-aliasing pre-filter, is  
actually computing the spline coefficients via a filtering-type  
operation. While a traditional anti-aliasing filter should blur an  
image (as a band-limiting step), the fact that the anti-aliasing pre- 
filter does not is of no concern since the filtered output isn't a  
band-limited set of pixels, but a set of coefficients for a band- 
limited spline.

The actual transform operators then use these coefficients to  
(properly) compute pixel values at different locations. I just  
assumed that the pre-filtering was broken (even on natural images  
with smooth variations) because images pre-filtered with the spline  
filter didn't look like traditional pre-filtered images ought.

IMO, ticket 213 should be closed as PEBCAK (I'll do that forthwith);  
further I'm sorry to have caused Peter to be further bothered about  
this non-issue.

Zach






On Mar 22, 2007, at 8:13 PM, James Turner wrote:

 By the way, ringing at sharp edges is an intrinsic feature of higher-
 order spline interpolation, right? I believe this kind of interpolant
 is really intended for smooth (band-limited) data. I'm not sure why
 the pre-filtering makes a difference though; I don't yet understand
 well enough what the pre-filter actually does.

 I'm not sure what people normally do in computer graphics, since I'm
 working more with natural band-limited images, but in the case of
 Stefan's rotate_artifacts example, wouldn't it be appropriate to
 use the 1st- or 0th-order spline instead of the 2nd order? If your
 real-life data are smooth enough, however, then in theory the
 ringing with higher orders should go away.

 Sorry if I'm stating the obvious and missing the real point! I just
 wanted to make sure this hasn't been overlooked... Likewise, sorry I
 didn't mention this before if it does turn out to be relevant. Let
 me know if you want references to explain what I said above.

 James.

 ___
 Numpy-discussion mailing list
 Numpy-discussion@scipy.org
 http://projects.scipy.org/mailman/listinfo/numpy-discussion

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-15 Thread James Turner
Hi Stefan,

Thanks for the suggestions!

 Is this related to
 
 http://projects.scipy.org/scipy/scipy/ticket/213
 
 in any way?

As far as I can see, the problems look different, but thanks for
the example of how to document this. I did confirm that your example
exhibits the same behaviour under numarray, in case that is useful.

 Code snippets to illustrate the problem would be welcome.

OK. I have had a go at producing a code snippet. I apologize that
this is based on numarray rather than numpy, since I'm using STScI
Python, but I think it should be very easy to convert if you have
numpy instead.

What I am doing is to transform overlapping input images onto a
common, larger grid and co-adding them. Although I'm using affine_
transform on 3D data from FITS images, the issue can be illustrated
using a simple 1D translation of a single 2D test array. The input
values are just [4., 3., 2., 1.] in each row. With a translation of
-0.1, the values should therefore be something like
[X, 3.1, 2.1, 1.1, X, X], where the Xs represent points outside the
original data range. What I actually get, however, is roughly
[X, 3.1, 2.1, 1.0, 1.9, X]. The 5th value of 1.9 contaminates the
co-added data in the final output array. Now I'm looking at this
element-by-element, I suppose the bad value of 1.9 is just a result
of extrapolating in order to preserve the original number of data
points, isn't it? Sorry I wasn't clear on that in my original post
-- but surely a blank value (as specified by cval) would be better?

I suppose I could work around this by blanking out the extrapolated
column after doing the affine_transform. I could calculate which is
the column to blank out based on the sense of the offset and the
input array dimensions. It seems pretty messy and inefficient though.
Another idea is to split the translation into integer and fractional
parts, keep the input and output array dimensions the same initially
and then copy the output into a larger array with integer offsets.
That is messy to keep track of though. Maybe a parameter could
instead be added to affine_transform that tells it to shrink the
number of elements instead of extrapolating? I'd be a bit out of my
depth trying to implement that though, even if the authors agree...
(maybe in a few months though).

Can anyone comment on whether this problem should be considered a
bug, or whether it's intended behaviour that I should work around?

The code snippet follows below. Thanks for your patience with
someone who isn't accustomed to posting questions like this
routinely :-).

James.

-

import numarray as N
import numarray.nd_image as ndi

# Create a 2D test pattern:
I = N.zeros((2,4),N.Float32)
I[:,:] = N.arange(4.0, 0.0, -1.0)

# Transformation parameters for a simple translation in 1D:
trmatrix = N.array([[1,0],[0,1]])
troffset = (0.0, -0.1)

# Apply the offset to the test pattern:
I_off1 = ndi.affine_transform(I, trmatrix, troffset, order=3, mode='constant',
  cval=-1.0, output_shape=(2,6))

I_off2 = ndi.affine_transform(I, trmatrix, troffset, order=3, mode='constant',
  cval=-1.0, output_shape=(2,6), prefilter=False)

# Compare the data before and after interpolation:
print I
print I_off1
print I_off2

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion


Re: [Numpy-discussion] nd_image.affine_transform edge effects

2007-03-14 Thread Stefan van der Walt
Hi James

On Fri, Mar 09, 2007 at 08:44:34PM -0300, James Turner wrote:
 Last year I wrote a program that uses the affine_transform()
 function in numarray to resample and co-add datacubes with WCS
 offsets in 3D. This function makes it relatively easy to handle
 N-D offsets and rotations with a smooth interpolant, which is
 exactly what I wanted. However, I am finding that the results
 contain edge effects which translate into artefacts in the
 middle of the final mosaic.

Is this related to

http://projects.scipy.org/scipy/scipy/ticket/213

in any way?

Code snippets to illustrate the problem would be welcome.

Regards
Stéfan
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion