Python skimage.morphology.watershed() Examples

def skimage_watershed_segmentation(mask, kernel=k_3x3, k=1):
    # mask = cv.dilate(mask, kernel, iterations=k)

    distance = ndimage.distance_transform_edt(mask)
    local_maxi = peak_local_max(distance, indices=False, footprint=kernel, labels=mask)

    markers = measure.label(local_maxi)
    labels_ws = watershed(-distance, markers, mask=mask)

    if labels_ws.max() < 2:
        return [mask], labels_ws

    res_masks = []
    for idx in range(1,  labels_ws.max() + 1):
        m = labels_ws == idx
        if m.sum() > 20:
            res_masks.append(m.astype(np.uint8))
    return res_masks, labels_ws 

def DynamicWatershedAlias(p_img, lamb, p_thresh = 0.5):
    """
    Applies our dynamic watershed to 2D prob/dist image.
    """
    b_img = (p_img > p_thresh) + 0
    Probs_inv = PrepareProb(p_img)


    Hrecons = HreconstructionErosion(Probs_inv, lamb)
    markers_Probs_inv = find_maxima(Hrecons, mask = b_img)
    markers_Probs_inv = label(markers_Probs_inv)
    ws_labels = watershed(Hrecons, markers_Probs_inv, mask=b_img)
    arrange_label = ArrangeLabel(ws_labels)
    wsl = generate_wsl(arrange_label)
    arrange_label[wsl > 0] = 0
    

    return arrange_label 

def overlay_cut_masks(images_dir, subdir_name, target_dir, cut_size=1):
    train_dir = os.path.join(images_dir, subdir_name)
    for mask_dirname in tqdm(glob.glob('{}/*/masks'.format(train_dir))):
        masks = []
        for ind, image_filepath in enumerate(glob.glob('{}/*'.format(mask_dirname))):
            image = np.asarray(Image.open(image_filepath))
            image = np.where(image > 0, ind + 1, 0)
            masks.append(image)
        labeled_masks = np.sum(masks, axis=0)
        overlayed_masks = np.where(labeled_masks, 1, 0)

        watershed_mask = watershed(overlayed_masks.astype(np.bool), labeled_masks, watershed_line=True)
        if watershed_mask.max() == watershed_mask.min():
            cut_masks = overlayed_masks
        else:
            borders = (watershed_mask == 0) & overlayed_masks
            selem = rectangle(cut_size, cut_size)
            dilated_borders = dilation(borders, selem=selem)
            cut_masks = np.where(dilated_borders, 0, overlayed_masks)

        target_filepath = '/'.join(mask_dirname.replace(images_dir, target_dir).split('/')[:-1]) + '.png'
        os.makedirs(os.path.dirname(target_filepath), exist_ok=True)
        imwrite(target_filepath, cut_masks) 

def _watershed(self, inraster, th_tree=15.):
        """ Simple implementation of a watershed tree crown delineation

        Parameters
        ----------
        inraster :   ndarray
                     raster of height values (e.g., CHM)
        th_tree :    float
                     minimum height of tree crown

        Returns
        -------
        ndarray
            raster of individual tree crowns
        """
        inraster_mask = inraster.copy()
        inraster_mask[inraster <= th_tree] = 0
        raster = inraster.copy()
        raster[np.isnan(raster)] = 0.
        labels = watershed(-raster, self.tree_markers, mask=inraster_mask)
        return labels 

def postprocess_victor(pred):
    av_pred = pred / 255.
    av_pred = av_pred[..., 2] * (1 - av_pred[..., 1])
    av_pred = 1 * (av_pred > 0.5)
    av_pred = av_pred.astype(np.uint8)

    y_pred = measure.label(av_pred, neighbors=8, background=0)
    props = measure.regionprops(y_pred)
    for i in range(len(props)):
        if props[i].area < 12:
            y_pred[y_pred == i + 1] = 0
    y_pred = measure.label(y_pred, neighbors=8, background=0)

    nucl_msk = (255 - pred[..., 2])
    nucl_msk = nucl_msk.astype('uint8')
    y_pred = watershed(nucl_msk, y_pred, mask=((pred[..., 2] > 80)), watershed_line=True)
    return y_pred

# test_dir = r'C:\dev\dsbowl\results_test\bowl_remap3\merged'
# borders_dir = r'C:\dev\dsbowl\results_test\bowl_remap_border2\merged' 

def label_mask(pred, main_threshold=0.3, seed_threshold=0.7, w_pixel_t=20, pixel_t=100):
    av_pred = pred / 255.
    av_pred = av_pred[..., 0] * (1 - av_pred[..., 2])
    av_pred = 1 * (av_pred > seed_threshold)
    av_pred = av_pred.astype(np.uint8)

    y_pred = measure.label(av_pred, neighbors=8, background=0)
    props = measure.regionprops(y_pred)
    for i in range(len(props)):
        if props[i].area < w_pixel_t:
            y_pred[y_pred == i + 1] = 0
    y_pred = measure.label(y_pred, neighbors=8, background=0)

    nucl_msk = (255 - pred[..., 0])
    nucl_msk = nucl_msk.astype('uint8')
    y_pred = watershed(nucl_msk, y_pred, mask=(pred[..., 0] > main_threshold * 255), watershed_line=True)

    props = measure.regionprops(y_pred)

    for i in range(len(props)):
        if props[i].area < pixel_t:
            y_pred[y_pred == i + 1] = 0
    y_pred = measure.label(y_pred, neighbors=8, background=0)
    return y_pred 

def watersplit(_probs, _points):
    points = _points.copy()

    points[points != 0] = np.arange(1, points.sum()+1)
    points = points.astype(float)

    probs = ndimage.black_tophat(_probs.copy(), 7)
    seg = watershed(probs, points)

    return find_boundaries(seg) 

def watershed(image, label=None):
    denoised = filters.rank.median(image, morphology.disk(2)) #过滤噪声
    #将梯度值低于10的作为开始标记点
    markers = filters.rank.gradient(denoised, morphology.disk(5)) < 10
    markers = ndi.label(markers)[0]

    gradient = filters.rank.gradient(denoised, morphology.disk(2)) #计算梯度
    labels =morphology.watershed(gradient, markers, mask=image) #基于梯度的分水岭算法

    fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(6, 6))
    axes = axes.ravel()
    ax0, ax1, ax2, ax3 = axes

    ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
    ax0.set_title("Original")
    # ax1.imshow(gradient, cmap=plt.cm.spectral, interpolation='nearest')
    ax1.imshow(gradient, cmap=plt.cm.gray, interpolation='nearest')
    ax1.set_title("Gradient")
    if label is not None:
        # ax2.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest')
        ax2.imshow(label, cmap=plt.cm.gray, interpolation='nearest')
    else:
        ax2.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest')
    ax2.set_title("Markers")
    ax3.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest')
    ax3.set_title("Segmented")

    for ax in axes:
        ax.axis('off')

    fig.tight_layout()
    plt.show() 

def breakup_region(component):
    distance = ndi.distance_transform_edt(component)
    skel = skeletonize(component)
    skeldist = distance*skel
    local_maxi = peak_local_max(skeldist, indices=False, footprint=disk(10))
    local_maxi=ndi.binary_closing(local_maxi,structure = disk(4),iterations = 2)
    markers = ndi.label(local_maxi)[0]
    labels = watershed(-distance, markers, mask=component)
    return(labels) 

def opencv_segmentation(mask, kernel=k_3x3, k=3):
    # noise removal
    opening = cv.morphologyEx(mask, cv.MORPH_OPEN, kernel, iterations=k)

    # sure background area
    sure_bg = cv.dilate(opening, kernel, iterations=k)

    # Finding sure foreground area
    dist_transform = cv.distanceTransform(opening,cv.DIST_L2, 5)
    ret, sure_fg = cv.threshold(dist_transform, 0.7*dist_transform.max(), 255, 0)

    # Finding unknown region
    sure_fg = np.uint8(sure_fg)
    unknown = cv.subtract(sure_bg, sure_fg)

    # Marker labelling
    ret, markers = cv.connectedComponents(sure_fg)

    # Add one to all labels so that sure background is not 0, but 1
    markers = markers + 1

    # Now, mark the region of unknown with zero
    markers[unknown > 0] = 0

    labels_ws = cv.watershed(cv.cvtColor(mask, cv.COLOR_GRAY2RGB), markers)

    if labels_ws.max() - 1 < 2:
        return [mask], labels_ws

    res_masks = []
    for idx in range(2,  labels_ws.max() + 1):
        m = labels_ws == idx
        if m.sum() > 5:
            m = cv.dilate(m.astype(np.uint8), kernel, iterations=1)
            res_masks.append(m)
    return res_masks, labels_ws 

def my_watershed(what, mask1, mask2):
    # markers = ndi.label(mask2, output=np.uint32)[0]
    # big_seeds = watershed(what, markers, mask=mask1, watershed_line=False)
    # m2 = mask1 - (big_seeds > 0)
    # mask2 = mask2 | m2

    markers = ndi.label(mask2, output=np.uint32)[0]
    labels = watershed(what, markers, mask=mask1, watershed_line=True)
    # labels = watershed(what, markers, mask=mask1, watershed_line=False)
    return labels 

def watershedding_2D(track,field_in,**kwargs):
    kwargs.pop('method',None)
    return segmentation_2D(track,field_in,method='watershed',**kwargs) 

def watershedding_3D(track,field_in,**kwargs):
    kwargs.pop('method',None)
    return segmentation_3D(track,field_in,method='watershed',**kwargs) 

def segmentation_2D(features,field,dxy,threshold=3e-3,target='maximum',level=None,method='watershed',max_distance=None):
    return segmentation(features,field,dxy,threshold=threshold,target=target,level=level,method=method,max_distance=max_distance) 

def my_watershed(what, mask1, mask2):
    markers = ndi.label(mask2, output=np.uint32)[0]
    labels = watershed(what, markers, mask=mask1, watershed_line=True)
    return labels 

def segmentation_3D(features,field,dxy,threshold=3e-3,target='maximum',level=None,method='watershed',max_distance=None):
    return segmentation(features,field,dxy,threshold=threshold,target=target,level=level,method=method,max_distance=max_distance) 

def create_separation(labels):
    tmp = dilation(labels > 0, square(12))
    tmp2 = watershed(tmp, labels, mask=tmp, watershed_line=True) > 0
    tmp = tmp ^ tmp2
    tmp = dilation(tmp, square(3))

    props = measure.regionprops(labels)

    msk1 = np.zeros_like(labels, dtype='bool')

    for y0 in range(labels.shape[0]):
        for x0 in range(labels.shape[1]):
            if not tmp[y0, x0]:
                continue
            if labels[y0, x0] == 0:
                sz = 5
            else:
                sz = 7
                if props[labels[y0, x0] - 1].area < 300:
                    sz = 5
                elif props[labels[y0, x0] - 1].area < 2000:
                    sz = 6
            uniq = np.unique(labels[max(0, y0 - sz):min(labels.shape[0], y0 + sz + 1),
                             max(0, x0 - sz):min(labels.shape[1], x0 + sz + 1)])
            if len(uniq[uniq > 0]) > 1:
                msk1[y0, x0] = True
    return msk1 

def proc_np_dist(pred):
    """
    Process Nuclei Prediction with Distance Map

    Args:
        pred: prediction output, assuming 
                channel 0 contain probability map of nuclei
                channel 1 containing the regressed distance map
    """
    blb_raw = pred[...,0]
    dst_raw = pred[...,1]

    blb = np.copy(blb_raw)
    blb[blb >  0.5] = 1
    blb[blb <= 0.5] = 0
    blb = measurements.label(blb)[0]
    blb = remove_small_objects(blb, min_size=10)
    blb[blb > 0] = 1   

    dst_raw[dst_raw < 0] = 0
    dst = np.copy(dst_raw)
    dst = dst * blb
    dst[dst  > 0.5] = 1
    dst[dst <= 0.5] = 0

    marker = dst.copy()
    marker = binary_fill_holes(marker) 
    marker = measurements.label(marker)[0]
    marker = remove_small_objects(marker, min_size=10)
    proced_pred = watershed(-dst_raw, marker, mask=blb)    
    return proced_pred

#### 

def overlay_masks_with_borders_json(images_dir, subdir_name, target_dir, borders_size=3, dilation_size=5):
    train_dir = os.path.join(images_dir, subdir_name)
    for mask_dirname in tqdm(glob.glob('{}/*/masks'.format(train_dir))):
        masks = []
        for ind, image_filepath in enumerate(glob.glob('{}/*'.format(mask_dirname))):
            image = np.asarray(Image.open(image_filepath))
            image = np.where(image > 0, ind + 1, 0)
            masks.append(image)
        labeled_masks = np.sum(masks, axis=0)
        overlayed_masks = np.where(labeled_masks, 1, 0)

        selem = rectangle(dilation_size, dilation_size)
        dilated_mask = dilation(overlayed_masks, selem=selem)
        watershed_mask = watershed((dilated_mask >= 0).astype(np.bool), labeled_masks, watershed_line=True)

        if watershed_mask.max() == watershed_mask.min():
            dilated_borders = np.zeros_like(overlayed_masks)
        else:
            borders = (watershed_mask == 0) & (dilated_mask > 0)
            selem = rectangle(borders_size, borders_size)
            dilated_borders = dilation(borders, selem=selem)

        nuclei = prepare_class_encoding(overlayed_masks)
        borders = prepare_class_encoding(dilated_borders)

        target_filepath = '/'.join(mask_dirname.replace(images_dir, target_dir).split('/')[:-1]) + '.json'
        os.makedirs(os.path.dirname(target_filepath), exist_ok=True)
        save_target_masks(target_filepath, nuclei, borders) 

def generate_wsl(ws):
    """
    Generates watershed line that correspond to areas of
    touching objects.
    """
    se = square(3)
    ero = ws.copy()
    ero[ero == 0] = ero.max() + 1
    ero = erosion(ero, se)
    ero[ws == 0] = 0

    grad = dilation(ws, se) - ero
    grad[ws == 0] = 0
    grad[grad > 0] = 255
    grad = grad.astype(np.uint8)
    return grad 

def watershed(masks, seeds, borders):
    seeds_detached = seeds * (1 - borders)
    markers = label(seeds_detached)
    labels = morph.watershed(masks, markers, mask=masks)
    return labels 

def overlay_masks_with_borders(images_dir, subdir_name, target_dir, borders_size=3, dilation_size=5):
    train_dir = os.path.join(images_dir, subdir_name)
    for mask_dirname in tqdm(glob.glob('{}/*/masks'.format(train_dir))):
        masks = []
        for ind, image_filepath in enumerate(glob.glob('{}/*'.format(mask_dirname))):
            image = np.asarray(Image.open(image_filepath))
            image = np.where(image > 0, ind + 1, 0)
            masks.append(image)
        labeled_masks = np.sum(masks, axis=0)
        overlayed_masks = np.where(labeled_masks, 1, 0)

        selem = rectangle(dilation_size, dilation_size)
        dilated_mask = dilation(overlayed_masks, selem=selem)
        watershed_mask = watershed((dilated_mask >= 0).astype(np.bool), labeled_masks, watershed_line=True)

        if watershed_mask.max() == watershed_mask.min():
            masks_with_borders = overlayed_masks
        else:
            borders = (watershed_mask == 0) & (dilated_mask > 0)
            selem = rectangle(borders_size, borders_size)
            dilated_borders = dilation(borders, selem=selem)
            masks_with_borders = np.where(dilated_borders, 2, overlayed_masks)

        target_filepath = '/'.join(mask_dirname.replace(images_dir, target_dir).split('/')[:-1]) + '.png'
        os.makedirs(os.path.dirname(target_filepath), exist_ok=True)
        imwrite(target_filepath, masks_with_borders) 

def detectCellShape(img, peaks, detectCellShapeParameter = None, threshold = None, save = None, verbose = False, 
                    subStack = None, out = sys.stdout, **parameter):
    """Find cell shapes as labeled image
    
    Arguments:
        img (array): image data
        peaks (array): point data of cell centers / seeds
        detectCellShape (dict):
            ============ =================== ===========================================================
            Name         Type                Descritption
            ============ =================== ===========================================================
            *threshold*  (float or None)     threshold to determine mask, pixel below this are background
                                             if None no mask is generated
            *save*       (tuple)             size of the box on which to perform the *method*
            *verbose*    (bool or int)       print / plot information about this step 
            ============ =================== ===========================================================
        verbose (bool): print progress info 
        out (object): object to write progress info to
        
    Returns:
        array: labeled image where each label indicates a cell 
    """    
    
    threshold = getParameter(detectCellShapeParameter, "threshold", threshold);
    save      = getParameter(detectCellShapeParameter, "save", save);    
    verbose   = getParameter(detectCellShapeParameter, "verbose", verbose);  
    
    if verbose:
        writeParameter(out = out, head = 'Cell shape detection:', threshold = threshold, save = save);    
    
    # extended maxima
    timer = Timer();
    
    if threshold is None:
        imgmask = None;
    else:
        imgmask = img > threshold;
        
    imgpeaks = voxelizePixel(peaks, dataSize = img.shape, weights = numpy.arange(1, peaks.shape[0]+1));
    
    imgws = watershed(-img, imgpeaks, mask = imgmask);
    #imgws = watershed_ift(-img.astype('uint16'), imgpeaks);
    #imgws[numpy.logical_not(imgmask)] = 0;
    
    if not save is None:
        writeSubStack(save, imgws.astype('int32'), subStack = subStack);
    
    
    if verbose > 1:
        #plotTiling(img)
        plotOverlayLabel(img * 0.01, imgws, alpha = False);
        #plotOverlayLabel(img, imgmax.astype('int64'), alpha = True)     
    
    if verbose:
        out.write(timer.elapsedTime(head = 'Cell Shape:') + '\n');
    
    return imgws 

def get_segmented_lungs(image):
    #Creation of the markers as shown above:
    marker_internal, marker_external, marker_watershed = generate_markers(image)
    
    #Creation of the Sobel-Gradient
    sobel_filtered_dx = ndimage.sobel(image, 1)
    sobel_filtered_dy = ndimage.sobel(image, 0)
    sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)
    sobel_gradient *= 255.0 / np.max(sobel_gradient)
    
    #Watershed algorithm
    watershed = morphology.watershed(sobel_gradient, marker_watershed)
    
    #Reducing the image created by the Watershed algorithm to its outline
    outline = ndimage.morphological_gradient(watershed, size=(3,3))
    outline = outline.astype(bool)
    
    #Performing Black-Tophat Morphology for reinclusion
    #Creation of the disk-kernel and increasing its size a bit
    blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],
                       [0, 1, 1, 1, 1, 1, 0],
                       [1, 1, 1, 1, 1, 1, 1],
                       [1, 1, 1, 1, 1, 1, 1],
                       [1, 1, 1, 1, 1, 1, 1],
                       [0, 1, 1, 1, 1, 1, 0],
                       [0, 0, 1, 1, 1, 0, 0]]
    #blackhat_struct = ndimage.iterate_structure(blackhat_struct, 8)
    blackhat_struct = ndimage.iterate_structure(blackhat_struct, 14) # <- retains more of the area, 12 works well. Changed to 14, 12 still excluded some parts.
    #Perform the Black-Hat
    outline += ndimage.black_tophat(outline, structure=blackhat_struct)
    
    #Use the internal marker and the Outline that was just created to generate the lungfilter
    lungfilter = np.bitwise_or(marker_internal, outline)
    #Close holes in the lungfilter
    #fill_holes is not used here, since in some slices the heart would be reincluded by accident
    lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)
    
    #Apply the lungfilter (note the filtered areas being assigned threshold_min HU)
    segmented = np.where(lungfilter == 1, image, threshold_min*np.ones(image.shape))
    
    #return segmented, lungfilter, outline, watershed, sobel_gradient, marker_internal, marker_external, marker_watershed
    return segmented 

def get_segmented_lungs(image):
    #Creation of the markers as shown above:
    marker_internal, marker_external, marker_watershed = generate_markers(image)
    
    #Creation of the Sobel-Gradient
    sobel_filtered_dx = ndimage.sobel(image, 1)
    sobel_filtered_dy = ndimage.sobel(image, 0)
    sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)
    sobel_gradient *= 255.0 / np.max(sobel_gradient)
    
    #Watershed algorithm
    watershed = morphology.watershed(sobel_gradient, marker_watershed)
    
    #Reducing the image created by the Watershed algorithm to its outline
    outline = ndimage.morphological_gradient(watershed, size=(3,3))
    outline = outline.astype(bool)
    
    #Performing Black-Tophat Morphology for reinclusion
    #Creation of the disk-kernel and increasing its size a bit
    blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],
                       [0, 1, 1, 1, 1, 1, 0],
                       [1, 1, 1, 1, 1, 1, 1],
                       [1, 1, 1, 1, 1, 1, 1],
                       [1, 1, 1, 1, 1, 1, 1],
                       [0, 1, 1, 1, 1, 1, 0],
                       [0, 0, 1, 1, 1, 0, 0]]
    #blackhat_struct = ndimage.iterate_structure(blackhat_struct, 8)
    blackhat_struct = ndimage.iterate_structure(blackhat_struct, 14) # <- retains more of the area, 12 works well. Changed to 14, 12 still excluded some parts.
    #Perform the Black-Hat
    outline += ndimage.black_tophat(outline, structure=blackhat_struct)
    
    #Use the internal marker and the Outline that was just created to generate the lungfilter
    lungfilter = np.bitwise_or(marker_internal, outline)
    #Close holes in the lungfilter
    #fill_holes is not used here, since in some slices the heart would be reincluded by accident
    lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)
    
    #Apply the lungfilter (note the filtered areas being assigned threshold_min HU)
    segmented = np.where(lungfilter == 1, image, threshold_min*np.ones(image.shape))
    
    #return segmented, lungfilter, outline, watershed, sobel_gradient, marker_internal, marker_external, marker_watershed
    return segmented 

def get_segmented_lungs(image):
    #Creation of the markers as shown above:
    marker_internal, marker_external, marker_watershed = generate_markers(image)
    
    #Creation of the Sobel-Gradient
    sobel_filtered_dx = ndimage.sobel(image, 1)
    sobel_filtered_dy = ndimage.sobel(image, 0)
    sobel_gradient = np.hypot(sobel_filtered_dx, sobel_filtered_dy)
    sobel_gradient *= 255.0 / np.max(sobel_gradient)
    
    #Watershed algorithm
    watershed = morphology.watershed(sobel_gradient, marker_watershed)
    
    #Reducing the image created by the Watershed algorithm to its outline
    outline = ndimage.morphological_gradient(watershed, size=(3,3))
    outline = outline.astype(bool)
    
    #Performing Black-Tophat Morphology for reinclusion
    #Creation of the disk-kernel and increasing its size a bit
    blackhat_struct = [[0, 0, 1, 1, 1, 0, 0],
                       [0, 1, 1, 1, 1, 1, 0],
                       [1, 1, 1, 1, 1, 1, 1],
                       [1, 1, 1, 1, 1, 1, 1],
                       [1, 1, 1, 1, 1, 1, 1],
                       [0, 1, 1, 1, 1, 1, 0],
                       [0, 0, 1, 1, 1, 0, 0]]
    #blackhat_struct = ndimage.iterate_structure(blackhat_struct, 8)
    blackhat_struct = ndimage.iterate_structure(blackhat_struct, 14) # <- retains more of the area, 12 works well. Changed to 14, 12 still excluded some parts.
    #Perform the Black-Hat
    outline += ndimage.black_tophat(outline, structure=blackhat_struct)
    
    #Use the internal marker and the Outline that was just created to generate the lungfilter
    lungfilter = np.bitwise_or(marker_internal, outline)
    #Close holes in the lungfilter
    #fill_holes is not used here, since in some slices the heart would be reincluded by accident
    lungfilter = ndimage.morphology.binary_closing(lungfilter, structure=np.ones((5,5)), iterations=3)
    
    #Apply the lungfilter (note the filtered areas being assigned threshold_min HU)
    segmented = np.where(lungfilter == 1, image, threshold_min*np.ones(image.shape))
    
    #return segmented, lungfilter, outline, watershed, sobel_gradient, marker_internal, marker_external, marker_watershed
    return segmented 

def watershed_segmentation(rgb_img, mask, distance=10):
    """Uses the watershed algorithm to detect boundary of objects. Needs a marker file which specifies area which is
       object (white), background (grey), unknown area (black).

    Inputs:
    rgb_img             = image to perform watershed on needs to be 3D (i.e. np.shape = x,y,z not np.shape = x,y)
    mask                = binary image, single channel, object in white and background black
    distance            = min_distance of local maximum

    Returns:
    analysis_images     = list of output images

    :param rgb_img: numpy.ndarray
    :param mask: numpy.ndarray
    :param distance: int
    :return analysis_images: list
    """
    params.device += 1

    # Store debug mode
    debug = params.debug
    params.debug = None

    dist_transform = cv2.distanceTransformWithLabels(mask, cv2.DIST_L2, maskSize=0)[0]

    localMax = peak_local_max(dist_transform, indices=False, min_distance=distance, labels=mask)

    markers = ndi.label(localMax, structure=np.ones((3, 3)))[0]
    dist_transform1 = -dist_transform
    labels = watershed(dist_transform1, markers, mask=mask)

    img1 = np.copy(rgb_img)

    for x in np.unique(labels):
        rand_color = color_palette(len(np.unique(labels)))
        img1[labels == x] = rand_color[x]

    img2 = apply_mask(img1, mask, 'black')

    joined = np.concatenate((img2, rgb_img), axis=1)

    estimated_object_count = len(np.unique(markers)) - 1

    # Reset debug mode
    params.debug = debug
    if params.debug == 'print':
        print_image(dist_transform, os.path.join(params.debug_outdir, str(params.device) + '_watershed_dist_img.png'))
        print_image(joined, os.path.join(params.debug_outdir, str(params.device) + '_watershed_img.png'))
    elif params.debug == 'plot':
        plot_image(dist_transform, cmap='gray')
        plot_image(joined)

    outputs.add_observation(variable='estimated_object_count', trait='estimated object count',
                            method='plantcv.plantcv.watershed', scale='none', datatype=int,
                            value=estimated_object_count, label='none')

    # Store images
    outputs.images.append([dist_transform, joined])

    return joined 

def _pharynx_orient(worm_img, min_blob_area):#, min_dist_btw_peaks=5):
    #%%
    
    blur = cv2.GaussianBlur(worm_img,(5,5),0) 
    
    #ret3,th3 = cv2.threshold(blur,0,255,cv2.THRESH_TOZERO+cv2.THRESH_OTSU)
    th, worm_mask = cv2.threshold(blur,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
    worm_cnt, cnt_area = binaryMask2Contour(worm_mask, min_blob_area=min_blob_area)
    
    worm_mask = np.zeros_like(worm_mask)
    cv2.drawContours(worm_mask, [worm_cnt.astype(np.int32)], 0, 1, -1)
    
    local_maxi = peak_local_max(blur,
                                indices=True, 
                                labels=worm_mask)
    
    #%%
    markers = np.zeros_like(worm_mask, dtype=np.uint8)
    kernel = np.ones((3,3),np.uint8)
    for x in local_maxi:
        markers[x[0], x[1]] = 1
    markers = cv2.dilate(markers,kernel,iterations = 1)
    markers = ndi.label(markers)[0]
    #strel = ndi.generate_binary_structure(3, 3)
    #markers = binary_dilation(markers, iterations=3)
    
    labels = watershed(-blur, markers, mask=worm_mask)
    props = regionprops(labels)
    
    #sort coordinates by area (the larger area is the head)
    props = sorted(props, key=lambda x: x.area, reverse=True)
    peaks_dict = {labels[x[0], x[1]]:x[::-1] for x in local_maxi}
    peaks_coords = np.array([peaks_dict[x.label] for x in props])
    
    if DEBUG:
        plt.figure()
        plt.subplot(1,3,1)
        plt.imshow(markers, cmap='gray', interpolation='none')
        
        plt.subplot(1,3,2)
        plt.imshow(labels)
        
        plt.subplot(1,3,3)
        plt.imshow(blur, cmap='gray', interpolation='none')
        
        for x,y in peaks_coords:
            plt.plot(x,y , 'or')
            
    if len(props) != 2:
        return np.full((2,2), np.nan) #invalid points return empty
        
    #%%
    return peaks_coords 

def create_mask(self,  labels):
        labels = measure.label(labels, neighbors=8, background=0)
        tmp = dilation(labels > 0, square(9))
        tmp2 = watershed(tmp, labels, mask=tmp, watershed_line=True) > 0
        tmp = tmp ^ tmp2
        tmp = dilation(tmp, square(7))
        msk = (255 * tmp).astype('uint8')

        props = measure.regionprops(labels)
        msk0 = 255 * (labels > 0)
        msk0 = msk0.astype('uint8')

        msk1 = np.zeros_like(labels, dtype='bool')

        max_area = np.max([p.area for p in props])

        for y0 in range(labels.shape[0]):
            for x0 in range(labels.shape[1]):
                if not tmp[y0, x0]:
                    continue
                if labels[y0, x0] == 0:
                    if max_area > 4000:
                        sz = 6
                    else:
                        sz = 3
                else:
                    sz = 3
                    if props[labels[y0, x0] - 1].area < 300:
                        sz = 1
                    elif props[labels[y0, x0] - 1].area < 2000:
                        sz = 2
                uniq = np.unique(labels[max(0, y0 - sz):min(labels.shape[0], y0 + sz + 1),
                                 max(0, x0 - sz):min(labels.shape[1], x0 + sz + 1)])
                if len(uniq[uniq > 0]) > 1:
                    msk1[y0, x0] = True
                    msk0[y0, x0] = 0

        msk1 = 255 * msk1
        msk1 = msk1.astype('uint8')

        msk2 = np.zeros_like(labels, dtype='uint8')
        msk = np.stack((msk0, msk1, msk2))
        msk = np.rollaxis(msk, 0, 3)
        return msk 

def create_mask(labels):
    labels = measure.label(labels, neighbors=8, background=0)
    tmp = dilation(labels > 0, square(9))    
    tmp2 = watershed(tmp, labels, mask=tmp, watershed_line=True) > 0
    tmp = tmp ^ tmp2
    tmp = dilation(tmp, square(7))
    msk = (255 * tmp).astype('uint8')
    
    props = measure.regionprops(labels)
    msk0 = 255 * (labels > 0)
    msk0 = msk0.astype('uint8')
    
    msk1 = np.zeros_like(labels, dtype='bool')
    
    max_area = np.max([p.area for p in props])
    
    for y0 in range(labels.shape[0]):
        for x0 in range(labels.shape[1]):
            if not tmp[y0, x0]:
                continue
            if labels[y0, x0] == 0:
                if max_area > 4000:
                    sz = 6
                else:
                    sz = 3
            else:
                sz = 3
                if props[labels[y0, x0] - 1].area < 300:
                    sz = 1
                elif props[labels[y0, x0] - 1].area < 2000:
                    sz = 2
            uniq = np.unique(labels[max(0, y0-sz):min(labels.shape[0], y0+sz+1), max(0, x0-sz):min(labels.shape[1], x0+sz+1)])
            if len(uniq[uniq > 0]) > 1:
                msk1[y0, x0] = True
                msk0[y0, x0] = 0
    
    msk1 = 255 * msk1
    msk1 = msk1.astype('uint8')
    
    msk2 = np.zeros_like(labels, dtype='uint8')
    msk = np.stack((msk0, msk1, msk2))
    msk = np.rollaxis(msk, 0, 3)
    return msk