Sun 04 Jul 2010 10:30:37 AM UTC, comment #3:
Sorry about the late reply (the last half year has been crazy for me...). To be honest I am not sure if I am the author of the code in question (I might be; I just don't remember), but I'll do a review anyway :-)
- ## If X and Y vectors produce a grid from them
+ ## Check dimensions of X and Y
if (isvector (X) && isvector (Y))
X = X(:).';
Y = Y(:);
Seems like a reasonable change to me (the new comment is more descriptive).
- ## If Xi and Yi are vectors of different orientation build a grid
- if ((rows (XI) == 1 && columns (YI) == 1)
- || (columns (XI) == 1 && rows (YI) == 1))
- ## Do nothing
+ ## Check dimensions of XI and YI
+ if (isvector (XI) && isvector (YI))
+ XI = XI(:).';
+ YI = YI(:);
This seems simpler then the current code. I am not sure if the
+ XI = XI(:).';
+ YI = YI(:);
part is really needed, but I guess it is better to play it safe, so this seems fine to me.
function b = isgriddata (X)
d1 = diff (X, 1, 1);
- d2 = diff (X, 1, 2);
- b = all (d1 (:) == 0) & all (d2 (:) == d2 (1));
+ b = all (d1 (:) == 0);
endfunction
This removes the check that the derivative of X is constant. I guess this is not a requirement for X to be gridded (consider e.g. the output of 'meshgrid (logspace (1, 10, 10))'). So, I think this change is fine as well.
In summary, I'd say check in the change.
|
Mon 19 Apr 2010 08:52:02 PM UTC, comment #1:
The attached patch fixes the fails of demo("interp2"). But now the cubic interpolation code within interp2 will never be used (as far as I can see). Could you have a look Søren (as you seem to be author of that code).
--- a/scripts/general/interp2.m Sun Apr 18 21:27:53 2010 +0200
+++ b/scripts/general/interp2.m Mon Apr 19 22:41:16 2010 +0200
@@ -318,7 +318,7 @@
else
- ## If X and Y vectors produce a grid from them
+ ## Check dimensions of X and Y
if (isvector (X) && isvector (Y))
X = X(:).';
Y = Y(:);
@@ -332,10 +332,10 @@
endif
endif
- ## If Xi and Yi are vectors of different orientation build a grid
- if ((rows (XI) == 1 && columns (YI) == 1)
- || (columns (XI) == 1 && rows (YI) == 1))
- ## Do nothing
+ ## Check dimensions of XI and YI
+ if (isvector (XI) && isvector (YI))
+ XI = XI(:).';
+ YI = YI(:);
elseif (! size_equal (XI, YI))
error ("XI and YI must be matrices of same size");
endif
@@ -409,8 +409,7 @@
function b = isgriddata (X)
d1 = diff (X, 1, 1);
- d2 = diff (X, 1, 2);
- b = all (d1 (:) == 0) & all (d2 (:) == d2 (1));
+ b = all (d1 (:) == 0);
endfunction
## Compute the bicubic interpolation coefficients
|
Sun 18 Apr 2010 09:03:51 PM UTC, original submission:
There seem to be several bugs in the cubic and spline interpolation code of interp2.
Reproduce with:
octave> demo interp2
[...]
interp2 example 7:
A=[13,-1,12;5,4,3;1,6,2];
x=[0,1,2]; y=[10,11,12];
xi=linspace(min(x),max(x),17);
yi=linspace(min(y),max(y),26)';
mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic'));
[x,y] = meshgrid(x,y);
hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
interp2 example 7: failed
mx_el_or: nonconformant arguments (op1 is 1x17, op2 is 26x1)
error: called from:
error: /home/thorsten/hg/octave/scripts/general/interp2.m at line 352, column 8
The offending line is
inside = !(XI < X (1) | XI > X (end) | YI < Y (1) | YI > Y (end));
To me it seems that this works only if xi and yi are meshgrid data as well. But this case is handled in the preceding conditional. Strange... Also strange: the bicubic function can handle the input of the above example, i.e.
mesh(xi,yi,bicubic(x,y,A,xi,yi));
works. So what is the purpose of the 2nd implementation of cubic interpolation in interp2.m?
next failure of the demos:
[...]
interp2 example 8:
[x,y,A] = peaks(10);
x = x(1,:)'; y = y(:,1);
xi=linspace(min(x),max(x),41);
yi=linspace(min(y),max(y),41)';
mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic'));
[x,y] = meshgrid(x,y);
hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
interp2 example 8: failed
interp2: input data must have `meshgrid' format
Here, the data generated by peaks are not accepted as meshgrid data due to a floating point problem. Also, the function isgriddata (defined in interp2.m) seems to require not only equal lines in the grid matrix but also equally spaced values. Is that intended (or really necessary?).
and the last failure in the demos:
interp2 example 9:
A=[13,-1,12;5,4,3;1,6,2];
x=[0,1,2]; y=[10,11,12];
xi=linspace(min(x),max(x),17);
yi=linspace(min(y),max(y),26)';
mesh(xi,yi,interp2(x,y,A,xi,yi,'spline'));
[x,y] = meshgrid(x,y);
hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
interp2 example 9: failed
interp2: input data must have `meshgrid' format
Looking at the sources, the _splinen_ function that actually calculates the interpolation is called only with the first row/column of the x and y variables, so why require meshgrids here?
regards
Thorsten
|