Work on plotting.

Author: tsalo

analysis/plot_aslprep_cbf.py

@@ -1,4 +1,4 @@
-"""Plot scalar maps from QSIRecon."""
+"""Plot CBF maps from ASLPrep."""
 
 import os
 from glob import glob
@@ -48,4 +48,3 @@
         fig.suptitle(title)
         plt.savefig(os.path.join(out_dir, f"{title.replace(' ', '_')}.png"))
         plt.close()
-

analysis/plot_rawdata.py

@@ -0,0 +1,240 @@
+#!/usr/bin/env python
+
+import shutil
+import subprocess
+from pathlib import Path
+
+import ants
+import imageio.v3 as iio
+import matplotlib.pyplot as plt
+import nibabel as nb
+import numpy as np
+from nilearn import image, plotting
+
+
+data_root = Path("/Users/mcieslak/projects/hbcd/pipeline_paper/data_examples")
+
+subid = "9645436710"
+sesid = "V02"
+
+# The dir-PA nifti is always the first
+raw_nifti = data_root / "assembly_bids" / f"sub-{subid}_ses-{sesid}_dir-PA_run-1_dwi.nii"
+processed_nifti = (
+    data_root
+    / "derivatives"
+    / "qsiprep"
+    / f"sub-{subid}_ses-{sesid}_space-ACPC_desc-preproc_dwi.nii"
+)
+raw_to_acpc_xfm = processed_nifti.parent / f"sub-{subid}_ses-{sesid}_from-raw_to-ACPC_rigid.mat"
+raw_mean_path = raw_nifti.parent / f"sub-{subid}_ses-{sesid}_dir-PA_run-1_dwi_mean.nii"
+processed_mean_path = (
+    processed_nifti.parent / f"sub-{subid}_ses-{sesid}_space-ACPC_desc-preproc_dwi_mean.nii"
+)
+
+# If there is not transform file, we need to run the registration
+if not Path(raw_to_acpc_xfm).exists():
+    # Load the nifti files
+    raw_img = nb.load(raw_nifti)
+    processed_img = nb.load(processed_nifti)
+
+    # Get the data and calculate mean of first 4 volumes
+    raw_data = raw_img.get_fdata()
+    processed_data = processed_img.get_fdata()
+
+    raw_mean = raw_data[..., :4].mean(axis=3)
+    processed_mean = processed_data[..., :4].mean(axis=3)
+
+    # Create new nifti images with the mean data
+    raw_mean_img = nb.Nifti1Image(raw_mean, raw_img.affine, raw_img.header)
+    processed_mean_img = nb.Nifti1Image(processed_mean, processed_img.affine, processed_img.header)
+
+    # Save the mean images
+    nb.save(raw_mean_img, raw_mean_path)
+    nb.save(processed_mean_img, processed_mean_path)
+
+    # Load the mean images directly with ANTs
+    raw_mean_ants = ants.image_read(str(raw_mean_path))
+    processed_mean_ants = ants.image_read(str(processed_mean_path))
+
+    # Perform rigid registration
+    registration = ants.registration(
+        fixed=processed_mean_ants, moving=raw_mean_ants, type_of_transform="Similarity"
+    )
+    shutil.copy(registration["fwdtransforms"][0], raw_to_acpc_xfm)
+    # Apply transform using lanczoswindowedsinc interpolation
+    registered_processed = ants.apply_transforms(
+        fixed=processed_mean_ants,
+        moving=raw_mean_ants,
+        transformlist=registration["fwdtransforms"],
+        interpolator="lanczosWindowedSinc",
+    )
+
+    # Save the registered image
+    registered_path = (
+        processed_nifti.parent / f"sub-{subid}_ses-{sesid}_space-ACPC_desc-preproc_dwi_mean.nii"
+    )
+    ants.image_write(registered_processed, str(registered_path))
+
+# use nilearn image.crop_img to get the cropped reference image
+ref_img_path = raw_nifti.parent / f"sub-{subid}_ses-{sesid}_space-ACPC_desc-cropped_dwi.nii"
+if not Path(ref_img_path).exists():
+    # Call 3dAutobox using subprocess
+    cmd = [
+        "3dAutobox",
+        "-prefix",
+        str(ref_img_path),
+        "-input",
+        str(processed_mean_path),
+        "-npad",
+        "5",
+    ]
+
+    process = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+
+    stdout, stderr = process.communicate()
+
+    if process.returncode != 0:
+        raise RuntimeError(f"3dAutobox failed with error: {stderr.decode()}")
+
+
+def resample_processed_into_raw(image_index):
+    """Select a 3d volume from raw_nifti and transform the corresponding
+    volume from processed_nifti so they can be plotted next to each other.
+
+    Parameters
+    ----------
+    image_index : int
+        The volume number to extract from the 4D datasets
+    """
+    # Load the reference image
+    ref_img = ants.image_read(str(ref_img_path))
+
+    # Load the image from image_index using nilearn, save as a single-volume nifti
+    # No need to transform the processed image, it is already in ACPC space
+    processed_vol_path = (
+        processed_nifti.parent
+        / f"sub-{subid}_ses-{sesid}_space-ACPC_desc-preproc_dwi_vol-{image_index}.nii"
+    )
+    processed_vol = image.index_img(str(processed_nifti), image_index)
+    processed_vol.to_filename(processed_vol_path)
+
+    # Transform the raw image at image_index to the ACPC space
+    raw_nii = image.index_img(str(raw_nifti), image_index)
+    raw_nii_path = raw_nifti.parent / f"sub-{subid}_ses-{sesid}_vol-{image_index}_raw.nii"
+    raw_nii.to_filename(raw_nii_path)
+    raw_ants = ants.image_read(str(raw_nii_path))
+    raw_vol = ants.apply_transforms(
+        fixed=ref_img,
+        moving=raw_ants,
+        transformlist=[str(raw_to_acpc_xfm)],
+        interpolator="lanczosWindowedSinc",
+    )
+    resampled_raw_nii_path = (
+        raw_nifti.parent / f"sub-{subid}_ses-{sesid}_vol-{image_index}_space-ACPC_raw.nii"
+    )
+    ants.image_write(raw_vol, str(resampled_raw_nii_path))
+
+    temp_files_to_delete = [
+        raw_nii_path,
+    ]
+
+    for file in temp_files_to_delete:
+        Path(file).unlink()
+
+    return resampled_raw_nii_path, processed_vol_path
+
+
+def make_figure(raw_nii_path, registered_nii_path, image_index, crop_proportion=0.1):
+    """Create a figure comparing raw and registered volumes.
+
+    Parameters
+    ----------
+    raw_nii_path : str or Path
+        Path to the raw nifti file
+    registered_nii_path : str or Path
+        Path to the registered nifti file
+    image_index : int
+        Volume number for labeling
+    crop_proportion : float, optional
+        Proportion of image edges to crop (default: 0.15 for 15%)
+    """
+    # Load the specific volumes
+    raw_nii = image.load_img(raw_nii_path)
+    registered_nii = image.load_img(registered_nii_path)
+
+    slices_to_plot = {
+        "x": [-10.0],
+        "y": [-25.0],
+        "z": [4.0],
+    }
+
+    # Calculate number of rows needed for the plots
+    n_planes = len(slices_to_plot)
+
+    # Create the figure
+    fig, axes = plt.subplots(n_planes, 2, figsize=(15, 6.5 * n_planes))
+    # Get the data range for this plane
+    raw_data = raw_nii.get_fdata()
+    reg_data = registered_nii.get_fdata()
+
+    # Calculate 2nd and 98th percentiles from combined data, clipped at 0
+    all_data = np.clip(np.concatenate([raw_data.ravel(), reg_data.ravel()]), 0, None)
+    vmin = np.percentile(all_data, 0.5)
+    vmax = np.percentile(all_data, 99.5)
+
+    # Plot each plane
+    for idx, (plane, coords) in enumerate(slices_to_plot.items()):
+        # Plot raw volume
+        plotting.plot_anat(
+            raw_nii,
+            display_mode=plane,
+            cut_coords=coords,
+            title="",
+            axes=axes[idx, 0],
+            draw_cross=False,
+            annotate=False,
+            vmin=vmin,
+            vmax=vmax,
+        )
+
+        # Plot registered volume
+        plotting.plot_anat(
+            registered_nii,
+            display_mode=plane,
+            cut_coords=coords,
+            title="",
+            axes=axes[idx, 1],
+            draw_cross=False,
+            annotate=False,
+            vmin=vmin,
+            vmax=vmax,
+        )
+
+    # Remove all whitespace between subplots and at edges
+    plt.subplots_adjust(hspace=0, wspace=0, left=0, right=1, bottom=0, top=1)
+
+    # Save the figure
+    fig_path = raw_nifti.parent / f"sub-{subid}_ses-{sesid}_vol-{image_index}_comparison.png"
+    plt.savefig(fig_path, dpi=300, bbox_inches="tight")
+    plt.close()
+
+    # Crop the saved figure
+    img = iio.imread(fig_path)
+
+    # Calculate crop dimensions based on the provided proportion
+    h, w = img.shape[:2]
+    crop_h = int(h * crop_proportion)
+    crop_w = int(w * crop_proportion)
+
+    # Crop the image
+    cropped = img[crop_h:-crop_h, crop_w:-crop_w]
+
+    # Save the cropped image
+    iio.imwrite(fig_path, cropped)
+
+
+vols_to_plot = [14, 15, 16, 17, 18, 19, 20, 21, 22]
+for vol in vols_to_plot:
+    print(f"Plotting volume {vol}")
+    raw_nii_path, registered_nii_path = resample_processed_into_raw(vol)
+    make_figure(raw_nii_path, registered_nii_path, vol)

analysis/plot_tedana_t2star.py

@@ -0,0 +1,46 @@
+"""Plot T2* maps from tedana."""
+
+import os
+from glob import glob
+
+import matplotlib.pyplot as plt
+from nilearn import image, plotting
+
+
+if __name__ == "__main__":
+    in_dir = "/cbica/projects/pafin/derivatives/fmriprep"
+    out_dir = "../figures"
+
+    patterns = {
+        "T2star": "sub-*/ses-1/func/*_space-MNI152NLin6Asym_res-2_T2starmap.nii.gz",
+    }
+    for title, pattern in patterns.items():
+        # Get all scalar maps
+        scalar_maps = sorted(glob(os.path.join(in_dir, pattern)))
+        scalar_maps = [f for f in scalar_maps if "PILOT" not in f]
+        print(f"{title}: {len(scalar_maps)}")
+
+        mean_img = image.mean_img(scalar_maps, copy_header=True)
+        sd_img = image.math_img("np.std(img, axis=3)", img=scalar_maps)
+
+        # Plot mean and SD
+        fig, axs = plt.subplots(2, 1, figsize=(10, 5))
+        plotting.plot_stat_map(
+            mean_img,
+            display_mode="z",
+            cut_coords=[-30, -15, 0, 15, 30, 45, 60],
+            title="Mean",
+            axes=axs[0],
+            figure=fig,
+        )
+        plotting.plot_stat_map(
+            sd_img,
+            display_mode="z",
+            cut_coords=[-30, -15, 0, 15, 30, 45, 60],
+            title="Standard Deviation",
+            axes=axs[1],
+            figure=fig,
+        )
+        fig.suptitle(title)
+        plt.savefig(os.path.join(out_dir, f"{title.replace(' ', '_')}.png"))
+        plt.close()