[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: why cell arrays are used in 'interp2'
|
From: |
Sergei Steshenko |
|
Subject: |
Re: why cell arrays are used in 'interp2' |
|
Date: |
Tue, 25 Dec 2012 10:33:04 -0800 (PST) |
----- Original Message -----
> From: Ben Abbott <address@hidden>
> To: Sergei Steshenko <address@hidden>
> Cc: Octave users list <address@hidden>
> Sent: Tuesday, December 25, 2012 7:59 PM
> Subject: Re: why cell arrays are used in 'interp2'
>
> On Dec 25, 2012, at 10:52 AM, Sergei Steshenko wrote:
>> ----- Original Message -----
>>> From: Ben Abbott <address@hidden>
>>> To: Sergei Steshenko <address@hidden>
>>> Cc: Octave users list <address@hidden>
>>> Sent: Tuesday, December 25, 2012 5:27 PM
>>> Subject: Re: why cell arrays are used in 'interp2'
>>>
>>> On Dec 24, 2012, at 10:25 PM, Sergei Steshenko wrote:
>>>
>>>> Hello,
>>>>
>>>> looking into
> 'octave-3.6.2/share/octave/3.6.2/m/general/interp2.m'
>>> file I see:
>>>>
>>>> "
>>>> 265 ## construct the cubic hermite base functions in x,
> y
>>>> 266
>>>> 267 ## formulas:
>>>> 268 ## b{1,1} = ( 2*t.^3 - 3*t.^2 + 1);
>>>> 269 ## b{2,1} = h.*( t.^3 - 2*t.^2 + t );
>>>> 270 ## b{1,2} = (-2*t.^3 + 3*t.^2 );
>>>> 271 ## b{2,2} = h.*( t.^3 - t.^2 );
>>>> 272
>>>> 273 ## optimized equivalents of the above:
>>>> 274 t1 = tx.^2;
>>>> 275 t2 = tx.*t1 - t1;
>>>> 276 xb{2,2} = hx.*t2;
>>>> 277 t1 = t2 - t1;
>>>> 278 xb{2,1} = hx.*(t1 + tx);
>>>> 279 t2 += t1;
>>>> 280 xb{1,2} = -t2;
>>>> 281 xb{1,1} = t2 + 1;
>>>> 282
>>>> 283 t1 = ty.^2;
>>>> 284 t2 = ty.*t1 - t1;
>>>> 285 yb{2,2} = hy.*t2;
>>>> 286 t1 = t2 - t1;
>>>> 287 yb{2,1} = hy.*(t1 + ty);
>>>> 288 t2 += t1;
>>>> 289 yb{1,2} = -t2;
>>>> 290 yb{1,1} = t2 + 1;
>>>> 291
>>>> 292 ZI = zeros (size (XI));
>>>> 293 for i = 1:2
>>>> 294 for j = 1:2
>>>> 295 zidx = sub2ind (size (Z), yidx+(j-1),
> xidx+(i-1));
>>>> 296 ZI += xb{1,i} .* yb{1,j} .* Z(zidx);
>>>> 297 ZI += xb{2,i} .* yb{1,j} .* DX(zidx);
>>>> 298 ZI += xb{1,i} .* yb{2,j} .* DY(zidx);
>>>> 299 ZI += xb{2,i} .* yb{2,j} .* DXY(zidx);
>>>> 300 endfor
>>>> 301 endfor
>>>>
>>>> ".
>>>>
>>>> It looks to me only numeric data is used as array elements, so why
> cell
>>> arrays and not regular 2d matrices ?
>>>>
>>>> Speed ? Space ? Something else I couldn't think of ?
>>>>
>>>>
>>>> Thanks,
>>>> Sergei.
>>>
>>> The cell-arrays contian 2D matrices whose sizes are [numel(yi),
> numel(xi)].
>>>
>>> Ben
>>
>> 4d matirces are also supported by Octave directly through m(t, u, v, x)
> notation, so are cell arrays to improve memory utilization ? Or speed ?
>>
>> Thanks,
>> Sergei.
>
> I don't know if there is speed advantage, but if the 4D array needs to be
> "squeezed" before the element-wise multiplication (.*), then the cell
> approach may be faster.
>
> If you (someone else?) would like to check, the xb and yb arrays can be
> changed
> to 4D arrays whose size is [numel(yi), numel(xi), 2, 2].
>
> Then the cell array indices "{m,n}" would be replaced by
> "(:,:,,m,n)". This will eliminate the need of the user to squeeze()
> the data before multiplying, but I don't know what is actually done in the
> c++ code.
>
> If you do ckeck and find it is faster, I'd be happy to push a changeset for
> you.
>
> Ben
>
I haven't even tried the function yet - decided to look into code in order to
verify whether I correctly understand the documentation.
When I saw the code I decided to ask - my (apparently wrong) assumption was
that Octave developers know much better than I and can justify design decisions
made by them.
The documentation, by the way, doesn't have examples - please see below.
Surprisingly, Matlab documentation _does_ have examples - see
http://matlab.izmiran.ru/help/techdoc/ref/interp2.html .
No, I won't offer a patch to implement what's already done by others. I don't
see it as pragmatic.
Regards,
Sergei.
"
octave:1> help interp2
`interp2' is a function from the file
/mnt/sdb8/sergei/AFSWD_debug/20121021/octave-3.6.2/share/octave/3.6.2/m/general/interp2.m
-- Function File: ZI = interp2 (X, Y, Z, XI, YI)
-- Function File: ZI = interp2 (Z, XI, YI)
-- Function File: ZI = interp2 (Z, N)
-- Function File: ZI = interp2 (..., METHOD)
-- Function File: ZI = interp2 (..., METHOD, EXTRAPVAL)
Two-dimensional interpolation. X, Y and Z describe a surface
function. If X and Y are vectors their length must correspondent
to the size of Z. X and Y must be monotonic. If they are
matrices they must have the `meshgrid' format.
`interp2 (X, Y, Z, XI, YI, ...)'
Returns a matrix corresponding to the points described by the
matrices XI, YI.
If the last argument is a string, the interpolation method can
be specified. The method can be 'linear', 'nearest' or
'cubic'. If it is omitted 'linear' interpolation is assumed.
`interp2 (Z, XI, YI)'
Assumes `X = 1:rows (Z)' and `Y = 1:columns (Z)'
`interp2 (Z, N)'
Interleaves the matrix Z n-times. If N is omitted a value of
`N = 1' is assumed.
The variable METHOD defines the method to use for the
interpolation. It can take one of the following values
'nearest'
Return the nearest neighbor.
'linear'
Linear interpolation from nearest neighbors.
'pchip'
Piecewise cubic Hermite interpolating polynomial.
'cubic'
Cubic interpolation from four nearest neighbors.
'spline'
Cubic spline interpolation--smooth first and second
derivatives throughout the curve.
If a scalar value EXTRAPVAL is defined as the final value, then
values outside the mesh as set to this value. Note that in this
case METHOD must be defined as well. If EXTRAPVAL is not defined
then NA is assumed.
See also: interp1
Additional help for built-in functions and operators is
available in the on-line version of the manual. Use the command
`doc <topic>' to search the manual index.
Help and information about Octave is also available on the WWW
at http://www.octave.org and via the address@hidden
mailing list.
octave:2>
"