Tags: algorithm, amimplementing, complex, following, interesting, interp2, interpolation, matlab, numbers, programming, requires, running, stepsfimage

Interpolation over complex numbers

On Programmer » Matlab

10,168 words with 9 Comments; publish: Tue, 29 Apr 2008 14:55:00 GMT; (200109.38, « »)

Hi,

I'm running into an interesting problem with interp2. I am

implementing an algorithm that requires the following steps:

fimage = FFT2(image)

POLARTRANSFORM(fimage)

...

CARTESIANTRANSFORM(fimage)

image2 = IFFT2(fimage)

What I've discovered is the process of polar/cartesian transforming

introduces serious errors into the phase component of the FFT'd image

but NOT into the magnitude component. I do the polar/cartesian polar

transforms using interp2 as follows:

function [pim, radius, theta] = c2p(im, N, rlower, rupper)

[rows, cols] = size(im);

cx = cols/2.0 + 0.5; % Add 0.5 because index starts at 1

cy = rows/2.0 + 0.5;

[theta, radius] = meshgrid(linspace(0, 2*pi, N),

logspace(log10(rlower), log10(rupper), N));

xi = radius.*cos(theta) + cx; % Locations in image to interpolate

data

yi = radius.*sin(theta) + cy; % from.

[x,y] = meshgrid([1:cols],[1:rows]);

% Older versions of interp2 may not handle complex numbers correctly!

% Need to interpolate the magnitude and phase separately to avoid

problems.

pim = interp2(x, y, abs(im), xi, yi, 'cubic').*exp(i*interp2(x, y,

angle(im), xi, yi, 'cubic'));

mask = isnan(pim);

pim(mask) = 0;

Is there a reason I am getting junk of the phase interpolation but

not for the magnitude?

All Comments

Leave a comment...

  • 9 Comments
    • "Elvis Dieguez" <edieguez.matlab.todaysummary.com.ieeeREMOVEME.org> wrote in message

      news:ef4cfe8.-1.matlab.todaysummary.com.webcrossing.raydaftYaTP...

      > Hi,

      > I'm running into an interesting problem with interp2. I am

      > implementing an algorithm that requires the following steps:

      > fimage = FFT2(image)

      > POLARTRANSFORM(fimage)

      > ...

      > CARTESIANTRANSFORM(fimage)

      > image2 = IFFT2(fimage)

      > What I've discovered is the process of polar/cartesian transforming

      > introduces serious errors into the phase component of the FFT'd image

      > but NOT into the magnitude component. I do the polar/cartesian polar

      > transforms using interp2 as follows:

      > function [pim, radius, theta] = c2p(im, N, rlower, rupper)

      > [rows, cols] = size(im);

      > cx = cols/2.0 + 0.5; % Add 0.5 because index starts at 1

      > cy = rows/2.0 + 0.5;

      > [theta, radius] = meshgrid(linspace(0, 2*pi, N),

      > logspace(log10(rlower), log10(rupper), N));

      > xi = radius.*cos(theta) + cx; % Locations in image to interpolate

      > data

      > yi = radius.*sin(theta) + cy; % from.

      > [x,y] = meshgrid([1:cols],[1:rows]);

      > % Older versions of interp2 may not handle complex numbers correctly!

      > % Need to interpolate the magnitude and phase separately to avoid

      > problems.

      > pim = interp2(x, y, abs(im), xi, yi, 'cubic').*exp(i*interp2(x, y,

      > angle(im), xi, yi, 'cubic'));

      > mask = isnan(pim);

      > pim(mask) = 0;

      > Is there a reason I am getting junk of the phase interpolation but

      > not for the magnitude?

      Without thinking too hard or looking at your code, my guess is a problem

      with phase unwrapping.

      #1; Tue, 29 Apr 2008 14:56:00 GMT
    • Elvis Dieguez wrote:

      > I'm running into an interesting problem with interp2. I am

      > implementing an algorithm that requires the following steps:

      > fimage = FFT2(image)

      > POLARTRANSFORM(fimage)

      > ...

      > CARTESIANTRANSFORM(fimage)

      > image2 = IFFT2(fimage)

      > What I've discovered is the process of polar/cartesian transforming

      > introduces serious errors into the phase component of the FFT'd

      > image

      > but NOT into the magnitude component. I do the polar/cartesian

      > polar

      > transforms using interp2 as follows:

      As Ken pointed out you run into a phase unwrapping problem. Phase

      information has a discontinuity at 2*pi which is not handled by

      imresize. For example: if two neighboring pixels have the angles 0

      and 2*pi imresize may use the average value for the output pixel.

      This obviously introduces a huge error. Instead of using magnitude

      and phase images as input to imresize you could use real and

      imaginary part images. These images can be resized savely.

      -Herbert

      #2; Tue, 29 Apr 2008 14:57:00 GMT
    • Herbert Ramoser wrote:

      >

      > Elvis Dieguez wrote:

      > transforming

      > As Ken pointed out you run into a phase unwrapping problem. Phase

      > information has a discontinuity at 2*pi which is not handled by

      > imresize. For example: if two neighboring pixels have the angles 0

      > and 2*pi imresize may use the average value for the output pixel.

      > This obviously introduces a huge error. Instead of using magnitude

      > and phase images as input to imresize you could use real and

      > imaginary part images. These images can be resized savely.

      Don't know much about 2D, but in 1D you get unwanted effects if you

      interpolate the real and imaginary parts rather than the magnitude

      and phase.

      #3; Tue, 29 Apr 2008 14:58:00 GMT
    • On Feb 7, 4:07 am, "Steve Amphlett"

      <Firstname.Lastname.matlab.todaysummary.com.where_I_work.com> wrote:

      > Don't know much about 2D, but in 1D you get unwanted effects if you

      > interpolate the real and imaginary parts rather than the magnitude

      > and phase.

      >

      Really?

      I don't believe it.

      But if you're right, then all the work by hundreds of researchers on

      developing global tide models using data from TOPEX/Poseidon

      oceanographic satellite is down the drain.

      Can you demonstrate a single case where interpolation on real and

      imaginary parts gives "unwanted" effects.

      Seems to me interpolation on real and imag parts is the ONLY valid way

      to interpolate in 1D or 2D without getting screwed up with phase, as

      mentioned by others.

      #4; Tue, 29 Apr 2008 14:59:00 GMT
    • NZTideMan wrote:

      >

      > On Feb 7, 4:07 am, "Steve Amphlett"

      > <Firstname.Lastname.matlab.todaysummary.com.where_I_work.com> wrote:

      > you

      > magnitude

      > Really?

      > I don't believe it.

      > But if you're right, then all the work by hundreds of researchers

      > on

      > developing global tide models using data from TOPEX/Poseidon

      > oceanographic satellite is down the drain.

      > Can you demonstrate a single case where interpolation on real and

      > imaginary parts gives "unwanted" effects.

      > Seems to me interpolation on real and imag parts is the ONLY valid

      > way

      > to interpolate in 1D or 2D without getting screwed up with phase,

      > as

      > mentioned by others.

      Ok, let's say you have two sine waves at very similar frequencies and

      amplitudes but significantly different phases and want to smoothly

      move from one to the other. If you interpolate the real and

      imaginary parts independently the amplitude will dip as you move from

      one point to the other, whereas you'd probably want the amplitude and

      phase to gradually adjust each in a monotonic way. What you really

      want in the complex plane is a spiral from one point to the other,

      not a straight line.

      #5; Tue, 29 Apr 2008 15:00:00 GMT
    • Steve Amphlett wrote:

      >

      > Ok, let's say you have two sine waves at very similar frequencies

      > and

      > amplitudes but significantly different phases and want to smoothly

      > move from one to the other. If you interpolate the real and

      > imaginary parts independently the amplitude will dip as you move

      > from

      > one point to the other, whereas you'd probably want the amplitude

      > and

      > phase to gradually adjust each in a monotonic way. What you really

      > want in the complex plane is a spiral from one point to the other,

      > not a straight line.

      ... I'm not advocating averaging or any other DSP calcs in general

      using amplitude and phase. Just highlighting a case where unusual

      interpolation requirements may lean toward it. As in other areas,

      one curve that's mathematically in the middle of two others may not

      be what's desired.

      #6; Tue, 29 Apr 2008 15:01:00 GMT
    • On Feb 7, 9:48 am, "Steve Amphlett" <m....matlab.todaysummary.com.home.now> wrote:

      > Steve Amphlett wrote:

      >

      > ... I'm not advocating averaging or any other DSP calcs in general

      > using amplitude and phase. Just highlighting a case where unusual

      > interpolation requirements may lean toward it. As in other areas,

      > one curve that's mathematically in the middle of two others may not

      > be what's desired.

      In your (pathological) case, I'd interpolate on real and imag on one

      signal, then the other and superpose.

      #7; Tue, 29 Apr 2008 15:02:00 GMT
    • NZTideMan wrote:

      >

      > In your (pathological) case, I'd interpolate on real and imag on

      > one

      > signal, then the other and superpose.

      Pathological? It comes up all the time in the simulation world. If

      you have noise predictions at a set of finite (engine) speeds and

      want to synthesise a slow transient through them, this is exactly

      what you have.

      #8; Tue, 29 Apr 2008 15:03:00 GMT
    • Herbert Ramoser wrote:

      I've tried doing the interpolation over real and imaginary parts and

      the errors are different but just as bad. I believe you are correct

      in diagnosing the problem (i.e. the phase discontinuity). I've tried

      using 'unwrap' but that doesn't solve the problem. Could it be that

      'unwrap' works column or row wise thereby still leaving

      discontinuities on the 2D surface?

      > As Ken pointed out you run into a phase unwrapping problem. Phase

      > information has a discontinuity at 2*pi which is not handled by

      > imresize. For example: if two neighboring pixels have the angles 0

      > and 2*pi imresize may use the average value for the output pixel.

      > This obviously introduces a huge error. Instead of using magnitude

      > and phase images as input to imresize you could use real and

      > imaginary part images. These images can be resized savely.

      > -Herbert

      #9; Tue, 29 Apr 2008 15:04:00 GMT