Python skimage.morphology.binary_dilation() Examples

def _expand_raster(raster, distance = (4,2)):
    try:
        from skimage import draw, morphology
    except:
        raise ImportError("""
            The fill function requires the module "scikit-image"
            to operate.  Please retry after installing scikit-image:

            $ pip install --upgrade scikit-image """)
    if distance[0] <= 0.5 and distance[1] <= 0.5: return raster

    num_pixels = np.array(np.ceil(distance), dtype = int)
    neighborhood = np.zeros((num_pixels[1]*2+1, num_pixels[0]*2+1), dtype=np.bool)
    rr, cc = draw.ellipse(num_pixels[1], num_pixels[0], distance[1]+0.5, distance[0]+0.5)
    neighborhood[rr, cc] = 1

    return morphology.binary_dilation(image = raster, selem=neighborhood) 

def basin(label_mask, wall):
    h,w = np.shape(label_mask)
    y, x = np.mgrid[0:h, 0:w]
    struct = generate_binary_structure(2,2)
    shifty, shiftx = np.mgrid[0:3, 0:3]
    shifty = (shifty-1).flatten()
    shiftx = (shiftx-1).flatten()

    for i in range(4):
        obdr = label_mask^binary_dilation(label_mask, struct)
        ibdr = label_mask^binary_erosion(label_mask, struct)
        yob, xob = y[obdr], x[obdr]        
        ynb, xnb = yob.reshape(-1,1)+shifty, xob.reshape(-1,1)+shiftx
        
        wallnb = np.min(map_coords(wall, (ynb, xnb))*(map_coords(ibdr, (ynb, xnb))==1)+\
                        5*(map_coords(ibdr, (ynb, xnb))!=1),1)
        keep = (wall[yob,xob]>wallnb)&(wallnb<=4)
        label_mask[yob[keep], xob[keep]]=True
        if np.sum(keep)==0:
            break
    return label_mask 

def get_connected_component_shape(layer, event):
        cords = np.round(layer.coordinates).astype(int)
        val = layer.get_value()
        if val is None:
            return
        if val != 0:
            data = layer.data
            binary = data == val
            if 'Shift' in event.modifiers:
                binary_new = binary_erosion(binary)
                data[binary] = 0
                data[binary_new] = val
            else:
                binary_new = binary_dilation(binary)
                data[binary_new] = val
            size = np.sum(binary_new)
            layer.data = data
            msg = f'clicked at {cords} on blob {val} which is now {size} pixels'
        else:
            msg = f'clicked at {cords} on background which is ignored'
        layer.status = msg
        print(msg) 

def compute_binary_mask_sprengel(spectrogram, threshold):
    """ Computes a binary mask for the spectrogram
    # Arguments
        spectrogram : a numpy array representation of a spectrogram (2-dim)
        threshold   : a threshold for times larger than the median
    # Returns
        binary_mask : the binary mask
    """
    # normalize to [0, 1)
    norm_spectrogram = normalize(spectrogram)

    # median clipping
    binary_image = median_clipping(norm_spectrogram, threshold)

    # erosion
    binary_image = morphology.binary_erosion(binary_image, selem=np.ones((4, 4)))

    # dilation
    binary_image = morphology.binary_dilation(binary_image, selem=np.ones((4, 4)))

    # extract mask
    mask = np.array([np.max(col) for col in binary_image.T])
    mask = smooth_mask(mask)

    return mask 

def get_image_area_to_sample(img):
  """
  calculate set g_c, which has two properties
  1) They represent background pixels 
  2) They are within a certain distance to the object
  :param img: Image that represents the object instance
  """

  #TODO: In the paper 'Deep Interactive Object Selection', they calculate g_c first based on the original object instead
  # of the dilated one.

  # Dilate the object by d_margin pixels to extend the object boundary
  img_area = np.copy(img)
  img_area = morphology.binary_dilation(img_area, morphology.diamond(D_MARGIN)).astype(np.uint8)

  g_c = np.logical_not(img_area).astype(int)
  g_c[np.where(distance_transform_edt(g_c) > D)] = 0

  return g_c 

def getEdgeEnhancedWeightMap(label, label_ids =[0,1,2,3], scale=1, edgescale=1, assign_equal_wt=False):
    shape = (0,)+label.shape[1:]
    weight_map = np.empty(shape, dtype='uint8')
    if assign_equal_wt:
        return np.ones_like(label)
    for i in range(label.shape[0]): 
        #Estimate weight maps:
        weights = np.ones(len(label_ids))
        slice_map = np.ones(label[i,:,:].shape)
        for _id in label_ids:
            class_frequency = np.sum(label[i,:,:] == label_ids[_id])
            if class_frequency:
                weights[label_ids.index(_id)] = scale*label[i,:,:].size/class_frequency
                slice_map[np.where(label[i,:,:]==label_ids.index(_id))] = weights[label_ids.index(_id)]
                edge = np.float32(morph.binary_dilation(
                    canny(np.float32(label[i,:,:]==label_ids.index(_id)),sigma=1), selem=selem))
                edge_frequency = np.sum(np.sum(edge==1.0))
                if edge_frequency:    
                    slice_map[np.where(edge==1.0)] += edgescale*label[i,:,:].size/edge_frequency
            # print (weights)
            # utils.imshow(edge, cmap='gray')
        # utils.imshow(weight_map, cmap='gray')
        weight_map = np.append(weight_map, np.expand_dims(slice_map, axis=0), axis=0)
    return np.float32(weight_map) 

def test_apply():
    input_mask_collection = binary_mask_collection_2d()
    output_mask_collection = input_mask_collection._apply(binary_dilation)

    assert input_mask_collection._pixel_ticks == output_mask_collection._pixel_ticks
    assert input_mask_collection._physical_ticks == output_mask_collection._physical_ticks
    assert input_mask_collection._log == output_mask_collection._log
    assert len(input_mask_collection) == len(output_mask_collection)

    region_0, region_1 = output_mask_collection.masks()

    assert region_0.name == '0'
    assert region_1.name == '1'

    temp = np.ones((2, 6), dtype=np.bool)
    assert np.array_equal(region_0, temp)
    temp = np.ones((3, 4), dtype=np.bool)
    temp[0, 0] = 0
    assert np.array_equal(region_1, temp)

    assert np.array_equal(region_0[Axes.Y.value], [0, 1])
    assert np.array_equal(region_0[Axes.X.value], [0, 1, 2, 3, 4, 5])

    assert np.array_equal(region_1[Axes.Y.value], [2, 3, 4])
    assert np.array_equal(region_1[Axes.X.value], [2, 3, 4, 5]) 

def _apply_dilation(self, snow_masks):
        """ Apply binary dilation for each mask in the series """
        if self.dilation_size:
            snow_masks = np.array([binary_dilation(mask, disk(self.dilation_size)) for mask in snow_masks])
        return snow_masks 

def _dilate_mask(mask, dilation_radius=5):
  """Damages a mask for a single object by dilating it in numpy.

  Args:
    mask: Boolean numpy array of shape(height, width, 1).
    dilation_radius: Integer, the radius of the used disk structure element.

  Returns:
    The dilated version of mask.
  """
  disk = morphology.disk(dilation_radius, dtype=np.bool)
  dilated_mask = morphology.binary_dilation(
      np.squeeze(mask, axis=2), selem=disk)[..., np.newaxis]
  return dilated_mask 

def modify_w_unet(base_label, preds, thres=0.25):
    base_label = base_label.copy()
    vals = np.unique(base_label[base_label>0])
    struct = generate_binary_structure(2,2)
    for nb_dilation in range(3):    
        for val in vals:
            label_mask = base_label==val
            base_label[binary_dilation(label_mask,struct)&(preds[:,:,0]>thres)&(base_label==0)]=val
    return base_label 

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 get_simple_eroded_dilated_mask(mask, erode_selem_size, dilate_selem_size, small_annotations_size):
    if mask.sum() > small_annotations_size**2:
        selem = rectangle(erode_selem_size, erode_selem_size)
        mask_ = binary_erosion(mask, selem=selem)
    else:
        selem = rectangle(dilate_selem_size, dilate_selem_size)
        mask_ = binary_dilation(mask, selem=selem)
    return mask_ 

def grow_mask(anat, aseg, ants_segs=None, ww=7, zval=2.0, bw=4):
    """
    Grow mask including pixels that have a high likelihood.

    GM tissue parameters are sampled in image patches of ``ww`` size.
    This is inspired on mindboggle's solution to the problem:
    https://github.com/nipy/mindboggle/blob/master/mindboggle/guts/segment.py#L1660

    """
    selem = sim.ball(bw)

    if ants_segs is None:
        ants_segs = np.zeros_like(aseg, dtype=np.uint8)

    aseg[aseg == 42] = 3  # Collapse both hemispheres
    gm = anat.copy()
    gm[aseg != 3] = 0

    refined = refine_aseg(aseg)
    newrefmask = sim.binary_dilation(refined, selem) - refined
    indices = np.argwhere(newrefmask > 0)
    for pixel in indices:
        # When ATROPOS identified the pixel as GM, set and carry on
        if ants_segs[tuple(pixel)] == 2:
            refined[tuple(pixel)] = 1
            continue

        window = gm[
            pixel[0] - ww:pixel[0] + ww,
            pixel[1] - ww:pixel[1] + ww,
            pixel[2] - ww:pixel[2] + ww,
        ]
        if np.any(window > 0):
            mu = window[window > 0].mean()
            sigma = max(window[window > 0].std(), 1.0e-5)
            zstat = abs(anat[tuple(pixel)] - mu) / sigma
            refined[tuple(pixel)] = int(zstat < zval)

    refined = sim.binary_opening(refined, selem)
    return refined 

def create_mask(aseg_data, dnum, enum):
    from skimage.morphology import binary_dilation, binary_erosion
    from skimage.measure import label
    print ("Creating dilated mask ...")

    # treat lateral orbital frontal and parsorbitalis special to avoid capturing too much of eye nerve
    lat_orb_front_mask = np.logical_or(aseg_data == 2012, aseg_data == 1012)
    parsorbitalis_mask = np.logical_or(aseg_data == 2019, aseg_data == 1019)
    frontal_mask = np.logical_or(lat_orb_front_mask, parsorbitalis_mask)
    print("Frontal region special treatment: ", format(np.sum(frontal_mask)))

    # reduce to binary
    datab = (aseg_data > 0)
    datab[frontal_mask] = 0
    # dilate and erode
    for x in range(dnum):
        datab = binary_dilation(datab, np.ones((3, 3, 3)))
    for x in range(enum):
        datab = binary_erosion(datab, np.ones((3, 3, 3)))

    # extract largest component
    labels = label(datab)
    assert (labels.max() != 0)  # assume at least 1 real connected component
    print("  Found {} connected component(s)!".format(labels.max()))

    if labels.max() > 1:
        print("  Selecting largest component!")
        datab = (labels == np.argmax(np.bincount(labels.flat)[1:]) + 1)

    # add frontal regions back to mask
    datab[frontal_mask] = 1

    # set mask
    aseg_data[~datab] = 0
    aseg_data[datab] = 1
    return aseg_data 

def sprengel_binary_mask_from_wave_file(filepath):
    fs, x = utils.read_wave_file(filepath)
    Sxx = sp.wave_to_amplitude_spectrogram(x, fs)
    Sxx_log = sp.wave_to_log_amplitude_spectrogram(x, fs)

    # plot spectrogram
    fig = plt.figure(1)
    subplot_image(Sxx_log, 411, "Spectrogram")

    Sxx = pp.normalize(Sxx)
    binary_image = pp.median_clipping(Sxx, 3.0)

    subplot_image(binary_image + 0, 412, "Median Clipping")

    binary_image = morphology.binary_erosion(binary_image, selem=np.ones((4, 4)))

    subplot_image(binary_image + 0, 413, "Erosion")

    binary_image = morphology.binary_dilation(binary_image, selem=np.ones((4, 4)))

    subplot_image(binary_image + 0, 414, "Dilation")

    mask = np.array([np.max(col) for col in binary_image.T])
    mask = morphology.binary_dilation(mask, np.ones(4))
    mask = morphology.binary_dilation(mask, np.ones(4))

    # plot_vector(mask, "Mask")

    fig.set_size_inches(10, 12)
    plt.tight_layout()
    fig.savefig(utils.get_basename_without_ext(filepath) + "_binary_mask.png", dpi=100) 

def smooth_mask(mask):
    """ Smooths a binary mask using 4x4 dilation
        # Arguments
            mask : the binary mask
        # Returns
            mask : a smoother binary mask
    """
    n_hood = np.ones(4)
    mask = morphology.binary_dilation(mask, n_hood)
    mask = morphology.binary_dilation(mask, n_hood)

    # type casting is a bitch
    return mask 

def compute_binary_mask_lasseck(spectrogram, threshold):
    # normalize to [0, 1)
    norm_spectrogram = normalize(spectrogram)

    # median clipping
    binary_image = median_clipping(norm_spectrogram, threshold)

    # closing binary image (dilation followed by erosion)
    binary_image = morphology.binary_closing(binary_image, selem=np.ones((4, 4)))

    # dialate binary image
    binary_image = morphology.binary_dilation(binary_image, selem=np.ones((4, 4)))

    # apply median filter
    binary_image = filters.median(binary_image, selem=np.ones((2, 2)))

    # remove small objects
    binary_image = morphology.remove_small_objects(binary_image, min_size=32, connectivity=1)

    mask = np.array([np.max(col) for col in binary_image.T])
    mask = smooth_mask(mask)

    return mask


# TODO: This method needs some real testing 

def outline_polygons(self, width=EDGE_WIDTH, color=LABEL_EDGE):
        from skimage.morphology import binary_dilation, disk
        im = np.asarray(self.image).copy()
        outset = binary_dilation(im == LABEL_POSITIVE, disk(width / 2))
        inset = binary_dilation(im != LABEL_POSITIVE, disk(width - width / 2))
        boundary = outset & inset
        im[boundary] = color
        self.image = Image.fromarray(im)
        self.artist = ImageDraw.Draw(self.image) 

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 create_shoreline_buffer(im_shape, georef, image_epsg, pixel_size, settings):
    """
    Creates a buffer around the reference shoreline. The size of the buffer is 
    given by settings['max_dist_ref'].

    KV WRL 2018

    Arguments:
    -----------
    im_shape: np.array
        size of the image (rows,columns)
    georef: np.array
        vector of 6 elements [Xtr, Xscale, Xshear, Ytr, Yshear, Yscale]
    image_epsg: int
        spatial reference system of the image from which the contours were extracted
    pixel_size: int
        size of the pixel in metres (15 for Landsat, 10 for Sentinel-2)
    settings: dict with the following keys
        'output_epsg': int
            output spatial reference system
        'reference_shoreline': np.array
            coordinates of the reference shoreline
        'max_dist_ref': int
            maximum distance from the reference shoreline in metres

    Returns:    
    -----------
    im_buffer: np.array
        binary image, True where the buffer is, False otherwise

    """
    # initialise the image buffer
    im_buffer = np.ones(im_shape).astype(bool)

    if 'reference_shoreline' in settings.keys():

        # convert reference shoreline to pixel coordinates
        ref_sl = settings['reference_shoreline']
        ref_sl_conv = SDS_tools.convert_epsg(ref_sl, settings['output_epsg'],image_epsg)[:,:-1]
        ref_sl_pix = SDS_tools.convert_world2pix(ref_sl_conv, georef)
        ref_sl_pix_rounded = np.round(ref_sl_pix).astype(int)

        # make sure that the pixel coordinates of the reference shoreline are inside the image
        idx_row = np.logical_and(ref_sl_pix_rounded[:,0] > 0, ref_sl_pix_rounded[:,0] < im_shape[1])
        idx_col = np.logical_and(ref_sl_pix_rounded[:,1] > 0, ref_sl_pix_rounded[:,1] < im_shape[0])
        idx_inside = np.logical_and(idx_row, idx_col)
        ref_sl_pix_rounded = ref_sl_pix_rounded[idx_inside,:]

        # create binary image of the reference shoreline (1 where the shoreline is 0 otherwise)
        im_binary = np.zeros(im_shape)
        for j in range(len(ref_sl_pix_rounded)):
            im_binary[ref_sl_pix_rounded[j,1], ref_sl_pix_rounded[j,0]] = 1
        im_binary = im_binary.astype(bool)

        # dilate the binary image to create a buffer around the reference shoreline
        max_dist_ref_pixels = np.ceil(settings['max_dist_ref']/pixel_size)
        se = morphology.disk(max_dist_ref_pixels)
        im_buffer = morphology.binary_dilation(im_binary, se)

    return im_buffer 

def davis_f_measure(foreground_mask, gt_mask, bound_th=0.008):
    """
    Compute mean,recall and decay from per-frame evaluation.
    Calculates precision/recall for boundaries between foreground_mask and
    gt_mask using morphological operators to speed it up.

    Arguments:
        foreground_mask (ndarray): binary segmentation image.
        gt_mask         (ndarray): binary annotated image.

    Returns:
        F (float): boundaries F-measure
        P (float): boundaries precision
        R (float): boundaries recall
    """
    assert np.atleast_3d(foreground_mask).shape[2] == 1

    bound_pix = bound_th if bound_th >= 1 else \
        np.ceil(bound_th * np.linalg.norm(foreground_mask.shape))

    # Get the pixel boundaries of both masks
    fg_boundary = seg2bmap(foreground_mask)
    gt_boundary = seg2bmap(gt_mask)

    fg_dil = binary_dilation(fg_boundary, disk(bound_pix))
    gt_dil = binary_dilation(gt_boundary, disk(bound_pix))

    # Get the intersection
    gt_match = gt_boundary * fg_dil
    fg_match = fg_boundary * gt_dil

    # Area of the intersection
    n_fg = np.sum(fg_boundary)
    n_gt = np.sum(gt_boundary)

    # % Compute precision and recall
    if n_fg == 0 and n_gt > 0:
        precision = 1
        recall = 0
    elif n_fg > 0 and n_gt == 0:
        precision = 0
        recall = 1
    elif n_fg == 0 and n_gt == 0:
        precision = 1
        recall = 1
    else:
        precision = np.sum(fg_match) / float(n_fg)
        recall = np.sum(gt_match) / float(n_gt)

    # Compute F measure
    if precision + recall == 0:
        F = 0
    else:
        F = 2 * precision * recall / (precision + recall)

    return F 

def merge_sinks(Label, Sinks, Radius=5):
    """
    Merges attraction basins obtained from gradient flow tracking using
    sink locations.

    Parameters
    ----------
    Segmentation : array_like
        Label image where positive values correspond to foreground pixels that
        share mutual sinks.
    Sinks : array_like
        N x 2 array containing the (x,y) locations of the tracking sinks. Each
        row is an (x,y) pair - in that order.
    Radius : float
        Radius used to merge sinks. Sinks closer than this radius to one
        another will have their regions of attraction merged.
        Default value = 5.

    Returns
    -------
    Merged : array_like
        Label image where attraction regions are merged.

    """

    # build seed image
    SeedImage = np.zeros(Label.shape)
    for i in range(Sinks.shape[0]):
        SeedImage[Sinks[i, 1], Sinks[i, 0]] = i+1

    # dilate sink image
    Dilated = mp.binary_dilation(SeedImage, mp.disk(Radius))

    # generate new labels for merged seeds, define memberships
    Labels = ms.label(Dilated)
    New = Labels[Sinks[:, 1].astype(np.int), Sinks[:, 0].astype(np.int)]

    # get unique list of seed clusters
    Unique = np.arange(1, New.max()+1)

    # generate new seed list
    Merged = np.zeros(Label.shape)

    # get pixel list for each sink object
    Props = ms.regionprops(Label.astype(np.int))

    # fill in new values
    for i in Unique:
        Indices = np.nonzero(New == i)[0]
        for j in Indices:
            Coords = Props[j].coords
            Merged[Coords[:, 0], Coords[:, 1]] = i

    return Merged 

def db_eval_boundary(foreground_mask, gt_mask, bound_th):
    """
    Compute mean,recall and decay from per-frame evaluation.
    Calculates precision/recall for boundaries between foreground_mask and
    gt_mask using morphological operators to speed it up.
    Arguments:
        foreground_mask (ndarray): binary segmentation image.
        gt_mask         (ndarray): binary annotated image.
    Returns:
        F (float): boundaries F-measure
        P (float): boundaries precision
        R (float): boundaries recall
    """
    assert np.atleast_3d(foreground_mask).shape[2] == 1

    bound_pix = bound_th if bound_th >= 1 else \
        np.ceil(bound_th * np.linalg.norm(foreground_mask.shape))

    # Get the pixel boundaries of both masks
    fg_boundary = seg2bmap(foreground_mask);
    gt_boundary = seg2bmap(gt_mask);

    fg_dil = binary_dilation(fg_boundary, disk(bound_pix))
    gt_dil = binary_dilation(gt_boundary, disk(bound_pix))

    # Get the intersection
    gt_match = gt_boundary * fg_dil
    fg_match = fg_boundary * gt_dil

    # Area of the intersection
    n_fg = np.sum(fg_boundary)
    n_gt = np.sum(gt_boundary)

    # % Compute precision and recall
    if n_fg == 0 and n_gt > 0:
        precision = 1
        recall = 0
    elif n_fg > 0 and n_gt == 0:
        precision = 0
        recall = 1
    elif n_fg == 0 and n_gt == 0:
        precision = 1
        recall = 1
    else:
        precision = np.sum(fg_match) / float(n_fg)
        recall = np.sum(gt_match) / float(n_gt)

    # Compute F measure
    if precision + recall == 0:
        F = 0
    else:
        F = 2 * precision * recall / (precision + recall);

    return F, precision, recall, np.sum(fg_match), n_fg, np.sum(gt_match), n_gt 

def overlay_skeleton_2d(image, skeleton, *,
                        image_cmap=None, color=(1, 0, 0), alpha=1,
                        dilate=0, axes=None):
    """Overlay the skeleton pixels on the input image.

    Parameters
    ----------
    image : array, shape (M, N[, 3])
        The input image. Can be grayscale or RGB.
    skeleton : array, shape (M, N)
        The input 1-pixel-wide skeleton.

    Other Parameters
    ----------------
    image_cmap : matplotlib colormap name or object, optional
        If the input image is grayscale, colormap it with this colormap.
        The default is grayscale.
    color : tuple of float in [0, 1], optional
        The RGB color for the skeleton pixels.
    alpha : float, optional
        Blend the skeleton pixels with the given alpha.
    dilate : int, optional
        Dilate the skeleton by this amount. This is useful when rendering
        large images where aliasing may cause some pixels of the skeleton
        not to be drawn.
    axes : matplotlib Axes
        The Axes on which to plot the image. If None, new ones are created.

    Returns
    -------
    axes : matplotlib Axes
        The Axis on which the image is drawn.
    """
    image = _normalise_image(image, image_cmap=image_cmap)
    skeleton = skeleton.astype(bool)
    if dilate > 0:
        selem = morphology.disk(dilate)
        skeleton = morphology.binary_dilation(skeleton, selem)
    if axes is None:
        fig, axes = plt.subplots()
    image[skeleton] = alpha * np.array(color) + (1 - alpha) * image[skeleton]
    axes.imshow(image)
    axes.axis('off')
    return axes 

def segment_active_contour(img, centers):
    """ segmentation using acive contours

    :param ndarray img: input image / segmentation
    :param [[int, int]] centers: position of centres / seeds
    :return (ndarray, [[int, int]]): resulting segmentation, updated centres
    """
    logging.debug('segment: active_contour...')
    # http://scikit-image.org/docs/dev/auto_examples/edges/plot_active_contours.html
    segm = np.zeros(img.shape[:2])
    img_smooth = ndimage.filters.gaussian_filter(img, 5)
    center_circles, _, _ = create_circle_center(img.shape[:2], centers)
    for i, snake in enumerate(center_circles):
        snake = segmentation.active_contour(img_smooth, snake.astype(float),
                                            alpha=0.015, beta=10, gamma=0.001,
                                            w_line=0.0, w_edge=1.0,
                                            max_px_move=1.0,
                                            max_iterations=2500,
                                            convergence=0.2)
        seg = np.zeros(segm.shape, dtype=bool)
        x, y = np.array(snake).transpose().tolist()
        # rr, cc = draw.polygon(x, y)
        seg[map(int, x), map(int, y)] = True
        seg = morphology.binary_dilation(seg, selem=morphology.disk(3))
        bb_area = int((max(x) - min(x)) * (max(y) - min(y)))
        logging.debug('bounding box area: %d', bb_area)
        seg = morphology.remove_small_holes(seg, min_size=bb_area)
        segm[seg] = i + 1
    return segm, centers, None