预白化MRI：SENSE重建后噪声协方差非单位矩阵问题排查

我在MATLAB中基于仿真体模数据开展SENSE MRI重建研究，向多线圈k空间数据添加含已知空间不变协方差矩阵Σ的复高斯相关噪声后，执行白化变换以将噪声协方差矩阵近似为单位矩阵。但重建后从图像背景块估计的噪声协方差矩阵Ψ̃并未接近单位矩阵，与预期不符。

已完成操作：

• 生成带已知协方差矩阵Σ的复高斯相关噪声并添加至k空间数据；
• 从图像空白区域计算噪声协方差矩阵Ψ；
• 通过Ψ的特征分解构造白化矩阵W = D⁻¹ᐟ² * V'；
• 利用块对角克罗内克积（W_block = kron(speye(N), W)）应用W；
• 使用白化后的编码矩阵Ã和测量向量Ỹ执行标准SENSE伪逆重建。

预期结果：

W * Psi * W' ≈ I
Psi_tilde ≈ I

但第二个条件未满足，需排查问题出在白化步骤还是后续重建流程。以下是完整代码：

% Just whiteing for phantom data
clear;  clc;

FOV = 256;
Ny = FOV;
Nx = FOV;
Nc = 8;
R=2;
alias_span = Ny/R;
data = load('modifiedshep.mat'); 
whos('-file', 'modifiedshep.mat')
raw = data.raw;
DATA = raw;
N_acs_y = 20;  
N_acs_x = 20;  
c_sens = estimate_sensitivity_acs(DATA, N_acs_y, N_acs_x);
magnitude = abs(DATA);                
phase = angle(DATA);                  
max_mag = max(magnitude(:));       
DATA = (magnitude / max_mag) .* exp(1j * phase);

%% 
sigma2 = 0.005^2;     % Variance of noise in each coil (diagonal elements)
rho = 0.99;            % Correlation coefficient between coils (0 = uncorrelated, 1 = full correlation)
corr = 1;             % Set to 1 to enable correlation

% Build covariance matrix [Nc x Nc]
if corr == 0
    Sigma = sigma2 * eye(Nc);   % No correlation
else
    Sigma = sigma2 * (eye(Nc) + rho * (ones(Nc) - eye(Nc)));  % With correlation
end

% Eigendecomposition for noise shaping
[q, a] = eig(Sigma);
t = q * sqrt(a);    % Transform: white -> correlated noise

% Generate white complex Gaussian noise [Ny x Nx x Nc]
noise_white = randn(Ny, Nx, Nc) + 1i * randn(Ny, Nx, Nc);  

% Apply spatially constant covariance per pixel
noise_corr = zeros(Ny, Nx, Nc);
for i = 1:Ny
    for j = 1:Nx
        noise_corr(i,j,:) = t * squeeze(noise_white(i,j,:));
    end
end

% Add noise to DATA (before undersampling!)
DATA = DATA + noise_corr;
idx = setdiff(1:Ny,1:R:Ny);
DATA(idx,:,:) = 0;

coil_images = ifft2c(DATA);
size(coil_images)

%% 
num_patches = 4;
patch_size =  20;
samples_per_patch = patch_size * patch_size;
total_samples = samples_per_patch * num_patches;
num_coils = size(coil_images, 3);

noise_matrix = zeros(total_samples, num_coils);

% Define corners
corners = {                  % top-left
    [1, patch_size, 1, patch_size];                      % top-left
    [1, patch_size, Ny-patch_size+1, Ny];                % top-right
    [Nx-patch_size+1, Nx, 1, patch_size];                 % bottom-left
    [Nx-patch_size+1, Nx, Ny-patch_size+1, Ny]           % bottom-right
};

% Extract patches and fill noise_matrix
for p = 1:num_patches
    r_start = corners{p}(1);
    r_end = corners{p}(2);
    c_start = corners{p}(3);
    c_end = corners{p}(4);
    
    patch = coil_images(r_start:r_end, c_start:c_end, :);       
    patch_reshaped = reshape(patch, [], num_coils);               
    
    idx_start = (p-1)*samples_per_patch + 1;
    idx_end = p*samples_per_patch;
    
    noise_matrix(idx_start:idx_end, :) = patch_reshaped;
end
noise_matrix = reshape(coil_images(r_start:r_end, c_start:c_end, :), [], Nc);
% Now compute covariance
Psi = (noise_matrix' * noise_matrix)/size(noise_matrix,1);
figure(1),
imagesc(abs(Psi));
colorbar;
title('Magnitude of \Psi');
xlabel('Coil #'); ylabel('Coil #');
%% 


[V, D] = eig(Psi);         % Stabilization constant
D = diag(D);              % Extract eigenvalues

% --- Step 3: Compute inverse sqrt of eigenvalues with stabilization ---
D_inv_sqrt = diag(1 ./ sqrt(D)); 

% --- Step 4: Whitening matrix ---
W =  D_inv_sqrt * V';  

check = W * Psi * W';
figure(2),
imagesc(abs(check));
colormap jet;
colorbar;
title('w*psi*W');
xlabel('Coil #'); ylabel('Coil #');
%% 
A = zeros(Nc, R, Ny/R, Nx);
for x = 1:Nx
    for y = 1:Ny/R
        for r = 1:R
            yy = y + (r-1)*alias_span;
            for l = 1:Nc
                A(l, r, y, x) = c_sens(yy, x, l);  % unwhitened
            end
        end
    end
end
Y = zeros(Nc,1,Ny/R,Nx);
for x=1:Nx
    for y=1:Ny/R
        for l=1:Nc
            Y(l,1,y,x) = coil_images(y,x,l);
        end
    end
end
N_pixels = (Ny / R) * Nx;   % Total number of unfolded pixels
A_big = zeros(Nc * N_pixels, R);
Y_big = zeros(Nc * N_pixels, 1);

count = 1;
for x = 1:Nx
    for y = 1:(Ny/R)
        A_block = A(:,:,y,x);  % [Nc × R]
        Y_block = Y(:,:,y,x);  % [Nc × 1]

        row_start = (count - 1) * Nc + 1;
        row_end   = count * Nc;

        A_big(row_start:row_end, :) = A_block;
        Y_big(row_start:row_end, :) = Y_block;

        count = count + 1;
    end
end
W_block = kron(speye(N_pixels), W); 
A_tilde_big = W_block * A_big;  % [Nc*N_pixels × R]
Y_tilde_big = W_block * Y_big; 
A_tilde = zeros(Nc, R, Ny/R, Nx);
Y_tilde = zeros(Nc, 1, Ny/R, Nx);

count = 1;
for x = 1:Nx
    for y = 1:(Ny/R)
        row_start = (count - 1) * Nc + 1;
        row_end   = count * Nc;

        A_tilde(:,:,y,x) = A_tilde_big(row_start:row_end, :);
        Y_tilde(:,:,y,x) = Y_tilde_big(row_start:row_end);

        count = count + 1;
    end
end
X = zeros(R, Ny/R, Nx);  % Store the reconstructed result

for x = 1:Nx
    for y = 1:(Ny/R)
        A_local = squeeze(A_tilde(:,:,y,x));   % [Nc x R]
        y_local = squeeze(Y_tilde(:,:,y,x));   % [Nc x 1]

        % Solve using pseudoinverse (more stable than manual inverse)
        x_local = pinv(A_local' * A_local) * (A_local' * y_local);

        % Alternatively: x_local = (A_local' \ (A_local' * y_local));  % if A is well-conditioned

        X(:,y,x) = x_local;  % Store result
    end
end

recon = zeros(Ny, Nx);  % Final unfolded image

alias_span = Ny / R;  % How far apart the aliased pixels were

for x = 1:Nx
    for y = 1:(Ny/R)
        for r = 1:R
            yy = y + (r-1)*alias_span;
            recon(yy, x) = X(r, y, x);
        end
    end
end

%% 

figure(4);
imagesc(abs(recon));
colormap gray;
axis image off;
title('Reconstructed Image (Magnitude)');

%% 
coil_images_white = zeros(Ny, Nx, Nc);  % Whitened coil images
alias_span = Ny / R;

for x = 1:Nx
    for y = 1:(Ny/R)
        for r = 1:R
            yy = y + (r-1)*alias_span;  % actual y-position in full image
            for c = 1:Nc
                coil_images_white(yy, x, c) = Y_tilde(c,1,y,x);
            end
        end
    end
end


%% 
num_patches = 4;
patch_size =  20;
samples_per_patch = patch_size * patch_size;
total_samples = samples_per_patch * num_patches;
num_coils = size(coil_images, 3);

noise_matrix = zeros(total_samples, num_coils);

% Define corners
corners = {                  % top-left
    [1, patch_size, 1, patch_size];                      % top-left
    [1, patch_size, Ny-patch_size+1, Ny];                % top-right
    [Nx-patch_size+1, Nx, 1, patch_size];                 % bottom-left
    [Nx-patch_size+1, Nx, Ny-patch_size+1, Ny]           % bottom-right
};

% Extract patches and fill noise_matrix
for p = 1:num_patches
    r_start = corners{p}(1);
    r_end = corners{p}(2);
    c_start = corners{p}(3);
    c_end = corners{p}(4);
    
    patch = coil_images_white(r_start:r_end, c_start:c_end, :);       
    patch_reshaped = reshape(patch, [], num_coils);               
    
    idx_start = (p-1)*samples_per_patch + 1;
    idx_end = p*samples_per_patch;
    
    noise_matrix_tilde(idx_start:idx_end, :) = patch_reshaped;
end

% Now compute covariance
Psi_tilde = (noise_matrix_tilde' * noise_matrix_tilde)/size(noise_matrix_tilde,1);
figure(5),
imagesc(abs(Psi_tilde));

colorbar; colormap jet;
title('Magnitude of \Psi');
xlabel('Coil #'); ylabel('Coil #');

%% 


function out = ifft2c(input)
    out = fftshift(ifft(ifftshift(input,1),[],1),1);
    out = fftshift(ifft(ifftshift(out,2),[],2),2);
end

function c_sens = estimate_sensitivity_acs(raw, N_acs_y, N_acs_x)

    [Ny, Nx, Nc] = size(raw);

    center_y = floor(Ny / 2) + 1;
    center_x = floor(Nx / 2) + 1;

    acs_start_y = center_y - floor(N_acs_y / 2);
    acs_end_y   = acs_start_y + N_acs_y - 1;
    acs_start_x = center_x - floor(N_acs_x / 2);
    acs_end_x   = acs_start_x + N_acs_x - 1;


    acs_kspace = raw(acs_start_y:acs_end_y, acs_start_x:acs_end_x, :);


    hamming_y = hamming(N_acs_y);
    hamming_x = hamming(N_acs_x);
    hamming_2D = hamming_y * hamming_x';
    for c = 1:Nc
        acs_kspace(:,:,c) = acs_kspace(:,:,c) .* hamming_2D;
    end

    acs_padded = zeros(Ny, Nx, Nc);
    acs_padded(acs_start_y:acs_end_y, acs_start_x:acs_end_x, :) = acs_kspace;


    coil_imgs = zeros(Ny, Nx, Nc);
    for c = 1:Nc
        coil_imgs(:,:,c) = ifftshift(ifft2(fftshift(acs_padded(:,:,c))));
    end


    sos_img = sqrt(sum(abs(coil_imgs).^2, 3));
    sos_img(sos_img == 0) = 1; 
    c_sens = zeros(Ny, Nx, Nc);
    for c = 1:Nc
        c_sens(:,:,c) = coil_imgs(:,:,c) ./ sos_img;
    end
end

核心问题排查及修正 #

1. 噪声协方差矩阵Ψ的样本提取错误
   循环填充4个背景块到noise_matrix后，被一行代码覆盖为仅最后一个块的样本：

   noise_matrix = reshape(coil_images(r_start:r_end, c_start:c_end, :), [], Nc);
   

   需删除该行代码，确保使用全部4个背景块的样本计算Ψ。

2. 白化后线圈图像构造逻辑错误
   Y_tilde是欠采样k空间对应的图像域数据，仅包含欠采样点，不能直接映射为全图像的白化线圈数据。正确做法是对原始全线圈图像逐像素应用白化矩阵W：

   % 替换原有coil_images_white构造代码
   coil_images_white = zeros(Ny, Nx, Nc);
   for i = 1:Ny
       for j = 1:Nx
           coil_images_white(i,j,:) = W * squeeze(coil_images(i,j,:));
       end
   end
   

3. Ψ̃计算变量未初始化
   计算noise_matrix_tilde前未初始化，需添加：

   noise_matrix_tilde = zeros(total_samples, num_coils);
   

4. 特征分解数值稳定性不足
   极小特征值会导致白化矩阵数值异常，需添加阈值截断：

   eps_val = 1e-6; % 可调整的数值阈值
   D(D < eps_val) = eps_val;
   D_inv_sqrt = diag(1 ./ sqrt(D));
   

内容的提问来源于stack exchange，提问作者maria