cs

cg.py

+'''
+Author: zhuwanjie 2268677665@qq.com
+Date: 2025-04-25 17:07:30
+LastEditors: zhuwanjie 2268677665@qq.com
+LastEditTime: 2025-04-27 09:50:41
+FilePath: \egg\cg.py
+Description: 这是默认设置,请设置`customMade`, 打开koroFileHeader查看配置 进行设置: https://github.com/OBKoro1/koro1FileHeader/wiki/%E9%85%8D%E7%BD%AE
+'''
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+from scipy.fftpack import fft2, ifft2
+
+
+def generate_sampling_mask(shape, sampling_rate):
+    """生成随机采样掩码"""
+    num_samples = int(sampling_rate * shape[0] * shape[1])
+    mask = np.zeros(shape, dtype=np.float32)
+    indices = np.random.choice(shape[0] * shape[1], num_samples, replace=False)
+    mask.flat[indices] = 1
+    return mask
+
+
+def cg_mri_reconstruction(kspace_sampled, mask, num_iters=200):
+    """
+    使用共轭梯度法进行MRI重建
+    :param kspace_sampled: 采样的k空间数据
+    :param mask: 采样掩码
+    :param num_iters: 迭代次数
+    :return: 重建的图像
+    """
+    # 初始化重建图像
+    img_reconstructed = ifft2(kspace_sampled)
+
+    def A(x):
+        # 正向算子：从图像到采样k空间
+        kspace = fft2(x)
+        return kspace * mask
+
+    def At(y):
+        # 伴随算子：从采样k空间到图像
+        return ifft2(y)
+
+    # 初始残差
+    r = kspace_sampled - A(img_reconstructed)
+    p = At(r)
+
+    epsilon = 1e-10  # 一个很小的正数，用于避免除零错误
+
+    for i in range(num_iters):
+        Ap = A(p)
+        denominator_alpha = np.sum(np.conj(Ap) * Ap)
+        if np.abs(denominator_alpha) < epsilon:
+            denominator_alpha = epsilon
+        alpha = np.sum(np.conj(r) * r) / denominator_alpha
+        img_reconstructed = img_reconstructed + alpha * p
+        r_new = r - alpha * Ap
+        denominator_beta = np.sum(np.conj(r) * r)
+        if np.abs(denominator_beta) < epsilon:
+            denominator_beta = epsilon
+        beta = np.sum(np.conj(r_new) * r_new) / denominator_beta
+        p = At(r_new) + beta * p
+        r = r_new
+
+    return np.abs(img_reconstructed)
+
+def main():
+    # 设置当前工作目录
+    current_dir = os.getcwd()
+
+    # 加载完整的k空间数据
+    kspace_path = os.path.join(current_dir, 'fid_20250417_103329.npy')
+    if not os.path.exists(kspace_path):
+        print(f"错误: 找不到文件 {kspace_path}")
+        return
+
+    try:
+        kspace_full = np.load(kspace_path)
+        print(f"成功加载k空间数据: {kspace_path}")
+    except Exception as e:
+        print(f"错误: 加载k空间数据时出错: {e}")
+        return
+
+    # 设置采样率
+    sampling_rate = 1  # 可以根据需要调整
+
+    # 生成采样掩码
+    mask1 = generate_sampling_mask(kspace_full.shape, sampling_rate)
+
+    # 采样k空间
+
+    # 共轭梯度法重建图像
+
+    # 遍历当前目录下的所有 .npy 文件，排除 'mrd_data.npy'
+    for filename in os.listdir(current_dir):
+        if filename.lower().endswith('.npy') and filename != 'fid_20250417_103329.npy':
+            mask_path = os.path.join(current_dir, filename)
+            try:
+                mask = np.load(mask_path)
+                print(f"处理文件: {mask_path}")
+            except Exception as e:
+                print(f"错误: 加载文件 {mask_path} 时出错: {e}")
+                continue
+
+            # 检查 kspace_full 和 mask 的形状是否匹配
+            if kspace_full.shape != mask.shape:
+                print(f"警告: 文件 {filename} 的形状与 k 空间数据不匹配，跳过。")
+                continue
+
+            # 采样 k 空间
+            kspace_sampled = kspace_full * mask1
+
+            # 重建图像
+            try:
+                img_reconstructed = cg_mri_reconstruction(kspace_sampled, mask1, num_iters=100)
+                print(f"成功重建图像: {filename}")
+            except Exception as e:
+                print(f"错误: 重建文件 {filename} 时出错: {e}")
+                continue
+
+            # 构造保存的 PNG 文件名
+            base_name = os.path.splitext(filename)[0]
+            png_filename = f"{base_name}_cg_reconstructed_fid_20250417_103329.png"
+            png_path = os.path.join(current_dir, png_filename)
+            #plt.imshow(img_reconstructed)
+            plt.imshow(img_reconstructed, cmap='gray')  # 设置为灰度图显示
+            plt.axis('off')  # 取消坐标轴
+            plt.tight_layout(pad=0)
+            plt.savefig(png_path, dpi=300)  # 保存为 PNG 文件
+
+
+if __name__ == "__main__":
+    main()

cs.py

+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import nnls
+from scipy.sparse import csr_matrix
+
+
+def read_mrd(file_path, header_size=512, data_shape=(384, 448), data_type=np.complex64):
+    """
+    读取 .mrd 文件并解析为复数浮点型数组。
+
+    :param file_path: .mrd 文件路径
+    :param header_size: 文件头的字节大小
+    :param data_shape: 数据的维度（例如 K 空间或图像空间的形状）
+    :param data_type: 数据类型（默认为 np.complex64,即单精度复数）
+    :return: 解析后的复数数组
+    """
+    # 将二进制数据解析为实数类型（float32 或 float64）
+    if data_type == np.complex64:
+        float_type = np.float32  # 单精度浮点数
+    elif data_type == np.complex128:
+        float_type = np.float64  # 双精度浮点数
+    else:
+        raise ValueError("Unsupported complex data type. Use np.complex64 or np.complex128.")
+    # 确定单个复数占用的字节数
+    bytes_per_sample = np.dtype(data_type).itemsize  # complex64 -> 8 bytes (4 实部 + 4 虚部)
+    float_type_size = np.dtype(float_type).itemsize  # 每个浮点数的字节大小
+
+    # 计算数据部分的字节数
+    total_data_points = np.prod(data_shape) * 2  # 每个复数由两个浮点数组成
+    data_bytes = total_data_points * float_type_size  # 数据部分的总字节数
+
+    with open(file_path, 'rb') as f:
+        # 跳过头部
+        f.seek(header_size)
+
+        # 读取数据部分
+        raw_data = f.read(data_bytes)
+        if len(raw_data) != data_bytes:
+            raise ValueError(f"Error: Expected {data_bytes} bytes of data, but got {len(raw_data)} bytes.")
+
+        # 解析为一维数组（实部和虚部交替存储）
+        real_imag_array = np.frombuffer(raw_data, dtype=float_type)
+
+        # 将实部和虚部合成为复数
+        complex_data = real_imag_array[::2] + 1j * real_imag_array[1::2]
+
+        # 重塑为指定维度
+        complex_data = complex_data.reshape(data_shape)
+    return complex_data
+
+
+def omp(A, y, sparsity_level):
+    """
+    正交匹配追踪算法实现
+    :param A: 测量矩阵（采样矩阵与稀疏矩阵的乘积）
+    :param y: 采样得到的测量信号
+    :param sparsity_level: 信号的稀疏度
+    :return: 恢复的稀疏信号
+    """
+    m, n = A.shape
+    r = y.copy()
+    omega = []
+    x = np.zeros(n)
+    for _ in range(sparsity_level):
+        correlations = np.abs(A.T @ r)
+        lambda_new = np.argmax(correlations)
+        omega.append(lambda_new)
+        A_omega = A[:, omega]
+        A_omega_dense = A_omega.toarray()
+        x_omega, _ = nnls(A_omega_dense, np.real(y))
+        x[omega] = x_omega
+        r = y - A @ x
+    return x
+
+
+# 示例调用
+file_path = r"D:\Users\Desktop\egg\Scan2.mrd"
+data_shape = (384, 448)  # 替换为实际的 K 空间形状
+header_size = 512  # 假设头部大小为 512 字节
+data_type = np.complex64  # 假设是单精度复数（complex64）
+
+# 读取复数数据
+mrd_data = read_mrd(file_path, header_size=header_size, data_shape=data_shape, data_type=data_type)
+np.save('mrd_data1.npy', mrd_data)
+
+# 计算复数数据的幅度（模）
+magnitude = np.abs(mrd_data)
+
+# 可视化原始 K 空间图像
+plt.imshow(magnitude, cmap='gray', aspect='auto')
+plt.title("Original K-Space Image")
+plt.colorbar(label="Normalized Intensity")
+plt.savefig("pic/originalK.png", dpi=300)  # 保存为 PNG 文件
+plt.show()
+
+# 预定义采样矩阵（手动生成稀疏随机矩阵）
+sampling_rate = 0.2  # 采样率
+n = data_shape[0] * data_shape[1]
+m = int(sampling_rate * n)
+density = 0.1  # 稀疏矩阵的非零元素密度
+nnz = int(m * n * density)  # 非零元素数量
+
+# 手动生成随机索引
+rows = np.random.randint(0, m, nnz)
+cols = np.random.randint(0, n, nnz)
+data = np.random.randn(nnz)
+
+# 创建稀疏矩阵
+Phi = csr_matrix((data, (rows, cols)), shape=(m, n))
+
+# 逐列计算稀疏矩阵的范数
+col_norms = np.sqrt(np.asarray(Phi.power(2).sum(axis=0))).flatten()
+# 避免除零错误
+col_norms[col_norms == 0] = 1
+# 逐列归一化稀疏矩阵
+Phi = csr_matrix((Phi.data / col_norms[Phi.indices], Phi.indices, Phi.indptr), shape=Phi.shape)
+
+# 完整的 k 空间数据向量化
+kspace_full = mrd_data.flatten()
+
+# 定义傅里叶变换操作函数
+def fft_operation(x):
+    x_reshaped = x.reshape(data_shape)
+    return np.fft.fft2(x_reshaped) / np.sqrt(n)
+
+# 压缩感知采样
+def phi_psi_product(x):
+    fft_x = fft_operation(x)
+    fft_x_flat = fft_x.flatten()
+    return Phi @ fft_x_flat
+
+y = phi_psi_product(kspace_full)
+
+# 稀疏度
+sparsity_level = 50
+
+# 定义 A 矩阵的作用（避免显式构建）
+def A_operator(x):
+    return phi_psi_product(x)
+
+# 通过迭代算法（OMP）恢复稀疏信号
+def modified_omp(A_operator, y, sparsity_level, n):
+    r = y.copy()
+    omega = []
+    x = np.zeros(n)
+    for _ in range(sparsity_level):
+        def correlation_operator(i):
+            basis = np.zeros(n)
+            basis[i] = 1
+            A_basis = A_operator(basis)
+            return np.abs(np.dot(A_basis.conj(), r))
+        correlations = np.array([correlation_operator(i) for i in range(n)])
+        lambda_new = np.argmax(correlations)
+        omega.append(lambda_new)
+        A_omega = np.array([A_operator(np.eye(n)[i]) for i in omega]).T
+        x_omega, _ = nnls(A_omega, np.real(y))
+        x[omega] = x_omega
+        Ax = np.array([A_operator(x[i] * np.eye(n)[i]) for i in range(n)]).sum(axis=0)
+        r = y - Ax
+    return x
+
+recovered_sparse_signal = modified_omp(A_operator, y, sparsity_level, n)
+
+# 从恢复的稀疏信号得到 k 空间信息
+recovered_k_space_vector = recovered_sparse_signal
+recovered_k_space = recovered_k_space_vector.reshape(data_shape)
+
+# 从 k 空间重建图像
+reconstructed_image = np.abs(np.fft.ifft2(recovered_k_space))
+
+# 可视化重建图像
+plt.imshow(reconstructed_image, cmap='gray')
+plt.colorbar()
+plt.title("Reconstructed Image (Compressed Sensing)")
+plt.savefig("pic/reconstructed_cs.png", dpi=300)  # 保存为 PNG 文件
+plt.show()
+

cs_Fista.py

+import numpy as np
+import os
+import matplotlib.pyplot as plt
+from scipy.fft import fft2, ifft2, fftshift
+
+from scipy.fftpack import dct, idct
+
+
+def dct2(a):
+    """二维离散余弦变换"""
+    return dct(dct(a.T, norm='ortho').T, norm='ortho')
+
+
+def idct2(a):
+    """二维离散余弦逆变换"""
+    return idct(idct(a.T, norm='ortho').T, norm='ortho')
+
+
+def fft2c(img):
+    """Centered 2D FFT."""
+    return np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(img)))
+
+
+def ifft2c(kspace):
+    """Centered 2D IFFT."""
+    return np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(kspace)))
+
+
+def soft_threshold(x, lam):
+    """软阈值操作"""
+    return np.sign(x) * np.maximum(np.abs(x) - lam, 0)
+
+
+def total_variation(x, alpha):
+    """计算图像的总变分"""
+    dx = np.diff(x, axis=0)
+    dy = np.diff(x, axis=1)
+    return alpha * (np.sum(np.abs(dx)) + np.sum(np.abs(dy)))
+
+
+def reconstruct_image(kspace_data):
+    im = fftshift(ifft2(kspace_data, (256, 256)), axes=0)
+    # 取模得到实值图像
+    magnitude_image = np.abs(im)
+    # 归一化图像
+    normalized_image = (magnitude_image - magnitude_image.min()) / (magnitude_image.max() - magnitude_image.min())
+    return normalized_image
+
+
+def fista_mri_reconstruction(kspace_sampled, mask, lam, num_iters=2000, step_size=0.05, tv_alpha=0.05):
+    """
+    使用快速迭代软阈值算法FISTA进行MRI重建并引入总变分约束
+    :param kspace_sampled: 采样的k空间数据
+    :param mask: 采样掩码
+    :param lam: 正则化参数
+    :param num_iters: 迭代次数
+    :param step_size: 步长
+    :param tv_alpha: 总变分约束强度
+    :return: 重建的图像
+    """
+    # 对k空间进行频域阈值处理
+    kspace_sampled_freq = fft2c(kspace_sampled)
+    kspace_sampled_freq_thresholded = soft_threshold(kspace_sampled_freq, 3)
+    kspace_sampled_thresholded = ifft2c(kspace_sampled_freq_thresholded)
+
+    # 初始重建（零填充），使用阈值处理后的k空间数据
+    img_reconstructed = reconstruct_image(kspace_sampled_thresholded)
+
+    # 显示初始化图像
+    plt.imshow(np.abs(img_reconstructed), cmap='gray')
+    plt.title('Initial Reconstructed MRI Image')
+    plt.show()
+    plt.imsave('initial_reconstructed_image.png', np.abs(img_reconstructed), cmap='gray')
+
+    y = img_reconstructed.copy()
+    t = 1
+
+    for i in range(num_iters):
+        # 将当前图像转换到稀疏域（使用二维离散余弦变换 DCT）
+        img_dct = dct2(y)
+        img_dct_thresholded = soft_threshold(img_dct, lam * step_size)
+
+        # 将稀疏系数转换回图像域
+        img_temp = idct2(img_dct_thresholded)
+
+        # 引入总变分约束
+        grad_tv = np.gradient(img_temp)
+        grad_tv = np.sqrt(grad_tv[0] ** 2 + grad_tv[1] ** 2)
+        img_temp = img_temp - step_size * tv_alpha * grad_tv
+
+        # 数据一致性步骤
+        kspace_reconstructed = fft2c(img_temp)
+        kspace_reconstructed = kspace_sampled * mask + kspace_reconstructed * (1 - mask)
+        img_new = reconstruct_image(kspace_reconstructed)
+
+        t_new = (1 + np.sqrt(1 + 4 * t ** 2)) / 2
+        y = img_new + ((t - 1) / t_new) * (img_new - img_reconstructed)
+
+        img_reconstructed = img_new
+        t = t_new
+
+        if (i + 1) % 10 == 0 or i == 0:
+            print(f"Iteration {i + 1}/{num_iters}")
+
+    return img_reconstructed
+
+
+def main(kspace_path, mask_path):
+    # 设置 FISTA 参数
+    lam = 0.05  # 正则化参数，可以根据需要调整
+    num_iters = 2000  # 迭代次数
+    step_size = 0.05  # 步长，可以根据需要调整
+    tv_alpha = 0.1  # 总变分约束强度，可以根据需要调整
+
+    kspace_full = np.load(kspace_path)
+    mask = np.load(mask_path)
+
+    # 使用完整k空间数据重建并显示图像
+    full_image = reconstruct_image(kspace_full)
+    plt.imshow(np.abs(full_image), cmap='gray')
+    plt.title('Full MRI Image from Full K-space')
+    plt.show()
+
+    # 检查kspace_full和mask的形状是否匹配
+    if kspace_full.shape != mask.shape:
+        print(f"警告: 文件 {mask_path} 的形状与k空间数据不匹配,程序退出。")
+        return
+
+    # 采样k空间
+    kspace_sampled = kspace_full * mask
+
+    # 重建图像
+    img_reconstructed = fista_mri_reconstruction(kspace_sampled, mask, lam, num_iters, step_size, tv_alpha)
+
+    # 构造保存的PNG文件名
+    base_name = os.path.splitext(os.path.basename(mask_path))[0]
+    current_dir = os.getcwd()
+    png_filename = f"{base_name}_frequency_fista_reconstructed.png"
+    png_path = os.path.join(current_dir, png_filename)
+
+    plt.imshow(img_reconstructed, cmap='gray')  # 设置为灰度图显示
+    plt.title('Reconstructed MRI Image with Frequency Domain Sparsity')
+    plt.show()
+    plt.imsave(png_path, np.abs(img_reconstructed), cmap='gray')
+
+
+if __name__ == "__main__":
+    # 在这里直接定义路径
+    kspace_path = r'D:\Users\Desktop\egg\fid_20250417_103329.npy'
+    mask_path = r'D:\Users\Desktop\egg\radial_mask.npy'
+    main(kspace_path, mask_path)

cs_ISTA.py

+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.fft import fft2, ifft2, fftshift, ifftshift
+import os
+import pywt
+from scipy.fftpack import dct, idct
+
+
+def fft2c(img):
+    """Centered 2D FFT."""
+    return np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(img)))
+
+
+def ifft2c(kspace):
+    """Centered 2D IFFT."""
+    return np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(kspace)))
+
+
+def generate_sampling_mask(shape, sampling_rate):
+    """生成随机采样掩码"""
+    num_samples = int(sampling_rate * shape[0] * shape[1])
+    mask = np.zeros(shape, dtype=np.float32)
+    indices = np.random.choice(shape[0] * shape[1], num_samples, replace=False)
+    mask.flat[indices] = 1
+    return mask
+
+
+def dct2(a):
+    """二维离散余弦变换"""
+    return dct(dct(a.T, norm='ortho').T, norm='ortho')
+
+
+def soft_threshold(x, lam):
+    """软阈值操作"""
+    return np.sign(x) * np.maximum(np.abs(x) - lam, 0)
+
+
+def idct2(a):
+    """二维离散余弦逆变换"""
+    return idct(idct(a.T, norm='ortho').T, norm='ortho')
+
+
+def reconstruct_image(kspace_data):
+    im = fftshift(ifft2(kspace_data, (256, 256)), axes=0)
+    # 取模得到实值图像
+    magnitude_image = np.abs(im)
+    # 归一化图像
+    normalized_image = (magnitude_image - magnitude_image.min()) / (magnitude_image.max() - magnitude_image.min())
+    return normalized_image
+
+
+def ista_mri_reconstruction(kspace_sampled, mask, lam, num_iters=100, step_size=1.0, wavelet='db4'):
+    # 初始重建（零填充）
+    img_reconstructed = reconstruct_image(kspace_sampled)
+
+    for i in range(num_iters):
+        # 将当前图像转换到稀疏域（使用离散小波变换 DWT）
+        coeffs = pywt.wavedec2(img_reconstructed, wavelet, level=3)
+        new_coeffs = []
+        for c in coeffs:
+            if isinstance(c, tuple):
+                c = tuple([soft_threshold(x, lam * step_size) for x in c])
+            else:
+                c = soft_threshold(c, lam * step_size)
+            new_coeffs.append(c)
+        # 将稀疏系数转换回图像域
+        img_temp = pywt.waverec2(new_coeffs, wavelet)
+        # 数据一致性步骤
+        kspace_reconstructed = fft2c(img_temp)
+        kspace_reconstructed = kspace_sampled * mask + kspace_reconstructed * (1 - mask)
+        img_reconstructed = reconstruct_image(kspace_reconstructed)
+
+        if (i + 1) % 10 == 0 or i == 0:
+            print(f"Iteration {i + 1}/{num_iters}")
+
+    return img_reconstructed
+
+
+def main(kspace_path, mask_path):
+    # 设置 FISTA 参数
+    lam = 0.05  # 正则化参数，可以根据需要调整
+    num_iters = 500  # 迭代次数
+    rho = 2.0  # 惩罚参数
+    step_size = 0.05  # 步长，可以根据需要调整
+    tv_alpha = 0.1  # 总变分约束强度，可以根据需要调整
+    wavelets = ['haar']
+
+    kspace_full = np.load(kspace_path)
+    mask = np.load(mask_path)
+
+    # 调整形状使其匹配
+    min_shape_0 = min(kspace_full.shape[0], mask.shape[0])
+    min_shape_1 = min(kspace_full.shape[1], mask.shape[1])
+    kspace_full = kspace_full[:min_shape_0, :min_shape_1]
+    mask = mask[:min_shape_0, :min_shape_1]
+
+    # 检查kspace_full和mask的形状是否匹配
+    if kspace_full.shape != mask.shape:
+        print(f"警告: 文件 {mask_path} 的形状与k空间数据不匹配,程序退出。")
+        return
+
+    # 采样k空间
+    kspace_sampled = kspace_full * mask
+
+    for wavelet in wavelets:
+        # 重建图像
+        img_reconstructed = ista_mri_reconstruction(kspace_sampled, mask, lam, num_iters, step_size)
+
+        # 构造保存的PNG文件名，包含小波基名称
+        base_name = os.path.splitext(os.path.basename(mask_path))[0]
+        current_dir = os.getcwd()
+        png_filename = f"{base_name}_{wavelet}_ista_reconstructed.png"
+        png_path = os.path.join(current_dir, png_filename)
+        plt.imshow(img_reconstructed, cmap='gray')  # 设置为灰度图显示
+        plt.title(f'Reconstructed MRI Image with {wavelet} wavelet')
+        plt.show()
+        plt.imsave(png_path, np.abs(img_reconstructed), cmap='gray')
+
+
+if __name__ == "__main__":
+    # 在这里直接定义路径
+    kspace_path = r'D:\Users\Desktop\egg\mrd_data1.npy'
+    mask_path = r'D:\Users\Desktop\egg\line_128.npy'
+    main(kspace_path, mask_path)

cs_admm.py

+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.fftpack import dct, idct
+import os
+
+
+def fft2c(img):
+    """Centered 2D FFT."""
+    return np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(img)))
+
+
+def ifft2c(kspace):
+    """Centered 2D IFFT."""
+    return np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(kspace)))
+
+
+def generate_sampling_mask(shape, sampling_rate):
+    """生成随机采样掩码"""
+    num_samples = int(sampling_rate * shape[0] * shape[1])
+    mask = np.zeros(shape, dtype=np.float32)
+    indices = np.random.choice(shape[0] * shape[1], num_samples, replace=False)
+    mask.flat[indices] = 1
+    return mask
+
+
+def dct2(a):
+    """二维离散余弦变换"""
+    return dct(dct(a.T, norm='ortho').T, norm='ortho')
+
+
+def idct2(a):
+    """二维离散余弦逆变换"""
+    return idct(idct(a.T, norm='ortho').T, norm='ortho')
+
+
+def soft_threshold(x, lam):
+    """软阈值操作"""
+    return np.sign(x) * np.maximum(np.abs(x) - lam, 0)
+
+
+def admm_mri_reconstruction(kspace_sampled, mask, lam, rho=1.0, num_iters=100):
+    """
+    使用交替方向乘子法进行 MRI 重建
+    :param kspace_sampled: 采样的 k 空间数据
+    :param mask: 采样掩码
+    :param lam: 正则化参数
+    :param rho: 惩罚参数
+    :param num_iters: 迭代次数
+    :return: 重建的图像
+    """
+    # 初始重建（零填充）
+    img_reconstructed = ifft2c(kspace_sampled)
+    z = dct2(img_reconstructed)
+    u = np.zeros_like(z)
+
+    for i in range(num_iters):
+        # 更新 x
+        kspace_x = kspace_sampled * mask + fft2c(ifft2c(idct2(z - u))) * (1 - mask)
+        img_reconstructed = ifft2c(kspace_x)
+
+        # 更新 z
+        sparse_coeff = dct2(img_reconstructed) + u
+        z = soft_threshold(sparse_coeff, lam / rho)
+
+        # 更新 u
+        u = u + (dct2(img_reconstructed) - z)
+
+        if (i + 1) % 10 == 0 or i == 0:
+            print(f"Iteration {i + 1}/{num_iters}")
+
+    return np.abs(img_reconstructed)
+
+
+def main():
+    # 设置当前工作目录
+    current_dir = os.getcwd()
+
+    # 加载完整的 k 空间数据
+    kspace_path = os.path.join(current_dir, 'fid_20250417_103329.npy')
+    if not os.path.exists(kspace_path):
+        print(f"错误: 找不到文件 {kspace_path}")
+        return
+
+    try:
+        kspace_full = np.load(kspace_path)
+        print(f"成功加载 k 空间数据: {kspace_path}")
+    except Exception as e:
+        print(f"错误: 加载 k 空间数据时出错: {e}")
+        return
+
+    # 设置 ADMM 参数
+    lam = 0.6  # 正则化参数，可以根据需要调整
+    rho = 2.0  # 惩罚参数
+    num_iters = 100 # 迭代次数
+
+    # 遍历当前目录下的所有 .npy 文件，排除 'mrd_data.npy'
+    for filename in os.listdir(current_dir):
+        if filename.lower().endswith('.npy') and filename != 'fid_20250417_103329.npy':
+            mask_path = os.path.join(current_dir, filename)
+            try:
+                mask = np.load(mask_path)
+                print(f"处理文件: {mask_path}")
+            except Exception as e:
+                print(f"错误: 加载文件 {mask_path} 时出错: {e}")
+                continue
+
+            # 检查 kspace_full 和 mask 的形状是否匹配
+            if kspace_full.shape != mask.shape:
+                print(f"警告: 文件 {filename} 的形状与 k 空间数据不匹配，跳过。")
+                continue
+
+            # 采样 k 空间
+            kspace_sampled = kspace_full * mask
+
+            # 重建图像
+            try:
+                img_reconstructed = admm_mri_reconstruction(kspace_sampled, mask, lam, rho, num_iters)
+                print(f"成功重建图像: {filename}")
+            except Exception as e:
+                print(f"错误: 使用 ADMM 重建文件 {filename} 时出错: {e}")
+                continue
+
+            # 构造保存的 PNG 文件名
+            base_name = os.path.splitext(filename)[0]
+            png_filename = f"{base_name}_ADMM_reconstructed_fid_20250417_103329.png"
+            png_path = os.path.join(current_dir, png_filename)
+            #plt.imshow(img_reconstructed)
+            plt.imshow(img_reconstructed, cmap='gray')  # 设置为灰度图显示
+            plt.axis('off')  # 取消坐标轴
+            plt.tight_layout(pad=0)
+            plt.savefig(png_path, dpi=300)  # 保存为 PNG 文件
+
+
+if __name__ == "__main__":
+    main()
+    
\ No newline at end of file

cs_bpnd.py

+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.fftpack import dct, idct
+import os
+from scipy.optimize import minimize
+
+
+def fft2c(img):
+    """Centered 2D FFT."""
+    return np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(img)))
+
+
+def ifft2c(kspace):
+    """Centered 2D IFFT."""
+    return np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(kspace)))
+
+
+def generate_sampling_mask(shape, sampling_rate):
+    """生成随机采样掩码"""
+    num_samples = int(sampling_rate * shape[0] * shape[1])
+    mask = np.zeros(shape, dtype=np.float32)
+    indices = np.random.choice(shape[0] * shape[1], num_samples, replace=False)
+    mask.flat[indices] = 1
+    return mask
+
+
+def dct2(a):
+    """二维离散余弦变换"""
+    return dct(dct(a.T, norm='ortho').T, norm='ortho')
+
+
+def idct2(a):
+    """二维离散余弦逆变换"""
+    return idct(idct(a.T, norm='ortho').T, norm='ortho')
+
+
+def bpdn_mri_reconstruction(kspace_sampled, mask, lam, num_iters=100):
+    """
+    使用基追踪去噪算法进行 MRI 重建
+    :param kspace_sampled: 采样的 k 空间数据
+    :param mask: 采样掩码
+    :param lam: 正则化参数
+    :param num_iters: 最大迭代次数
+    :return: 重建的图像
+    """
+    # 初始猜测
+    img_initial = ifft2c(kspace_sampled)
+
+    def objective(x):
+        """目标函数：数据拟合项 + 正则化项"""
+        img = x.reshape(mask.shape)
+        kspace_reconstructed = fft2c(img)
+        data_fidelity = np.linalg.norm(kspace_sampled * mask - kspace_reconstructed * mask) ** 2
+        sparse_coeff = dct2(img)
+        regularization = lam * np.linalg.norm(sparse_coeff.flatten(), 1)
+        return data_fidelity + regularization
+
+    # 优化求解
+    result = minimize(objective, img_initial.flatten(), method='L-BFGS-B', options={'maxiter': num_iters})
+    img_reconstructed = result.x.reshape(mask.shape)
+    return np.abs(img_reconstructed)
+
+
+# 下采样函数
+def downsample(data, factor):
+    return data[::factor, ::factor]
+
+
+def main():
+    # 设置当前工作目录
+    current_dir = os.getcwd()
+
+    # 加载完整的 k 空间数据
+    kspace_path = os.path.join(current_dir, 'mrd_data1.npy')
+    if not os.path.exists(kspace_path):
+        print(f"错误: 找不到文件 {kspace_path}")
+        return
+
+    try:
+        kspace_full = np.load(kspace_path)
+        # 进一步增大下采样因子，可根据实际情况调整
+        downsample_factor = 4
+        kspace_full = downsample(kspace_full, downsample_factor)
+        print(f"成功加载并下采样 k 空间数据: {kspace_path}")
+    except Exception as e:
+        print(f"错误: 加载 k 空间数据时出错: {e}")
+        return
+
+    # 设置 BPDN 参数
+    lam = 0.1  # 正则化参数，可以根据需要调整
+    num_iters = 100  # 最大迭代次数
+
+    # 遍历当前目录下的所有 .npy 文件，排除 'mrd_data.npy'
+    for filename in os.listdir(current_dir):
+        if filename.lower().endswith('.npy') and filename != 'mrd_data.npy':
+            mask_path = os.path.join(current_dir, filename)
+            try:
+                mask = np.load(mask_path)
+                # 对掩码也进行下采样
+                mask = downsample(mask, downsample_factor)
+                print(f"处理文件: {mask_path}")
+            except Exception as e:
+                print(f"错误: 加载文件 {mask_path} 时出错: {e}")
+                continue
+
+            # 检查 kspace_full 和 mask 的形状是否匹配
+            if kspace_full.shape != mask.shape:
+                print(f"警告: 文件 {filename} 的形状与 k 空间数据不匹配，跳过。")
+                continue
+
+            # 采样 k 空间
+            kspace_sampled = kspace_full * mask
+
+            # 重建图像
+            try:
+                img_reconstructed = bpdn_mri_reconstruction(kspace_sampled, mask, lam, num_iters)
+                print(f"成功重建图像: {filename}")
+            except Exception as e:
+                print(f"错误: 使用 BPDN 重建文件 {filename} 时出错: {e}")
+                continue
+
+            # 构造保存的 PNG 文件名
+            base_name = os.path.splitext(filename)[0]
+            png_filename = f"{base_name}_BPDN_reconstructed.png"
+            png_path = os.path.join(current_dir, png_filename)
+            plt.imshow(img_reconstructed)
+            plt.axis('off')  # 取消坐标轴
+            plt.tight_layout(pad=0)
+            plt.savefig(png_path, dpi=300)  # 保存为 PNG 文件
+
+
+if __name__ == "__main__":
+    main()
+    
+    
\ No newline at end of file

cs_cgm.py

+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.fftpack import dct, idct
+import os
+
+
+def fft2c(img):
+    """Centered 2D FFT."""
+    return np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(img)))
+
+
+def ifft2c(kspace):
+    """Centered 2D IFFT."""
+    return np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(kspace)))
+
+
+def generate_sampling_mask(shape, sampling_rate):
+    """生成随机采样掩码"""
+    num_samples = int(sampling_rate * shape[0] * shape[1])
+    mask = np.zeros(shape, dtype=np.float32)
+    indices = np.random.choice(shape[0] * shape[1], num_samples, replace=False)
+    mask.flat[indices] = 1
+    return mask
+
+
+def dct2(a):
+    """二维离散余弦变换"""
+    return dct(dct(a.T, norm='ortho').T, norm='ortho')
+
+
+def idct2(a):
+    """二维离散余弦逆变换"""
+    return idct(idct(a.T, norm='ortho').T, norm='ortho')
+
+
+def cg_mri_reconstruction(kspace_sampled, mask, num_iters=100):
+    """
+    使用共轭梯度法进行 MRI 重建
+    :param kspace_sampled: 采样的 k 空间数据
+    :param mask: 采样掩码
+    :param num_iters: 迭代次数
+    :return: 重建的图像
+    """
+    # 初始重建（零填充）
+    img_reconstructed = ifft2c(kspace_sampled)
+
+    def A_op(x):
+        """定义线性算子 A"""
+        return fft2c(ifft2c(x) * mask)
+
+    r = kspace_sampled - A_op(img_reconstructed)
+    p = r
+    for i in range(num_iters):
+        Ap = A_op(p)
+        alpha = np.sum(np.conj(r) * r) / np.sum(np.conj(p) * Ap)
+        img_reconstructed = img_reconstructed + alpha * p
+        r_new = r - alpha * Ap
+        beta = np.sum(np.conj(r_new) * r_new) / np.sum(np.conj(r) * r)
+        p = r_new + beta * p
+        r = r_new
+
+        if (i + 1) % 10 == 0 or i == 0:
+            print(f"Iteration {i + 1}/{num_iters}")
+
+    return np.abs(img_reconstructed)
+
+
+def main():
+    # 设置当前工作目录
+    current_dir = os.getcwd()
+
+    # 加载完整的 k 空间数据
+    kspace_path = os.path.join(current_dir, 'mrd_data1.npy')
+    if not os.path.exists(kspace_path):
+        print(f"错误: 找不到文件 {kspace_path}")
+        return
+
+    try:
+        kspace_full = np.load(kspace_path)
+        print(f"成功加载 k 空间数据: {kspace_path}")
+    except Exception as e:
+        print(f"错误: 加载 k 空间数据时出错: {e}")
+        return
+
+    # 设置共轭梯度法参数
+    num_iters = 100  # 迭代次数
+
+    # 遍历当前目录下的所有 .npy 文件，排除 'mrd_data.npy'
+    for filename in os.listdir(current_dir):
+        if filename.lower().endswith('.npy') and filename != 'mrd_data.npy':
+            mask_path = os.path.join(current_dir, filename)
+            try:
+                mask = np.load(mask_path)
+                print(f"处理文件: {mask_path}")
+            except Exception as e:
+                print(f"错误: 加载文件 {mask_path} 时出错: {e}")
+                continue
+
+            # 检查 kspace_full 和 mask 的形状是否匹配
+            if kspace_full.shape != mask.shape:
+                print(f"警告: 文件 {filename} 的形状与 k 空间数据不匹配，跳过。")
+                continue
+
+            # 采样 k 空间
+            kspace_sampled = kspace_full * mask
+
+            # 重建图像
+            try:
+                img_reconstructed = cg_mri_reconstruction(kspace_sampled, mask, num_iters)
+                print(f"成功重建图像: {filename}")
+            except Exception as e:
+                print(f"错误: 使用共轭梯度法重建文件 {filename} 时出错: {e}")
+                continue
+
+            # 构造保存的 PNG 文件名
+            base_name = os.path.splitext(filename)[0]
+            png_filename = f"{base_name}_CMG_reconstructed.png"
+            png_path = os.path.join(current_dir, png_filename)
+            plt.imshow(img_reconstructed)
+            plt.axis('off')  # 取消坐标轴
+            plt.tight_layout(pad=0)
+            plt.savefig(png_path, dpi=300)  # 保存为 PNG 文件
+
+
+if __name__ == "__main__":
+    main()
\ No newline at end of file

cs_ista_k.py

+import numpy as np
+import matplotlib.pyplot as plt
+import pywt
+from scipy.fft import fft2, ifft2, fftshift, ifftshift
+import os
+
+
+def fft2c(img):
+    """Centered 2D FFT."""
+    return np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(img)))
+
+
+def ifft2c(kspace):
+    """Centered 2D IFFT."""
+    return np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(kspace)))
+
+
+def soft_threshold(x, lam):
+    """软阈值操作"""
+    return np.sign(x) * np.maximum(np.abs(x) - lam, 0)
+
+
+def total_variation(x, alpha):
+    """计算图像的总变分"""
+    dx = np.diff(x, axis=0)
+    dy = np.diff(x, axis=1)
+    return alpha * (np.sum(np.abs(dx)) + np.sum(np.abs(dy)))
+
+
+def reconstruct_image(kspace_data):
+    im = fftshift(ifft2(kspace_data, (256, 256)), axes=0)
+    # 取模得到实值图像
+    magnitude_image = np.abs(im)
+    # 归一化图像
+    normalized_image = (magnitude_image - magnitude_image.min()) / (magnitude_image.max() - magnitude_image.min())
+    return normalized_image
+
+
+def ista_mri_reconstruction(kspace_sampled, mask, lam, num_iters=2000, step_size=0.05, tv_alpha=0.05):
+    """
+    使用迭代收缩阈值算法（ISTA）进行MRI重建，并引入总变分约束
+    :param kspace_sampled: 采样的k空间数据
+    :param mask: 采样掩码
+    :param lam: 正则化参数
+    :param num_iters: 迭代次数
+    :param step_size: 步长
+    :param tv_alpha: 总变分约束强度
+    :return: 重建的图像
+    """
+    # 初始重建（零填充）
+    img_reconstructed = reconstruct_image(kspace_sampled)
+
+    # 显示初始化图像
+    plt.imshow(np.abs(img_reconstructed), cmap='gray')
+    plt.title('Initial Reconstructed MRI Image')
+    plt.show()
+    plt.imsave('initial_reconstructed_image.png', np.abs(img_reconstructed), cmap='gray')
+
+    for i in range(num_iters):
+        # 数据一致性步骤
+        kspace_reconstructed = fft2c(img_reconstructed)
+        kspace_reconstructed = kspace_sampled * mask + kspace_reconstructed * (1 - mask)
+        img_temp = reconstruct_image(kspace_reconstructed)
+
+        # 将当前图像转换到稀疏域（使用离散小波变换 DWT）
+        coeffs = pywt.wavedec2(img_temp, 'db4', level=3)
+        new_coeffs = []
+        for c in coeffs:
+            if isinstance(c, tuple):
+                c = tuple([soft_threshold(x, lam * step_size) for x in c])
+            else:
+                c = soft_threshold(c, lam * step_size)
+            new_coeffs.append(c)
+
+        # 将稀疏系数转换回图像域
+        img_reconstructed = pywt.waverec2(new_coeffs, 'db4')
+
+        # 引入总变分约束
+        grad_tv = np.gradient(img_reconstructed)
+        grad_tv = np.sqrt(grad_tv[0] ** 2 + grad_tv[1] ** 2)
+        img_reconstructed = img_reconstructed - step_size * tv_alpha * grad_tv
+
+        if (i + 1) % 10 == 0 or i == 0:
+            print(f"Iteration {i + 1}/{num_iters}")
+
+    return img_reconstructed
+
+
+def main():
+    # 设置当前工作目录
+    current_dir = os.getcwd()
+
+    # 加载完整的k空间数据
+    kspace_path = os.path.join(current_dir, 'fid_20250418_163142.npy')
+    if not os.path.exists(kspace_path):
+        print(f"错误: 找不到文件 {kspace_path}")
+        return
+
+    try:
+        kspace_full = np.load(kspace_path)
+        print(f"成功加载k空间数据: {kspace_path}")
+    except Exception as e:
+        print(f"错误: 加载k空间数据时出错: {e}")
+        return
+
+    # 设置 ISTA 参数
+    lam = 0.05  # 正则化参数，可以根据需要调整
+    num_iters = 10000  # 迭代次数
+    step_size = 0.05  # 步长，可以根据需要调整
+    tv_alpha = 0.1  # 总变分约束强度，可以根据需要调整
+
+    # 遍历当前目录下的所有.npy文件，排除'mrd_data.npy'
+    for filename in os.listdir(current_dir):
+        if filename.lower().endswith('.npy') and filename!= 'fid_20250418_163142.npy':
+            mask_path = os.path.join(current_dir, filename)
+            try:
+                mask = np.load(mask_path)
+                print(f"处理文件: {mask_path}")
+            except Exception as e:
+                print(f"错误: 加载文件 {mask_path} 时出错: {e}")
+                continue
+
+            # 检查kspace_full和mask的形状是否匹配
+            if kspace_full.shape!= mask.shape:
+                print(f"警告: 文件 {filename} 的形状与k空间数据不匹配，跳过。")
+                continue
+
+            # 采样k空间
+            kspace_sampled = kspace_full * mask
+
+            # 重建图像
+            try:
+                img_reconstructed = ista_mri_reconstruction(kspace_sampled, mask, lam, num_iters, step_size, tv_alpha)
+                print(f"成功重建图像: {filename}")
+            except Exception as e:
+                print(f"错误: 重建文件 {filename} 时出错: {e}")
+                continue
+
+            # 构造保存的PNG文件名
+            base_name = os.path.splitext(filename)[0]
+            png_filename = f"{base_name}_ista_reconstructed_fid_20250417_103329.png"
+            png_path = os.path.join(current_dir, png_filename)
+            plt.imshow(img_reconstructed, cmap='gray')  # 设置为灰度图显示
+            plt.title('Reconstructed MRI Image')
+            plt.show()
+            plt.imsave(png_path, np.abs(img_reconstructed), cmap='gray')
+
+
+if __name__ == "__main__":
+    main()