Three-dimensional projection of a three-dimensional image on a two-dimensional plane

Question

I have a three-dimensional shade of gray, corresponding to ultrasound data. In Matlab, this three-dimensional volume is simply a three-dimensional matrix of MxNxP. The structure that interests me is not oriented along the z axis, but along the already known coordinate system (x'y'z '). What I have so far is similar to the following figure depicting the source (xyz) and local coordinate systems (x'y'z '):

I want to get a two-dimensional projection of this volume (that is, an image) through a specific plane in 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 easily achieved. those. if the volume in Matlab is equal to V, then:

projection = sum(V,3); 

Thus, the projection can be calculated as the sum of the 3rd size of the array. However, with a change in orientation, the problem becomes more complicated.

I looked at the radon transform (2D, which applies only to two-dimensional images, not to volumes), and also considers orthographic projections, but at the moment I do not know what to do!

Thanks for any advice!

Answer 1

New solution attempt:

After the tutorial http://blogs.mathworks.com/steve/2006/08/17/spatial-transformations-three-dimensional-rotation/ and making minor changes, I might have something that could help you. Keep in mind that I have little or no experience with bulk data in MATLAB, so the implementation is pretty hacky.

In the code below, I use tformarray () to rotate the structure in space. First, the data is centered and then rotated using rotationmat3D to create a spatial transformation before the data is returned to its original position.

Like I never used tformarray before, I carried datapoints that went beyond a certain area after rotation, just filling the data matrix (NxMxP) with zeros around. If someone knows a better way, let us know :)

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 a confirmation that the data rotates as expected and builds a closed surface when the data passes the value “1”. With blob2 you should know that you can design using simple sums.

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

Answer 2

Assuming that you have access to the coordinate base R = [x 'y' z '], and that these vectors are orthonormal, you can simply extract the representation in this basis by multiplying your data by a 3x3 R matrix, where x', y ' , z 'are column vectors.

With the data stored in D (Nx3), you can get a view with R by multiplying it by: Dmarked = D * R;

and now D = Dmarked * inv (R), so moving back and forth is unstable.

The following code may help to see the conversion. Here I create a synthetic dataset, rotate it, and then rotate it back. Fulfilling the sum (DR (:, 3)) will then be your sum over z '

 %#Create synthetic dataset N1 = 250; r1 = 1; dr1 = 0.1; dz1 = 0; mu1 = [0;0]; Sigma1 = eye(2); theta1 = 0 + (2*pi).*rand(N1,1); rRand1 = normrnd(r1,dr1,1,N1); rZ1 = rand(N1,1)*dz1+1; D = [([rZ1*0 rZ1*0] + repmat(rRand1',1,2)).*[sin(theta1) cos(theta1)] rZ1]; %Create roation matrix rot = pi/8; R = rotationmat3D(rot,[0 1 0]); % R = 0.9239 0 0.3827 % 0 1.0000 0 % -0.3827 0 0.9239 Rinv = inv(R); %Rotate data DR = D*R; %#Visaulize data f1 = figure(1);clf subplot(1,3,1); plot3(DR(:,1),DR(:,2),DR(:,3),'.');title('Your data') subplot(1,3,2); plot3(DR*Rinv(:,1),DR*Rinv(:,2),DR*Rinv(:,3),'.r'); view([0.5 0.5 0.2]);title('Representation using your [xmarked ymarked zmarked]'); subplot(1,3,3); plot3(D(:,1),D(:,2),D(:,3),'.'); view([0.5 0.5 0.2]);title('Original data before rotation'); 

Answer 3

If you have two normalized 3x1 vectors x2 and y2 corresponding to your local coordinate system (x 'and y').

Then, for the position P its local coordinate will be xP=P'x2 and yP=P'*y2 .

So, you can try to design your volume using the drive:

 [xyz]=ndgrid(1:M,1:N,1:P); posP=[x(:) y(:) z(:)]; xP=round(posP*x2); yP=round(posP*y2); xP=xP+min(xP(:))+1; yP=yP+min(yP(:))+1; V2=accumarray([xP(:),yP(:)],V(:)); 

If you provide your details, I will check it.