Re: why cell arrays are used in 'interp2'

[Ben Abbott]

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