PCA in 3D for 3 major axes

Question

I have an anatomical 3D image (B), which is an indexed image i, j, k:

B(1,1,1)=0 %for example. 

The file contains only 0 and 1 sec.

I can visualize it correctly with isosurface: isosurface(B);

I would like to cut off the volume at some coordinate in j (it differs for each volume).

The problem is that the volume is tilted vertically, perhaps it is 45%, so the cut will not correspond to the anatomical volume.

I would like to get a new orthogonal coordinate system for the data, so my plane in the j coordinate will more accurately cut off the anatomical volume.

I was told to do this using the PCA, but I don't know how to do it, and reading the help pages did not help. Any direction would be welcome!

EDIT: I follow the recommendations, and now I have a new volume with a zero position, but I think that the axes do not follow the anatomical image as they should. These are preliminary and postal images:

This is the code I used to create the images (I took the code from the answer and the idea from the comments):

 clear all; close all; clc; hippo3d = MRIread('lh_hippo.mgz'); vol = hippo3d.vol; [IJK] = size(vol); figure; isosurface(vol); % customize view and color-mapping of original volume daspect([1,1,1]) axis tight vis3d; camlight; lighting gouraud camproj perspective colormap(flipud(jet(16))); caxis([0 1]); colorbar xlabel x; ylabel y; zlabel z box on % create the 2D data matrix a = 0; for i=1:K for j=1:J for k=1:I a = a + 1; x(a) = i; y(a) = j; z(a) = k; v(a) = vol(k, j, i); end end end [COEFF SCORE] = princomp([x; y; z; v]'); % check that we get exactly the same image when going back figure; atzera = reshape(v, I, J, K); isosurface(atzera); % customize view and color-mapping for the check image daspect([1,1,1]) axis tight vis3d; camlight; lighting gouraud camproj perspective colormap(flipud(jet(16))); caxis([0 1]); colorbar xlabel x; ylabel y; zlabel z box on % Convert all columns from SCORE xx = reshape(SCORE(:,1), I, J, K); yy = reshape(SCORE(:,2), I, J, K); zz = reshape(SCORE(:,3), I, J, K); vv = reshape(SCORE(:,4), I, J, K); % prepare figure %clf figure; set(gcf, 'Renderer','zbuffer') % render isosurface at level=0.5 as a wire-frame isoval = 0.5; [pf,pv] = isosurface(xx, yy, zz, vv, isoval); p = patch('Faces',pf, 'Vertices',pv, 'FaceColor','none', 'EdgeColor',[0.5 1 0.5]); % customize view and color-mapping daspect([1,1,1]) axis tight vis3d;view(-45,35); camlight; lighting gouraud camproj perspective colormap(flipud(jet(16))); caxis([0 1]); colorbar xlabel x; ylabel y; zlabel z box on 

Can someone give a clue what might happen? It seems that the problem might be the reshape command, is it possible that I am canceling the previously done work?

Answer 1

Not sure about the PCA, but here is an example showing how to visualize 3D scalar volume data, and volume reduction on an inclined plane (without axis aligned). The code is inspired by this demo in the MATLAB documentation.

 % volume data [x,y,z,v] = flow(); vv = double(v < -3.2); % threshold to get volume with 0/1 values vv = smooth3(vv); % smooth data to get nicer visualization :) xmn = min(x(:)); xmx = max(x(:)); ymn = min(y(:)); ymx = max(y(:)); zmn = min(z(:)); zmx = max(z(:)); % let create a slicing plane at an angle=45 about x-axis, % get its coordinates, then immediately delete it n = 50; h = surface(linspace(xmn,xmx,n), linspace(ymn,ymx,n), zeros(n)); rotate(h, [-1 0 0], -45) xd = get(h, 'XData'); yd = get(h, 'YData'); zd = get(h, 'ZData'); delete(h) % prepare figure clf set(gcf, 'Renderer','zbuffer') % render isosurface at level=0.5 as a wire-frame isoval = 0.5; [pf,pv] = isosurface(x, y, z, vv, isoval); p = patch('Faces',pf, 'Vertices',pv, ... 'FaceColor','none', 'EdgeColor',[0.5 1 0.5]); isonormals(x, y, z, vv, p) % draw a slice through the volume at the rotated plane we created hold on h = slice(x, y, z, vv, xd, yd, zd); set(h, 'FaceColor','interp', 'EdgeColor','none') % draw slices at axis planes h = slice(x, y, z, vv, xmx, [], []); set(h, 'FaceColor','interp', 'EdgeColor','none') h = slice(x, y, z, vv, [], ymx, []); set(h, 'FaceColor','interp', 'EdgeColor','none') h = slice(x, y, z, vv, [], [], zmn); set(h, 'FaceColor','interp', 'EdgeColor','none') % customize view and color-mapping daspect([1,1,1]) axis tight vis3d; view(-45,35); camlight; lighting gouraud camproj perspective colormap(flipud(jet(16))); caxis([0 1]); colorbar xlabel x; ylabel y; zlabel z box on 

Below is a result showing an isosurface made in the form of a wire frame, in addition to the cut planes oriented along the axes, and one rotated. We see that the volume space inside the isosurface has values ​​equal to 1, and 0 on the outside

Answer 2

I do not think PCA solves your problem. If you apply the PCA to your data, this will give you three new axes. These axes are called main components (PCs). They have the property that the first PC has the greatest dispersion when data is projected onto it. The second PC should also have the greatest dispersion when data is projected onto it, taking into account that it should be orthogonal to the first, the third PC is similar.

Now the question arises, when you project your data into a new coordinate system (defined by a PC), will it correspond to the anatomical volume? Not necessarily and most likely will not. The right axes for your data do not have a PCA optimization goal.

Answer 3

Sorry, I tried to answer @ Tevfik-Aytekin, but the answer is too long.

Hope this answer is helpful for someone:

Hi @Tevfik, thanks for your reply. I have been struggling for several days with the same problem, and I think I realized it right now.

I think the difference in what you say is that I work with coordinates. When I execute the PCA on my coordinates, I get a 3x3 transformation matrix for my data (the COEFF matrix, which is unitary and orthogonal, it's just a rotation matrix), so I know that I get exactly the same volume during the conversion.

These are the following steps:

• I had (I, J, K), 3D volume.
• According to @werner's suggestion, I changed it to 4 column matrix (x, y, z, v), size (I * J * K, 4).
• The coordinates (x, y, z) are excluded when v == 0 and v too. So now the size (initial volume, 3). Only coordinates with a value of 1, so the length of the vector is the volume of the figure.
• Run PCA to get COEFF and SCORE.
• COEFF - 3x3 matrix. It is unitary and orthogonal, it is a rotation matrix for my volume data.
• I did the editing on the PCA1 axis (i.e. it deletes al values ​​when COEFF (1) is greater than some value). That was my problem, I wanted to cut the volume perpendicular to the main axis.
• That was enough for me, the exact coordinates give me what I wanted. But:
• I didn’t have to come back, because I just need the remaining volume, but
• To return, I just had to restore the original coordinates. So first I converted the remaining coordinates using SCORE * COEFF '.
• Then I created the matrix again (I * J * K, 4), making v column = 1 only when the transformed coordinates coincided with the new matrix (with the message, option 'row').
• I created the indexed volume back using reshape (v, I, J, K).
• If I visualize the new volume back, it is cut perpendicular to the main geometric axes of the figure, as I needed.

Please, I would really like to hear any comments or suggestions on this method.