diff --git a/content/index.ipynb b/content/index.ipynb
index fc6dcf1..911d3d4 100644
--- a/content/index.ipynb
+++ b/content/index.ipynb
@@ -240,6 +240,406 @@
     "os.makedirs(path_results, exist_ok=True)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 3     |     Overview of processing pipeline\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Run segmentation on GRE scan\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide_input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "\n",
+    "if notebook != 'figures':\n",
+    "    for subject in subjects:\n",
+    "        os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "        # The \"CoV\" RF shimming scenario was chosen as the segmentation baseline due to the more homogeneous signal intensity in the I-->S direction, which results in a better segmentation peformance in the C7-T2 region\n",
+    "        fname_manual_seg = os.path.join(path_labels, subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-SC_seg.nii.gz\")\n",
+    "        if os.path.exists(fname_manual_seg):\n",
+    "            # Manual segmentation already exists. Copy it to local folder\n",
+    "            print(f\"{subject}: Manual segmentation found\\n\")\n",
+    "            shutil.copyfile(fname_manual_seg, f\"{subject}_acq-CoV_T2starw_seg.nii.gz\")\n",
+    "            # Generate QC report to make sure the manual segmentation is correct\n",
+    "            !sct_qc -i {subject}_acq-CoV_T2starw.nii.gz -s {subject}_acq-CoV_T2starw_seg.nii.gz -p sct_deepseg_sc -qc {path_qc} -qc-subject {subject}\n",
+    "        else:\n",
+    "            # Manual segmentation does not exist. Run automatic segmentation.\n",
+    "            print(f\"{subject}: Manual segmentation not found\")\n",
+    "            !sct_deepseg_sc -i {subject}_acq-CoV_T2starw.nii.gz -c t2s -qc {path_qc}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Copy CSF masks from the derivatives folder.\n",
+    "\n",
+    "For more details about how these masks were created, see: [https://github.com/shimming-toolbox/rf-shimming-7t/issues/67](https://github.com/shimming-toolbox/rf-shimming-7t/issues/67)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide_input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "if notebook != 'figures':\n",
+    "    for subject in subjects:\n",
+    "        os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "        fname_manual_seg = os.path.join(path_labels, subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz\")\n",
+    "        if os.path.exists(fname_manual_seg):\n",
+    "            # Manual segmentation already exists. Copy it to local folder\n",
+    "            print(f\"{subject}: Manual segmentation found\\n\")\n",
+    "            shutil.copyfile(fname_manual_seg, f\"{subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz\")\n",
+    "            # Generate QC report to make sure the manual segmentation is correct\n",
+    "            !sct_qc -i {subject}_acq-CoV_T2starw.nii.gz -s {subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz -p sct_deepseg_sc -qc {path_qc} -qc-subject {subject}\n",
+    "        else:\n",
+    "            # Manual segmentation does not exist. Panic!\n",
+    "            print(f\"{subject}: Manual segmentation not found 😱\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Crop GRE scan for faster processing and better registration results\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide_input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "if notebook != 'figures':\n",
+    "    for subject in subjects:\n",
+    "        os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "        !sct_crop_image -i {subject}_acq-CoV_T2starw.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop.nii.gz\n",
+    "        !sct_crop_image -i {subject}_acq-CoV_T2starw_seg.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop_seg.nii.gz\n",
+    "        !sct_crop_image -i {subject}_acq-CoV_T2starw_label-CSF_seg.nii.gz -m {subject}_acq-CoV_T2starw_seg.nii.gz -dilate 20x20x0 -o {subject}_acq-CoV_T2starw_crop_label-CSF_seg.nii.gz"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Label vertebrae on GRE scan"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide_input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "\n",
+    "# Given the low resolution of the GRE scan, the automatic detection of C2-C3 disc is unreliable. Therefore we need to use the manual disc labels that are part of the dataset.\n",
+    "if notebook != 'figures':\n",
+    "    for subject in subjects:\n",
+    "        os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "        fname_label_discs = os.path.join(path_data, \"derivatives\", \"labels\", subject, \"anat\", f\"{subject}_acq-CoV_T2starw_label-discs_dseg.nii.gz\")\n",
+    "        !sct_label_utils -i {subject}_acq-CoV_T2starw_crop_seg.nii.gz -disc {fname_label_discs} -o {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz\n",
+    "        # Generate QC report to assess labeled segmentation\n",
+    "        !sct_qc -i {subject}_acq-CoV_T2starw_crop.nii.gz -s {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -p sct_label_vertebrae -qc {path_qc} -qc-subject {subject}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Register *_T2starw to CoV_T2starw ⏳"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide_input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "if notebook != 'figures':\n",
+    "    for subject in subjects:\n",
+    "        os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "        for shim_mode in shim_modes:\n",
+    "            # Don't do it for CoV_T2starw\n",
+    "            if shim_mode != 'CoV':\n",
+    "                !sct_register_multimodal -i {subject}_acq-{shim_mode}_T2starw.nii.gz -d {subject}_acq-CoV_T2starw_crop.nii.gz -dseg {subject}_acq-CoV_T2starw_crop_seg.nii.gz -param step=1,type=im,algo=dl -qc {path_qc}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide_input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "\n",
+    "if notebook != 'figures':\n",
+    "    for subject in subjects:\n",
+    "        os.chdir(os.path.join(path_data, subject, \"anat\"))\n",
+    "        for shim_mode in shim_modes:\n",
+    "            # Shim methods are registered to the CoV T2starw scan, so we need to use the added suffix to identify them\n",
+    "            if shim_mode == 'CoV':\n",
+    "                file_suffix = 'crop'\n",
+    "            else:\n",
+    "                file_suffix = 'reg'\n",
+    "            fname_result_sc = os.path.join(path_results, f\"{subject}_acq-{shim_mode}_T2starw_label-SC.csv\")\n",
+    "            !sct_extract_metric -i {subject}_acq-{shim_mode}_T2starw_{file_suffix}.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 -perslice 1 -o {fname_result_sc}\n",
+    "            fname_result_csf = os.path.join(path_results, f\"{subject}_acq-{shim_mode}_T2starw_label-CSF.csv\")\n",
+    "            !sct_extract_metric -i {subject}_acq-{shim_mode}_T2starw_{file_suffix}.nii.gz -f {subject}_acq-CoV_T2starw_crop_label-CSF_seg.nii.gz -method wa -vert 3:9 -vertfile {subject}_acq-CoV_T2starw_crop_seg_labeled.nii.gz -perslice 1 -o {fname_result_csf}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Make figure of CSF/SC signal ratio from T2starw scan"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide_input",
+     "report_output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# Make figure of CSF/SC signal ratio from T2starw scan\n",
+    "\n",
+    "# Go back to root data folder\n",
+    "os.chdir(path_data)\n",
+    "\n",
+    "def smooth_data(data, window_size=50):\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 C3 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])\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",
+    "# Vertebral level labels\n",
+    "vertebral_levels = [\"C3\", \"C4\", \"C5\", \"C6\", \"C7\", \"T1\", \"T2\"]\n",
+    "label_positions = normalized_vector[:-1] + np.diff(normalized_vector) / 2\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",
+    "line_width = 3\n",
+    "\n",
+    "# Check if axes is an array or a single object\n",
+    "if n_rows == 1:\n",
+    "    axes = [axes]\n",
+    "\n",
+    "# Data storage for statistics\n",
+    "data_stats = []\n",
+    "\n",
+    "# Iterate over each subject and create a subplot\n",
+    "for i, subject in enumerate(subjects):\n",
+    "    ax = axes[i]\n",
+    "\n",
+    "    for shim_mode in shim_modes:\n",
+    "        # Initialize list to collect data for this shim method\n",
+    "        method_data = []\n",
+    "\n",
+    "        # Get signal in SC\n",
+    "        file_csv = os.path.join(path_results, f\"{subject}_acq-{shim_mode}_T2starw_label-SC.csv\")\n",
+    "        df = pd.read_csv(file_csv)\n",
+    "        data_sc = df['WA()']\n",
+    "        data_sc_smoothed = smooth_data(data_sc)\n",
+    "\n",
+    "        # Get signal in CSF\n",
+    "        file_csv = os.path.join(path_results, f\"{subject}_acq-{shim_mode}_T2starw_label-CSF.csv\")\n",
+    "        df = pd.read_csv(file_csv)\n",
+    "        data_csf = df['WA()']\n",
+    "        data_csf_smoothed = smooth_data(data_csf)\n",
+    "        \n",
+    "        # Compute ratio\n",
+    "        data_sc_csf_ratio = data_csf_smoothed / data_sc_smoothed\n",
+    "\n",
+    "        # Normalize the x-axis to a 1-0 scale for each subject (to go from superior-inferior direction)\n",
+    "        x_subject = np.linspace(1, 0, len(data_sc_csf_ratio))\n",
+    "\n",
+    "        # Interpolate to the fixed grid\n",
+    "        interp_func = interp1d(x_subject, data_sc_csf_ratio, kind='linear', bounds_error=False, fill_value='extrapolate')\n",
+    "        resampled_data = interp_func(x_grid)\n",
+    "\n",
+    "        method_data.append(resampled_data)\n",
+    "\n",
+    "        # If there's data for this shim method, plot it and compute stats\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_mode}\", linewidth=line_width)\n",
+    "        \n",
+    "            # Compute stats on the non-resampled data (to avoid interpolation errors)\n",
+    "            mean_data = np.mean(data_sc_csf_ratio)\n",
+    "            sd_data = np.std(data_sc_csf_ratio)\n",
+    "            data_stats.append([subject, shim_mode, mean_data, sd_data])\n",
+    "\n",
+    "    # Set x-ticks and labels for the bottom subplot\n",
+    "    if i == n_rows - 1:  # Check if it's the last subplot\n",
+    "        # Create a secondary axis for vertebral level labels\n",
+    "        secax = ax.secondary_xaxis('bottom')\n",
+    "        secax.set_xticks(label_positions)\n",
+    "        secax.set_xticklabels(vertebral_levels, fontsize=font_size-4)\n",
+    "        secax.tick_params(axis='x', which='major', length=0)  # Hide tick marks\n",
+    "\n",
+    "    ax.set_xticks(custom_xticks)\n",
+    "    ax.set_xticklabels([''] * len(custom_xticks))\n",
+    "    ax.tick_params(axis='x', which='major', length=0)  # Hide tick marks for the primary axis\n",
+    "\n",
+    "    ax.set_title(f'{subject}', fontsize=font_size)\n",
+    "    ax.set_ylabel('CSF/Cord T2starw signal ratio', 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.savefig(os.path.join(path_results, 'fig_gre_csf-sc_ratio.png'), dpi=300, format='png')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Perform statistics"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide_input",
+     "report_output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# Perform statistics\n",
+    "\n",
+    "# Convert to DataFrame and save to CSV\n",
+    "df_stats = pd.DataFrame(data_stats, columns=['Subject', 'Shim_Mode', 'Average', 'Standard_Deviation'])\n",
+    "df_stats.to_csv(os.path.join(path_results, 'stats_t2starw.csv'), index=False)\n",
+    "\n",
+    "# Compute statistics across subjects\n",
+    "grouped_means = df_stats.groupby('Shim_Mode').agg({'Average': ['mean', 'std'], 'Standard_Deviation': ['mean', 'std']})\n",
+    "\n",
+    "# Format mean ± standard deviation\n",
+    "grouped_means['Average_formatted'] = grouped_means['Average']['mean'].map(\"{:.2f}\".format) + \" ± \" + grouped_means['Average']['std'].map(\"{:.2f}\".format)\n",
+    "grouped_means['Standard_Deviation_formatted'] = grouped_means['Standard_Deviation']['mean'].map(\"{:.2f}\".format) + \" ± \" + grouped_means['Standard_Deviation']['std'].map(\"{:.2f}\".format)\n",
+    "\n",
+    "# Drop multi-level index and only keep formatted columns\n",
+    "grouped_means = grouped_means.drop(columns=['Average', 'Standard_Deviation'])\n",
+    "grouped_means.columns = ['Average', 'Standard_Deviation']  # Rename columns for clarity\n",
+    "grouped_means.reset_index().to_csv(os.path.join(path_results, 'stats_t2starw_average_across_subjects.csv'), index=False)\n",
+    "\n",
+    "# Perform Repeated Measures ANOVA\n",
+    "aovrm = AnovaRM(df_stats, 'Average', 'Subject', within=['Shim_Mode'])\n",
+    "res = aovrm.fit()\n",
+    "print(res.summary())\n",
+    "\n",
+    "# Perform Post Hoc paired t-tests\n",
+    "\n",
+    "# Filter out subjects that don't have all shim modes\n",
+    "shim_modes = df_stats['Shim_Mode'].unique()\n",
+    "valid_subjects = df_stats.groupby('Subject').filter(lambda x: all(mode in x['Shim_Mode'].values for mode in shim_modes))\n",
+    "\n",
+    "# Pairwise comparisons\n",
+    "pairs = list(itertools.combinations(shim_modes, 2))\n",
+    "p_values = []\n",
+    "\n",
+    "for pair in pairs:\n",
+    "    data1 = valid_subjects[valid_subjects['Shim_Mode'] == pair[0]]['Average']\n",
+    "    data2 = valid_subjects[valid_subjects['Shim_Mode'] == pair[1]]['Average']\n",
+    "    _, p_value = ttest_rel(data1, data2)\n",
+    "    p_values.append(p_value)\n",
+    "\n",
+    "# Function for Benjamini-Hochberg FDR correction\n",
+    "def benjamini_hochberg(p_values):\n",
+    "    n = len(p_values)\n",
+    "    sorted_p_values = np.sort(p_values)\n",
+    "    sorted_index = np.argsort(p_values)\n",
+    "    adjusted_p_values = np.zeros(n)\n",
+    "\n",
+    "    for i, p in enumerate(sorted_p_values):\n",
+    "        adjusted_p_values[sorted_index[i]] = min(p * n / (i + 1), 1)\n",
+    "\n",
+    "    return adjusted_p_values\n",
+    "\n",
+    "# Applying Benjamini-Hochberg FDR correction\n",
+    "adjusted_p_values = benjamini_hochberg(p_values)\n",
+    "\n",
+    "# Creating a summary DataFrame\n",
+    "comparison_results = pd.DataFrame({'Pair': pairs, 'P-Value': p_values, 'Adjusted P-Value': adjusted_p_values})\n",
+    "\n",
+    "# Save to CSV\n",
+    "comparison_results.to_csv(os.path.join(path_results, 'stats_t2starw_paired_ttest.csv'), index=False)\n",
+    "\n",
+    "print(f\"Paired t-tests:\\n{comparison_results}\")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},