Re: FIXED (again): INTERP2 problem

[Pierre Baldensperger]

Paul Kienzle wrote:

> I'm happy to support XI and YI arbitrary, however according to the docs
> on the web, interp2(z,[1:5],[1:6]) should return a 5x6 matrix. Do you 
> have
> reason to believe otherwise?


Well, yes, sorry... Which web docs are you referring to?
Matlab behaves differently depending whether YI is a row or column vector:

- interp2(Z,[1:5],[1:6]) should not work
- interp2(Z,[1:5],[1:6]') should indeed auto-grid and return a 5x6 matrix

I have no Matlab handy to check that, but that is how I understood the
Matlab help for interp2, and that is also IMO the way it should ultimately
behave. Why would we forbid anybody to perform a 2D-to-1D interpolation
with that function (that is: entering a 1D-list of XI and YI and returning
a 1D vector ZI)?

In the present case, I'm concerned with limiting the field of application
of that routine, just because we want to avoid "lazy" users a (_trivial_)
preliminary call to meshgrid.

So if you really want to restore this "auto-meshgridding" behaviour, I'd
suggest you trigger it the "Matlab way" (or at least the way I understand
it ;-), that is:

- if XI and YI have the same size and orientation, do _not_ meshgrid
- if XI and YI have the same orientation but different sizes, _error_
- if XI and YI have different orientations, auto-meshgrid

Of course, that may cause programs relying on the "auto-meshgridding" to
fail when they provide XI and YI with the same orientation (even with a
"silent" failure when they have the same size).

Personally, I always feel uncomfortable with routines that "think" for
the user (same for programs, e.g. MSWord...), so I would entirely remove
that "auto-meshgridding" thing. Does it provide any performance or
memory gain?

Here's a test for non-gridded XI and YI based on my initial 2D polar to
rectangular interpolation example (simplify if necessary):
-------------------------------------------------
angles=0:5:45;
radius=3:10;
[meshang,meshrad]=meshgrid(angles,radius);
funcz=cos(0.2*(meshang+meshrad));
xmin=min(min(meshrad.*cos(meshang*pi/180)));
xmax=max(max(meshrad.*cos(meshang*pi/180)));
ymin=min(min(meshrad.*sin(meshang*pi/180)));
ymax=max(max(meshrad.*sin(meshang*pi/180)));
vecx=xmin:xmax;
vecy=ymin:ymax;
[meshvx,meshvy]=meshgrid(vecx,vecy);
vecang=atan2(meshvy,meshvx)*180/pi;
vecrad=sqrt(meshvx.^2+meshvy.^2);
vecz=interp2(meshang,meshrad,funcz,vecang,vecrad);
vecz(isnan(vecz))=0;
ref=[0 0.80973 0.67773 0.51871 0.33902 0.1458 -0.05322 -0.25012;\
    0 -0.41921 -0.81121 -0.91341 -0.93802 -0.87589 -0.87448 -0.90443;\
    0 0.47898 0.90645 0.53622 0.17884 -0.10336 -0.23031 -0.31921;\
    0 -0.88756 -0.32277 0.54954 0.87754 0.86503 0.85058 0.67038;\
    0 0 -0.78262 -0.76296 -0.22033 0.31549 0.66296 0.81102;\
    0 0 0 -0.59707 -0.93724 -0.77413 -0.316 0;\
    0 0 0 0 -0.35388 -0.88415 0 0;\
    0 0 0 0 0 -0.08131 0 0];
assert(max(max(abs(ref-vecz))) < 1.e-5);
-------------------------------------------------

Thank you for your help and explanations,
-Pierre.




-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------