Fast medical imaging cropping

Medical imaging cropping here is a SimpleITK how-to: take a large NIfTI, cut it along z into smaller chunks with RegionOfInterest, and write each chunk so spacing, origin, and direction stay honest. It is not a generic preprocessing explainer. It is not a PYCAD product.

If you meant change voxel spacing → resampling in medical imaging. If you meant the whole preprocess pipeline → what is image preprocessing. If you meant intensity scale → image normalization. If you meant tabular rows (impute / IQR / one-hot) → data preprocessing for machine learning. If you meant MONAI Spacingd / CropForegroundd for tumor seg → preprocessing 3D volumes in MONAI.

3351 already points at this URL for “crop (same grid, smaller FOV).” The job is the cut, not a new spacing.

Why you chunk a volume

A CT can be 400–600 slices. Training a 3D net on the whole array is often the wrong first move: the z-length is not the same from patient to patient, and a convolution on 500 slices will not fit a single GPU. Cutting along z into fixed-depth chunks is one way to get a batch that actually loads. You keep the original xy and the original spacing. You are not resampling.

sitk.RegionOfInterest copies the metadata. A numpy slice that you write back as a new image will invent an origin unless you set it. That is why this snippet goes through SimpleITK instead of array[z:z+n].

The cut

import os
import SimpleITK as sitk
import matplotlib.pyplot as plt

def load_and_process_image_in_chunks(filepath, chunk_size, output_dir):
    os.makedirs(output_dir, exist_ok=True)
    image = sitk.ReadImage(filepath)
    size = image.GetSize()
    print("direction", image.GetDirection())
    print("spacing", image.GetSpacing())
    print("origin", image.GetOrigin())

    for z in range(0, size[2], chunk_size):
        region_size = [size[0], size[1], min(chunk_size, size[2] - z)]
        region_index = [0, 0, z]
        region = sitk.RegionOfInterest(image, region_size, region_index)
        chunk_filename = os.path.join(output_dir, f"chunk_{z}.nii")
        sitk.WriteImage(region, chunk_filename)
        yield region, chunk_filename

def visualize_image(filepath):
    image = sitk.ReadImage(filepath)
    arr = sitk.GetArrayViewFromImage(image)
    for i in range(arr.shape[0]):
        plt.imshow(arr[i, :, :], cmap="gray")
        plt.title(f"slice {i + 1}/{arr.shape[0]}")
        plt.axis("off")
        plt.show()

filepath = "path_to_large_image.nii"
for chunk, name in load_and_process_image_in_chunks(filepath, 50, "output_chunks"):
    print("saved", name)
visualize_image("output_chunks/chunk_0.nii")

SimpleITK indexes are (x, y, z). The numpy view is (z, y, x). Mix those up and the chunk is a slab in the wrong direction. The last chunk is shorter when size[2] is not a multiple of chunk_size — that is what the min is for.

What the metadata is doing

• Spacing. Millimetres between voxel centres. Unchanged by this crop. Changing it is resampling.
• Direction. How the axes sit in the patient. Unchanged.
• Origin. Physical coordinates of index (0, 0, 0). Each chunk gets a new origin at its first voxel. SimpleITK sets that. A raw array write will not.

Install is pip install SimpleITK matplotlib. This is a 2023 PYCAD Team snippet, not a product.

What this page is not

• Not resampling, normalisation, or “what is preprocessing.” Those are the 680 keepers linked above.
• Not MONAI CropForegroundd / Spacingd. That compose is 3050.
• Not a PYCAD crop app or a “fast cropping” SKU.
• Not a join-us / YouTube / services footer. Dropped.

If the crop has to run inside a clinic viewer on a DICOM series, that is the imaging piece. Case studies.