Compute angle histogram in parallel on dataset using fourier transfor line fitting. Differ the threshold for the FT, and save the images so that we can compare.
In [2]:
from IPython import parallel
rc = parallel.Client()
In [3]:
len(rc)
Out[3]:
IPython magic, run on all instances
In [4]:
%px print 'hello'
Block magic
In [5]:
%%px
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
plt.ioff()
print(plt.isinteractive())
Load function which computes angle from image to workers
In [6]:
%%px
def angle_ft_line_fit(img, threshold=0.999, debug=False):
"""
Calculate preferred orientation in image with a line fit in FT.
Parameters
----------
threshold : float
Percentage of pixels to include in FT for calculating
threshold. 0.999 * 512**2 = 262 pixels
Returns
-------
float
Angle
"""
from skimage.exposure import cumulative_distribution
from scipy.stats import linregress
from numpy import abs, log, arctan, where, poly1d, pi
from numpy.fft import fftshift, fft2
from matplotlib import pyplot
plt = pyplot
# FT power spectrum
F = abs(fftshift(fft2(img)))
# do not calculate log(0)
F[F!=0], F[F==0] = log(F[F!=0]), log(F[F!=0].min())
# threshold
cdf = cumulative_distribution(F)
limit = where(cdf[0] > threshold)[0].min()
threshold_value = cdf[1][limit]
F = F > threshold_value
# points
y,x = where(F)
# cases
dx = abs(x.max()-x.min())
dy = abs(y.max()-y.min())
if dx==0:
# we have a vertical line
angle = 90
b = [0, 1]
# solve y=mx+c by least-squares regression
elif dx < dy:
# linregress is imprecise for dx < dy => swap x,y
m,c,r,pv,err = linregress(y,x)
b = (1/m, -c/m)
# calculate angle (assume line goes through center)
angle = (90 - arctan(b[0]) / pi * 180) % 180
else:
m,c,r,pv,err = linregress(y,x)
b = (m,c)
angle = (90 - arctan(b[0]) / pi * 180) % 180
# show image, FT and fit
if debug:
f, ax = plt.subplots(ncols=2, figsize=(8,4))
ax[0].imshow(img)
ax[1].imshow(F)
# add calculated line
# polynomial generator
p = poly1d(b)
height, width = img.shape
if angle != 90:
line = ([0, width], [p(0), p(width)])
else:
line = ([width//2, width//2], [0,height])
ax[1].plot(*line)
ax[1].set_title('ang: {:3.0f} r:{:0.2} err:{:0.2}'
.format(angle,r,err))
ax[1].set_xlim(0,width)
ax[1].set_ylim(height,0)
return angle
print('angle_ft_line_fit defined')
Function to compute histogram of angles (each angle correspond to a slice of the image). This function should be mapped in parallel.
In [26]:
def angle_histogram(arg):
# work around for me not knowing how to dview.map multiple arguments
threshold, filename = arg
from itertools import product
from skimage.filter import threshold_otsu
from skimage.io import imread
from numpy import zeros
img = imread(filename)
bs = 100 # approx block size
iy, ix = img.shape
by, bx = iy//bs, ix//bs # blocks
bsy, bsx = iy//by, ix//bx # block size, spread spare pixels
h = zeros(180) # histogram
for j,i in product(range(by), range(bx)):
x,y = j*bsx, i*bsy # pos
temp_img = img[y:y+bsy, x:x+bsx]
mean = temp_img.mean()
# small image
if temp_img.shape[0] < 50 or temp_img.shape[1] < 50:
continue
# emtpy image
if mean == 0:
continue
# threshold below noise-threshold
ot = threshold_otsu(temp_img)
if (ot < 1):
continue
angle = angle_ft_line_fit(temp_img, threshold=threshold)
angle = int(angle)
h[angle] += 1
return h, threshold, filename
Create a list with all combinations of image and threshold
In [39]:
from numpy import linspace
from itertools import product
files = !ls /notebooks/TFY4500/unstained/u*.png
thresholds = linspace(0.9, 1, num=11, endpoint=False)
combinations = list(product(thresholds, files))
combinations[0]
Out[39]:
Do one of the images in parallel (for testing)
In [40]:
dview = rc.direct_view()
# run in parallel
res = dview.apply(angle_histogram, combinations[0])
In [41]:
from time import sleep
while res.ready() == False:
sleep(1)
res.elapsed
Out[41]:
In [42]:
res[0] # one of the results (what angle_histogram returns)
Out[42]:
Do them all in parallel
In [43]:
res = dview.map(angle_histogram, combinations)
Wait for it to complete
In [47]:
from __future__ import print_function
from time import sleep
while res.ready() == False:
print('.', end='')
sleep(10)
res.elapsed
Out[47]:
In [19]:
# to stop it, interupt the above and run (or use *Clusters* tab)
#res.abort()
Check if we got all computations back
In [50]:
len(res.result) == len(combinations)
Out[50]:
Save results in pngs. As my docker-notebook instance don't cope well with matplotlib in parallel (eats up memory), do this sequentially so that we can find out what uses the memory.
In [65]:
# save results
import pickle
with open('unstained/histograms.pickle', mode='wb') as f:
pickle.dump(res.result, f)
In [77]:
%%writefile plot_memleak.py
import pickle
import matplotlib
import matplotlib.pyplot as plt
# valid backends are:
# pgf, cairo, macosx, cocoaagg, gdk, ps, gtkagg, nbagg, gtk, qt5agg,
# template, emf, gtk3cairo, gtk3agg, wx, qt4agg, tkagg, agg, svg,
# gtkcairo, wxagg, webagg, pdf
matplotlib.use('ps')
from skimage.io import imread
with open('unstained/histograms.pickle', 'r') as f:
results = pickle.load(f)
@profile
def saveplot(histogram, threshold, filename):
"Makes png of histogram alongside the image."
img = imread(filename)
# make plot
fig, axs = plt.subplots(ncols=2, figsize=(16,8))
axs[0].imshow(img)
axs[1].plot(histogram)
# save it
fn = filename.replace('ed/u', 'ed/angles-ft-line-fit-u')
fn = fn.replace('.png', '-{:1.3f}'.format(threshold) + '.png')
fig.savefig(fn)
# try to clear up memory
plt.close()
del img
import gc
gc.collect()
print('saved to {}'.format(fn))
# do a small loop
for i in range(5):
saveplot(*results[i])
In [60]:
# interrupt -> garbage collect
import gc
gc.collect()
Out[60]:
In [104]:
import matplotlib
import matplotlib.pyplot as plt
# valid backends are:
# pgf, cairo, macosx, cocoaagg, gdk, ps, gtkagg, nbagg, gtk, qt5agg,
# template, emf, gtk3cairo, gtk3agg, wx, qt4agg, tkagg, agg, svg,
# gtkcairo, wxagg, webagg, pdf
matplotlib.use('ps')
from skimage.io import imread
def saveplot(histogram, threshold, filename):
"Makes png of histogram alongside the image."
img = imread(filename)
# make plot
fig, axs = plt.subplots(ncols=2, figsize=(16,8))
axs[0].imshow(img)
axs[1].plot(histogram)
# save it
fn = filename.replace('ed/u', 'ed/angles-ft-line-fit-u')
fn = fn.replace('.png', '-{:1.3f}'.format(threshold) + '.png')
fig.savefig(fn)
# try to clear up memory
plt.clf()
plt.close()
#del img
#del fig
#print('gc: {}'.format(gc.collect()))
print('saved to {}'.format(fn))
In [105]:
# do a small loop
import gc
for i in range(5):
saveplot(*res[i])
print(gc.collect())
Garbage collect after running 'plt.clf()' and 'plt.close()' inside 'savefig()' seems to do the trick.
In [106]:
for r in res:
saveplot(*r)
gc.collect()
Lets include some of them here, to show off
In [ ]:
from IPython.display import Image
In [113]:
Image('unstained/angles-ft-line-fit-u13v1ch1z0-0.991.png')
Out[113]:
In [114]:
Image('unstained/angles-ft-line-fit-u13v0ch1z0-0.991.png')
Out[114]:
In [115]:
Image('unstained/angles-ft-line-fit-u0v0ch1z0-0.900.png')
Out[115]:
In [116]:
Image('unstained/angles-ft-line-fit-u0v1ch1z0-0.900.png')
Out[116]:
In [117]:
Image('unstained/angles-ft-line-fit-u11v0ch1z0-0.991.png')
Out[117]:
So having a small number of pixels to fit to in fourier spectrum seems to work well