Solving a non linear equation using the Least squares method

3 views (last 30 days)
This concern MRI data and diffusion tensor imaging
I am trying to estimate D within each voxel using the data in data1.mat and i am not sure if i have done it correctly
Any help would appreciated
syms d_xx d_yy d_zz d_xy d_xz d_yz real
D = [d_xx, d_xy, d_xz; d_xy, d_yy, d_yz; d_xz, d_yz, d_zz]
syms x y z real
g = [x; y; z]
sympref('FloatingPointOutput',false);
syms b S_0
S = S_0*exp(-b*g'*D*g)
log(S)
logS = expand(1/2* g'*D*g)
Load data
load data1.mat
b_val = 1000;
xvec = g(:,1);
yvec = g(:,2);
zvec = g(:,3);
Svec = reshape(S, [], 64);
S0vec = reshape(S0, [], 1);
Construct A and b
A = [-xvec.^2/2, -xvec.*yvec, -xvec.*zvec, -yvec.*2/2, -yvec.*zvec, -zvec.^2/2];
b = log(Svec./S0vec);
c = zeros(6, size(S0vec, 1));
for i = 1:size(S0vec, 1)
c(:,i) = A \ b(i,:)';
end
% Tensor for each voxel
D = zeros(3, 3, size(S0vec, 1));
for i = 1:size(S0vec, 1)
D(:,:,i) = [c(1,i), c(2,i), c(3,i); c(2,i), c(4,i), c(5,i); c(3,i), c(5,i), c(6,i)];
end
D;

Answers (0)

Categories

Find more on Polymers in Help Center and File Exchange

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!