3d grayscale volume projection onto 2D plane

Question

I have a 3-D grayscale volume corresponding to ultrasound data. In Matlab this 3-D volume is simply a 3-D matrix of MxNxP. The structure I'm interested in is not oriented along the z axis, but along a local coordinate system already known (x'y'z'). What I have up to this point is something like the figure shown below, depicting the original (xyz) and the local coordinate systems (x'y'z'):

I want to obtain the 2-D projection of this volume (i.e. an image) through a specific plane on the local coordinate system, say at z' = z0. How can I do this?

If the volume was oriented along the z axis this projection could be readily achieved. i.e. if the volume, in Matlab, is V, then:

projection = sum(V,3);

thus, the projection can be computed just as the sum along the 3rd dimension of the array. However with a change of orientation the problem becomes more complicated.

I've been looking at radon transform (2D, that applies only to 2-D images and not volumes) and also been considering ortographic projections, but at this point I'm clueless as to what to do!

Thanks for any advice!

Answer 1

New attempt at solution:

Following the tutorial http://blogs.mathworks.com/steve/2006/08/17/spatial-transformations-three-dimensional-rotation/ and making some small changes, I might have something which could help you. Bear in mind, I have little or no experience with volumetric data in MATLAB, so the implementation is quite hacky.

In the below code I use tformarray() to rotate the structure in space. First, the data is centered, then rotated using rotationmat3D to produce the spacial transformation, before the data is moved back to its original position.

As I have never used tformarray before, I handeled datapoints falling outside the defined region after rotation by simply padding the data matrix (NxMxP) with zeros all around. If anyone know a better way, please let us know :)

The code:

    %Synthetic dataset, 25x50x25
blob = flow();

%Pad to allow for rotations in space. Bad solution, 
%something better might be possible to better understanding
%of tformarray()
blob = padarray(blob,size(blob));

f1 = figure(1);clf;
s1=subplot(1,2,1);
p = patch(isosurface(blob,1));
set(p, 'FaceColor', 'red', 'EdgeColor', 'none');
daspect([1 1 1]);
view([1 1 1])
camlight
lighting gouraud

%Calculate center
blob_center = (size(blob) + 1) / 2;

%Translate to origin transformation
T1 = [1 0 0 0
    0 1 0 0
    0 0 1 0
    -blob_center 1];

%Rotation around [0 0 1]
rot = -pi/3;
Rot = rotationmat3D(rot,[0 1 1]);
T2 = [ 1  0  0   0
       0  1  0   0
       0  0  1   0
       0  0  0   1];
T2(1:3,1:3) = Rot;   

%Translation back
T3 = [1 0 0 0
    0 1 0 0
    0 0 1 0
    blob_center 1];

%Total transform
T = T1 * T2 * T3;

%See http://blogs.mathworks.com/steve/2006/08/17/spatial-transformations-three-dimensional-rotation/
tform = maketform('affine', T);
R = makeresampler('linear', 'fill');
TDIMS_A = [1 2 3];
TDIMS_B = [1 2 3];
TSIZE_B = size(blob);
TMAP_B = [];
F = 0;
blob2 = ...
tformarray(blob, tform, R, TDIMS_A, TDIMS_B, TSIZE_B, TMAP_B, F);

s2=subplot(1,2,2);
p2 = patch(isosurface(blob2,1));
set(p2, 'FaceColor', 'red', 'EdgeColor', 'none');
daspect([1 1 1]);
view([1 1 1])
camlight
lighting gouraud

The arbitrary visualization below is just to confirm that the data is rotated as expected, plotting a closed surface when the data passed the value '1'. With blob2, you should know be able to project by using simple sums.

figure(2)
subplot(1,2,1);imagesc(sum(blob,3));
subplot(1,2,2);imagesc(sum(blob2,3));