# -*- coding: utf-8 -*-
from __future__ import division
from past.utils import old_div
import cv2
import numpy as np
from ..lib.arrayops import (normsigmoid, normalize, Bandpass, Bandstop,
findminima, findmaxima, find_near, smooth, getOtsuThresh, convexityRatio,
filterFactory, brightness, background, thresh_biggestCnt, contours2mask,
findContours)
def _getBrightAlpha(backgray, foregray, window=None):
"""
Get alpha transparency for merging foreground to
background gray image according to brightness.
(This is a test and not intended for practical use)
:param backgray: background image.
:param foregray: foreground image.
:param window: window used to customizing alfa. It can be a binary or alpha mask,
values go from 0 for transparency to any value where the maximum is visible
i.e a window with all the same values does nothing. A binary mask can be used,
where 0 is transparent and 1 is visible. If not window is given alfa is not
altered and the intended alpha is returned.
:return: alfa mask
"""
# this method was obtained for a particular problem, change to an
# automated one
backmask = normalize(normsigmoid(backgray, 10, 180) + normsigmoid(backgray, 3.14, 192) +
normsigmoid(backgray, -3.14, 45))
foremask = normalize(normsigmoid(foregray, -1, 242)
* normsigmoid(foregray, 3.14, 50))
foremask = normalize(foremask * backmask)
foremask[foremask > 0.9] = 2.0
ksize = (21, 21)
foremask = normalize(cv2.blur(foremask, ksize))
if window is not None:
foremask *= normalize(window) # ensures that window is normilized to 1
return foremask
[docs]def get_beta_params_hist(P):
"""
Automatically find parameters for bright alpha masks
using a histogram analysis method.
:param P: gray image
:return: beta1 for minimum valley left of body, beta2 for brightest valley right of body
where the body starts at the tallest peak in the histogram.
"""
hist_P, bins = np.histogram(P.flatten(), 256, [0, 256])
window = "hamming" # 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
window_len = 51
hist_PS = smooth(hist_P, window_len, window=window, correct=True)
otsu = getOtsuThresh(hist_P) # Otsu value
minima = findminima(hist=hist_PS) # minima values
maxima = findmaxima(hist=hist_PS) # maxima values
if minima != [] and maxima != []:
data_body = find_near(maxima, thresh=otsu, side="right")
beta1 = find_near(minima, thresh=data_body, side="left")
beta2 = find_near(minima, thresh=data_body, side="right")
else: # case where histogram is uniform and does not have minima or maxima values
beta1 = beta2 = otsu
if beta2 == beta1:
beta2 += 1 # prevents overlapping
return beta1, beta2
[docs]def get_bright_alpha(backgray, foregray, window=None):
"""
Get alpha transparency for merging foreground to
background gray image according to brightness.
:param backgray: background image. (as float)
:param foregray: foreground image. (as float)
:param window: window used to customizing alfa. It can be a binary or alpha mask,
values go from 0 for transparency to any value where the maximum is visible
i.e a window with all the same values does nothing. A binary mask can be used,
where 0 is transparent and 1 is visible. If not window is given alfa is not
altered and the intended alpha is returned.
:return: alfa mask
"""
backmask = Bandpass(3, *get_beta_params_hist(backgray)
)(backgray) # beta1 = 50, beta2 = 190
foremask = Bandpass(3, *get_beta_params_hist(foregray)
)(foregray) # beta1 = 50, beta2 = 220
foremask = normalize(foremask * backmask)
if window is not None:
foremask *= normalize(window) # ensures that window is normilized to 1
return foremask
[docs]def get_beta_params_Otsu(P):
"""
Automatically find parameters for alpha masks using Otsu threshold value.
:param P: gray image
:return: beta1 for minimum histogram value, beta2 for Otsu value
"""
# process histogram for uint8 gray image
if P.any(): # if array is not empty
hist, bins = np.histogram(P.flatten(), 256)
# get Otsu thresh as beta2
beta2 = bins[getOtsuThresh(hist)]
return np.min(P), beta2 # beta1, beta2
return 0.0, 1.0
[docs]def get_layered_alpha(back, fore):
"""
Get bright alpha mask (using Otsu method)
:param back: BGR background image
:param fore: BGR foreground image
:return: alpha mask
"""
# find retinal area and its alpha
mask_back, alpha_back = retinal_mask(back, biggest=True, addalpha=True)
mask_fore, alpha_fore = retinal_mask(fore, biggest=True, addalpha=True)
# convert uint8 to float
backgray = brightness(back).astype(float)
foregray = brightness(fore).astype(float)
# scale from 0-1 to 0-255
backm = alpha_back * 255
forem = alpha_fore * 255
# get alpha masks fro background and foreground
backmask = Bandstop(
3, *get_beta_params_Otsu(backm[mask_back.astype(bool)]))(backgray)
foremask = Bandpass(
3, *get_beta_params_Otsu(forem[mask_fore.astype(bool)]))(foregray)
# merge masks
alphamask = normalize(foremask * backmask * (old_div(backm, 255.)))
return alphamask
[docs]def retina_markers_thresh(P):
"""
Retinal markers thresholds to find background,
retinal area and optic disc with flares based
in the histogram.
:param P: gray image
:return: min,b1,b2,max
where:
black background < min
b1 > retina < b2
flares > max
"""
# calculate histogram
#hist_P = histogram(P)[0]
hist_P, bins = np.histogram(P.flatten(), 256, [0, 256])
# filter histogram to smooth it.
window = "hamming" # 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
window_len = 51
# it smooths and expands the histogram
hist_PS = smooth(hist_P, window_len, window=window, correct=True)
#from scipy.signal import savgol_filter
#hist_PS = savgol_filter(hist_P, window_len, polyorder = 3)
otsu = getOtsuThresh(hist_P) # Otsu value
minima = findminima(hist=hist_PS) # minima values
maxima = findmaxima(hist=hist_PS) # maxima values
# INITIAL DATA
data_min_left = P.min() # initial minimum value
data_min = find_near(maxima, thresh=data_min_left)
#data_min_right = find_near(minima, thresh=data_min, side="right")
data_max = P.max()
data_max_left = find_near(minima, thresh=data_max, side="left")
# SECONDARY DATA
data_body = find_near(maxima, thresh=otsu, side="right")
data_body_left = find_near(minima, thresh=data_body, side="left")
#data_body_right = find_near(minima,thresh=data_body,side="right")
return data_min, data_body_left, data_body, data_max_left
[docs]def find_optic_disc_watershed(img, P):
"""
Find optic disk in image using a watershed method.
:param img: BGR image
:param P: gray image
:return: optic_disc, Crs, markers, watershed
"""
assert img is not None and P is not None
#fore = cv2.cvtColor(P,cv2.COLOR_GRAY2BGR)
# create watershed
data_min, data_body_left, data_body, data_max_left = retina_markers_thresh(
P)
watershed = np.zeros_like(P).astype("int32")
mk_back, mk_body, mk_flare = 1, 2, 3
# background FIXMED use P.min() and aproaching to local maxima
watershed[P <= data_min] = mk_back
watershed[np.bitwise_and(P > data_body_left, P <
data_body)] = mk_body # main body
# find bright objects
flares_thresh = (P >= data_max_left).astype(np.uint8)
contours, hierarchy = findContours(
flares_thresh.copy(), cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
mks_flare = []
for i, cnt in enumerate(contours):
index = mk_flare + i
mks_flare.append(index)
# Flares. this can be used approaching
watershed[contours2mask(cnt, shape=watershed.shape)] = index
# to local maxima, but brightest areas are almost
# always saturated so no need to use it
markers = watershed.copy()
# apply watershed to watershed
# FIXME perhaps the function should be cv2.floodFill?
cv2.watershed(img, watershed)
# different classification algorithms could be used using watershed
contours_flares = []
for mk_flare in mks_flare:
brightest = np.uint8(watershed == mk_flare)
contours, hierarchy = findContours(
brightest, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
contours_flares.extend(contours)
Crs = [(i, j) for i, j in [(convexityRatio(cnt), cnt)
for cnt in contours_flares] if i != 0] # convexity ratios
Crs.sort(reverse=True, key=lambda x: x[0])
candidate = Crs[-1]
ellipse = cv2.fitEllipse(candidate[1])
optic_disc = np.zeros_like(P)
cv2.ellipse(optic_disc, ellipse, 1, -1) # get elliptical ROI
return optic_disc, Crs, markers, watershed
[docs]def layeredfloods(img, gray=None, backmask=None, step=1, connectivity=4, weight=False):
"""
Create an alpha mask from an image using a weighted layered flooding algorithm,
:param img: BGR image
:param gray: Gray image
:param backmask: background mask
:param step: step to increase upDiff in the floodFill algorithm. If weight is True
step also increases the weight of the layers.
:param connectivity: pixel connectivity of 4 or 8 to use in the floodFill algorithm
:param weight: Increase progressively the weight of the layers using the step parameter.
:return: alpha mask
"""
if gray is None:
if len(img.shape) > 2:
gray = brightness(img)
else:
gray = img
# initialization of flags for flooding
flags = connectivity | cv2.FLOODFILL_FIXED_RANGE | cv2.FLOODFILL_MASK_ONLY
h, w = gray.shape # get image shape
# pad mask to add 1 layer of pixels
# this is needed in the flood operation
mask = np.zeros((h + 2, w + 2), np.uint8)
# create the background mask
if backmask is None:
backmask = background(gray)
mask_background = np.zeros((h + 2, w + 2), np.uint8)
mask_background[1:-1, 1:-1] = backmask
# get area of background mask
area_background = np.sum(mask_background)
# low and high limits for flooding
lo, hi = 255, 0
# first seed point of the minimum color
xs, ys = np.where(gray == np.min(gray))
seed_pt = xs[0], ys[0]
# flood adding overal layers starting from background
area = 0
hi = step
all = np.zeros_like(mask, np.float32)
while hi <= 255:
mask[:] = 0
cv2.floodFill(img, mask, seed_pt, (255, 255, 255),
(lo,) * 3, (hi,) * 3, flags)
if weight and area >= area_background:
all += mask.astype("float") * hi
elif not weight and area >= area_background:
all += mask.astype("float")
else:
area = np.sum(np.bitwise_and(mask, mask_background))
hi += step
all = all[1:-1, 1:-1]
all /= all.max()
return all
[docs]def retinal_mask(img, biggest=False, addalpha=False):
"""
Obtain the mask of the retinal area in an image.
For a simpler and lightweight algorithm see :func:`retinal_mask_watershed`.
:param img: BGR or gray image
:param biggest: True to return only biggest object
:param addalpha: True to add additional alpha mask parameter
:return: if addalpha:
binary mask, alpha mask
else:
binary mask
"""
# pad image to add 1 layer of pixels
# this is used to correctly flood all the background
if len(img.shape) > 2:
h, w, c = img.shape
imgP = np.zeros((h + 2, w + 2, 3), np.uint8)
imgP[1:-1, 1:-1, :] = img
P = brightness(imgP)
else:
h, w = img.shape
imgP = np.zeros((h + 2, w + 2), np.uint8)
imgP[1:-1, 1:-1] = img
P = imgP
# get flooded alpha mask
mask_alpha = 1 - layeredfloods(imgP, gray=P)[1:-1, 1:-1]
hist, bins = np.histogram(mask_alpha.flatten(), 256)
thresh = bins[getOtsuThresh(hist)]
mask_binary = (mask_alpha > thresh).astype(np.uint8)
if biggest: # process to give only the biggest object
a = thresh_biggestCnt(mask_binary)
mask_binary = np.zeros_like(mask_binary)
cv2.drawContours(mask_binary, [a], 0, 1, -1) # get current fillet ROI
if addalpha:
#if biggest: mask_alpha[mask_binary==0] = 0
return mask_binary, mask_alpha
return mask_binary
[docs]def retinal_mask_watershed(img, parameters=(10, 30, None), addMarkers=False):
"""
Quick and simple watershed method to obtain the mask of the retinal area in an image.
For a more robust algorithm see :func:`retinal_mask`.
:param img: BGR or gray image
:param parameters: tuple of parameters to pass to :func:`filterFactory`
:param addMarkers: True to add additional Marker mask. It contains 0 for unknown
areas, 1 for background and 2 for retinal area.
:return: if addMarkers:
binary mask, Markers mask
else:
binary mask
"""
if parameters is not None:
myfilter = filterFactory(*parameters) # alfa,beta1,beta2
img = (myfilter(img.astype("float")) *
255).astype("uint8") # *fore.astype("float")
P = brightness(img) # get scaled image brightness
thresh, sure_bg = cv2.threshold(
P, 0, 1, cv2.THRESH_BINARY + cv2.THRESH_OTSU) # obtain over threshold
thresh, sure_fg = cv2.threshold(P, thresh + 10, 1, cv2.THRESH_BINARY)
markers = np.ones_like(sure_fg).astype("int32") # make background markers
markers[sure_bg == 1] = 0 # mark unknown markers
markers[sure_fg == 1] = 2 # mark sure object markers
cv2.watershed(img, markers) # get watershed on markers
thresh, mask = cv2.threshold(markers.astype(
"uint8"), 1, 1, cv2.THRESH_BINARY) # get binary image of contour
if addMarkers:
return mask, thresh
return mask
