MATLAB Central Newsreader - Voxel loop optimisation

Voxel loop optimisation

Sven — Thu, 27 Mar 2008 02:13:02 -0400

Hi there, I've got a 3d matrix and I'm trying to select a
subset of elements to create a 2d matrix. At the moment I'm
doing this in a loop, and was wondering if anyone had a
trick they could teach me to avoid this loop becoming a
bottleneck.

Imagine a 3x3x3 rubix cube of elements:
>> rc = rand(3,3,3);

I can select, say, a plane through the 3rd dimension as follows:

>> myPlanar2d = squeeze(rc(2,:,:));

But what if I want to move along the 3rd dimension of my
rubix cube and select, say, the [1 3 2] rows instead of just
the second row as above? At the moment I'm doing this in a
loop as follows:

idxs = [1 3 2];
myIndexed2d = zeros(3,3);
for i = 1:length(idxs)
myIndexed2d(i,:) = rc(idxs(i),:,i);
end

Is there a shortcut I can take that avoids my loop?

Thanks,
Sven.

Re: Voxel loop optimisation

Phil Goddard — Thu, 27 Mar 2008 06:05:10 -0400

You could use ARRAYFUN in conjunction with an appropriate
function handle, i.e.

fh = @(idx,page) rc(idx,:,page)';
myIndexed2d = cell2mat(arrayfun(fh,idxs,1:length
(idxs),'UniformOutput',false))';

However, if you're looking for speed, then in newer
versions of MATLAB (at least since since R2006a) the loop
is likely to be faster for the purely numeric data that
you're manipulating.

Phil.

Re: Voxel loop optimisation

Peter Boettcher — Thu, 27 Mar 2008 14:39:07 -0400

"Sven " <sven.holcombe@gmail.deleteme.com> writes:

> Hi there, I've got a 3d matrix and I'm trying to select a
> subset of elements to create a 2d matrix. At the moment I'm
> doing this in a loop, and was wondering if anyone had a
> trick they could teach me to avoid this loop becoming a
> bottleneck.
>
> Imagine a 3x3x3 rubix cube of elements:
>>> rc = rand(3,3,3);
>
> I can select, say, a plane through the 3rd dimension as follows:
>
>>> myPlanar2d = squeeze(rc(2,:,:));
>
> But what if I want to move along the 3rd dimension of my
> rubix cube and select, say, the [1 3 2] rows instead of just
> the second row as above? At the moment I'm doing this in a
> loop as follows:
>
> idxs = [1 3 2];
> myIndexed2d = zeros(3,3);
> for i = 1:length(idxs)
> myIndexed2d(i,:) = rc(idxs(i),:,i);
> end
>
> Is there a shortcut I can take that avoids my loop?

You can use sub2ind. Build 3 3x3 matrices, where the first contains the
i indices desired for each of the elements, the second the j indices, etc.

lin_ind = sub2ind(size(rc), repmat(idxs.', 1, 3), repmat(1:3, 3, 1), repmat((1:3).', 1, 3));

myIndexed2d = rc(lin_ind);

I leave the timing to you. The loop will probably be faster, though
sub2ind might win if you need to do this repeatedly, and precompute one
or more of the index matrices.

meshgrid or ndgrid can also help generate the "easy" matrices. For
instance, the second two matrices above (j and k) can be generated like:

[j k] = ndgrid(1:3);

-Peter

Re: Voxel loop optimisation

Sven — Tue, 01 Apr 2008 00:35:04 -0400

"Phil Goddard" <philgoddardNOSPAM@telus.net> wrote in
message <fsfdem$ait$1@fred.mathworks.com>...
> You could use ARRAYFUN in conjunction with an appropriate
> function handle, i.e.
>
> fh = @(idx,page) rc(idx,:,page)';
> myIndexed2d = cell2mat(arrayfun(fh,idxs,1:length
> (idxs),'UniformOutput',false))';
>
> However, if you're looking for speed, then in newer
> versions of MATLAB (at least since since R2006a) the loop
> is likely to be faster for the purely numeric data that
> you're manipulating.
>
> Phil.

Thanks Phil and Peter.
The real voxel data I'm working with is much larger than
3x3x3, so the sub2ind method ends up taking too much memory
to build the index arrays. The arrayfun method is new to me
so I'll try this and check its speed compared to the loop.

Thanks again,
Sven.