Extract subset of time-curves from a 4D volume

Dear all,

I have a 4D volume and I can easily extract the temporal curve for each voxel.

I would like to extract only the curves for the voxels that correspond to a certain condition.

For example: V is the 4D volume and R is a 3D mask volume. I would like to extract only the curves of the points in which R>0 but without using a for cicle.

Normally I would write:

for k=1:size(V,3)
    for i=1:size(V,1)
        for j =1:size(V,2)
            if (R(i,j,k)>0)
                y = squeeze(V(i,j,k,:));
                a(i,:) = y;
                ......
        end
    end
end

I tried to use indexing but I can only use it with 3D volumes.

Is there an efficient solution (I mean fast) to do this?

Thanks a lot

[Turlough Hughes]

I'm not sure what you're data should look like in the end. Can you explain a bit more about the curves you are trying to produce. Are these supposed to be 2d plots?

[Domenico]

Yes, they are 1D functions.

The volume I'm speaking of is a 4D temporal volume in which each 3D volume represents a different time point.

[Turlough Hughes]

You could do something like this.

V2 = V; % make a copy of V
V2(repmat(R<=0,1,1,1,size(V2,4)))=NaN; % Remove values where R <= 0

My feeling is that you should maintain the shape of V, so I opted to insert NaN's at coordinates where your condition R>0 is false, ie where R<=0 is true. Is that the sort of answer you're looking for?

[Domenico]

Thanks a lot.

In this way you obtain a sort of "masking" of the 4D volume, right? From that masking, however, I would like to extract each single temporal 1D curve, i.e. for every survived pixel I would like to obtain its curve (without using for)

Answer 1

[Turlough Hughes]

I can reduce this to a single for loop which iterates through the timesteps but i'm not sure about doing this entirely without loop. Here, I concatentate the data of a given timestep into a column vector and insert it into results. I took the liberty of including the x,y and z indices in columns 1,2 and 3, respectively, so each timestep is added sequentially thereafter.

Rmask = R>0; % Binary mask
[Ix, Iy, Iz] = ind2sub(size(Rmask),find(Rmask)); % get x,y,z indices for mask
result = [Ix, Iy, Iz, zeros(numel(Ix),size(V,4))]; % preallocate space

for t = 1:size(V,4)
    Vtemp = V(:,:,:,t); % grab one timestep
    result(:,t+3) = Vtemp(Rmask); % Convert to column vector and insert into results
end

[Domenico]

Thanks a lot

I will try and let you know ASAP!

[Turlough Hughes]

Actually, its probably better to store intensities seperately, so that would be as follows if you wish to do so (I also just modified the name for result to being Vr):

Rmask = R>0; % Binary mask

% Get indices and store in Ixyz
[Ix, Iy, Iz] = ind2sub(size(Rmask),find(Rmask)); % get x,y,z indices for mask
Ixyz = [Ix Iy Iz];

Vr = [zeros(size(Ixyz,1),size(V,4))]; % preallocate space for V within Rmask
for t = 1:size(V,4)
    Vtemp = V(:,:,:,t); % grab one timestep
    Vr(:,t) = Vtemp(Rmask); % Convert to column vector and store
end

[Domenico]

Thanks a lot, it works!

I supposed the existence of a Matlab trick to avoid the for cicle, however this reduction works well.

[Turlough Hughes]

Actually, there is a nice way to do it and without for loops

Firstly get your indices.

[Ix, Iy, Iz] = ind2sub(size(Rmask),find(Rmask)); % get x,y,z indices for mask
Ixyz = [Ix Iy Iz];

Secondly I vectorize V and place it in Vr, then for the first time step you can get your corresponding logical index by R(:)>0, however for all timesteps you would need to use the repmat function to repeat the index according to the number of timesteps. Finally I reshape the data so it has the same dimensions as in previous examples, i.e. each column is a time step and each row is a position in x,y and z.

Vr = V(:);
Vr = Vr(repmat(R(:)>0,[size(V,4),1]));
Vr = reshape(Vr,nnz(R>0),size(V,4));

[Turlough Hughes]

I haven't tried but I expect this will be appreciably faster than my previous attempt. Also, when you are happy with an answer can you accept it.

[Domenico]

Ok, only a question: what do you do with Ixyz? You obtain it but you don't use it for Vr. Is right?

[Turlough Hughes]

Ixyz is the corresponding index in V. The first row of Vr shows values changing over time at a single point within V, the first row of Ixyz inidicates the corresponding xyz indices within V which can usually be related to an xyz position. I just assumed you would need it.