Code covered by the BSD License

Highlights fromMatlab code for my Graduate Thesis

from Matlab code for my Graduate Thesis by Troy
Numerically solves the diffusion equations as it pertains to medical imaging.

vvmodeuler3Dneumann(varargin)
```function [w1, w2, w3] = vvmodeuler3Dneumann(varargin)
% crankNicolson: uses Crank-Nicolson algorithm to approximate the solution
%   to the parabolic PDE:
%       u_{t}(x,t) - alpha^2 u_{xx}(x,t) = 0, 0<x<L, 0<t<T
%   subject to the boundary conditions
%       u(0,t) = u(L,t) = 0, 0<t<T
%   and the initial conditions
%       u(x,0) = f(x), 0<=x<=L
%
% arguments:
%   L (vector 1by3) - upper bound of spatial (x) variable
%       (Default L = [1,1,1])
%   T (scalar) - upper bound of time (t) variable
%       (Default T = 1)
%   alpha (scalar) - square root of coefficient of u_{xx} term
%       (Default alpha = 1)
%   m (vector 1by3) - number of discrete spatial intervals
%       (Default m = [100,100,100])
%   n (scalar) - number of discrete time intervals
%       (Default n = 100)
%
%   w (m+1 x m+1 x n) - approximation to u(x,t) at discrete space/time positions
%

% author: Troy J. Winkstern
% email: tjw8191@rit.edu
% date: 30 Jan 2011

% parse input arguments
[L,T,alpha,m,n] = parseInputs(varargin{:});
L_x = L(1);
L_y = L(2);
L_z = L(3);
m_x = m(1);
m_y = m(2);
m_z = m(3);

% initialize h, k, lambda, and w
h_x = L_x/m_x;
h_y = L_y/m_y;
h_z = L_z/m_z;
k = T./n;
lambda_x = (alpha.^2).*k./(h_x.^2);
lambda_y = (alpha.^2).*k./(h_y.^2);
lambda_z = (alpha.^2).*k./(h_z.^2);
w1 = zeros(2*m_y+1,2*m_x+1,2*m_z+1,n);
w2 = zeros(2*m_y+1,2*m_x+1,2*m_z+1,n);
w3 = zeros(2*m_y+1,2*m_x+1,2*m_z+1,n);

% initialize w1
for i=1:m_x+1
for j=1:m_y+1
for k=1:m_z+1

w1(j,i,k,1) = cos(pi.*h_x.*(i-1)./L_x).*cos(pi.*h_y.*(j-1)./L_y)...
.*cos(pi.*h_z.*(k-1)./L_z);

end
end
end

% initialize w2
for i=1:m_x+1
for j=1:m_y+1
for k=1:m_z+1

w1(j,i,k,1) = cos(pi.*h_x.*(i-1)./L_x).*cos(pi.*h_y.*(j-1)./L_y)...
.*cos(pi.*h_z.*(k-1)./L_z);

end
end
end

% initialize w3
for i=1:m_x+1
for j=1:m_y+1
for k=1:m_z+1

w1(j,i,k,1) = cos(pi.*h_x.*(i-1)./L_x).*cos(pi.*h_y.*(j-1)./L_y)...
.*cos(pi.*h_z.*(k-1)./L_z);

end
end
end

for t=1:(n-1)
f1 = laplacian(w1(:,:,:,t),h_y,h_x,h_z);
y1 = w1(:,:,:,t) + (k/2)*f1;

f2 = laplacian(w2(:,:,:,t),h_y,h_x,h_z);
y2 = w2(:,:,:,t) + (k/2)*f2;

f3 = laplacian(w3(:,:,:,t),h_y,h_x,h_z);
y3 = w3(:,:,:,t) + (k/2)*f3;

% went out to length 2m
X1 = [sqrt(1./(2*m_z)).*real(fft(y1, 2*m_y,1)):zeros(1,201,1)];

X2 = [sqrt(1./(2*m_z)).*real(fft(y2, 2*m_y,1)):zeros(1,201,1)];

X3 = [sqrt(1./(2*m_z)).*real(fft(y3, 2*m_y,1)):zeros(1,201,1)];

Y1 = [sqrt(1./(2*m_z)).*real(fft(X1, 2*m_x,2)):zeros(201,1,1)];

Y2 = [sqrt(1./(2*m_z)).*real(fft(X2, 2*m_x,2)):zeros(201,1,1)];

Y3 = [sqrt(1./(2*m_z)).*real(fft(X3, 2*m_x,2)):zeros(201,1,1)];

V1 = [sqrt(1./(2*m_z)).*real(fft(Y1, 2*m_z,3)):zeros(1,1,201)];

V2 = [sqrt(1./(2*m_z)).*real(fft(Y2, 2*m_z,3)):zeros(1,1,201)];

V3 = [sqrt(1./(2*m_z)).*real(fft(Y3, 2*m_z,3)):zeros(1,1,201)];

% angles divided by m
a_x = repmat((pi./m_x).*(0:m_x),[m_y+1,1,m_z+1]);
a_y = repmat((pi./m_y).*(0:m_y).',[1,m_x+1,m_z+1]);
a_z = repmat(reshape((pi./m_z).*(0:m_z),[1,1,m_z+1]),[m_y+1,m_x+1,1]);

% changed c to accomodate for 3-D
c1 = 1 + ((lambda_x)*(lambda_y)*(lambda_z)*(h_x)^2*(h_y)^2*(h_z)^2) - ...
(1/2)*((lambda_x)*cos(a_x)*(h_x)^2).*((lambda_y)*cos(a_y)*(h_y)^2) ...
.*((lambda_z)*cos(a_z)*(h_z)^2) - (1/2)*(((lambda_x)*cos(a_x)*(h_x)^2).^2)...
.*(((lambda_y)*cos(a_y)*(h_y)^2).^2).*(((lambda_z)*cos(a_z)*(h_z)^2).^2);

c2 = 1 + ((lambda_x)*(lambda_y)*(lambda_z)*(h_x)^2*(h_y)^2*(h_z)^2) - ...
(1/2)*((lambda_x)*cos(a_x)*(h_x)^2).*((lambda_y)*cos(a_y)*(h_y)^2) ...
.*((lambda_z)*cos(a_z)*(h_z)^2) - (1/2)*(((lambda_x)*cos(a_x)*(h_x)^2).^2)...
.*(((lambda_y)*cos(a_y)*(h_y)^2).^2).*(((lambda_z)*cos(a_z)*(h_z)^2).^2);

c3 = 1 + ((lambda_x)*(lambda_y)*(lambda_z)*(h_x)^2*(h_y)^2*(h_z)^2) - ...
(1/2)*((lambda_x)*cos(a_x)*(h_x)^2).*((lambda_y)*cos(a_y)*(h_y)^2) ...
.*((lambda_z)*cos(a_z)*(h_z)^2) - (1/2)*(((lambda_x)*cos(a_x)*(h_x)^2).^2)...
.*(((lambda_y)*cos(a_y)*(h_y)^2).^2).*(((lambda_z)*cos(a_z)*(h_z)^2).^2);

% went out to length 2m
Z1 = sqrt(1./(2*m_y)).*real(fft(V1./c1, 2*m_y,1));
Z1 = Z1([1:2*m_y 1],:,:);

Z2 = sqrt(1./(2*m_y)).*real(fft(V2./c2, 2*m_y,1));
Z2 = Z2([1:2*m_y 1],:,:);

Z3 = sqrt(1./(2*m_y)).*real(fft(V3./c3, 2*m_y,1));
Z3 = Z3([1:2*m_y 1],:,:);

R1 = sqrt(1./(2*m_x)).*real(fft(Z1, 2*m_x,2));
R1 = R1(:,[1:2*m_x 1],:);

R2 = sqrt(1./(2*m_x)).*real(fft(Z2, 2*m_x,2));
R2 = R2(:,[1:2*m_x 1],:);

R3 = sqrt(1./(2*m_x)).*real(fft(Z3, 2*m_x,2));
R3 = R3(:,[1:2*m_x 1],:);

S1 = sqrt(2./m_z).*imag(fft(R1, 2*m_z,3));
S1 = S1(:,:,[1:2*m_z 1]);

S2 = sqrt(2./m_z).*imag(fft(R2, 2*m_z,3));
S2 = S2(:,:,[1:2*m_z 1]);

S3 = sqrt(2./m_z).*imag(fft(R3, 2*m_z,3));
S3 = S3(:,:,[1:2*m_z 1]);

w1(:,:,:,t+1) = S1;
w2(:,:,:,t+1) = S2;
w3(:,:,:,t+1) = S3;
end

%crop w1,w2,w3
w1 = w1(1:m_y+1,1:m_x+1,1:m_z+1,:);
w2 = w2(1:m_y+1,1:m_x+1,1:m_z+1,:);
w3 = w3(1:m_y+1,1:m_x+1,1:m_z+1,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction parseInputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L,T,alpha,m,n] = parseInputs(varargin)

% check number of input arguments
nargs = length(varargin);
error(nargchk(0,5,nargs));

% get/set n
if nargs<5
n = [];
else
n = varargin{5};
end
if isempty(n)
n = 10;
end

% check n has valid value
if (numel(n)>1) || (n<1) || ~isequal(round(n),n)
error('n must be a positive integer.');
end

% get/set m
if nargs<4
m = [];
else
m = varargin{4};
end
if isempty(m)
m = [10,10,10];
end

% check m has valid value
if (numel(m)~=3) || any(m<1) || any(~isequal(round(m),m))
error('m must be a vector of positive integers.');
end

% get/set alpha
if nargs<3
alpha = [];
else
alpha = varargin{3};
end
if isempty(alpha)
alpha = 1;
end

% check alpha has valid value
if numel(alpha)>1
error('alpha must be a real scalar.');
end

% get/set T
if nargs<2
T = [];
else
T = varargin{2};
end
if isempty(T)
T = 1;
end

% check T has valid value
if (numel(T)>1) || (T<=0)
error('T must be a positive real scalar.');
end

% get/set L
if nargs<1
L = [];
else
L = varargin{1};
end
if isempty(L)
L = [1,1,1];
end

% check L has valid value
if (numel(L)~=3) || any(L<=0)
error('L must be a positive vector with 3 elements.');
end

% warning if L not integer
if any(~isequal(round(L),L))
warning('L should be an vector of integers so that the inital condition takes on the value of 0 for t = 0 and t = T.');
end

%another subfunction
function y = laplacian(x,h_y,h_x,h_z)
[m,n,p] = size(x);
y = zeros(m,n,p);

%compute the 8 corners
y(1,1,1) = (2*x(1,2,1)-2*x(1,1,1))/(h_x.^2)+(2*x(2,1,1)-2*x(1,1,1))/(h_y.^2)...
+(2*x(1,1,2)-2*x(1,1,1))/(h_z.^2);
y(m,1,1) = (2*x(m,2,1)-2*x(m,1,1))/(h_x.^2)+(2*x(m-1,1,1)-2*x(m,1,1))/(h_y.^2)...
+(2*x(m,1,2)-2*x(m,1,1))/(h_z.^2);
y(1,n,1) = (2*x(1,n-1,1)-2*x(1,n,1))/(h_x.^2)+(2*x(2,n,1)-2*x(1,n,1))/(h_y.^2)...
+(2*x(1,n,2)-2*x(1,n,1))/(h_z.^2);
y(1,1,p) = (2*x(1,2,p)-2*x(1,1,p))/(h_x.^2)+(2*x(2,1,p)-2*x(1,1,p))/(h_y.^2)...
+(2*x(1,1,p-1)-2*x(1,1,p))/(h_z.^2);
y(m,n,1) = (2*x(m,n-1,1)-2*x(m,n,1))/(h_x.^2)+(2*x(m-1,n,1)-2*x(m,n,1))/(h_y.^2)...
+(2*x(m,n,2)-2*x(m,n,1))/(h_z.^2);
y(1,n,p) = (2*x(1,n-1,p)-2*x(1,n,p))/(h_x.^2)+(2*x(2,n,p)-2*x(1,n,p))/(h_y.^2)...
+(2*x(1,n,p-1)-2*x(1,n,p))/(h_z.^2);
y(m,1,p) = (2*x(m,2,p)-2*x(m,1,p))/(h_x.^2)+(2*x(m-1,1,p)-2*x(m,1,p))/(h_y.^2)...
+(2*x(m,1,p-1)-2*x(m,1,p))/(h_z.^2);
y(m,n,p) = (2*x(m,n-1,p)-2*x(m,n,p))/(h_x.^2)+(2*x(m-1,n,p)-2*x(m,n,p))/(h_y.^2)...
+(2*x(m,n,p-1)-2*x(m,n,p))/(h_z.^2);

%compute the 12 edges
for i = 2:n-1
y(1,i,1) = (x(1,i-1,1)+x(1,i+1,1)-2*x(1,i,1))/(h_x.^2) + ...
(2*x(2,i,1)-2*x(1,i,1))/(h_y.^2) + (2*x(1,i,2)-2*x(1,i,1))/(h_z.^2);
end

for i = 2:n-1
y(m,i,1) = (x(m,i-1,1)+x(m,i+1,1)-2*x(m,i,1))/(h_x.^2) + ...
(2*x(m-1,i+1,1)-2*x(m,i,1))/(h_y.^2) + (2*x(m,i,2)-2*x(m,i,1))/(h_z.^2);
end

for i = 2:n-1
y(1,i,p) = (x(1,i-1,p)+x(1,i+1,p)-2*x(1,i,p))/(h_x.^2) + ...
(2*x(2,i+1,p)-2*x(1,i,p))/(h_y.^2) + (2*x(1,i,p-1)-2*x(1,i,p))/(h_z.^2);
end

for i = 2:n-1
y(m,i,p) = (x(m,i-1,p)+x(m,i+1,p)-2*x(m,i,p))/(h_x.^2) + ...
(2*x(m-1,i+1,p)-2*x(m,i,p))/(h_y.^2) + (2*x(m,i,p-1)-2*x(m,i,p))/(h_z.^2);
end

for j = 2:m-1
y(j,1,1) = (x(j-1,1,1)+x(j+1,1,1)-2*x(j,1,1))/(h_y.^2) + ...
(2*x(j,2,1)-2*x(j,1,1))/(h_x.^2) + (2*x(j,1,2)-2*x(j,1,1))/(h_z.^2);
end

for j = 2:m-1
y(j,n,1) = (x(j-1,n,1)+x(j+1,n,1)-2*x(j,n,1))/(h_y.^2) + ...
(2*x(j,n-1,1)-2*x(j,n,1))/(h_x.^2) + (2*x(j,n,2)-2*x(j,n,1))/(h_z.^2);
end

for j = 2:m-1
y(j,1,p) = (x(j-1,1,p)+x(j+1,1,p)-2*x(j,1,p))/(h_y.^2) + ...
(2*x(j,2,p)-2*x(j,1,p))/(h_x.^2) + (2*x(j,1,p-1)-2*x(j,1,p))/(h_z.^2);
end

for j = 2:m-1
y(j,n,p) = (x(j-1,n,p)+x(j+1,n,p)-2*x(j,n,p))/(h_y.^2) + ...
(2*x(j,n-1,p)-2*x(j,n,p))/(h_x.^2) + (2*x(j,n,p-1)-2*x(j,n,p))/(h_z.^2);
end

for k = 2:p-1
y(1,1,k) = (x(m,1,k-1)+x(m,1,k+1)-2*x(m,1,k))/(h_z.^2) + ...
(2*x(m,2,k)-2*x(m,1,k))/(h_x.^2) + (2*x(m-1,1,k)-2*x(m,1,k))/(h_y.^2);
end

for k = 2:p-1
y(m,1,k) = (x(1,1,k-1)+x(1,1,k+1)-2*x(1,1,k))/(h_z.^2) + ...
(2*x(1,2,k)-2*x(1,1,k))/(h_x.^2) + (2*x(2,1,k)-2*x(1,1,k))/(h_y.^2);
end

for k = 2:p-1
y(1,n,k) = (x(1,n,k-1)+x(1,n,k+1)-2*x(1,n,k))/(h_z.^2) + ...
(2*x(1,n-1,k)-2*x(1,n,k))/(h_x.^2) + (2*x(2,n,k)-2*x(1,n,k))/(h_y.^2);
end

for k = 2:p-1
y(m,n,k) = (x(m,n,k-1)+x(m,n,k+1)-2*x(m,n,k))/(h_z.^2) + ...
(2*x(m,n-1,k)-2*x(m,n,k))/(h_x.^2) + (2*x(m-1,n,k)-2*x(m,n,k))/(h_y.^2);
end

%compute the inside
for i = 2:n-1
for j = 2:m-1
for k = 2:p-1

y(j,i,k) = (x(j-1,i,k)+x(j+1,i,k)-2*x(j,i,k))/(h_y.^2) + ...
(x(j,i-1,k)+x(j,i+1,k)-2*x(j,i,k))/(h_x.^2) + ...
(x(j,i,k-1)+x(j,i,k+1)-2*x(j,i,k))/(h_z.^2);

end
end
end```