MOLLI Simulator

发表于 2023-01-11 更新于 2024-11-07 分类于 MRI 评论数：本文字数： 3.2k 阅读时长 ≈ 5 分钟

这篇是准备复现 MOLLI 的基本算法，基于 Bloch Simulator. 这个算法看起来已经不是特别老了，但我暂时也没去借鉴别人的……这里先用 Matlab 进行最基本的模拟。我想过程应该也不是特别复杂？就先做一下试试看吧。

参考文献：

[1] Messroghli, D. R., Radjenovic, A., Kozerke, S., Higgins, D. M., Sivananthan, M. U., & Ridgway, J. P. (2004). Modified Look-Locker inversion recovery (MOLLI) for high-resolutionT1 mapping of the heart. Magnetic Resonance in Medicine, 52(1), 141–146.

[2] Shao, J., Rapacchi, S., Nguyen, K.-L., & Hu, P. (2016). Myocardial T1 mapping at 3.0 tesla using an inversion recovery spoiled gradient echo readout and bloch equation simulation with slice profile correction (BLESSPC) T1 estimation algorithm: FLASH-MOLLI with BLESSPC T1 Mapping. Journal of Magnetic Resonance Imaging, 43(2), 414–425. https://doi.org/10.1002/jmri.24999

SSFP

这篇文章首先是提到自己是利用 SSFP 的序列来进行图像采集的。

We acquired images in late diastole using a single-shot steady-state free-precession (SSFP) technique, combined with sensitivity encoding to achieve a data acquisition window of <200 ms duration.

这个东西我们之前做过了……虽然那个版本的考虑有点太简单了。

function [Msig,Mss] = sssignal(flip,T1,T2,TE,TR,dfreq)

Rflip = yrot(flip);
[Atr,Btr] = freeprecess(TR-TE,T1,T2,dfreq);
[Ate,Bte] = freeprecess(TE,T1,T2,dfreq);
Mss = inv(eye(3)-Ate*Rflip*Atr) * (Ate*Rflip*Btr+Bte);
Msig = Mss(1)+1i*Mss(2);

程序就是这么个程序。这里的 FA 应该是180，因为是 inversion pulse.

问题分析

模拟信号

在我们之前的 MRI simulation 中，我们都没有考虑任何的梯度磁场。因为我们这里的 dfreq 是默认一维的（只和不同的 z 坐标有关）。理想的情况是：假设一个空间，在这个空间中每一个点都具有一个关于它所在位置所得到的信息，然后在打开各个梯度线圈的时候分别去计算它们的的情况。但是就我们之前的讨论看下来，我们的计算全部都是只考虑的，包括 MOLLI 本身也只考虑的问题。

我们到底是否需要把整个 Sequence 的情况都考虑进去呢？我其实是不太明白的……应该是不用的吧。

有其他文章提到，每次打 inversion pulse 的时候，会存在一个误差的情况（信号损失？等问题）……但这种问题可以稍后再讨论。

总之可以先模拟只存在的情况。

但是存在更大的问题：在模拟信号中，我们的源数据从哪来？直至目前，我们的数据全部都是假定，然后在给定的序列下进行模拟，观察其的情况。但是实际操作的时候，我们使用的是人体在MRI机器下的扫描，势必要问扫描出来的数据中是否能够给出每一个 Slice 下每一个 pixel / voxel 中的每一个质子的？我们得到的数据都是图片信息，那么图片信息撑死也就每个 pixel 有一个灰度值。必须要有三通道的图片信息（譬如 256*256*3）的图片才能够解决这个问题。

T1 mapping

我通过一个例子 T1 Mapping (stanford.edu) 发现，T1 mapping 和我们的模拟信号根本不是一回事。比如在他的源代码中，他分为两个部分：模拟部分与实际数据应用的部分。但是这两个部分是完全分割开的：模拟了半天，实际扫描MRI数据的部分中并没有用到任何模拟的地方。

注：在这个例子中，我发现它的最优化过程中直接将弛豫后的矢量模作为最优化目标函数：直接取这个模最大的那个数，扩充这个数左右区间（比如最大值为280，则去寻找270-290区间，最小刻度为0.1），并进一步寻找最大模。这个目标优化的策略有道理吗？

数据模拟

首先，在仅有的情况下，我们默认所有的（这里默认每个 pixel 上对应一个）都是数值向上的，即只有的值，并且这个。

然后在 inverse pulse 进来之后，相当于形成一个的旋转，所以调用 yrot(math.pi)可以得到一个旋转矩阵。

旋转之后进一步按照下面这个式子进行弛豫：

$M_{z^\prime}\left(t\right) = M_z^0\left(1 - e^{- t / T_1}\right) + M_{z^\prime}\left(0_+\right)e^{- t / T_1}$

弛豫的时间范围就是在各个之间进行。

具体来说，问题有很多：

在这里的，，都是指什么？是哪些 pulse、哪些量之间的时间间隔？
BLESSPC[2] 中考虑到的 subslice 是有考虑各自不同的 flip angle 的。但是我认知的 flip angle 好像只是 pulse 中的那个量呀？
在 T1 mapping 中取最大模的优化目标有意义吗？有没有什么依据？
我们似乎需要对着下面这个式子进行拟合：
$y = A - Be^{-t/T_1^\star}$
那这个式子是怎么来的？为什么要按照这个式子去拟合？
就算我们默认这个式子是已经得到的，那我们是要如何拟合？我们已经知道有不同的，不同的，但是拟合的标准是什么？
是否需要考虑 inhomogeneity？即假设是非稳定的，加入的量？（我不认为这里和这个问题有什么关系，所以这个问题不一定有意义）

关于 3 中提到的式子， [1] 似乎提到了一些解释，但不是很清楚：

y denotes signal intensity, corresponds to the apparent, modified in an LL experiment.

那我的理解姑且是这样的：因为通过数据模拟，我们可以从一张图片上模拟得到某段时间内所有的的值。然后根据这个，我们拟合出 $，，$ ，最后根据下面这个式子去算：

$T_1 = T_1^\star\left(\left(B/A\right) - 1\right)$

如果只是对每个 pixel 进行估计，那么 $，，$ 就都是数字。

2023-4-5 updated

到目前为止，我（在同伴的激励下）完成了这里基本所有的内容。但是由于本文只是单纯地讲到 MOLLI 的操作细节，这里就不介绍数据读取的梯度实现。上面所考虑到的这些都只是基于 inversion pulse 所带来的数据拟合。这些数据如果只是直接从仿真中读取的花，拟合结果是相当好的。但如果是从仿真后的 kspace -> ifft 重建读取，那结果可能就要差很多。
我目前所做的内容：

使用FLASH (Gradient Spoiling) + MOLLI 读取仿真的直接结果，拟合后获得不同值的仿真误差百分比（）——大概是这个范围？
使用b-SSFP + MOLLI 做仿真+梯度模拟+数据读取与建立，获取 k-space 图像后进行 ifft 变换，对获得的数据进行 pixel-wise 的拟合。（目前正在缓慢的程序运行中）

实验目前还在进行中。不知道有没有时间能够把之前做过的这些内容做一个经验+技术总结…如果可以的话，把遇到的问题全部记录一下。