Python skimage.filters 模块，threshold_otsu() 实例源码

def convert_new(fname, target_size):
    print('Processing image: %s' % fname)
    img = Image.open(fname)
    blurred = img.filter(ImageFilter.BLUR)
    ba = np.array(blurred)
    ba_gray = rgb2gray(ba)
    val = filters.threshold_otsu(ba_gray)
    # foreground = (ba_gray > val).astype(np.uint8)
    foreground = closing(ba_gray > val, square(3))
    # kernel = morphology.rectangle(5, 5)
    # foreground = morphology.binary_dilation(foreground, kernel)
    labels = measure.label(foreground)
    properties = measure.regionprops(labels)
    properties = sorted(properties, key=lambda p: p.area, reverse=True)
    # draw_top_regions(properties, 3)
    # return ba
    bbox = properties[0].bbox
    bbox = (bbox[1], bbox[0], bbox[3], bbox[2])
    cropped = img.crop(bbox)
    resized = cropped.resize([target_size, target_size])
    return np.array(resized)

def convert_new_regions(fname, target_size):
    print('Processing image: %s' % fname)
    img = Image.open(fname)
    blurred = img.filter(ImageFilter.BLUR)
    ba = np.array(blurred)
    ba_gray = rgb2gray(ba)
    val = filters.threshold_otsu(ba_gray)
    # foreground = (ba_gray > val).astype(np.uint8)
    foreground = closing(ba_gray > val, square(3))
    # kernel = morphology.rectangle(5, 5)
    # foreground = morphology.binary_dilation(foreground, kernel)
    labels = measure.label(foreground)
    properties = measure.regionprops(labels)
    properties = sorted(properties, key=lambda p: p.area, reverse=True)
    draw_top_regions(properties, 3)
    return ba

def convert(fname, target_size):
    # print('Processing image: %s' % fname)
    img = Image.open(fname)
    blurred = img.filter(ImageFilter.BLUR)
    ba = np.array(blurred)
    ba_gray = rgb2gray(ba)
    val = filters.threshold_otsu(ba_gray)
    # foreground = (ba_gray > val).astype(np.uint8)
    foreground = closing(ba_gray > val, square(3))
    # kernel = morphology.rectangle(5, 5)
    # foreground = morphology.binary_dilation(foreground, kernel)
    labels = measure.label(foreground)
    properties = measure.regionprops(labels)
    properties = sorted(properties, key=lambda p: p.area, reverse=True)
    # draw_top_regions(properties, 3)
    # return ba
    bbox = properties[0].bbox
    bbox = (bbox[1], bbox[0], bbox[3], bbox[2])
    cropped = img.crop(bbox)
    resized = cropped.resize([target_size, target_size])
    return resized

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 otsu_thresholding(im_float):

    """
        apply otsu thresholding on the whole slide

        im_float = input image as float type 

        return otsued image with gray, otsu threshold value 
    """

    print("threshold_otsu\n")
    threshold_global_Otsu = threshold_otsu(im_float)

    #thresholding 
    im_bool= (im_float > threshold_global_Otsu)
    im_int = im_bool.astype(float)
    print im_int*255
    return im_int*255, threshold_global_Otsu

def otsu_thresholding(im_float):

    """
        apply otsu thresholding on the whole slide

        im_float = input image as float type 

        return otsued image with gray, otsu threshold value 
    """

    print("threshold_otsu\n")
    threshold_global_Otsu = threshold_otsu(im_float)

    #thresholding 
    im_bool= (im_float > threshold_global_Otsu)
    im_int = im_bool.astype(float)
    print im_int*255
    return im_int*255, threshold_global_Otsu

def get_real_images(paths):
    real_images = []
    for path in paths:
        # Calculate a threshold to do image binarization, all colors at every pixel will be translated to number 0(white) or 1(black)
        camera = io.imread(path)
        val = filters.threshold_otsu(camera)
        result = (camera < val)*1.0
        real_images.append(result)
    np_images = numpy.array(real_images)
    np_images = np_images.reshape(np_images.shape[0], np_images.shape[1] * np_images.shape[2])
    return np_images

def to_binary_otsu(img, invert=False):
    if img.dtype == np.bool:
        img = np.array(img, 'uint8')
    if img.max() == img.min():
        if img.min() == 1:
            return np.array(img * 255, 'uint8')
        else:
            return np.array(img, 'uint8')
    else:
        t = threshold_otsu(img)
    img[img <= t] = 255 if invert else 0
    img[img > t] = 0 if invert else 255
    return np.array(img, 'uint8')

def otsu_threshold(im):
    """
    :param im:
    :return:
    """
    val = filters.threshold_otsu(im.get_data())
    return set_new_data(im, im.get_data() * (im.get_data() > val))

def binarisation(src_image):
    if len(src_image.shape) == 3:
        image = (src_image.sum(axis=2) / 3).astype('ubyte')
    else:
        image = src_image
    thresh = threshold_otsu(image)
    binary = (image > thresh).astype('ubyte')
    binary1 = 1 - binary
    im = 255 - np.multiply(255 - image, binary1)
    block_size = 35
    binary = threshold_adaptive(image, block_size, offset=20)
    binary = binary.astype('ubyte')
    return binary

def binarisation(src_image):
    if len(src_image.shape) == 3:
        image = (src_image.sum(axis=2) / 3).astype('ubyte')
    else:
        image = src_image
    thresh = threshold_otsu(image)
    binary = (image > thresh).astype('ubyte')
    binary1 = 1 - binary
    im = 255 - np.multiply(255 - image, binary1)
    block_size = 35
    binary = threshold_adaptive(image, block_size, offset=20)
    binary = binary.astype('ubyte')
    return binary

# 100

def binarisation(src_image):
    if len(src_image.shape) == 3:
        image = (src_image.sum(axis=2) / 3).astype('ubyte')
    else:
        image = src_image
    thresh = threshold_otsu(image)
    binary = (image > thresh).astype('ubyte')
    binary1 = 1 - binary
    im = 255 - np.multiply(255 - image, binary1)
    block_size = 35
    binary = threshold_adaptive(im, block_size, offset=20)
    binary = binary.astype('ubyte')
    return binary

def get_candidate_symbol_regions(image, text_regions, updated_width, updated_height):
    img = skimage.io.imread(image.name)[:, :, :3]
    if not (updated_height == len(img) and updated_width == len(img[0])):
        img = skimage.transform.resize(img, (updated_height, updated_width))

    symbol_regions = dict()
    for x, y, w, h in text_regions:
        text_region_image = img[y: y + h, x: x + w]
        text_region_image_width = len(text_region_image[0])
        text_region_image_height = len(text_region_image)

        text_region_gray_image = skimage.color.rgb2gray(text_region_image)
        text_region_binary_image = image <= threshold_otsu(text_region_gray_image)

        temp = TemporaryFile(".png")
        skimage.io.imsave(temp.name, text_region_binary_image)
        text_region_binary_image = skimage.io.imread(temp.name)

        text_region_blurred_image = gaussian_filter(text_region_binary_image, sigma=3.5)
        text_region_blobs = text_region_blurred_image > text_region_blurred_image.mean()

        text_region_labels = skimage.morphology.label(text_region_blobs, neighbors=4)

        symbol_blobs = ndimage.find_objects(text_region_labels)
        candidate_symbol_regions = set()

        for c1, c2 in symbol_blobs:
            if (c2.stop - c2.start) * c1.stop - c1.start > (text_region_image.shape[0] * text_region_image.shape[1]) * (0.026):
                if (c2.stop - c2.start) * c1.stop - c1.start < (text_region_image.shape[0] * text_region_image.shape[1]) * (0.90):
                    candidate_symbol_regions.add(
                        (c2.start, c1.start, c2.stop - c2.start, c1.stop - c1.start))

        symbol_regions[str((x, y, w, h))] = dict()
        symbol_regions[str((x, y, w, h))]["image"] = text_region_image
        symbol_regions[str((x, y, w, h))]["regions"] = candidate_symbol_regions
        symbol_regions[str((x, y, w, h))]["width"] = text_region_image_width
        symbol_regions[str((x, y, w, h))]["height"] = text_region_image_height

    return symbol_regions

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 find_ellipse(image, mode='ellipse_aligned', min_length=24):
    """ Find bright ellipse contours on a black background.

    This routine thresholds the image (using the Otsu threshold), finds the
    longest contour and fits an ellipse to this contour.

    Parameters
    ----------
    image : ndarray, 2d
    mode : {'ellipse', 'ellipse_aligned', 'circle'}
        'ellipse' or None finds an arbitrary ellipse (default)
        'circle' finds a circle
        'ellipse_aligned' finds an ellipse with its axes aligned along [x y] axes
    min_length : number, optional
        minimum length of the ellipse contour, in pixels. Default 24.

    Returns
    -------
    yr, xr, yc, xc when dimension order was y, x (most common)
    xr, yr, xc, yc when dimension order was x, y
    """
    assert image.ndim == 2
    # Threshold the image
    thresh = threshold_otsu(image)
    binary = image > thresh

    # Find the contours of 0.5 value. For a thresholded ellipse contour, this
    # likely finds 2 contours: the inner and the outer.
    contours = find_contours(binary, 0.5, fully_connected='high')
    if len(contours) == 0:
        raise ValueError('No contours found')

    # Eliminate short contours
    contours = [c for c in contours if len(c) >= min_length]

    # fit circles to the rest, keep the one with lowest residual deviation
    result = [np.nan] * 4
    residual = None
    for c in contours:
        try:
            (xr, yr), (xc, yc), _ = fit_ellipse(c.T, mode=mode)
            if np.any(np.isnan([xr, yr, xc, yc])):
                continue
            x, y = c.T
            r = np.sum((((xc - x)/xr)**2 + ((yc - y)/yr)**2 - 1)**2)/len(c)
            if residual is None or r < residual:
                result = xr, yr, xc, yc
                residual = r
        except np.linalg.LinAlgError:
            pass

    return result

def threshold(image, *, sigma=0., radius=0, offset=0.,
              method='sauvola', smooth_method='Gaussian'):
    """Use scikit-image filters to "intelligently" threshold an image.

    Parameters
    ----------
    image : array, shape (M, N, ...[, 3])
        Input image, conformant with scikit-image data type
        specification [1]_.
    sigma : float, optional
        If positive, use Gaussian filtering to smooth the image before
        thresholding.
    radius : int, optional
        If given, use local median thresholding instead of global.
    offset : float, optional
        If given, reduce the threshold by this amount. Higher values
        result in fewer pixels above the threshold.
    method: {'sauvola', 'niblack', 'median'}
        Which method to use for thresholding. Sauvola is 100x faster, but
        median might be more accurate.
    smooth_method: {'Gaussian', 'TV', 'NL'}
        Which method to use for smoothing. Choose from Gaussian smoothing,
        total variation denoising, and non-local means denoising.

    Returns
    -------
    thresholded : image of bool, same shape as `image`
        The thresholded image.

    References
    ----------
    .. [1] http://scikit-image.org/docs/dev/user_guide/data_types.html
    """
    if sigma > 0:
        if smooth_method.lower() == 'gaussian':
            image = filters.gaussian(image, sigma=sigma)
        elif smooth_method.lower() == 'tv':
            image = restoration.denoise_tv_bregman(image, weight=sigma)
        elif smooth_method.lower() == 'nl':
            image = restoration.denoise_nl_means(image,
                                             patch_size=round(2 * sigma))
    if radius == 0:
        t = filters.threshold_otsu(image) + offset
    else:
        if method == 'median':
            footprint = hyperball(image.ndim, radius=radius)
            t = ndi.median_filter(image, footprint=footprint) + offset
        elif method == 'sauvola':
            w = 2 * radius + 1
            t = threshold_sauvola(image, window_size=w, k=offset)
        elif method == 'niblack':
            w = 2 * radius + 1
            t = threshold_niblack(image, window_size=w, k=offset)
        else:
            raise ValueError('Unknown method %s. Valid methods are median,'
                             'niblack, and sauvola.' % method)
    thresholded = image > t
    return thresholded

def skeleton_from_image_discrete(image, sigma = 1, absolute_threshold = None, threshold_factor = 0.95, with_head_tail = False, verbose = False):
  """Detect skeleton points of wormshape in image

  Arguments:
    image (array): the image to detect venterline of worm from
    sigma (float or None): width of Gaussian smoothing on image, if None use raw image
    absolute_threshold (float or None): if set use this as the threshold, if None the threshold is set via Otsu
    threshold_level (float): in case the threshold is determined by Otsu multiply by this factor
    verbose (bool): plot results

  Returns:
    array (nx2): unsorted skeleton points
  """
  ### smooth image
  if sigma is not None:
    imgs = filters.gaussian_filter(np.asarray(image, float), sigma);
  else:
    imgs = image;

  ### get worm foreground
  if absolute_threshold is not None:
    level = absolute_threshold;
  else:
    level = threshold_factor * threshold_otsu(imgs);
  imgth = imgs < level;  

  ### skeletonize
  skel = skeletonize(imgth);
  y,x = np.where(skel);

  if verbose:
    plt.imshow(imgth, interpolation = 'none');
    plt.scatter(x,y,c = 'k', s = 40);

  if with_head_tail:
    # find end points:
    adj = skeleton_to_adjacency(skel);
    nhs = np.array([len(v) for v in adj.values()]);
    ht = np.where(nhs == 1)[0];

    if verbose:
      xy = np.vstack([x,y]).T;
      if ht.shape[0] > 0:
        xyht = xy[ht];
        plt.scatter(xyht[:,0], xyht[:,1], s = 60, c = 'r');

    return np.vstack([x,y]).T, ht;
  else:
    return np.vstack([x,y]).T;

def main():
    """Load image, calculate optimal threshold, binarize, plot."""
    # load image
    img = data.coins()
    height, width = img.shape
    nb_pixels = height * width

    # precalculate some values for speedup
    # average pixel value
    g_avg = np.average(img)
    # P(pixel-value), i.e. #pixels-with-value / #all-pixels
    p_g = [0] * 256
    for g in range(0, 256):
        p_g[g] = np.sum(img == g) / nb_pixels

    # Otsu method
    # calculations are based on standard formulas
    q_best = None
    threshold_best = None
    img_bin_best = None
    # iterate over all possible thresholds
    for t in range(1, 255):
        img_bin = np.zeros(img.shape)
        img_bin[img >= t] = 1

        p1 = np.sum(img_bin) / nb_pixels
        p0 = 1 - p1

        g0 = np.average(img[img_bin == 0]) if np.sum(img[img_bin == 0]) > 0 else 0
        g1 = np.average(img[img_bin == 1]) if np.sum(img[img_bin == 1]) > 0 else 0

        var0 = sum([(g-g0)**2 * p_g[g] for g in range(0, t+1)])
        var1 = sum([(g-g1)**2 * p_g[g] for g in range(t+1, 256)])

        var_between = p0 * (g0 - g_avg)**2 + p1 * (g1 - g_avg)**2
        var_inner = p0 * var0**2 + p1 * var1**2

        # q is the relation of variance between classes and variance within classes
        q = var_between / var_inner if var_inner > 0 else 0

        print(t, p0, p1, g0, g1, g_avg, var_between, var_inner, q)
        if q_best is None or q_best < q:
            q_best = q
            threshold_best = t
            img_bin_best = img <= t

    # ground truth, based on scikit-image
    gt_tresh = skifilters.threshold_otsu(img)
    ground_truth = img <= gt_tresh

    # plot
    util.plot_images_grayscale(
        [img, img_bin_best, ground_truth],
        ["Image", "Otsu", "Otsu (Ground Truth)"]
    )