|
| 1 | +# Selecting a region of interest (ROI) |
| 2 | + |
| 3 | +```python |
| 4 | +from pathlib import Path |
| 5 | +import logging |
| 6 | +import numpy as np |
| 7 | +import scipy.ndimage as ndi |
| 8 | +import matplotlib.pyplot as plt |
| 9 | +import nibabel as nib |
| 10 | +import mritk |
| 11 | +``` |
| 12 | + |
| 13 | +First we make sure that the necessary data is downloaded and present in a folder called "gonzo" in the same directory as this script |
| 14 | + |
| 15 | +```python |
| 16 | +mri_data_dir = Path("gonzo") |
| 17 | +``` |
| 18 | + |
| 19 | +If you haven't downloaded he gonzo dataset, you can do so with the following command: |
| 20 | +mritk datasets download gonzo |
| 21 | +```shell |
| 22 | + mritk datasets download gonzo -o gonzo |
| 23 | +``` |
| 24 | + |
| 25 | +If you want to learn more about the gonzo dataset, you can check the command |
| 26 | +```shell |
| 27 | +mritk datasets info gonzo |
| 28 | +``` |
| 29 | +which will point you do the url https://doi.org/10.5281/zenodo.14266867. |
| 30 | + |
| 31 | + |
| 32 | +```python |
| 33 | +# We now load the full T1-weighted image |
| 34 | +t1_path = mri_data_dir / "mri-dataset/mri_dataset/sub-01/ses-01/anat/sub-01_ses-01_T1w.nii.gz" |
| 35 | +t1_data = mritk.MRIData.from_file(t1_path) |
| 36 | +``` |
| 37 | + |
| 38 | +We define a new grid of points in physical space that we want to extract from the original image. |
| 39 | +This grid is defined by the ranges of x, y, and z coordinates and the number of points along each axis. |
| 40 | + |
| 41 | +```python |
| 42 | +xs = np.linspace(0, 70, 80) |
| 43 | +ys = np.linspace(0, 20, 20) |
| 44 | +zs = np.linspace(-40, 90, 110) |
| 45 | +``` |
| 46 | + |
| 47 | +```python |
| 48 | +# Create a 3D grid of points as one long vector |
| 49 | +grid_x, grid_y, grid_z = np.meshgrid(xs, ys, zs, indexing="ij") |
| 50 | +grid_points = np.vstack([grid_x.ravel(), grid_y.ravel(), grid_z.ravel()]).T |
| 51 | +``` |
| 52 | + |
| 53 | +```python |
| 54 | +# We also define a new affine transformation for the extracted piece, which in this case is just the identity matrix. This means that the extracted piece will be in the same physical space as the original image. |
| 55 | +new_affine = np.eye(4) |
| 56 | +new_affine[0, 0] = xs[1] - xs[0] |
| 57 | +new_affine[1, 1] = ys[1] - ys[0] |
| 58 | +new_affine[2, 2] = zs[1] - zs[0] |
| 59 | +new_affine[0, 3] = xs[0] |
| 60 | +new_affine[1, 3] = ys[0] |
| 61 | +new_affine[2, 3] = zs[0] |
| 62 | +``` |
| 63 | + |
| 64 | +We now compute the corresponding voxel indices in the original image for each point in our grid using the affine transformation of the original image. |
| 65 | + |
| 66 | +```python |
| 67 | +vi = mritk.data.physical_to_voxel_indices(grid_points, affine=t1_data.affine) # Shape: (N, 3) |
| 68 | +``` |
| 69 | + |
| 70 | +Finally, we extract the values at the specified voxel indices and reshape them back to the original grid shape. |
| 71 | + |
| 72 | +```python |
| 73 | +v = t1_data.data[tuple(vi.T)] |
| 74 | +v = v.reshape(grid_x.shape) |
| 75 | +``` |
| 76 | + |
| 77 | +```python |
| 78 | +# We save the extracted values as a new NIfTI file with the same affine as the original image. |
| 79 | +piece_t1_data = mritk.MRIData(data=v, affine=new_affine) |
| 80 | +piece_t1_data.save("piece_T1.nii.gz") |
| 81 | +``` |
| 82 | + |
| 83 | +We can now visualize the extracted piece of the T1 image using napari or any other NIfTI viewer. The extracted piece will correspond to the specified grid of points in physical space, allowing us to focus on a specific region of interest (ROI) within the original image. |
| 84 | +We can do this using the command: |
| 85 | +```shell |
| 86 | +mritk napari piece_T1.nii.gz gonzo/mri-dataset/mri_dataset/sub-01/ses-01/anat/sub-01_ses-01_T1w.nii.gz |
| 87 | +``` |
| 88 | + |
| 89 | + |
| 90 | +We can now try to do the same thing for the concentration files, which are stored in a different folder and have a different naming convention. We will loop through all the concentration files, extract the values at the specified voxel indices, |
| 91 | +and save them as new NIfTI files with the same affine as the original image. |
| 92 | + |
| 93 | +```python |
| 94 | +concentration_files = list( |
| 95 | + sorted( |
| 96 | + (mri_data_dir / "mri-processed/mri_processed_data/sub-01/concentrations").glob("sub-01_ses-0*_concentration.nii.gz") |
| 97 | + ) |
| 98 | + ) |
| 99 | +``` |
| 100 | + |
| 101 | +```python |
| 102 | +vs = [] |
| 103 | +``` |
| 104 | + |
| 105 | +```python |
| 106 | +for i, path in enumerate(concentration_files): |
| 107 | + print(path) |
| 108 | + conc_data = mritk.MRIData.from_file(path) |
| 109 | + vi = mritk.data.physical_to_voxel_indices(grid_points, affine=conc_data.affine) # Shape: (N, 3) |
| 110 | + v = conc_data.data[tuple(vi.T)] # Extract values at the specified voxel indices |
| 111 | + v = v.reshape(grid_x.shape) |
| 112 | + vs.append(v) |
| 113 | +``` |
| 114 | + |
| 115 | +```python |
| 116 | +piece_data = mritk.MRIData(data=np.stack(vs), affine=new_affine) |
| 117 | +piece_data.save(f"piece_conc.nii.gz") |
| 118 | +``` |
| 119 | + |
| 120 | +We can visualize the extracted concentration files in napari as well, using the command: |
| 121 | +```shell |
| 122 | +mritk napari piece_conc.nii.gz piece_T1.nii.gz "gonzo/mri-dataset/mri_dataset/sub-01/ses-01/anat/sub-01_ses-01_T1w.nii.gz" |
| 123 | +``` |
| 124 | + |
0 commit comments