Merge pull request #89 from shimming-toolbox/jca/82-b1-value-in-plots

Display B1+ efficiency value inside the figure

Author: dpapp86

data_processing.ipynb

@@ -439,6 +439,7 @@
     "\n",
     "# Adjust the layout so labels and titles do not overlap\n",
     "plt.tight_layout()\n",
+    "plt.savefig(os.path.join(path_results, 'fig_gre_csf-sc_ratio.png'), dpi=300, format='png')\n",
     "plt.show()"
    ]
   },
@@ -734,6 +735,7 @@
     "\n",
     "# Adjust the layout so labels and titles do not overlap\n",
     "plt.tight_layout()\n",
+    "plt.savefig(os.path.join(path_results, 'fig_b1plus.png'), dpi=300, format='png')\n",
     "plt.show()"
    ]
   },
@@ -843,6 +845,14 @@
     "\n",
     "# Select subject to show\n",
     "subject = 'sub-05'\n",
+    "os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
+    "file_mask = f\"{subject}_acq-anatCP_TB1TFL_mask-shimming.nii.gz\"\n",
+    "\n",
+    "# Load the statistics from the CSV file\n",
+    "stats_df = pd.read_csv(os.path.join(path_results, 'stats_b1plus.csv'))\n",
+    "\n",
+    "# Filter for the specific subject\n",
+    "subject_stats = stats_df[stats_df['Subject'] == subject]\n",
     "\n",
     "# Defining crop limits for resulting figure\n",
     "xmin = 20\n",
@@ -859,55 +869,58 @@
     "font_size = 14\n",
     "axes=axes.flatten()\n",
     "\n",
-    "# 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",
+    "# First, plot the anatomical image with an overlay of the mask\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",
+    "# Load data\n",
+    "CP_anat=nib.load(f\"{subject}_acq-anatCP_TB1TFL.nii.gz\")\n",
+    "CP_SC=nib.load(file_mask)\n",
+    "CP_nTpV=nib.load(f\"{subject}_acq-CP_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[0]\n",
+    "splot.imshow((CP_anat.T), cmap='gray', origin='lower',vmin=0,vmax=2000)#, interpolation='spline36')\n",
+    "splot.set_title('Anat', size=font_size)\n",
+    "splot.axis('off')\n",
+    "\n",
+    "# Then, plot each B1+ map, with an overlay of the mean and CV inside the cord\n",
+    "for i,shim_mode in enumerate(shim_modes):\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",
+    "    # Find the statistics for the current shim mode\n",
+    "    shim_stats = subject_stats[subject_stats['Shim_Mode'] == shim_mode]\n",
+    "    if not shim_stats.empty:\n",
+    "        mean_val = shim_stats.iloc[0]['Average']\n",
+    "        std_val = shim_stats.iloc[0]['Standard_Deviation']\n",
+    "        cv = std_val / mean_val * 100  # Coefficient of variation in percentage\n",
+    "        annotation_text = f\"{mean_val:.2f} nT/V\\n{cv:.2f}%\"\n",
+    "        splot.annotate(annotation_text, (0.05, 0.95), xycoords='axes fraction', \n",
+    "                       fontsize=10, color='white', \n",
+    "                       verticalalignment='top', horizontalalignment='left')\n",
     "\n",
     "plt.tight_layout()\n",
     "plt.subplots_adjust(wspace=0.1, hspace=0, right=0.9)\n",
@@ -923,6 +936,7 @@
     "# 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.savefig(os.path.join(path_results, 'fig_b1plus_map.png'), dpi=300, format='png')\n",
     "plt.show\n"
    ]
   },