def stk_to_rois(stk, threshold, min_size, max_window=8, downscale_factor=2):
    thresholded_stk = stk > threshold
    thresholded_stk = remove_small_objects(thresholded_stk, min_size)
    distance = ndi.distance_transform_edt(thresholded_stk)
    cropped_stk = stk.copy()
    cropped_stk[np.logical_not(thresholded_stk)] = 0
    combined_stk = cropped_stk + distance/distance.max()
    local_max = peak_local_max(combined_stk, indices=False, 
                               footprint=np.ones((max_window, max_window)), 
                               labels=thresholded_stk)
    markers = ndi.label(local_max)[0]
    labels = watershed(-combined_stk, markers, mask=thresholded_stk)
    new_markers = markers.copy()
    for i in set(labels.flatten()):
        if i == 0: continue
        if np.sum(labels==i) < min_size:
            new_markers[markers==i] = 0
    labels = watershed(-combined_stk, new_markers, mask=thresholded_stk)
    labels_set = set(labels.flatten())
    rois = []
    for label in labels_set:
        if label == 0: continue
        if np.sum((labels==label).astype(int)) < min_size: continue
        nroi = np.zeros((stk.shape[0], stk.shape[1]))
        cx,cy = np.where(labels==label)
        cx,cy = int(cx.mean()), int(cy.mean())
        x,y = np.ogrid[0:nroi.shape[0], 0:nroi.shape[1]]
        r = 4
        mask =  (cx-x)**2 + (cy-y)**2 <= r*r
        nroi[mask] = 1
        #nroi[labels==label] = 1
        rois.append(zoom(nroi, downscale_factor, order=0))
    rois = np.array(rois)
    return rois, thresholded_stk, labels

def do_watershed(image, markers,  tfile, shape, bstruct, algorithm, mg_size, use_ww_wl, wl, ww, q):
    mask = np.memmap(tfile, shape=shape, dtype='uint8', mode='r+')
                
    if use_ww_wl:
        if algorithm == 'Watershed':
            tmp_image = ndimage.morphological_gradient(
                           get_LUT_value(image, ww, wl).astype('uint16'),
                           mg_size)
            tmp_mask = watershed(tmp_image, markers.astype('int16'), bstruct)
        else:
            tmp_image = get_LUT_value(image, ww, wl).astype('uint16')
            #tmp_image = ndimage.gaussian_filter(tmp_image, self.config.mg_size)
            #tmp_image = ndimage.morphological_gradient(
                           #get_LUT_value(image, ww, wl).astype('uint16'),
                           #self.config.mg_size)
            tmp_mask = watershed_ift(tmp_image, markers.astype('int16'), bstruct)
    else:
        if algorithm == 'Watershed':
            tmp_image = ndimage.morphological_gradient((image - image.min()).astype('uint16'), mg_size)
            tmp_mask = watershed(tmp_image, markers.astype('int16'), bstruct)
        else:
            tmp_image = (image - image.min()).astype('uint16')
            #tmp_image = ndimage.gaussian_filter(tmp_image, self.config.mg_size)
            #tmp_image = ndimage.morphological_gradient((image - image.min()).astype('uint16'), self.config.mg_size)
            tmp_mask = watershed_ift(tmp_image, markers.astype('int8'), bstruct)
    mask[:] = tmp_mask
    mask.flush()
    q.put(1)

    def analyse(self, **kwargs):
        image_object = kwargs['image']

        if image_object is None:
            raise RuntimeError()

        # Read the image
        image = cv2.imread(self.image_utils.getOutputFilename(image_object.id))

        if image is None:
            print('File not found')
            return

        # Work on green channel
        gray = image[:, :, 1]

        # Apply otsu thresholding
        thresh = filters.threshold_otsu(gray)
        gray[gray < thresh] = 0

        # Apply histogram equalization
        gray = exposure.equalize_adapthist(gray) * 255

        # Create elevation map
        elevation_map = filters.sobel(gray)

        gray = gray.astype(int)

        # Create cell markers
        markers = numpy.zeros_like(gray)
        markers[gray < 100] = 2  # seen as white in plot
        markers[gray > 150] = 1  # seen as black in plot

        # Segment with watershed using elevation map
        segmentation = morphology.watershed(elevation_map, markers)
        segmentation = ndi.binary_fill_holes(segmentation - 1)
        # labeled_image, n = ndi.label(segmentation)

        # Watershed with distance transform
        kernel = numpy.ones((5, 5), numpy.uint8)

        distance = ndi.distance_transform_edt(segmentation)
        distance2 = cv2.erode(distance, kernel)
        distance2 = cv2.dilate(distance2, kernel)
        local_max = peak_local_max(distance2, num_peaks=1, indices=False, labels=segmentation)
        markers2 = ndi.label(local_max)[0]
        labels = morphology.watershed(-distance2, markers2, mask=segmentation)

        # Extract regions (caching signifies more memory use)
        regions = regionprops(labels, cache=True)

        # Filter out big wrong regions
        regions = [region for region in regions if region.area < 2000]

        # Set result
        result = str(len(regions))

        return result

def LargestWatershedRegion(shapes,dims,skipBias): 
    L=len(shapes)-skipBias     
    shapes=shapes.reshape((-1,) + dims[1:]) 
    D=len(dims)
    num_peaks=4
#    structure=np.ones(tuple(3*np.ones((np.ndim(shapes)-1,1))))
    for ll in range(L): 
        temp=shapes[ll]
        local_maxi = peak_local_max(gaussian_filter(temp,[1]*(D-1)), exclude_border=False, indices=False, num_peaks=num_peaks)
        markers,junk = label(local_maxi)
        nonzero_mask=temp>0
        if np.sum(nonzero_mask)>(3**3)*num_peaks:
            labels = watershed(-temp, markers, mask=nonzero_mask) #watershed regions
            ind = 1
            temp2 = np.copy(temp)
            temp2[labels!=1]=0
            total_intensity = sum(temp2.reshape(-1,))
            for kk in range(2,labels.max()+1):
                temp2 = np.copy(temp)
                temp2[labels!=kk]=0
                total_intensity2 = sum(temp2.reshape(-1,))
                if total_intensity2>total_intensity:
                    ind = kk
                    total_intensity=total_intensity2
            temp[labels!=ind]=0
            shapes[ll]=temp
    shapes=shapes.reshape((len(shapes),-1))
    return shapes

    def segment(self, src):
        '''
            pychron: preprocessing cv.Mat
        '''
#        image = pychron.ndarray[:]
#         image = asarray(pychron)
        image = src[:]
        if self.use_adaptive_threshold:
#            block_size = 25
            markers = threshold_adaptive(image, self.block_size)

            n = markers[:].astype('uint8')
            n[markers == True] = 255
            n[markers == False] = 1
            markers = n

        else:
            markers = zeros_like(image)
            markers[image < self.threshold_low] = 1
            markers[image > self.threshold_high] = 255

        elmap = sobel(image, mask=image)
        wsrc = watershed(elmap, markers, mask=image)

#         wsrc = wsrc.astype('uint8')
        return invert(wsrc)

    def segment(self, data, peaks):
        """Perform a watershed segmentation based on local maxima."""
        markers = np.zeros_like(data)
        markers[tuple(peaks.T)] = np.arange(len(peaks)) + 1
        seg = morphology.watershed(-data, markers, mask=data > 0)

        return seg, markers

def detectOpticDisc(image):
    kernel = octagon(10, 10)
    thresh = threshold_otsu(image[:,:,1])
    binary = image > thresh
    print binary.dtype
    luminance = convertToHLS(image)[:,:,2]
    t = threshold_otsu(luminance)
    t = erosion(luminance, kernel)
    
    
    labels = segmentation.slic(image[:,:,1], n_segments = 3)
    out = color.label2rgb(labels, image[:,:,1], kind='avg')
    skio.imshow(out)
    
    x, y = computeCentroid(t)
    print x, y
    rows, cols, _ = image.shape
    p1 = closing(image[:,:,1],kernel)
    p2 = opening(p1, kernel)
    p3 = reconstruction(p2, p1, 'dilation')
    p3 = p3.astype(np.uint8)
    #g = dilation(p3, kernel)-erosion(p3, kernel)
    #g = rank.gradient(p3, disk(5))
    g = cv2.morphologyEx(p3, cv2.MORPH_GRADIENT, kernel)
    #markers = rank.gradient(p3, disk(5)) < 10
    markers = drawCircle(rows, cols, x, y, 85)
    #markers = ndimage.label(markers)[0]
    #skio.imshow(markers)
    g = g.astype(np.uint8)
    #g = cv2.cvtColor(g, cv2.COLOR_GRAY2RGB)
    w = watershed(g, markers)
    print np.max(w), np.min(w)
    w = w.astype(np.uint8)
    #skio.imshow(w)
    return w

def segment(image, thresh):
    #preprocess image
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    
    #perform euclidean distance transform
    distances = ndimage.distance_transform_edt(thresh)
    localMax = peak_local_max(distances, indices = False, min_distance = 3, labels = thresh)

    #perform connected component analysis on local peaks
    markers = ndimage.label(localMax, structure = np.ones((3, 3)))[0]
    labels = watershed(-distances, markers, mask = thresh)

    #loop over labels returned from watershed to mark them
    for label in np.unique(labels):
        if label == 0:
            continue

        mask = np.zeros(gray.shape, dtype="uint8")
        mask[labels == label] = 255

        #find contours in mask and choose biggest contour by area
        contours = cv2.findContours(mask.copy(), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)[-2]
        contour = max(contours, key = cv2.contourArea)

        #draw circle around max size contour
        ((x, y), r) = cv2.minEnclosingCircle(contour)
        cv2.circle(image, (int(x), int(y)), int(r), (0, 255, 0), 2)
    
    #show final image
    cv2.imshow("Output", image)
    return len(np.unique(labels) - 1)

def watershed(image):
    hsv_image = color.rgb2hsv(image)

    low_res_image = rescale(hsv_image[:, :, 0], SCALE)
    local_mean = mean(low_res_image, disk(50))
    local_minimum_flat = np.argmin(local_mean)
    local_minimum = np.multiply(np.unravel_index(local_minimum_flat, low_res_image.shape), round(1 / SCALE))

    certain_bone_pixels = np.full_like(hsv_image[:, :, 0], False, bool)
    certain_bone_pixels[
    local_minimum[0] - INITIAL_WINDOW_SIZE/2:local_minimum[0]+INITIAL_WINDOW_SIZE/2,
    local_minimum[1] - INITIAL_WINDOW_SIZE/2:local_minimum[1]+INITIAL_WINDOW_SIZE/2
    ] = True

    certain_non_bone_pixels = np.full_like(hsv_image[:, :, 0], False, bool)
    certain_non_bone_pixels[0:BORDER_SIZE, :] = True
    certain_non_bone_pixels[-BORDER_SIZE:-1, :] = True
    certain_non_bone_pixels[:, 0:BORDER_SIZE] = True
    certain_non_bone_pixels[:, -BORDER_SIZE:-1] = True

    smoothed_hsv = median(hsv_image[:, :, 0], disk(50))
    threshold = MU * np.median(smoothed_hsv[certain_bone_pixels])

    possible_bones = np.zeros_like(hsv_image[:, :, 0])
    possible_bones[smoothed_hsv < threshold] = 1

    markers = np.zeros_like(possible_bones)
    markers[certain_bone_pixels] = 1
    markers[certain_non_bone_pixels] = 2

    labels = morphology.watershed(-possible_bones, markers)

    return labels

    def expand_watershed(self, pubsub_evt):
        markers = self.matrix
        image = self.viewer.slice_.matrix
        self.viewer.slice_.do_threshold_to_all_slices()
        mask = self.viewer.slice_.current_mask.matrix[1:, 1:, 1:]
        ww = self.viewer.slice_.window_width
        wl = self.viewer.slice_.window_level
        if BRUSH_BACKGROUND in markers and BRUSH_FOREGROUND in markers:
            tmp_image = ndimage.morphological_gradient(get_LUT_value(image, ww, wl).astype('uint16'), self.mg_size)
            tmp_mask = watershed(tmp_image, markers)

            if self.viewer.overwrite_mask:
                mask[:] = 0
                mask[tmp_mask == 1] = 253
            else:
                mask[(tmp_mask==2) & ((mask == 0) | (mask == 2) | (mask == 253))] = 2
                mask[(tmp_mask==1) & ((mask == 0) | (mask == 2) | (mask == 253))] = 253

            #mask[:] = tmp_mask
            self.viewer.slice_.current_mask.matrix[0] = 1
            self.viewer.slice_.current_mask.matrix[:, 0, :] = 1
            self.viewer.slice_.current_mask.matrix[:, :, 0] = 1

            self.viewer.slice_.discard_all_buffers()
            self.viewer.slice_.current_mask.clear_history()
            Publisher.sendMessage('Reload actual slice')

    def segment(self, src):
        image = src.ndarray[:]
        if self.use_adaptive_threshold:
            block_size = 25
            markers = threshold_adaptive(image, block_size) * 255
            markers = invert(markers)

        else:
            markers = zeros_like(image)
            markers[image < self.threshold_low] = 1
            markers[image > self.threshold_high] = 255

        elmap = sobel(image, mask=image)
        wsrc = watershed(elmap, markers, mask=image)

#        elmap = ndimage.distance_transform_edt(image)
#        local_maxi = is_local_maximum(elmap, image,
#                                      ones((3, 3))
#                                      )
#        markers = ndimage.label(local_maxi)[0]
#        wsrc = watershed(-elmap, markers, mask=image)
#        fwsrc = ndimage.binary_fill_holes(out)
#        return wsrc
        if self.use_inverted_image:
            out = invert(wsrc)
        else:
            out = wsrc

#        time.sleep(1)
#        do_later(lambda:self.show_image(image, -elmap, out))
        return out

def testSkimage():
    img = Image.open('../img/1.png')
    img = np.array(img)
    imggray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    # (thresh, imgbw) = cv2.threshold(imggray, 128, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)

    # canny detector
    # from skimage.feature import canny
    # edges = canny(imggray/ 255.)
    from scipy import ndimage as ndi
    # fill_imgbw = ndi.binary_fill_holes(edges)
    # label_objects, nb_labels = ndi.label(fill_imgbw)
    # sizes = np.bincount(label_objects.ravel())
    # mask_sizes = sizes > 20
    # mask_sizes[0] = 0
    # cleaned_imgbw = mask_sizes[label_objects]

    markers = np.zeros_like(imggray)
    markers[imggray < 120] = 1
    markers[imggray > 150] = 2

    from skimage.filters import sobel
    elevation_map = sobel(imggray)
    from skimage.morphology import watershed
    segmentation = watershed(elevation_map, markers)

    # from skimage.color import label2rgb
    # segmentation = ndi.binary_fill_holes(segmentation - 10)
    # labeled_coins, _ = ndi.label(segmentation)
    # image_label_overlay = label2rgb(labeled_coins, image=imggray)
    plt.imshow(segmentation, cmap='gray')
    plt.show()
    return

def black_background(image, kernel):

    shifted = cv2.pyrMeanShiftFiltering(image, 10, 39)
    gray = cv2.cvtColor(shifted, cv2.COLOR_BGR2GRAY)
    thresh = cv2.threshold(gray, 0, 255,
                           cv2.THRESH_BINARY | cv2.THRESH_OTSU)[1]
    D = ndimage.distance_transform_edt(thresh)
    localMax = peak_local_max(D, indices=False, min_distance=10,
                              labels=thresh)

    # perform a connected component analysis on the local peaks,
    # using 8-connectivity, then appy the Watershed algorithm
    markers = ndimage.label(localMax, structure=np.ones((3, 3)))[0]
    labels = watershed(-D, markers, mask=thresh)
    # create a mask
    mask2 = np.zeros(gray.shape, dtype="uint8")
    #  loop over the unique labels returned by the Watershed  algorithm for
    for label in np.unique(labels):
        # if the label is zero, we are examining the 'background' so simply ignore it
        if label == 0:
            continue
        # otherwise, allocate memory for the label region and draw
        # it on the mask
        mask2[labels == label] = 255
    return mask2

def segmentation(file_name):
    data_x, data_y, data_z = get_data(file_name)
    shape_x = len(np.unique(data_x))
    shape_y = len(np.unique(data_y))
    X = data_x.reshape(shape_x, shape_y)
    Y = data_y.reshape(shape_x, shape_y)
    Z = data_z.reshape(shape_x, shape_y)

    markers = np.zeros_like(Z)
    markers[Z < 0.15] = 1
    markers[Z > 0.3] = 2
    elevation_map = roberts(Z)
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3), sharex=True, sharey=True)
    # ax.imshow(Z)
    # ax.imshow(elevation_map, cmap=plt.cm.jet, interpolation='nearest')
    segmentation = watershed(elevation_map, markers)
    ax2.imshow(segmentation, interpolation='nearest')
    # ax.axis('off')
    # ax.set_title('segmentation')
    segmentation = ndi.binary_fill_holes(segmentation - 1)
    labeled_coins, _ = ndi.label(segmentation)
    ax1.imshow(Z, cmap=plt.cm.gray, interpolation='nearest')
    ax1.contour(segmentation, [0.5], linewidths=1.2, colors='y')
    ax1.axis('off')
    ax1.set_adjustable('box-forced')

    plt.show()

 def run(self, workspace):
     labeled_nuclei = workspace.object_set.get_objects(self.primary_objects.value).get_segmented()
     cell_image = workspace.image_set.get_image(self.image_name.value).pixel_data[:,:]
     image_collection = []
     cell_treshold = otsu(cell_image, min_threshold=0, max_threshold=1)
     
     cell_binary = (cell_image >= cell_treshold)
     cell_distance = scipym.distance_transform_edt(cell_binary).astype(np.uint16)
     cell_labeled = skm.watershed(-cell_distance, labeled_nuclei, mask=cell_binary)
     
      
     #
     #fil hall and filter on syze the object in cell_labeled
     #
     cell_labeled = self.filter_on_border(cell_labeled)
     cell_labeled = fill_labeled_holes(cell_labeled)
 
     objects = cellprofiler.objects.Objects()
     objects.segmented = cell_labeled
     objects.parent_image = cell_image
     
     workspace.object_set.add_objects(objects, self.object_name.value)        
     image_collection.append((cell_image, "Original"))
     image_collection.append((cell_labeled, "Labelized image"))
     workspace.display_data.image_collection = image_collection

    def __call__(self, image, window_size=10, threshold=0, fill_holes=True,
                 outline_smoothing=2, remove_borderobjects=True, size_min=1,
                 *args, **kw):

        thresh = threshold_adaptive(image, block_size=window_size,
                                    offset=-1*threshold)

        if outline_smoothing >= 1:
            thresh = outlineSmoothing(thresh, outline_smoothing)

        thresh = remove_small_objects(thresh, size_min)

        seeds = ndi.label(clear_border(~thresh))[0]
        thresh = ndi.binary_fill_holes(thresh)
        smask = seeds.astype(bool)

        # object don't touch border after outline smoothing
        if remove_borderobjects:
            thresh = clear_border(thresh)

        img = np.zeros(thresh.shape)
        img[~smask] = 1
        edt = ndi.morphology.distance_transform_edt(img)
        edt -= ndi.morphology.distance_transform_edt(seeds)

        labels = watershed(edt, seeds)
        labels[smask] = 0
        labels[~thresh] = 0

        return labels

def Image_ws_tranche(image):
    
    laser = Detect_laser(image)
    laser_tranche = tranche_image(laser,60)
    
    image_g = skimage.color.rgb2gray(image)
    image_g = image_g * laser_tranche
    
    image_med = rank2.median((image_g*255).astype('uint8'),disk(8))
    
    image_clahe = exposure.equalize_adapthist(image_med, clip_limit=0.03)
    image_clahe_stretch = exposure.rescale_intensity(image_clahe, out_range=(0, 256))

    image_grad = rank2.gradient(image_clahe_stretch,disk(3))
    
    image_grad_mark = image_grad<20
    image_grad_forws = rank2.gradient(image_clahe_stretch,disk(1))
    
    image_grad_mark_closed = closing(image_grad_mark,disk(1))
    
    Labelised = (skimage.measure.label(image_grad_mark_closed,8,0))+1
    Watersheded  = watershed(image_grad_forws,Labelised)
    
    cooc = coocurence_liste(Watersheded,laser,3)
    
    x,y = compte_occurences(cooc)
    
    return x,y

def seeded_watershed(boundary, seed_threshold = 0, seed_size = 5, mask=None):
    """Extract seeds from boundary prediction and runs seeded watershed.

    Args:
        boundary (3D numpy array) = boundary predictions
        seed_threshold (int) = Add seeds where boundary prob is <= threshold
        seed_size (int) = seeds must be >= seed size
        mask (3D numpy array) = true to watershed, false to ignore
    Returns:
        3d watershed
    """
    
    from skimage import morphology as skmorph
    from numpy import bincount 

    # get seeds
    from scipy.ndimage import label as label2
    seeds = label2(boundary <= seed_threshold, output=numpy.uint32)[0]

    # remove small seeds
    if seed_size > 0:
        component_sizes = bincount(seeds.ravel())
        small_components = component_sizes < seed_size 
        small_locations = small_components[seeds]
        seeds[small_locations] = 0

    # mask out background (don't have to 0 out seeds since)
    supervoxels = skmorph.watershed(boundary, seeds,
            None, None, mask)

    return supervoxels

    def segment(self, image):
        """
        """
        # image = src[:]
        if self.use_adaptive_threshold:
            bs = self.blocksize
            if not bs % 2:
                bs += 1

            markers = threshold_adaptive(image, bs)

            # n = markers[:].astype('uint8')
            n = markers.astype('uint8')
            # n[markers] = 255
            # n[invert(markers)] = 1
            # markers = n
            return n
        else:
            markers = zeros_like(image)
            # print('image',image.max(), image.min())
            # print('le', image<self.threshold_low)
            # print('ge', image>self.threshold_high)
            markers[image <= self.threshold_low] = 1
            markers[image > self.threshold_high] = 255

        #elmap = sobel(image, mask=image)
        elmap = canny(image, sigma=1)
        wsrc = watershed(elmap, markers, mask=image)

        return invert(wsrc)

def isolate_single_seed(stack_path, output_path):
    """Load stack then isolate seed."""

    sigma = 10
    iterations = 1

    raw_stack = Image3D.from_path(stack_path)

    print raw_stack.shape

    smoothed = gaussian_filter(raw_stack, sigma).view(Image3D)

    edges = sobel_magnitude_nd(smoothed).view(Image3D)

    #edges.save('edges')

    labels = np.zeros(raw_stack.shape)
    cx, cy, cz = map(lambda x: x/2, raw_stack.shape)
    labels[cx, cy, cz] = 1

    threshold = threshold_otsu(smoothed)
    thresholded = smoothed > threshold

    #thresholded.view(Image3D).save('thresh')

    segmentation = watershed(edges, markers=labels, mask=thresholded)

    #segmentation.view(Image3D).save('seg')

    dilated = binary_dilation(segmentation, iterations=iterations)

    isolate = np.multiply(raw_stack, dilated)

    isolate.view(Image3D).save(output_path)