Principal eigenvectors, eigenvalues and cross product of diffusion ellipsoid successfully calculated.

A = size(MT2);
m = A(1);
n = A(2);
p = A(3);
q = 9;
s=3;
matr=zeros(m,n,p,9);
MT3eval = zeros(m,n,p,3);
MT3evec = zeros(m,n,p,9);

%We're going to manually alter the small number of points which have
%negative eigenvalues (~1% of points). We do this by replacing them by the mean of their
%neighbours. If the program gets stuck at this point, then a point has a
%negative eigenvalue AND the mean of it's neighbours also has a negative
%eigenvalue. This point may have to be manually changed.
counter = 2;
%while counter > 1
    %counter = 1
    matr = [MT3(ii,jj,kk,1),MT3(ii,jj,kk,2),MT3(ii,jj,kk,3); MT3(ii,jj,kk,2),MT3(ii,jj,kk,4),MT3(ii,jj,kk,5);MT3(ii,jj,kk,3),MT3(ii,jj,kk,5),MT3(ii,jj,kk,6)];
    for ii = 1:m
        for jj = 1:n
            for kk = 1:p
                for ll=1:q
               
                [u v] = eig(matr);
                %evec(ll) = v(ll)
                MT3evec(ii,jj,kk,ll)=v(ll);
                MT3eval(ii,jj,kk,1)=u(1);
                MT3eval(ii,jj,kk,2)=u(5);
                MT3eval(ii,jj,kk,3)=u(9);
               
        end
            end
        end
       
end
to get the eval matrix and evec matrix

then

%%Cross product of diffusion tensor principle eigenvectors
% use variables from runnin 'diffusion cross' as all scans are the same
% size
Cross13=zeros(m,n,p,9);
Cross23=zeros(m,n,p,9);

Principle1_1=zeros(m,n,p,3);
Principle1_2=zeros(m,n,p,3);
Principle1_3=zeros(m,n,p,3);
Principle3_1=zeros(m,n,p,3);
Principle3_2=zeros(m,n,p,3);
Principle3_3=zeros(m,n,p,3);

Principle1_1(:,:,:,1)=MT1evec(:,:,:,1).*MT1eval(:,:,:,1);
Principle1_2(:,:,:,1)=MT1evec(:,:,:,2).*MT1eval(:,:,:,2);
Principle1_3(:,:,:,1)=MT1evec(:,:,:,3).*MT1eval(:,:,:,3);
Principle1_1(:,:,:,2)=MT1evec(:,:,:,4).*MT1eval(:,:,:,1);
Principle1_2(:,:,:,2)=MT1evec(:,:,:,5).*MT1eval(:,:,:,2);
Principle1_3(:,:,:,2)=MT1evec(:,:,:,6).*MT1eval(:,:,:,3);
Principle1_1(:,:,:,3)=MT1evec(:,:,:,7).*MT1eval(:,:,:,1);
Principle1_2(:,:,:,3)=MT1evec(:,:,:,8).*MT1eval(:,:,:,2);
Principle1_3(:,:,:,3)=MT1evec(:,:,:,9).*MT1eval(:,:,:,3);

Principle2_1(:,:,:,1)=MT2evec(:,:,:,1).*MT2eval(:,:,:,1);
Principle2_2(:,:,:,1)=MT2evec(:,:,:,2).*MT2eval(:,:,:,2);
Principle2_3(:,:,:,1)=MT2evec(:,:,:,3).*MT2eval(:,:,:,3);
Principle2_1(:,:,:,2)=MT2evec(:,:,:,4).*MT2eval(:,:,:,1);
Principle2_2(:,:,:,2)=MT2evec(:,:,:,5).*MT2eval(:,:,:,2);
Principle2_3(:,:,:,2)=MT2evec(:,:,:,6).*MT2eval(:,:,:,3);
Principle2_1(:,:,:,3)=MT2evec(:,:,:,7).*MT2eval(:,:,:,1);
Principle2_2(:,:,:,3)=MT2evec(:,:,:,8).*MT2eval(:,:,:,2);
Principle2_3(:,:,:,3)=MT2evec(:,:,:,9).*MT2eval(:,:,:,3);

Principle3_1(:,:,:,1)=MT3evec(:,:,:,1).*MT3eval(:,:,:,1);
Principle3_1(:,:,:,1)=MT3evec(:,:,:,2).*MT3eval(:,:,:,2);
Principle3_1(:,:,:,1)=MT3evec(:,:,:,3).*MT3eval(:,:,:,3);
Principle3_1(:,:,:,2)=MT3evec(:,:,:,4).*MT3eval(:,:,:,1);
Principle3_1(:,:,:,2)=MT3evec(:,:,:,5).*MT3eval(:,:,:,2);
Principle3_1(:,:,:,2)=MT3evec(:,:,:,6).*MT3eval(:,:,:,3);
Principle3_1(:,:,:,3)=MT3evec(:,:,:,7).*MT3eval(:,:,:,1);
Principle3_1(:,:,:,3)=MT3evec(:,:,:,8).*MT3eval(:,:,:,2);
Principle3_1(:,:,:,3)=MT3evec(:,:,:,9).*MT3eval(:,:,:,3);

cross131=cross(Principle1_1,Principle3_1,4);
cross231=cross(Principle1_1,Principle3_1,4);
cross132=cross(Principle1_2,Principle3_2,4);
cross232=cross(Principle1_2,Principle3_2,4);
cross133=cross(Principle1_3,Principle3_3,4);
cross233=cross(Principle1_3,Principle3_3,4);

Cross13(:,:,:,1)=cross131(:,:,:,1);
Cross13(:,:,:,2)=cross132(:,:,:,1);
Cross13(:,:,:,3)=cross133(:,:,:,1);
Cross13(:,:,:,4)=cross131(:,:,:,2);
Cross13(:,:,:,5)=cross132(:,:,:,2);
Cross13(:,:,:,6)=cross133(:,:,:,2);
Cross13(:,:,:,7)=cross131(:,:,:,3);
Cross13(:,:,:,8)=cross132(:,:,:,3);
Cross13(:,:,:,9)=cross133(:,:,:,3);

Cross23(:,:,:,1)=cross231(:,:,:,1);
Cross23(:,:,:,2)=cross232(:,:,:,1);
Cross23(:,:,:,3)=cross233(:,:,:,1);
Cross23(:,:,:,4)=cross231(:,:,:,2);
Cross23(:,:,:,5)=cross232(:,:,:,2);
Cross23(:,:,:,6)=cross233(:,:,:,2);
Cross23(:,:,:,7)=cross231(:,:,:,3);
Cross23(:,:,:,8)=cross232(:,:,:,3);
Cross23(:,:,:,9)=cross233(:,:,:,3);

to form the cross products of each of the three vectors forming the principle axes of the diffusion tensor ellipsoid.  Long hand and ugly script, but it works (I think)

Now to visualise them, perhaps a 3D quiver plot? 

Slices is a bit cumbersome in this perspective, unless the ellipsoid is projected and just the x and y values used? 

Will think on best methods of visualising and work on 7T data set, which is still 'holey'

Principle eigenvectors and cross product of mouse difusion tensors.

Previous post had a wrong method, the values used were not the principle eigenvectors. 

Diffusion matrix is upper symmetric,
[Dxx Dxy Dzy
  0      Dyy Dyz
 0       0      Dzz]

=

[Dxx Dxy Dzy
 Dyx Dyy Dyz
 Dzx Dzy Dzz]
             

so when w take the eigenvectors and values in MATLAB with

matr[u v]=eig(D)

we get a 3x3 matrix of eigenvalues u
and a 3x3 matrix of eigenvectors v

u is a diagonal matrix and the elements are the eigenvalues

v is a matrix of 3, 3x1 column vectors which are the principle eigenvectors of D

To take the cross product, we should take the cross product of each principle eigenvector with the corresponding eigenvector extracted from the diffusion matrix at the next time point.

So to get the matrix of the cross products, we will need to take the cross product three times for each pair and form this into a new matrix, unless there is a quicker command.

Mouse Tensor data

Extracted the diffusion tensor successfully from the mouse data today using MrTrix3

Commands:

navigate to folder with scans in, used APQ_0299

cd scan1
ls
## In each folder scan1, scan2, scan3, there is a file EPIDW_N4.nii.gz which is the diffusion weighted scan in a .nii.gz format, the mask has been created by Nyoman and is there as mask.nii the gradient file is the same for all the scans and is encoding.txt which is located in the parent folder of all the scans, two folders back.

inside folder /scan1:
dwi2tensor EPIDW_N4.nii.gz EPItensor1.nii -mask mask.nii -grad ../../encoding.txt
will create the tensor file as a .nii

Repeat this step inside scan2 folder and scan3 folder, changing the output filename as desired.  I used EPItensor, EPItensor2, EPItensor3.

Check in mrview to see if there are no weird islands, artifacts or missing bits.

Then load in MATLAB using:

MouseTens1 = load_untouch_nii('EPItensor1.nii')

which will create the structure MouseTens1 with substructures hdr (header info) and variables ext (extension info), untouch, and img.  img is the matrix with the actual diffusion data/

MouseTens1.hdr.dime.dim gives the size of the array.  4 (time), 64,48,76 (x,y,z), 6 (diffusion matrix values), 1,1,1


MT1=MouseTens1.img;

will create a matrix size 64 x 48 x 76 x 6 which is our diffusion matrix

we should then be able to compare these matrices directly.


Principle eigenvectors:

The 4th dimension of the matrix MT1,MT2, MT3 contains the six components of the diagonal diffusion matrix.  Components 1, 4, 6 are Dxx, Dyy, Dzz, the magnitudes of the principle eigenvectors in the x,y,z directions. 

We extract these into three new matrixes:
MT1ev=MT1(:,:,:,[1 4 6])
etc.
we can then take the cross product via
Cross12=cross(MT1ev,MT2ev);
etc
and compare the cross product matrices. 
What is the best numerical value to use as a marker of change?