Why does my function that should modify cells from an input return exactly the input every time?

44 次查看（过去 30 天）
Skyler
Skyler 2024-12-7，23:49
评论： Stephen23 2024-12-8，6:51
I am attempting to model diffusion for an MRI project. In doing so, I created a function called diffuseCells() that takes several inputs and should modify the contents of the cells according to the rules in the functions. When I work through an individual example in the Command Window, everything seems to work. However, when I run the function it returns a cell structure that is exactly the same as the input structure. Any help would be greatly appreciated - I have worked on this single bug for hours and can't find the error. At this point I'm more or less banging my head against the wall. The function is below, and the entire script is below that.
The offending function:
%Inputs: cell array of protons, phantom, Gx(t),Gy(t), b0, dt,t
%Output: cell array of protons for t+dt
function t_cell = diffuseCells(p_cell,p, G1, G2,b0,dt,t)
%Calculate gradients
Gxc = G1(t)*linspace(round(-1*size(p_cell,2)/2), round(size(p_cell,2)/2), size(p_cell,2));
Gyc = G2(t)*linspace(round(-1*size(p_cell,1)/2), round(size(p_cell,1)/2), size(p_cell,1));
%Initialize temp cell storage for randomwalk calc
t_cell = cell(size(p_cell)); % Initialize
for i = 1:numel(p_cell)
t_cell{i} = p_cell{i}; % Copy each element
end
%t_cell = p_cell;
%at each timepoint iterate through cells and update spins and phases,
%then update position according to random walk
for i=1:size(p_cell,2) % x coord
for j=size(p_cell,1) % y coord
a = cell2mat(p_cell(i,j));
%update phase1 = mod(1,phase0 + spin (Hz) * dt (uS)),1), spin=4257 Hz/G * (b0 + Gx +Gy) (G) =
%gamma*field strength
a(1,:) = 4257*(b0+Gxc(i)+Gyc(j));
a(2,:) = mod(a(2,:)+a(1,:)*dt/(10^6),1);
%Storage for a after random walk
a_r=a;
%update position according to random walk
for k=1:size(a,2)
current = a(:,k);
[dx,dy] = randomwalk(p(i,j)); %use diffusion coefficients stored in phantom as n for random walk
if(p(i+dx,j+dy) == 0) %If moving out of object, stop
dx=0;
dy=0;
end
%move proton if it needs to be moved
if(dx ~= 0 || dy ~= 0)
%place in new cell
%t_cell{i+dx,j+dy} = [];
if(i+dx <= size(p_cell,2) && i+dx >= 1 && y+dy <= size(p_cell,1) && y+dy >=1)
t_cell{i+dx,j+dy} = [cell2mat(p_cell(i+dx,j+dy)),current];
end
%remove from old cell
a_r(:,k) = [];
end
end
%update cell with new proton entries
%t_cell{i,j} = [];
t_cell{i,j} = a_r;
end
end
end
Entire script:
%%
%define phantom, convert to cells
pt = phantom('Modified Shepp-Logan', 256);
%imagesc(p)
%Process phantom to have pixel values as diffusion coefficients
%unique(pt) = -0.0000,0,0.1,0.2,0.3,0.4,1 so multiply by ten to get
%diffusion coefficients
%NOTE: divide output by ten when displaying phantom to get good contrast
p = round(pt*10);
imagesc(pt/10)
p_c = makeCells(p);
%Define gradients and dimensions for k space traversal
FOVx = 25.6;FOVy = 19.2;dx = .4; dy = 2.4;
[G1,G2,G3,dt,Npts,Gx,Gy,kx,ky] = makeEPI(FOVx,FOVy,dx,dy);
b0 = 3*10000; %b0 = 3T * 10,000G/T
%Loop through time and simulate diffusion
%for t=1:size(Gx,1)
% nextCell = diffuseCells(p_c,p,Gx,Gy,b0,dt,t);
% p_c = nextCell;
%end
%%
%troubleshooting area
%[dx, dy] = randomwalk(0)
nextCell = diffuseCells(p_c,p,Gx,Gy,b0,dt,1);
%[dx,dy] = randomwalk(p(i,j));
%%
%Inputs: cell array of protons, phantom, Gx(t),Gy(t), b0, dt,t
%Output: cell array of protons for t+dt
function t_cell = diffuseCells(p_cell,p, G1, G2,b0,dt,t)
%Calculate gradients
Gxc = G1(t)*linspace(round(-1*size(p_cell,2)/2), round(size(p_cell,2)/2), size(p_cell,2));
Gyc = G2(t)*linspace(round(-1*size(p_cell,1)/2), round(size(p_cell,1)/2), size(p_cell,1));
%Initialize temp cell storage for randomwalk calc
t_cell = cell(size(p_cell)); % Initialize
for i = 1:numel(p_cell)
t_cell{i} = p_cell{i}; % Copy each element
end
%t_cell = p_cell;
%at each timepoint iterate through cells and update spins and phases,
%then update position according to random walk
for i=1:size(p_cell,2) % x coord
for j=size(p_cell,1) % y coord
a = cell2mat(p_cell(i,j));
%update phase1 = mod(1,phase0 + spin (Hz) * dt (uS)),1), spin=4257 Hz/G * (b0 + Gx +Gy) (G) =
%gamma*field strength
a(1,:) = 4257*(b0+Gxc(i)+Gyc(j));
a(2,:) = mod(a(2,:)+a(1,:)*dt/(10^6),1);
%Storage for a after random walk
a_r=a;
%update position according to random walk
for k=1:size(a,2)
current = a(:,k);
[dx,dy] = randomwalk(p(i,j)); %use diffusion coefficients stored in phantom as n for random walk
if(p(i+dx,j+dy) == 0) %If moving out of object, stop
dx=0;
dy=0;
end
%move proton if it needs to be moved
if(dx ~= 0 || dy ~= 0)
%place in new cell
%t_cell{i+dx,j+dy} = [];
if(i+dx <= size(p_cell,2) && i+dx >= 1 && y+dy <= size(p_cell,1) && y+dy >=1)
t_cell{i+dx,j+dy} = [cell2mat(p_cell(i+dx,j+dy)),current];
end
%remove from old cell
a_r(:,k) = [];
end
end
%update cell with new proton entries
%t_cell{i,j} = [];
t_cell{i,j} = a_r;
end
end
end
function [G1,G2,G3,dt,Npts,gx,gy,kx,ky] = makeEPI(FOVx,FOVy,dx,dy)
% Inputs
% FOVx = x-axis field of view in cm
% FOVy = y-axis field of view in cm
% dx = x-axis resolution in cm
% dy = y-axis resolution in cm
%
%
% Outputs (note: all amplitudes are positive)
% G1: G1 gradient amplitude in G/cm;
% G2: G2 gradient amplitude in G/cm;
% G3: G3 gradient amplitude in G/cm;
% dt: readout ADC sampling rate in microseconds;
% Npts: number of points in the gradient vector (sampled at dTS = 10 usec);
% gx: Npts x 1 column vector; x-gradient sampled at dTS;
% gy: Npts x 1 column vector; y-gradient sampled at dTS;
% kx = (Npts+1) x 1 column vector with kx values; sampled at dTS;
% ky = (Npts+1) x 1 column vector with ky values; sampled at dTS;
%
gam = 4257; %gamma/(2*pi) in Hz/G
G1dephase = 320; % width of readout dephaser in microseconds.
G1readout = 640; % width of each readout gradient in microseconds.
G2time = 10; % width of each phase encode blip in microseconds.
G3time = 320; % width of phase encode dephaser in microseconds.
dTS = G2time; % sampling rate for waveforms;
% Inputs
% FOVx = x-axis field of view in cm
% FOVy = y-axis field of view in cm
% dx = x-axis resolution in cm
% dy = y-axis resolution in cm
%
%
% Outputs (note: all amplitudes are positive)
% G1: G1 gradient amplitude in G/cm;
% G2: G2 gradient amplitude in G/cm;
% G3: G3 gradient amplitude in G/cm;
% dt: readout ADC sampling rate in microseconds;
% Npts: number of points in the gradient vector (sampled at dTS = 10 usec);
% gx: Npts x 1 column vector; x-gradient sampled at dTS;
% gy: Npts x 1 column vector; y-gradient sampled at dTS;
% kx = (Npts+1) x 1 column vector with kx values; sampled at dTS;
% ky = (Npts+1) x 1 column vector with ky values; sampled at dTS;
%
% Notes:
% 1) assume the gradients begin with t = 0 at the start of the dephasers.
% 2) The first element of each gradient vector describes the gradient for time interval [0 dTS]
% 3) the first element of each k-space vector should be zero -- i.e. start at
% the center of k-space.
% 4) You may assume that the gradients can change instantaneously.
gam = 4257; %gamma/(2*pi) in Hz/G
G1dephase = 320; % width of readout dephaser in microseconds.
G1readout = 640; % width of each readout gradient in microseconds.
G2time = 10; % width of each phase encode blip in microseconds.
G3time = 320; % width of phase encode dephaser in microseconds.
dTS = G2time; % sampling rate for waveforms;
% PUT YOUR CODE HERE; replace NaN or 1 values with correct values.
% Calculating gradient amplitudes (in G/cm)
G1 = 1/(dx*gam*G1readout*10^-6); % Gx = 1/(FOVx*gam*G1readout*10^-6)
G2 = 1/(FOVy*gam*G2time*10^-6); % G2: phase encode gradient amplitude, typically for phase encoding (y-direction)
G3 = 1/(2*dy*gam*G3time*10^-6);
%1/(dy*gam*G3time*10^-6); %depends on dy
%dy/(gam*G3time*10^-6);
%gam*2*pi/(FOVy*dy);
%(pi/2)/(gam*G3time);
%(pi/2)/(gam*dy*2*pi*G3time*10^-6);
%(pi/2)/(gam*dy*2*pi*G3time*10^-6); %90 degree shift
%1/(FOVy*gam*G3time*10^-6); % G3: phase encode dephaser, how to make? shift to pi/2?
% Time intervals and Npts calculation
%wx = 1/dx;
%wy = 1/dy;
%dt = G2time;
dt = 1/(gam*G1*FOVx)*10^6;
%1/(2*max(wx,wy))*10^6; % us
%Npts = round((G1dephase+(FOVy/dy)*(G1readout+G2time-1))/dt);
Nx = FOVx/dx;
Ny = FOVy/dy;
Npts = round((G1dephase+Ny*G1readout+(Ny-1)*G2time)/dTS);
%FOVx/dx
%(FOVx/dx)*(FOVy/dy);
%Npts = (4*(G1dephase+G1readout)+G1dephase)/dTs
gx = zeros(Npts, 1); % x-gradient
gy = zeros(Npts, 1); % y-gradient
kx = zeros(Npts + 1, 1); % k-space values for x-direction
ky = zeros(Npts + 1, 1); % k-space values for y-direction
Nd = round(G1dephase/dTS); %dephase
Nr = round(G1readout/dTS); %readout
Nc = Nd+1; %current
Nb = round(G2time/dTS); %blip
%Initialize dephase
gx(1:Nd) = -G1;
gy(1:round(G3time/dTS))=-G3;
for i=1:Ny
j = mod(i,2);
if(j==1)
gx(Nc:Nc+Nr-1)=G1;
else
gx(Nc:Nc+Nr-1)=-G1;
end
if(i<Ny)
gy(Nc+Nr:Nc+Nr+Nb-1)=G2;
end
Nc = Nc + Nr + Nb;
end
gx(G1dephase/dt+1)=G1;
% Calculate k-space vectors (kx, ky) based on the gradient waveforms
for i = 2:Npts + 1
kx(i) = kx(i-1) + gam * gx(i-1) * (dt * 1e-6); % kx calculation in Hz/cm
ky(i) = ky(i-1) + gam * gy(i-1) * (dt * 1e-6); % ky calculation in Hz/cm
end
figure; clf;
nexttile;
dTS = 10;
totalTime = Npts*dTS
taxis = [0:dTS:totalTime];
plot(taxis(2:end),gx);
title("Gx vs Time")
xlabel("Time (uS)")
ylabel("Gradient amplitude (G/cm)")
legend("Gx")
nexttile;grid;
plot(taxis(2:end),gy);
title("Gy vs Time")
xlabel("Time (uS)")
ylabel("Gradient amplitude (G/cm)")
legend("Gx","Gy")
nexttile;grid;
plot(kx,ky,'-o');
title("Kx vs Ky")
xlabel("Kx (1/cm)")
ylabel("Ky (1/cm)")
grid;
end
%Inputs: phantom with pixel values as diffusion coefficients (must be
%integers)
%Output: cell array with 3 protons per pixel
function p_cell = makeCells(p)
p_cell = cell(size(p));
for i=1:size(p,1)
for j=1:size(p,2)
%initialize spins as 0, phases as 0
p_cell{i,j} = [1,0,3;0,0,0];
end
end
end
%Inputs: # of steps to take
%Output: change in x and change in y position
function [dx, dy] = randomwalk(n)
dx = 0;
dy = 0;
for i=1:n
%Uniform random variate [0 1]
var = rand();
%choose one of 8 directions randomly
if(0 <= var && var <= 1/8)
%d='UL';
dx = dx-1;
dy = dy+1;
elseif(1/8 < var && var <= 1/4)
%d='U';
dy = dy+1;
elseif(1/4 < var && var <= 3/8)
%d='UR';
dy = dy+1;
dx = dx+1;
elseif(3/8 < var && var <= 1/2)
%d='R';
dx = dx+1;
elseif(1/2 < var && var <= 5/8)
%d='DR';
dy = dy-1;
dx = dx+1;
elseif(5/8 < var && var <= 3/4)
%d='D';
dy = dy-1;
elseif(3/4< var && var <= 7/8)
%d='DL';
dy = dy-1;
dx = dx-1;
else
%d='UL';
dy = dy+1;
dx = dx-1;
end
end
end

请先登录，再回答此问题。

采纳的回答

Stephen23
Stephen23 2024-12-8，3:48
编辑：Stephen23 2024-12-8，4:01
This line:
for j=size(p_cell,1) % y coord
should be
for j=1:size(p_cell,1) % y coord
% ^^
After that you will need to fix some y's which should be j's. And some indexing.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by