Python SimpleITK.Image() Examples

def resample_images(
            image: sitk.Image,
            transforms: List[sitk.Euler3DTransform],
            interpolation: Interpolation,
            ) -> List[sitk.Image]:
        floating = reference = image
        default_value = np.float64(sitk.GetArrayViewFromImage(image).min())
        transforms = transforms[1:]  # first is identity
        images = [image]  # first is identity
        for transform in transforms:
            resampler = sitk.ResampleImageFilter()
            resampler.SetInterpolator(get_sitk_interpolator(interpolation))
            resampler.SetReferenceImage(reference)
            resampler.SetOutputPixelType(sitk.sitkFloat32)
            resampler.SetDefaultPixelValue(default_value)
            resampler.SetTransform(transform)
            resampled = resampler.Execute(floating)
            images.append(resampled)
        return images 

def gen_warp_field(shape, sigma=0.1, grid_size=3):
    '''Generate an spline warp field'''
    import SimpleITK as sitk
    # Initialize bspline transform
    args = shape+(sitk.sitkFloat32,)
    ref_image = sitk.Image(*args)
    tx = sitk.BSplineTransformInitializer(ref_image, [grid_size, grid_size])

    # Initialize shift in control points:
    # mesh size = number of control points - spline order
    p = sigma * np.random.randn(grid_size+3, grid_size+3, 2)

    # Anchor the edges of the image
    p[:, 0, :] = 0
    p[:, -1:, :] = 0
    p[0, :, :] = 0
    p[-1:, :, :] = 0

    # Set bspline transform parameters to the above shifts
    tx.SetParameters(p.flatten())

    # Compute deformation field
    displacement_filter = sitk.TransformToDisplacementFieldFilter()
    displacement_filter.SetReferenceImage(ref_image)
    displacement_field = displacement_filter.Execute(tx)

    return displacement_field 

def __init__(self,
                 srr_mask_sitk,
                 suffix="_mask",
                 sigma=2,
                 lower=0.5,
                 upper=100,
                 ):

        if not isinstance(srr_mask_sitk, sitk.Image):
            raise ValueError("Input must be of type sitk.Image")

        self._srr_mask_sitk = srr_mask_sitk
        self._suffix = suffix
        self._sigma = sigma
        self._lower = lower
        self._upper = upper

        self._mask_sitk = None
        self._mask = None 

def get_rigid_transforms(
            self,
            degrees_params: np.ndarray,
            translation_params: np.ndarray,
            image: sitk.Image,
            ) -> List[sitk.Euler3DTransform]:
        center_ijk = np.array(image.GetSize()) / 2
        center_lps = image.TransformContinuousIndexToPhysicalPoint(center_ijk)
        identity = np.eye(4)
        matrices = [identity]
        for degrees, translation in zip(degrees_params, translation_params):
            radians = np.radians(degrees).tolist()
            motion = sitk.Euler3DTransform()
            motion.SetCenter(center_lps)
            motion.SetRotation(*radians)
            motion.SetTranslation(translation.tolist())
            motion_matrix = self.transform_to_matrix(motion)
            matrices.append(motion_matrix)
        transforms = [self.matrix_to_transform(m) for m in matrices]
        return transforms 

def add_artifact(
            self,
            image: sitk.Image,
            transforms: List[sitk.Euler3DTransform],
            times: np.ndarray,
            interpolation: Interpolation,
            ):
        images = self.resample_images(image, transforms, interpolation)
        arrays = [sitk.GetArrayViewFromImage(im) for im in images]
        arrays = [array.transpose() for array in arrays]  # ITK to NumPy
        spectra = [self.fourier_transform(array) for array in arrays]
        self.sort_spectra(spectra, times)
        result_spectrum = np.empty_like(spectra[0])
        last_index = result_spectrum.shape[2]
        indices = (last_index * times).astype(int).tolist()
        indices.append(last_index)
        ini = 0
        for spectrum, fin in zip(spectra, indices):
            result_spectrum[..., ini:fin] = spectrum[..., ini:fin]
            ini = fin
        result_image = np.real(self.inv_fourier_transform(result_spectrum))
        return result_image.astype(np.float32) 

def get_brain_image_sitk(self):
        if self._sitk_brain_image is None:
            raise ValueError("Brain was not asked for. "
                             "Do not set option '-n' and run again.")

        return self._sitk_brain_image

    ##
    #       Get computed brain image mask
    # \date       2016-10-12 14:33:53+0100
    #
    # \param      self  The object
    #
    # \return     The brain mask as sitk.Image object
    # 

def get_input_image_sitk(self):
        if self._sitk is None:
            raise ValueError("Input image was not read yet.")

        return sitk.Image(self._sitk)

    ##
    # Gets the brain masked stack.
    # \date       2018-01-18 00:44:49+0000
    #
    # \param      self            The object
    # \param      filename        The filename
    # \param      extract_slices  Extract slices of stack; boolean
    #
    # \return     Returns image as Stack object holding obtained brain mask
    # 

def sitk_to_nib(image: sitk.Image) -> Tuple[np.ndarray, np.ndarray]:
    data = sitk.GetArrayFromImage(image).transpose()
    spacing = np.array(image.GetSpacing())
    direction = np.array(image.GetDirection())
    origin = image.GetOrigin()
    if len(direction) == 9:
        rotation = direction.reshape(3, 3)
    elif len(direction) == 4:  # ignore first dimension if 2D (1, 1, H, W)
        rotation_2d = direction.reshape(2, 2)
        rotation = np.eye(3)
        rotation[1:3, 1:3] = rotation_2d
        spacing = 1, *spacing
        origin = 0, *origin
    rotation = np.dot(FLIP_XY, rotation)
    rotation_zoom = rotation * spacing
    translation = np.dot(FLIP_XY, origin)
    affine = np.eye(4)
    affine[:3, :3] = rotation_zoom
    affine[:3, 3] = translation
    return data, affine 

def from_stack(cls,
                   stack,
                   compute_brain_image=False,
                   compute_brain_mask=True,
                   compute_skull_image=False,
                   dir_tmp=os.path.join(DIR_TMP, "BrainExtractionTool")):

        self = cls(compute_brain_image=compute_brain_image,
                   compute_brain_mask=compute_brain_mask,
                   compute_skull_image=compute_skull_image,
                   dir_tmp=dir_tmp)

        self._stack = stack
        self._sitk = sitk.Image(stack.sitk)

        return self

    ##
    #       Sets the sitk image for brain stripping
    # \date       2016-10-12 15:46:20+0100
    #
    # \param      self        The object
    # \param      sitk_image  The sitk image as sitk.Image object
    #
    # 

def _add_slice_contribution(slice, coverage_sitk):

        #
        slice_sitk = sitk.Image(slice.sitk)
        spacing = np.array(slice_sitk.GetSpacing())
        spacing[-1] = slice.get_slice_thickness()
        slice_sitk.SetSpacing(spacing)

        contrib_nda = sitk.GetArrayFromImage(slice_sitk)
        contrib_nda[:] = 1
        contrib_sitk = sitk.GetImageFromArray(contrib_nda)
        contrib_sitk.CopyInformation(slice_sitk)

        coverage_sitk += sitk.Resample(
            contrib_sitk,
            coverage_sitk,
            sitk.Euler3DTransform(),
            sitk.sitkNearestNeighbor,
            0,
            coverage_sitk.GetPixelIDValue(),
        )

        return coverage_sitk 

def verify_same_geometry(img_1: sitk.Image, img_2: sitk.Image):
    ori1, spacing1, direction1, size1 = img_1.GetOrigin(), img_1.GetSpacing(), img_1.GetDirection(), img_1.GetSize()
    ori2, spacing2, direction2, size2 = img_2.GetOrigin(), img_2.GetSpacing(), img_2.GetDirection(), img_2.GetSize()

    same_ori = np.all(np.isclose(ori1, ori2))
    if not same_ori: print("the origin does not match between the images")
    same_spac = np.all(np.isclose(spacing1, spacing2))
    if not same_spac: print("the spacing does not match between the images")
    same_dir = np.all(np.isclose(direction1, direction2))
    if not same_dir: print("the direction does not match between the images")
    same_size = np.all(np.isclose(size1, size2))
    if not same_size: print("the size does not match between the images")

    if same_ori and same_spac and same_dir and same_size:
        return True
    else:
        return False 

def __init__(self, stacks, reconstruction_sitk):
        self._stacks = stacks
        self._reconstruction_sitk = reconstruction_sitk

        self._coverage_sitk = None

    ##
    # Gets the slice coverage as Image. The (integer) intensity values reflect
    # the number of slices that have contributed to this particular voxel.
    # \date       2019-02-23 21:18:10+0000
    #
    # \param      self  The object
    #
    # \return     Slice coverage as sitk.Image uint8 image.
    # 

def show_slices(self):
        sitkh.plot_stack_of_slices(
            self.sitk, cmap="Greys_r", title=self.get_filename())

    # Write information of Stack to HDD to given directory:
    #  - sitk.Image object as entire volume
    #  - each single slice with its associated spatial transformation (optional)
    #  \param[in] directory string specifying where the output will be written to (default="/tmp/")
    #  \param[in] filename string specifying the filename. If not given the assigned one within Stack will be chosen.
    #  \param[in] write_slices boolean indicating whether each Slice of the stack shall be written (default=False) 

def get_coverage_sitk(self):
        if self._coverage_sitk is None:
            raise RuntimeError("Execute 'run' first")
        return sitk.Image(self._coverage_sitk)

    ##
    # Compute slice coverage
    # \date       2019-02-23 21:19:55+0000
    #
    # \param      self  The object
    # 

def run(self):

        # create zero image
        coverage_sitk = sitk.Image(self._reconstruction_sitk) * 0

        for i, stack in enumerate(self._stacks):
            print("Slices of stack %d/%d ... " % (i + 1, len(self._stacks)))

            # Add each individual slice contribution
            for slice in stack.get_slices():
                coverage_sitk = self._add_slice_contribution(
                    slice, coverage_sitk)

        # Cast to unsigned integer
        self._coverage_sitk = sitk.Cast(coverage_sitk, sitk.sitkUInt8)

    ##
    # Adds a slice contribution.
    # \date       2019-02-23 21:27:12+0000
    #
    # \param      slice          Slice as sl.Slice object
    # \param      coverage_sitk  sitk.Image reflecting the current iteration of
    #                            slice coverage
    #
    # \return     Updated slice contribution, sitk.Image
    # 

def from_filename(cls,
                      dir_input,
                      filename,
                      compute_brain_image=False,
                      compute_brain_mask=True,
                      compute_skull_image=False,
                      dir_tmp=os.path.join(DIR_TMP, "BrainExtractionTool")):

        self = cls(compute_brain_image=compute_brain_image,
                   compute_brain_mask=compute_brain_mask,
                   compute_skull_image=compute_skull_image,
                   dir_tmp=dir_tmp)
        self._sitk = sitkh.read_nifti_image_sitk(
            os.path.join(dir_input, "%s.nii.gz" % filename),
            sitk.sitkFloat64)

        return self

    ##
    # Initialize brain stripping class based on given sitk.Image object
    # \date       2016-10-12 12:18:35+0100
    #
    # \param      cls                  The cls
    # \param      sitk_image           The sitk image
    # \param      compute_brain_image  Boolean flag for computing brain image
    # \param      compute_brain_mask   Boolean flag for computing brain image
    #                                  mask
    # \param      compute_skull_image  Boolean flag for computing skull mask
    # \param      dir_tmp              Directory where temporary results are
    #                                  written to, string
    #
    # \return     object
    # 

def set_input_image_sitk(self, sitk_image):
        self._sitk = sitk.Image(sitk_image)

    ##
    #       Set flag of whether or not to compute the brain image
    # \date       2016-10-12 12:35:46+0100
    #
    # \param      self                 The object
    # \param      compute_brain_image  Boolean flag
    # 

def set_bet_options(self, bet_options):
        self._bet_options = bet_options

    ##
    #       Gets the input image
    # \date       2016-10-12 14:41:05+0100
    #
    # \param      self  The object
    #
    # \return     The input image as sitk.Image object
    # 

def get_mask_sitk(self):
        return sitk.Image(self._mask_sitk) 

def get_processed_mask_sitk(self):
        return sitk.Image(self._mask_sitk) 

def _extract_slices(self, slice_thickness, slice_numbers=None):

        slices = [None] * self._N_slices

        if slice_numbers is None:
            slice_numbers = range(0, self._N_slices)

        if len(slice_numbers) != self._N_slices:
            raise ValueError(
                "slice_numbers must correspond to the number of slices "
                "of the image volume")

        # Extract slices and add masks
        for i in range(0, self._N_slices):
            slices[i] = sl.Slice.from_sitk_image(
                slice_sitk=self.sitk[:, :, i:i + 1],
                filename=self._filename,
                slice_number=slice_numbers[i],
                slice_sitk_mask=self.sitk_mask[:, :, i:i + 1],
                slice_thickness=slice_thickness,
            )

        return slices

    # Create a binary mask consisting of ones
    #  \return binary_mask as sitk.Image object consisting of ones 

def myshow3d(img, xslices=[], yslices=[], zslices=[], title=None, margin=0.05, dpi=80):
    size = img.GetSize()
    img_xslices = [img[s,:,:] for s in xslices]
    img_yslices = [img[:,s,:] for s in yslices]
    img_zslices = [img[:,:,s] for s in zslices]
  
    maxlen = max(len(img_xslices), len(img_yslices), len(img_zslices))
  
      
    img_null = sitk.Image([0,0], img.GetPixelIDValue(), img.GetNumberOfComponentsPerPixel())
  
    img_slices = []
    d = 0
  
    if len(img_xslices):
	img_slices += img_xslices + [img_null]*(maxlen-len(img_xslices))
	d += 1
      
    if len(img_yslices):
	img_slices += img_yslices + [img_null]*(maxlen-len(img_yslices))
	d += 1
    
    if len(img_zslices):
	img_slices += img_zslices + [img_null]*(maxlen-len(img_zslices))
	d +=1
  
    if maxlen != 0:
	if img.GetNumberOfComponentsPerPixel() == 1:
	    img = sitk.Tile(img_slices, [maxlen,d])
	#TODO check in code to get Tile Filter working with VectorImages
	else:
	    img_comps = []
	    for i in range(0,img.GetNumberOfComponentsPerPixel()):
		img_slices_c = [sitk.VectorIndexSelectionCast(s, i) for s in img_slices]
		img_comps.append(sitk.Tile(img_slices_c, [maxlen,d]))
	    img = sitk.Compose(img_comps)
	    
    
    myshow(img, title, margin, dpi) 

def copy_geometry(image: sitk.Image, ref: sitk.Image):
    image.SetOrigin(ref.GetOrigin())
    image.SetDirection(ref.GetDirection())
    image.SetSpacing(ref.GetSpacing())
    return image 

def apply_transform_to_file(
        input_path: TypePath,
        transform,  # : Transform seems to create a circular import (TODO)
        output_path: TypePath,
        type: str = INTENSITY,
        ):
    from . import Image, ImagesDataset, Subject
    subject = Subject(image=Image(input_path, type))
    dataset = ImagesDataset([subject], transform=transform)
    transformed = dataset[0]
    dataset.save_sample(transformed, dict(image=output_path)) 

def _get_subject_from_tensor(tensor: torch.Tensor) -> Subject:
        subject_dict = {}
        for channel_index, channel_tensor in enumerate(tensor):
            name = f'channel_{channel_index}'
            image = Image(tensor=channel_tensor, type=INTENSITY)
            subject_dict[name] = image
        subject = Subject(subject_dict)
        dataset = ImagesDataset([subject])
        sample = dataset[0]
        return sample 

def crop(self, index_ini, index_fin):
        new_origin = nib.affines.apply_affine(self.affine, index_ini)
        new_affine = self.affine.copy()
        new_affine[:3, 3] = new_origin
        i0, j0, k0 = index_ini
        i1, j1, k1 = index_fin
        patch = self.data[0, i0:i1, j0:j1, k0:k1].clone()
        return Image(tensor=patch, affine=new_affine, type=self.type) 

def as_sitk(self) -> sitk.Image:
        return nib_to_sitk(self[DATA][0], self[AFFINE]) 

def get_bspline_transform(
            image: sitk.Image,
            num_control_points: Tuple[int, int, int],
            coarse_field: np.ndarray,
            ) -> sitk.BSplineTransformInitializer:
        mesh_shape = [n - SPLINE_ORDER for n in num_control_points]
        bspline_transform = sitk.BSplineTransformInitializer(image, mesh_shape)
        parameters = coarse_field.flatten(order='F').tolist()
        bspline_transform.SetParameters(parameters)
        return bspline_transform 

def sitk_to_nib(image: sitk.Image):
        return sitk_to_nib(image) 

def run_segmentation_propagation(self):

        if self._stack is None or self._template is None:
            raise ValueError("Specify stack and template first")

        # Choose interpolator
        try:
            interpolator_str = self._interpolator
            interpolator = eval("sitk.sitk" + interpolator_str)
        except:
            raise ValueError("Error: interpolator is not known")

        self._stack_sitk = sitk.Image(self._stack.sitk)

        # Register stack to template
        if self._registration_method is not None:
            self._registration_method.set_fixed(self._template)
            self._registration_method.set_moving(self._stack)
            self._registration_method.use_fixed_mask(self._use_template_mask)
            self._registration_method.run()

            self._registration_transform_sitk = self._registration_method.get_registration_transform_sitk()
            self._registration_transform_sitk = eval(
                "sitk." + self._registration_transform_sitk.GetName() + "(self._registration_transform_sitk.GetInverse())")

            self._stack_sitk = sitkh.get_transformed_sitk_image(
                self._stack_sitk, self._registration_transform_sitk)

        # Propagate mask
        self._stack_sitk_mask = sitk.Resample(self._template.sitk_mask, self._stack_sitk, sitk.Euler3DTransform(
        ), interpolator, 0, self._template.sitk_mask.GetPixelIDValue())

        # Dilate mask
        if self._dilation_radius > 0:

            stack_mask_morpher = stmorph.StackMaskMorphologicalOperations.from_sitk_mask(
                mask_sitk=self._stack_sitk_mask,
                dilation_radius=self._dilation_radius,
                dilation_kernel=self._dilation_kernel,
                use_dilation_in_plane_only=self._use_dilation_in_plane_only,
            )
            stack_mask_morpher.run_dilation()
            self._stack_sitk_mask = stack_mask_morpher.get_processed_mask_sitk()

            # sitkh.show_sitk_image(self._stack_sitk_mask)