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? 

1D surgery case

One dimensional surgery simulation has finally been successfully coded, taking out a segment and leaving that as space with a boundary as the simulation, we also tried stitching the function bak together, so there are two cases, in small tumours, the appropriate case would be stiching together, as the brain biomechanics would allow that degree of deformation, but in larger resections, the former case is more biologically appropriate. 

Code:

%%1D-Fisher equation

%define constants
D = 2.5e-6;     %Mean Cancer-cell Diffusion (mm^2/s)
rho = 4e-7;         %Mean Cancer-cell growth rate (cells/second/voxel)
k = 2.4e5;          %Maximum number of cells per unit volume (cells/mm^3)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          



%% Define Domain

%The size of DFull is 128*128*60*6. We choose n, m and p so that at full
%resolution (r=1) the size of D is 130*130*62*6, and we add a single layer
%of empty cells around the brain. This is just some insurance, because if
%the brain is touching the boundary of the domain, errors can occur.
%based on DTI data, 2.5mm voxels

                               %xpoints
xrange = [2.5:2.5:130*2.5];                  %All these dimensions are in mm, and come straight from the MRI data
dx = 2.5;      %space step size
xmax=max(xrange);
%right boundary
%m = floor(130);                       %ypoints
%p = floor(62);                        %zpoints
tsize =200;                              %tsteps.  Unlike the space steps, this is variable. The stability condition requires it varies.
trange = [0, (2e8)];                   %The time dimension in seconds
tmax=max(trange);
dt = (trange(2)-trange(1))/tsize;          %time step size
tsteps = 1:dt:trange(2);
%stability
%% Initial Conditions

% Define Initial Condition as a function of x. Remember the units
% of c are cells per mm^3

% ic = @(x) 1*exp(-((x-mean(xrange))^2/1));
C=zeros(1,132); %initialise C
N = length(C);
C(62:67)=k/4; %I.C. cells here

CC = zeros(1,N); %post surgery function
%% Solve equation 1D FTCS scheme
  %Explicitly step through time with anisotropic diffusion
  %Equation dc/ct=D(\del^2)c + rho*c(1 -c/k), the Fisher KPP equation
  %becomes in FTCS scheme
  % c(i,k+1) = rc(i-1,k) + (1-2r)c(i,k) + rc(i+1,k)
  % where r = Ddt/(dx)^2
  %so set
  r=((D*dt)/(dx)^2)
   % Define C as a vector
    %dCdt = (-2*D/dx^2).*C;
 %C = C + dCdt*dt;
 %C = C + dt*rho.*abs(C).*(1-abs(C)/k);

 C(1) = 0;
C(132) = 0;
Cnew1 = zeros(1,N);
Cnew2 = zeros(1,N);
%C = C + dCdt*dt;
   

 figure
 hold off
 ConcSumUncut=zeros(1,tsize);

for t=1:tsize
    CC=C;
   
 


       for i = 3:(N-2)
             C(i) = r*CC(i+1) + (1-2*r)*CC(i) + r*CC(i-1);
       end
          
  
       %Also add a cell growth term - the Fisher Equation.
    %Note the absolute values are just back-ups to ensure stability. We
    %know that the explicit diffusion equation is numerically stable, as
    %small disturbances will die down. But with an extra growth term, small
    %disturbances may lead to negative numbers, which are then inflated.
    %This is just a safety-measure.
            C = C + dt*rho.*abs(C).*(1-abs(C)/k);
            plot(CC)  % ./k to normalise
            axis([0 132 0 k])
            pause(0.1)
            disp(t)
            %while t <= 30
                %continue
              

%   if t == 30  %then do surgery
                  
%is=find(C>0.2*k,1)
%for i=1:N
   % if C(i) < 0.2*k
        %CC(i) = C(i);
   % end
%end
%plot(N,CC)
%legend('Tumour', 'Resected tumour')


%for i = 1:(is-1)
     Cnew1(i+((N/2)-is)+1) = CC(i);
   
%end
%for i = (N-is):N
    Cnew2(i-(N/2-is)) = CC(i);

%end
%C = Cnew1 + Cnew2;



 %feed this back into main loop
    ConcSumUncut(t)=sum(C);
    end
          

            
 %plot(Csurg)
 

figure
plot(ConcSumUncut)
%end
%%Surgery function
%cuts C for C>is and stiches back together the ends

%%%One Dimensional tumour resection simulation.  Boundary at cut edges.
%%post-resection.
%%1D-Fisher equation

%define constants
D = 2.5e-6;     %Mean Cancer-cell Diffusion (mm^2/s)
rho = 4e-7;         %Mean Cancer-cell growth rate (cells/second/voxel)
k = 2.4e5;          %Maximum number of cells per unit volume (cells/mm^3)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          



%% Define Domain

%The size of DFull is 128*128*60*6. We choose n, m and p so that at full
%resolution (r=1) the size of D is 130*130*62*6, and we add a single layer
%of empty cells around the brain. This is just some insurance, because if
%the brain is touching the boundary of the domain, errors can occur.
%based on DTI data, 2.5mm voxels

                               %xpoints
xrange = [2.5:2.5:130*2.5];                  %All these dimensions are in mm, and come straight from the MRI data
dx = 2.5;      %space step size
xmax=max(xrange);
%right boundary
%m = floor(130);                       %ypoints
%p = floor(62);                        %zpoints
tsize =200;                              %tsteps.  Unlike the space steps, this is variable. The stability condition requires it varies.
trange = [0, (2e8)];                   %The time dimension in seconds
tmax=max(trange);
dt = (trange(2)-trange(1))/tsize;          %time step size
tsteps = 1:dt:trange(2);
%stability

%% Initial Conditions

% Define Initial Condition as a function of x. Remember the units
% of c are cells per mm^3

% ic = @(x) 1*exp(-((x-mean(xrange))^2/1));
C=zeros(1,132); %initialise C
N = length(C);
C(62:67)=k/4; %I.C. cells here

CC = zeros(1,N); %post surgery function

%% Solve equation 1D FTCS scheme
  %Explicitly step through time with anisotropic diffusion
  %Equation dc/ct=D(\del^2)c + rho*c(1 -c/k), the Fisher KPP equation
  %becomes in FTCS scheme
  % c(i,k+1) = rc(i-1,k) + (1-2r)c(i,k) + rc(i+1,k)
  % where r = Ddt/(dx)^2
  %so set
  r=((D*dt)/(dx)^2)
   % Define C as a vector
    %dCdt = (-2*D/dx^2).*C;
 %C = C + dCdt*dt;
 %C = C + dt*rho.*abs(C).*(1-abs(C)/k);

 C(1) = 0;
C(132) = 0;

   

 figure
 hold off
 ConcSum=zeros(1,tsize);

for t=1:30
    CC=C;
   
 


       for i = 2:(N-1)
             C(i) = r*CC(i+1) + (1-2*r)*CC(i) + r*CC(i-1);
       end
          
  
       %Also add a cell growth term - the Fisher Equation.
    %Note the absolute values are just back-ups to ensure stability. We
    %know that the explicit diffusion equation is numerically stable, as
    %small disturbances will die down. But with an extra growth term, small
    %disturbances may lead to negative numbers, which are then inflated.
    %This is just a safety-measure.
            C = C + dt*rho.*abs(C).*(1-abs(C)/k);
         % C = C + dt*0.001*rho.*abs(C).*(1-abs(C)/k).*(abs((C) - 0.2*k));
            plot(CC)  % ./k to normalise
            axis([0 132 0 k])
            pause(0.1)
            disp(t)
            %while t <= 30
                %continue
              

   if t == 30  %then do surgery
                  
CC = zeros(1,N);
is=find(C>0.2*k,1);
for i=1:N
    if C(i) < 0.2*k
        CC(i) = C(i);
    end
end
C = CC;
plot(CC)
%C(is+1) = C(is);
%C(N-is-1) = C(N-is);
 %feed this back into main loop
    
    end
          

ConcSum(t)=sum(C);
            
 %plot(Csurg)
 
end
               if  i == is-1;
             C(i) = (1-r)*CC(i) + r*CC(i-1);
               if   i ==(N-is+1);
               C(i) = r*CC(i+1) + (1-r)*CC(i);
               end
               end
              
for t=31:tsize
     CC=C;
   
 
 for i=2:(is-2)
       C(i) = r*CC(i+1) + (1-2*r)*CC(i) + r*CC(i-1);
          
         
 end

 for i=(N-is):N-1
     
       C(i) = r*CC(i+1) + (1-2*r)*CC(i) + r*CC(i-1);
     
      
 end
   %for i = (N-is+1):(N-1)
      % C(i) = r*CC(i+1) + (1-2*r)*CC(i) + r*CC(i-1);
  


       %Also add a cell growth term - the Fisher Equation.
    %Note the absolute values are just back-ups to ensure stability. We
    %know that the explicit diffusion equation is numerically stable, as
    %small disturbances will die down. But with an extra growth term, small
    %disturbances may lead to negative numbers, which are then inflated.
    %This is just a safety-measure.
            C = C + dt*rho.*abs(C).*(1-(abs(C)/k));
           %Cubic growth term?
           %C = C + dt*0.001*rho.*abs(C).*(1-abs(C)/k).*abs((C-(0.02*k)./C));
            plot(CC)  % ./k to normalise
            axis([0 132 0 (k)])
            pause(0.1)
            disp(t)
            ConcSum(t) = sum(C);
end
figure
plot(ConcSum)
%end
%%Surgery function
%cuts C for C>is and stiches back together the ends

Summary of current issues:

Working on segmenting T1 and T2 images from Simon to export to MATLAB and try to quantify difference in tumour spread for boundary.  

1. Need a different tool?  Check is FSL or MRTrix will do it, but I think they will only make a mask and brain extraction for diffusion data?  
2. analysis of the data in MATLAB, converting scalar field to meaning
3. differentiating between tumour and (say) bone of other dense material with same image value.
4. Maybe using FreeSurfer and/or 3D Slicer to do segmentatin and extraction.

Growing tumours:

Have found optimal values for stable solutions generally seem to exist at around Stability condition =  0.02 or thereabouts. Can simulate tumours quite effectively when parameters are right.

Simulation issues:  

1. Need to confirm how anisotropic the in vivo tumours are, scaling anisotropy appears to make things 'jagged' but still more like a spiky ball as opposed to strictly following any structures.  Investigate this further and visually compare with images.

Almost totally sure that the problem is too many negative eigenvalues in the tensor.  An analysis showed that around 8% of the eigenvalues of the tensor are negative (using tensor MouseTensorNyFiles, from Nyoman's data). 

Will check the other tensors and think on how to fix it.

 

Isotropic case for mouse brain did not yield any answers, although some tumours grew better especially when the number of time steps was increased. Tumours still stopped growing and concentration dropped to zero or negative.The upper video is of one of the best tumours, will run with same parameters to see if it can become any bigger.

A good example of an instability.   Will check data around this time that instability appears.