Python skimage.filters.threshold_otsu() Examples

def __call__(self, image, test):
        image_data = numpy.asarray(image, dtype=numpy.float32)[:, :, :3]
        rgb_image_data = image_data.transpose(2, 0, 1)
        lab_image_data = rgb2lab(image_data / 255).transpose(2, 0, 1).astype(numpy.float32)
        luminous_image_data = lab_image_data[0].astype(numpy.uint8)

        try:
            th = threshold_otsu(luminous_image_data)
        except:
            import traceback
            print(traceback.format_exc())
            th = 0

        linedrawing = (luminous_image_data > th).astype(numpy.float32)
        linedrawing = numpy.expand_dims(linedrawing, axis=0)

        return lab_image_data, linedrawing, rgb_image_data 

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 project_object_edge(img, dimension):
    """ scale the image, binarise with Othu and project to one dimension

    :param ndarray img:
    :param int dimension: select dimension for projection
    :return list(float):

    >>> img = np.zeros((20, 10, 3))
    >>> img[2:6, 1:7, :] = 1
    >>> img[10:17, 4:6, :] = 1
    >>> project_object_edge(img, 0).tolist()  # doctest: +NORMALIZE_WHITESPACE
    [0.0, 0.0, 0.7, 0.7, 0.7, 0.7, 0.0, 0.0, 0.0, 0.0,
     0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.0, 0.0, 0.0]
    """
    assert dimension in (0, 1), 'not supported dimension %i' % dimension
    assert img.ndim == 3, 'unsupported image shape %r' % img.shape
    img_gray = np.mean(img, axis=-1)
    img_gray = GaussianBlur(img_gray, (5, 5), 0)
    p_low, p_high = np.percentile(img_gray, (1, 95))
    img_gray = rescale_intensity(img_gray, in_range=(p_low, p_high))
    img_bin = img_gray > threshold_otsu(img_gray)
    img_edge = np.mean(img_bin, axis=1 - dimension)
    return img_edge 

def predict(obj, ctx):
    """
    A prediction function.

    Args:
        obj (bytearray): a list of object (or a single object) to predict on
        ctx (dict): a model training context

    Returns:
        a list with predictions

    """
    binary_image = Image.open(BytesIO(obj)).convert('L')
    img = resize(np.array(binary_image), (28, 28), preserve_range=True)
    threshold = threshold_otsu(img)
    img_colors_inverted = list(map(lambda el: 0 if el > threshold else int(255 - el), img.flatten()))
    img_transformed = ctx['scaler'].transform([img_colors_inverted])
    df = pd.DataFrame(img_transformed)  # convert to DataFrame (pandas)
    return ctx['model'].predict(df).tolist() 

def remove_background(img: np.ndarray) -> np.ndarray:
        """ Remove noise using OTSU's method.

        Parameters
        ----------
        img : np.ndarray
            The image to be processed

        Returns
        -------
        np.ndarray
            The image with background removed
        """

        img = img.astype(np.uint8)
        # Binarize the image using OTSU's algorithm. This is used to find the center
        # of mass of the image, and find the threshold to remove background noise
        threshold = filters.threshold_otsu(img)

        # Remove noise - anything higher than the threshold. Note that the image is still grayscale
        img[img > threshold] = 255

        return img 

def run(args):
    logging.basicConfig(level=logging.INFO)

    slide = openslide.OpenSlide(args.wsi_path)

    # note the shape of img_RGB is the transpose of slide.level_dimensions
    img_RGB = np.transpose(np.array(slide.read_region((0, 0),
                           args.level,
                           slide.level_dimensions[args.level]).convert('RGB')),
                           axes=[1, 0, 2])

    img_HSV = rgb2hsv(img_RGB)

    background_R = img_RGB[:, :, 0] > threshold_otsu(img_RGB[:, :, 0])
    background_G = img_RGB[:, :, 1] > threshold_otsu(img_RGB[:, :, 1])
    background_B = img_RGB[:, :, 2] > threshold_otsu(img_RGB[:, :, 2])
    tissue_RGB = np.logical_not(background_R & background_G & background_B)
    tissue_S = img_HSV[:, :, 1] > threshold_otsu(img_HSV[:, :, 1])
    min_R = img_RGB[:, :, 0] > args.RGB_min
    min_G = img_RGB[:, :, 1] > args.RGB_min
    min_B = img_RGB[:, :, 2] > args.RGB_min

    tissue_mask = tissue_S & tissue_RGB & min_R & min_G & min_B

    np.save(args.npy_path, tissue_mask) 

def threshold_segmentation(img):
    # calculate the overview level size and retrieve the image
    img_hsv = img.convert('HSV')
    img_hsv_np = np.array(img_hsv)

    # dilate image and then threshold the image
    schannel = img_hsv_np[:, :, 1]
    mask = np.zeros(schannel.shape)

    schannel = dilation(schannel, star(3))
    schannel = ndimage.gaussian_filter(schannel, sigma=(5, 5), order=0)
    threshold_global = threshold_otsu(schannel)

    mask[schannel > threshold_global] = FOREGROUND
    mask[schannel <= threshold_global] = BACKGROUND

    return mask 

def MaxBodyBox(input):
    Otsu=filters.threshold_otsu(input[input.shape[0]//2])
    Seg=np.zeros(input.shape)
    Seg[input>=Otsu]=255
    Seg=Seg.astype(np.int)
    ConnectMap=label(Seg, connectivity= 2)
    Props = regionprops(ConnectMap)
    Area=np.zeros([len(Props)])
    Area=[]
    Bbox=[]
    for j in range(len(Props)):
        Area.append(Props[j]['area'])
        Bbox.append(Props[j]['bbox'])
    Area=np.array(Area)
    Bbox=np.array(Bbox)
    argsort=np.argsort(Area)
    Area=Area[argsort]
    Bbox=Bbox[argsort]
    Area=Area[::-1]
    Bbox=Bbox[::-1,:]
    MaximumBbox=Bbox[0]
    return Otsu,MaximumBbox 

def otsu(*args, **kwargs):
        raise NotImplementedError("Install skimage!")


#####################################################################
## utils for inspecting model layers, replacing layers, ...
##################################################################### 

def ssmdt(dt, ssmiter):
    dt = ssm(dt, anisotropic=True, iterations=ssmiter)
    dt[dt < filters.threshold_otsu(dt)] = 0
    dt = skfmm.distance(dt, dx=5e-2)
    dt = skfmm.distance(np.logical_not(dt), dx=5e-3)
    dt[dt > 0.04] = 0.04
    dt = dt.max() - dt
    dt[dt <= 0.038] = 0
    return dt 

def resize(img, dest_width, dest_height):
    """
    Resize an image with given destination dimensions.
    @Author: Appcell
    @param img: the image to be resized
    @dest_width: width of image after resizing
    @dest_height: height of image after resizing
    @return: a numpy.ndarray object of resized image.
    """
    return cv2.resize(img, (dest_width, dest_height))


# def contrast_adjust_log(img, gain):
#     """
#     Perform logarithm correction to increase contrast
#     @Author: Rigel
#     @param img: the image to be corrected
#     @gain: constant multiplier
#     @return: a numpy.ndarray object of corrected image.
#     """
#     return adjust_log(img, gain)


# def binary_otsu(img):
#     """
#     Perform binarization with otsu algorithm
#     @Author: Rigel
#     @param img: the image to be corrected
#     @gain: constant multiplier
#     @return: a numpy.ndarray boolean object of binary image.
#     """
#     global_thresh = threshold_otsu(img)
#     return img > global_thresh 

def get_regions(slice_or_arr, fill_holes=False, clear_borders=True, threshold='otsu'):
    """Get the skimage regions of a black & white image."""
    if threshold == 'otsu':
        thresmeth = filters.threshold_otsu
    elif threshold == 'mean':
        thresmeth = np.mean
    if isinstance(slice_or_arr, Slice):
        edges = filters.scharr(slice_or_arr.image.array.astype(np.float))
        center = slice_or_arr.image.center
    elif isinstance(slice_or_arr, np.ndarray):
        edges = filters.scharr(slice_or_arr.astype(np.float))
        center = (int(edges.shape[1]/2), int(edges.shape[0]/2))
    edges = filters.gaussian(edges, sigma=1)
    if isinstance(slice_or_arr, Slice):
        box_size = 100/slice_or_arr.mm_per_pixel
        thres_img = edges[int(center.y-box_size):int(center.y+box_size),
                          int(center.x-box_size):int(center.x+box_size)]
        thres = thresmeth(thres_img)
    else:
        thres = thresmeth(edges)
    bw = edges > thres
    if clear_borders:
        segmentation.clear_border(bw, buffer_size=int(max(bw.shape)/50), in_place=True)
    if fill_holes:
        bw = ndimage.binary_fill_holes(bw)
    labeled_arr, num_roi = measure.label(bw, return_num=True)
    regionprops = measure.regionprops(labeled_arr, edges)
    return labeled_arr, regionprops, num_roi 

def Normalization(Image):
    Spacing=Image.GetSpacing()
    Origin = Image.GetOrigin()
    Direction = Image.GetDirection()
    Array=sitk.GetArrayFromImage(Image)
    Array_new=Array.copy()
    Array_new+=np.min(Array_new)
    Array_new=Array_new[Array_new.shape[0]/2-5:Array_new.shape[0]/2+5]
    Mask=Array_new.copy()
    for i in range(Array_new.shape[0]):
        otsu=filters.threshold_otsu(Array_new[i])
        Mask[i][Array_new[i]<0.5*otsu]=0
        Mask[i][Array_new[i]>=0.5*otsu]=1
    MaskSave=sitk.GetImageFromArray(Mask)
    MaskSave=sitk.BinaryDilate(MaskSave,10)
    MaskSave=sitk.BinaryErode(MaskSave,10)
    Mask=sitk.GetArrayFromImage(MaskSave)

    Avg=np.average(Array[Array_new.shape[0]/2-5:Array_new.shape[0]/2+5],weights=Mask)
    Std=np.sqrt(np.average(abs(Array[Array_new.shape[0]/2-5:Array_new.shape[0]/2+5] - Avg)**2,weights=Mask))
    Array=(Array.astype(np.float32)-Avg)/Std
    Array[Array>UpperBound]=UpperBound
    Array[Array<LowerBound]=LowerBound 
    Array=((Array.astype(np.float64)-np.min(Array))/(np.max(Array)-np.min(Array))*255).astype(np.uint8)
    Image=sitk.GetImageFromArray(Array)
    Image.SetDirection(Direction)
    Image.SetOrigin(Origin)
    Image.SetSpacing(Spacing)
    return Image,MaskSave 

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 evaluate(self,img):
    img = cv2.medianBlur(img,5)
    if self.auto_thresh:
      self.thresh = threshold_otsu(img)
    if self.white_spots:
      bw = (img > self.thresh).astype(np.uint8)
    else:
      bw = (img <= self.thresh).astype(np.uint8)
    #cv2.imshow(self.name,bw*255)
    #cv2.waitKey(5)
    if not .1* img.size < np.count_nonzero(bw) < .8*img.size:
      print("reevaluating threshold!!")
      print("Ratio:",np.count_nonzero(bw)/img.size)
      print("old:",self.thresh)
      self.thresh = threshold_otsu(img)
      print("new:",self.thresh)
    m = cv2.moments(bw)
    r = {}
    try:
      r['x'] = m['m10']/m['m00']
      r['y'] = m['m01']/m['m00']
    except ZeroDivisionError:
      return -1
    x,y,w,h = cv2.boundingRect(bw)
    #if (h,w) == img.shape:
    #  return -1
    r['bbox'] = y,x,y+h,x+w
    return r 

def _otsu(data: torch.Tensor, **kwargs):
    """
    Use Otsu's method, which assumes a GMM with 2 components
    but uses some heuristic to maximize the variance differences.
    """
    h = 2 ** int(1 + math.log2(data.shape[0]) / 2)
    fake_img = data.view(h, -1).cpu().numpy()
    return otsu(fake_img, h) 

def autothreshold(data, threshold_type = 'otsu', z = 2.3264, set_max_ceiling = True):
	if threshold_type.endswith('_p'):
		data = data[data>0]
	else:
		data = data[data!=0]
	if (threshold_type == 'otsu') or (threshold_type == 'otsu_p'):
		lthres = filters.threshold_otsu(data)
		uthres = data[data>lthres].mean() + (z*data[data>lthres].std())
		# Otsu N (1979) A threshold selection method from gray-level histograms. IEEE Trans. Sys., Man., Cyber. 9: 62-66.
	elif (threshold_type == 'li')  or (threshold_type == 'li_p'):
		lthres = filters.threshold_li(data)
		uthres = data[data>lthres].mean() + (z*data[data>lthres].std())
		# Li C.H. and Lee C.K. (1993) Minimum Cross Entropy Thresholding Pattern Recognition, 26(4): 617-625
	elif (threshold_type == 'yen') or (threshold_type == 'yen_p'):
		lthres = filters.threshold_yen(data)
		uthres = data[data>lthres].mean() + (z*data[data>lthres].std())
		# Yen J.C., Chang F.J., and Chang S. (1995) A New Criterion for Automatic Multilevel Thresholding IEEE Trans. on Image Processing, 4(3): 370-378.
	elif threshold_type == 'zscore_p':
		lthres = data.mean() - (z*data.std())
		uthres = data.mean() + (z*data.std())
		if lthres < 0:
			lthres = 0.001
	else:
		lthres = data.mean() - (z*data.std())
		uthres = data.mean() + (z*data.std())
	if set_max_ceiling: # for the rare case when uthres is larger than the max value
		if uthres > data.max():
			uthres = data.max()
	return lthres, uthres 

def threshold_mask(raw_output, threshold):
    if threshold == 0:
        try:
            threshold = filters.threshold_otsu(raw_output)
        except:
            threshold = 0.5

    print('\tThreshold: {}'.format(threshold))

    raw_output[raw_output > threshold] = 1
    raw_output[raw_output < 1] = 0

    all_labels = measure.label(raw_output)
    props = measure.regionprops(all_labels)
    props.sort(key=lambda x: x.area, reverse=True)
    thresholded_mask = np.zeros(raw_output.shape)

    if len(props) >= 2:
        if props[0].area / props[1].area > 5:  # if the largest is way larger than the second largest
            thresholded_mask[all_labels == props[0].label] = 1  # only turn on the largest component
        else:
            thresholded_mask[all_labels == props[0].label] = 1  # turn on two largest components
            thresholded_mask[all_labels == props[1].label] = 1
    elif len(props):
        thresholded_mask[all_labels == props[0].label] = 1

    thresholded_mask = scipy.ndimage.morphology.binary_fill_holes(thresholded_mask).astype(np.uint8)

    return thresholded_mask 

def _thresh_bw(pix_valid):
    # calculate otsu_threshold as lower limit. Otsu understimates the threshold.
    try:
        otsu_thresh = skf.threshold_otsu(pix_valid)
    except:
        return np.nan

    # calculate the histogram
    pix_hist = np.bincount(pix_valid)

    # the higher limit is the most frequent value in the distribution
    # (background)
    largest_peak = np.argmax(pix_hist)
    if otsu_thresh < largest_peak and otsu_thresh + 2 < len(pix_hist) - 1:
        # smooth the histogram to find a better threshold
        pix_hist = np.convolve(pix_hist, np.ones(3), 'same')
        cumhist = np.cumsum(pix_hist)

        xx = np.arange(otsu_thresh, cumhist.size)
        try:
            # the threshold is calculated as the first pixel level above the otsu threshold 
            # at which there would be larger increase in the object area.
            hist_ratio = pix_hist[xx] / cumhist[xx]
            thresh = np.where(
                (hist_ratio[3:] - hist_ratio[:-3]) > 0)[0][0] + otsu_thresh
        except IndexError:
            thresh = np.argmin(
                pix_hist[
                    otsu_thresh:largest_peak]) + otsu_thresh
    else:
        # if otsu is larger than the maximum peak keep otsu threshold
        thresh = otsu_thresh
    return thresh 

def __otsu_bin__(img,invert):
    # with scikit-image
    thresh = threshold_otsu(img)
    # threshed_image = img > thresh

    threshed_image = np.zeros(img.shape,np.uint8)
    threshed_image[np.where(img > thresh)] = 255

    # with open cv
    # _,threshed_image = cv2.threshold(img,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)

    return threshed_image,thresh 

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 otsu(data, nbins):
    from skimage.filters import threshold_otsu
    thresh = threshold_otsu(data, nbins)
    return data > thresh 

def remove_noise(img):
    th = filters.threshold_otsu(img)
    bin_img = img > th
    regions, num = ndimage.label(bin_img)
    areas = []
    for i in range(num):
        areas.append(np.sum(regions == i+1))
    img[regions != np.argmax(areas)+1] = 0
    return img 

def compCentroid_detect1(fcn, savefolder):
	data_dict = sio.loadmat(fcn)
	f = matlab_style_gauss2D((10,10),0.25)
	A = cv2.filter2D(data_dict['A'], -1, f)

	level = threshold_otsu(A) #otsu threshold of image
	bw = A > level #binary image
	L,num = label(bw,8,return_num=True) #label  the segmented blobs
	#import pdb;pdb.set_trace()
	plot_x = np.zeros((num, 1)) # location of centroid
	plot_y = np.zeros((num, 1))

	sum_x = np.zeros((num, 1))
	sum_y = np.zeros((num, 1))
	area = np.zeros((num, 1))
	score = np.zeros((num, 1))

	height,width = bw.shape[0], bw.shape[1]
	for i in range(height):
		for j in range(width):
			if L[i,j] != 0:
				N = L[i,j]
				sum_x[N-1] = sum_x[N-1]+i*A[i,j]
				sum_y[N-1] = sum_y[N-1]+j*A[i,j]
				area[N-1] = area[N-1] + 1
				score[N-1] = score[N-1] + A[i,j]

	plot_x = np.around(sum_x*1.0/score)
	plot_y = np.around(sum_y*1.0/score)
	score = score*1.0/area
	centroid = np.zeros((num,2))
	for row in range(num):
		centroid[row,0] = plot_x[row,0]
		centroid[row,1] = plot_y[row,0]
	#centroid = np.mat(centroid)
	savefile = savefolder + fcn[-9:]
	sio.savemat(savefile,{'centroid':centroid, 'area':area, 'score':score}) 

def remove_noise(img):
    th = filters.threshold_otsu(img)
    bin_img = img > th
    regions, num = ndimage.label(bin_img)
    areas = []
    for i in range(num):
        areas.append(np.sum(regions == i+1))
    img[regions != np.argmax(areas)+1] = 0
    return img 

def __call__(self, img_small):
        m = morphology.square(self.square_size)
        img_th = morphology.black_tophat(img_small, m)
        img_sob = abs(filters.sobel_v(img_th))
        img_closed = morphology.closing(img_sob, m)
        threshold = filters.threshold_otsu(img_closed)
        return img_closed > threshold 

def fallback(self,img):
    """Called when the spots are lost"""
    if self.safe_mode or self.fallback_mode:
      if self.fallback_mode:
        self.pipe.send("Fallback failed")
      else:
        self.pipe.send("[safe mode] Could not compute barycenter")
      self.fallback_mode = False
      return -1
    self.fallback_mode = True
    print("Loosing spot! Trying to reevaluate threshold...")
    self.thresh = threshold_otsu(img)
    return self.evaluate(img) 

def compute_multivarian_otsu(features):
    """ compute otsu individually over each sample dimension
    WARNING: this compute only localy  and since it does compare all
    combinations of orienting the asign for tight cases it may not decide

    :param ndarray features:
    :return list(bool):

    >>> np.random.seed(0)
    >>> fts = np.row_stack([np.random.random((5, 3)) - 1,
    ...                     np.random.random((5, 3)) + 1])
    >>> fts[:, 1] = - fts[:, 1]
    >>> compute_multivarian_otsu(fts).astype(int)
    array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1])
    """
    ys = np.zeros(features.shape)
    for i in range(features.shape[-1]):
        thr = filters.threshold_otsu(features[:, i])
        asign = features[:, i] > thr
        if i > 0:
            m = np.mean(ys[:, :i], axis=1)
            d1 = np.mean(np.abs(asign - m))
            d2 = np.mean(np.abs(~asign - m))
            # check if for this dimension it wount be better to swap it
            if d2 < d1:
                asign = ~asign
        ys[:, i] = asign
    y = np.mean(ys, axis=1) > 0.5
    return y


# def estim_class_model_gmm(features, nb_classes, init='kmeans'):
#     """ from all features estimate Gaussian Mixture Model and assuming
#     each cluster is a single class compute probability that each feature
#     belongs to each class
#
#     :param [[float]] features: list of features per segment
#     :param int nb_classes: number of classes
#     :return [[float]]: probabilities that each feature belongs to each class
#     """
#     logging.debug('estimate GMM for all given features %s and %i component',
#                   repr(features.shape), nb_classes)
#     # http://scikit-learn.org/stable/modules/generated/sklearn.mixture.GMM.html
#     mm = mixture.GMM(n_components=nb_classes, covariance_type='full', n_iter=999)
#     if init == 'kmeans':
#         # http://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html
#         kmeans = cluster.KMeans(n_clusters=nb_classes, init='k-means++', n_jobs=-1)
#         y = kmeans.fit_predict(features)
#         mm.fit(features, y)
#     else:
#         mm.fit(features)
#     logging.info('compute probability of each feature to all component')
#     return mm 

def find_wl_contours1(im_ndwi, cloud_mask, im_ref_buffer):
    """
    Traditional method for shoreline detection using a global threshold.
    Finds the water line by thresholding the Normalized Difference Water Index 
    and applying the Marching Squares Algorithm to contour the iso-value 
    corresponding to the threshold.

    KV WRL 2018

    Arguments:
    -----------
    im_ndwi: np.ndarray
        Image (2D) with the NDWI (water index)
    cloud_mask: np.ndarray
        2D cloud mask with True where cloud pixels are
    im_ref_buffer: np.array
        Binary image containing a buffer around the reference shoreline

    Returns:    
    -----------
    contours_wl: list of np.arrays
        contains the coordinates of the contour lines

    """

    # reshape image to vector
    vec_ndwi = im_ndwi.reshape(im_ndwi.shape[0] * im_ndwi.shape[1])
    vec_mask = cloud_mask.reshape(cloud_mask.shape[0] * cloud_mask.shape[1])
    vec = vec_ndwi[~vec_mask]
    # apply otsu's threshold
    vec = vec[~np.isnan(vec)]
    t_otsu = filters.threshold_otsu(vec)
    # use Marching Squares algorithm to detect contours on ndwi image
    im_ndwi_buffer = np.copy(im_ndwi)
    im_ndwi_buffer[~im_ref_buffer] = np.nan
    contours = measure.find_contours(im_ndwi_buffer, t_otsu)

    # remove contours that contain NaNs (due to cloud pixels in the contour)
    contours_nonans = []
    for k in range(len(contours)):
        if np.any(np.isnan(contours[k])):
            index_nan = np.where(np.isnan(contours[k]))[0]
            contours_temp = np.delete(contours[k], index_nan, axis=0)
            if len(contours_temp) > 1:
                contours_nonans.append(contours_temp)
        else:
            contours_nonans.append(contours[k])
    contours = contours_nonans

    return contours