ENH: Use the attenuation distance to plot the DWI signal

Use the attenuation distance to plot the DWI signal.
nipreps · Dec 22, 2024 · 7047f9f · 7047f9f
1 parent 28a8dfe
commit 7047f9f
Showing 1 changed file with 93 additions and 7 deletions.
diff --git a/docs/notebooks/dwi_simulated_gp.ipynb b/docs/notebooks/dwi_simulated_gp.ipynb
@@ -65,9 +65,7 @@
     "        gtab, mevals, S0=S0, angles=angles, fractions=fractions, snr=snr\n",
     "    )\n",
     "\n",
-    "    grad = np.vstack([gtab.bvecs.T, gtab.bvals])\n",
-    "\n",
-    "    return signal, sticks, grad"
+    "    return signal, sticks, gtab"
    ],
    "id": "a0f5bab019855954",
    "outputs": [],
@@ -84,12 +82,96 @@
    "cell_type": "code",
    "source": [
     "hsph_dirs = 90\n",
-    "signal, sticks, grad = create_multitensor_dmri_signal(hsph_dirs)"
+    "signal, sticks, gtab = create_multitensor_dmri_signal(hsph_dirs)"
    ],
    "id": "e9545781fe5cf3b8",
    "outputs": [],
    "execution_count": null
   },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": "Define the method that computes the signal attenuation.",
+   "id": "40b8ba2f818a49f8"
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "def diffusion_distance(bvecs, Q):\n",
+    "    # Borrowed from https://github.com/arokem/osmosis/blob/72d848ddc172ce6a72baba4c0d4cf3675db1a381/osmosis/tensor.py#L233-L248\n",
+    "    # sphADC = np.dot(np.dot(bvecs.T, np.matrix(Q).getI()), bvecs)\n",
+    "    sphADC = np.dot(np.dot(bvecs.T, np.linalg.inv(Q)), bvecs)\n",
+    "    dist = np.diag(1 / np.sqrt(sphADC))\n",
+    "\n",
+    "    return dist"
+   ],
+   "id": "d40db09fa860c25e",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": "Perform the DTI fit on the signal and compute the signal attenuation along each diffusion-encoding gradient direction.",
+   "id": "25f32e883f520b86"
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "import dipy.reconst.dti as dti\n",
+    "\n",
+    "dti_model = dti.TensorModel(gtab)\n",
+    "dti_fit = dti_model.fit(signal)\n",
+    "\n",
+    "distance = diffusion_distance(gtab.bvecs.T, dti_fit.quadratic_form)"
+   ],
+   "id": "ba6f132d088a1c0a",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "markdown",
+   "source": "We now define a function to plot the DWI signal.",
+   "id": "e9c2295d90cfbfc1"
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "from matplotlib import pyplot as plt \n",
+    "%matplotlib inline\n",
+    "\n",
+    "def plot_dwi(dwi_data, bvecs, b0s_mask, voxel_idx, b0_idx, distance):\n",
+    "\n",
+    "    s0_voxel = dwi_data[voxel_idx[0], voxel_idx[1], voxel_idx[2], b0_idx]\n",
+    "    s_voxel = dwi_data[voxel_idx[0], voxel_idx[1], voxel_idx[2], ~b0s_mask]\n",
+    "\n",
+    "    # Scale gradient vector values with the distance\n",
+    "    # s_normal = s_voxel/s0_voxel\n",
+    "    # scaled_vectors = bvecs[:, ~b0s_mask] / s_normal[np.newaxis, :]\n",
+    "    scaled_vectors = bvecs[:, ~b0s_mask] / distance[~b0s_mask]\n",
+    "\n",
+    "    # Plot the DWI signal as a 3D point cloud\n",
+    "    fig = plt.figure()\n",
+    "    ax = fig.add_subplot(111, projection=\"3d\")\n",
+    "\n",
+    "    _ = ax.scatter(scaled_vectors[0, :], scaled_vectors[1, :], scaled_vectors[2, :], c=\"blue\", marker=\"*\")\n",
+    "\n",
+    "    # Set labels\n",
+    "    ax.set_xlabel(\"X\")\n",
+    "    ax.set_ylabel(\"Y\")\n",
+    "    ax.set_zlabel(\"Z\")\n",
+    "    ax.set_title(\"Normalized dMRI signal\")\n",
+    "\n",
+    "    return fig"
+   ],
+   "id": "1ed5c04d4b49d8d9",
+   "outputs": [],
+   "execution_count": null
+  },
   {
    "metadata": {},
    "cell_type": "markdown",
@@ -99,13 +181,17 @@
   {
    "metadata": {},
    "cell_type": "code",
-   "execution_count": null,
    "source": [
     "voxel_idx = [0, 0, 0]\n",
-    "dwi_data = signal[np.newaxis, np.newaxis, np.newaxis, :]"
+    "dwi_data = signal[np.newaxis, np.newaxis, np.newaxis, :]\n",
+    "# There is only one b0 volume in the simulated signal\n",
+    "b0_idx = 0\n",
+    "_ = plot_dwi(dwi_data, gtab.bvecs.T, gtab.b0s_mask, voxel_idx, b0_idx, distance)\n",
+    "plt.show()"
    ],
    "id": "c07f103d9bd347cc",
-   "outputs": []
+   "outputs": [],
+   "execution_count": null
   },
   {
    "metadata": {},