Added in B1+ map visualisation. (#72)

* Added in B1+ map visualisation. Added B1+ map visualisation: Solves #44 Adding the color bar messes with the vertical spacing a bit, but I have left it as is for the sake of generating it all in one figure. Everything else has been left untocuched * Changed title We should not refer to a specific number in the paper (in case the number change) Also, I prefer to refer to B1+ maps as in the paper (while we call the procedure RF shimming) * Cosmetic changes * Fixed wrong use of logical test vs. assignment * Show the mask that was used for RF shimming Instead of the SC mask. See: #72 (review) * Do not threshold the mask * Better color dynamic * Added legend on colorbar * Show mask edge for better clarity, refactoring Separated the code that generates the mask into another cell so the figure generation cell runs much faster * adjusted colorbar, font size * Show more S-I extent --------- Co-authored-by: jcohenadad <[email protected]>
shimming-toolbox · Jan 23, 2024 · 04276c0 · 04276c0
1 parent bf4cf2f
commit 04276c0
Showing 1 changed file with 110 additions and 3 deletions.
diff --git a/data_processing.ipynb b/data_processing.ipynb
@@ -689,13 +689,120 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "a916a6f7",
+   "id": "1dde4c00",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Visualize RF maps to generate Figure 2\n",
+    "# Prepare RF shimming mask for figure\n",
+    "\n",
+    "# Select subject to show\n",
+    "subject = 'sub-05'\n",
+    "os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
+    "\n",
+    "# Create RF shimming mask from the segmentation (re-create the live RF shimming procedure)\n",
+    "file_mask = f\"{subject}_acq-anatCP_TB1TFL_mask-shimming.nii.gz\"\n",
+    "!sct_maths -i {subject}_acq-anatCP_TB1TFL_seg_labeled.nii.gz -thr 3 -uthr 9 -o {file_mask}\n",
+    "!sct_maths -i {file_mask} -bin 1 -o {file_mask}\n",
+    "!sct_create_mask -i {subject}_acq-anatCP_TB1TFL.nii.gz -p centerline,{file_mask} -size 28mm -f cylinder -o {file_mask}\n",
+    "\n",
+    "# Dilate mask\n",
+    "!sct_maths -i {file_mask} -dilate 2 -shape disk -dim 2 -o {subject}_acq-anatCP_TB1TFL_mask-shimming_dil.nii.gz\n",
+    "# Subtract dilated mask with mask to get the edge\n",
+    "!sct_maths -i {subject}_acq-anatCP_TB1TFL_mask-shimming_dil.nii.gz -sub {file_mask} -o {file_mask}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "93d72ec3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create figure of B1+ maps\n",
+    "\n",
+    "# Select subject to show\n",
+    "subject = 'sub-05'\n",
+    "\n",
+    "# Defining crop limits for resulting figure\n",
+    "xmin = 20\n",
+    "xmax = 110\n",
+    "ymin = 20\n",
+    "ymax = 150\n",
+    "\n",
+    "# Defining dynamic range\n",
+    "dynmin = 5 \n",
+    "dynmax = 30\n",
+    "\n",
+    "# Create a figure with multiple subplots\n",
+    "fig, axes = plt.subplots(2, 4, figsize=(8, 6))\n",
+    "font_size = 14\n",
+    "axes=axes.flatten()\n",
     "\n",
-    "# TODO"
+    "# Plotter loop (avoiding the generation of an ungly 4D data structure)\n",
+    "for i,shim_mode in enumerate(shim_modes):\n",
+    "    if i==0: # Grabbing the CP mode anat to use for displaying the SC and cropping noise for the B1+ map\n",
+    "        \n",
+    "        # Load data\n",
+    "        CP_anat=nib.load(f\"{subject}_acq-anat{shim_mode}_TB1TFL.nii.gz\")\n",
+    "        CP_SC=nib.load(file_mask)\n",
+    "        CP_nTpV=nib.load(f\"{subject}_acq-{shim_mode}_TB1map.nii.gz\")\n",
+    "        \n",
+    "        # Defining mask based on the magnitude image intensity threshold\n",
+    "        cslice=CP_anat.shape[2] // 2 -2 #shows the SC seg best\n",
+    "        threshold=300\n",
+    "        mask=CP_anat.get_fdata()[xmin:xmax,ymin:ymax, cslice]\n",
+    "        mask=np.where(mask > threshold, 1, 0)\n",
+    "        \n",
+    "        # Cropping anat, SC, B1+\n",
+    "        CP_anat=CP_anat.get_fdata()[xmin:xmax,ymin:ymax, cslice]\n",
+    "        CP_anat=CP_anat*mask\n",
+    "        CP_SC=CP_SC.get_fdata()[xmin:xmax,ymin:ymax, cslice]\n",
+    "        CP_nTpV=CP_nTpV.get_fdata()[xmin:xmax,ymin:ymax, cslice]\n",
+    "        CP_nTpV=CP_nTpV*mask\n",
+    "        \n",
+    "        # All opacity overalys look ugly: workaround, set the anat slice to a max value where the segmentation exists\n",
+    "        CP_anat[CP_SC>0.5]=4000;\n",
+    "        \n",
+    "        # Plotting anat overlayed with SC\n",
+    "        splot=axes[i]\n",
+    "        splot.imshow((CP_anat.T), cmap='gray', origin='lower',vmin=0,vmax=2000)#, interpolation='spline36')\n",
+    "        splot.set_title('SC overlay', size=font_size)\n",
+    "        splot.axis('off')\n",
+    "        \n",
+    "        # Plotting the B1+ map\n",
+    "        splot=axes[i+1]\n",
+    "        splot.imshow((CP_nTpV.T), cmap='viridis', origin='lower',vmin=dynmin,vmax=dynmax)#, interpolation='spline36')\n",
+    "        splot.set_title(shim_mode, size=font_size)\n",
+    "        splot.axis('off')\n",
+    "\n",
+    "        \n",
+    "    else:\n",
+    "        # Load data\n",
+    "        B1map=nib.load(f\"{subject}_acq-{shim_mode}_TB1map.nii.gz\")\n",
+    "        B1map=B1map.get_fdata()[xmin:xmax,ymin:ymax, cslice]\n",
+    "        B1map=B1map*mask\n",
+    "        \n",
+    "        # Plot\n",
+    "        splot=axes[i+1]\n",
+    "        im = splot.imshow((B1map.T), cmap='viridis', origin='lower',vmin=dynmin,vmax=dynmax)#, interpolation='spline36')\n",
+    "        splot.set_title(shim_mode, size=font_size)\n",
+    "        splot.axis('off')\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.subplots_adjust(wspace=0.1, hspace=0, right=0.9)\n",
+    "\n",
+    "# Colorbar\n",
+    "# Assume that the colorbar should start at the bottom of the lower row of subplots and\n",
+    "# extend to the top of the upper row of subplots\n",
+    "cbar_bottom = 0.06  # This might need adjustment\n",
+    "cbar_height = 0.83  # This represents the total height of both rows of subplots\n",
+    "cbar_ax = fig.add_axes([0.93, cbar_bottom, 0.03, cbar_height])\n",
+    "cbar = plt.colorbar(im, cax=cbar_ax)\n",
+    "\n",
+    "# cbar_ax = fig.add_axes([0.95, 0.5, 0.04, 0.4])\n",
+    "# cbar = plt.colorbar(im, cax=cbar_ax)\n",
+    "cbar_ax.set_title('nT/V', size=12)\n",
+    "plt.show\n"
    ]
   }
  ],