Skip to content

Commit

Permalink
Add MPRAGE processing and visualisation (#46)
Browse files Browse the repository at this point in the history 
Fixes #37
Fixes #38
Fixes #39
Fixes #43 
Fixes #45
  • Loading branch information
@dpapp86
dpapp86 authored Jan 16, 2024
1 parent 9cddd3e commit 30753b9
Showing 1 changed file with 329 additions and 5 deletions.
334 changes: 329 additions & 5 deletions data_processing.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,11 @@
"- Apply the computed warping field to bring the segmentation and vertebral levels to the B1 map\n",
"- Convert the B1 map to nT/V units\n",
"- Extract the B1 map value within the spinal cord\n",
"- Segment the spinal cord on the CoV MPRAGE scan\n",
"- Label the vertebral levels using automatic labeling\n",
"- Register the CP mode MPRAGE to the CoV scan\n",
"- Warp segmentation and labeling to the CP mode scan\n",
"- Visualize the CP mode and CoV mode MPRAGE scan, calculate coefficient of variation within the cord for both, plot signal intensity\n",
"\n",
">Slow processes are indicated with the emoji ⏳"
]
Expand Down Expand Up @@ -166,7 +171,8 @@
"path_data = os.getcwd()\n",
"print(f\"path_data: {path_data}\")\n",
"path_qc = os.path.join(path_data, \"qc\")\n",
"shim_modes = [\"CP\", \"CoV\", \"patient\", \"phase\", \"SAReff\", \"target\", \"volume\"]\n",
"shim_modes = [\"CP\", \"patient\", \"volume\", \"phase\", \"CoV\", \"target\", \"SAReff\"]\n",
"shim_modes_MPRAGE = [\"CP\", \"CoV\"]\n",
"# shim_modes = [\"CP\", \"CoV\"] # debugging\n",
"print(f\"shim_modes: {shim_modes}\")\n",
"subjects = sorted(glob.glob(\"sub-*\"))\n",
Expand All @@ -185,6 +191,9 @@
"outputs": [],
"source": [
"# Run segmentation on GRE scan\n",
"# ℹ️ The \"CoV reduction\" RF shimming scenario was chosen as the segmentation baseline due to the more \n",
"# homogeneous signal intensity in the I-->S direction, which results in a better segmentation peformance\n",
"# in the C7-T2 region\n",
"\n",
"for subject in subjects:\n",
" os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
Expand Down Expand Up @@ -259,12 +268,12 @@
"metadata": {},
"outputs": [],
"source": [
"# Extract the signal intensity on the GRE scan within the spinal cord between levels C1 and T2 (included)\n",
"# Extract the signal intensity on the GRE scan within the spinal cord between levels C3 and T2 (included), which correspond to the region where RF shimming was prescribed\n",
"\n",
"for subject in subjects:\n",
" # TODO: loop across other shim methods \n",
" os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
" !sct_extract_metric -i {subject}_acq-CoV_T2starw_crop.nii.gz -f {subject}_acq-CoV_T2starw_crop_seg.nii.gz -method wa -vert 1:9 -vertfile {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -append 1 -perslice 1 -o gre_CoV.csv\n"
" !sct_extract_metric -i {subject}_acq-CoV_T2starw_crop.nii.gz -f {subject}_acq-CoV_T2starw_crop_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -append 1 -perslice 1 -o gre_CoV.csv\n"
]
},
{
Expand Down Expand Up @@ -359,12 +368,12 @@
"metadata": {},
"outputs": [],
"source": [
"# Extract B1+ value along the spinal cord between levels C1 and T2 (included)\n",
"# Extract B1+ value along the spinal cord between levels C3 and T2 (included)\n",
"\n",
"for subject in subjects:\n",
" os.chdir(os.path.join(path_data, subject, \"fmap\"))\n",
" for shim_mode in shim_modes:\n",
" !sct_extract_metric -i {subject}_acq-{shim_mode}_TB1map.nii.gz -f {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz -method wa -vert 1:9 -vertfile {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz -perslice 1 -o TB1map_{shim_mode}.csv"
" !sct_extract_metric -i {subject}_acq-{shim_mode}_TB1map.nii.gz -f {subject}_acq-anat{shim_mode}_TB1TFL_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-anat{shim_mode}_TB1TFL_seg_labeled.nii.gz -perslice 1 -o TB1map_{shim_mode}.csv"
]
},
{
Expand Down Expand Up @@ -480,6 +489,8 @@
" all_data = []\n",
" for file_path in glob.glob(os.path.join(path_data, subject, \"fmap\", f\"TB1map_*{shim_mode}*.csv\")):\n",
" df = pd.read_csv(file_path)\n",
" # Drop rows where VertLevel is 1 or 2 (because RF shimming was not performed at these levels)\n",
" df = df[df['VertLevel'].isin(['1', '2']) == False]\n",
" wa_data = df['WA()']\n",
" all_data.extend(wa_data)\n",
"\n",
Expand Down Expand Up @@ -512,6 +523,319 @@
" posthoc_result = pairwise_tukeyhsd(df_summary['Average'], df_summary['Shim_Mode'])\n",
" print(\"Posthoc Tukey HSD Result:\\n\", posthoc_result)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a916a6f7",
"metadata": {},
"outputs": [],
"source": [
"# Visualize RF maps to generate Figure 2\n",
"\n",
"# TODO"
]
},
{
"cell_type": "markdown",
"id": "cc07e3e6",
"metadata": {},
"source": [
"# MPRAGE PROCESSING BLOCK\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "49a99b31",
"metadata": {},
"outputs": [],
"source": [
"# Run segmentation on MPRAGE CoV scan\n",
"\n",
"for subject in subjects:\n",
" os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
" !sct_deepseg_sc -i {subject}_acq-CoV_T1w.nii.gz -c t1 -qc {path_qc}"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "663bb98c",
"metadata": {},
"outputs": [],
"source": [
"#Crop the scans for easier and faster coregistration\n",
"#Using a more permissive bounding box than GRE\n",
"\n",
"for subject in subjects:\n",
" os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
" !sct_crop_image -i {subject}_acq-CoV_T1w.nii.gz -m {subject}_acq-CoV_T1w_seg.nii.gz -dilate 34x34x0 -o {subject}_acq-CoV_T1w_crop.nii.gz\n",
" !sct_crop_image -i {subject}_acq-CoV_T1w_seg.nii.gz -m {subject}_acq-CoV_T1w_seg.nii.gz -dilate 34x34x0 -o {subject}_acq-CoV_T1w_crop_seg.nii.gz"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "75071e07",
"metadata": {},
"outputs": [],
"source": [
"# Vertebral labeling on the CoV MPRAGE scan\n",
"\n",
"for subject in subjects:\n",
" os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
" !sct_label_vertebrae -i {subject}_acq-CoV_T1w_crop.nii.gz -s {subject}_acq-CoV_T1w_crop_seg.nii.gz -c t1 -qc {path_qc}"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "436a2728",
"metadata": {},
"outputs": [],
"source": [
"#Coregistering the CP mode to the CoV case\n",
"#Seems to be working good enough\n",
"for subject in subjects:\n",
" os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
" !sct_register_multimodal -i {subject}_acq-CP_T1w.nii.gz -d {subject}_acq-CoV_T1w_crop.nii.gz -dseg {subject}_acq-CoV_T1w_crop_seg.nii.gz -qc {path_qc}"
]
},
{
"cell_type": "markdown",
"id": "a2074892",
"metadata": {},
"source": [
"## Verify QC report (MPRAGE segmentation and coregistration)\n",
"\n",
"\n",
"Open the QC report located under `ds004906/qc/index.html`. Make sure the registration are correct before resuming the analysis."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0ffcfb28",
"metadata": {},
"outputs": [],
"source": [
"# UNNECESSARY, but keeping code it in case we want to register the GRE scans to the T1w scans.\n",
"# https://github.com/shimming-toolbox/rf-shimming-7t/issues/49\n",
"\n",
"# Warping spinal cord segmentation and vertebral level to the CP mode scan for each subject\n",
"\n",
"# for subject in subjects:\n",
"# os.chdir(os.path.join(path_data, subject, \"anat\")) \n",
"# !sct_apply_transfo -i {subject}_acq-CoV_T1w_crop_seg.nii.gz -d {subject}_acq-CP_T1w.nii.gz -w warp_{subject}_acq-CoV_T1w_crop2{subject}_acq-CP_T1w.nii.gz -x linear -o {subject}_acq-CP_T1w_seg_reg.nii.gz\n",
"# !sct_apply_transfo -i {subject}_acq-CoV_T1w_crop_seg_labeled.nii.gz -d {subject}_acq-CP_T1w.nii.gz -w warp_{subject}_acq-CoV_T1w_crop2{subject}_acq-CP_T1w.nii.gz -x linear -o {subject}_acq-CP_T1w_seg_labeled_reg.nii.gz"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "75c43da0",
"metadata": {},
"outputs": [],
"source": [
"# Extract the signal intensity on the MPRAGE scan within the spinal cord between levels C3 and T2 (included)\n",
"\n",
"for subject in subjects: \n",
" os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
" #CP\n",
" !sct_extract_metric -i {subject}_acq-CP_T1w_reg.nii.gz -f {subject}_acq-CoV_T1w_crop_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T1w_crop_seg_labeled.nii.gz -append 1 -perslice 1 -o {subject}_acq-CP_T1w.csv\n",
" #CoV\n",
" !sct_extract_metric -i {subject}_acq-CoV_T1w_crop.nii.gz -f {subject}_acq-CoV_T1w_crop_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T1w_crop_seg_labeled.nii.gz -append 1 -perslice 1 -o {subject}_acq-CoV_T1w.csv"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "63e6181a",
"metadata": {},
"outputs": [],
"source": [
"#Flattening sagittaly\n",
"for subject in subjects:\n",
" os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
" #flatten \n",
" !sct_flatten_sagittal -i {subject}_acq-CoV_T1w_crop.nii.gz -s {subject}_acq-CoV_T1w_crop_seg.nii.gz\n",
" !sct_flatten_sagittal -i {subject}_acq-CP_T1w_reg.nii.gz -s {subject}_acq-CoV_T1w_crop_seg.nii.gz"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "77d24f64",
"metadata": {},
"outputs": [],
"source": [
"# Display images\n",
"\n",
"for subject in subjects:\n",
" os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
" CoV_flatten = os.path.join(path_data, subject, \"anat\", f\"{subject}_acq-CoV_T1w_crop_flatten.nii.gz\")\n",
" CoV_flatten_data = nib.load(CoV_flatten).get_fdata()\n",
" CP_flatten = os.path.join(path_data, subject, \"anat\", f\"{subject}_acq-CP_T1w_reg_flatten.nii.gz\")\n",
" CP_flatten_data=nib.load(CP_flatten).get_fdata()\n",
"\n",
" # Data manipulation to get the central slice\n",
" # Already cropped inf-sup and right-left, no need to crop more\n",
" CoV_flatten_data_cslice=CoV_flatten_data[CoV_flatten_data.shape[0] // 2,:,:]\n",
" CP_flatten_data_cslice=CP_flatten_data[CP_flatten_data.shape[0] // 2,:,:]\n",
" #And now plot side by side\n",
" plt.figure(figsize=(3, 4))\n",
" vmin = -1 \n",
" vmax = -0.70\n",
"\n",
" plt.subplot(1, 2, 1)\n",
" plt.imshow(CP_flatten_data_cslice.T, cmap='gray', origin='lower',vmin=vmin,vmax=vmax) \n",
" plt.title('CP mode',color='blue')\n",
" plt.axis('off')\n",
"\n",
" plt.subplot(1, 2, 2)\n",
" plt.imshow(CoV_flatten_data_cslice.T, cmap='gray', origin='lower',vmin=vmin, vmax=vmax)\n",
" plt.title('CoV reduction',color='orange')\n",
" plt.axis('off')\n",
"\n",
" # Add the subject title at the top\n",
" plt.suptitle(subject, y=0.9) # Adjust 'y' as needed for position\n",
" plt.tight_layout()\n",
" plt.subplots_adjust(wspace=0.1, hspace=0)\n",
" plt.show() "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cd2a3c9c",
"metadata": {},
"outputs": [],
"source": [
"# Generate signal intensity curves\n",
"# Make figure of B1+ values along the spinal cord across shim methods\n",
"\n",
"# Go back to root data folder\n",
"os.chdir(os.path.join(path_data))\n",
"\n",
"def smooth_data(data, window_size=11):\n",
" \"\"\" Apply a simple moving average to smooth the data. \"\"\"\n",
" return uniform_filter1d(data, size=window_size, mode='nearest')\n",
"\n",
"# Fixed grid for x-axis\n",
"x_grid = np.linspace(0, 1, 100)\n",
"\n",
"# z-slices corresponding to levels C1 to T2 on the PAM50 template. These will be used to scale the x-label of each subject.\n",
"original_vector =np.array([907, 870, 833, 800, 769, 735, 692, 646]) # np.array([984, 938, 907, 870, 833, 800, 769, 735, 692, 646])\n",
"\n",
"# Normalize the PAM50 z-slice numbers to the 1-0 range (to show inferior-superior instead of superior-inferior)\n",
"min_val = original_vector.min()\n",
"max_val = original_vector.max()\n",
"normalized_vector = 1 - ((original_vector - min_val) / (max_val - min_val))\n",
"\n",
"# Use this normalized vector as x-ticks\n",
"custom_xticks = normalized_vector\n",
"\n",
"# Number of subjects determines the number of rows in the subplot\n",
"n_rows = len(subjects)\n",
"\n",
"# Create a figure with multiple subplots\n",
"fig, axes = plt.subplots(n_rows, 1, figsize=(10, 6 * n_rows))\n",
"font_size = 18\n",
"\n",
"# Check if axes is an array or a single object\n",
"if n_rows == 1:\n",
" axes = [axes]\n",
"\n",
"# Iterate over each subject and create a subplot\n",
"for i, subject in enumerate(subjects):\n",
" ax = axes[i]\n",
" \n",
" os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
"\n",
" for shim_method in shim_modes_MPRAGE:\n",
" # Initialize list to collect data for this shim method\n",
" method_data = []\n",
"\n",
" for file_path in glob.glob(f\"{subject}_acq-{shim_method}_T1w.csv\"):\n",
" df = pd.read_csv(file_path)\n",
" wa_data = df['WA()']\n",
"\n",
" # Normalize the x-axis to a 0-1 scale for each subject\n",
" x_subject = np.linspace(0, 1, len(wa_data))\n",
"\n",
" # Interpolate to the fixed grid\n",
" interp_func = interp1d(x_subject, wa_data, kind='linear', bounds_error=False, fill_value='extrapolate')\n",
" resampled_data = interp_func(x_grid)\n",
"\n",
" # Apply smoothing\n",
" smoothed_data = smooth_data(resampled_data)\n",
"\n",
" method_data.append(smoothed_data)\n",
"\n",
" # If there's data for this shim method, plot it\n",
" if method_data:\n",
" # Plotting each file's data separately\n",
" for resampled_data in method_data:\n",
" ax.plot(x_grid, resampled_data, label=f\"{shim_method}\")\n",
"\n",
" # Set custom x-ticks\n",
" ax.set_xticks(custom_xticks)\n",
" ax.set_xticklabels([''] * len(original_vector))\n",
"\n",
" ax.set_title(f'{subject}', fontsize=font_size)\n",
" ax.set_ylabel('Signal intensity [a.u.]', fontsize=font_size)\n",
" ax.tick_params(axis='y', which='major', labelsize=font_size-4)\n",
"\n",
" # Add legend only to the first subplot\n",
" if i == 0:\n",
" ax.legend(fontsize=font_size)\n",
"\n",
" ax.grid(True)\n",
"\n",
"# Adjust the layout so labels and titles do not overlap\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1a690f03",
"metadata": {},
"outputs": [],
"source": [
"# Create mean and CoV \n",
"\n",
"# Go back to root data folder\n",
"os.chdir(path_data)\n",
"\n",
"# Data storage\n",
"data_summary = []\n",
"\n",
"# Compute mean and SD of MPRAGE sigint for each subject and each shim mode\n",
"for subject in subjects:\n",
" for shim_mode in shim_modes_MPRAGE:\n",
" all_data = []\n",
" for file_path in glob.glob(os.path.join(path_data, subject, \"anat\", f\"{subject}_acq-{shim_mode}_T1w.csv\")):\n",
" # print(file_path)\n",
" df = pd.read_csv(file_path)\n",
" wa_data = df['WA()']\n",
" all_data.extend(wa_data)\n",
"\n",
" if all_data:\n",
" mean_data = np.mean(all_data)\n",
" sd_data = np.std(all_data)\n",
" CoV_data=np.round((sd_data/mean_data)*100,2)\n",
" data_summary.append([subject, shim_mode, mean_data, sd_data, CoV_data])\n",
"\n",
"# Convert to DataFrame and save to CSV\n",
"df_summary = pd.DataFrame(data_summary, columns=['Subject', 'Shim_Mode', 'Average', 'Standard_Deviation', 'CoV'])\n",
"df_summary.to_csv(os.path.join(path_results, 'subject_shim_mode_summary.csv'), index=False)\n",
"\n",
"df_summary_onlyCoV=df_summary[['Subject', 'Shim_Mode', 'CoV']]\n",
"print(df_summary_onlyCoV)\n"
]
}
],
"metadata": {
Expand Down

0 comments on commit 30753b9

Please sign in to comment.