PCA on a 3D image to obtain 3 principal axes

Question

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

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

The file contains only 0s and 1s.

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

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

The problem is that the volume is tilted vertically, it maybe has 45% degrees, so the cut will not be following the anatomical volume.

I would like to obtain a new orthogonal coordinate system for the data, so my plane in coordinate j would cut the anatomical volume in a more accurate way.

I've been told to do it with PCA, but I don't have a clue how to do it, and reading the help pages haven't been of help. Any direction will be welcome!

EDIT: I have been following the recommendations, and now I got a new volume, zero-centered, but I think that axes don't follow the anatomical image as they should. These are the pre and post images:

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

clear all; close all; clc;
hippo3d = MRIread('lh_hippo.mgz');
vol = hippo3d.vol;

[I J K] = 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 anybody provide a hint what might be happening? It seems that the problem might be the reshape command, Is it possible that I am canceling out the job previously done?

Answer 1

Not sure about PCA, but here is an example showing how to visualize a 3D scalar volume data, and cutting the volume at a tilted plane (non-axis aligned). 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 the result showing the isosurface rendered as wire-frame, in addition to slicing planes both axes-aligned and one rotated. We can see that the volume space on the inside of the isosurface has values equal to 1, and 0 on the outside

Answer 2

I don't think that PCA solves your problem. If you apply PCA to your data it will give you three new axes. These axes are called principal components (PCs). They have the property that the first PC has the largest variance when the data is projected on it. The second PC must also has the largest variance when data is projected on it subject to the constraint that it should be orthogonal to the first, the third PC is similar.

Now the question is when you project your data into the new coordinate system (defined by the PCs) will it match the anatomical volume? Not necessarily and most probably will not. The right axes for your data do not have the optimization objective of PCA.

Answer 3

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

Hopefully this answer will be useful for somebody:

Hi @Tevfik, thanks for your answer. I've struggling for days with this same problem, and I think I got it right now.

I think that the difference respect to what you are saying is that I am working with coordinates. When I perform PCA over my coordinates, I get a 3x3 transformation matrix for my data (COEFF matrix, which is unitary and orthogonal, it is just a rotation matrix), so I know that I get exactly the same volume when transformed.

These are the steps I followed:

• I had a (I,J,K), 3D volume.
• As per @werner suggestion, I changed it to a 4 column matrix (x,y,z,v), size (I*J*K, 4).
• Eliminated the coordinates (x,y,z) when v == 0, and v too. So right now, size is (original volume, 3). Only the coordinates with value 1, so the length of the vector is the volume of the figure.
• Perform PCA to obtain COEFF and SCORE.
• COEFF is a 3x3 matrix. It is unitary and orthogonal, it is a rotation matrix for my volume data.
• I did the editing in the PCA1 axis (i.e. delete al values when COEFF(1) is bigger than some-value). This was my problem, I wanted to cut the volume perpendicular to the main axis.
• This was enough for me, the reamining coordinates are giving me the volume I wanted. But:
• I didn't need to go back, as I just needed the remaining volume, but
• In order to go back, I just had to reconstruct the original coordinates. So first I transformed the remaining coordinates with SCORE*COEFF'.
• Then I created again a (I*J*K, 4) matrix, making the v column = 1 only when the transformed coordinates matched the new matrix (with ismember, 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, just as I needed.

Please, I would really like to hear any comment or suggestion on the method.