0

我正在学习 MRI 并行成像中的图像重建。我正在 Maltab 中编写关于论文“并行磁共振成像的通用方法（DK Sodickson，2001）”的混合重建程序。我的代码仅适用于欠采样 M=2，但不适用于 M=4。我不知道为什么。你可以帮帮我吗？

clear;clc
Nc=8;FOV=256;Nx=FOV;Ny=FOV;
ph=phantom('modified shepp-logan',FOV);
c_sens=zeros(Ny,Nx,Nc);
for n=1:Nc
for i=1:Ny
    for j=1:Nx
        mean=(FOV/(Nc+1))*n;
        var=30;
        a = 1/(var*(2*pi)^(0.5));
        b = (i-mean)^2;
        d = 2*((var)^2);
        k = (-1)*(b/d);
        c_sens(i,j,n)=a*exp(k);
    end
end
end
c_img=zeros(Ny,Nx,Nc);
for n=1:Nc
c_img(:,:,n) = ph.*c_sens(:,:,n); 
end
raw=fftshift(fft2(fftshift(c_img)));
%-----------------------------------------------------------------
M=4;
sampling_location = (1:M:Ny); % Ny/2+1 is sampled.
reduced_k_space = zeros(Ny/M, Nx, Nc);
for n = 1:Nc
reduced_k_space(:, :, n)=raw(sampling_location,:, n);
end
block_size=2;
total_block = Nc*block_size;        %total number of subblocks
exps = zeros(block_size, Ny);
y = -128:127;
for m = 1:block_size
exps(m,:) = exp(-2*pi*1i*(m-1)*y/Ny);
end
S_comp = zeros(Nx,Ny);
for x=1:Nx
C = [];
for n=1:Nc
    C = [C; c_sens(:,x,n)'];    
end
    
B = [];
tmp_b = [];

%for ky = 0:M:Ny
for ky = 0:M:(Ny/M) - 1
    for i=1:Nc
        temp_B = C(i,:) .* exp(-2*pi*1i*ky*y/Ny); % coil sensitivity * gradient modulation 
        tmp_b = [tmp_b; temp_B];
    end
    B = [B; tmp_b];
    tmp_b=[];
end

for h=1:total_block/2
    b_comb(h,:) = B(h,:);
end 

% Weights
weights = exps*pinv(b_comb);

permuted_raw = permute(reduced_k_space, [3, 1, 2]);
reshaped_raw = reshape(permuted_raw, Nc*Nx/M, Ny);
        
for n=1:(Ny/block_size)
    S_comp(block_size*(n-1)+1 : block_size*n, x) = weights * reshaped_raw((total_block/2)*(n-1) + 1 : (total_block/2)*n, x);
end
end

% final reconstructed image
reconimg = ifftshift(ifft2(ifftshift(S_comp)));
figure(1)
imshow(abs(reconimg),[])
4

0 回答 0