Python SimpleITK.GetArrayViewFromImage() 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 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 add_artifact(
            self,
            image: sitk.Image,
            spikes_positions: np.ndarray,
            intensity_factor: float,
            ):
        array = sitk.GetArrayViewFromImage(image).transpose()
        spectrum = self.fourier_transform(array)
        shape = np.array(spectrum.shape)
        mid_shape = shape // 2
        indices = np.floor(spikes_positions * shape).astype(int)
        for index in indices:
            diff = index - mid_shape
            i, j, k = mid_shape + diff
            spectrum[i, j, k] = spectrum.max() * intensity_factor
            # If we wanted to add a pure cosine, we should add spikes to both
            # sides of k-space. However, having only one is a better
            # representation og the actual cause of the artifact in real
            # scans.
            #i, j, k = mid_shape - diff
            #spectrum[i, j, k] = spectrum.max() * intensity_factor
        result = np.real(self.inv_fourier_transform(spectrum))
        return result.astype(np.float32) 

def get_borders_mean(image, filter_otsu=True):
    # pylint: disable=bad-whitespace
    array = sitk.GetArrayViewFromImage(image)
    borders_tuple = (
        array[ 0,  :,  :],
        array[-1,  :,  :],
        array[ :,  0,  :],
        array[ :, -1,  :],
        array[ :,  :,  0],
        array[ :,  :, -1],
    )
    borders_flat = np.hstack([border.ravel() for border in borders_tuple])
    if not filter_otsu:
        return borders_flat.mean()
    borders_reshaped = borders_flat.reshape(1, 1, -1)
    borders_image = sitk.GetImageFromArray(borders_reshaped)
    otsu = sitk.OtsuThresholdImageFilter()
    otsu.Execute(borders_image)
    threshold = otsu.GetThreshold()
    values = borders_flat[borders_flat < threshold]
    if values.any():
        default_value = values.mean()
    else:
        default_value = borders_flat.mean()
    return default_value 

def my_show(*args):
    for count, src in enumerate(args):
        nda = sitk.GetArrayViewFromImage(src)
        plt.subplot(1, len(args), count + 1)
        plt.imshow(nda)
    plt.show() 

def sitk_to_np_no_copy(image_sitk):
    return sitk.GetArrayViewFromImage(image_sitk) 

def sitk_to_np(image_sitk, type=None):
    if type is None:
        return sitk.GetArrayFromImage(image_sitk)
    else:
        return sitk.GetArrayViewFromImage(image_sitk).astype(type) 

def surface_distance(label_image_0, label_image_1):
    # code adapted from https://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/Python_html/34_Segmentation_Evaluation.html
    try:
        # calculate distances on label contours
        reference_distance_map = sitk.SignedMaurerDistanceMap(label_image_1, squaredDistance=False, useImageSpacing=True)
        reference_distance_map_arr = sitk.GetArrayViewFromImage(reference_distance_map)
        reference_surface = sitk.LabelContour(label_image_1)
        reference_surface_arr = sitk.GetArrayViewFromImage(reference_surface)

        segmented_distance_map = sitk.SignedMaurerDistanceMap(label_image_0, squaredDistance=False, useImageSpacing=True)
        segmented_distance_map_arr = sitk.GetArrayViewFromImage(segmented_distance_map)
        segmented_surface = sitk.LabelContour(label_image_0)
        segmented_surface_arr = sitk.GetArrayViewFromImage(segmented_surface)

        seg2ref_distances = np.abs(reference_distance_map_arr[segmented_surface_arr == 1])
        ref2seg_distances = np.abs(segmented_distance_map_arr[reference_surface_arr == 1])

        all_surface_distances = np.concatenate([seg2ref_distances, ref2seg_distances])

        # # Multiply the binary surface segmentations with the distance maps. The resulting distance
        # # maps contain non-zero values only on the surface (they can also contain zero on the surface)
        # seg2ref_distance_map = reference_distance_map * sitk.Cast(segmented_surface, sitk.sitkFloat32)
        # ref2seg_distance_map = segmented_distance_map * sitk.Cast(reference_surface, sitk.sitkFloat32)
        #
        # statistics_image_filter = sitk.StatisticsImageFilter()
        # # Get the number of pixels in the reference surface by counting all pixels that are 1.
        # statistics_image_filter.Execute(reference_surface)
        # num_reference_surface_pixels = int(statistics_image_filter.GetSum())
        # # Get the number of pixels in the reference surface by counting all pixels that are 1.
        # statistics_image_filter.Execute(segmented_surface)
        # num_segmented_surface_pixels = int(statistics_image_filter.GetSum())
        #
        # # Get all non-zero distances and then add zero distances if required.
        # seg2ref_distance_map_arr = sitk.GetArrayViewFromImage(seg2ref_distance_map)
        # seg2ref_distances = list(seg2ref_distance_map_arr[seg2ref_distance_map_arr != 0])
        # seg2ref_distances = seg2ref_distances + list(np.zeros(num_segmented_surface_pixels - len(seg2ref_distances)))
        # ref2seg_distance_map_arr = sitk.GetArrayViewFromImage(ref2seg_distance_map)
        # ref2seg_distances = list(ref2seg_distance_map_arr[ref2seg_distance_map_arr != 0])
        # ref2seg_distances = ref2seg_distances + list(np.zeros(num_reference_surface_pixels - len(ref2seg_distances)))
        #
        # all_surface_distances = seg2ref_distances + ref2seg_distances

        current_mean_surface_distance = np.mean(all_surface_distances)
        current_median_surface_distance = np.median(all_surface_distances)
        current_std_surface_distance = np.std(all_surface_distances)
        current_max_surface_distance = np.max(all_surface_distances)
    except:
        current_mean_surface_distance = np.nan
        current_median_surface_distance = np.nan
        current_std_surface_distance = np.nan
        current_max_surface_distance = np.nan
        pass

    return current_mean_surface_distance, current_median_surface_distance, current_std_surface_distance, current_max_surface_distance 

def to_itk_image(image_like):

    if isinstance(image_like, itk.Image):
        return image_like

    if is_arraylike(image_like):
        array = np.asarray(image_like)
        case_use_view = array.flags['OWNDATA']
        if have_dask and isinstance(image_like, dask.array.core.Array):
            case_use_view = False
        array = np.ascontiguousarray(array)
        if case_use_view:
            image_from_array = itk.image_view_from_array(array)
        else:
            image_from_array = itk.image_from_array(array)
        return image_from_array
    elif have_vtk and isinstance(image_like, vtk.vtkImageData):
        from vtk.util import numpy_support as vtk_numpy_support
        array = vtk_numpy_support.vtk_to_numpy(
            image_like.GetPointData().GetScalars())
        array.shape = tuple(image_like.GetDimensions())[::-1]
        image_from_array = itk.image_view_from_array(array)
        image_from_array.SetSpacing(image_like.GetSpacing())
        image_from_array.SetOrigin(image_like.GetOrigin())
        return image_from_array
    elif have_simpleitk and isinstance(image_like, sitk.Image):
        array = sitk.GetArrayViewFromImage(image_like)
        image_from_array = itk.image_view_from_array(array)
        image_from_array.SetSpacing(image_like.GetSpacing())
        image_from_array.SetOrigin(image_like.GetOrigin())
        direction = image_like.GetDirection()
        npdirection = np.asarray(direction)
        npdirection = np.reshape(npdirection, (-1, 3))
        itkdirection = itk.matrix_from_array(npdirection)
        image_from_array.SetDirection(itkdirection)
        return image_from_array
    elif have_imagej:
        import imglyb
        if isinstance(image_like,
                      imglyb.util.ReferenceGuardingRandomAccessibleInterval):
            array = imglyb.to_numpy(image_like)
            image_from_array = itk.image_view_from_array(array)
            return image_from_array
    elif isinstance(image_like, itk.ProcessObject):
        return itk.output(image_like)

    return None