Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
207 lines (167 sloc) 8.71 KB
import numpy as np
import scipy as sp
from scipy import ndimage as ndi
from skimage import color, exposure, feature
from skimage.filters import threshold_otsu
from skimage.morphology import remove_small_objects, disk, square, opening, dilation
from skimage.measure import find_contours, approximate_polygon
from skimage import segmentation
from skimage.future import graph
from skimage import transform as tf
from skimage.draw import polygon
def find_largest(img, threshold = 'opening'):
"""Function to find largest bright object in an image.
Workflow:
- Convert to binary
- Find and retain only largest object in binary image
- Fill holes
- Apply opening and dilation to remove minutiae
Current options to convert to binary:
- 'opening': applies global Otsu thresholding after a minor morphological opening (default)
- 'global': applies global Otsu thresholding
- 'adaptive': applies adaptive threshold after gaussian smoothing
- 'tips': eliminates both ~black and ~nearly white backgrounds
- 'rag': applies RAG thresholding, followed by global Otsu. Slowest method.
Coming soon:
- 'textural': eliminates backgrounds of (any) constant colour based on GLCM or binary patterns
- 'bilateral': applies adaptive threshold, but after bilateral smoothing"""
# Gray scale conversion and contrast stretching
gry = color.rgb2gray(img);
p2, p95 = np.percentile(gry, (2, 95))
if threshold == 'rag':
labels1 = segmentation.slic(img, compactness=50, n_segments=500)
out1 = color.label2rgb(labels1, img, kind='avg')
thresh = min(np.shape(out1)[:1])/10 # threshold = 1/10th of the smallest side of the image
# with a large value we get few colours
g = graph.rag_mean_color(img, labels1)
labels2 = graph.cut_threshold(labels1, g, thresh)
out2 = color.label2rgb(labels2, img, kind='avg')
out2 = color.rgb2gray(out2);
global_thresh = threshold_otsu(out2)
binary = out2 < global_thresh
elif threshold == 'tips':
rescale = exposure.rescale_intensity(gry, in_range=(p2, p95))
binary = np.logical_and(color.rgb2gray(rescale) > 0.05, color.rgb2gray(rescale) < 0.95)
elif threshold == 'adaptive':
gry = sp.ndimage.filters.gaussian_filter(color.rgb2gray(gry),7)
rescale = exposure.rescale_intensity(gry, in_range=(p2, p95))
thresh_sauvola = threshold_sauvola(rescale, window_size=99)
binary = rescale < thresh_sauvola
elif threshold == 'global':
rescale = exposure.rescale_intensity(gry, in_range=(p2, p95))
global_thresh = threshold_otsu(rescale)
binary = rescale < global_thresh
else:
rescale = 1-opening(1-(exposure.rescale_intensity(gry, in_range=(p2, p95))), disk(2))
global_thresh = threshold_otsu(rescale)
binary = rescale < global_thresh
# Detect largest bright element in the binary image. Making the assumption it would be the map.
# Eliminate everything else (text, colorbar, holes, ...).
# Label all white objects (made up of ones)
label_objects, nb_labels = ndi.label(binary) # ndimage.label actually labels 0 (background) as 0 and then
# labels every nonzero object as 1, 2, ... n.
# Calculate every labeled object's size.
# np.bincount ignores whether input is an image or another type of array.
# It just calculates the binary sizes, including for the 0 (background).
sizes = np.bincount(label_objects.ravel())
sizes[0] = 0 # This sets the size of the background to 0 so that if it happened to be larger than
# the largest white object it would not matter
# Keep only largest object
binary_objects = remove_small_objects(binary, max(sizes))
# Remove holes from it (black regions inside white object)
binary_holes = ndi.morphology.binary_fill_holes(binary_objects)
# This may help remove minutiae on the outside (tick marks and tick labels)
binary_mask = opening(binary_holes, disk(7))
return binary_mask
def ordered(points):
"""Function to sort corners based on angle from centroid.
Modified from: http://stackoverflow.com/a/31235064/1034648"""
x = points[:,0]
y = points[:,1]
cx = np.mean(x)
cy = np.mean(y)
a = np.arctan2(y - cy, x - cx)
order = a.ravel().argsort()
x = x[order]
y = y[order]
return np.vstack([x,y])
def rectify_seismic(img, binary_mask):
"""Function to warp to a rectangle the area in the input img defined by binary_mask
It returns the warped area as an image"""
# Find mask contour, approximate it with a quadrilateral, find and sort corners
contour = np.squeeze(find_contours(binary_mask, 0))
coords = approximate_polygon(contour, tolerance=50)
# sort the corners with the exception of the last one (repetition of first corner)
sortedCoords = ordered(coords[:-1]).T
# Define size of output image based on largest width and height in the input
w1 = np.sqrt(((sortedCoords[0, 1]-sortedCoords[3, 1])**2)+((sortedCoords[0, 0]-sortedCoords[3, 0])**2))
w2 = np.sqrt(((sortedCoords[1, 1]-sortedCoords[2, 1])**2)+((sortedCoords[1, 0]-sortedCoords[2, 0])**2))
h1 = np.sqrt(((sortedCoords[0, 1]-sortedCoords[1, 1])**2)+((sortedCoords[0, 0]-sortedCoords[1, 0])**2))
h2 = np.sqrt(((sortedCoords[3, 1]-sortedCoords[2, 1])**2)+((sortedCoords[3, 0]-sortedCoords[2, 0])**2))
w = max(int(w1), int(w2))
h = max(int(h1), int(h2))
# Define rectangular destination coordinates (homologous points) for warping
dst = np.array([[0, 0],
[h-1, 0],
[h-1, w-1],
[0, w-1]], dtype = 'float32')
# Estimate warping transform, apply to input image (mask portion), and output
dst[:,[0,1]] = dst[:,[1,0]]
sortedCoords[:,[0,1]] = sortedCoords[:,[1,0]]
tform = tf.ProjectiveTransform()
tform.estimate(dst,sortedCoords)
warped = tf.warp(img, tform, output_shape=(h-1, w-1))
return warped
def auto_canny(img, sigma = 0.33):
"""Zero-parameter, automatic Canny edge detection using scikit-image.
Original function from pyimagesearch: Zero-parameter, automatic Canny edge with with Python and OpenCV
www.pyimagesearch.com/2015/04/06/zero-parameter-automatic-canny-edge-detection-with-python-and-opencv"""
# compute the median of the single channel pixel intensities
v = np.median(img)
# apply automatic Canny edge detection using the computed median
lower = float(max(0.0, (1.0 - sigma) * v))
upper = float(min(1.0, (1.0 + sigma) * v))
edged = feature.canny(img, sigma, lower, upper)
# return the edged image
return edged
# Function to test if 4 points form a rectangle. Adapted to Python from:
# http://stackoverflow.com/questions/2303278/find-if-4-points-on-a-plane-form-a-rectangle/2304031#2304031
def is_rectangle(points):
x = points[:,0]
y = points[:,1]
cx = np.mean(x)
cy = np.mean(y)
dd1 = abs(cx-x[0]) + abs(cy-y[0])
dd2 = abs(cx-x[1]) + abs(cy-y[1])
dd3 = abs(cx-x[2]) + abs(cy-y[2])
dd4 = abs(cx-x[3]) + abs(cy-y[3])
return abs(dd1-dd2)<=10 and abs(dd1-dd3)<=10 and abs(dd1-dd4)<=10
# A better alternative may be to use angles between the semi-diagonals connecting centroid with 4 corners (sorted).
def remove_annotations(img, g = None):
"""Function to remove rectangular areas with annotations from an image.
Workflow:
- Blur a bit and find canny edges
- Find contours of edges and approximate
- Retain only approximated contours (rectangles) based on two conditions:
having 4 sides and having 4 equal semi-diagonals
- Convert retained contours to polygons
- Add to a binary mask
- Use mask to remove the rectangular areas in the input"""
if g is not None:
gry = sp.signal.medfilt2d(color.rgb2gray(img), g)
else:
gry = color.rgb2gray(img)
#
cn = auto_canny(color.rgb2gray(gry)) # canny edges with automatic parameters based on median pixel intensity
cny_dilated = dilation(cn, disk(2.5)) # dilate edges
mask = np.zeros(np.shape(cny_dilated), dtype=np.uint8)
for contour in find_contours(cny_dilated, 0):
coords = approximate_polygon(contour, tolerance=30)
if len(coords) == 5 and is_rectangle(ordered(coords[:-1]).T):
rr,cc = polygon(ordered(coords[:-1])[0], ordered(coords[:-1])[1])
mask[rr,cc]=1
dmask = dilation(mask, square(10))
image = img.copy()
for layer in range(image.shape[-1]):
image[np.where(dmask)] = np.NaN
return image
You can’t perform that action at this time.