scipy.ndimage
, and for the EllipseFit which should rarely of real need)For this notebook:
Your comments are welcome.
Acknowledgements would be highly appreciated; for academic citation please use E. Afik, Robust and highly performant ring detection algorithm for 3d particle tracking using 2d microscope imaging. Sci. Rep. 5, 13584; doi: 10.1038/srep13584 (2015).
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import os
## load (and compile if necessary) the ridge_directed_ring_detector
import pyximport; pyximport.install()
import ridge_directed_ring_detector as ring_detector
#ring_detector.__file__
## a function to plot the rings
pi = np.pi
colourmap = plt.jet
txt_colour = 'r'
ring_colour = 'y'
def plot_rings(xyr, minrad=1, ring_colour = 'g', ring_style=[None,'dotted','dashed'][0], lw=1, txt_colour='b', title=None):
'''
given a 3-column array of (x,y) for ring centres & a radius, plots the
rings
'''
#import matplotlib.pyplot as plt
for c in xyr:
if c[-1]>minrad:
circ = plt.Circle(c[:2][::-1],c[-1],ec=ring_colour,fill=False, lw=lw, ls=ring_style)
plt.gca().add_patch(circ)
i,j,r = c
if txt_colour:
plt.text(j,i,'.', fontsize=20, color=txt_colour)
plt.draw()
ridge_hough = ring_detector.RidgeHoughTransform()
default_params = ridge_hough.params.copy()
ridge_hough.params['sigma'] = 0.8
ridge_hough.params['curv_thresh'] = -50
ridge_hough.params['circle_thresh'] = 2 * np.pi * .33
ridge_hough.params['Rmin'] = 4
ridge_hough.params['Rmax'] = 50
ridge_hough.params['vote_thresh'] = 4
ridge_hough.params['dr'] = 3
ridge_hough.params['eccentricity'] = 0
experimental_params = ridge_hough.params.copy()
# print 'default img analysis parameters:'
# for key in sorted(default_params.keys()):
# print '\t %s = %s' % (key, default_params[key])
print '\nimg analysis parameters in the experiment:'
for key in sorted(experimental_params.keys()):
print '\t %s = %s' % (key, experimental_params[key])
img analysis parameters in the experiment: Rmax = 50 Rmin = 4 circle_thresh = 2.07345115137 curv_thresh = -50 dr = 3 eccentricity = 0 ksize = 5 sigma = 0.8 vote_thresh = 4
## Load picture
imgName = '../data/unprocessed_fig_46.png'
img = plt.imread(imgName )
if img.max()<=1.:
img*=255./img.max()
IMG = img.astype(np.float32)
%%time
ridge_hough = ring_detector.RidgeHoughTransform(IMG)
ridge_hough.img_preprocess()
#ridge_hough.debugging_rings_detection()
ridge_hough.rings_detection()
#detected_rings = ridge_hough.output['rings']
rings = ridge_hough.output['rings_subpxl'] # keep for later (sub-pxl mask example)
print 'NB: first run may be several times slower and not indicative.\n'
NB: first run may be several times slower and not indicative. CPU times: user 183 ms, sys: 3.89 ms, total: 187 ms Wall time: 186 ms
fig,ax = plt.subplots(figsize=(15,10))
im = ax.matshow(IMG, cmap=plt.cm.bone_r, vmax=80)
# fig.colorbar(im)
plot_rings(rings, ring_colour='r')
plt.grid(True)
print 'found {} rings'.format(rings.shape[0])
found 58 rings
if True:
import pandas as pd
#pd.set_option('display.max_rows', 10)
rings_df = pd.DataFrame(rings[rings.T[-1].argsort()[::-1]],
columns=['row','col','rad'],
index=range(1,rings.shape[0]+1),
)
rings_df
row | col | rad | |
---|---|---|---|
1 | 408.176910 | 461.251709 | 48.545418 |
2 | 330.256012 | 401.067688 | 48.238937 |
3 | 343.312378 | 213.903900 | 46.409515 |
4 | 256.172424 | 305.647461 | 46.350906 |
5 | 287.775879 | 612.518372 | 43.204762 |
6 | 393.092865 | 838.347351 | 43.097740 |
7 | 466.228119 | 178.406967 | 42.978962 |
8 | 708.943970 | 249.834793 | 42.313602 |
9 | 357.993530 | 760.224548 | 42.171597 |
10 | 341.070374 | 849.519592 | 41.789619 |
11 | 611.496704 | 52.919159 | 40.357243 |
12 | 647.907166 | 167.828720 | 39.501373 |
13 | 237.966568 | 336.721283 | 37.410336 |
14 | 188.510330 | 767.723755 | 36.317535 |
15 | 292.771240 | 270.796783 | 36.015724 |
16 | 451.596619 | 232.172852 | 34.557400 |
17 | 437.062927 | 156.936401 | 34.303852 |
18 | 278.976196 | 465.001923 | 32.312885 |
19 | 224.318817 | 297.148956 | 31.863239 |
20 | 328.458252 | 657.893372 | 31.599180 |
21 | 391.790588 | 850.987610 | 31.531713 |
22 | 547.764832 | 122.137726 | 30.778458 |
23 | 484.692566 | 96.152039 | 30.183065 |
24 | 25.397516 | 860.769287 | 29.092058 |
25 | 405.008453 | 337.241180 | 28.975523 |
26 | 363.025116 | 546.884460 | 27.798862 |
27 | 507.947754 | 94.895638 | 26.934837 |
28 | 209.471985 | 836.348816 | 26.920893 |
29 | 117.612564 | 747.031677 | 26.076532 |
30 | 228.922379 | 386.000885 | 24.882486 |
31 | 392.876770 | 788.330750 | 24.357092 |
32 | 314.935852 | 781.981567 | 23.164736 |
33 | 293.951538 | 818.902039 | 22.348557 |
34 | 375.199554 | 120.219872 | 22.097008 |
35 | 554.976013 | 186.061630 | 22.074947 |
36 | 457.207977 | 639.623535 | 19.026182 |
37 | 311.242310 | 282.900635 | 16.326797 |
38 | 311.807770 | 329.090118 | 16.182253 |
39 | 279.010162 | 261.195801 | 15.344191 |
40 | 647.061646 | 164.963791 | 15.176590 |
41 | 647.785217 | 230.112228 | 14.739820 |
42 | 662.276855 | 256.523102 | 14.501662 |
43 | 663.716797 | 165.119324 | 13.543608 |
44 | 356.987885 | 824.924072 | 13.377859 |
45 | 254.525543 | 408.055634 | 13.233380 |
46 | 334.132080 | 262.668884 | 13.000826 |
47 | 420.898499 | 703.607361 | 12.900146 |
48 | 658.893433 | 178.048523 | 12.288940 |
49 | 666.953125 | 60.880520 | 12.070793 |
50 | 290.352051 | 447.239746 | 11.642728 |
51 | 272.655518 | 186.926315 | 10.405554 |
52 | 110.185440 | 848.866516 | 9.536131 |
53 | 347.125000 | 329.544128 | 9.361567 |
54 | 417.960297 | 782.177917 | 9.085448 |
55 | 455.712128 | 712.429138 | 7.903854 |
56 | 353.116516 | 660.972229 | 7.545264 |
57 | 73.845337 | 899.971313 | 7.310769 |
58 | 143.483170 | 782.976501 | 4.520035 |
58 rows × 3 columns
Least principal eigenvalue of the Hessian takes roughly 4% of the execusion time, therefore
optimisation attempts should focus on the steps following the Least principal curvature calculation
%%timeit
ridge_hough = ring_detector.RidgeHoughTransform(IMG)
ridge_hough.img_preprocess()
# ridge_hough.debugging_rings_detection()
# ridge_hough.rings_detection()
# # detected_rings = ridge_hough.output['rings']
# rings = ridge_hough.output['rings_subpxl'] # keep for later (sub-pxl mask example)
100 loops, best of 3: 7.92 ms per loop
from skimage import data
from scipy import ndimage, misc
image = data.coins()[0:95, 70:370]
%%time
## as the algorithm is currently meant for ridges, convert the image to a gradient magnitude image
IMG = image.astype(np.float32)
sigma = .5
d_mag = ndimage.gaussian_gradient_magnitude(IMG, sigma)
IMG = d_mag
ridge_hough = ring_detector.RidgeHoughTransform(IMG)
ridge_hough.params['sigma'] = 3
ridge_hough.params['Rmin'] = 15
ridge_hough.params['Rmax'] = 30
ridge_hough.params['curv_thresh'] = -35
ridge_hough.params['circle_thresh'] = .4*2*pi
ridge_hough.params['vote_thresh'] = 3
ridge_hough.params['dr'] = 3
ridge_hough.img_preprocess()
#ridge_hough.debugging_rings_detection()
ridge_hough.rings_detection()
#detected_rings = ridge_hough.output['rings']
CHT_rings = ridge_hough.output['rings_subpxl'] # keep for later (sub-pxl mask example)
print 'NB: first run may be several times slower and not indicative.\n'
NB: first run may be several times slower and not indicative. CPU times: user 8.18 ms, sys: 4 µs, total: 8.18 ms Wall time: 7.95 ms
plt.matshow(image, cmap=plt.cm.gray)
plot_rings(CHT_rings, ring_colour='r')
plt.grid(True)
print CHT_rings
[[ 52.33530807 206.2250824 18.84978867] [ 56.56129456 30.30374527 18.23406982] [ 50.78614807 85.13285065 22.39951324] [ 50.84337997 145.36039734 22.25479317] [ 43.76795578 264.86791992 27.9619236 ]]
NB this algorithm does not rely on the user specifying how many circles are in the image
if False:
image = data.coins()[0:95, 70:370].astype(np.float32)
sigma = .5
di,dj = [ndimage.gaussian_filter(image,sigma,order=der) for der in ((0,1),(1,0))]
d_mag = np.sqrt(di**2+dj**2)
plt.matshow(d_mag)
plt.colorbar()
print 'gradient magnitude: \n\tmean = {} \n\tstd = {}'.format(d_mag.mean(), d_mag.std())
gradient magnitude: mean = 6.95160913467 std = 11.1133260727
plt.matshow(np.asarray(ridge_hough.deriv['principal_curv']))
plt.colorbar()
#A = [plot(circ[1],circ[0],'o') for circ in CHT_rings]
plot_rings(CHT_rings)
plt.grid(True)
#print CHT_rings