In [ ]:
import javabridge
import bioformats
from matplotlib.widgets import RectangleSelector, Button
import numpy as np
import matplotlib.pyplot as plt
from skimage.transform import rescale, pyramid_gaussian, resize
%matplotlib osx
In [ ]:
# load image
path="/Users/david/Desktop/wing_overview.nd2"
ds = 2
layers = 3
In [ ]:
javabridge.start_vm(class_path=bioformats.JARS, run_headless=True)
img = bioformats.load_image(path, rescale=False)
pyr = [p for p in pyramid_gaussian(img, max_layer= layers, downscale = ds)]
In [ ]:
%matplotlib
layer_label = 3
In [ ]:
def add_selection_callback(ev):
add_selection_callback.accu.append((select_callback.x1, select_callback.x2, select_callback.y1, select_callback.y2))
print('added selection')
add_selection_callback.accu = []
def select_callback(eclick, erelease):
select_callback.x1, select_callback.y1 = eclick.xdata, eclick.ydata
select_callback.x2, select_callback.y2 = erelease.xdata, erelease.ydata
print('selected ' + str((select_callback.x1, select_callback.x2, select_callback.y1, select_callback.y2)))
select_callback.x1, select_callback.x2, select_callback.y1, select_callback.y2 = 0, 0, 0, 0
plt.imshow(pyr[layer_label])
rs = RectangleSelector(plt.gca(), select_callback,
drawtype='box', useblit=True,
button=[1, 3], # don't use middle button
minspanx=5, minspany=5,
spancoords='pixels',
interactive=True)
# add selection button
axnext = plt.axes([0.45, -0.002, 0.1, 0.075])
but = Button(axnext, "ADD" )
but.on_clicked(add_selection_callback)
plt.show()
In [ ]:
%matplotlib inline
layer_to_use = 2
cuts = []
for x1, x2, y1, y2 in add_selection_callback.accu:
# cut from full resolution
x1 *= ds**(layer_label - layer_to_use)
y1 *= ds**(layer_label - layer_to_use)
x2 *= ds**(layer_label - layer_to_use)
y2 *= ds**(layer_label - layer_to_use)
x1 = int(x1)
y1 = int(y1)
x2 = int(x2)
y2 = int(y2)
#plt.figure()
#plt.imshow(pyr[0][y1:y2, x1:x2])
cuts.append(pyr[layer_to_use][y1:y2, x1:x2])
In [ ]:
print('selection ({}, {}) -> ({}, {})'.format(x1, y1, x2, y2))
plt.imshow(pyr[layer_to_use][y1:y2, x1:x2])
In [ ]:
from skimage.feature import ORB, match_descriptors, plot_matches
from skimage.transform import AffineTransform
from skimage.measure import ransac
from sklearn.cluster import MeanShift, estimate_bandwidth, DBSCAN
plt.rcParams['figure.figsize'] = [12.0,12.0]
idx1 = 0
idx2 = 8
offx = 200
offy = 67
sizy = 256
sizx = 256
img1 = cuts[idx1]
#img1 = resize(cuts[idx1], (sizy, sizx))
#img2 = resize(cuts[idx2], (sizy, sizx))
#img2 = pyr[layer_to_use][offy:(offy + sizy), offx:(offx+sizx)]
img2 = pyr[layer_to_use][1000:5000, 1000:4000]
#img2 = pyr[layer_to_use]
orb1 = ORB(n_keypoints=100, harris_k=0.14, fast_threshold=0.03, fast_n=12)
orb = ORB(n_keypoints=2500, harris_k=0.14, fast_threshold=0.03, fast_n=12)
orb1.detect_and_extract(img1)
desc1 = orb1.descriptors
key1 = orb1.keypoints
orb.detect_and_extract(img2)
desc2 = orb.descriptors
key2 = orb.keypoints
plt.figure()
m = match_descriptors(desc1, desc2)
plot_matches(plt.gca(), img1, img2, key1, key2, m)
#bandwidth = estimate_bandwidth(key2, quantile=0.2, n_samples=500)
#ms = MeanShift(bandwidth=bandwidth, bin_seeding=True, cluster_all=True)
#ms.fit(key2)
ms = DBSCAN(eps=40, min_samples=50)
ms.fit(key2)
labels = ms.labels_
#cluster_centers = ms.cluster_centers_
print(len(set(labels)))
cents = []
for lab in [l for l in set(labels) if l >= 0]:
cp = key2[labels==lab]
if len(cp)>100:
continue
cent = np.apply_along_axis(np.mean, 0, cp)
cents.append(cent)
plt.figure()
m = match_descriptors(desc1, desc2[labels==lab])
plot_matches(plt.gca(), img1, img2, key1, key2[labels==lab], m)
d1 = key1[m[:,0],:]
d2 = key2[m[:,1],:]
if (m.shape[0] < minsamples):
continue
model, inlier = ransac((d1, d2), AffineTransform, minsamples, 10)
nInlier = sum(inlier)
inlierRatio = sum(inlier)/len(inlier)
mat = model.params
print(mat)
print(mat.dot(np.array([0,0,1])))
print(np.array([0,img1.shape[1],1]).dot(mat))
print(np.array([img1.shape[0],0,1]).dot(mat))
print(np.array([img1.shape[0],img1.shape[1],1]).dot(mat))
print(labels)
print(cents)
plt.figure()
plt.imshow(img2, cmap='gray')
plt.scatter([y[1] for y in cents], [x[0] for x in cents], s=150, marker='x', c='red')
In [ ]:
cp = key2[labels==lab]
np.apply_along_axis(np.mean, 0, cp)
In [ ]:
minsamples = 5
d1 = key1[m[:,0],:]
d2 = key2[m[:,1],:]
model, inlier = ransac((d1, d2), AffineTransform, minsamples, 20)
sum(inlier) > minsamples
In [ ]:
cut = 256
minInlier = 16
minInlierRatio = .4
print(pyr[layer_to_use].shape)
for x in np.arange(0, pyr[layer_to_use].shape[1] - cut - 2, cut/2):
for y in np.arange(0, pyr[layer_to_use].shape[0] - cut - 2, cut/2):
img2 = pyr[layer_to_use][y:(y+cut-1), x:(x+cut-1)]
img2 = resize(img2, (sizy, sizx))
try:
orb.detect_and_extract(img2)
desc2 = orb.descriptors
key2 = orb.keypoints
m = match_descriptors(desc1, desc2, cross_check=True)
d1 = key1[m[:,0],:]
d2 = key2[m[:,1],:]
model, inlier = ransac((d1, d2), AffineTransform, minsamples, 10)
nInlier = sum(inlier)
inlierRatio = sum(inlier)/len(inlier)
if (nInlier > minInlier and inlierRatio > minInlierRatio):
plt.figure()
plt.imshow(img2)
print('nInlier={}, ratio={}'.format(nInlier, inlierRatio))
except Exception as e:
print(None)
In [ ]:
from skimage.filters import threshold_adaptive
from skimage.morphology import binary_closing, ball, disk
img1 = pyr[layer_to_use]
thrd = threshold_adaptive(img1, 201)
In [ ]:
from skimage.morphology import binary_closing, binary_opening, ball, disk
from skimage.measure import regionprops, label
thrd2 = np.bitwise_not(thrd)
thrd2 = binary_closing(thrd2, selem=disk(3))
rp = regionprops(label(thrd2))[0]
In [ ]:
In [ ]:
ls = [r.label for r in regionprops(label(thrd2)) if r.area>20000 and
r.area<80000 and r.solidity > .6]
#and r.area/r.filled_area > .9
In [ ]:
res = np.zeros(thrd.shape)
l = label(thrd2)
for li in ls:
res += (l == li)
In [ ]:
from skimage.morphology import remove_small_holes, watershed, binary_erosion
from skimage.feature import peak_local_max
from scipy import ndimage as ndi
plt.rcParams['figure.figsize'] = [12.0,12.0]
r2 = remove_small_holes(res.astype(np.bool), 5000)
r2 = binary_erosion(r2, selem=disk(11))
plt.imshow(label(r2))
for r in regionprops(label(r2)):
print(np.array(r.bbox) * (ds**layer_to_use))
for p in pyr:
print(p.shape)
In [ ]:
from skimage.morphology import remove_small_holes, watershed, binary_erosion, remove_small_objects
from skimage.feature import peak_local_max
from scipy import ndimage as ndi
from skimage.morphology import binary_closing, binary_opening, ball, disk
from skimage.measure import regionprops, label
from skimage.filters import threshold_adaptive, threshold_otsu, threshold_local
from skimage.morphology import binary_closing, ball, disk, binary_opening
import javabridge
import bioformats
from matplotlib.widgets import RectangleSelector, Button
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from skimage.transform import rescale, pyramid_gaussian, resize
from skimage.color import label2rgb
from simple_detection import *
%matplotlib inline
plt.rcParams['figure.figsize'] = [12.0,12.0]
In [ ]:
def bbox_pix2unit(bbox, start, pixsize):
res = np.array(bbox, dtype=float).reshape((2,2)) * np.array(pixsize, dtype=float) + np.array(start, dtype=float)
return res.reshape((4,))
def aspect(bbox):
(ymin, xmin, ymax, xmax) = bbox
exy = ymax - ymin
exx = xmax - xmin
return (exy / exx) if (exx > exy) else (exx / exy)
def detect_wings_simple(img, start, pixsize,
ds=2, layers=2, thresh_window=351,
minarea=20000, maxarea=80000, minsolidity=.6,
minaspect=.3, plot=False, threshold_fun=None):
# downsample
pyr = [p for p in pyramid_gaussian(img, max_layer= layers, downscale = ds)]
img_ds = pyr[layers]
# adaptive threshold
if threshold_fun is None:
thrd = img_ds > threshold_local(img_ds, thresh_window)
else:
thrd = img_ds > threshold_local(img_ds, thresh_window, method='generic', param=threshold_fun)
# clean a bit
thrd = np.bitwise_not(thrd)
thrd = binary_opening(thrd, selem=disk(4))
# filter objs
ls = [r.label for r in regionprops(label(thrd)) if r.area>minarea and
r.area<maxarea and r.solidity>minsolidity and aspect(r.bbox) > minaspect]
# filtered binary
res = np.zeros(thrd.shape)
l = label(thrd)
for li in ls:
res += (l == li)
# more cleaning, plus some erosion to separate touching wings
r2 = remove_small_holes(res.astype(np.bool), 25000)
r2 = binary_erosion(r2, selem=disk(3))
if plot:
image_label_overlay = label2rgb(label(r2), image=img_ds)
fig, ax = plt.subplots(figsize=(12, 12))
ax.imshow(image_label_overlay)
# get bboxes
bboxes = []
for r in regionprops(label(r2)):
if r.area < (minarea * .8 ):
continue
bboxes.append(np.array(r.bbox) * (ds**layers))
if plot:
minr, minc, maxr, maxc = r.bbox
rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr,
fill=False, edgecolor='red', linewidth=2)
ax.add_patch(rect)
# pixels to units
return [bbox_pix2unit(b, start, pixsize) for b in bboxes]
def scale_bbox(bbox, expand_factor = .15):
(ymin, xmin, ymax, xmax) = tuple(bbox)
yrange = ymax - ymin
xrange = xmax - xmin
return (ymin - yrange * expand_factor / 2., xmin - xrange * expand_factor / 2.,
ymax + yrange * expand_factor / 2., xmax + xrange * expand_factor / 2.)
def read_bf(path):
javabridge.start_vm(class_path=bioformats.JARS, run_headless=True)
img = bioformats.load_image(path, rescale=False)
return img
In [ ]:
img = read_bf('/Volumes/CALM_4TB/Overviews/NG_Overview_020.nd2')
In [ ]:
import logging
logging.basicConfig(format='%(asctime)s - %(levelname)s: %(message)s', level=logging.DEBUG)
detect_wings_simple(img, [-20726,53331], [1.61,1.61], [1,-1], plot=True)
In [ ]:
import os
basedir = '/Volumes/CALM_4TB/Overviews/'
files = [os.path.join(basedir, f) for f in os.listdir(basedir)]
for file in files:
img = read_bf(file)
detect_wings_simple(img, [0,0], [0,0], plot=True)
In [ ]:
detect_wings_simple(img, [0,0], [0,0], plot=True, layers=2)
In [ ]:
pyr = [p for p in pyramid_gaussian(img, max_layer= 4, downscale = 2)]
In [ ]:
img_ds = pyr[2]
# adaptive threshold
thrd = img_ds > threshold_local(img_ds, 301)
#thrd = img_ds > threshold_local(img_ds, 51, method='generic', param=threshold_otsu)
# clean a bit
thrd = np.bitwise_not(thrd)
thrd = binary_opening(thrd, selem=disk(7))
# filter objs
ls = [r.label for r in regionprops(label(thrd)) if r.area>20000 and
r.area<60000 and r.solidity>.6 and aspect(r.bbox) > .25]
# filtered binary
res = np.zeros(thrd.shape)
l = label(thrd)
for li in ls:
res += (l == li)
# more cleaning, plus some erosion to separate touching wings
r2 = remove_small_holes(res.astype(np.bool), 5000)
r2 = binary_erosion(r2, selem=disk(11))
In [ ]:
print(ls)
plt.imshow(remove_small_holes(res.astype(np.bool), 50000))