Luna Radiology

cli

coregister_volumes

coregister_volumes(input_itk_volume, input_itk_geometry, resample_pixel_spacing, output_dir, order, save_npy)

Resamples and co-registeres all volumes to occupy the same physical coordinates of a reference geometry (given as a itk_volume) and desired voxel size

Parameters:

Name	Type	Description	Default
input_itk_volume	str	path to itk compatible image volume (.mhd, .nrrd, .nii, etc.)	required
input_itk_geometry	str	path to itk compatible image volume (.mhd, .nrrd, .nii, etc.) to use as a reference geometry	required
output_dir	str	output/working directory	required
resample_pixel_spacing	float	voxel size in mm	required
order	int	interpolation order [0-5]	required
save_npy(bool)		whether to also save a numpy file representing the volume	required

Returns:

Name	Type	Description
dict		metadata about function call

Source code in src/luna/radiology/cli/coregister_volumes.py

def coregister_volumes(
    input_itk_volume: str,
    input_itk_geometry: str,
    resample_pixel_spacing: float,
    output_dir: str,
    order: int,
    save_npy: bool,
):
    """Resamples and co-registeres all volumes to occupy the same physical coordinates of a reference geometry (given as a itk_volume) and desired voxel size

    Args:
        input_itk_volume (str): path to itk compatible image volume (.mhd, .nrrd, .nii, etc.)
        input_itk_geometry (str): path to itk compatible image volume (.mhd, .nrrd, .nii, etc.) to use as a reference geometry
        output_dir (str): output/working directory
        resample_pixel_spacing (float): voxel size in mm
        order (int): interpolation order [0-5]
        save_npy(bool): whether to also save a numpy file representing the volume

    Returns:
        dict: metadata about function call
    """
    d_properties = {}

    resample_pixel_spacing = np.full((3), resample_pixel_spacing)

    image_class_object_volume = read_itk_image(
        input_itk_volume, modality=str(Path(input_itk_volume).stem)
    )
    image_class_object_geometry = read_itk_image(
        input_itk_geometry, modality=str(Path(input_itk_geometry).stem)
    )

    _ = interpolate(
        image_class_object_volume,
        resample_pixel_spacing,
        reference_geometry=image_class_object_geometry,
        order=order,
    )

    image_file = image_class_object_volume.export(file_path=output_dir)
    d_properties["itk_volume"] = image_file

    if save_npy:
        np.save(image_file + ".npy", image_class_object_volume.get_voxel_grid())
        d_properties["npy_volume"] = image_file + ".npy"

    # d_properties['volume_size']    = list (image_class_object_volume.size.flatten())
    # d_properties['volume_spacing'] = list (image_class_object_volume.spacing.flatten())
    # d_properties['volume_percentiles'] = list (np.percentile(image_class_object_volume.get_voxel_grid(), np.linspace(0, 100, 11)).flatten())

    output_file = f"{output_dir}/volume_registered_features.parquet"
    pd.DataFrame([d_properties]).to_parquet(output_file)

    d_properties["feature_data"] = output_file

    return d_properties

interpolate(image, resample_spacing, reference_geometry, order=3)

Run interplation

Parameters:

Name	Type	Description	Default
image	imageClass	mirp image class object	required
resample_spacing	ndarray	spacing of resample as a 1D vector	required
reference_geometry	imageClass	output/working directory	required
order	int	interpolation order [0-5]	3

Returns:

Name	Type	Description
imageClass		mirp image class object after resample

Source code in src/luna/radiology/cli/coregister_volumes.py

def interpolate(image, resample_spacing, reference_geometry, order=3):
    """Run interplation

    Args:
        image (imageClass): mirp image class object
        resample_spacing (np.ndarray): spacing of resample as a 1D vector
        reference_geometry (imageClass): output/working directory
        order (int): interpolation order [0-5]

    Returns:
        imageClass: mirp image class object after resample
    """

    image_size = image.size
    image_spacing = image.spacing
    image_voxels = image.voxel_grid
    image_origin = image.origin
    logger.info("Image size=%s, spacing=%s" % (image_size, image_spacing))

    reference_origin = reference_geometry.origin

    reference_offset = (reference_origin - image_origin) / image_spacing
    logger.info("Reference offset=%s" % (reference_offset))

    grid_spacing = resample_spacing / image_spacing

    grid_size = np.ceil(
        np.multiply(
            reference_geometry.size, reference_geometry.spacing / resample_spacing
        )
    )

    logger.info("Grid spacing=%s, size=%s" % (grid_spacing, grid_size))

    map_z, map_y, map_x = np.mgrid[: grid_size[0], : grid_size[1], : grid_size[2]]

    map_z = map_z * grid_spacing[0] + reference_offset[0]
    map_y = map_y * grid_spacing[1] + reference_offset[1]
    map_x = map_x * grid_spacing[2] + reference_offset[2]

    logger.info("Z=%s, Y=%s, Z=%s" % (map_z.shape, map_y.shape, map_x.shape))

    resampled_image = scipy.ndimage.map_coordinates(
        input=image_voxels.astype(np.float32),
        coordinates=np.array([map_z, map_y, map_x], dtype=np.float32),
        order=order,
        mode="nearest",
    )

    print(resampled_image.shape)
    logger.info("Resampled size=%s" % (list(resampled_image.shape)))

    image.set_spacing(resample_spacing)
    image.set_origin(image_origin + (reference_origin - image_origin))
    image.set_voxel_grid(voxel_grid=resampled_image)

    logger.info("New origin=%s" % (image.origin))

    return resampled_image

dicom_to_itk

calculate_normalization(dcms)

Calculates the SUV normalization (g/BQ)

Parameters:

Name	Type	Description	Default
dcms	list[str]	list of dicom paths	required

Returns:

Name	Type	Description
float		Normalization in (g/BQ)

Source code in src/luna/radiology/cli/dicom_to_itk.py

def calculate_normalization(dcms):
    """Calculates the SUV normalization (g/BQ)

    Args:
        dcms (list[str]): list of dicom paths

    Returns:
        float: Normalization in (g/BQ)
    """
    logger.info("About to convert PET volume to SUV units")

    ds = dcmread(dcms[0])

    dose = float(ds.RadiopharmaceuticalInformationSequence[0].RadionuclideTotalDose)
    weight = float(ds.PatientWeight) * 1000

    logger.info(f"Radionuclide dose={dose}, patient weight={weight}")

    norm = weight / dose

    return norm

check_delay_correction(dcms)

Ensures all dicom images were delay corrected to their Aquisition TIme

Parameters:

Name	Type	Description	Default
dcms	list[str]	list of dicom paths	required

Source code in src/luna/radiology/cli/dicom_to_itk.py

def check_delay_correction(dcms):
    """Ensures all dicom images were delay corrected to their Aquisition TIme

    Args:
        dcms (list[str]): list of dicom paths
    """
    for dcm in dcms:
        ds = dcmread(dcm)
        if not ds.DecayCorrection == "START":
            logger.error(
                "check_delay_correction - FAILED - Cannot process a PET volume not constructed of 'START' time delay corrected slices!"
            )
            raise RuntimeError(
                "Cannot process a PET volume not constructed of 'START' time delay corrected slices!"
            )

    logger.info(
        "check_delay_correction - PASSED - All slices were decay corrected to their START AquisitionTime"
    )

check_pet(dcms)

Ensures dicom directory contains PET images

Parameters:

Name	Type	Description	Default
dcms	list[str]	list of dicom paths	required

Returns:

Name	Type	Description
bool		True if all dicoms are PT modality

Source code in src/luna/radiology/cli/dicom_to_itk.py

def check_pet(dcms):
    """Ensures dicom directory contains PET images

    Args:
        dcms (list[str]): list of dicom paths

    Returns:
        bool: True if all dicoms are PT modality
    """
    for dcm in dcms:
        ds = dcmread(dcm)
        if not ds.Modality == "PT":
            logger.warning(
                "check_pet - FAILED - Trying to apply PT corrections to non-PT image, gracefully skipping!"
            )
            return False

    logger.info("check_pet - PASSED - These are PT images")

    return True

convert_pet_volume_to_suv(dicom_urlpath, input_volume)

Renormalizes a PET ITK volume to SUVs, saves a new volume in-place

Parameters:

Name	Type	Description	Default
dicom_urlpath	str	path to matching dicom series within a folder	required
input_volume	str	path to PT volume	required

Source code in src/luna/radiology/cli/dicom_to_itk.py

def convert_pet_volume_to_suv(dicom_urlpath, input_volume):
    """Renormalizes a PET ITK volume to SUVs, saves a new volume in-place

    Args:
        dicom_urlpath (str): path to matching dicom series within a folder
        input_volume (str): path to PT volume
    """
    dcms = list(Path(dicom_urlpath).glob("*.dcm"))

    if check_pet(dcms):
        check_delay_correction(dcms)

        norm = calculate_normalization(dcms)
        logger.info(f"Calculated SUV normalization factor={norm}")

        image, header = medpy.io.load(input_volume)

        image = image * norm
        logger.info(f"SUM SUV={image.sum()}")

        medpy.io.save(image, input_volume, hdr=header)

        logger.info(f"Saved normalized volume: {input_volume}")

dicom_to_itk(dicom_urlpath, output_urlpath, itk_image_type, itk_c_type, convert_to_suv=False)

Generate an ITK compatible image from a dicom series/folder

Parameters:

Name	Type	Description	Default
dicom_urlpath	str	path to dicom series within a folder	required
output_urlpath	str	output/working directory	required
itk_image_type	str	ITK volume image type to output (mhd, nrrd, nii, etc.)	required
itk_c_type	str	pixel (C) type for pixels, e.g. float or unsigned short	required

Returns:

Name	Type	Description
dict		metadata about function call

Source code in src/luna/radiology/cli/dicom_to_itk.py

def dicom_to_itk(
    dicom_urlpath, output_urlpath, itk_image_type, itk_c_type, convert_to_suv=False
):
    """Generate an ITK compatible image from a dicom series/folder

    Args:
        dicom_urlpath (str): path to dicom series within a folder
        output_urlpath (str): output/working directory
        itk_image_type (str): ITK volume image type to output (mhd, nrrd, nii, etc.)
        itk_c_type (str): pixel (C) type for pixels, e.g. float or unsigned short

    Returns:
        dict: metadata about function call
    """
    PixelType = itk.ctype(itk_c_type)
    ImageType = itk.Image[PixelType, 3]

    namesGenerator = itk.GDCMSeriesFileNames.New()
    namesGenerator.SetUseSeriesDetails(True)
    namesGenerator.AddSeriesRestriction("0008|0021")
    namesGenerator.SetGlobalWarningDisplay(False)
    namesGenerator.SetDirectory(dicom_urlpath)

    seriesUIDs = namesGenerator.GetSeriesUIDs()
    num_dicoms = len(seriesUIDs)

    if num_dicoms < 1:
        logger.warning("No DICOMs in: " + dicom_urlpath)
        return None

    logger.info(
        "The directory {} contains {} DICOM Series".format(
            dicom_urlpath, str(num_dicoms)
        )
    )

    n_slices = 0
    volume = {}
    for uid in seriesUIDs:
        logger.info("Reading: " + uid)
        fileNames = namesGenerator.GetFileNames(uid)
        if len(fileNames) < 1:
            continue

        n_slices = len(fileNames)

        reader = itk.ImageSeriesReader[ImageType].New()
        dicomIO = itk.GDCMImageIO.New()
        reader.SetImageIO(dicomIO)
        reader.SetFileNames(fileNames)
        reader.ForceOrthogonalDirectionOff()

        writer = itk.ImageFileWriter[ImageType].New()

        outFileName = os.path.join(
            output_urlpath, uid + "_volumetric_image." + itk_image_type
        )
        writer.SetFileName(outFileName)
        writer.UseCompressionOn()
        writer.SetInput(reader.GetOutput())
        logger.info("Writing: " + outFileName)
        writer.Update()

        img, _ = medpy.io.load(outFileName)
        volume[outFileName] = np.prod(img.shape)

    if convert_to_suv:
        convert_pet_volume_to_suv(dicom_urlpath, outFileName)

    path = next(Path(dicom_urlpath).glob("*.dcm"))
    ds = dcmread(path)

    # If there are multiple seriesUIDs in a single DICOM dir, return
    # the largest one by volume in the output properties
    outFileName = max(volume, key=volume.get)

    # Prepare metadata and commit
    properties = {
        "itk_volume": outFileName,
        "num_slices": n_slices,
        "segment_keys": {
            "radiology_patient_name": str(ds.PatientName),
            "radiology_accession_number": str(ds.AccessionNumber),
            "radiology_series_instance_uuid": str(ds.SeriesInstanceUID),
            "radiology_series_number": str(ds.SeriesNumber),
            "radiology_modality": str(ds.Modality),
        },
    }

    return properties

extract_radiomics

extract_radiomics_multiple_labels(itk_volume_urlpath, itk_labels_urlpath, lesion_indices, pyradiomics_config, output_urlpath, check_geometry_strict=False, enable_all_filters=False)

Extract radiomics using pyradiomics, with additional checking of the input geometry match, resulting in a feature csv file

Parameters:

Name	Type	Description	Default
itk_volume_urlpath	str	path to itk compatible image volume (.mhd, .nrrd, .nii, etc.)	required
itk_labels_urlpath	str	path to itk compatible label volume (.mha)	required
output_urlpath	str	output/working directory	required
lesion_indices	List[int]	list of lesion indicies (label values) to process	required
pyradiomics_config	dict	keyword arguments to pass to pyradiomics featureextractor	required
check_geometry_strict	str	enforce strict match of the input geometries (otherwise, possible mismatch is silently ignored)	False
enable_all_filters	bool	turns on all image filters	False

Returns:

Name	Type	Description
dict		metadata about function call

Source code in src/luna/radiology/cli/extract_radiomics.py

def extract_radiomics_multiple_labels(
    itk_volume_urlpath: str,
    itk_labels_urlpath: str,
    lesion_indices: List[int],
    pyradiomics_config: dict,
    output_urlpath: str,
    check_geometry_strict: bool = False,
    enable_all_filters: bool = False,
):
    """Extract radiomics using pyradiomics, with additional checking of the input geometry match, resulting in a feature csv file

    Args:
        itk_volume_urlpath (str): path to itk compatible image volume (.mhd, .nrrd, .nii, etc.)
        itk_labels_urlpath (str): path to itk compatible label volume (.mha)
        output_urlpath (str): output/working directory
        lesion_indices (List[int]): list of lesion indicies (label values) to process
        pyradiomics_config (dict): keyword arguments to pass to pyradiomics featureextractor
        check_geometry_strict (str): enforce strict match of the input geometries (otherwise, possible mismatch is silently ignored)
        enable_all_filters (bool): turns on all image filters

    Returns:
        dict: metadata about function call
    """

    if type(lesion_indices) == int:
        lesion_indices = [lesion_indices]

    image, image_header = medpy.io.load(itk_volume_urlpath)

    if Path(itk_labels_urlpath).is_dir():
        label_path_list = [str(path) for path in Path(itk_labels_urlpath).glob("*")]
    elif Path(itk_labels_urlpath).is_file():
        label_path_list = [itk_labels_urlpath]
    else:
        raise RuntimeError("Issue with detecting label format")

    available_labels = set()
    for label_path in label_path_list:
        label, label_header = medpy.io.load(label_path)

        available_labels.update

        available_labels.update(np.unique(label))

        logger.info(f"Checking {label_path}")

        if (
            check_geometry_strict
            and not image_header.get_voxel_spacing() == label_header.get_voxel_spacing()
        ):
            raise RuntimeError(
                f"Voxel spacing mismatch, image.spacing={image_header.get_voxel_spacing()}, label.spacing={label_header.get_voxel_spacing()}"
            )

        if not image.shape == label.shape:
            raise RuntimeError(
                f"Shape mismatch: image.shape={image.shape}, label.shape={label.shape}"
            )

    df_result = pd.DataFrame()

    for lesion_index in available_labels.intersection(lesion_indices):
        extractor = featureextractor.RadiomicsFeatureExtractor(
            label=lesion_index, **pyradiomics_config
        )

        if enable_all_filters:
            extractor.enableAllImageTypes()

        for label_path in label_path_list:
            result = extractor.execute(itk_volume_urlpath, label_path)

            result["lesion_index"] = lesion_index

            df_result = pd.concat([df_result, pd.Series(result).to_frame()], axis=1)

    output_filename = os.path.join(output_urlpath, "radiomics.csv")

    df_result.T.to_csv(output_filename, index=False)

    logger.info(df_result.T)

    properties = {
        "feature_csv": output_filename,
        "lesion_indices": lesion_indices,
    }

    return properties

extract_voxels

extract_voxels(input_itk_volume, output_dir)

Save a numpy file from a given ITK volume

Parameters:

Name	Type	Description	Default
input_itk_volume	str	path to itk compatible image volume (.mhd, .nrrd, .nii, etc.)	required
output_dir	str	output/working directory	required

Returns:

Name	Type	Description
dict		metadata about function call

Source code in src/luna/radiology/cli/extract_voxels.py

def extract_voxels(input_itk_volume, output_dir):
    """Save a numpy file from a given ITK volume

    Args:
        input_itk_volume (str): path to itk compatible image volume (.mhd, .nrrd, .nii, etc.)
        output_dir (str): output/working directory

    Returns:
        dict: metadata about function call
    """
    file_stem = Path(input_itk_volume).stem
    # file_ext = Path(input_itk_volume).suffix

    outFileName = os.path.join(output_dir, file_stem + ".npy")

    image, header = medpy.io.load(input_itk_volume)

    np.save(outFileName, image)

    logger.info(f"Extracted voxels of shape {image.shape}")

    # Prepare metadata and commit
    properties = {
        "npy_volume": outFileName,
    }

    return properties

generate_threshold_mask

generate_threshold_mask(input_itk_volume, threshold, area_closing_radius, expansion_radius, save_npy, output_dir)

Use: https://scikit-image.org/docs/dev/api/skimage.segmentation.html#skimage.segmentation.expand_labels And: https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.distance_transform_edt.html

Source code in src/luna/radiology/cli/generate_threshold_mask.py

def generate_threshold_mask(
    input_itk_volume,
    threshold,
    area_closing_radius,
    expansion_radius,
    save_npy,
    output_dir,
):
    """
    Use: https://scikit-image.org/docs/dev/api/skimage.segmentation.html#skimage.segmentation.expand_labels
    And: https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.distance_transform_edt.html
    """
    d_properties = {}

    image_class_object_volume = read_itk_image(input_itk_volume)

    voxel_size = image_class_object_volume.spacing[0]
    voxel_grid = image_class_object_volume.get_voxel_grid()

    voxel_grid = np.where(voxel_grid > threshold, 1, 0).astype(bool)
    logger.info(f"Initial mask sum = {voxel_grid.sum()}")

    logger.info("Applying area closing....")
    closing_surface_area = 12 * area_closing_radius**2 / voxel_size
    voxel_grid = remove_small_holes(voxel_grid, closing_surface_area)
    logger.info(
        f"area_closing(SA={closing_surface_area}) mask sum = {voxel_grid.sum()}"
    )

    voxel_edt = distance_transform_edt(-voxel_grid.astype(int) + 1, sampling=voxel_size)
    voxel_mask = np.where(voxel_edt < expansion_radius, 1, 0)

    image_class_object_volume.set_voxel_grid(voxel_grid=voxel_edt)

    image_file = image_class_object_volume.export(file_path=output_dir)
    d_properties["itk_labels"] = image_file

    if save_npy:
        np.save(image_file + ".npy", voxel_mask.astype(np.uint8))
        d_properties["npy_labels"] = image_file + ".npy"
        np.save(image_file + ".edt.npy", voxel_edt.astype(np.float32))
        d_properties["npy_edt_labels"] = image_file + ".edt.npy"

    d_properties["n_mask_voxels"] = voxel_grid.sum()

    output_file = f"{output_dir}/volume_registered_features.parquet"
    pd.DataFrame([d_properties]).to_parquet(output_file)

    d_properties["segment_keys"] = {"radiology_submodality": "DERIVED"}

    d_properties["feature_data"] = output_file

    return d_properties

match_metadata

match_metadata(dicom_tree_folder, input_itk_labels, output_dir)

Generate an ITK compatible image from a dicom series/folder

Parameters:

Name	Type	Description	Default
dicom_tree_folder	str	path to root/tree dicom folder that may contained multiple dicom series	required
input_itk_labels	str	path to itk compatible label volume (.mha)	required
output_dir	str	output/working directory	required

Returns:

Name	Type	Description
dict		metadata about function call

Source code in src/luna/radiology/cli/match_metadata.py

def match_metadata(dicom_tree_folder: str, input_itk_labels: str, output_dir: str):
    """Generate an ITK compatible image from a dicom series/folder

    Args:
        dicom_tree_folder (str): path to root/tree dicom folder that may contained multiple dicom series
        input_itk_labels (str): path to itk compatible label volume (.mha)
        output_dir (str): output/working directory

    Returns:
        dict: metadata about function call
    """
    dicom_folders = set()
    for dicom in Path(dicom_tree_folder).rglob("*.dcm"):
        dicom_folders.add(dicom.parent)

    label, _ = medpy.io.load(input_itk_labels)
    found_dicom_paths = set()
    for dicom_folder in dicom_folders:
        n_slices_dcm = len(list((dicom_folder).glob("*.dcm")))
        if label.shape[2] == n_slices_dcm:
            found_dicom_paths.add(dicom_folder)

    logger.info(found_dicom_paths)
    if not len(found_dicom_paths) > 0:
        raise RuntimeError("Could not find matching scans!")

    matched = None
    for found_dicom_path in found_dicom_paths:
        path = next(found_dicom_path.glob("*.dcm"))
        ds = dcmread(path)
        if not ds.Modality == "CT":
            continue
        print(
            "Matched: z=",
            label.shape[2],
            found_dicom_path,
            ds.PatientName,
            ds.AccessionNumber,
            ds.SeriesInstanceUID,
            ds.StudyDescription,
            ds.SeriesDescription,
            ds.SliceThickness,
        )
        matched = found_dicom_path

    if matched is None:
        raise RuntimeError("Could not find matching CT!")

    properties = {
        "dicom_folder": str(matched),
        "zdim": label.shape[2],
    }

    return properties

randomize_contours

randomize_contours(input_itk_volume, input_itk_labels, resample_pixel_spacing, resample_smoothing_beta, output_dir)

Generate randomize contours given and image, label after resampling using MIRP processing library

First, we interpolate the inputs to a isotropic spacing, then get the supervoxel pertubations using the MIRP methods, as defined in:

Zwanenburg, A., Leger, S., Agolli, L. et al. Assessing robustness of radiomic features by image perturbation. Sci Rep 9, 614 (2019). https://doi.org/10.1038/s41598-018-36938-4

Parameters:

Name	Type	Description	Default
input_itk_volume	str	path to itk compatible image volume (.mhd, .nrrd, .nii, etc.)	required
input_itk_labels	str	path to itk compatible label volume (.mha)	required
output_dir	str	output/working directory	required
resample_pixel_spacing	float	voxel size in mm	required
resample_smoothing_beta	float	smoothing beta of gaussian filter	required

Returns:

Name	Type	Description
dict		metadata about function call

Source code in src/luna/radiology/cli/randomize_contours.py

def randomize_contours(
    input_itk_volume: str,
    input_itk_labels: str,
    resample_pixel_spacing: float,
    resample_smoothing_beta: float,
    output_dir: str,
):
    """Generate randomize contours given and image, label after resampling using MIRP processing library

    First, we interpolate the inputs to a isotropic spacing, then get the supervoxel pertubations using the MIRP methods,
    as defined in:

    Zwanenburg, A., Leger, S., Agolli, L. et al. Assessing robustness of radiomic features by image perturbation. Sci Rep 9, 614 (2019). https://doi.org/10.1038/s41598-018-36938-4

    Args:
        input_itk_volume (str): path to itk compatible image volume (.mhd, .nrrd, .nii, etc.)
        input_itk_labels (str): path to itk compatible label volume (.mha)
        output_dir (str): output/working directory
        resample_pixel_spacing (float): voxel size in mm
        resample_smoothing_beta (float): smoothing beta of gaussian filter

    Returns:
        dict: metadata about function call
    """
    logger.info("Hello, processing %s, %s", input_itk_volume, input_itk_labels)
    settings = Settings()

    print(settings)

    resample_pixel_spacing = np.full((3), resample_pixel_spacing)

    settings.img_interpolate.new_spacing = resample_pixel_spacing
    settings.roi_interpolate.new_spacing = resample_pixel_spacing
    settings.img_interpolate.smoothing_beta = resample_smoothing_beta

    # Read
    image_class_object = read_itk_image(input_itk_volume, "CT")
    roi_class_object_list = read_itk_segmentation(input_itk_labels)

    # Crop for faster interpolation
    image_class_object, roi_class_object_list = crop_image(
        img_obj=image_class_object,
        roi_list=roi_class_object_list,
        boundary=50.0,
        z_only=True,
    )

    # Interpolation
    image_class_object = interpolate_image(
        img_obj=image_class_object, settings=settings
    )
    roi_class_object_list = interpolate_roi(
        img_obj=image_class_object, roi_list=roi_class_object_list, settings=settings
    )

    # Export
    image_file = image_class_object.export(file_path=f"{output_dir}/main_image")

    # ROI processing
    roi_class_object = combine_all_rois(
        roi_list=roi_class_object_list, settings=settings
    )
    label_file = roi_class_object.export(
        img_obj=image_class_object, file_path=f"{output_dir}/main_label"
    )

    roi_class_object_list, svx_class_object_list = randomise_roi_contours(
        img_obj=image_class_object, roi_list=roi_class_object_list, settings=settings
    )

    roi_supervoxels = combine_all_rois(
        roi_list=svx_class_object_list, settings=settings
    )
    voxels_file = roi_supervoxels.export(
        img_obj=image_class_object, file_path=f"{output_dir}/supervoxels"
    )

    for roi in combine_pertubation_rois(
        roi_list=roi_class_object_list, settings=settings
    ):
        if "COMBINED" in roi.name:
            roi.export(
                img_obj=image_class_object, file_path=f"{output_dir}/pertubations"
            )

    print(image_file, label_file)

    # Construct return dicts
    properties = {
        "itk_volume": image_file,
        "itk_labels": label_file,
        "mirp_pertubations": f"{output_dir}/pertubations",
        "mirp_supervoxels": voxels_file,
    }

    return properties

window_volume

window_volume(input_itk_volume, output_dir, low_level, high_level)

Applies a window function (clipping) to an input itk volume, outputs windowed volume

Parameters:

Name	Type	Description	Default
input_itk_volume	str	path to itk compatible image volume (.mhd, .nrrd, .nii, etc.)	required
output_dir	str	output/working directory	required
low_level	float	lower bound of clipping operation	required
high_level	float	higher bound of clipping operation	required

Returns:

Name	Type	Description
dict		metadata about function call

Source code in src/luna/radiology/cli/window_volume.py

def window_volume(
    input_itk_volume: str, output_dir: str, low_level: float, high_level: float
):
    """Applies a window function (clipping) to an input itk volume, outputs windowed volume

    Args:
        input_itk_volume (str): path to itk compatible image volume (.mhd, .nrrd, .nii, etc.)
        output_dir (str): output/working directory
        low_level (float): lower bound of clipping operation
        high_level (float): higher bound of clipping operation

    Returns:
        dict: metadata about function call
    """
    file_stem = Path(input_itk_volume).stem
    file_ext = Path(input_itk_volume).suffix

    outFileName = os.path.join(output_dir, file_stem + ".windowed" + file_ext)

    logger.info("Applying window [%s,%s]", low_level, high_level)

    image, header = medpy.io.load(input_itk_volume)
    image = np.clip(image, low_level, high_level)
    medpy.io.save(image, outFileName, header)
    # Prepare metadata and commit
    properties = {
        "itk_volume": outFileName,
    }

    return properties

mirp

Created on April 27, 2021

@author: pashaa@mskcc.org

imageClass

USED -> IMAGE CLASS

ImageClass

Source code in src/luna/radiology/mirp/imageClass.py

class ImageClass:
    # Class for image volumes

    def __init__(
        self,
        voxel_grid,
        origin,
        spacing,
        orientation,
        modality=None,
        spat_transform="base",
        no_image=False,
        metadata=None,
        slice_table=None,
    ):
        # Set details regarding voxel orientation and such
        self.origin = np.array(origin)
        self.orientation = np.array(orientation)

        self.spat_transform = (
            spat_transform  # Signifies whether the current image is a base image or not
        )
        self.slice_table = slice_table

        # The spacing, the affine matrix and its inverse are set using the set_spacing method.
        self.spacing = None
        self.m_affine = None
        self.m_affine_inv = None

        # Set voxel spacing. This also set the affine matrix and its inverse.
        self.set_spacing(new_spacing=np.array(spacing))

        # Image name
        self.name = None

        # Initialise voxel grid dependent parameters
        self.isEncoded_voxel_grid = None
        self.voxel_grid = None
        self.size = None
        self.dtype_name = None

        # Interpolation settings
        self.interpolated = False
        self.interpolation_algorithm = None

        # Bin settings
        self.binned = False
        self.bin_width = None

        # Discretisation settings
        self.discretised = False
        self.discretisation_algorithm = None
        self.discretisation_settings = None

        # Noise addition parameters
        self.noise = -1.0
        self.noise_iter = 0

        # Translation parameters
        self.transl_fraction_x = 0.0
        self.transl_fraction_y = 0.0
        self.transl_fraction_z = 0.0

        # Rotation parameters
        self.rotation_angle = 0.0

        # Set voxel grid and image
        if not no_image:
            self.is_missing = False
            self.set_voxel_grid(voxel_grid=voxel_grid)
        else:
            self.is_missing = True

        # Set metadata and a list of update tags
        self.metadata: Union[FileDataset, None] = metadata
        self.as_parametric_map = False

        # Image modality
        if modality is None and metadata is not None:
            # Set imaging modality using metadata
            self.modality = self.get_metadata(
                tag=(0x0008, 0x0060), tag_type="str"
            )  # Imaging modality
        elif modality is None:
            self.modality = "GENERIC"
        else:
            self.modality = modality

        # Normalisation flags.
        self.is_normalised = False

    def copy(self, drop_image=False):
        # Creates a new copy of the image object
        img_copy = copy.deepcopy(self)

        if drop_image:
            img_copy.drop_image()

        return img_copy

    def show(self, img_slice):
        import pylab

        if self.is_missing:
            return

        pylab.imshow(self.get_voxel_grid()[img_slice, :, :], cmap=pylab.cm.bone)
        pylab.show()

    def set_origin(self, new_origin):
        self.origin = new_origin

    def set_spacing(self, new_spacing):
        # Update spacing
        self.spacing: np.ndarray = new_spacing

        # Recompute the affine matrices
        m_affine = np.zeros((3, 3), dtype=np.float)

        # z-coordinates
        m_affine[:, 0] = self.spacing[0] * np.array(
            [self.orientation[0], self.orientation[1], self.orientation[2]]
        )

        # y-coordinates
        m_affine[:, 1] = self.spacing[1] * np.array(
            [self.orientation[3], self.orientation[4], self.orientation[5]]
        )

        # x-coordinates
        m_affine[:, 2] = self.spacing[2] * np.array(
            [self.orientation[6], self.orientation[7], self.orientation[8]]
        )

        self.m_affine = m_affine
        self.m_affine_inv = np.linalg.inv(self.m_affine)

    def set_voxel_grid(self, voxel_grid):
        """Sets voxel grid"""

        # Determine size
        self.size = np.array(voxel_grid.shape)
        self.dtype_name = voxel_grid.dtype.name

        # Encode voxel grid
        self.encode_voxel_grid(voxel_grid=voxel_grid)

    # Return None for missing images
    def get_voxel_grid(self) -> np.ndarray:
        """Return the voxel grid as a ndarray"""
        if self.is_missing:
            return None

        if self.isEncoded_voxel_grid:
            # Decode voxel grid (typically roi)
            decoded_voxel = np.zeros(np.prod(self.size), dtype=np.bool)

            # Check if the voxel grid contains values
            if self.voxel_grid is not None:
                decode_zip = copy.deepcopy(self.voxel_grid)

                for ii, jj in decode_zip:
                    decoded_voxel[ii : jj + 1] = True

            # Shape into correct form
            decoded_voxel = decoded_voxel.reshape(self.size)

            return decoded_voxel
        else:
            return self.voxel_grid

    def encode_voxel_grid(self, voxel_grid):
        """Performs run length encoding of the voxel grid"""

        # Determine whether the voxel grid should be encoded (only True for boolean data types; typically roi)
        if self.dtype_name == "bool":
            # Run length encoding for "True"
            rle_end = np.array(
                np.append(
                    np.where(voxel_grid.ravel()[1:] != voxel_grid.ravel()[:-1]),
                    np.prod(self.size) - 1,
                )
            )
            rle_start = np.cumsum(np.append(0, np.diff(np.append(-1, rle_end))))[:-1]
            rle_val = voxel_grid.ravel()[rle_start]

            # Check whether the voxel grid is empty (consists of 0s)
            if np.all(~rle_val):
                self.voxel_grid = None
                self.isEncoded_voxel_grid = True
            else:
                # Select only True values entries for further compression
                rle_start = rle_start[rle_val]
                rle_end = rle_end[rle_val]

                # Create zip
                self.voxel_grid = zip(rle_start, rle_end)
                self.isEncoded_voxel_grid = True
        else:
            self.voxel_grid = voxel_grid
            self.isEncoded_voxel_grid = False

    def decode_voxel_grid(self):
        """Performs run length decoding of the voxel grid and converts it to a numpy array"""
        if self.dtype_name == "bool" and self.isEncoded_voxel_grid:
            decoded_voxel = np.zeros(np.prod(self.size), dtype=np.bool)

            # Check if the voxel grid contains values
            if self.voxel_grid is not None:
                decode_zip = copy.deepcopy(self.voxel_grid)
                for ii, jj in decode_zip:
                    decoded_voxel[ii : jj + 1] = True

            # Set shape to original grid
            decoded_voxel = decoded_voxel.reshape(self.size)

            # Update self.voxel_grid and isEncoded_voxel_grid tags
            self.voxel_grid = decoded_voxel
            self.isEncoded_voxel_grid = False

    def decimate(self, by_slice):
        """
        Decimates image voxel grid by removing every second element
        :param by_slice:
        :return:
        """

        # Skip for missing images
        if self.is_missing:
            return

        # Get the voxel grid
        img_voxel_grid = self.get_voxel_grid()

        # Update the voxel grid
        if by_slice:
            # Drop every second pixel
            img_voxel_grid = img_voxel_grid[
                :, slice(None, None, 2), slice(None, None, 2)
            ]

            # Update voxel spacing
            self.spacing[[1, 2]] *= 2.0
        else:
            # Drop every second voxel
            img_voxel_grid = img_voxel_grid[
                slice(None, None, 2), slice(None, None, 2), slice(None, None, 2)
            ]

            # Update voxel spacing
            self.spacing *= 2.0

        # Update voxel grid. This also updates the size attribute.
        self.set_voxel_grid(voxel_grid=img_voxel_grid)

    def interpolate(self, by_slice, settings):
        """Performs interpolation of the image volume"""
        from luna.radiology.mirp.imageProcess import (  # aauker: Circular import
            gaussian_preprocess_filter,
            interpolate_to_new_grid,
        )

        # Skip for missing images
        if self.is_missing:
            return

        # Local interpolation constants
        if None not in settings.img_interpolate.new_spacing:
            iso_spacing = settings.img_interpolate.new_spacing[0]
            new_spacing = np.array(
                [iso_spacing, iso_spacing, iso_spacing]
            )  # Desired spacing in mm
        elif type(settings.img_interpolate.new_non_iso_spacing) in [list, tuple]:
            if None not in settings.img_interpolate.new_non_iso_spacing:
                non_iso_spacing = settings.img_interpolate.new_non_iso_spacing
                new_spacing = np.array(non_iso_spacing)
            else:
                new_spacing = self.spacing
        else:
            new_spacing = self.spacing

        print(f"Interpolating main image to {new_spacing}")

        # Read additional details
        order = (
            settings.img_interpolate.spline_order
        )  # Order of multidimensional spline filter (0=nearest neighbours, 1=linear, 3=cubic)
        interpolate_flag = (
            settings.img_interpolate.interpolate
        )  # Whether to interpolate or not

        # Set spacing for interpolation across slices to the original spacing in case interpolation is only conducted within the slice
        if by_slice:
            new_spacing[0] = self.spacing[0]

        # Image translation
        translate_z = 0  # settings.vol_adapt.translate_z[0]
        translate_y = 0  # settings.vol_adapt.translate_y[0]
        translate_x = 0  # settings.vol_adapt.translate_x[0]

        # Convert to [0.0, 1.0] range
        translate_x = translate_x - np.floor(translate_x)
        translate_y = translate_y - np.floor(translate_y)
        translate_z = translate_z - np.floor(translate_z)
        trans_vec = np.array([translate_z, translate_y, translate_x])

        # Add translation fractions
        self.transl_fraction_x = translate_x
        self.transl_fraction_y = translate_y
        self.transl_fraction_z = translate_z

        # Skip if translation in both directions is 0.0
        if (
            translate_x == 0.0
            and translate_y == 0.0
            and translate_z == 0.0
            and not interpolate_flag
        ):
            return None

        # Check if pre-processing is required
        if settings.img_interpolate.anti_aliasing:
            self.set_voxel_grid(
                voxel_grid=gaussian_preprocess_filter(
                    orig_vox=self.get_voxel_grid(),
                    orig_spacing=self.spacing,
                    sample_spacing=new_spacing,
                    param_beta=settings.img_interpolate.smoothing_beta,
                    mode="nearest",
                    by_slice=by_slice,
                )
            )

        # Interpolate image and positioning
        (
            self.size,
            sample_spacing,
            upd_voxel_grid,
            grid_origin,
        ) = interpolate_to_new_grid(
            orig_dim=self.size,
            orig_spacing=self.spacing,
            orig_vox=self.get_voxel_grid(),
            sample_spacing=new_spacing,
            translation=trans_vec,
            order=order,
            mode="nearest",
            align_to_center=True,
        )

        # Update origin before spacing, because computing the origin requires the original affine matrix.
        self.origin = self.origin + np.dot(self.m_affine, np.transpose(grid_origin))

        # Update spacing and affine matrix.
        self.set_spacing(sample_spacing)

        # Round intensities in case of modalities with inherently discretised intensities
        if (self.modality == "CT") and (self.spat_transform == "base"):
            upd_voxel_grid = np.round(upd_voxel_grid)
        elif (self.modality == "PT") and (self.spat_transform == "base"):
            upd_voxel_grid[upd_voxel_grid < 0.0] = 0.0

        # Set interpolation
        self.interpolated = True

        # Set interpolation algorithm
        if order == 0:
            self.interpolation_algorithm = "nnb"
        elif order == 1:
            self.interpolation_algorithm = "lin"
        elif order > 1:
            self.interpolation_algorithm = "si" + str(order)

        if settings.img_interpolate.bin:
            self.binned = True
            self.bin_width = settings.img_interpolate.bin_width
            self.bins = np.arange(-1000, 1000, settings.img_interpolate.bin_width)
            upd_voxel_grid = np.digitize(upd_voxel_grid, self.bins).astype(np.float)

        # Set voxel grid
        self.set_voxel_grid(voxel_grid=upd_voxel_grid)

    def add_noise(self, noise_level, noise_iter):
        """Adds Gaussian noise to the image volume
        noise_level: standard deviation of image noise present"""

        # Add noise iteration number
        self.noise_iter = noise_iter

        # Skip for missing images
        if self.is_missing:
            return

        # Skip for invalid noise levels
        if noise_level is None:
            return
        if np.isnan(noise_level) or noise_level < 0.0:
            return

        # Add Gaussian noise to image
        voxel_grid = self.get_voxel_grid()
        voxel_grid += np.random.normal(loc=0.0, scale=noise_level, size=self.size)

        # Check for corrections due to image modality
        if self.spat_transform == "base":
            # Round CT values to the nearest integer
            if self.modality == "CT":
                voxel_grid = np.round(a=voxel_grid, decimals=0)

            # Set minimum PT to 0.0
            if self.modality == "PT":
                voxel_grid[voxel_grid < 0.0] = 0.0

        # Set noise level in image
        self.noise = noise_level

        self.set_voxel_grid(voxel_grid=voxel_grid)

    def saturate(self, intensity_range, fill_value=None):
        """
        Saturate image intensities using an intensity range
        :param intensity_range: range of intensity values
        :param fill_value: fill value for out-of-range intensities. If None, the upper and lower ranges are used
        :return:
        """
        # Skip for missing images
        if self.is_missing:
            return

        intensity_range = np.array(copy.deepcopy(intensity_range))

        if np.any(~np.isnan(intensity_range)):
            # Get voxel grid
            voxel_grid = self.get_voxel_grid()

            # Lower boundary
            if not np.isnan(intensity_range[0]):
                if fill_value is None:
                    voxel_grid[voxel_grid < intensity_range[0]] = intensity_range[0]
                else:
                    voxel_grid[voxel_grid < intensity_range[0]] = fill_value[0]

            # Upper boundary
            if not np.isnan(intensity_range[1]):
                if fill_value is None:
                    voxel_grid[voxel_grid > intensity_range[1]] = intensity_range[1]
                else:
                    voxel_grid[voxel_grid > intensity_range[1]] = fill_value[1]

            # Set the updated voxel grid
            self.set_voxel_grid(voxel_grid=voxel_grid)

    def normalise_intensities(
        self, norm_method="none", intensity_range=None, saturation_range=None, mask=None
    ):
        """
        Normalises image intensities
        :param norm_method: string defining the normalisation method. Should be one of "none", "range", "standardisation"
        :param intensity_range: range of intensities for normalisation
        :return:
        """

        # Skip for missing images
        if self.is_missing:
            return

        if intensity_range is None:
            intensity_range = [np.nan, np.nan]

        if mask is None:
            mask = np.ones(self.size, dtype=np.bool)
        else:
            mask = mask.astype(np.bool)

        if np.sum(mask) == 0:
            mask = np.ones(self.size, dtype=np.bool)

        if saturation_range is None:
            saturation_range = [np.nan, np.nan]

        if norm_method == "none":
            return

        elif norm_method == "range":
            # Normalisation to [0, 1] range using fixed intensities.

            # Get voxel grid
            voxel_grid = self.get_voxel_grid()

            # Find maximum and minimum intensities
            if np.isnan(intensity_range[0]):
                min_int = np.min(voxel_grid[mask])
            else:
                min_int = intensity_range[0]

            if np.isnan(intensity_range[1]):
                max_int = np.max(voxel_grid[mask])
            else:
                max_int = intensity_range[1]

            # Normalise by range
            if not max_int == min_int:
                voxel_grid = (voxel_grid - min_int) / (max_int - min_int)
            else:
                voxel_grid = voxel_grid - min_int

            # Update the voxel grid
            self.set_voxel_grid(voxel_grid=voxel_grid)

            self.is_normalised = True

        elif norm_method == "relative_range":
            # Normalisation to [0, 1]-ish range using relative intensities.

            # Get voxel grid
            voxel_grid = self.get_voxel_grid()

            min_int_rel = 0.0
            if not np.isnan(intensity_range[0]):
                min_int_rel = intensity_range[0]

            max_int_rel = 1.0
            if not np.isnan(intensity_range[1]):
                max_int_rel = intensity_range[1]

            # Compute minimum and maximum intensities.
            value_range = [np.min(voxel_grid[mask]), np.max(voxel_grid[mask])]
            min_int = value_range[0] + min_int_rel * (value_range[1] - value_range[0])
            max_int = value_range[0] + max_int_rel * (value_range[1] - value_range[0])

            # Normalise by range
            if not max_int == min_int:
                voxel_grid = (voxel_grid - min_int) / (max_int - min_int)
            else:
                voxel_grid = voxel_grid - min_int

            # Update the voxel grid
            self.set_voxel_grid(voxel_grid=voxel_grid)

            self.is_normalised = True

        elif norm_method == "quantile_range":
            # Normalisation to [0, 1]-ish range based on quantiles.

            # Get voxel grid
            voxel_grid = self.get_voxel_grid()

            min_quantile = 0.0
            if not np.isnan(intensity_range[0]):
                min_quantile = intensity_range[0]

            max_quantile = 1.0
            if not np.isnan(intensity_range[1]):
                max_quantile = intensity_range[1]

            # Compute quantiles from voxel grid.
            min_int = np.quantile(voxel_grid[mask], q=min_quantile)
            max_int = np.quantile(voxel_grid[mask], q=max_quantile)

            # Normalise by range
            if not max_int == min_int:
                voxel_grid = (voxel_grid - min_int) / (max_int - min_int)
            else:
                voxel_grid = voxel_grid - min_int

            # Update the voxel grid
            self.set_voxel_grid(voxel_grid=voxel_grid)

            self.is_normalised = True

        elif norm_method == "standardisation":
            # Normalisation to mean 0 and standard deviation 1.

            # Get voxel grid
            voxel_grid = self.get_voxel_grid()

            # Determine mean and standard deviation of the voxel intensities
            mean_int = np.mean(voxel_grid[mask])
            sd_int = np.std(voxel_grid[mask])

            # Protect against invariance.
            if sd_int == 0.0:
                sd_int = 1.0

            # Normalise
            voxel_grid = (voxel_grid - mean_int) / sd_int

            # Update the voxel grid
            self.set_voxel_grid(voxel_grid=voxel_grid)

            self.is_normalised = True
        else:
            raise ValueError(
                f"{norm_method} is not a valid method for normalising intensity values."
            )

        self.saturate(intensity_range=saturation_range)

    def rotate(self, angle):
        """Rotate volume along z-axis."""

        # Skip for missing images
        if self.is_missing:
            return

        import scipy.ndimage as ndi

        from luna.radiology.mirp.featureSets.volumeMorphology import get_rotation_matrix

        # Find actual output size of x-y plane
        new_z_dim = np.asmatrix([self.size[0], 0.0, 0.0]) * get_rotation_matrix(
            np.radians(angle), dim=3, rot_axis=0
        )
        new_y_dim = np.asmatrix([0.0, self.size[1], 0.0]) * get_rotation_matrix(
            np.radians(angle), dim=3, rot_axis=0
        )
        new_x_dim = np.asmatrix([0.0, 0.0, self.size[2]]) * get_rotation_matrix(
            np.radians(angle), dim=3, rot_axis=0
        )
        new_dim_flt = np.squeeze(
            np.array(np.abs(new_z_dim))
            + np.array(np.abs(new_y_dim) + np.abs(new_x_dim))
        )

        # Get voxel grid
        voxel_grid = self.get_voxel_grid()

        # Rotate voxels along angle in the y-x plane and find truncated output size
        voxel_grid = ndi.rotate(
            voxel_grid.astype(np.float32),
            angle=angle,
            axes=(1, 2),
            reshape=True,
            order=1,
            mode="nearest",
        )
        new_dim_int = np.array(np.shape(voxel_grid)) * 1.0

        if (self.modality == "CT") and (self.spat_transform == "base"):
            voxel_grid = np.round(voxel_grid)

        # Update spacing
        self.spacing *= new_dim_int / new_dim_flt

        # Set rotation angle
        self.rotation_angle = angle

        # Update voxel grid with rotated voxels
        self.set_voxel_grid(voxel_grid=voxel_grid)

    def crop(
        self,
        ind_ext_z=None,
        ind_ext_y=None,
        ind_ext_x=None,
        xy_only=False,
        z_only=False,
    ):
        """ "Crop image to the provided map extent."""

        # Skip for missing images
        if self.is_missing:
            return

        # Determine corresponding voxel indices
        max_ind = np.ceil(
            np.array((np.max(ind_ext_z), np.max(ind_ext_y), np.max(ind_ext_x)))
        ).astype(np.int)
        min_ind = np.floor(
            np.array((np.min(ind_ext_z), np.min(ind_ext_y), np.min(ind_ext_x)))
        ).astype(np.int)

        # Set bounding indices
        max_bound_ind = np.minimum(max_ind, self.size).astype(np.int)
        min_bound_ind = np.maximum(min_ind, np.array([0, 0, 0])).astype(np.int)

        # Get voxel grid
        voxel_grid = self.get_voxel_grid()

        # Create corresponding image volumes by slicing original volume
        if z_only:
            voxel_grid = voxel_grid[min_bound_ind[0] : max_bound_ind[0] + 1, :, :]
            min_bound_ind[1] = 0
            min_bound_ind[2] = 0
        elif xy_only:
            voxel_grid = voxel_grid[
                :,
                min_bound_ind[1] : max_bound_ind[1] + 1,
                min_bound_ind[2] : max_bound_ind[2] + 1,
            ]
            min_bound_ind[0] = 0
            max_bound_ind[0] = self.size[0].astype(np.int)
        else:
            voxel_grid = voxel_grid[
                min_bound_ind[0] : max_bound_ind[0] + 1,
                min_bound_ind[1] : max_bound_ind[1] + 1,
                min_bound_ind[2] : max_bound_ind[2] + 1,
            ]

        # Update origin and z-slice position
        self.origin = self.origin + np.dot(self.m_affine, np.transpose(min_bound_ind))

        # Update voxel grid
        self.set_voxel_grid(voxel_grid=voxel_grid)

    def crop_to_size(self, center, crop_size, xy_only=False):
        """Crop images to the exact size"""

        # Skip for missing images
        if self.is_missing:
            return

        # Make local copy
        crop_size = np.array(copy.deepcopy(crop_size))

        # Determine the new grid origin in the original index space. Only the dimensions with a number are updated
        grid_origin = np.round(center - crop_size / 2.0).astype(np.int)

        # Update grid origin and crop_size for the remainder of the calculation
        grid_origin[np.isnan(crop_size)] = 0
        crop_size[np.isnan(crop_size)] = self.size[np.isnan(crop_size)]

        # Determine coordinates of the box that can be copied in the original space
        max_ind_orig = grid_origin + crop_size
        min_ind_orig = grid_origin

        # Update coordinates based on boundaries in the original images
        max_ind_orig = np.minimum(max_ind_orig, self.size).astype(np.int)
        min_ind_orig = np.maximum(min_ind_orig, [0, 0, 0]).astype(np.int)

        # Determine coordinates where this box should land, i.e. perform the coordinate transformation to grid index space.
        max_ind_grid = max_ind_orig - grid_origin
        min_ind_grid = min_ind_orig - grid_origin

        # Create an empty voxel_grid to copy to
        cropped_grid = np.full(crop_size.astype(np.int), fill_value=np.nan)

        # Get slice of voxel grid
        voxel_grid = self.get_voxel_grid()[
            min_ind_orig[0] : max_ind_orig[0],
            min_ind_orig[1] : max_ind_orig[1],
            min_ind_orig[2] : max_ind_orig[2],
        ]

        # Put the voxel grid slice into the cropped grid
        cropped_grid[
            min_ind_grid[0] : max_ind_grid[0],
            min_ind_grid[1] : max_ind_grid[1],
            min_ind_grid[2] : max_ind_grid[2],
        ] = voxel_grid

        # Replace any remaining NaN values in the grid by the lowest intensity in voxel_grid
        cropped_grid[np.isnan(cropped_grid)] = np.min(voxel_grid)

        # Restore the original dtype in case it got lost
        cropped_grid = cropped_grid.astype(voxel_grid.dtype)

        # Update origin
        self.origin = self.origin + np.dot(self.m_affine, np.transpose(grid_origin))

        # Set voxel grid
        self.set_voxel_grid(voxel_grid=cropped_grid)

    def set_spatial_transform(self, transform_method: str):
        if transform_method == "base":
            self.spat_transform = "base"

        else:
            self.spat_transform = transform_method
            self.as_parametric_map = True

    def compute_diagnostic_features(self, append_str: str = ""):
        """Creates diagnostic features for the image stack"""

        # Set feature names
        feat_names = [
            "img_dim_x",
            "img_dim_y",
            "img_dim_z",
            "vox_dim_x",
            "vox_dim_y",
            "vox_dim_z",
            "mean_int",
            "min_int",
            "max_int",
        ]

        # Generate an initial table
        feature_table = pd.DataFrame(
            np.full(shape=(1, len(feat_names)), fill_value=np.nan)
        )
        feature_table.columns = feat_names

        if not self.is_missing:
            # Update columns with actual values
            feature_table["img_dim_x"] = self.size[2]
            feature_table["img_dim_y"] = self.size[1]
            feature_table["img_dim_z"] = self.size[0]
            feature_table["vox_dim_x"] = self.spacing[2]
            feature_table["vox_dim_y"] = self.spacing[1]
            feature_table["vox_dim_z"] = self.spacing[0]
            feature_table["mean_int"] = np.mean(self.get_voxel_grid())
            feature_table["min_int"] = np.min(self.get_voxel_grid())
            feature_table["max_int"] = np.max(self.get_voxel_grid())

            # Update column names
        feature_table.columns = [
            "_".join(["diag", feature, append_str]).rstrip("_")
            for feature in feature_table.columns
        ]

        return feature_table

    def export(self, file_path):
        """
        Exports the image to the requested directory
        :param file_path: directory to write image to
        :return:
        """

        # Skip if the image is missing
        if self.is_missing:
            return

        # Construct file name
        file_name = self.get_export_descriptor() + ".nii"

        # Export image to file
        return self.write(file_path=file_path, file_name=file_name)

    def get_export_descriptor(self):
        """
        Generates an image descriptor based on parameters of the image
        :return:
        """

        descr_list = ["volumes"]

        # Image name and spatial transformation
        if self.name is not None:
            descr_list += [self.name]
        if self.modality is not None:
            descr_list += [self.modality]

        if self.interpolated:
            # Interpolation
            descr_list += [
                self.interpolation_algorithm,
                "x",
                str(self.spacing[2])[:5],
                "y",
                str(self.spacing[1])[:5],
                "z",
                str(self.spacing[0])[:5],
            ]
        if self.binned:
            descr_list += ["bin", str(self.bin_width)]

        if self.rotation_angle != 0.0:
            # Rotation angle
            descr_list += ["rot", str(self.rotation_angle)[:5]]

        if not (
            self.transl_fraction_x == 0.0
            and self.transl_fraction_y == 0.0
            and self.transl_fraction_z == 0.0
        ):
            # Translation fraction
            descr_list += [
                "trans",
                "x",
                str(self.transl_fraction_x)[:5],
                "y",
                str(self.transl_fraction_y)[:5],
                "z",
                str(self.transl_fraction_z)[:5],
            ]

        if self.noise != -1.0:
            # Noise level
            descr_list += ["noise", str(self.noise)[:5], "iter", str(self.noise_iter)]

        if not self.spat_transform == "base":
            descr_list += [self.spat_transform]

        return "_".join(descr_list)

    def convert2sitk(self):
        """Converts image object back to Simple ITK
        This step may precede writing to file."""

        import SimpleITK as sitk

        sitk.ProcessObject.SetGlobalDefaultNumberOfThreads(4)

        # Skip if the image is missing
        if self.is_missing:
            return None

        # Get image data type and set a valid data type that can be read by simple itk
        vox_dtype = self.dtype_name
        if vox_dtype in [
            "float16",
            "float32",
            "float64",
            "float80",
            "float96",
            "float128",
        ]:
            cast_type = np.float32
        elif vox_dtype in ["bool", "uint8"]:
            cast_type = np.uint8
        else:
            raise RuntimeError(f"Invalid ITK value {vox_dtype}")

        # Convert image voxels
        sitk_img = sitk.GetImageFromArray(
            self.get_voxel_grid().astype(cast_type), isVector=False
        )
        sitk_img.SetOrigin(self.origin[::-1])
        sitk_img.SetSpacing(self.spacing[::-1])
        sitk_img.SetDirection(self.orientation[::-1])

        return sitk_img

    def write(self, file_path, file_name):
        """Writes the image to a file"""
        import os

        import SimpleITK as sitk

        sitk.ProcessObject.SetGlobalDefaultNumberOfThreads(4)

        # Check if path exists
        file_path = os.path.normpath(file_path)
        if not os.path.exists(file_path):
            os.makedirs(file_path)

        # Add file and file name
        file_path = os.path.join(file_path, file_name)

        # Convert image to simple itk format
        sitk_img = self.convert2sitk()

        # Write to file using simple itk, and the image is not missing
        if sitk_img is not None:
            print(f"Writing to {file_path}")
            sitk.WriteImage(sitk_img, file_path, True)

        return file_path

    def write_dicom(self, file_path, file_name, bit_depth=16):
        if self.metadata is None:
            raise ValueError(
                "Image slice cannot be written to DICOM as metadata is missing."
            )

        # Set pixeldata
        self.set_pixel_data(bit_depth=bit_depth)

        # Remove 0x0002 group data.
        for filemeta_tag in list(self.metadata.group_dataset(0x0002).keys()):
            del self.metadata[filemeta_tag]

        # Save to file
        self.metadata.save_as(
            filename=os.path.join(file_path, file_name), write_like_original=False
        )

    def write_dicom_series(self, file_path, file_name="", bit_depth=16):
        if self.metadata is None:
            raise ValueError(
                "Image cannot be written as a DICOM series as metadata is missing."
            )

        # Check if the write folder exists
        if not os.path.isdir(file_path):
            if os.path.isfile(file_path):
                # Check if the write folder is a file.
                raise IOError(
                    f"{file_path} is an existing file, not a directory. No DICOM images were exported."
                )
            else:
                os.makedirs(file_path, exist_ok=True)

        if self.as_parametric_map:
            self._convert_to_parametric_map_iod()

        # Provide a new series UID
        self.set_metadata(tag=(0x0020, 0x000E), value=create_new_uid(dcm=self.metadata))

        if self.modality == "PT":
            # Update the number of slices attribute.
            self.set_metadata(tag=(0x0054, 0x0081), value=self.size[0])

        # Export per slice
        for ii in np.arange(self.size[0]):
            # Obtain the slice
            slice_obj: ImageClass = self.get_slices(slice_number=ii)[0]

            # Generate the file name
            slice_file_name = file_name + f"{ii:06d}" + ".dcm"

            # Set instance number
            slice_obj.set_metadata(tag=(0x0020, 0x0013), value=ii)

            # Update the SOP instance UID to avoid collisions
            slice_obj.set_metadata(
                tag=(0x0008, 0x0018), value=create_new_uid(dcm=slice_obj.metadata)
            )

            # Update instance creation date and time
            slice_obj.set_metadata(tag=(0x0008, 0x0012), value=time.strftime("%Y%m%d"))
            slice_obj.set_metadata(tag=(0x0008, 0x0013), value=time.strftime("%H%M%S"))

            # Export
            slice_obj.write_dicom(
                file_path=file_path, file_name=slice_file_name, bit_depth=bit_depth
            )

    def get_slices(self, slice_number=None):
        img_obj_list = []

        # Create a copy of the current object
        base_img_obj = self.copy()

        # Remove attributes that need to be set
        base_img_obj.isEncoded_voxel_grid = None
        base_img_obj.voxel_grid = None
        base_img_obj.size = None
        base_img_obj.dtype_name = None

        if slice_number is None:
            voxel_grid = self.get_voxel_grid()

            # Iterate over slices
            for ii in np.arange(self.size[0]):
                slice_img_obj = copy.deepcopy(base_img_obj)

                # Update origin and slice position
                slice_img_obj.origin += np.dot(self.m_affine, np.array([ii, 0, 0]))

                # Update name
                if slice_img_obj.name is not None:
                    slice_img_obj.name += "_slice_" + str(ii)

                # Copy voxel grid without dropping dimensions.
                if voxel_grid is not None:
                    slice_img_obj.set_voxel_grid(
                        voxel_grid=voxel_grid[ii : ii + 1, :, :]
                    )

                # Add to list
                img_obj_list += [slice_img_obj]
        else:
            # Extract a single slice
            slice_img_obj = copy.deepcopy(base_img_obj)

            # Update origin and slice position
            slice_img_obj.origin += np.dot(
                self.m_affine, np.array([slice_number, 0, 0])
            )

            # Update name
            if slice_img_obj.name is not None:
                slice_img_obj.name += "_slice_" + str(slice_number)

            # Copy voxel grid without dropping dimensions.
            if not self.is_missing:
                slice_img_obj.set_voxel_grid(
                    voxel_grid=self.get_voxel_grid()[
                        slice_number : slice_number + 1, :, :
                    ]
                )

            # Add to list
            img_obj_list += [slice_img_obj]

        return img_obj_list

    def drop_image(self):
        """Drops image, e.g. to free up memory. We don't set the is_m"""
        self.isEncoded_voxel_grid = None
        self.voxel_grid = None

    def get_metadata(self, tag, tag_type, default=None):
        # Do not attempt to read the metadata if no metadata is present.
        if self.metadata is None:
            return

        return get_pydicom_meta_tag(
            dcm_seq=self.metadata, tag=tag, tag_type=tag_type, default=default
        )

    def set_metadata(self, tag, value, force_vr=None):
        # Do not update the metadata if no metadata is present.
        if self.metadata is None:
            return None

        set_pydicom_meta_tag(
            dcm_seq=self.metadata, tag=tag, value=value, force_vr=force_vr
        )

    def has_metadata(self, tag):
        if self.metadata is None:
            return None

        else:
            return get_pydicom_meta_tag(dcm_seq=self.metadata, tag=tag, test_tag=True)

    def delete_metadata(self, tag):
        if self.metadata is None:
            return None

        elif self.has_metadata(tag):
            del self.metadata[tag]

        else:
            pass

    def cond_set_metadata(self, tag, value, force_vr=None):
        if self.set_metadata is None:
            return None

        elif not self.has_metadata(tag):
            self.set_metadata(tag=tag, value=value, force_vr=force_vr)

        else:
            pass

    def update_image_plane_metadata(self):
        # Update pixel spacing, image orientation, image position, slice thickness, slice location, rows and columns based on the image object

        # Do not update the metadata if no metadata is present.
        if self.metadata is None:
            return

        # Pixel spacing
        pixel_spacing = [self.spacing[2], self.spacing[1]]
        self.set_metadata(tag=(0x0028, 0x0030), value=pixel_spacing)

        # Image orientation
        image_orientation = self.orientation[::-1][:6]
        self.set_metadata(tag=(0x0020, 0x0037), value=image_orientation)

        # Image position
        self.set_metadata(tag=(0x0020, 0x0032), value=self.origin[::-1])

        # Slice thickness
        self.set_metadata(tag=(0x0018, 0x0050), value=self.spacing[0])

        # Slice location
        self.set_metadata(tag=(0x0020, 0x1041), value=self.origin[0])

        # Rows (y)
        self.set_metadata(tag=(0x0028, 0x0010), value=self.size[1])

        # Columns (x)
        self.set_metadata(tag=(0x0028, 0x0011), value=self.size[2])

    def set_pixel_data(self, bit_depth=16):
        # Important tags to update:

        if self.metadata is None:
            return

        if self.size[0] > 1:
            warnings.warn(
                "Cannot set pixel data for image with more than one slice.", UserWarning
            )
            return

        # Set samples per pixel
        self.set_metadata(tag=(0x0028, 0x0002), value=1)

        # Set photometric interpretation
        self.set_metadata(tag=(0x0028, 0x0004), value="MONOCHROME2")

        # Remove the Pixel Data Provider URL attribute
        self.delete_metadata(tag=(0x0028, 0x7FE0))

        # Determine how pixel data are stored.
        if self.as_parametric_map:
            self._set_pixel_data_float(bit_depth=16)

        else:
            self._set_pixel_data_int(bit_depth=16)

    def _set_pixel_data_int(self, bit_depth):
        # Set dtype for the image
        if bit_depth == 8:
            pixel_type = np.int8
        elif bit_depth == 16:
            pixel_type = np.int16
        elif bit_depth == 32:
            pixel_type = np.int32
        elif bit_depth == 64:
            pixel_type = np.int64
        else:
            raise ValueError(
                f"Bit depth of DICOM images should be one of 8, 16 (default), 32, or 64., Found: {bit_depth}"
            )

        # Update metadata related to the image data
        self.update_image_plane_metadata()

        pixel_grid = np.squeeze(self.get_voxel_grid(), axis=0)

        # Always write 16-bit data
        self.set_metadata(tag=(0x0028, 0x0100), value=bit_depth)  # Bits allocated
        self.set_metadata(tag=(0x0028, 0x0101), value=bit_depth)  # Bits stored
        self.set_metadata(tag=(0x0028, 0x0102), value=bit_depth - 1)  # High-bit
        self.set_metadata(
            tag=(0x0028, 0x0103), value=1
        )  # Pixel representation (we assume signed integers)

        # Standard settings for lowest and highest pixel value
        if self.modality == "CT":
            rescale_intercept = 0.0
            rescale_slope = 1.0

        elif self.modality == "PT":
            rescale_intercept = 0.0
            rescale_slope = np.max(pixel_grid) / (2 ** (bit_depth - 1) - 1)

        elif self.modality == "MR":
            rescale_intercept = 0.0
            rescale_slope = 1.0

        else:
            raise TypeError(f"Unknown modality {self.modality}")

        # Inverse rescaling prior to dicom storage
        pixel_grid = (pixel_grid - rescale_intercept) / rescale_slope

        # Convert back to int16
        pixel_grid = pixel_grid.astype(pixel_type)

        # Store to PixelData
        self.set_metadata(
            tag=(0x7FE0, 0x0010), value=pixel_grid.tobytes(), force_vr="OW"
        )

        # Delete other PixelData containers
        self.delete_metadata(tag=(0x7FE0, 0x0008))  # Float pixel data
        self.delete_metadata(tag=(0x7FE0, 0x0009))  # Double float pixel data

        # Update rescale intercept and slope (in case of CT and PET only)
        if self.modality in ["CT", "PT"]:
            self.set_metadata(tag=(0x0028, 0x1052), value=rescale_intercept)
            self.set_metadata(tag=(0x0028, 0x1053), value=rescale_slope)

        # Remove elements of the VOI LUT module
        self.delete_metadata(tag=(0x0028, 0x3010))  # VOI LUT sequence
        self.delete_metadata(tag=(0x0028, 0x1050))  # Window center
        self.delete_metadata(tag=(0x0028, 0x1051))  # Window width
        self.delete_metadata(
            tag=(0x0028, 0x1055)
        )  # Window center and width explanation
        self.delete_metadata(tag=(0x0028, 0x1056))  # VOI LUT function

        # Update smallest and largest image pixel value. Cannot set more than 16 bits due to limitations of the
        # tag.
        if bit_depth <= 16:
            self.set_metadata(
                tag=(0x0028, 0x0106), value=np.min(pixel_grid), force_vr="SS"
            )
            self.set_metadata(
                tag=(0x0028, 0x0107), value=np.max(pixel_grid), force_vr="SS"
            )

    def _set_pixel_data_float(self, bit_depth):
        if not self.as_parametric_map:
            raise ValueError(
                "Floating point representation in DICOM is only supported by parametric maps, but the image in MIRP is not marked for conversion of the metadata to"
                "parametric maps."
            )

        # Set dtype for the image
        if bit_depth == 16:
            pixel_type = np.int16
        elif bit_depth == 32:
            pixel_type = np.float32
        elif bit_depth == 64:
            pixel_type = np.float64
        else:
            raise ValueError(
                f"Bit depth of floating point DICOM images should be 16, 32 (default), or 64. Found: {bit_depth}"
            )

        # Set the number of allocated bits
        self.set_metadata(tag=(0x0028, 0x0100), value=bit_depth)  # Bits allocated

        # Update metadata related to the image data
        self.update_image_plane_metadata()

        # Get the pixel grid of the slice
        pixel_grid = np.squeeze(self.get_voxel_grid(), axis=0)

        # Define rescale intercept and slope
        if bit_depth == 16:
            # DICOM ranges
            dcm_range = float(2**bit_depth - 1)
            dcm_min = -float(2 ** (bit_depth - 1))
            dcm_max = float(2 ** (bit_depth - 1)) - 1.0

            # Data ranges
            value_min = np.min(pixel_grid)
            value_max = np.max(pixel_grid)
            value_range = np.max(pixel_grid) - np.min(pixel_grid)

            # Rescale slope and intercept
            rescale_slope = value_range / dcm_range
            rescale_intercept = (value_min * dcm_max - value_max * dcm_min) / dcm_range

            # Inverse rescaling prior to dicom storage
            pixel_grid = (
                pixel_grid * (dcm_range) - value_min * dcm_max + value_max * dcm_min
            )
            pixel_grid *= 1.0 / value_range

            # Round pixel values for safety
            pixel_grid = np.round(pixel_grid)
        else:
            rescale_intercept = 0.0
            rescale_slope = 1.0

        # Cast to the right data type
        pixel_grid = pixel_grid.astype(pixel_type)

        # Store pixel data to the right attribute
        if bit_depth == 16:
            # Store to Pixel Data attribute
            self.set_metadata(
                tag=(0x7FE0, 0x0010), value=pixel_grid.tobytes(), force_vr="OW"
            )

            # Set Image Pixel module-specific tags
            self.set_metadata(tag=(0x0028, 0x0101), value=bit_depth)  # Bits stored
            self.set_metadata(tag=(0x0028, 0x0102), value=bit_depth - 1)  # High-bit
            self.set_metadata(
                tag=(0x0028, 0x0103), value=1
            )  # Pixel representation (we assume signed integers)

            # Update smallest and largest pixel value
            self.set_metadata(
                tag=(0x0028, 0x0106), value=np.min(pixel_grid), force_vr="SS"
            )
            self.set_metadata(
                tag=(0x0028, 0x0107), value=np.max(pixel_grid), force_vr="SS"
            )

            # Delete other PixelData containers
            self.delete_metadata(tag=(0x7FE0, 0x0008))  # Float pixel data
            self.delete_metadata(tag=(0x7FE0, 0x0009))  # Double float pixel data

        elif bit_depth == 32:
            # Store to Float Pixel Data attribute
            self.set_metadata(
                tag=(0x7FE0, 0x0008), value=np.ravel(pixel_grid).tolist(), force_vr="OF"
            )

            # Delete other PixelData containers
            self.delete_metadata(tag=(0x7FE0, 0x0009))  # Double float pixel data
            self.delete_metadata(tag=(0x7FE0, 0x0010))  # Integer type data

        elif bit_depth == 64:
            # Store to Double Float Pixel Data attribute
            self.set_metadata(
                tag=(0x7FE0, 0x0009), value=np.ravel(pixel_grid).tolist(), force_vr="OD"
            )

            # Delete other PixelData containers
            self.delete_metadata(tag=(0x7FE0, 0x0008))  # Float pixel data
            self.delete_metadata(tag=(0x7FE0, 0x0010))  # Integer type data

        # Reset rescale intercept and slope attributes to default values for parametric maps.
        self.set_metadata(
            tag=(0x0028, 0x1052), value=rescale_intercept
        )  # Rescale intercept
        self.set_metadata(tag=(0x0028, 0x1053), value=rescale_slope)  # Rescale slope

        # Remove elements of the VOI LUT module
        self.delete_metadata(tag=(0x0028, 0x3010))  # VOI LUT sequence
        self.delete_metadata(tag=(0x0028, 0x1050))  # Window center
        self.delete_metadata(tag=(0x0028, 0x1051))  # Window width
        self.delete_metadata(
            tag=(0x0028, 0x1055)
        )  # Window center and width explanation
        self.delete_metadata(tag=(0x0028, 0x1056))  # VOI LUT function

        # Number of frames
        self.set_metadata(tag=(0x0028, 0x0008), value=1)

    def _convert_to_parametric_map_iod(self):
        if self.metadata is None:
            return None

        # Create a copy of the metadata.
        old_dcm: FileDataset = copy.deepcopy(self.metadata)

        # Update the SOP class to that of a parametric map image
        self.set_metadata(tag=(0x0008, 0x0016), value="1.2.840.10008.5.1.4.1.1.30")

        # Update the image type attribute
        image_type = self.get_metadata(
            tag=(0x0008, 0x0008), tag_type="mult_str", default=[]
        )
        image_type = [image_type[ii] if ii < len(image_type) else "" for ii in range(4)]
        image_type[0] = "DERIVED"
        image_type[1] = "PRIMARY"
        image_type[2] = image_type[2] if not image_type[2] == "" else "STATIC"
        image_type[3] = "MIXED" if self.spat_transform == "base" else "FILTERED"

        self.set_metadata(tag=(0x0008, 0x0008), value=image_type)

        # Parametric Map Image module attributes that may be missing.
        self.cond_set_metadata(
            tag=(0x2050, 0x0020), value="IDENTITY"
        )  # Presentation LUT shape
        self.cond_set_metadata(
            tag=(0x0018, 0x9004), value="RESEARCH"
        )  # Content qualification
        self.cond_set_metadata(tag=(0x0028, 0x0301), value="NO")  # Burned-in Annotation
        self.cond_set_metadata(
            tag=(0x0028, 0x0302), value="YES"
        )  # Recognisable facial features
        self.cond_set_metadata(
            tag=(0x0070, 0x0080),
            value=self.get_export_descriptor().upper().strip()[:15],
        )  # Content label
        self.cond_set_metadata(
            tag=(0x0070, 0x0081), value=self.get_export_descriptor()[:63]
        )  # Content description
        self.cond_set_metadata(tag=(0x0070, 0x0084), value="Doe^John")

        # Set the source instance sequence
        source_instance_list = []
        for reference_instance_sop_uid in self.slice_table.sop_instance_uid:
            ref_inst = Dataset()
            set_pydicom_meta_tag(
                dcm_seq=ref_inst,
                tag=(0x0008, 0x1150),
                value=get_pydicom_meta_tag(
                    dcm_seq=old_dcm, tag=(0x0008, 0x0016), tag_type="str"
                ),
            )
            set_pydicom_meta_tag(
                dcm_seq=ref_inst, tag=(0x0008, 0x1155), value=reference_instance_sop_uid
            )

            source_instance_list += [ref_inst]

        self.set_metadata(tag=(0x0008, 0x2112), value=Sequence(source_instance_list))

        # Attributes from the enhanced general equipment module may be missing.
        self.cond_set_metadata(tag=(0x0008, 0x0070), value="unknown")  # Manufacturer
        self.cond_set_metadata(tag=(0x0008, 0x1090), value="unknown")  # Model name
        self.cond_set_metadata(
            tag=(0x0018, 0x1000), value="unknown"
        )  # Device Serial Number
        self.set_metadata(tag=(0x0018, 0x1020), value="MIRP " + get_version())

        # Items from multi-frame function groups may be missing. We currently only use a single frame.
        self.set_metadata(
            tag=(0x5200, 0x9229), value=Sequence()
        )  # Shared functional groups sequence
        self.set_metadata(
            tag=(0x5200, 0x9230), value=Sequence()
        )  # Per-frame functional groups sequence

        # Multi-frame Dimension module

        # Dimension organisation sequence. We copy the frame of reference as UID.
        dim_org_seq_elem = Dataset()
        set_pydicom_meta_tag(
            dim_org_seq_elem,
            tag=(0x0020, 0x9164),
            value=self.get_metadata(tag=(0x0020, 0x0052), tag_type="str"),
        )  # Dimension organisation UID
        self.set_metadata(tag=(0x0020, 0x9221), value=Sequence([dim_org_seq_elem]))

        # Dimension Index sequence. We point to the instance number.
        dim_index_seq_elem = Dataset()
        set_pydicom_meta_tag(
            dim_index_seq_elem, tag=(0x0020, 0x9165), value=(0x0020, 0x0013)
        )  # Dimension index pointer
        set_pydicom_meta_tag(
            dim_index_seq_elem,
            tag=(0x0020, 0x9164),
            value=self.get_metadata(tag=(0x00200052), tag_type="str"),
        )  # Dimension organisation UID
        self.set_metadata(tag=(0x0020, 0x9222), value=Sequence([dim_index_seq_elem]))

add_noise(noise_level, noise_iter)

Adds Gaussian noise to the image volume noise_level: standard deviation of image noise present

Source code in src/luna/radiology/mirp/imageClass.py

def add_noise(self, noise_level, noise_iter):
    """Adds Gaussian noise to the image volume
    noise_level: standard deviation of image noise present"""

    # Add noise iteration number
    self.noise_iter = noise_iter

    # Skip for missing images
    if self.is_missing:
        return

    # Skip for invalid noise levels
    if noise_level is None:
        return
    if np.isnan(noise_level) or noise_level < 0.0:
        return

    # Add Gaussian noise to image
    voxel_grid = self.get_voxel_grid()
    voxel_grid += np.random.normal(loc=0.0, scale=noise_level, size=self.size)

    # Check for corrections due to image modality
    if self.spat_transform == "base":
        # Round CT values to the nearest integer
        if self.modality == "CT":
            voxel_grid = np.round(a=voxel_grid, decimals=0)

        # Set minimum PT to 0.0
        if self.modality == "PT":
            voxel_grid[voxel_grid < 0.0] = 0.0

    # Set noise level in image
    self.noise = noise_level

    self.set_voxel_grid(voxel_grid=voxel_grid)

compute_diagnostic_features(append_str='')

Creates diagnostic features for the image stack

Source code in src/luna/radiology/mirp/imageClass.py

def compute_diagnostic_features(self, append_str: str = ""):
    """Creates diagnostic features for the image stack"""

    # Set feature names
    feat_names = [
        "img_dim_x",
        "img_dim_y",
        "img_dim_z",
        "vox_dim_x",
        "vox_dim_y",
        "vox_dim_z",
        "mean_int",
        "min_int",
        "max_int",
    ]

    # Generate an initial table
    feature_table = pd.DataFrame(
        np.full(shape=(1, len(feat_names)), fill_value=np.nan)
    )
    feature_table.columns = feat_names

    if not self.is_missing:
        # Update columns with actual values
        feature_table["img_dim_x"] = self.size[2]
        feature_table["img_dim_y"] = self.size[1]
        feature_table["img_dim_z"] = self.size[0]
        feature_table["vox_dim_x"] = self.spacing[2]
        feature_table["vox_dim_y"] = self.spacing[1]
        feature_table["vox_dim_z"] = self.spacing[0]
        feature_table["mean_int"] = np.mean(self.get_voxel_grid())
        feature_table["min_int"] = np.min(self.get_voxel_grid())
        feature_table["max_int"] = np.max(self.get_voxel_grid())

        # Update column names
    feature_table.columns = [
        "_".join(["diag", feature, append_str]).rstrip("_")
        for feature in feature_table.columns
    ]

    return feature_table

convert2sitk()

Converts image object back to Simple ITK This step may precede writing to file.

Source code in src/luna/radiology/mirp/imageClass.py

def convert2sitk(self):
    """Converts image object back to Simple ITK
    This step may precede writing to file."""

    import SimpleITK as sitk

    sitk.ProcessObject.SetGlobalDefaultNumberOfThreads(4)

    # Skip if the image is missing
    if self.is_missing:
        return None

    # Get image data type and set a valid data type that can be read by simple itk
    vox_dtype = self.dtype_name
    if vox_dtype in [
        "float16",
        "float32",
        "float64",
        "float80",
        "float96",
        "float128",
    ]:
        cast_type = np.float32
    elif vox_dtype in ["bool", "uint8"]:
        cast_type = np.uint8
    else:
        raise RuntimeError(f"Invalid ITK value {vox_dtype}")

    # Convert image voxels
    sitk_img = sitk.GetImageFromArray(
        self.get_voxel_grid().astype(cast_type), isVector=False
    )
    sitk_img.SetOrigin(self.origin[::-1])
    sitk_img.SetSpacing(self.spacing[::-1])
    sitk_img.SetDirection(self.orientation[::-1])

    return sitk_img

crop(ind_ext_z=None, ind_ext_y=None, ind_ext_x=None, xy_only=False, z_only=False)

"Crop image to the provided map extent.

Source code in src/luna/radiology/mirp/imageClass.py

def crop(
    self,
    ind_ext_z=None,
    ind_ext_y=None,
    ind_ext_x=None,
    xy_only=False,
    z_only=False,
):
    """ "Crop image to the provided map extent."""

    # Skip for missing images
    if self.is_missing:
        return

    # Determine corresponding voxel indices
    max_ind = np.ceil(
        np.array((np.max(ind_ext_z), np.max(ind_ext_y), np.max(ind_ext_x)))
    ).astype(np.int)
    min_ind = np.floor(
        np.array((np.min(ind_ext_z), np.min(ind_ext_y), np.min(ind_ext_x)))
    ).astype(np.int)

    # Set bounding indices
    max_bound_ind = np.minimum(max_ind, self.size).astype(np.int)
    min_bound_ind = np.maximum(min_ind, np.array([0, 0, 0])).astype(np.int)

    # Get voxel grid
    voxel_grid = self.get_voxel_grid()

    # Create corresponding image volumes by slicing original volume
    if z_only:
        voxel_grid = voxel_grid[min_bound_ind[0] : max_bound_ind[0] + 1, :, :]
        min_bound_ind[1] = 0
        min_bound_ind[2] = 0
    elif xy_only:
        voxel_grid = voxel_grid[
            :,
            min_bound_ind[1] : max_bound_ind[1] + 1,
            min_bound_ind[2] : max_bound_ind[2] + 1,
        ]
        min_bound_ind[0] = 0
        max_bound_ind[0] = self.size[0].astype(np.int)
    else:
        voxel_grid = voxel_grid[
            min_bound_ind[0] : max_bound_ind[0] + 1,
            min_bound_ind[1] : max_bound_ind[1] + 1,
            min_bound_ind[2] : max_bound_ind[2] + 1,
        ]

    # Update origin and z-slice position
    self.origin = self.origin + np.dot(self.m_affine, np.transpose(min_bound_ind))

    # Update voxel grid
    self.set_voxel_grid(voxel_grid=voxel_grid)

crop_to_size(center, crop_size, xy_only=False)

Crop images to the exact size

Source code in src/luna/radiology/mirp/imageClass.py

def crop_to_size(self, center, crop_size, xy_only=False):
    """Crop images to the exact size"""

    # Skip for missing images
    if self.is_missing:
        return

    # Make local copy
    crop_size = np.array(copy.deepcopy(crop_size))

    # Determine the new grid origin in the original index space. Only the dimensions with a number are updated
    grid_origin = np.round(center - crop_size / 2.0).astype(np.int)

    # Update grid origin and crop_size for the remainder of the calculation
    grid_origin[np.isnan(crop_size)] = 0
    crop_size[np.isnan(crop_size)] = self.size[np.isnan(crop_size)]

    # Determine coordinates of the box that can be copied in the original space
    max_ind_orig = grid_origin + crop_size
    min_ind_orig = grid_origin

    # Update coordinates based on boundaries in the original images
    max_ind_orig = np.minimum(max_ind_orig, self.size).astype(np.int)
    min_ind_orig = np.maximum(min_ind_orig, [0, 0, 0]).astype(np.int)

    # Determine coordinates where this box should land, i.e. perform the coordinate transformation to grid index space.
    max_ind_grid = max_ind_orig - grid_origin
    min_ind_grid = min_ind_orig - grid_origin

    # Create an empty voxel_grid to copy to
    cropped_grid = np.full(crop_size.astype(np.int), fill_value=np.nan)

    # Get slice of voxel grid
    voxel_grid = self.get_voxel_grid()[
        min_ind_orig[0] : max_ind_orig[0],
        min_ind_orig[1] : max_ind_orig[1],
        min_ind_orig[2] : max_ind_orig[2],
    ]

    # Put the voxel grid slice into the cropped grid
    cropped_grid[
        min_ind_grid[0] : max_ind_grid[0],
        min_ind_grid[1] : max_ind_grid[1],
        min_ind_grid[2] : max_ind_grid[2],
    ] = voxel_grid

    # Replace any remaining NaN values in the grid by the lowest intensity in voxel_grid
    cropped_grid[np.isnan(cropped_grid)] = np.min(voxel_grid)

    # Restore the original dtype in case it got lost
    cropped_grid = cropped_grid.astype(voxel_grid.dtype)

    # Update origin
    self.origin = self.origin + np.dot(self.m_affine, np.transpose(grid_origin))

    # Set voxel grid
    self.set_voxel_grid(voxel_grid=cropped_grid)

decimate(by_slice)

Decimates image voxel grid by removing every second element :param by_slice: :return:

Source code in src/luna/radiology/mirp/imageClass.py

def decimate(self, by_slice):
    """
    Decimates image voxel grid by removing every second element
    :param by_slice:
    :return:
    """

    # Skip for missing images
    if self.is_missing:
        return

    # Get the voxel grid
    img_voxel_grid = self.get_voxel_grid()

    # Update the voxel grid
    if by_slice:
        # Drop every second pixel
        img_voxel_grid = img_voxel_grid[
            :, slice(None, None, 2), slice(None, None, 2)
        ]

        # Update voxel spacing
        self.spacing[[1, 2]] *= 2.0
    else:
        # Drop every second voxel
        img_voxel_grid = img_voxel_grid[
            slice(None, None, 2), slice(None, None, 2), slice(None, None, 2)
        ]

        # Update voxel spacing
        self.spacing *= 2.0

    # Update voxel grid. This also updates the size attribute.
    self.set_voxel_grid(voxel_grid=img_voxel_grid)

decode_voxel_grid()

Performs run length decoding of the voxel grid and converts it to a numpy array

Source code in src/luna/radiology/mirp/imageClass.py

def decode_voxel_grid(self):
    """Performs run length decoding of the voxel grid and converts it to a numpy array"""
    if self.dtype_name == "bool" and self.isEncoded_voxel_grid:
        decoded_voxel = np.zeros(np.prod(self.size), dtype=np.bool)

        # Check if the voxel grid contains values
        if self.voxel_grid is not None:
            decode_zip = copy.deepcopy(self.voxel_grid)
            for ii, jj in decode_zip:
                decoded_voxel[ii : jj + 1] = True

        # Set shape to original grid
        decoded_voxel = decoded_voxel.reshape(self.size)

        # Update self.voxel_grid and isEncoded_voxel_grid tags
        self.voxel_grid = decoded_voxel
        self.isEncoded_voxel_grid = False

drop_image()

Drops image, e.g. to free up memory. We don't set the is_m

Source code in src/luna/radiology/mirp/imageClass.py

def drop_image(self):
    """Drops image, e.g. to free up memory. We don't set the is_m"""
    self.isEncoded_voxel_grid = None
    self.voxel_grid = None

encode_voxel_grid(voxel_grid)

Performs run length encoding of the voxel grid

Source code in src/luna/radiology/mirp/imageClass.py

def encode_voxel_grid(self, voxel_grid):
    """Performs run length encoding of the voxel grid"""

    # Determine whether the voxel grid should be encoded (only True for boolean data types; typically roi)
    if self.dtype_name == "bool":
        # Run length encoding for "True"
        rle_end = np.array(
            np.append(
                np.where(voxel_grid.ravel()[1:] != voxel_grid.ravel()[:-1]),
                np.prod(self.size) - 1,
            )
        )
        rle_start = np.cumsum(np.append(0, np.diff(np.append(-1, rle_end))))[:-1]
        rle_val = voxel_grid.ravel()[rle_start]

        # Check whether the voxel grid is empty (consists of 0s)
        if np.all(~rle_val):
            self.voxel_grid = None
            self.isEncoded_voxel_grid = True
        else:
            # Select only True values entries for further compression
            rle_start = rle_start[rle_val]
            rle_end = rle_end[rle_val]

            # Create zip
            self.voxel_grid = zip(rle_start, rle_end)
            self.isEncoded_voxel_grid = True
    else:
        self.voxel_grid = voxel_grid
        self.isEncoded_voxel_grid = False

export(file_path)

Exports the image to the requested directory :param file_path: directory to write image to :return:

Source code in src/luna/radiology/mirp/imageClass.py

def export(self, file_path):
    """
    Exports the image to the requested directory
    :param file_path: directory to write image to
    :return:
    """

    # Skip if the image is missing
    if self.is_missing:
        return

    # Construct file name
    file_name = self.get_export_descriptor() + ".nii"

    # Export image to file
    return self.write(file_path=file_path, file_name=file_name)

get_export_descriptor()

Generates an image descriptor based on parameters of the image :return:

Source code in src/luna/radiology/mirp/imageClass.py

def get_export_descriptor(self):
    """
    Generates an image descriptor based on parameters of the image
    :return:
    """

    descr_list = ["volumes"]

    # Image name and spatial transformation
    if self.name is not None:
        descr_list += [self.name]
    if self.modality is not None:
        descr_list += [self.modality]

    if self.interpolated:
        # Interpolation
        descr_list += [
            self.interpolation_algorithm,
            "x",
            str(self.spacing[2])[:5],
            "y",
            str(self.spacing[1])[:5],
            "z",
            str(self.spacing[0])[:5],
        ]
    if self.binned:
        descr_list += ["bin", str(self.bin_width)]

    if self.rotation_angle != 0.0:
        # Rotation angle
        descr_list += ["rot", str(self.rotation_angle)[:5]]

    if not (
        self.transl_fraction_x == 0.0
        and self.transl_fraction_y == 0.0
        and self.transl_fraction_z == 0.0
    ):
        # Translation fraction
        descr_list += [
            "trans",
            "x",
            str(self.transl_fraction_x)[:5],
            "y",
            str(self.transl_fraction_y)[:5],
            "z",
            str(self.transl_fraction_z)[:5],
        ]

    if self.noise != -1.0:
        # Noise level
        descr_list += ["noise", str(self.noise)[:5], "iter", str(self.noise_iter)]

    if not self.spat_transform == "base":
        descr_list += [self.spat_transform]

    return "_".join(descr_list)

get_voxel_grid()

Return the voxel grid as a ndarray

Source code in src/luna/radiology/mirp/imageClass.py

def get_voxel_grid(self) -> np.ndarray:
    """Return the voxel grid as a ndarray"""
    if self.is_missing:
        return None

    if self.isEncoded_voxel_grid:
        # Decode voxel grid (typically roi)
        decoded_voxel = np.zeros(np.prod(self.size), dtype=np.bool)

        # Check if the voxel grid contains values
        if self.voxel_grid is not None:
            decode_zip = copy.deepcopy(self.voxel_grid)

            for ii, jj in decode_zip:
                decoded_voxel[ii : jj + 1] = True

        # Shape into correct form
        decoded_voxel = decoded_voxel.reshape(self.size)

        return decoded_voxel
    else:
        return self.voxel_grid

interpolate(by_slice, settings)

Performs interpolation of the image volume

Source code in src/luna/radiology/mirp/imageClass.py

def interpolate(self, by_slice, settings):
    """Performs interpolation of the image volume"""
    from luna.radiology.mirp.imageProcess import (  # aauker: Circular import
        gaussian_preprocess_filter,
        interpolate_to_new_grid,
    )

    # Skip for missing images
    if self.is_missing:
        return

    # Local interpolation constants
    if None not in settings.img_interpolate.new_spacing:
        iso_spacing = settings.img_interpolate.new_spacing[0]
        new_spacing = np.array(
            [iso_spacing, iso_spacing, iso_spacing]
        )  # Desired spacing in mm
    elif type(settings.img_interpolate.new_non_iso_spacing) in [list, tuple]:
        if None not in settings.img_interpolate.new_non_iso_spacing:
            non_iso_spacing = settings.img_interpolate.new_non_iso_spacing
            new_spacing = np.array(non_iso_spacing)
        else:
            new_spacing = self.spacing
    else:
        new_spacing = self.spacing

    print(f"Interpolating main image to {new_spacing}")

    # Read additional details
    order = (
        settings.img_interpolate.spline_order
    )  # Order of multidimensional spline filter (0=nearest neighbours, 1=linear, 3=cubic)
    interpolate_flag = (
        settings.img_interpolate.interpolate
    )  # Whether to interpolate or not

    # Set spacing for interpolation across slices to the original spacing in case interpolation is only conducted within the slice
    if by_slice:
        new_spacing[0] = self.spacing[0]

    # Image translation
    translate_z = 0  # settings.vol_adapt.translate_z[0]
    translate_y = 0  # settings.vol_adapt.translate_y[0]
    translate_x = 0  # settings.vol_adapt.translate_x[0]

    # Convert to [0.0, 1.0] range
    translate_x = translate_x - np.floor(translate_x)
    translate_y = translate_y - np.floor(translate_y)
    translate_z = translate_z - np.floor(translate_z)
    trans_vec = np.array([translate_z, translate_y, translate_x])

    # Add translation fractions
    self.transl_fraction_x = translate_x
    self.transl_fraction_y = translate_y
    self.transl_fraction_z = translate_z

    # Skip if translation in both directions is 0.0
    if (
        translate_x == 0.0
        and translate_y == 0.0
        and translate_z == 0.0
        and not interpolate_flag
    ):
        return None

    # Check if pre-processing is required
    if settings.img_interpolate.anti_aliasing:
        self.set_voxel_grid(
            voxel_grid=gaussian_preprocess_filter(
                orig_vox=self.get_voxel_grid(),
                orig_spacing=self.spacing,
                sample_spacing=new_spacing,
                param_beta=settings.img_interpolate.smoothing_beta,
                mode="nearest",
                by_slice=by_slice,
            )
        )

    # Interpolate image and positioning
    (
        self.size,
        sample_spacing,
        upd_voxel_grid,
        grid_origin,
    ) = interpolate_to_new_grid(
        orig_dim=self.size,
        orig_spacing=self.spacing,
        orig_vox=self.get_voxel_grid(),
        sample_spacing=new_spacing,
        translation=trans_vec,
        order=order,
        mode="nearest",
        align_to_center=True,
    )

    # Update origin before spacing, because computing the origin requires the original affine matrix.
    self.origin = self.origin + np.dot(self.m_affine, np.transpose(grid_origin))

    # Update spacing and affine matrix.
    self.set_spacing(sample_spacing)

    # Round intensities in case of modalities with inherently discretised intensities
    if (self.modality == "CT") and (self.spat_transform == "base"):
        upd_voxel_grid = np.round(upd_voxel_grid)
    elif (self.modality == "PT") and (self.spat_transform == "base"):
        upd_voxel_grid[upd_voxel_grid < 0.0] = 0.0

    # Set interpolation
    self.interpolated = True

    # Set interpolation algorithm
    if order == 0:
        self.interpolation_algorithm = "nnb"
    elif order == 1:
        self.interpolation_algorithm = "lin"
    elif order > 1:
        self.interpolation_algorithm = "si" + str(order)

    if settings.img_interpolate.bin:
        self.binned = True
        self.bin_width = settings.img_interpolate.bin_width
        self.bins = np.arange(-1000, 1000, settings.img_interpolate.bin_width)
        upd_voxel_grid = np.digitize(upd_voxel_grid, self.bins).astype(np.float)

    # Set voxel grid
    self.set_voxel_grid(voxel_grid=upd_voxel_grid)

normalise_intensities(norm_method='none', intensity_range=None, saturation_range=None, mask=None)

Normalises image intensities :param norm_method: string defining the normalisation method. Should be one of "none", "range", "standardisation" :param intensity_range: range of intensities for normalisation :return:

Source code in src/luna/radiology/mirp/imageClass.py

def normalise_intensities(
    self, norm_method="none", intensity_range=None, saturation_range=None, mask=None
):
    """
    Normalises image intensities
    :param norm_method: string defining the normalisation method. Should be one of "none", "range", "standardisation"
    :param intensity_range: range of intensities for normalisation
    :return:
    """

    # Skip for missing images
    if self.is_missing:
        return

    if intensity_range is None:
        intensity_range = [np.nan, np.nan]

    if mask is None:
        mask = np.ones(self.size, dtype=np.bool)
    else:
        mask = mask.astype(np.bool)

    if np.sum(mask) == 0:
        mask = np.ones(self.size, dtype=np.bool)

    if saturation_range is None:
        saturation_range = [np.nan, np.nan]

    if norm_method == "none":
        return

    elif norm_method == "range":
        # Normalisation to [0, 1] range using fixed intensities.

        # Get voxel grid
        voxel_grid = self.get_voxel_grid()

        # Find maximum and minimum intensities
        if np.isnan(intensity_range[0]):
            min_int = np.min(voxel_grid[mask])
        else:
            min_int = intensity_range[0]

        if np.isnan(intensity_range[1]):
            max_int = np.max(voxel_grid[mask])
        else:
            max_int = intensity_range[1]

        # Normalise by range
        if not max_int == min_int:
            voxel_grid = (voxel_grid - min_int) / (max_int - min_int)
        else:
            voxel_grid = voxel_grid - min_int

        # Update the voxel grid
        self.set_voxel_grid(voxel_grid=voxel_grid)

        self.is_normalised = True

    elif norm_method == "relative_range":
        # Normalisation to [0, 1]-ish range using relative intensities.

        # Get voxel grid
        voxel_grid = self.get_voxel_grid()

        min_int_rel = 0.0
        if not np.isnan(intensity_range[0]):
            min_int_rel = intensity_range[0]

        max_int_rel = 1.0
        if not np.isnan(intensity_range[1]):
            max_int_rel = intensity_range[1]

        # Compute minimum and maximum intensities.
        value_range = [np.min(voxel_grid[mask]), np.max(voxel_grid[mask])]
        min_int = value_range[0] + min_int_rel * (value_range[1] - value_range[0])
        max_int = value_range[0] + max_int_rel * (value_range[1] - value_range[0])

        # Normalise by range
        if not max_int == min_int:
            voxel_grid = (voxel_grid - min_int) / (max_int - min_int)
        else:
            voxel_grid = voxel_grid - min_int

        # Update the voxel grid
        self.set_voxel_grid(voxel_grid=voxel_grid)

        self.is_normalised = True

    elif norm_method == "quantile_range":
        # Normalisation to [0, 1]-ish range based on quantiles.

        # Get voxel grid
        voxel_grid = self.get_voxel_grid()

        min_quantile = 0.0
        if not np.isnan(intensity_range[0]):
            min_quantile = intensity_range[0]

        max_quantile = 1.0
        if not np.isnan(intensity_range[1]):
            max_quantile = intensity_range[1]

        # Compute quantiles from voxel grid.
        min_int = np.quantile(voxel_grid[mask], q=min_quantile)
        max_int = np.quantile(voxel_grid[mask], q=max_quantile)

        # Normalise by range
        if not max_int == min_int:
            voxel_grid = (voxel_grid - min_int) / (max_int - min_int)
        else:
            voxel_grid = voxel_grid - min_int

        # Update the voxel grid
        self.set_voxel_grid(voxel_grid=voxel_grid)

        self.is_normalised = True

    elif norm_method == "standardisation":
        # Normalisation to mean 0 and standard deviation 1.

        # Get voxel grid
        voxel_grid = self.get_voxel_grid()

        # Determine mean and standard deviation of the voxel intensities
        mean_int = np.mean(voxel_grid[mask])
        sd_int = np.std(voxel_grid[mask])

        # Protect against invariance.
        if sd_int == 0.0:
            sd_int = 1.0

        # Normalise
        voxel_grid = (voxel_grid - mean_int) / sd_int

        # Update the voxel grid
        self.set_voxel_grid(voxel_grid=voxel_grid)

        self.is_normalised = True
    else:
        raise ValueError(
            f"{norm_method} is not a valid method for normalising intensity values."
        )

    self.saturate(intensity_range=saturation_range)

rotate(angle)

Rotate volume along z-axis.

Source code in src/luna/radiology/mirp/imageClass.py

def rotate(self, angle):
    """Rotate volume along z-axis."""

    # Skip for missing images
    if self.is_missing:
        return

    import scipy.ndimage as ndi

    from luna.radiology.mirp.featureSets.volumeMorphology import get_rotation_matrix

    # Find actual output size of x-y plane
    new_z_dim = np.asmatrix([self.size[0], 0.0, 0.0]) * get_rotation_matrix(
        np.radians(angle), dim=3, rot_axis=0
    )
    new_y_dim = np.asmatrix([0.0, self.size[1], 0.0]) * get_rotation_matrix(
        np.radians(angle), dim=3, rot_axis=0
    )
    new_x_dim = np.asmatrix([0.0, 0.0, self.size[2]]) * get_rotation_matrix(
        np.radians(angle), dim=3, rot_axis=0
    )
    new_dim_flt = np.squeeze(
        np.array(np.abs(new_z_dim))
        + np.array(np.abs(new_y_dim) + np.abs(new_x_dim))
    )

    # Get voxel grid
    voxel_grid = self.get_voxel_grid()

    # Rotate voxels along angle in the y-x plane and find truncated output size
    voxel_grid = ndi.rotate(
        voxel_grid.astype(np.float32),
        angle=angle,
        axes=(1, 2),
        reshape=True,
        order=1,
        mode="nearest",
    )
    new_dim_int = np.array(np.shape(voxel_grid)) * 1.0

    if (self.modality == "CT") and (self.spat_transform == "base"):
        voxel_grid = np.round(voxel_grid)

    # Update spacing
    self.spacing *= new_dim_int / new_dim_flt

    # Set rotation angle
    self.rotation_angle = angle

    # Update voxel grid with rotated voxels
    self.set_voxel_grid(voxel_grid=voxel_grid)

saturate(intensity_range, fill_value=None)

Saturate image intensities using an intensity range :param intensity_range: range of intensity values :param fill_value: fill value for out-of-range intensities. If None, the upper and lower ranges are used :return:

Source code in src/luna/radiology/mirp/imageClass.py

def saturate(self, intensity_range, fill_value=None):
    """
    Saturate image intensities using an intensity range
    :param intensity_range: range of intensity values
    :param fill_value: fill value for out-of-range intensities. If None, the upper and lower ranges are used
    :return:
    """
    # Skip for missing images
    if self.is_missing:
        return

    intensity_range = np.array(copy.deepcopy(intensity_range))

    if np.any(~np.isnan(intensity_range)):
        # Get voxel grid
        voxel_grid = self.get_voxel_grid()

        # Lower boundary
        if not np.isnan(intensity_range[0]):
            if fill_value is None:
                voxel_grid[voxel_grid < intensity_range[0]] = intensity_range[0]
            else:
                voxel_grid[voxel_grid < intensity_range[0]] = fill_value[0]

        # Upper boundary
        if not np.isnan(intensity_range[1]):
            if fill_value is None:
                voxel_grid[voxel_grid > intensity_range[1]] = intensity_range[1]
            else:
                voxel_grid[voxel_grid > intensity_range[1]] = fill_value[1]

        # Set the updated voxel grid
        self.set_voxel_grid(voxel_grid=voxel_grid)

set_voxel_grid(voxel_grid)

Sets voxel grid

Source code in src/luna/radiology/mirp/imageClass.py

def set_voxel_grid(self, voxel_grid):
    """Sets voxel grid"""

    # Determine size
    self.size = np.array(voxel_grid.shape)
    self.dtype_name = voxel_grid.dtype.name

    # Encode voxel grid
    self.encode_voxel_grid(voxel_grid=voxel_grid)

write(file_path, file_name)

Writes the image to a file

Source code in src/luna/radiology/mirp/imageClass.py

def write(self, file_path, file_name):
    """Writes the image to a file"""
    import os

    import SimpleITK as sitk

    sitk.ProcessObject.SetGlobalDefaultNumberOfThreads(4)

    # Check if path exists
    file_path = os.path.normpath(file_path)
    if not os.path.exists(file_path):
        os.makedirs(file_path)

    # Add file and file name
    file_path = os.path.join(file_path, file_name)

    # Convert image to simple itk format
    sitk_img = self.convert2sitk()

    # Write to file using simple itk, and the image is not missing
    if sitk_img is not None:
        print(f"Writing to {file_path}")
        sitk.WriteImage(sitk_img, file_path, True)

    return file_path

imageMetaData

parse_image_correction(dcm_seq, tag, correction_abbr, image_corrections)

Parses image correction information :param dcm_seq: pydicom object :param tag: tag for image correction :param correction_abbr: abbreviation for the specific correction in the image_corrections list :param image_corrections: list of image corrections from DICOM tag (0028,0051) :return: whether an image correction was performed ("YES", "NO" or "")

This allows determining image correction with a single line call

Source code in src/luna/radiology/mirp/imageMetaData.py

def parse_image_correction(dcm_seq, tag, correction_abbr, image_corrections):
    """
    Parses image correction information
    :param dcm_seq: pydicom object
    :param tag: tag for image correction
    :param correction_abbr: abbreviation for the specific correction in the image_corrections list
    :param image_corrections: list of image corrections from DICOM tag (0028,0051)
    :return: whether an image correction was performed ("YES", "NO" or "")

    This allows determining image correction with a single line call
    """

    is_corrected = get_pydicom_meta_tag(
        dcm_seq=dcm_seq, tag=tag, tag_type="str", default=None
    )
    if is_corrected is None and correction_abbr in image_corrections:
        is_corrected = "YES"
    elif is_corrected is None and correction_abbr not in image_corrections:
        is_corrected = "NO"
    else:
        is_corrected = ""

    return is_corrected

imagePerturbations

adapt_roi_size(roi_list, settings)

Adapt roi size by growing or shrinking the roi

Source code in src/luna/radiology/mirp/imagePerturbations.py

def adapt_roi_size(roi_list, settings):
    """Adapt roi size by growing or shrinking the roi"""

    # Adapt roi size by shrinking or increasing the roi
    new_roi_list = []

    # Get the adaptation size and type. Rois with adapt_size > 0.0 are dilated. Rois with adapt_size < 0.0 are eroded.
    # The type determines whether the roi is grown/shrunk with by certain distance ("distance") or to a certain volume fraction ("fraction")
    adapt_size_list = settings.vol_adapt.roi_adapt_size
    adapt_type = settings.vol_adapt.roi_adapt_type

    # Iterate over roi objects in the roi list and adaptation sizes
    for roi_obj in roi_list:
        for adapt_size in adapt_size_list:
            if adapt_size > 0.0 and adapt_type == "distance":
                new_roi_obj = roi_obj.copy()
                new_roi_obj.dilate(by_slice=settings.general.by_slice, dist=adapt_size)

                # Update name and adaptation size
                new_roi_obj.name += "_grow" + str(adapt_size)
                new_roi_obj.adapt_size = adapt_size

                # Add to roi list
                new_roi_list += [new_roi_obj]

            elif adapt_size < 0.0 and adapt_type == "distance":
                new_roi_obj = roi_obj.copy()
                new_roi_obj.erode(
                    by_slice=settings.general.by_slice,
                    dist=adapt_size,
                    eroded_vol_fract=settings.vol_adapt.eroded_vol_fract,
                )

                # Update name and adaptation size
                new_roi_obj.name += "_shrink" + str(np.abs(adapt_size))
                new_roi_obj.adapt_size = adapt_size

                # Add to roi list
                new_roi_list += [new_roi_obj]

            elif adapt_type == "fraction" and not adapt_size == 0.0:
                new_roi_obj = roi_obj.copy()
                new_roi_obj.adapt_volume(
                    by_slice=settings.general.by_slice, vol_grow_fract=adapt_size
                )

                # Update name and adaptation size
                if adapt_size > 0:
                    new_roi_obj.name += "_grow" + str(adapt_size)
                else:
                    new_roi_obj.name += "_shrink" + str(np.abs(adapt_size))
                new_roi_obj.adapt_size = adapt_size

                # Add to roi list
                new_roi_list += [new_roi_obj]

            else:
                new_roi_list += [roi_obj]

    # Check for non-updated rois
    roi_names = extract_roi_names(new_roi_list)
    uniq_roi_names, uniq_index, uniq_counts = np.unique(
        np.asarray(roi_names), return_index=True, return_counts=True
    )
    if np.size(uniq_index) != len(roi_names):
        uniq_roi_list = [new_roi_list[ii] for ii in uniq_index]
    else:
        uniq_roi_list = new_roi_list

    # Return expanded roi list
    return uniq_roi_list

randomise_roi_contours(img_obj, roi_list, settings)

Use SLIC to randomise the roi based on supervoxels

Source code in src/luna/radiology/mirp/imagePerturbations.py

def randomise_roi_contours(img_obj, roi_list, settings):
    """Use SLIC to randomise the roi based on supervoxels"""

    # Check whether randomisation should take place
    if not settings.vol_adapt.randomise_roi:
        return roi_list

    from scipy.ndimage import binary_closing

    from luna.radiology.mirp.utilities import world_to_index

    new_roi_list = []
    svx_roi_list = []

    # Iterate over roi objects
    for roi_ind in np.arange(0, len(roi_list)):
        print(f">>> Processing ROI with label [{roi_list[roi_ind].label_value}]")

        # Resect image to speed up segmentation process
        res_img_obj, res_roi_obj = crop_image(
            img_obj=img_obj, roi_obj=roi_list[roi_ind], boundary=5.0, z_only=False
        )

        print(f"Res_roi_obj shape = {res_roi_obj.roi.size}")
        # Calculate statistics on post-processed, cropped ROI
        res_roi_obj.calculate_roi_statistics(img_obj=res_img_obj, tag="postprocess")

        # tumor_volume     = res_roi_obj.roi.get_voxel_grid().sum() * np.prod(img_obj.spacing)
        # tumor_volume_1up = binary_dilation(res_roi_obj.roi.get_voxel_grid()).sum() * np.prod(img_obj.spacing)
        # tumor_surface_area = tumor_volume_1up-tumor_volume
        # print ("Volume, Differential Volume: ", tumor_volume, tumor_surface_area)

        # min_n_voxels = np.max([20.0, 250.0 / np.prod(res_img_obj.spacing)])
        # segment_guess =  int(np.prod(res_img_obj.size) / min_n_voxels)
        # print ("Starting guess: ", segment_guess)

        # for n_segments in np.linspace(segment_guess, segment_guess*5, 50):
        #     # Get supervoxels
        #     n_segments = int(n_segments)

        img_segments = get_supervoxels(
            img_obj=res_img_obj, roi_obj=res_roi_obj, settings=settings, n_segments=None
        )

        # Determine overlap of supervoxels with contour
        overlap_indices, overlap_fract, overlap_size = get_supervoxel_overlap(
            roi_obj=res_roi_obj, img_segments=img_segments
        )

        # Set the highest overlap to 1.0 to ensure selection of at least 1 supervoxel
        # aauker: aka, highest overlapping supervoxel is always included
        overlap_fract[np.argmax(overlap_fract)] = 1.0

        # Include supervoxels with 80% coverage and exclude those with less then 20% coverage
        a = 0.80
        b = 0.20

        overlap_fract[overlap_fract > a] = 1.0
        overlap_fract[overlap_fract < b] = 0.0

        candidate_indices = overlap_indices[
            np.logical_and(overlap_fract > 0.0, overlap_fract < 1.0)
        ]
        candidate_segments = np.where(
            np.isin(img_segments, candidate_indices), img_segments, 0
        )

        average_segment_size = (
            np.prod(img_obj.spacing)
            * np.where(candidate_segments > 0, 1, 0).sum()
            / len(candidate_indices)
        )

        print(f"Average segment size: {average_segment_size}")

        # if average_segment_size < 250: break

        # break # Use initial guess...for now

        print("Candidate segments: ", len(candidate_indices))

        # Determine grid indices of the resected grid with respect to the original image grid
        grid_origin = world_to_index(
            coord=res_img_obj.origin,
            origin=img_obj.origin,
            spacing=img_obj.spacing,
            affine=img_obj.m_affine,
        )
        grid_origin = grid_origin.astype(np.int)

        # Iteratively create randomised regions of interest
        for ii in np.arange(settings.vol_adapt.roi_random_rep):
            # Draw random numbers between 0.0 and 1.0
            # Use numpy rng for deterministic behavior (updated)
            random_incl = np.random.default_rng(ii).random(size=len(overlap_fract))

            # Select those segments where the random number is less than the overlap fraction - i.e. the fraction is the
            # probability of selecting the supervoxel
            incl_segments = overlap_indices[np.less(random_incl, overlap_fract)]

            # Replace randomised contour in original roi voxel space
            roi_vox = np.zeros(shape=roi_list[roi_ind].roi.size, dtype=np.bool)
            roi_vox[
                grid_origin[0] : grid_origin[0] + res_roi_obj.roi.size[0],
                grid_origin[1] : grid_origin[1] + res_roi_obj.roi.size[1],
                grid_origin[2] : grid_origin[2] + res_roi_obj.roi.size[2],
            ] = np.reshape(
                np.in1d(np.ravel(img_segments), incl_segments), res_roi_obj.roi.size
            )

            # Apply binary closing to close gaps
            roi_vox = binary_closing(input=roi_vox)

            # Update voxels in original roi, adapt name and set randomisation id

            repl_roi = roi_list[roi_ind].copy()
            repl_roi.roi.set_voxel_grid(
                voxel_grid=roi_vox
            )  # Replace copied original contour with randomised contour
            repl_roi.name += "_svx_" + str(ii)  # Adapt roi name
            repl_roi.svx_randomisation_id = ii  # Update randomisation id
            new_roi_list += [repl_roi]

        # Update voxels in original roi, adapt name and set randomisation id
        # Replace randomised contour in original roi voxel space
        roi_vox = np.zeros(shape=roi_list[roi_ind].roi.size, dtype=np.uint8)
        roi_vox[
            grid_origin[0] : grid_origin[0] + res_roi_obj.roi.size[0],
            grid_origin[1] : grid_origin[1] + res_roi_obj.roi.size[1],
            grid_origin[2] : grid_origin[2] + res_roi_obj.roi.size[2],
        ] = candidate_segments

        repl_roi = roi_list[roi_ind].copy()
        repl_roi.roi.set_voxel_grid(
            voxel_grid=roi_vox
        )  # Replace copied original contour with randomised contour
        repl_roi.name += "_SUPERVOXEL"  # Adapt roi name
        repl_roi.svx_randomisation_id = -1  # Update randomisation id
        svx_roi_list += [repl_roi]

    return new_roi_list, svx_roi_list

rotate_image(img_obj, settings=None, rot_angle=None, roi_list=None)

Rotation of image and rois

Source code in src/luna/radiology/mirp/imagePerturbations.py

def rotate_image(img_obj, settings=None, rot_angle=None, roi_list=None):
    """Rotation of image and rois"""

    if settings is not None:
        rot_angle = settings.vol_adapt.rot_angles
    elif rot_angle is None:
        logging.error(
            "No rotation angles were provided. A single rotation angle is expected."
        )

    if len(rot_angle) > 1:
        logging.warning(
            "Multiple rotation angles were provided. Only the first is selected."
        )

    if type(rot_angle) is list:
        rot_angle = rot_angle[0]

    if rot_angle in [0.0, 360.0]:
        return img_obj, roi_list

    # Rotate rois
    if roi_list is not None:
        for ii in np.arange(0, len(roi_list)):
            roi_list[ii].rotate(angle=rot_angle, img_obj=img_obj)

    # Rotate image object
    img_obj.rotate(angle=rot_angle)

    return img_obj, roi_list

imageProcess

crop_image(img_obj, roi_list=None, roi_obj=None, boundary=0.0, z_only=False)

The function is used to slice a subsection of the image so that further processing is facilitated in terms of memory and computational requirements.

Source code in src/luna/radiology/mirp/imageProcess.py

def crop_image(img_obj, roi_list=None, roi_obj=None, boundary=0.0, z_only=False):
    """The function is used to slice a subsection of the image so that further processing is facilitated in terms of
    memory and computational requirements."""

    ####################################################################################################################
    # Initial steps
    ####################################################################################################################

    # Temporarily parse roi_obj to list, if roi_obj is provided and not roi_list. This is done for easier code maintenance.
    if roi_list is None:
        roi_list = [roi_obj]
        return_roi_obj = True
    else:
        return_roi_obj = False

    ####################################################################################################################
    # Determine region of interest bounding box
    ####################################################################################################################
    roi_ext_x = []
    roi_ext_y = []
    roi_ext_z = []

    # Determine extent of all rois
    for roi_obj in roi_list:
        # Skip if the ROI is missing
        if roi_obj.roi is None:
            continue

        z_ind, y_ind, x_ind = np.where(roi_obj.roi.get_voxel_grid() > 0.0)

        # Skip if the ROI is empty
        if len(z_ind) == 0 or len(y_ind) == 0 or len(x_ind) == 0:
            continue

        roi_ext_z += [np.min(z_ind), np.max(z_ind)]
        roi_ext_y += [np.min(y_ind), np.max(y_ind)]
        roi_ext_x += [np.min(x_ind), np.max(x_ind)]

    # Check if the combined ROIs are empty
    if not (len(roi_ext_z) == 0 or len(roi_ext_y) == 0 or len(roi_ext_x) == 0):
        # Express mm boundary in voxels.
        boundary = np.ceil(boundary / img_obj.spacing).astype(np.int)

        # Concatenate extents for rois and add boundary to generate map extent
        ind_ext_z = np.array(
            [np.min(roi_ext_z) - boundary[0], np.max(roi_ext_z) + boundary[0]]
        )
        ind_ext_y = np.array(
            [np.min(roi_ext_y) - boundary[1], np.max(roi_ext_y) + boundary[1]]
        )
        ind_ext_x = np.array(
            [np.min(roi_ext_x) - boundary[2], np.max(roi_ext_x) + boundary[2]]
        )

        ####################################################################################################################
        # Resect image based on roi extent
        ####################################################################################################################

        img_res = img_obj.copy()
        img_res.crop(
            ind_ext_z=ind_ext_z, ind_ext_y=ind_ext_y, ind_ext_x=ind_ext_x, z_only=z_only
        )

        ####################################################################################################################
        # Resect rois based on roi extent
        ####################################################################################################################

        # Copy roi objects before resection
        roi_res_list = [roi_res_obj.copy() for roi_res_obj in roi_list]

        # Resect in place
        [
            roi_res_obj.crop(
                ind_ext_z=ind_ext_z,
                ind_ext_y=ind_ext_y,
                ind_ext_x=ind_ext_x,
                z_only=z_only,
            )
            for roi_res_obj in roi_res_list
        ]

    else:
        # This happens if all rois are empty - only copies of the original image object and the roi are returned
        img_res = img_obj.copy()
        roi_res_list = [roi_res_obj.copy() for roi_res_obj in roi_list]

    ####################################################################################################################
    # Return to calling function
    ####################################################################################################################

    if return_roi_obj:
        return img_res, roi_res_list[0]
    else:
        return img_res, roi_res_list

get_supervoxel_overlap(roi_obj, img_segments, mask=None)

Determines overlap of supervoxels with other the region of interest

Source code in src/luna/radiology/mirp/imageProcess.py

def get_supervoxel_overlap(roi_obj, img_segments, mask=None):
    """Determines overlap of supervoxels with other the region of interest"""

    # Return None in case image segments and/or ROI are missing
    if img_segments is None or roi_obj.roi is None:
        return None, None, None

    # Check segments overlapping with the current contour
    if mask == "morphological" and roi_obj.roi_morphology is not None:
        overlap_segment_labels, overlap_size = np.unique(
            np.multiply(img_segments, roi_obj.roi_morphology.get_voxel_grid()),
            return_counts=True,
        )
    elif mask == "intensity" and roi_obj.roi_intensity is not None:
        overlap_segment_labels, overlap_size = np.unique(
            np.multiply(img_segments, roi_obj.roi_intensity.get_voxel_grid()),
            return_counts=True,
        )
    else:
        overlap_segment_labels, overlap_size = np.unique(
            np.multiply(img_segments, roi_obj.roi.get_voxel_grid()), return_counts=True
        )

    # Find super voxels with non-zero overlap with the roi
    overlap_size = overlap_size[overlap_segment_labels > 0]
    overlap_segment_labels = overlap_segment_labels[overlap_segment_labels > 0]

    if len(overlap_size) == 0:
        raise RuntimeError(
            "No valid supervoxels found, this can happen if the entire grid recieved label 0 from slic, did you window out your tumor?"
        )

    # Check the actual size of the segments overlapping with the current contour
    segment_size = list(
        map(lambda x: np.sum([img_segments == x]), overlap_segment_labels)
    )

    # Calculate the fraction of overlap
    overlap_frac = overlap_size / segment_size

    return overlap_segment_labels, overlap_frac, overlap_size

get_supervoxels(img_obj, roi_obj, settings, n_segments=None)

Extracts supervoxels from an image

Source code in src/luna/radiology/mirp/imageProcess.py

def get_supervoxels(img_obj, roi_obj, settings, n_segments=None):
    """Extracts supervoxels from an image"""

    # Check if image and/or roi exist, and skip otherwise
    if img_obj.is_missing or roi_obj.roi is None:
        print(
            "You tried to get supervoxel labels, however didn't pass in an image, so be careful!"
        )
        return None

    # Get image object grid aauker: I don't think this is neccessary anymore
    img_voxel_grid = img_obj.get_voxel_grid()
    roi_bool_mask = roi_obj.roi.get_voxel_grid().astype(np.bool)

    outside = binary_dilation(
        roi_bool_mask, iterations=int(10 // np.min(img_obj.spacing))
    )

    min_n_voxels = np.max(
        [20.0, 100.0 / np.prod(img_obj.spacing)]
    )  # Even smaller super voxels....
    segment_guess = int(np.sum(outside) / min_n_voxels)
    print("Starting guess: ", segment_guess)

    # inside  = binary_erosion(bool_mask,  iterations=10)
    # ring_mask = np.logical_and( outside, np.invert(inside))

    # Calculate roi pixel grid, low-level pixels, high-level pixels
    roi_level_pixels = np.ma.array(img_voxel_grid, mask=np.invert(roi_bool_mask))
    low_level, high_level = (
        roi_level_pixels.min() - 3 * roi_level_pixels.std(),
        roi_level_pixels.max() + 3 * roi_level_pixels.std(),
    )
    img_voxel_grid = np.clip(img_voxel_grid, low_level, high_level)

    # Convert to float with range [0.0, 1.0]
    img_voxel_grid = (img_voxel_grid - np.min(img_voxel_grid)) / np.ptp(img_voxel_grid)

    # Just also convert to floats
    img_voxel_grid = img_voxel_grid.astype(np.float)

    # Create a slic segmentation of the image stack
    img_segments = slic(
        image=img_voxel_grid,
        n_segments=segment_guess,
        spacing=img_obj.spacing,
        mask=outside,
        max_iter=50,
        sigma=1.0,
        compactness=0.04,
        multichannel=False,
        convert2lab=False,
        enforce_connectivity=True,
        start_label=1,
    )

    return img_segments

interpolate_to_new_grid(orig_dim, orig_spacing, orig_vox, sample_dim=None, sample_spacing=None, grid_origin=None, translation=np.array([0.0, 0.0, 0.0]), order=1, mode='nearest', align_to_center=True, processor='scipy')

Resamples input grid and returns the output grid. :param orig_dim: dimensions of the input grid :param orig_origin: origin (in world coordinates) of the input grid :param orig_spacing: spacing (in world measures) of the input grid :param orig_vox: input grid :param sample_dim: desired output size (determined within the function if None) :param sample_origin: desired output origin (in world coordinates; determined within the function if None) :param sample_spacing: desired sample spacing (in world measures; should be provided if sample_dim or sample_origin is None) :param translation: a translation vector that is used to shift the interpolation grid (in voxel measures) :param order: interpolation spline order (0=nnb, 1=linear, 2=order 2 spline, 3=cubic splice, max 5). :param mode: describes how to handle extrapolation beyond input grid. :param align_to_center: whether the input and output grids should be aligned by their centers (True) or their origins (False) :param processor: which function to use for interpolation: "scipy" for scipy.ndimage.map_coordinates and "sitk" for SimpleITK.ResampleImageFilter :return:

Source code in src/luna/radiology/mirp/imageProcess.py

def interpolate_to_new_grid(
    orig_dim,
    orig_spacing,
    orig_vox,
    sample_dim=None,
    sample_spacing=None,
    grid_origin=None,
    translation=np.array([0.0, 0.0, 0.0]),
    order=1,
    mode="nearest",
    align_to_center=True,
    processor="scipy",
):
    """
    Resamples input grid and returns the output grid.
    :param orig_dim: dimensions of the input grid
    :param orig_origin: origin (in world coordinates) of the input grid
    :param orig_spacing: spacing (in world measures) of the input grid
    :param orig_vox: input grid
    :param sample_dim: desired output size (determined within the function if None)
    :param sample_origin: desired output origin (in world coordinates; determined within the function if None)
    :param sample_spacing: desired sample spacing (in world measures; should be provided if sample_dim or sample_origin is None)
    :param translation: a translation vector that is used to shift the interpolation grid (in voxel measures)
    :param order: interpolation spline order (0=nnb, 1=linear, 2=order 2 spline, 3=cubic splice, max 5).
    :param mode: describes how to handle extrapolation beyond input grid.
    :param align_to_center: whether the input and output grids should be aligned by their centers (True) or their origins (False)
    :param processor: which function to use for interpolation: "scipy" for scipy.ndimage.map_coordinates and "sitk" for SimpleITK.ResampleImageFilter
    :return:
    """

    # Check if sample spacing is provided
    if sample_dim is None and sample_spacing is None:
        logging.error("Sample spacing is required for interpolation, but not provided.")

    # If no sample spacing is provided, assume original spacing. Note that for most purposes sample spacing should be provided
    if sample_spacing is None:
        sample_spacing = orig_spacing

    # Set sample spacing and orig_spacing to float
    sample_spacing = sample_spacing.astype(np.float)
    orig_spacing = orig_spacing.astype(np.float)

    # If no sample dimensions are provided, assume that the user wants to sample the original grid
    if sample_dim is None:
        sample_dim = np.ceil(np.multiply(orig_dim, orig_spacing / sample_spacing))

    # Set grid spacing (i.e. a fractional spacing in input voxel dimensions)
    grid_spacing = sample_spacing / orig_spacing

    # Set grid origin, if not provided previously
    if grid_origin is None:
        if align_to_center:
            grid_origin = (
                0.5 * (np.array(orig_dim) - 1.0)
                - 0.5 * (np.array(sample_dim) - 1.0) * grid_spacing
            )

        else:
            grid_origin = np.array([0.0, 0.0, 0.0])

        # Update with translation vector
        grid_origin += translation * grid_spacing

    if processor == "scipy":
        import scipy.ndimage as ndi

        # Convert sample_spacing and sample_origin to normalised original spacing (where voxel distance is 1 in each direction)
        # This is required for the use of ndi.map_coordinates, which uses the original grid as reference.
        # Generate interpolation map grid
        map_z, map_y, map_x = np.mgrid[
            : sample_dim[0], : sample_dim[1], : sample_dim[2]
        ]

        # Transform map to normalised original space
        map_z = map_z * grid_spacing[0] + grid_origin[0]
        map_z = map_z.astype(np.float32)
        map_y = map_y * grid_spacing[1] + grid_origin[1]
        map_y = map_y.astype(np.float32)
        map_x = map_x * grid_spacing[2] + grid_origin[2]
        map_x = map_x.astype(np.float32)

        # Interpolate orig_vox on interpolation grid
        map_vox = ndi.map_coordinates(
            input=orig_vox.astype(np.float32),
            coordinates=np.array([map_z, map_y, map_x], dtype=np.float32),
            order=order,
            mode=mode,
        )

    elif processor == "sitk":
        import SimpleITK as sitk

        sitk.ProcessObject.SetGlobalDefaultNumberOfThreads(4)

        # Convert input voxel grid to sitk image. Note that SimpleITK expects x,y,z ordering, while we use z,y,
        # x ordering. Hence origins, spacings and sizes are inverted for both input image (sitk_orig_img) and
        # ResampleImageFilter objects.
        sitk_orig_img = sitk.GetImageFromArray(
            orig_vox.astype(np.float32), isVector=False
        )
        sitk_orig_img.SetOrigin(np.array([0.0, 0.0, 0.0]))
        sitk_orig_img.SetSpacing(np.array([1.0, 1.0, 1.0]))

        interpolator = sitk.ResampleImageFilter()

        # Set interpolator algorithm; SimpleITK has more interpolators, but for now use the older scheme for scipy.
        if order == 0:
            interpolator.SetInterpolator(sitk.sitkNearestNeighbor)
        elif order == 1:
            interpolator.SetInterpolator(sitk.sitkLinear)
        elif order == 2:
            interpolator.SetInterpolator(sitk.sitkBSplineResamplerOrder2)
        elif order == 3:
            interpolator.SetInterpolator(sitk.sitkBSpline)

        # Set output origin and output spacing
        interpolator.SetOutputOrigin(grid_origin[::-1])
        interpolator.SetOutputSpacing(grid_spacing[::-1])
        interpolator.SetSize(sample_dim[::-1].astype(int).tolist())

        map_vox = sitk.GetArrayFromImage(interpolator.Execute(sitk_orig_img))
    else:
        raise ValueError('The selected processor should be one of "scipy" or "sitk"')

    # Return interpolated grid and spatial coordinates
    return sample_dim, sample_spacing, map_vox, grid_origin

transform_images(img_obj, roi_list, settings, compute_features=False, extract_images=False, file_path=None)

Performs image transformations and calculates features. :param img_obj: image object :param roi_list: list of region of interest objects :param settings: configuration settings :param compute_features: flag to enable feature computation :param extract_images: flag to enable image exports :param file_path: path for image exports :return: list of features computed in the transformed image

Source code in src/luna/radiology/mirp/imageProcess.py

def transform_images(
    img_obj,
    roi_list,
    settings,
    compute_features=False,
    extract_images=False,
    file_path=None,
):
    """
    Performs image transformations and calculates features.
    :param img_obj: image object
    :param roi_list: list of region of interest objects
    :param settings: configuration settings
    :param compute_features: flag to enable feature computation
    :param extract_images: flag to enable image exports
    :param file_path: path for image exports
    :return: list of features computed in the transformed image
    """

    # Empty list for storing features
    feat_list = []

    # Check if image transformation is required
    if not settings.img_transform.perform_img_transform:
        return feat_list

    # Get spatial filters to apply
    spatial_filter = settings.img_transform.spatial_filters

    # Iterate over spatial filters
    for curr_filter in spatial_filter:
        if curr_filter == "wavelet":
            # Wavelet filters
            from mirp.imageFilters.waveletFilter import WaveletFilter

            filter_obj = WaveletFilter(settings=settings)
            feat_list += filter_obj.apply_transformation(
                img_obj=img_obj,
                roi_list=roi_list,
                settings=settings,
                compute_features=compute_features,
                extract_images=extract_images,
                file_path=file_path,
            )

        elif curr_filter == "laplacian_of_gaussian":
            # Laplacian of Gaussian filters
            from mirp.imageFilters.laplacianOfGaussian import LaplacianOfGaussianFilter

            filter_obj = LaplacianOfGaussianFilter(settings=settings)
            feat_list += filter_obj.apply_transformation(
                img_obj=img_obj,
                roi_list=roi_list,
                settings=settings,
                compute_features=compute_features,
                extract_images=extract_images,
                file_path=file_path,
            )

        elif curr_filter == "laws":
            # Laws' kernels
            from mirp.imageFilters.lawsFilter import LawsFilter

            filter_obj = LawsFilter(settings=settings)
            feat_list += filter_obj.apply_transformation(
                img_obj=img_obj,
                roi_list=roi_list,
                settings=settings,
                compute_features=compute_features,
                extract_images=extract_images,
                file_path=file_path,
            )

        elif curr_filter == "mean":
            # Mean / uniform filter
            from mirp.imageFilters.meanFilter import MeanFilter

            filter_obj = MeanFilter(settings=settings)
            feat_list += filter_obj.apply_transformation(
                img_obj=img_obj,
                roi_list=roi_list,
                settings=settings,
                compute_features=compute_features,
                extract_images=extract_images,
                file_path=file_path,
            )

        else:
            raise ValueError(
                f"{curr_filter} is not implemented as a spatial filter. Please use one of wavelet, laplacian_of_gaussian, mean or laws."
            )

    return feat_list

imageReaders

read_itk_image(path_to_itk_file, modality=None)

This takes a itk volume in, and spits out this "ImageClass thing

Source code in src/luna/radiology/mirp/imageReaders.py

def read_itk_image(path_to_itk_file, modality=None):
    """This takes a itk volume in, and spits out this "ImageClass thing"""

    print("Loading ", path_to_itk_file)

    # Load the image
    sitk_img = sitk.ReadImage(path_to_itk_file)

    # Import the image volume
    voxel_grid = sitk.GetArrayFromImage(sitk_img).astype(np.float)

    # Determine origin, spacing, and orientation
    image_origin = np.array(sitk_img.GetOrigin())[::-1]
    image_spacing = np.array(sitk_img.GetSpacing())[::-1]
    image_orientation = np.array(sitk_img.GetDirection())[::-1]

    # Create an ImageClass object from the input image.
    image_obj = ImageClass(
        voxel_grid=voxel_grid,
        origin=image_origin,
        spacing=image_spacing,
        orientation=image_orientation,
        modality=modality,
        spat_transform="base",
        no_image=False,
    )

    return image_obj

roiClass

RoiClass

Source code in src/luna/radiology/mirp/roiClass.py

class RoiClass:
    # Class for regions of interest

    def __init__(
        self,
        name,
        contour,
        roi_mask=None,
        g_range=np.array([np.nan, np.nan]),
        incl_threshold=0.5,
        metadata=None,
        label_value=1,
    ):
        self.name = name
        self.label_value = int(label_value)

        # Fixed parameters
        self.g_range = g_range  # Range of allowed grey level intensities
        self.incl_threshold = incl_threshold  # Threshold for partial volume effect
        self.adapt_size = 0.0  # Shrinkage and growth of roi
        self.svx_randomisation_id = (
            -1
        )  # Randomisation id for supervoxel roi randomisation

        # ROI masks
        self.roi: Union[
            ImageClass, None
        ] = roi_mask  # Union of intensity and morphology masks
        self.roi_intensity: Union[ImageClass, None] = None  # Intensity mask of the ROI
        self.roi_morphology: Union[
            ImageClass, None
        ] = None  # Morphological mask of the ROI

        # Diagnostics features
        self.diagnostic_list = []
        self.metadata: FileDataset = metadata
        self.stats = {}

    def calculate_roi_statistics(self, img_obj, tag="none"):
        masked_voxel_array = np.ma.array(
            img_obj.get_voxel_grid(), mask=~self.roi.get_voxel_grid()
        )  # Invert binary roi grid
        skew = float(mstats.skew(masked_voxel_array, axis=None))
        #  Compute skew factors
        self.stats[tag] = {
            "roi_level_mu": masked_voxel_array.mean(),
            "roi_level_std": masked_voxel_array.std(),
            "roi_level_range": masked_voxel_array.ptp(),
            "roi_level_min": masked_voxel_array.min(),
            "roi_level_max": masked_voxel_array.max(),
            "roi_level_skew": skew,
        }

    def get_roi_statistics(self, tag="none"):
        return self.stats[tag]

    def copy(self):
        # Creates a new copy of the roi
        return copy.deepcopy(self)

    def create_mask_from_contours(
        self,
        img_obj,
        draw_method="ray_cast",
        disconnected_segments="keep_as_is",
        settings=None,
    ):
        # Creates an image based on provided contours

        # Skip if image object is empty
        if img_obj.is_missing:
            self.roi = None
            return

        if settings is not None:
            disconnected_segments = settings.general.divide_disconnected_roi

        # Create an empty roi volume
        roi_mask = np.zeros(img_obj.size, dtype=np.bool)

        # Iterate over contours to fill out the mask
        for contour in self.contour:
            # Multiple methods are implemented. All methods return a slice_list (containing slice numbers (z)) and a mask list, which contain boolean masks for respective
            # slices. This are then inserted at the specified slice positions, using an OR operation. This operation is required to avoid overwriting different slices.

            # Ray casting method to draw segmentation map based on polygon contour
            if draw_method == "ray_cast":
                slice_list, mask_list = contour.contour_to_grid_ray_cast(
                    img_obj=img_obj
                )
                for ii in np.arange(len(slice_list)):
                    slice_id = slice_list[ii]
                    roi_mask[slice_id, :, :] = np.logical_or(
                        roi_mask[slice_id, :, :], mask_list[ii]
                    )

        if disconnected_segments == "keep_largest":
            # Check if the created roi mask consists of multiple, separate segments, and keep only the largest.
            import skimage.measure

            # Label regions
            roi_label_mask, n_regions = skimage.measure.label(
                input=roi_mask, connectivity=2, return_num=True
            )

            # Determine size of regions
            roi_sizes = np.zeros(n_regions)
            for ii in np.arange(start=0, stop=n_regions):
                roi_sizes[ii] = np.sum(roi_label_mask == ii + 1)

            # Select largest region
            roi_mask = roi_label_mask == np.argmax(roi_sizes) + 1

        # Store roi as image object
        self.roi = ImageClass(
            voxel_grid=roi_mask,
            origin=img_obj.origin,
            spacing=img_obj.spacing,
            orientation=img_obj.orientation,
        )

        # Remove contour information
        self.contour = None

    def decimate(self, by_slice):
        """
        Decimates the roi
        :param by_slice: boolean, 2D (True) or 3D (False)
        :return:
        """

        # Resect masks
        if self.roi is not None:
            self.roi.decimate(by_slice=by_slice)
        if self.roi_intensity is not None:
            self.roi_intensity.decimate(by_slice=by_slice)
        if self.roi_morphology is not None:
            self.roi_morphology.decimate(by_slice=by_slice)

    def interpolate(self, img_obj, settings):
        print("Resampling ROI ", self.roi.spacing, " -> ", img_obj.spacing)
        # Skip if image and/or is missing
        if img_obj is None or self.roi is None:
            return

        if settings.img_interpolate.anti_aliasing:
            self.roi.set_voxel_grid(
                voxel_grid=gaussian_preprocess_filter(
                    orig_vox=self.roi.get_voxel_grid(),
                    orig_spacing=self.roi.spacing,
                    sample_spacing=img_obj.spacing,
                    param_beta=settings.img_interpolate.smoothing_beta,
                    mode="nearest",
                    by_slice=settings.general.by_slice,
                )
            )

        # Register with image
        self.register(img_obj=img_obj, settings=settings)

        # Binarise
        self.binarise_mask()

    def register(self, img_obj: ImageClass, settings, apply_to_self=True):
        """Register roi with image
        Do not apply threshold until after interpolation"""

        if apply_to_self is False:
            roi_copy = self.copy()
            roi_copy.register(img_obj=img_obj, apply_to_self=True)
            return roi_copy

        # Skip if image and/or is missing
        if img_obj is None or self.roi is None:
            return

        # Check whether registration is required
        registration_required = False

        # Mismatch in grid dimension
        if np.any([np.abs(np.array(self.roi.size) - np.array(img_obj.size)) > 0.0]):
            registration_required = True

        # Mismatch in origin
        if np.any([np.abs(self.roi.origin - img_obj.origin) > 0.0]):
            registration_required = True

        # Mismatch in spacing
        if np.any([np.abs(self.roi.spacing - img_obj.spacing) > 0.0]):
            registration_required = True

        if not np.allclose(self.roi.orientation, img_obj.orientation):
            raise ValueError(
                "Cannot register segmentation and image object due to different alignments. "
                "Please use an external programme to transfer segmentation to the image."
            )

        if registration_required:
            # Register roi to image; this transforms the roi grid into
            (
                self.roi.size,
                sample_spacing,
                voxel_grid,
                grid_origin,
            ) = interpolate_to_new_grid(
                orig_dim=self.roi.size,
                orig_spacing=self.roi.spacing,
                orig_vox=self.roi.get_voxel_grid(),
                sample_dim=img_obj.size,
                sample_spacing=img_obj.spacing,
                grid_origin=np.dot(
                    self.roi.m_affine_inv,
                    np.transpose(img_obj.origin - self.roi.origin),
                ),
                order=settings.roi_interpolate.spline_order,
                mode="nearest",
                align_to_center=False,
            )

            # Update origin before spacing, because computing the origin requires the original affine matrix.
            self.roi.origin = self.roi.origin + np.dot(
                self.roi.m_affine, np.transpose(grid_origin)
            )

            # Update spacing and affine matrix.
            self.roi.set_spacing(sample_spacing)

            # Update voxel grid
            self.roi.set_voxel_grid(voxel_grid=voxel_grid)

    def binarise_mask(self):
        if self.roi is None:
            return

        if not self.roi.dtype_name == "bool":
            self.roi.set_voxel_grid(
                voxel_grid=np.around(self.roi.get_voxel_grid(), 6)
                >= np.around(self.incl_threshold, 6)
            )

    def generate_masks(self):
        """ "Generate roi intensity and morphology masks"""

        if self.roi is None:
            self.roi_intensity = None
            self.roi_morphology = None
        else:
            self.roi_intensity = self.roi.copy()
            self.roi_morphology = self.roi.copy()

    def update_roi(self):
        """Update region of interest based on intensity and morphological masks"""

        if (
            self.roi is None
            or self.roi_intensity is None
            or self.roi_morphology is None
        ):
            return

        self.roi.set_voxel_grid(
            voxel_grid=np.logical_or(
                self.roi_intensity.get_voxel_grid(),
                self.roi_morphology.get_voxel_grid(),
            )
        )

    def crop(
        self,
        ind_ext_z=None,
        ind_ext_y=None,
        ind_ext_x=None,
        xy_only=False,
        z_only=False,
    ):
        """ "Resects roi"""

        # Resect masks
        if self.roi is not None:
            self.roi.crop(
                ind_ext_z=ind_ext_z,
                ind_ext_y=ind_ext_y,
                ind_ext_x=ind_ext_x,
                xy_only=xy_only,
                z_only=z_only,
            )

        if self.roi_intensity is not None:
            self.roi_intensity.crop(
                ind_ext_z=ind_ext_z,
                ind_ext_y=ind_ext_y,
                ind_ext_x=ind_ext_x,
                xy_only=xy_only,
                z_only=z_only,
            )

        if self.roi_morphology is not None:
            self.roi_morphology.crop(
                ind_ext_z=ind_ext_z,
                ind_ext_y=ind_ext_y,
                ind_ext_x=ind_ext_x,
                xy_only=xy_only,
                z_only=z_only,
            )

    def crop_to_size(self, center, crop_size, xy_only=False):
        """ "Crops roi to a pre-defined size"""

        # Crop masks to size
        if self.roi is not None:
            self.roi.crop_to_size(center=center, crop_size=crop_size, xy_only=xy_only)
        if self.roi_intensity is not None:
            self.roi_intensity.crop_to_size(
                center=center, crop_size=crop_size, xy_only=xy_only
            )
        if self.roi_morphology is not None:
            self.roi_morphology.crop_to_size(
                center=center, crop_size=crop_size, xy_only=xy_only
            )

    def resegmentise_mask(self, img_obj, by_slice, method, settings):
        # Resegmentation of the roi map based on grey level values

        from skimage.measure import label
        from skimage.morphology import remove_small_holes

        # Skip if required voxel grids are missing
        if (
            img_obj.is_missing
            or self.roi_intensity is None
            or self.roi_morphology is None
        ):
            return

        ################################################################################################################
        # Resegmentation that affects both intensity and morphological maps
        ################################################################################################################

        # Initialise range
        updated_range = np.array([np.nan, np.nan])

        if bool(set(method).intersection(["threshold", "range"])):
            # Filter out voxels with intensity outside prescribed range

            # Local constant
            g_thresh = settings.roi_resegment.g_thresh  # Threshold values

            # Upper threshold
            if not np.isnan(g_thresh[1]):
                updated_range[1] = copy.deepcopy(g_thresh[1])

            # Lower threshold
            if not np.isnan(g_thresh[0]):
                updated_range[0] = copy.deepcopy(g_thresh[0])

            # Set the threshold values as g_range
            self.g_range = g_thresh

        if bool(set(method).intersection(["sigma", "outlier"])):
            # Remove voxels with outlier intensities

            # Local constant
            sigma = settings.roi_resegment.sigma
            img_voxel_grid = img_obj.get_voxel_grid()
            roi_voxel_grid = self.roi_intensity.get_voxel_grid()

            # Check if the voxel grid is not empty
            if np.any(roi_voxel_grid):
                # Calculate mean and standard deviation of intensities in roi
                mean_int = np.mean(img_voxel_grid[roi_voxel_grid])
                sd_int = np.std(img_voxel_grid[roi_voxel_grid])

                if not np.isnan(updated_range[0]):
                    updated_range[0] = np.max(
                        [updated_range[0], mean_int - sigma * sd_int]
                    )
                else:
                    updated_range[0] = mean_int - sigma * sd_int

                if not np.isnan(updated_range[1]):
                    updated_range[1] = np.min(
                        [updated_range[1], mean_int + sigma * sd_int]
                    )
                else:
                    updated_range[1] = mean_int + sigma * sd_int

        if not np.isnan(updated_range[0]) or not np.isnan(updated_range[1]):
            # Update intensity mask
            roi_voxel_grid = self.roi_intensity.get_voxel_grid()

            if not np.isnan(updated_range[0]):
                roi_voxel_grid = np.logical_and(
                    (img_obj.get_voxel_grid() >= updated_range[0]), roi_voxel_grid
                )

            if not np.isnan(updated_range[1]):
                roi_voxel_grid = np.logical_and(
                    (img_obj.get_voxel_grid() <= updated_range[1]), roi_voxel_grid
                )

            # Set roi voxel volume
            self.roi_intensity.set_voxel_grid(voxel_grid=roi_voxel_grid)

        ################################################################################################################
        # Resegmentation that affects only morphological maps
        ################################################################################################################
        if bool(set(method).intersection("close_volume")):
            # Close internal volumes

            # Read minimal volume required
            max_fill_volume = settings.roi_resegment.max_fill_volume

            # Get voxel grid of the roi morphological mask
            roi_voxel_grid = self.roi_morphology.get_voxel_grid()

            # Determine fill volume (in voxels); if max_fill_volume is less than 0.0, fill all holes
            if max_fill_volume < 0.0:
                fill_volume = np.prod(np.array(self.roi_morphology.size)) + 1.0
            else:
                fill_volume = (
                    np.floor(max_fill_volume / np.prod(self.roi_morphology.spacing))
                    + 1.0
                )

            # If the maximum fill volume is smaller than the minimal size of a hole
            if fill_volume < 1.0:
                return None

            # Label all non-roi voxels and get label corresponding to voxels outside of the roi
            non_roi_label = label(
                np.pad(roi_voxel_grid, 1, mode="constant", constant_values=0),
                background=1,
                connectivity=3,
            )
            outside_label = non_roi_label[0, 0, 0]

            # Crop non-roi labels and determine non-roi voxels outside of the mask
            non_roi_label = non_roi_label[1:-1, 1:-1, 1:-1]
            vox_outside = non_roi_label == outside_label

            # Determine mask of voxels which are not internal holes
            vox_not_internal = np.logical_or(roi_voxel_grid, vox_outside)

            # Check if there are any holes, otherwise continue
            if not np.any(~vox_not_internal):
                return None

            if by_slice:
                # 2D approach to filling holes

                for ii in np.arange(0, self.roi_morphology.size[0]):
                    # Skip operations on slides that do not contain voxels in the mask or no holes in the slice
                    if not np.any(roi_voxel_grid[ii, :, :]):
                        continue
                    if not (np.any(~vox_not_internal[ii, :, :])):
                        continue

                    # Fill holes up to fill_volume in voxel number
                    vox_filled = remove_small_holes(
                        vox_not_internal[ii, :, :],
                        min_size=np.int(fill_volume),
                        connectivity=2,
                    )

                    # Update mask by removing outside voxels from the mask
                    roi_voxel_grid[ii, :, :] = np.squeeze(
                        np.logical_and(vox_filled, ~vox_outside[ii, :, :])
                    )
            else:
                # 3D approach to filling holes

                # Fill holes up to fill_volume in voxel number
                vox_filled = remove_small_holes(
                    vox_not_internal, min_size=np.int(fill_volume), connectivity=3
                )

                # Update mask by removing outside voxels from the mask
                roi_voxel_grid = np.logical_and(vox_filled, ~vox_outside)

            # Update voxel grid
            self.roi_morphology.set_voxel_grid(voxel_grid=roi_voxel_grid)

        if bool(set(method).intersection("remove_disconnected")):
            # Remove disconnected voxels

            # Discover prior disconnected volumes from the roi voxel grid
            vox_disconnected = label(
                self.roi.get_voxel_grid(), background=0, connectivity=3
            )
            vox_disconnected_labels = np.unique(vox_disconnected)

            # Set up an empty morphological masks
            upd_vox_mask = np.full(
                shape=self.roi_morphology.size, fill_value=False, dtype=np.bool
            )

            # Get the minimum volume fraction for inclusion as voxels
            min_vol_fract = settings.roi_resegment.min_vol_fract

            # Iterate over disconnected labels
            for curr_volume_label in vox_disconnected_labels:
                # Skip background
                if vox_disconnected_labels == 0:
                    continue

                # Mask only current volume, skip if empty
                curr_mask = np.logical_and(
                    self.roi_morphology.get_voxel_grid(),
                    vox_disconnected == curr_volume_label,
                )
                if not np.any(curr_mask):
                    continue

                # Find fully disconnected voxels groups and count them
                vox_mask = label(curr_mask, background=0, connectivity=3)
                vox_mask_labels, vox_label_count = np.unique(
                    vox_mask, return_counts=True
                )

                # Filter out the background counts
                valid_label_id = np.nonzero(vox_mask_labels)
                vox_mask_labels = vox_mask_labels[valid_label_id]
                vox_label_count = vox_label_count[valid_label_id]

                # Normalise to maximum
                vox_label_count = vox_label_count / np.max(vox_label_count)

                # Select labels fulfilling the minimal size
                vox_mask_labels = vox_mask_labels[vox_label_count >= min_vol_fract]

                for vox_mask_label_id in vox_mask_labels:
                    upd_vox_mask += vox_mask == vox_mask_label_id

                # Update morphological voxel grid
                self.roi_morphology.set_voxel_grid(voxel_grid=upd_vox_mask > 0)

    def is_empty(self):
        """Checks whether the roi or one of its masks is empty"""

        # Original roi object
        if self.roi is not None:
            n_roi_voxels = np.int(np.sum(self.roi.get_voxel_grid()))
            if n_roi_voxels == 0:
                return True

        # Roi intensity mask
        if self.roi_intensity is not None:
            n_roi_int_voxels = np.int(np.sum(self.roi_intensity.get_voxel_grid()))
            if n_roi_int_voxels == 0:
                return True

        # Roi morphological mask
        if self.roi_morphology is not None:
            n_roi_morph_voxels = np.int(np.sum(self.roi_morphology.get_voxel_grid()))
            if n_roi_morph_voxels == 0:
                return True

        # If none of the above rois was empty, return false
        return False

    def rotate(self, angle, img_obj):
        """Rotates roi in the y-x plane"""

        # Register with image prior to rotation
        self.register(img_obj=img_obj)

        # Rotate roi
        self.roi.rotate(angle)

    def dilate(self, by_slice, dist=None, vox_dist=None):
        import scipy.ndimage as ndi

        # Skip if the roi does not exist
        if self.roi is None:
            return

        # Check dtype of the roi voxel grid and binarise if necessary
        if not self.roi.dtype_name == "bool":
            self.binarise_mask()

        # Check if any distance is provided for dilation
        if vox_dist is None and dist is None:
            raise RuntimeError("No dilation distance provided.")

        # Check whether voxel are isometric
        if by_slice:
            spacing = self.roi.spacing[[1, 2]]
        else:
            spacing = self.roi.spacing

        if np.any(spacing - np.max(spacing) != 0.0):
            raise RuntimeWarning(
                "Non-uniform voxel spacing was detected. Roi dilation requires uniform voxel spacing."
            )

        # Derive filter extension and distance
        if dist is not None:
            base_ext: int = np.max([np.floor(dist / np.max(spacing)).astype(np.int), 0])
        else:
            base_ext: int = np.int(vox_dist)
            dist = vox_dist * np.max(spacing)

        # Check if an actual extension is required.
        if base_ext > 0:
            # Create displacement map
            df_base = pd.DataFrame(
                {
                    "x": rep(
                        x=np.arange(-base_ext, base_ext + 1),
                        each=(2 * base_ext + 1) * (2 * base_ext + 1),
                        times=1,
                    ),
                    "y": rep(
                        x=np.arange(-base_ext, base_ext + 1),
                        each=2 * base_ext + 1,
                        times=2 * base_ext + 1,
                    ),
                    "z": rep(
                        x=np.arange(-base_ext, base_ext + 1),
                        each=1,
                        times=(2 * base_ext + 1) * (2 * base_ext + 1),
                    ),
                }
            )

            # Calculate distances for displacement map
            df_base["dist"] = np.sqrt(
                np.sum(
                    np.multiply(
                        df_base.loc[:, ("z", "y", "x")].values, self.roi.spacing
                    )
                    ** 2.0,
                    axis=1,
                )
            )

            # Identify elements in range
            if by_slice:
                df_base["in_range"] = np.logical_and(
                    df_base.dist <= dist, df_base.z == 0
                )
            else:
                df_base["in_range"] = df_base.dist <= dist

            # Update voxel coordinates to start at [0,0,0]
            df_base.loc[:, ["x", "y", "z"]] -= df_base.loc[0, ["x", "y", "z"]]

            # Generate geometric filter structure
            geom_struct = np.zeros(
                shape=(
                    np.max(df_base.z) + 1,
                    np.max(df_base.y) + 1,
                    np.max(df_base.x) + 1,
                ),
                dtype=np.bool,
            )
            geom_struct[
                df_base.z.astype(np.int),
                df_base.y.astype(np.int),
                df_base.x.astype(np.int),
            ] = df_base.in_range

            # Dilate roi mask amd store voxel grid
            self.roi.set_voxel_grid(
                voxel_grid=ndi.binary_dilation(
                    self.roi.get_voxel_grid(), structure=geom_struct, iterations=1
                )
            )

        else:
            print(
                "No dilation: distance %s is too small compared to voxel spacing %s.",
                str(dist),
                str(np.max(spacing)),
            )

    def adapt_volume(self, by_slice, vol_grow_fract=None):
        import scipy.ndimage

        # Skip if the roi does not exist
        if self.roi is None:
            return

        # Check dtype of the roi voxel grid and binarise if necessary
        if not self.roi.dtype_name == "bool":
            print("Converting roi to boolean before dilation.")
            self.binarise_mask()

        # Check if any distance is provided for dilation
        if vol_grow_fract is None:
            raise RuntimeError("No dilation volume fraction was provided.")

        # Check whether voxels are isometric
        if by_slice:
            spacing = self.roi.spacing[[1, 2]]
        else:
            spacing = self.roi.spacing

        if np.any(spacing - np.max(spacing) != 0.0):
            raise RuntimeWarning(
                "Non-uniform voxel spacing was detected. Roi volume adaptation requires uniform voxel spacing."
            )

        # Set geometrical structure
        geom_struct = scipy.ndimage.generate_binary_structure(3, 1)
        if by_slice:
            geom_struct[(0, 2), :, :] = False  # Set structures in different slices to 0

        if not vol_grow_fract == 0.0 and not self.is_empty():
            # Determine original volume
            previous_roi = self.roi.get_voxel_grid()
            orig_volume = np.sum(previous_roi)

            # Iteratively grow or shrink the volume. The loop terminates through break statements
            while True:
                if vol_grow_fract > 0.0:
                    updated_roi = scipy.ndimage.binary_dilation(
                        previous_roi, structure=geom_struct, iterations=1
                    )
                else:
                    updated_roi = scipy.ndimage.binary_erosion(
                        previous_roi, structure=geom_struct, iterations=1
                    )

                new_volume = np.sum(updated_roi)
                if new_volume == 0:
                    break

                if (
                    vol_grow_fract > 0.0
                    and new_volume / orig_volume - 1.0 >= vol_grow_fract
                ):
                    break

                if (
                    vol_grow_fract < 0.0
                    and new_volume / orig_volume - 1.0 <= vol_grow_fract
                ):
                    break

                # Replace previous roi by the updated roi
                previous_roi = updated_roi

            # Randomly add/remove border voxels until desired growth/shrinkage is achieved
            if not new_volume / orig_volume - 1.0 == vol_grow_fract:
                additional_vox = np.abs(
                    int(
                        np.floor(
                            orig_volume * (1.0 + vol_grow_fract) - np.sum(previous_roi)
                        )
                    )
                )
                if additional_vox > 0:
                    border_voxel_ind = np.array(
                        np.where(np.logical_xor(previous_roi, updated_roi))
                    )
                    select_ind = np.random.choice(
                        a=border_voxel_ind.shape[1], size=additional_vox, replace=False
                    )
                    border_voxel_ind = border_voxel_ind[:, select_ind]
                    if vol_grow_fract > 0.0:
                        previous_roi[
                            border_voxel_ind[0, :],
                            border_voxel_ind[1, :],
                            border_voxel_ind[2, :],
                        ] = True
                    else:
                        previous_roi[
                            border_voxel_ind[0, :],
                            border_voxel_ind[1, :],
                            border_voxel_ind[2, :],
                        ] = False

            # Set the new roi
            self.roi.set_voxel_grid(voxel_grid=previous_roi)

    def erode(self, by_slice, eroded_vol_fract=0.8, dist=None, vox_dist=None):
        """ " Erosion of the roi segment"""

        import scipy.ndimage as ndi

        # Skip if the roi does not exist
        if self.roi is None:
            return

        # Check whether voxels are booleans
        if not self.roi.dtype_name == "bool":
            print("Converting roi to boolean before erosion.")
            self.binarise_mask()

        # Check if any distance is provided for dilation
        if vox_dist is None and dist is None:
            raise RuntimeError("No erosion distance provided.")

        # Check whether voxel are isometric
        if by_slice:
            spacing = self.roi.spacing[[1, 2]]
        else:
            spacing = self.roi.spacing

        if np.any(spacing - np.max(spacing) != 0.0):
            raise RuntimeWarning(
                "Non-uniform voxel spacing was detected. Roi erosion requires uniform voxel spacing."
            )

        # Set geometrical structure
        geom_struct = ndi.generate_binary_structure(3, 1)
        if by_slice:
            geom_struct[(0, 2), :, :] = False  # Set structures in different slices to 0

        # Set number of erosion steps
        if vox_dist is None:
            erode_steps = np.max(
                [np.round(np.abs(dist) / np.max(spacing)).astype(np.int), 0]
            )
        else:
            erode_steps = np.abs(vox_dist.astype(np.int))
            dist = vox_dist * np.max(spacing)

        if erode_steps > 0:
            # Determine initial volume
            voxels_prev = self.roi.get_voxel_grid()
            vol_init = np.sum(voxels_prev)

            # Iterate over erosion steps
            for step in np.arange(0, erode_steps):
                # Perform erosion
                voxels_upd = ndi.binary_erosion(
                    voxels_prev, structure=geom_struct, iterations=1
                )

                # Calculate volume of the eroded volume
                vol_curr = np.sum(voxels_upd)

                # Stop erosion if the volume shrinks below 80 percent of the original volume due to erosion and return
                # return voxels from the previous erosion step.
                if vol_curr * 1.0 / vol_init < eroded_vol_fract:
                    voxels_upd = voxels_prev
                    break
                else:
                    voxels_prev = voxels_upd

            # Set updated voxels
            self.roi.set_voxel_grid(voxel_grid=voxels_upd)
        else:
            print(
                "No erosion: distance %s is too small compared to voxel spacing %s.",
                str(dist),
                str(np.max(spacing)),
            )

    def decode_voxel_grid(self):
        """Converts run length encoded grids to conventional volumes"""

        # Decode main ROI object
        if self.roi is not None:
            self.roi.decode_voxel_grid()

        # Decode intensity and morphological masks
        if self.roi_intensity is not None:
            self.roi_intensity.decode_voxel_grid()
        if self.roi_morphology is not None:
            self.roi_morphology.decode_voxel_grid()

    def as_pandas_dataframe(
        self,
        img_obj,
        intensity_mask=False,
        morphology_mask=False,
        distance_map=False,
        by_slice=False,
    ):
        """Converts the image and roi voxel grids to a pandas dataframe for further processing"""

        # Return None if the image and/or ROI are missing
        if img_obj.is_missing or self.roi is None:
            return None

        # Check if the masks exist and assign if not
        if intensity_mask and self.roi_intensity is None:
            self.roi_intensity = self.roi.copy()
        if (morphology_mask or distance_map) and self.roi_morphology is None:
            self.roi_morphology = self.roi.copy()

        # Create table from test object
        img_dims = img_obj.size
        index_id = np.arange(start=0, stop=np.prod(img_dims))
        coords = np.unravel_index(indices=index_id, dims=img_dims)
        df_img = pd.DataFrame(
            {
                "index_id": index_id,
                "g": np.ravel(img_obj.get_voxel_grid()),
                "x": coords[2],
                "y": coords[1],
                "z": coords[0],
            }
        )

        if intensity_mask:
            df_img["roi_int_mask"] = np.ravel(
                self.roi_intensity.get_voxel_grid()
            ).astype(np.bool)
        if morphology_mask:
            df_img["roi_morph_mask"] = np.ravel(
                self.roi_morphology.get_voxel_grid()
            ).astype(np.bool)
        if distance_map:
            # Calculate distance by sequential border erosion
            from scipy.ndimage import binary_erosion, generate_binary_structure

            # Set up distance map and morphological voxel grid
            dist_map = np.zeros(img_dims)
            morph_voxel_grid = self.roi_morphology.get_voxel_grid()

            if by_slice:
                # Distances are determined in 2D
                binary_struct = generate_binary_structure(rank=2, connectivity=1)

                # Iterate over slices
                for ii in np.arange(0, img_dims[0]):
                    # Calculate distance by sequential border erosion
                    roi_eroded = np.squeeze(morph_voxel_grid[ii, :, :])

                    # Iterate distance from border
                    while np.sum(roi_eroded) > 0:
                        roi_eroded = binary_erosion(roi_eroded, structure=binary_struct)
                        dist_map[ii, :, :] += roi_eroded * 1

            else:
                # Distances are determined in 3D
                binary_struct = generate_binary_structure(rank=3, connectivity=1)

                # Copy of roi morphology mask
                roi_eroded = copy.deepcopy(morph_voxel_grid)

                # Incrementally erode the morphological mask
                while np.sum(roi_eroded) > 0:
                    roi_eroded = binary_erosion(roi_eroded, structure=binary_struct)
                    dist_map += roi_eroded * 1

            # Update distance from border, as minimum distance is 1
            dist_map[morph_voxel_grid] += 1

            # Add distance map to table
            df_img["border_distance"] = np.ravel(dist_map).astype(np.int32)

        return df_img

    def compute_diagnostic_features(self, img_obj, append_str=""):
        """Creates diagnostic features for the ROI"""

        # Set feature names
        feat_names = [
            "int_map_dim_x",
            "int_map_dim_y",
            "int_map_dim_z",
            "int_bb_dim_x",
            "int_bb_dim_y",
            "int_bb_dim_z",
            "int_vox_dim_x",
            "int_vox_dim_y",
            "int_vox_dim_z",
            "int_vox_count",
            "int_mean_int",
            "int_min_int",
            "int_max_int",
            "mrp_map_dim_x",
            "mrp_map_dim_y",
            "mrp_map_dim_z",
            "mrp_bb_dim_x",
            "mrp_bb_dim_y",
            "mrp_bb_dim_z",
            "mrp_vox_dim_x",
            "mrp_vox_dim_y",
            "mrp_vox_dim_z",
            "mrp_vox_count",
            "mrp_mean_int",
            "mrp_min_int",
            "mrp_max_int",
        ]

        # Create pandas dataframe with one row and feature columns
        df = pd.DataFrame(np.full(shape=(1, len(feat_names)), fill_value=np.nan))
        df.columns = feat_names

        # Skip further analysis if the image and/or roi are missing
        if img_obj.is_missing or self.roi is None:
            return df

        # Register with image on function call
        roi_copy = self.register(img_obj, apply_to_self=False)

        # Binarise (if required)
        roi_copy.binarise_mask()

        # Make copies of intensity and morphological masks (if required)
        if roi_copy.roi_intensity is None:
            roi_copy.roi_intensity = roi_copy.roi
        if roi_copy.roi_morphology is None:
            roi_copy.roi_morphology = roi_copy.roi

        # Get image and roi voxel grids
        img_voxel_grid = img_obj.get_voxel_grid()
        int_voxel_grid = roi_copy.roi_intensity.get_voxel_grid()
        mrp_voxel_grid = roi_copy.roi_morphology.get_voxel_grid()

        # Compute bounding boxes
        int_bounding_box_dim = np.squeeze(
            np.diff(roi_copy.get_bounding_box(roi_voxel_grid=int_voxel_grid), axis=0)
            + 1
        )
        mrp_bounding_box_dim = np.squeeze(
            np.diff(roi_copy.get_bounding_box(roi_voxel_grid=mrp_voxel_grid), axis=0)
            + 1
        )

        # Set intensity mask features
        df["int_map_dim_x"] = roi_copy.roi_intensity.size[2]
        df["int_map_dim_y"] = roi_copy.roi_intensity.size[1]
        df["int_map_dim_z"] = roi_copy.roi_intensity.size[0]
        df["int_bb_dim_x"] = int_bounding_box_dim[2]
        df["int_bb_dim_y"] = int_bounding_box_dim[1]
        df["int_bb_dim_z"] = int_bounding_box_dim[0]
        df["int_vox_dim_x"] = roi_copy.roi_intensity.spacing[2]
        df["int_vox_dim_y"] = roi_copy.roi_intensity.spacing[1]
        df["int_vox_dim_z"] = roi_copy.roi_intensity.spacing[0]
        df["int_vox_count"] = np.sum(int_voxel_grid)
        df["int_mean_int"] = np.mean(img_voxel_grid[int_voxel_grid])
        df["int_min_int"] = np.min(img_voxel_grid[int_voxel_grid])
        df["int_max_int"] = np.max(img_voxel_grid[int_voxel_grid])

        # Set morphological mask features
        df["mrp_map_dim_x"] = roi_copy.roi_morphology.size[2]
        df["mrp_map_dim_y"] = roi_copy.roi_morphology.size[1]
        df["mrp_map_dim_z"] = roi_copy.roi_morphology.size[0]
        df["mrp_bb_dim_x"] = mrp_bounding_box_dim[2]
        df["mrp_bb_dim_y"] = mrp_bounding_box_dim[1]
        df["mrp_bb_dim_z"] = mrp_bounding_box_dim[0]
        df["mrp_vox_dim_x"] = roi_copy.roi_morphology.spacing[2]
        df["mrp_vox_dim_y"] = roi_copy.roi_morphology.spacing[1]
        df["mrp_vox_dim_z"] = roi_copy.roi_morphology.spacing[0]
        df["mrp_vox_count"] = np.sum(mrp_voxel_grid)
        df["mrp_mean_int"] = np.mean(img_voxel_grid[mrp_voxel_grid])
        df["mrp_min_int"] = np.min(img_voxel_grid[mrp_voxel_grid])
        df["mrp_max_int"] = np.max(img_voxel_grid[mrp_voxel_grid])

        # Update column names
        df.columns = [
            "_".join(["diag", feature, append_str]).strip("_") for feature in df.columns
        ]

        del roi_copy

        self.diagnostic_list += [df]

    def get_bounding_box(self, roi_voxel_grid):
        # Calculates coordinates of ROI bounding box
        z_ind, y_ind, x_ind = np.where(roi_voxel_grid)
        max_ind = np.array((np.max(z_ind), np.max(y_ind), np.max(x_ind)))
        min_ind = np.array((np.min(z_ind), np.min(y_ind), np.min(x_ind)))
        del z_ind, y_ind, x_ind

        return min_ind, max_ind

    def get_center_slice(self):
        """Identify location of the central slice in the roi"""

        # Return a NaN if no roi is present
        if self.roi is None:
            return np.nan

        # Determine indices of voxels included in the roi
        z_ind, y_ind, x_ind = np.where(self.roi.get_voxel_grid())
        z_center = (np.max(z_ind) + np.min(z_ind)) // 2

        return z_center

    def get_all_slices(self):
        """Identify location of all slices in the roi"""

        # Return NaN in case the roi is missing
        if self.roi is None:
            return np.array([np.nan])

        z_ind, y_ind, x_ind = np.where(self.roi.get_voxel_grid())

        return np.unique(z_ind)

    def export(self, img_obj, file_path):
        """
        Export roi to file
        :param img_obj:
        :param file_path:
        :return:
        """

        roi_str_components = [img_obj.get_export_descriptor()]
        roi_str_components += [self.get_export_descriptor()]

        # Write morphological and intensity roi
        if self.roi_morphology is not None and self.roi_intensity is not None:
            self.roi_morphology.write(
                file_path=file_path,
                file_name="_".join(roi_str_components + ["morph.nii.gz"]),
            )
            self.roi_intensity.write(
                file_path=file_path,
                file_name="_".join(roi_str_components + ["int.nii.gz"]),
            )

        elif self.roi is not None:
            return self.roi.write(
                file_path=file_path,
                file_name="_".join(roi_str_components + ["roi.nii"]),
            )

        else:
            return

    def get_export_descriptor(self):
        """
        Generates an export string for identifying a file
        :return: export string
        """
        descr_list = []

        if self.adapt_size != 0.0:
            # Volume adaptation
            descr_list += ["vol", str(self.adapt_size)]
        if self.svx_randomisation_id != -1:
            # Contour randomisation
            descr_list += ["svx", str(self.svx_randomisation_id)]

        descr_list += [self.name]

        return "_".join(descr_list)

    def get_slices(self, slice_number=None):
        # Extract roi objects for each slice

        roi_obj_list = []

        # Create a copy of the current object
        base_roi_obj = self.copy()

        # Remove attributes that need to be set
        base_roi_obj.roi = None
        base_roi_obj.roi_intensity = None
        base_roi_obj.roi_morphology = None

        if slice_number is None:
            # Extract mask for each slice.  Copy the base roi object.

            if self.roi is not None:
                roi_slices = self.roi.get_slices()

            if self.roi_intensity is not None:
                roi_int_slices = self.roi_intensity.get_slices()
            else:
                roi_int_slices = None

            if self.roi_morphology is not None:
                roi_morph_slices = self.roi_morphology.get_slices()
            else:
                roi_morph_slices = None

            # Add masks to a roi object for each slice
            for ii in np.arange(self.roi.size[0]):
                slice_roi_obj = copy.deepcopy(base_roi_obj)

                if self.roi is not None:
                    slice_roi_obj.roi = roi_slices[ii]
                if self.roi_intensity is not None:
                    slice_roi_obj.roi_intensity = roi_int_slices[ii]
                if self.roi_morphology is not None:
                    slice_roi_obj.roi_morphology = roi_morph_slices[ii]

                # Add to list
                roi_obj_list += [slice_roi_obj]
        else:
            # Extract a single slice. Copy the base roi object.
            slice_roi_obj = copy.deepcopy(base_roi_obj)

            # Add the mask for the requested slice
            if self.roi is not None:
                slice_roi_obj.roi = self.roi.get_slices(slice_number=slice_number)
            if self.roi_intensity is not None:
                slice_roi_obj.roi_intensity = self.roi_intensity.get_slices(
                    slice_number=slice_number
                )
            if self.roi_morphology is not None:
                slice_roi_obj.roi_morphology = self.roi_morphology.get_slices(
                    slice_number=slice_number
                )

            # Add to list
            roi_obj_list += [slice_roi_obj]

        return roi_obj_list

    def write_dicom(self, file_path, file_name="RS.dcm"):
        import os

        if self.metadata is None:
            return None

        # Check if the write folder exists
        if not os.path.isdir(file_path):
            if os.path.isfile(file_path):
                # Check if the write folder is a file.
                raise IOError(
                    f"{file_path} is an existing file, not a directory. No DICOM images were exported."
                )
            else:
                os.makedirs(file_path, exist_ok=True)

        self.metadata.save_as(
            filename=os.path.join(file_path, file_name), write_like_original=False
        )

    def get_metadata(self, tag, tag_type, default=None):
        # Do not attempt to read the metadata if no metadata is present.
        if self.metadata is None:
            return

        return get_pydicom_meta_tag(
            dcm_seq=self.metadata, tag=tag, tag_type=tag_type, default=default
        )

    def set_metadata(self, tag, value, force_vr=None):
        # Do not update the metadata if no metadata is present.
        if self.metadata is None:
            return None

        set_pydicom_meta_tag(
            dcm_seq=self.metadata, tag=tag, value=value, force_vr=force_vr
        )

    def has_metadata(self, tag):
        if self.metadata is None:
            return None

        else:
            return get_pydicom_meta_tag(dcm_seq=self.metadata, tag=tag, test_tag=True)

    def rename(self, new):
        if self.metadata is not None:
            # Obtain the old name
            old = self.name

            if not self.has_metadata(tag=(0x3006, 0x0020)):
                raise ValueError(
                    "The DICOM metaheader does not contain a Structure Set ROI sequence."
                )

            # Iterate over roi elements in the roi sequence
            for ii, roi_element in enumerate(self.metadata[0x3006, 0x0020]):
                # Find ROI name that matches the old name
                if (
                    get_pydicom_meta_tag(
                        dcm_seq=roi_element, tag=(0x3006, 0x0026), tag_type="str"
                    )
                    == old
                ):
                    set_pydicom_meta_tag(
                        dcm_seq=roi_element, tag=(0x3006, 0x0026), value=new
                    )

            # Assign a new name
            self.name = new
        else:
            # Assign a new name
            self.name = new

as_pandas_dataframe(img_obj, intensity_mask=False, morphology_mask=False, distance_map=False, by_slice=False)

Converts the image and roi voxel grids to a pandas dataframe for further processing

Source code in src/luna/radiology/mirp/roiClass.py

def as_pandas_dataframe(
    self,
    img_obj,
    intensity_mask=False,
    morphology_mask=False,
    distance_map=False,
    by_slice=False,
):
    """Converts the image and roi voxel grids to a pandas dataframe for further processing"""

    # Return None if the image and/or ROI are missing
    if img_obj.is_missing or self.roi is None:
        return None

    # Check if the masks exist and assign if not
    if intensity_mask and self.roi_intensity is None:
        self.roi_intensity = self.roi.copy()
    if (morphology_mask or distance_map) and self.roi_morphology is None:
        self.roi_morphology = self.roi.copy()

    # Create table from test object
    img_dims = img_obj.size
    index_id = np.arange(start=0, stop=np.prod(img_dims))
    coords = np.unravel_index(indices=index_id, dims=img_dims)
    df_img = pd.DataFrame(
        {
            "index_id": index_id,
            "g": np.ravel(img_obj.get_voxel_grid()),
            "x": coords[2],
            "y": coords[1],
            "z": coords[0],
        }
    )

    if intensity_mask:
        df_img["roi_int_mask"] = np.ravel(
            self.roi_intensity.get_voxel_grid()
        ).astype(np.bool)
    if morphology_mask:
        df_img["roi_morph_mask"] = np.ravel(
            self.roi_morphology.get_voxel_grid()
        ).astype(np.bool)
    if distance_map:
        # Calculate distance by sequential border erosion
        from scipy.ndimage import binary_erosion, generate_binary_structure

        # Set up distance map and morphological voxel grid
        dist_map = np.zeros(img_dims)
        morph_voxel_grid = self.roi_morphology.get_voxel_grid()

        if by_slice:
            # Distances are determined in 2D
            binary_struct = generate_binary_structure(rank=2, connectivity=1)

            # Iterate over slices
            for ii in np.arange(0, img_dims[0]):
                # Calculate distance by sequential border erosion
                roi_eroded = np.squeeze(morph_voxel_grid[ii, :, :])

                # Iterate distance from border
                while np.sum(roi_eroded) > 0:
                    roi_eroded = binary_erosion(roi_eroded, structure=binary_struct)
                    dist_map[ii, :, :] += roi_eroded * 1

        else:
            # Distances are determined in 3D
            binary_struct = generate_binary_structure(rank=3, connectivity=1)

            # Copy of roi morphology mask
            roi_eroded = copy.deepcopy(morph_voxel_grid)

            # Incrementally erode the morphological mask
            while np.sum(roi_eroded) > 0:
                roi_eroded = binary_erosion(roi_eroded, structure=binary_struct)
                dist_map += roi_eroded * 1

        # Update distance from border, as minimum distance is 1
        dist_map[morph_voxel_grid] += 1

        # Add distance map to table
        df_img["border_distance"] = np.ravel(dist_map).astype(np.int32)

    return df_img

compute_diagnostic_features(img_obj, append_str='')

Creates diagnostic features for the ROI

Source code in src/luna/radiology/mirp/roiClass.py

def compute_diagnostic_features(self, img_obj, append_str=""):
    """Creates diagnostic features for the ROI"""

    # Set feature names
    feat_names = [
        "int_map_dim_x",
        "int_map_dim_y",
        "int_map_dim_z",
        "int_bb_dim_x",
        "int_bb_dim_y",
        "int_bb_dim_z",
        "int_vox_dim_x",
        "int_vox_dim_y",
        "int_vox_dim_z",
        "int_vox_count",
        "int_mean_int",
        "int_min_int",
        "int_max_int",
        "mrp_map_dim_x",
        "mrp_map_dim_y",
        "mrp_map_dim_z",
        "mrp_bb_dim_x",
        "mrp_bb_dim_y",
        "mrp_bb_dim_z",
        "mrp_vox_dim_x",
        "mrp_vox_dim_y",
        "mrp_vox_dim_z",
        "mrp_vox_count",
        "mrp_mean_int",
        "mrp_min_int",
        "mrp_max_int",
    ]

    # Create pandas dataframe with one row and feature columns
    df = pd.DataFrame(np.full(shape=(1, len(feat_names)), fill_value=np.nan))
    df.columns = feat_names

    # Skip further analysis if the image and/or roi are missing
    if img_obj.is_missing or self.roi is None:
        return df

    # Register with image on function call
    roi_copy = self.register(img_obj, apply_to_self=False)

    # Binarise (if required)
    roi_copy.binarise_mask()

    # Make copies of intensity and morphological masks (if required)
    if roi_copy.roi_intensity is None:
        roi_copy.roi_intensity = roi_copy.roi
    if roi_copy.roi_morphology is None:
        roi_copy.roi_morphology = roi_copy.roi

    # Get image and roi voxel grids
    img_voxel_grid = img_obj.get_voxel_grid()
    int_voxel_grid = roi_copy.roi_intensity.get_voxel_grid()
    mrp_voxel_grid = roi_copy.roi_morphology.get_voxel_grid()

    # Compute bounding boxes
    int_bounding_box_dim = np.squeeze(
        np.diff(roi_copy.get_bounding_box(roi_voxel_grid=int_voxel_grid), axis=0)
        + 1
    )
    mrp_bounding_box_dim = np.squeeze(
        np.diff(roi_copy.get_bounding_box(roi_voxel_grid=mrp_voxel_grid), axis=0)
        + 1
    )

    # Set intensity mask features
    df["int_map_dim_x"] = roi_copy.roi_intensity.size[2]
    df["int_map_dim_y"] = roi_copy.roi_intensity.size[1]
    df["int_map_dim_z"] = roi_copy.roi_intensity.size[0]
    df["int_bb_dim_x"] = int_bounding_box_dim[2]
    df["int_bb_dim_y"] = int_bounding_box_dim[1]
    df["int_bb_dim_z"] = int_bounding_box_dim[0]
    df["int_vox_dim_x"] = roi_copy.roi_intensity.spacing[2]
    df["int_vox_dim_y"] = roi_copy.roi_intensity.spacing[1]
    df["int_vox_dim_z"] = roi_copy.roi_intensity.spacing[0]
    df["int_vox_count"] = np.sum(int_voxel_grid)
    df["int_mean_int"] = np.mean(img_voxel_grid[int_voxel_grid])
    df["int_min_int"] = np.min(img_voxel_grid[int_voxel_grid])
    df["int_max_int"] = np.max(img_voxel_grid[int_voxel_grid])

    # Set morphological mask features
    df["mrp_map_dim_x"] = roi_copy.roi_morphology.size[2]
    df["mrp_map_dim_y"] = roi_copy.roi_morphology.size[1]
    df["mrp_map_dim_z"] = roi_copy.roi_morphology.size[0]
    df["mrp_bb_dim_x"] = mrp_bounding_box_dim[2]
    df["mrp_bb_dim_y"] = mrp_bounding_box_dim[1]
    df["mrp_bb_dim_z"] = mrp_bounding_box_dim[0]
    df["mrp_vox_dim_x"] = roi_copy.roi_morphology.spacing[2]
    df["mrp_vox_dim_y"] = roi_copy.roi_morphology.spacing[1]
    df["mrp_vox_dim_z"] = roi_copy.roi_morphology.spacing[0]
    df["mrp_vox_count"] = np.sum(mrp_voxel_grid)
    df["mrp_mean_int"] = np.mean(img_voxel_grid[mrp_voxel_grid])
    df["mrp_min_int"] = np.min(img_voxel_grid[mrp_voxel_grid])
    df["mrp_max_int"] = np.max(img_voxel_grid[mrp_voxel_grid])

    # Update column names
    df.columns = [
        "_".join(["diag", feature, append_str]).strip("_") for feature in df.columns
    ]

    del roi_copy

    self.diagnostic_list += [df]

crop(ind_ext_z=None, ind_ext_y=None, ind_ext_x=None, xy_only=False, z_only=False)

"Resects roi

Source code in src/luna/radiology/mirp/roiClass.py

def crop(
    self,
    ind_ext_z=None,
    ind_ext_y=None,
    ind_ext_x=None,
    xy_only=False,
    z_only=False,
):
    """ "Resects roi"""

    # Resect masks
    if self.roi is not None:
        self.roi.crop(
            ind_ext_z=ind_ext_z,
            ind_ext_y=ind_ext_y,
            ind_ext_x=ind_ext_x,
            xy_only=xy_only,
            z_only=z_only,
        )

    if self.roi_intensity is not None:
        self.roi_intensity.crop(
            ind_ext_z=ind_ext_z,
            ind_ext_y=ind_ext_y,
            ind_ext_x=ind_ext_x,
            xy_only=xy_only,
            z_only=z_only,
        )

    if self.roi_morphology is not None:
        self.roi_morphology.crop(
            ind_ext_z=ind_ext_z,
            ind_ext_y=ind_ext_y,
            ind_ext_x=ind_ext_x,
            xy_only=xy_only,
            z_only=z_only,
        )

crop_to_size(center, crop_size, xy_only=False)

"Crops roi to a pre-defined size

Source code in src/luna/radiology/mirp/roiClass.py

def crop_to_size(self, center, crop_size, xy_only=False):
    """ "Crops roi to a pre-defined size"""

    # Crop masks to size
    if self.roi is not None:
        self.roi.crop_to_size(center=center, crop_size=crop_size, xy_only=xy_only)
    if self.roi_intensity is not None:
        self.roi_intensity.crop_to_size(
            center=center, crop_size=crop_size, xy_only=xy_only
        )
    if self.roi_morphology is not None:
        self.roi_morphology.crop_to_size(
            center=center, crop_size=crop_size, xy_only=xy_only
        )

decimate(by_slice)

Decimates the roi :param by_slice: boolean, 2D (True) or 3D (False) :return:

Source code in src/luna/radiology/mirp/roiClass.py

def decimate(self, by_slice):
    """
    Decimates the roi
    :param by_slice: boolean, 2D (True) or 3D (False)
    :return:
    """

    # Resect masks
    if self.roi is not None:
        self.roi.decimate(by_slice=by_slice)
    if self.roi_intensity is not None:
        self.roi_intensity.decimate(by_slice=by_slice)
    if self.roi_morphology is not None:
        self.roi_morphology.decimate(by_slice=by_slice)

decode_voxel_grid()

Converts run length encoded grids to conventional volumes

Source code in src/luna/radiology/mirp/roiClass.py

def decode_voxel_grid(self):
    """Converts run length encoded grids to conventional volumes"""

    # Decode main ROI object
    if self.roi is not None:
        self.roi.decode_voxel_grid()

    # Decode intensity and morphological masks
    if self.roi_intensity is not None:
        self.roi_intensity.decode_voxel_grid()
    if self.roi_morphology is not None:
        self.roi_morphology.decode_voxel_grid()

erode(by_slice, eroded_vol_fract=0.8, dist=None, vox_dist=None)

" Erosion of the roi segment

Source code in src/luna/radiology/mirp/roiClass.py

def erode(self, by_slice, eroded_vol_fract=0.8, dist=None, vox_dist=None):
    """ " Erosion of the roi segment"""

    import scipy.ndimage as ndi

    # Skip if the roi does not exist
    if self.roi is None:
        return

    # Check whether voxels are booleans
    if not self.roi.dtype_name == "bool":
        print("Converting roi to boolean before erosion.")
        self.binarise_mask()

    # Check if any distance is provided for dilation
    if vox_dist is None and dist is None:
        raise RuntimeError("No erosion distance provided.")

    # Check whether voxel are isometric
    if by_slice:
        spacing = self.roi.spacing[[1, 2]]
    else:
        spacing = self.roi.spacing

    if np.any(spacing - np.max(spacing) != 0.0):
        raise RuntimeWarning(
            "Non-uniform voxel spacing was detected. Roi erosion requires uniform voxel spacing."
        )

    # Set geometrical structure
    geom_struct = ndi.generate_binary_structure(3, 1)
    if by_slice:
        geom_struct[(0, 2), :, :] = False  # Set structures in different slices to 0

    # Set number of erosion steps
    if vox_dist is None:
        erode_steps = np.max(
            [np.round(np.abs(dist) / np.max(spacing)).astype(np.int), 0]
        )
    else:
        erode_steps = np.abs(vox_dist.astype(np.int))
        dist = vox_dist * np.max(spacing)

    if erode_steps > 0:
        # Determine initial volume
        voxels_prev = self.roi.get_voxel_grid()
        vol_init = np.sum(voxels_prev)

        # Iterate over erosion steps
        for step in np.arange(0, erode_steps):
            # Perform erosion
            voxels_upd = ndi.binary_erosion(
                voxels_prev, structure=geom_struct, iterations=1
            )

            # Calculate volume of the eroded volume
            vol_curr = np.sum(voxels_upd)

            # Stop erosion if the volume shrinks below 80 percent of the original volume due to erosion and return
            # return voxels from the previous erosion step.
            if vol_curr * 1.0 / vol_init < eroded_vol_fract:
                voxels_upd = voxels_prev
                break
            else:
                voxels_prev = voxels_upd

        # Set updated voxels
        self.roi.set_voxel_grid(voxel_grid=voxels_upd)
    else:
        print(
            "No erosion: distance %s is too small compared to voxel spacing %s.",
            str(dist),
            str(np.max(spacing)),
        )

export(img_obj, file_path)

Export roi to file :param img_obj: :param file_path: :return:

Source code in src/luna/radiology/mirp/roiClass.py

def export(self, img_obj, file_path):
    """
    Export roi to file
    :param img_obj:
    :param file_path:
    :return:
    """

    roi_str_components = [img_obj.get_export_descriptor()]
    roi_str_components += [self.get_export_descriptor()]

    # Write morphological and intensity roi
    if self.roi_morphology is not None and self.roi_intensity is not None:
        self.roi_morphology.write(
            file_path=file_path,
            file_name="_".join(roi_str_components + ["morph.nii.gz"]),
        )
        self.roi_intensity.write(
            file_path=file_path,
            file_name="_".join(roi_str_components + ["int.nii.gz"]),
        )

    elif self.roi is not None:
        return self.roi.write(
            file_path=file_path,
            file_name="_".join(roi_str_components + ["roi.nii"]),
        )

    else:
        return

generate_masks()

"Generate roi intensity and morphology masks

Source code in src/luna/radiology/mirp/roiClass.py

def generate_masks(self):
    """ "Generate roi intensity and morphology masks"""

    if self.roi is None:
        self.roi_intensity = None
        self.roi_morphology = None
    else:
        self.roi_intensity = self.roi.copy()
        self.roi_morphology = self.roi.copy()

get_all_slices()

Identify location of all slices in the roi

Source code in src/luna/radiology/mirp/roiClass.py

def get_all_slices(self):
    """Identify location of all slices in the roi"""

    # Return NaN in case the roi is missing
    if self.roi is None:
        return np.array([np.nan])

    z_ind, y_ind, x_ind = np.where(self.roi.get_voxel_grid())

    return np.unique(z_ind)

get_center_slice()

Identify location of the central slice in the roi

Source code in src/luna/radiology/mirp/roiClass.py

def get_center_slice(self):
    """Identify location of the central slice in the roi"""

    # Return a NaN if no roi is present
    if self.roi is None:
        return np.nan

    # Determine indices of voxels included in the roi
    z_ind, y_ind, x_ind = np.where(self.roi.get_voxel_grid())
    z_center = (np.max(z_ind) + np.min(z_ind)) // 2

    return z_center

get_export_descriptor()

Generates an export string for identifying a file :return: export string

Source code in src/luna/radiology/mirp/roiClass.py

def get_export_descriptor(self):
    """
    Generates an export string for identifying a file
    :return: export string
    """
    descr_list = []

    if self.adapt_size != 0.0:
        # Volume adaptation
        descr_list += ["vol", str(self.adapt_size)]
    if self.svx_randomisation_id != -1:
        # Contour randomisation
        descr_list += ["svx", str(self.svx_randomisation_id)]

    descr_list += [self.name]

    return "_".join(descr_list)

is_empty()

Checks whether the roi or one of its masks is empty

Source code in src/luna/radiology/mirp/roiClass.py

def is_empty(self):
    """Checks whether the roi or one of its masks is empty"""

    # Original roi object
    if self.roi is not None:
        n_roi_voxels = np.int(np.sum(self.roi.get_voxel_grid()))
        if n_roi_voxels == 0:
            return True

    # Roi intensity mask
    if self.roi_intensity is not None:
        n_roi_int_voxels = np.int(np.sum(self.roi_intensity.get_voxel_grid()))
        if n_roi_int_voxels == 0:
            return True

    # Roi morphological mask
    if self.roi_morphology is not None:
        n_roi_morph_voxels = np.int(np.sum(self.roi_morphology.get_voxel_grid()))
        if n_roi_morph_voxels == 0:
            return True

    # If none of the above rois was empty, return false
    return False

register(img_obj, settings, apply_to_self=True)

Register roi with image Do not apply threshold until after interpolation

Source code in src/luna/radiology/mirp/roiClass.py

def register(self, img_obj: ImageClass, settings, apply_to_self=True):
    """Register roi with image
    Do not apply threshold until after interpolation"""

    if apply_to_self is False:
        roi_copy = self.copy()
        roi_copy.register(img_obj=img_obj, apply_to_self=True)
        return roi_copy

    # Skip if image and/or is missing
    if img_obj is None or self.roi is None:
        return

    # Check whether registration is required
    registration_required = False

    # Mismatch in grid dimension
    if np.any([np.abs(np.array(self.roi.size) - np.array(img_obj.size)) > 0.0]):
        registration_required = True

    # Mismatch in origin
    if np.any([np.abs(self.roi.origin - img_obj.origin) > 0.0]):
        registration_required = True

    # Mismatch in spacing
    if np.any([np.abs(self.roi.spacing - img_obj.spacing) > 0.0]):
        registration_required = True

    if not np.allclose(self.roi.orientation, img_obj.orientation):
        raise ValueError(
            "Cannot register segmentation and image object due to different alignments. "
            "Please use an external programme to transfer segmentation to the image."
        )

    if registration_required:
        # Register roi to image; this transforms the roi grid into
        (
            self.roi.size,
            sample_spacing,
            voxel_grid,
            grid_origin,
        ) = interpolate_to_new_grid(
            orig_dim=self.roi.size,
            orig_spacing=self.roi.spacing,
            orig_vox=self.roi.get_voxel_grid(),
            sample_dim=img_obj.size,
            sample_spacing=img_obj.spacing,
            grid_origin=np.dot(
                self.roi.m_affine_inv,
                np.transpose(img_obj.origin - self.roi.origin),
            ),
            order=settings.roi_interpolate.spline_order,
            mode="nearest",
            align_to_center=False,
        )

        # Update origin before spacing, because computing the origin requires the original affine matrix.
        self.roi.origin = self.roi.origin + np.dot(
            self.roi.m_affine, np.transpose(grid_origin)
        )

        # Update spacing and affine matrix.
        self.roi.set_spacing(sample_spacing)

        # Update voxel grid
        self.roi.set_voxel_grid(voxel_grid=voxel_grid)

rotate(angle, img_obj)

Rotates roi in the y-x plane

Source code in src/luna/radiology/mirp/roiClass.py

def rotate(self, angle, img_obj):
    """Rotates roi in the y-x plane"""

    # Register with image prior to rotation
    self.register(img_obj=img_obj)

    # Rotate roi
    self.roi.rotate(angle)

update_roi()

Update region of interest based on intensity and morphological masks

Source code in src/luna/radiology/mirp/roiClass.py

def update_roi(self):
    """Update region of interest based on intensity and morphological masks"""

    if (
        self.roi is None
        or self.roi_intensity is None
        or self.roi_morphology is None
    ):
        return

    self.roi.set_voxel_grid(
        voxel_grid=np.logical_or(
            self.roi_intensity.get_voxel_grid(),
            self.roi_morphology.get_voxel_grid(),
        )
    )

rep(x, each=1, times=1)

"This functions replicates the R rep function for tiling and repeating vectors

Source code in src/luna/radiology/mirp/roiClass.py

def rep(x, each=1, times=1):
    """ "This functions replicates the R rep function for tiling and repeating vectors"""
    each = int(each)
    times = int(times)

    if each > 1:
        x = np.repeat(x, repeats=each)

    if times > 1:
        x = np.tile(x, reps=times)

    return x

utilities

check_string(input_string)

Updates string characters that may lead to errors when written.

Source code in src/luna/radiology/mirp/utilities.py

def check_string(input_string):
    """Updates string characters that may lead to errors when written."""
    input_string = input_string.replace(" ", "_")
    input_string = input_string.replace(",", "_")
    input_string = input_string.replace(";", "_")
    input_string = input_string.replace(":", "_")
    input_string = input_string.replace('"', "_")
    input_string = input_string.replace("=", "_equal_")
    input_string = input_string.replace(">", "_greater_")
    input_string = input_string.replace("<", "_smaller_")
    input_string = input_string.replace("&", "_and_")
    input_string = input_string.replace("|", "_or_")

    return input_string

extract_roi_names(roi_list)

Extract the names of the regions of interest in roi_list :param roi_list: List of roi objects :return: List of roi names

Source code in src/luna/radiology/mirp/utilities.py

def extract_roi_names(roi_list):
    """
    Extract the names of the regions of interest in roi_list
    :param roi_list: List of roi objects
    :return: List of roi names
    """
    return [roi.name for roi in roi_list]

index_to_world(index, origin, spacing)

"Translates index to world coordinates

Source code in src/luna/radiology/mirp/utilities.py

def index_to_world(index, origin, spacing):
    """ "Translates index to world coordinates"""
    return origin + index * spacing

makedirs_check(path)

Checks if the given path is an existing directory structure, otherwise creates it.

Source code in src/luna/radiology/mirp/utilities.py

def makedirs_check(path):
    """
    Checks if the given path is an existing directory
    structure, otherwise creates it.
    """

    if not os.path.isdir(path):
        os.makedirs(path)

world_to_index(coord, origin, spacing, affine=None)

"Translates world coordinates to index

Source code in src/luna/radiology/mirp/utilities.py

def world_to_index(coord, origin, spacing, affine=None):

    """ "Translates world coordinates to index"""
    if not np.array_equal(spacing, np.diag(affine)):
        print("Warning, affine does not match spacing!", spacing, affine)
    if affine is None:
        return (coord - origin) / spacing
    else:
        return (coord - origin) / np.diag(affine)