Path: news.mathworks.com!newsfeed-00.mathworks.com!panix!bloom-beacon.mit.edu!llnews!53ab2750!not-for-mail
From: Peter Boettcher <boettcher@ll.mit.edu>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Voxel loop optimisation
References: <fsevre$cg7$1@fred.mathworks.com>
Message-ID: <muytzisjkok.fsf@G99-Boettcher.llan.ll.mit.edu>
Organization: MIT Lincoln Laboratory
User-Agent: Gnus/5.110006 (No Gnus v0.6) Emacs/23.0.0 (gnu/linux)
Cancel-Lock: sha1:SBWHfO8cua3NVDVQrnn1fehHvJE=
MIME-Version: 1.0
Content-Type: text/plain; charset=us-ascii
Lines: 46
Date: Thu, 27 Mar 2008 10:39:07 -0400
NNTP-Posting-Host: 155.34.163.114
X-Complaints-To: news@ll.mit.edu
X-Trace: llnews 1206628161 155.34.163.114 (Thu, 27 Mar 2008 10:29:21 EDT)
NNTP-Posting-Date: Thu, 27 Mar 2008 10:29:21 EDT
Xref: news.mathworks.com comp.soft-sys.matlab:459479


"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