Added spin-echo T2star illustrations and Some bug fixes

Contrast/t2star_spinecho_illustration.m

@@ -0,0 +1,96 @@
+write_animations = 1;
+spinecho_flag = 1;
+% Simulations of off-resonance and T2* decay
+
+B0 = 1.5e3; % 1500 mT = 1.5 T
+% using units of mT and ms for the bloch_rotate function
+T2 = 80; % ms
+T1 = Inf; % neglect T1
+% create range of off-resonance in a voxel
+Nvoxel = 8^2;
+off_resonance_ranges = [0 0.07 0.25]*1e-6; % off resonance range, in ppm
+off_resonance_description = {'no' 'mild' 'severe'};
+
+for I = 1:3
+
+    imall = [];
+    off_resonance_range = off_resonance_ranges(I); % off resonance range, in ppm
+
+dBZ = linspace(-off_resonance_range,off_resonance_range, Nvoxel); % assumes a uniform distribution
+dBZ = (rand(Nvoxel) -0.5 )* off_resonance_range;  % Uniform distribution
+dBZ = randn(Nvoxel) * off_resonance_range;  % Gaussian distribution
+
+
+% start after RF excitation
+Mstart = [1,0,0].'; 
+
+Tmax = T2;
+dt = .5; % ms
+N = Tmax/dt;
+t = [1:N]*dt;
+Mall = zeros(3,N,Nvoxel);
+Mall(:,1,:) = repmat(Mstart, [1 1 Nvoxel]);
+for INvoxel = 1:Nvoxel
+    Bvoxel = [0 0 dBZ(INvoxel)*B0]; % off-resonance, in rotating frame
+    for It = 1:N-1
+        if spinecho_flag && It == round(N/3)
+            Mall(2,It,INvoxel) = -Mall(2,It,INvoxel); % 180-pulse to invert My
+        end
+        Mall(:,It+1,INvoxel) = bloch_relax(  bloch_rotate(Mall(:,It,INvoxel),dt,Bvoxel), dt, 1, T1, T2 );
+    end
+end
+
+Mall_voxel = sum(Mall, 3)/Nvoxel;
+
+[x y] = meshgrid(1:sqrt(Nvoxel));
+Splot = 0.5;
+fs = figure(1);
+%%
+for It = 1:N
+    Mxy_plot = Mall(1,It,:) + 1i*Mall(2,It,:);
+    Mxy_plot = reshape(Mxy_plot, [sqrt(Nvoxel) sqrt(Nvoxel)]);
+
+    subplot(1,3,1)
+    quiver(x-real(Mxy_plot)/2,y-imag(Mxy_plot)/2,...
+            real(Mxy_plot), imag(Mxy_plot), ...
+            0)
+    axis([0 9 0 9]), axis square
+    title('M_{XY} within a voxel')
+
+    subplot(1,3,2)
+    Mxy_voxel = sum(sum(Mxy_plot))/length(Mxy_plot(:));
+    quiver(-real(Mxy_voxel)/2,-imag(Mxy_voxel)/2,...
+            real(Mxy_voxel), imag(Mxy_voxel), ...
+            0)
+    axis([-0.5 0.5 -0.5 0.5]), axis square
+    title('M_{XY} summed across a voxel ')
+
+    subplot(1,3,3)
+    plot(t(1:It), Mall_voxel(1,1:It), t(1:It), Mall_voxel(2,1:It),'--')
+    xlim([0 max(t)])
+    ylim([-0.1 1])
+    title('Voxel signal')
+    xlabel('time [ms]')
+
+    drawnow
+
+    currFrame = getframe(fs);
+    if It == 1 && isempty(imall)
+        im = frame2im(currFrame);
+        [imind,map] = rgb2ind(im,256,'nodither');
+    end
+    imall = cat(4, imall, rgb2ind(frame2im(currFrame),map,'nodither'));
+
+end
+
+
+if write_animations
+    if spinecho_flag
+        root_name = 't2star-spinecho_';
+    else
+        root_name = 't2star_';
+    end
+    imwrite(imall,map,[root_name off_resonance_description{I} '_off-resonance.gif'],'DelayTime',0,'LoopCount',Inf) %g443800
+end
+
+end

Notebooks/Accelerated Imaging Methods.ipynb

@@ -45,21 +45,6 @@
     "    * Reconstruct an image from undersampled raw data\n"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Formulating the Reconstruction Problem\n",
-    "\n",
-    "For image reconstruction, it is helpful to reformulate the MRI Signal Equation, which we can pose for MRI as a linear system and using standard mathematical notation describing systems as:\n",
-    "\n",
-    "$$ \\mathbf{y} = \\mathbf{Ex} + \\mathbf{n} $$\n",
-    "\n",
-    "where $\\mathbf{y}$ is the acquired data, $\\mathbf{E}$ is the encoding matrix, $\\mathbf{x}$ is the spatial distribution of the transverse magnetization (e.g. image), and $\\mathbf{n}$ is noise\n",
-    "\n",
-    "The encoding matrix, $E$, must include a discrete Fourier Transform matrix, representing our data is the Fourier Transform of the transverse magnetization, evaluated at k-space locations.  It should also include coil sensitivity profiles in order to support parallel imaging formulations."
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -70,9 +55,7 @@
     "\n",
     "$$ \\mathbf{y} = \\mathbf{Ex} + \\mathbf{n} $$\n",
     "\n",
-    "In 2D Cartesian\n",
-    "\n",
-    "vectorized data and images\n",
+    "where $\\mathbf{y}$ is the acquired data, $\\mathbf{E}$ is the encoding matrix, $\\mathbf{x}$ is the spatial distribution of the transverse magnetization (e.g. image), and $\\mathbf{n}$ is noise.  In this formulation, the image is vectorized.  For example, 2D FT sampled data and the corresponding 2D image would be converted to:\n",
     "\n",
     "$$ y = \\left[ \n",
     "    \\begin{array}{c}\n",
@@ -100,9 +83,7 @@
     "\\right]$$\n",
     "\n",
     "\n",
-    "where $\\mathbf{y}$ is the acquired data, $\\mathbf{E}$ is the encoding matrix, $\\mathbf{x}$ is the spatial distribution of the transverse magnetization (e.g. image), and $\\mathbf{n}$ is noise\n",
-    "\n",
-    "The encoding matrix, $E$, must include a discrete Fourier Transform matrix, representing our data is the Fourier Transform of the transverse magnetization, evaluated at k-space locations.  It should also include coil sensitivity profiles in order to support parallel imaging formulations."
+    "The encoding matrix, $E$, at a minimum includes a discrete Fourier Transform matrix, $\\mathbf{F}$, representing our data is the Fourier Transform of the transverse magnetization, evaluated at k-space locations. "
    ]
   },
   {
@@ -111,9 +92,19 @@
    "source": [
     "## Parallel Imaging\n",
     "\n",
-    "For parallel imaging (PI), we need to consider the coil sensitivity profiles, $\\mathbf{C}$, into encoding matrix along with a Fourier Transform encoding matrix, $\\mathbf{F}$:\n",
+    "For parallel imaging (PI), we need to consider the coil sensitivity profiles, $\\mathbf{C}_i$, for each RF coil into encoding matrix along with a Fourier Transform encoding matrix, $\\mathbf{F}$, as well as a k-space sampling operator, $\\mathbf{S}$, for the measurements from each RF coil, $\\mathbf{y}_i$:\n",
+    "\n",
+    "$$\\mathbf{y}_i = \\mathbf{E}_i \\mathbf{x} + \\mathbf{n}_i = \\mathbf{S} \\mathbf{F} \\mathbf{C_i} \\mathbf{x} + \\mathbf{n}_i $$\n",
+    "\n",
+    "Here we can return to our original formulation by concatenating the coil dimension, for example as:\n",
+    "\n",
+    "$$\\mathbf{y} = [\\mathbf{y}_1 \\ \\mathbf{y}_2 \\ldots \\mathbf{y}_N]$$ \n",
+    "\n",
+    "$$\\mathbf{E} = [\\mathbf{E}_1 \\ \\mathbf{E}_2 \\ldots \\mathbf{E}_N]$$ \n",
+    "$$\\mathbf{y} = [\\mathbf{n}_1 \\ \\mathbf{n}_2 \\ldots \\mathbf{n}_N]$$ \n",
+    "\n",
+    "$$ \\mathbf{y} = \\mathbf{Ex} + \\mathbf{n} $$\n",
     "\n",
-    "$$\\mathbf{E} = \\mathbf{F} \\mathbf{C}$$\n",
     "\n",
     "### Image-space Methods\n",
     "Image-space parallel imaging methods (e.g. SENSE) can be formulated as the following optimization problem\n",
@@ -180,27 +171,40 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Neural Network/Deep Learning Reconstructions\n",
+    "## Machine Learning Reconstructions\n",
     "\n",
-    "Convolutional neural networks (CNNs), which form the backbone of deep learning (DL), can be used to convert k-space to image data from a subset of data samples as well.  Conceptually, these methods can be trained to learn how to incorporate coil sensitivity information like parallel imaging and typical image sparsity patterns like compressed sensing.  Since they learn from real-world data, they can learn information that is the most relevant to MRI data and thus have been shown to support higher acceleration factors. \n",
+    "Machine learning using convolutional neural networks (CNNs), which form the backbone of deep learning (DL), can be used to convert k-space to image data from a subset of data samples as well.  Conceptually, these methods can be trained to learn how to incorporate coil sensitivity information like parallel imaging and typical image sparsity patterns like compressed sensing.  Since they learn from real-world data, they can learn information that is the most relevant to MRI data and thus have been shown to support higher acceleration factors. \n",
     "\n",
-    "Perhaps the most popular class of DL reconstruction methods for MRI are so-called \"un-rolled\" networks.  This term confirms \n",
+    "These methods can often be considered general solutions to a regularized least-squares objective\n",
     "\n",
-    "Compressed Sensing is formulated as the following optimization problem\n",
+    "$$\\hat{x}_{reg} = \\arg \\min_\\mathbf{x} \\frac{1}{2} \\| \\mathbf{y} - \\mathbf{Ex} \\|^2_2 + \\mathcal{R}(\\mathbf{x}) $$\n",
     "\n",
-    "$$\\hat{x}_{CS} = \\arg \\min_\\mathbf{x} \\frac{1}{2} \\| \\mathbf{y} = \\mathbf{Ex} \\|^2_2 +\\tau \\| \\mathbf{Wx} \\|_1 $$\n",
+    "where $\\mathcal{R}(\\cdot)$ can be a parallel imaging or compressed sensing regularizer as above, but can also be implicitly implemented via machine learning techniques.  This means that the machine learning reconstruction is attempting to constrain or regularize the reconstruction.\n",
     "\n",
-    "which includes a data consistency term where the data multiplied by the encoding matrix must match the reconstructed image, and a regularization term that enforces that the image is sparse in some other domain through the sparsifying transform, $\\mathbf{W}$.\n",
+    "Perhaps the most popular class of DL reconstruction methods for MRI are so-called \"un-rolled\" networks.  This term comes from the fact these networks are designed to mimic iterations of classical algorithms to solve optimization problems, and each iteration has been \"un-rolled\" into a series of cascaded neural networks.  Some popular methods include MoDL and Variational Networks.\n",
     "\n",
-    "Popular sparsifying transforms include total variation (TV), total generalized variation (TGV), and wavelets.\n",
     "\n",
-    "The k-space sampling patterns used for these methods typically use pseudo-random undersampling with a variable density that preferentially increases the number of samples near the center of k-space\n",
+    "### SNR in ML Reconstructions\n",
+    "It is difficult to define SNR when using machine learning reconstruction methods, as they inherently perform some denoising when constraining the reconstruction.\n",
     "\n",
     "\n",
-    "### SNR in Compressed Sensing\n",
+    "### Artifacts with ML Reconstructions\n",
+    "\n",
+    "A challenge with machine learning reconstructions is that artifacts can be difficult to identify.  They can result in an over-smoothed appearance.  These methods can also overfit to the learned anatomy and data, meaning features might be erased or hallucinated, although using physics-based techniques provide some assurances of data-consistency.\n",
     "\n",
-    "### Artifacts with Compressed Sensing\n"
+    "### Reference\n",
+    "\n",
+    "For a recent, comprehensive reference on thes methods I reccommend:\n",
+    "\n",
+    "Hammernik K, Küstner T, Yaman B, Huang Z, Rueckert D, Knoll F, Akçakaya M. Physics-Driven Deep Learning for Computational Magnetic Resonance Imaging: Combining physics and machine learning for improved medical imaging. IEEE Signal Process Mag. 2023 Jan;40(1):98-114. doi: 10.1109/msp.2022.3215288. Epub 2023 Jan 2. PMID: 37304755; PMCID: PMC10249732.\n",
+    "\n",
+    "https://arxiv.org/pdf/2203.12215.pdf\n"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": []
   }
  ],
  "metadata": {

Notebooks/Artifacts.ipynb

@@ -2,8 +2,12 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 2,
-   "metadata": {},
+   "execution_count": null,
+   "metadata": {
+    "vscode": {
+     "languageId": "octave"
+    }
+   },
    "outputs": [
     {
      "name": "stdout",
@@ -26,7 +30,7 @@
    "source": [
     "# Artifacts\n",
     "\n",
-    "This notebook includes a simulation of various MRI artifacts, as well as a high-level [Artifact Comparison](#Artifact-Comparison) below.  The wikipedia entry https://en.wikipedia.org/wiki/MRI_artifact is also very comprehensive.\n",
+    "This notebook includes a simulation of various MRI artifacts, as well as a high-level Artifact Comparison below.  The wikipedia entry https://en.wikipedia.org/wiki/MRI_artifact is also very comprehensive.\n",
     "\n",
     "## Learning Goals\n",
     "\n",
@@ -64,8 +68,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
-   "metadata": {},
+   "execution_count": null,
+   "metadata": {
+    "vscode": {
+     "languageId": "octave"
+    }
+   },
    "outputs": [
     {
      "name": "stdout",