Python skimage.morphology.disk() Examples

def dilate(data, size, shape, dim=None):
    """
    Dilate data using ball structuring element
    :param data: Image or numpy array: 2d or 3d array
    :param size: int: If shape={'square', 'cube'}: Corresponds to the length of an edge (size=1 has no effect).
    If shape={'disk', 'ball'}: Corresponds to the radius, not including the center element (size=0 has no effect).
    :param shape: {'square', 'cube', 'disk', 'ball'}
    :param dim: {0, 1, 2}: Dimension of the array which 2D structural element will be orthogonal to. For example, if
    you wish to apply a 2D disk kernel in the X-Y plane, leaving Z unaffected, parameters will be: shape=disk, dim=2.
    :return: numpy array: data dilated
    """
    if isinstance(data, Image):
        im_out = data.copy()
        im_out.data = dilate(data.data, size, shape, dim)
        return im_out
    else:
        return dilation(data, selem=_get_selem(shape, size, dim), out=None) 

def addNoiseAndGray(surf):
    # https://stackoverflow.com/questions/34673424/how-to-get-numpy-array-of-rgb-colors-from-pygame-surface
    imgdata = pygame.surfarray.array3d(surf)
    imgdata = imgdata.swapaxes(0, 1)
    # print('imgdata shape %s' % imgdata.shape)  # shall be IMG_HEIGHT * IMG_WIDTH
    imgdata2 = noise_generator('s&p', imgdata)

    img2 = Image.fromarray(np.uint8(imgdata2))
    # img2.save('/home/zhichyu/Downloads/2sp.jpg')
    grayscale2 = ImageOps.grayscale(img2)
    # grayscale2.save('/home/zhichyu/Downloads/2bw2.jpg')
    # return grayscale2

    array = np.asarray(np.uint8(grayscale2))
    # print('array.shape %s' % array.shape)
    selem = disk(random.randint(0, 1))
    eroded = erosion(array, selem)
    return eroded 

def masked(img, gt, mask, alpha=1):
    """Returns image with GT lung field outlined with red, predicted lung field
    filled with blue."""
    rows, cols = img.shape[:2]
    color_mask = np.zeros((rows, cols, 3))
    boundary = morphology.dilation(gt, morphology.disk(3)) ^ gt
    color_mask[mask == 1] = [0, 0, 1]
    color_mask[boundary == 1] = [1, 0, 0]
    
    img_hsv = color.rgb2hsv(img)
    color_mask_hsv = color.rgb2hsv(color_mask)

    img_hsv[..., 0] = color_mask_hsv[..., 0]
    img_hsv[..., 1] = color_mask_hsv[..., 1] * alpha

    img_masked = color.hsv2rgb(img_hsv)
    return img_masked 

def _get_blob_mask(ROI_image, thresh, thresh_block_size, is_light_background, analysis_type):
    # get binary image, 
    if is_light_background:
        ## apply a median filter to reduce rough edges / sharpen the boundary btw worm and background
        ROI_image_th = cv2.medianBlur(ROI_image, 3)
        ROI_mask = ROI_image_th < thresh
    else:
        if analysis_type == "PHARYNX":
            # for fluorescent pharynx labeled images, refine the threshold with a local otsu (http://scikit-image.org/docs/dev/auto_examples/plot_local_otsu.html)
            # this compensates for local variations in brightness in high density regions, when many worms are close to each other
            ROI_rank_otsu = skf.rank.otsu(ROI_image, skm.disk(thresh_block_size))
            ROI_mask = (ROI_image>ROI_rank_otsu)
            # as a local threshold introcudes artifacts at the edge of the mask, also use a global threshold to cut these out
            ROI_mask &= (ROI_image>=thresh)
        else:
            # this case applies for example to worms where the whole body is fluorecently labeled
            ROI_image_th = cv2.medianBlur(ROI_image, 3)
            ROI_mask = ROI_image_th >= thresh
        
    ROI_mask &= (ROI_image != 0)
    ROI_mask = ROI_mask.astype(np.uint8)

    return ROI_mask, thresh # returning thresh here seems redundant, as it isn't actually changed 

def masked(img, gt, mask, alpha=1):
    """Returns image with GT lung field outlined with red, predicted lung field
    filled with blue."""
    rows, cols = img.shape[:2]
    color_mask = np.zeros((rows, cols, 3))
    boundary = morphology.dilation(gt, morphology.disk(3)) ^ gt
    color_mask[mask == 1] = [0, 0, 1]
    color_mask[boundary == 1] = [1, 0, 0]
    
    img_hsv = color.rgb2hsv(img)
    color_mask_hsv = color.rgb2hsv(color_mask)

    img_hsv[..., 0] = color_mask_hsv[..., 0]
    img_hsv[..., 1] = color_mask_hsv[..., 1] * alpha

    img_masked = color.hsv2rgb(img_hsv)
    return img_masked 

def get_dark_mask(full_data):
    #get darker objects that are unlikely to be worm
    if full_data.shape[0] < 2:
        #nothing to do here returning
        return np.zeros((full_data.shape[1], full_data.shape[2]), np.uint8)
    
    #this mask shoulnd't contain many worms
    img_h = cv2.medianBlur(np.max(full_data, axis=0), 5)
    #this mask is likely to contain a lot of worms
    img_l = cv2.medianBlur(np.min(full_data, axis=0), 5)
    
    #this is the difference (the tagged pixels should be mostly worms)
    img_del = img_h-img_l
    th_d = threshold_otsu(img_del)
    
    #this is the maximum of the minimum pixels of the worms...
    th = np.max(img_l[img_del>th_d])
    #this is what a darkish mask should look like
    dark_mask = cv2.dilate((img_h<th).astype(np.uint8), disk(11))
    
    return dark_mask 

def erode(data, size, shape, dim=None):
    """
    Dilate data using ball structuring element
    :param data: Image or numpy array: 2d or 3d array
    :param size: int: If shape={'square', 'cube'}: Corresponds to the length of an edge (size=1 has no effect).
    If shape={'disk', 'ball'}: Corresponds to the radius, not including the center element (size=0 has no effect).
    :param shape: {'square', 'cube', 'disk', 'ball'}
    :param dim: {0, 1, 2}: Dimension of the array which 2D structural element will be orthogonal to. For example, if
    you wish to apply a 2D disk kernel in the X-Y plane, leaving Z unaffected, parameters will be: shape=disk, dim=2.
    :return: numpy array: data dilated
    """
    if isinstance(data, Image):
        im_out = data.copy()
        im_out.data = erode(data.data, size, shape, dim)
        return im_out
    else:
        return erosion(data, selem=_get_selem(shape, size, dim), out=None) 

def get_diffraction_test_image(self, dtype=np.float32):
        image_x, image_y = self.image_x, self.image_y
        cx, cy = image_x / 2, image_y / 2
        image = np.zeros((image_y, image_x), dtype=np.float32)
        iterator = zip(self._x_list, self._y_list, self._intensity_list)
        for x, y, i in iterator:
            if self.diff_intensity_reduction is not False:
                dr = np.hypot(x - cx, y - cy)
                i = self._get_diff_intensity_reduction(dr, i)
            image[y, x] = i
        disk = morphology.disk(self.disk_r, dtype=dtype)
        image = convolve2d(image, disk, mode="same")
        if self.rotation != 0:
            image = rotate(image, self.rotation, reshape=False)
        if self.blur != 0:
            image = gaussian_filter(image, self.blur)
        if self._background_lorentz_width is not False:
            image += self._get_background_lorentz()
        if self.intensity_noise is not False:
            noise = np.random.random((image_y, image_x)) * self.intensity_noise
            image += noise
        return image 

def local_max_roll(fm, k0, k1, diff):
    max_ls = []
    for ksize in range(k0, k1):
        selem = disk(ksize)
        fm_max = local_max(fm, selem, diff)
        max_ls.append(fm_max)
    lmax = (np.stack(max_ls, axis=0).sum(axis=0) > 0).astype(np.uint8)
    nlab, max_lab = cv2.connectedComponents(lmax)
    max_res = np.zeros_like(lmax)
    for lb in range(1, nlab):
        area = max_lab == lb
        if np.sum(area) > 1:
            crds = tuple(int(np.median(c)) for c in np.where(area))
            max_res[crds] = 1
        else:
            max_res[np.where(area)] = 1
    return max_res 

def dilation(x, radius=3):
    """ Return greyscale morphological dilation of an image,
    see `skimage.morphology.dilation <http://scikit-image.org/docs/dev/api/skimage.morphology.html#skimage.morphology.dilation>`_.

    Parameters
    -----------
    x : 2D array image.
    radius : int for the radius of mask.
    """
    from skimage.morphology import disk, dilation
    mask = disk(radius)
    x = dilation(x, selem=mask)
    return x




## Sequence 

def dilation(x, radius=3):
    """ Return greyscale morphological dilation of an image,
    see `skimage.morphology.dilation <http://scikit-image.org/docs/dev/api/skimage.morphology.html#skimage.morphology.dilation>`_.

    Parameters
    -----------
    x : 2D array image.
    radius : int for the radius of mask.
    """
    from skimage.morphology import disk, dilation
    mask = disk(radius)
    x = dilation(x, selem=mask)
    return x




## Sequence 

def finalProcessingSpur(s, params):
    logging.info(f"{s['filename']} - \tfinalProcessingSpur")
    disk_radius = int(params.get("disk_radius", "25"))
    selem = disk(disk_radius)
    mask = s["img_mask_use"]
    mask_opened = binary_opening(mask, selem)
    mask_spur = ~mask_opened & mask

    io.imsave(s["outdir"] + os.sep + s["filename"] + "_spur.png", img_as_ubyte(mask_spur))

    prev_mask = s["img_mask_use"]
    s["img_mask_use"] = mask_opened

    s.addToPrintList("spur_pixels",
                     printMaskHelper(params.get("mask_statistics", s["mask_statistics"]), prev_mask, s["img_mask_use"]))

    if len(s["img_mask_use"].nonzero()[0]) == 0:  # add warning in case the final tissue is empty
        logging.warning(
            f"{s['filename']} - After BasicModule.finalProcessingSpur NO tissue remains detectable! Downstream modules likely to be incorrect/fail")
        s["warnings"].append(
            f"After BasicModule.finalProcessingSpur NO tissue remains detectable! Downstream modules likely to be incorrect/fail") 

def ps_disk(radius):
    r"""
    Creates circular disk structuring element for morphological operations

    Parameters
    ----------
    radius : float or int
        The desired radius of the structuring element

    Returns
    -------
    strel : 2D-array
        A 2D numpy bool array of the structring element
    """
    rad = int(np.ceil(radius))
    other = np.ones((2 * rad + 1, 2 * rad + 1), dtype=bool)
    other[rad, rad] = False
    disk = spim.distance_transform_edt(other) < radius
    return disk 

def test_correct_disk_x_y_and_radius_random(self):
        x, y, px, py = 56, 48, 4, 5
        x, y = randint(45, 55, size=(py, px)), randint(45, 55, size=(py, px))
        r = randint(20, 40, size=(py, px))
        s = mdtd.generate_4d_data(
            probe_size_x=px,
            probe_size_y=py,
            image_size_x=120,
            image_size_y=100,
            disk_x=x,
            disk_y=y,
            disk_r=5,
            disk_I=20,
            ring_x=x,
            ring_y=y,
            ring_r=r,
            ring_I=5,
            blur=True,
            downscale=False,
        )
        s_com = s.center_of_mass()
        s_r = s.radial_average(centre_x=s_com.inav[0].data, centre_y=s_com.inav[1].data)
        s_r = s_r.isig[15:]  # Do not include the disk
        r -= 15  # Need to shift the radius, due to not including the disk
        assert (s_r.data.argmax(axis=-1) == r).all() 

def pair_to_rgb(gen_img, tar_img, background='black', use_dilation=False, disk_value=2):
    """
    args:
        gen_img: (ndarray) in [img_rows, img_cols], dytpe=unit8
        tar_img: (ndarray) in [img_rows, img_cols], dytpe=unit8
        background: (str) ['black', 'white']
    return:
        rgb_img: red -> false positive;
                 green -> true positive;
                 blue -> false positive;
    """
    # enhance outline border
    if use_dilation:
        gen_img = dilation(gen_img, disk(disk_value))
        tar_img = dilation(tar_img, disk(disk_value))

    if background == "black":
        # saving rgb results
        rgb_img = np.zeros((gen_img.shape[0], gen_img.shape[1], 3), np.uint8)
        # assign false negative as red channel
        rgb_img[:, :, 0][np.logical_and(gen_img == 255, tar_img == 0)] = 255
        # assign true positive as green channel
        rgb_img[:, :, 1][np.logical_and(gen_img == 255, tar_img == 255)] = 255
        # assign false positive as blue channel
        rgb_img[:, :, 2][np.logical_and(gen_img == 0, tar_img == 255)] = 255
    else:
        # saving rgb results
        rgb_img = np.ones(
            (gen_img.shape[0], gen_img.shape[1], 3), np.uint8) * 255
        # assign false negative as red channel
        rgb_img[:, :, 1][np.logical_and(gen_img == 255, tar_img == 0)] = 0
        rgb_img[:, :, 2][np.logical_and(gen_img == 255, tar_img == 0)] = 0
        # assign true positive as green channel
        rgb_img[:, :, 0][np.logical_and(gen_img == 255, tar_img == 255)] = 0
        rgb_img[:, :, 2][np.logical_and(gen_img == 255, tar_img == 255)] = 0
        # assign false positive as blue channel
        rgb_img[:, :, 0][np.logical_and(gen_img == 0, tar_img == 255)] = 0
        rgb_img[:, :, 1][np.logical_and(gen_img == 0, tar_img == 255)] = 0
    return rgb_img 

def trim_small_clusters(im, size=1):
    r"""
    Remove isolated voxels or clusters smaller than a given size

    Parameters
    ----------
    im : ND-array
        The binary image from which voxels are to be removed
    size : scalar
        The threshold size of clusters to trim.  As clusters with this many
        voxels or fewer will be trimmed.  The default is 1 so only single
        voxels are removed.

    Returns
    -------
    im : ND-image
        A copy of ``im`` with clusters of voxels smaller than the given
        ``size`` removed.

    """
    if im.dims == 2:
        strel = disk(1)
    elif im.ndims == 3:
        strel = ball(1)
    else:
        raise Exception('Only 2D or 3D images are accepted')
    filtered_array = np.copy(im)
    labels, N = spim.label(filtered_array, structure=strel)
    id_sizes = np.array(spim.sum(im, labels, range(N + 1)))
    area_mask = (id_sizes <= size)
    filtered_array[area_mask[labels]] = 0
    return filtered_array 

def compute_median(img, params):
    median_disk_size = int(params.get("median_disk_size", 3))
    return median(rgb2gray(img), selem=disk(median_disk_size))[:, :, None] 

def getIntensityThresholdOtsu(s, params):
    logging.info(f"{s['filename']} - \tLightDarkModule.getIntensityThresholdOtsu")
    name = "otsu"
    local = strtobool(params.get("local", "False"))
    radius = float(params.get("radius", 15))
    selem = disk(radius)

    img = s.getImgThumb(s["image_work_size"])
    img = color.rgb2gray(img)

    if local:
        thresh = rank.otsu(img, selem)
        name += "local"
    else:
        thresh = threshold_otsu(img)

    map = img < thresh

    s["img_mask_" + name] = map > 0
    if strtobool(params.get("invert", "False")):
        s["img_mask_" + name] = ~s["img_mask_" + name]

    io.imsave(s["outdir"] + os.sep + s["filename"] + "_" + name + ".png", img_as_ubyte(s["img_mask_" + name]))

    prev_mask = s["img_mask_use"]
    s["img_mask_use"] = s["img_mask_use"] & s["img_mask_" + name]

    s.addToPrintList(name,
                     printMaskHelper(params.get("mask_statistics", s["mask_statistics"]), prev_mask, s["img_mask_use"]))

    if len(s["img_mask_use"].nonzero()[0]) == 0:  # add warning in case the final tissue is empty
        logging.warning(f"{s['filename']} - After LightDarkModule.getIntensityThresholdOtsu:{name} NO tissue remains "
                        f"detectable! Downstream modules likely to be incorrect/fail")
        s["warnings"].append(f"After LightDarkModule.getIntensityThresholdOtsu:{name} NO tissue remains detectable! "
                             f"Downstream modules likely to be incorrect/fail")

    return 

def binary_dilation(x, radius=3):
    """ Return fast binary morphological dilation of an image.
    see `skimage.morphology.binary_dilation <http://scikit-image.org/docs/dev/api/skimage.morphology.html#skimage.morphology.binary_dilation>`_.

    Parameters
    -----------
    x : 2D array image.
    radius : int for the radius of mask.
    """
    from skimage.morphology import disk, binary_dilation
    mask = disk(radius)
    x = binary_dilation(image, selem=mask)
    return x 

def test_array_different_dimensions(self, nav_dims):
        shape = list(np.random.randint(2, 6, size=nav_dims))
        shape.extend([50, 50])
        chunks = [1] * nav_dims
        chunks.extend([25, 25])
        dask_array = da.random.random(size=shape, chunks=chunks)
        binary_image = sm.disk(5)
        match_array_dask = dt._template_match_with_binary_image(
            dask_array, binary_image=binary_image
        )
        assert len(dask_array.shape) == nav_dims + 2
        assert dask_array.shape == match_array_dask.shape
        match_array = match_array_dask.compute()
        assert dask_array.shape == match_array.shape 

def test_morphology_fft_dilate_2D(self):
        im = self.im[:, :, 50]
        truth = spim.binary_dilation(im, structure=disk(3))
        test = ps.tools.fftmorphology(im, strel=disk(3), mode='dilation')
        assert np.all(truth == test) 

def watershed_mask(
            self,
            stain_thresh: Number,
            markers: BinaryMaskCollection,
            disk_size: Optional[int],
    ) -> BinaryMaskCollection:
        """Create a watershed mask that is the union of the spot intensities above stain_thresh and
        a marker image generated from nuclei

        Parameters
        ----------
        stain_thresh : Number
            threshold to apply to the stain image
        markers : BinaryMaskCollection
            markers image generated from nuclei
        disk_size : Optional[int]
            if provided, execute a morphological opening operation over the thresholded stain image

        Returns
        -------
        BinaryMaskCollection :
            watershed mask
        """
        thresholded_stain = ThresholdBinarize(stain_thresh).run(self.stain)
        markers_and_stain = Merge.SimpleMerge().run([thresholded_stain, markers])
        watershed_mask = Filter.Reduce(
            "logical_or",
            lambda shape: np.zeros(shape=shape, dtype=np.bool)
        ).run(markers_and_stain)
        if disk_size is not None:
            disk_img = disk(disk_size)
            watershed_mask = Filter.Map(
                "morphology.binary_open",
                disk_img,
                module=FunctionSource.skimage
            ).run(watershed_mask)

        return watershed_mask 

def filter_nuclei(self, nuclei_thresh: float, disk_size: Optional[int]) -> BinaryMaskCollection:
        """Binarize the nuclei image using a thresholded binarizer and perform morphological binary
        opening.

        Parameters
        ----------
        nuclei_thresh : float
            Threshold the nuclei image at this value
        disk_size : int
            if passed, execute a binary opening of the filtered image

        Returns
        -------
        BinaryMaskCollection :
            mask collection with one mask, which is
        """
        nuclei_binarized = ThresholdBinarize(nuclei_thresh).run(self.nuclei_mp_scaled)
        if disk_size is not None:
            disk_img = disk(disk_size)
            nuclei_binarized = Filter.Map(
                "morphology.binary_open",
                disk_img,
                module=FunctionSource.skimage
            ).run(nuclei_binarized)

        # should only produce one binary mask.
        assert len(nuclei_binarized) == 1
        return nuclei_binarized 

def _white_tophat(self, image: xr.DataArray) -> xr.DataArray:
        if self.is_volume:
            structuring_element = ball(self.masking_radius)
        else:
            structuring_element = disk(self.masking_radius)
        return white_tophat(image, selem=structuring_element) 

def get_markers(m_b, c):
    # threshold
    c_thresh = threshold_otsu(c)
    c_b = c > c_thresh

    mk_ = np.where(c_b, 0, m_b)

    area, radius = mean_blob_size(m_b)
    struct_size = int(0.25 * radius)
    struct_el = morph.disk(struct_size)
    m_padded = pad_mask(mk_, pad=struct_size)
    m_padded = morph.erosion(m_padded, selem=struct_el)
    mk_ = crop_mask(m_padded, crop=struct_size)
    mk_, _ = ndi.label(mk_)
    return mk_ 

def clean_mask(m, c):
    # threshold
    m_thresh = threshold_otsu(m)
    c_thresh = threshold_otsu(c)
    m_b = m > m_thresh
    c_b = c > c_thresh

    # combine contours and masks and fill the cells
    m_ = np.where(m_b | c_b, 1, 0)
    m_ = ndi.binary_fill_holes(m_)

    # close what wasn't closed before
    area, radius = mean_blob_size(m_b)
    struct_size = int(1.25 * radius)
    struct_el = morph.disk(struct_size)
    m_padded = pad_mask(m_, pad=struct_size)
    m_padded = morph.binary_closing(m_padded, selem=struct_el)
    m_ = crop_mask(m_padded, crop=struct_size)

    # open to cut the real cells from the artifacts
    area, radius = mean_blob_size(m_b)
    struct_size = int(0.75 * radius)
    struct_el = morph.disk(struct_size)
    m_ = np.where(c_b & (~m_b), 0, m_)
    m_padded = pad_mask(m_, pad=struct_size)
    m_padded = morph.binary_opening(m_padded, selem=struct_el)
    m_ = crop_mask(m_padded, crop=struct_size)

    # join the connected cells with what we had at the beginning
    m_ = np.where(m_b | m_, 1, 0)
    m_ = ndi.binary_fill_holes(m_)

    # drop all the cells that weren't present at least in 25% of area in the initial mask
    m_ = drop_artifacts(m_, m_b, min_coverage=0.25)

    return m_ 

def _get_intensities(z, vectors, radius=1):
    """Basic intensity integration routine, takes the maximum value at the
    given vector positions with the number of pixels given by `radius`.

    Parameters
    ----------
    vectors : DiffractionVectors
        Vectors to the locations of the spots to be
        integrated.
    radius: int,
        Number of pixels within which to find the largest maximum

    Returns
    -------
    intensities : np.array
        List of extracted intensities
    """
    i, j = np.array(vectors.data).astype(int).T

    if radius > 1:
        footprint = morphology.disk(radius)
        filtered = ndi.maximum_filter(z, footprint=footprint)
        intensities = filtered[j, i].reshape(-1, 1)  # note that the indices are flipped
    else:
        intensities = z[j, i].reshape(-1, 1)  # note that the indices are flipped

    return np.array(intensities) 

def test_correct_disk_x_y_and_radius_random(self):
        x, y, px, py = 56, 48, 4, 5
        x, y = randint(45, 55, size=(py, px)), randint(45, 55, size=(py, px))
        r = randint(20, 40, size=(py, px))
        s = mdtd.generate_4d_data(
            probe_size_x=px,
            probe_size_y=py,
            image_size_x=120,
            image_size_y=100,
            disk_x=x,
            disk_y=y,
            disk_r=5,
            disk_I=20,
            ring_x=x,
            ring_y=y,
            ring_r=r,
            ring_I=5,
            blur=True,
            downscale=False,
        )
        dask_array = da.from_array(s.data, chunks=(4, 4, 50, 50))
        s = LazyDiffraction2D(dask_array)
        s_com = s.center_of_mass()
        s_r = s.radial_average(centre_x=s_com.inav[0].data, centre_y=s_com.inav[1].data)
        s_r = s_r.isig[15:]  # Do not include the disk
        r -= 15  # Need to shift the radius, due to not including the disk
        assert (s_r.data.argmax(axis=-1) == r).all() 

def test_1d_dask_array_error(self):
        binary_image = sm.disk(5)
        dask_array = da.random.random(size=50, chunks=10)
        with pytest.raises(ValueError):
            dt._template_match_with_binary_image(dask_array, binary_image=binary_image) 

def make_boundaries(label, thickness=None):
        """
        Input is an image label, output is a numpy array mask encoding the boundaries of the objects
        Extract pixels at the true boundary by dilation - erosion of label.
        Don't just pick the void label as it is not exclusive to the boundaries.
        """
        assert(thickness is not None)
        import skimage.morphology as skm
        void = 255
        mask = np.logical_and(label > 0, label != void)[0]
        selem = skm.disk(thickness)
        boundaries = np.logical_xor(skm.dilation(mask, selem),
                                    skm.erosion(mask, selem))
        return boundaries