crop.py

"""
Copyright 2020 Division of Medical Image Computing, German Cancer Research Center (DKFZ), Heidelberg, Germany

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
"""

import os
import shutil
import pickle
import numpy as np

from loguru import logger
from multiprocessing.pool import Pool
from pathlib import Path
from typing import List, Tuple, Sequence
from scipy.ndimage import binary_fill_holes

from nndet.io.paths import get_case_id_from_path
from nndet.io.load import load_case_from_list


def create_nonzero_mask(data: np.ndarray) -> np.ndarray:
    """
    Create a nonzero mask from data

    Args:
        data (np.ndarray): input data [C, X, Y, Z]

    Returns:
        np.ndarray: binary mask on nonzero regions [X, Y, Z]
    """
    assert len(data.shape) == 4 or len(data.shape) == 3, \
        "data must have shape (C, X, Y, Z) or shape (C, X, Y)"
    nonzero_mask = np.max(data != 0, axis=0)
    nonzero_mask = binary_fill_holes(nonzero_mask.astype(bool))
    return nonzero_mask


def get_bbox_from_mask(mask: np.ndarray, outside_value: int = 0) -> List[Tuple]:
    """
    Create a bounding box from a mask

    Args:
        mask (np.ndarray): mask [X, Y, Z]
        outside_value (int): background value

    Returns:
        np.ndarray: [(dim0_min, dim0_max), (dim1_min, dim1_max), (dim2_min, dim2_max))
    """
    mask_voxel_coords = (mask != outside_value).nonzero()

    min0idx = int(np.min(mask_voxel_coords[0]))
    max0idx = int(np.max(mask_voxel_coords[0])) + 1

    min1idx = int(np.min(mask_voxel_coords[1]))
    max1idx = int(np.max(mask_voxel_coords[1])) + 1
    idx = [(min0idx, max0idx), (min1idx, max1idx)]

    if len(mask_voxel_coords) == 3:
        min2idx = int(np.min(mask_voxel_coords[2]))
        max2idx = int(np.max(mask_voxel_coords[2])) + 1
        idx.append((min2idx, max2idx))
    return idx


def crop_to_bbox_no_channels(image, bbox: Sequence[Sequence[int]]):
    """
    Crops image to bounding box (in spatial dimensions)

    Args:
        image (arraylike): 2d or 3d array
        bbox (Sequence[Sequence[int]]): bounding box coordinated in an interleaved fashion
            (e.g. (x1, x2), (y1, y2), (z1, z2))

    Returns:
        arraylike: cropped array
    """
    resizer = tuple([slice(_dim[0], _dim[1]) for _dim in bbox])
    return image[resizer]


def crop_to_bbox(data: np.ndarray, bbox: Sequence[Sequence[int]]):
    """
    Crops image to bounding box (performed per channel)

    Args:
        data (np.ndarray): 3d or 4d array [C, X, Y, (Z)]
        bbox (Sequence[Sequence[int]]): bounding box coordinated in an interleaved fashion
            (e.g. (x1, x2), (y1, y2), (z1, z2))

    Returns:
        np.ndarray: cropped array
    """
    cropped_data = []
    for c in range(data.shape[0]):
        cropped = crop_to_bbox_no_channels(data[c], bbox)
        cropped_data.append(cropped)
    data = np.stack(cropped_data)
    return data


def crop_to_nonzero(data, seg=None, nonzero_label=-1):
    """
    Crop data to nonzero region of data

    Args:
        data (np.ndarray): data to crop
        seg (np.ndarray): segmentation
        nonzero_label (int): nonzero label is written into segmentation map
            where only background was found

    Returns:
        np.ndarray: cropped data
        np.ndarray: cropped and filled (with nonzero_label) segmentation
        List[Tuple[int]]: bounding box of nonzero region
    """
    nonzero_mask = create_nonzero_mask(data)
    bbox = get_bbox_from_mask(nonzero_mask, 0)

    data = crop_to_bbox(data, bbox)
    seg = crop_to_bbox(seg, bbox)
    nonzero_mask = crop_to_bbox_no_channels(nonzero_mask, bbox)[None]

    if seg is not None:
        seg[(seg == 0) & (nonzero_mask == 0)] = nonzero_label
    else:
        nonzero_mask = nonzero_mask.astype(np.int32)
        nonzero_mask[nonzero_mask == 0] = nonzero_label
        nonzero_mask[nonzero_mask > 0] = 0
        seg = nonzero_mask
    return data, seg, bbox


class ImageCropper(object):
    def __init__(self, num_processes: int, output_dir: Path = None):
        """
        Helper class to crop images to non zero region (must hold for all modalities)
        In the case of BRaTS and ISLES data this results in a significant reduction in image size

        Args:
            num_processes (int): number of processes to use for cropping
            output_dir (Path): path to output directory
        """
        self.output_dir = Path(output_dir) if output_dir is not None else None
        self.num_processes = num_processes
        self.maybe_init_output_dir()

    def maybe_init_output_dir(self):
        if self.output_dir is not None:
            (self.output_dir / "imagesTr").mkdir(parents=True, exist_ok=True)
            (self.output_dir / "labelsTr").mkdir(parents=True, exist_ok=True)

    def run_cropping(self,
                     case_files: List[List[Path]],
                     overwrite_existing: bool = False,
                     output_dir: Path = None,
                     copy_gt_data: bool = True,
                     ):
        """
        Crops data to non zero region and saves them into output_dir
        Optional: also copies ground truth data

        Args:
            case_files (List[List[Path]]): list with all cases in the structure [Case[Case Files]];
                where case files are sorted to corresponding modalities (last file is the label file)
            overwrite_existing (bool): overwrite existing crops
            output_dir (Path): path to output directory
            copy_gt_data (bool): copies ground truth data to output directory
        """
        if output_dir is not None:
            self.output_dir = Path(output_dir)
            self.maybe_init_output_dir()

        if copy_gt_data:
            self.copy_gt_data(case_files)

        list_of_args = []
        for _i, case in enumerate(case_files):
            case_id = get_case_id_from_path(str(case[0]))
            assert not case_id.endswith(".gz") and not case_id.endswith(".nii")
            list_of_args.append((case, case_id, overwrite_existing))

        if self.num_processes == 0:
            for a in list_of_args:
                self.process_data(*a)
        else:
            with Pool(processes=self.num_processes) as p:
                p.starmap(self.process_data, list_of_args)

    def copy_gt_data(self, case_files: List[List[Path]]):
        """
        Copy ground truth to output directory
        """
        output_dir_gt = self.output_dir / "labelsTr"
        if output_dir_gt.is_dir():
            shutil.rmtree(output_dir_gt)
        source_dir_gt = case_files[0][-1].parent
        shutil.copytree(source_dir_gt, output_dir_gt)

    def process_data(self, case: List[Path], case_id: str, overwrite_existing: bool = False):
        """
        Extract nonzero region from all cases and create a single array where segmentation
        is located in the last channel and save as npz (saved in key `data`)
        Additional properties per case are saved inside a pkl file

        Args:
            case (List[Path]): list of paths to data and label (label is always at the last position
                and data is sorted after modalities)
            case_id (str): case identifier
            overwrite_existing (bool): overwrite existing data
        """
        try:
            logger.info(f"Processing case {case_id}")
            npz_exists = (self.output_dir / "imagesTr" / f"{case_id}.npz").is_file()
            pkl_exists = (self.output_dir / "imagesTr" / f"{case_id}.pkl").is_file()

            if not (npz_exists and pkl_exists) or overwrite_existing:
                data, seg, properties = self.load_crop_from_list_of_files(case[:-1], case[-1])

                all_data = np.vstack((data, seg))
                np.savez_compressed(self.output_dir / "imagesTr" / f"{case_id}.npz", data=all_data)
                with open(self.output_dir / "imagesTr" / f"{case_id}.pkl", 'wb') as f:
                    pickle.dump(properties, f)
            else:
                logger.warning(f"Case {case_id} already exists and overwrite is deactivated")
        except Exception as e:
            logger.info(f"exception in: {case_id}: {e}")
            raise e

    @staticmethod
    def load_crop_from_list_of_files(data_files: List[Path], seg_file: Path = None):
        """
        Load and crop form list of files

        Args:
            data_files (List[Path]): paths to data files
            seg_file (Path): pth to segmentation

        Returns:
            np.ndarray: cropped data
            np.ndarray: cropped (and filled segmentation: -1 where no forground exists) label
            dict: additional properties
                `original_size_of_raw_data`: original shape of data (correctly reordered)
                `original_spacing`: original spacing (correctly reordered)
                `list_of_data_files`: paths of data files
                `seg_file`: path to label file
                `itk_origin`: origin in world coordinates
                `itk_spacing`: spacing in world coordinates
                `itk_direction`: direction in world coordinates
                `crop_bbox`: List[Tuple[int]] cropped bounding box
                `classes`: present classes in segmentation
                `size_after_cropping`: size after cropping
        """
        data, seg, properties = load_case_from_list(data_files, seg_file)
        return ImageCropper.crop(data, properties, seg)

    @staticmethod
    def crop(data: np.ndarray, properties: dict, seg: np.ndarray = None):
        """
        Crop data and segmentation to non zero region

        Args:
            data (np.ndarray): data to crop [C, X, Y, Z]
            properties (dict): additional properties
            seg (np.ndarray): segmentation [1, X, Y, Z]

        Returns:
            data (np.ndarray): data to crop [C, X, Y, Z]
            seg (np.ndarray): segmentation [1, X, Y, Z]
            properties (dict): newly added properties
                `crop_bbox`: List[Tuple[int]] cropped bounding box
                `classes`: present classes in segmentation
                `size_after_cropping`: size after cropping
        """
        shape_before = data.shape
        data, seg, bbox = crop_to_nonzero(data, seg, nonzero_label=-1)
        shape_after = data.shape
        logger.info(f"Shape before crop {shape_before}; after crop {shape_after}; "
                    f"spacing {np.array(properties['original_spacing'])}")

        properties["crop_bbox"] = bbox
        properties['classes'] = np.unique(seg)
        seg[seg < -1] = 0
        properties["size_after_cropping"] = data[0].shape
        return data, seg, properties