In [ ]:
from skimage import io, exposure, transform, img_as_ubyte
In [346]:
files = !ls unstained/*tif
for f in files:
img = io.imread(f)
rimg = filters.gaussian_filter(img, 3)
rimg = transform.resize(rimg, (img.shape[0]//4, img.shape[1]//4))
rimg = exposure.equalize_adapthist(rimg, clip_limit=0.05)
rimg = rimg / rimg.max()
rimg = img_as_ubyte(rimg)
f = f.replace('tif', 'png')
io.imsave(f, img)
io.imsave(f.replace('/u', '/view_u'), rimg)
In [347]:
!rm unstained/*tif
In [342]:
# for rescaling view images
files = !ls unstained/u*png
for f in files:
#f = files[0]
img = io.imread(f)
img = filters.gaussian_filter(img, 3)
rimg = transform.resize(img, (img.shape[0]//4, img.shape[1]//4))
#rimg = filters.gaussian_filter(rimg, 0.8)
rimg = exposure.equalize_adapthist(rimg, clip_limit=0.05)
rimg = rimg / rimg.max()
rimg = img_as_ubyte(rimg)
#imshow(rimg);
f = f.replace('/u', '/view_u')
io.imsave(f, rimg)
In [264]:
!rm unstained/*tif
Function to compute angle with FT line fit
In [265]:
%pylab inline
from skimage.exposure import cumulative_distribution
from scipy.stats import linregress
from skimage.filters import threshold_otsu
from matplotlib import ticker
plt.rcParams['image.cmap'] = 'gray_r'
plt.rcParams['image.interpolation'] = 'none'
In [294]:
def angle_ft_line_fit(img, debug=False):
"""
Calculate preferred orientation in image with a line fit in FT.
Parameters
----------
debug : bool
Wheter to show a graph of image + FT line fit.
Returns
-------
float
Angle
"""
# FT power spectrum
F = np.abs(fftshift(fft2(img)))
# do not calculate log(0)
F = np.log(F.clip(0))
# threshold
number_of_pixels = img.shape[0]//2
threshold = 1 - number_of_pixels/img.size
# keep top value pixels
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
b = [0, 1]
r = 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)
else:
m,c,r,pv,err = linregress(x,y)
b = (m,c)
# calculate angle (assume line goes through center)
# as our coordinate system have y positive down - use negative m
angle = int(arctan(-b[0]) / pi * 180 + 90) % 180
# show image, FT and fit
if debug:
f, ax = subplots(ncols=2,figsize=(8,8))
# add calculated line
# polynomial generator
p = poly1d(b)
height, width = img.shape
if angle != 90:
line = ([0, width], [p(0), p(width)])
ang0_rad = -angle/180*pi
dy0 = tan(ang0_rad)
line0 = ([0, width], [height//2-width//2*dy0, height//2+width//2*dy0])
else:
line = ([width//2, width//2], [0,height])
line0 = ([0, width], [height//2, height//2])
im0 = ax[0].imshow(img)
im1 = ax[1].imshow(F)
cb0 = plt.colorbar(im0, ax=ax[0], fraction=0.045)
cb1 = plt.colorbar(im1, ax=ax[1], fraction=0.045)
tick_locator = ticker.MaxNLocator(nbins=5)
cb0.locator = tick_locator
cb0.update_ticks()
#binary image
tick_locator = ticker.MaxNLocator(nbins=1)
cb1.locator = tick_locator
cb1.update_ticks()
ax[1].plot(line[0], line[1], 'g', linewidth=2)
ax[0].plot(line0[0], line0[1], 'b', linewidth=2)
ax[1].plot(line0[0], line0[1], 'b', linewidth=2)
ax[0].set_title('mean:{:0.0f}'.format(img.mean()))
ax[1].set_title('angle:{:0.0f}, r$^2$:{:0.2}'.format(angle,r**2))
ax[0].set_xlim(0,width)
ax[0].set_ylim(height,0)
ax[1].set_xlim(30,70)
ax[1].set_ylim(70,30)
ax[0].get_xaxis().set_ticks([])
ax[0].get_yaxis().set_ticks([])
ax[1].get_xaxis().set_ticks([])
ax[1].get_yaxis().set_ticks([])
plt.show()
return angle, r**2
Create three line fit of the same image with different threshold, for comparison in report
In [271]:
from skimage.io import imread
from itertools import product
img = imread('unstained/u13v0ch0z0.png') # image which has strong response
# blocks
bs = 100 # approx block size
iy, ix = img.shape
by, bx = iy//bs, ix//bs # number of blocks
bsy, bsx = iy//by, ix//bx # block size, spread spare pixels
# per block
n = 0
for j,i in product(range(by), range(bx)):
x,y = j*bsx, i*bsy # pos
img_block = img[y:y+bsy, x:x+bsx]
mean = img_block.mean()
# emtpy image
if mean < 2:
continue
# pick some blocks
if n >= 25:
break
n+=1
angle_ft_line_fit(img_block, debug=True)
In [273]:
from skimage.io import imread
from skimage import exposure
from itertools import product
from sklearn.externals import joblib
from sklearn.datasets.base import Bunch
In [282]:
def histogram_ft_line(filename):
"Draw a histogram of angles based on angle_ft_line_fit()"
img = imread(filename) # image which has strong response
# blocks
bs = 100 # approx block size
iy, ix = img.shape
by, bx = iy//bs, ix//bs # number of blocks
bsy, bsx = iy//by, ix//bx # block size, spread spare pixels
# per block
h = []
for i,j in product(range(by), range(bx)):
y,x = i*bsy, j*bsx # pos
img_block = img[y:y+bsy, x:x+bsx]
mean = img_block.mean()
# emtpy image
if mean < 0.1:
continue
b = Bunch()
b.angle, b.r = angle_ft_line_fit(img_block)
h.append(b)
return h
In [348]:
files = !ls unstained/u*0.png
histograms_ft_line = []
for f in files:
b = Bunch()
b.f = f
b.h = histogram_ft_line(f)
histograms_ft_line.append(b)
joblib.dump(histograms_ft_line, 'unstained/hist_ft_line.pkl')
Out[348]:
In [349]:
for hist in histograms_ft_line:
img_name = hist.f.replace('/u', '/view_u')
img = io.imread(img_name)
h = np.zeros(180)
for angle in hist.h:
h[angle.angle] += 1
figure(figsize=(11,4))
subplot(121)
imshow(img)
yticks([])
xticks([])
title('(a) Image')
colorbar()
subplot(122)
plot(h)
title('(b) Histogram of angles')
plt.show()
In [350]:
for hist in histograms_ft_line:
img_name = hist.f.replace('/u', '/view_u')
img = io.imread(img_name)
h = np.zeros(180)
for angle in hist.h:
if angle.r > 0.3:
h[angle.angle] += 1
figure(figsize=(11,4))
subplot(121)
imshow(img)
yticks([])
xticks([])
title('(a) Image')
colorbar()
subplot(122)
plot(h)
title('(b) Histogram of angles')
plt.show()
In [276]:
import pyfftw
fft2 = pyfftw.interfaces.numpy_fft.fft2
from os.path import basename
In [351]:
files = !ls unstained/u*png
#files = ['unstained/u3v1ch0z0.png']
for f in files:
i = io.imread(f)
F = fft2(i)
F = np.abs(F)
F = fftshift(F)
F = np.log(F)
F = F/F.max()
io.imsave(f.replace('d/u','d/fourier_u'), F)
In [279]:
from skimage import io, filters
from scipy.optimize import curve_fit
import pyfftw
fft2 = pyfftw.interfaces.numpy_fft.fft2
def ft_radial_1d_signature(filename):
"""
Compute radial signature of FT where values are above a threshold.
"""
F = io.imread(filename)
# threshold
F[F<F.max()*0.65] = 0
# positions
x,y = np.where(F!=0)
# center of mass
# should be center, but its easy to calculate
cm = y.mean(), x.mean()
# make (0,0) to cm
y = y - cm[0]
x = x - cm[1]
angles = np.zeros(180)
n = 0
for point in zip(y,x):
# only upper half, as its symetric
if point[0] < 0:
continue
# special case in range 0-180
if point[0] == 0 and point[1] < 0:
continue
# center of mass
if point == (0,0):
continue
angle = floor(arctan2(*point) / pi * 180)
#angles[angle] += F[point[0]-cm[0], point[1]-cm[1]]
angles[angle] += 1
return angles
In [352]:
files = !ls unstained/fourier*png
histograms_ft_radial = []
for f in files:
b = Bunch()
b.f = f
b.h = ft_radial_1d_signature(f)
histograms_ft_radial.append(b)
joblib.dump(histograms_ft_radial, 'unstained/hist_ft_radial.pkl')
Out[352]:
In [354]:
for hist in histograms_ft_radial:
img_name = hist.f.replace('fourier_', 'view_')
img = io.imread(img_name)
figure(figsize=(11,4))
subplot(121)
imshow(img)
yticks([])
xticks([])
title('(a) Image')
colorbar()
subplot(122)
plot(hist.h)
title('(b) Histogram of angles')
plt.show()
In [359]:
def hog(filename):
img = io.imread(filename).astype(int)
dx = np.zeros_like(img)
dy = np.zeros_like(img)
angles = np.zeros_like(img)
dy[1:-1,:] = img[2:,:] - img[0:-2,:]
dx[:,1:-1] = img[:,2:] - img[:,0:-2]
angles = floor(arctan2(dy,dx)/pi*180 + 90) % 180
angles = angles.astype(int)
dy = dy[1:-1, 1:-1]
dx = dx[1:-1, 1:-1]
mag = sqrt(dx**2 + dy**2)
mag = mag/mag.max()/2 # normalize
angles = angles[1:-1, 1:-1]
h = np.zeros(180, dtype=float)
for i in range(180):
mask = angles == i
h[i] += mag[mask].sum()
return h
In [360]:
files = !ls unstained/u*png
histograms_gradient = []
for f in files:
b = Bunch()
b.f = f
b.h = hog(f)
histograms_gradient.append(b)
joblib.dump(histograms_gradient, 'unstained/hist_gradient.pkl')
Out[360]:
In [361]:
for hist in histograms_gradient:
img_name = hist.f.replace('/u', '/view_u')
img = io.imread(img_name)
figure(figsize=(11,4))
subplot(121)
imshow(img)
yticks([])
xticks([])
title('(a) Image')
colorbar()
subplot(122)
plot(hist.h)
title('(b) Histogram of angles')
plt.show()
In [69]:
import csv
tumor_cells = ([],[])
with open('tumor_cells.csv', 'r') as csvfile:
rows = csv.reader(csvfile, delimiter=',', quotechar='|')
for n, row in enumerate(rows):
if n < 3: continue
if len(row):
tumor_cells[0].append(int(row[0]))
tumor_cells[1].append(int(row[1]))
In [68]:
tumor_stroma = ([],[])
with open('tumor_stroma.csv', 'r') as csvfile:
rows = csv.reader(csvfile, delimiter=',', quotechar='|')
for n, row in enumerate(rows):
if n < 3: continue
if len(row):
tumor_stroma[0].append(int(row[0]))
tumor_stroma[1].append(int(row[1]))
In [108]:
tumor_env = ([],[])
with open('tumor_environment.csv', 'r') as csvfile:
rows = csv.reader(csvfile, delimiter=',', quotechar='|')
for n, row in enumerate(rows):
if n < 3: continue
if len(row):
tumor_env[0].append(int(row[0]))
tumor_env[1].append(int(row[1]))
In [112]:
tumor_collagen = ([],[])
with open('tumor_collagen.csv', 'r') as csvfile:
rows = csv.reader(csvfile, delimiter=',', quotechar='|')
for n, row in enumerate(rows):
if n < 3: continue
if len(row):
tumor_collagen[0].append(int(row[0]))
tumor_collagen[1].append(int(row[1]))
In [120]:
breast_cancer = ([],[])
with open('breast_cancer.csv', 'r') as csvfile:
rows = csv.reader(csvfile, delimiter=',', quotechar='|')
for n, row in enumerate(rows):
if n < 3: continue
if len(row):
breast_cancer[0].append(int(row[0]))
breast_cancer[1].append(int(row[1]))
In [122]:
%pylab inline
matplotlib.rcParams.update({'font.size': 16})
figure(figsize=(12,8))
#plt.plot(*breast_cancer, label='breast cancer', linestyle='', marker='o')
plt.plot(*tumor_cells, label='tumor cells', linestyle='', marker='o')
plt.plot(*tumor_env, label='tumor environment', linestyle='', marker='o')
plt.plot(*tumor_collagen, label='tumor collagen', linestyle='', marker='o')
plt.plot(*tumor_stroma, label='tumor stroma', linestyle='', marker='o')
xlim(1950,2014)
yscale('log')
xlabel('Year')
ylabel('Number of published articles in pubmed')
legend(loc='upper left', title='Search terms', fontsize=14, numpoints=1)
Out[122]:
In [10]:
%pylab inline
In [11]:
from skimage.io import imread
from skimage.exposure import equalize_hist, histogram, adjust_log, adjust_gamma, equalize_adapthist
from skimage import img_as_float, img_as_ubyte, img_as_int
In [12]:
matplotlib.rcParams.update({'font.size': 16})
img = imread('mp.png').astype(np.uint16)
eimg = equalize_adapthist(img, clip_limit=0.025)
eimg = img_as_ubyte(eimg)
h = histogram(img)
he = histogram(eimg)
In [18]:
fig,axs = subplots(nrows=2, ncols=2, figsize=(18,12))
fig.subplots_adjust(wspace=10, hspace=5)
fig.tight_layout(h_pad=3)
im0 = axs[0,0].imshow(img, cmap='gray_r')
im1 = axs[1,0].imshow(eimg, cmap='gray_r')
plt.colorbar(im0, ax=axs[0,0])
plt.colorbar(im1, ax=axs[1,0])
axs[0,1].plot(h[1], h[0], '.')
axs[1,1].plot(he[1], he[0], '.')
axs[0,0].set_title('(a) Original image', y=1.08)
axs[1,0].set_title('(c) Equalized image', y=1.08)
axs[0,1].set_title('(b) Histogram of levels original')
axs[1,1].set_title('(d) Histogram of levels equalized')
axs[0,1].set_ylabel('Number of pixels')
axs[1,1].set_ylabel('Number of pixels')
axs[0,1].set_xlabel('Intensity')
axs[1,1].set_xlabel('Intensity')
axs[0,1].ticklabel_format(style='sci', axis='y', scilimits=(0,0))
axs[1,1].ticklabel_format(style='sci', axis='y', scilimits=(0,0))
axs[0,1].set_yscale('log')
axs[1,1].set_yscale('log')
axs[0,0].xaxis.tick_top()
axs[1,0].xaxis.tick_top()
axs[0,1].set_xlim(0,255)
axs[1,1].set_xlim(0,255);
In [16]:
for i in range(10):
s = (img == i).sum()
print(s)
In [165]:
img[img==0] = 3
img = img - 3
In [185]:
for i in range(10):
print( (img == i).sum() )
In [189]:
eeimg = equalize_hist(img)#, mask=(img > 0))
eeimg = img_as_ubyte(eeimg)
h = histogram(img)
he = histogram(eeimg)
In [190]:
fig,axs = subplots(nrows=2, ncols=2, figsize=(18,12))
fig.subplots_adjust(wspace=10, hspace=5)
fig.tight_layout(h_pad=3)
im0 = axs[0,0].imshow(img, cmap='gray_r')
im1 = axs[1,0].imshow(eeimg, cmap='gray_r')
plt.colorbar(im0, ax=axs[0,0])
plt.colorbar(im1, ax=axs[1,0])
axs[0,1].plot(h[1], h[0], '.')
axs[1,1].plot(he[1], he[0], '.')
axs[0,0].set_title('(a) Original image')
axs[1,0].set_title('(c) Equalized image')
axs[0,1].set_title('(b) Histogram of levels original')
axs[1,1].set_title('(d) Histogram of levels equalized')
axs[0,1].set_ylabel('Number of pixels')
axs[1,1].set_ylabel('Number of pixels')
axs[0,1].set_xlabel('Intensity')
axs[1,1].set_xlabel('Intensity')
axs[0,1].ticklabel_format(style='sci', axis='y', scilimits=(0,0))
axs[1,1].ticklabel_format(style='sci', axis='y', scilimits=(0,0))
axs[0,1].set_yscale('log')
axs[1,1].set_yscale('log')
axs[0,1].set_xlim(0,255)
axs[1,1].set_xlim(0,255);
In [17]:
%pylab inline
#from skimage.io import imread
import matplotlib.gridspec as gridspec
plt.rcParams['image.interpolation'] = 'none'
plt.rcParams['image.cmap'] = 'gray'
figsize(4,4)
In [18]:
size = 256
img = np.zeros((size,size), dtype=np.uint8)
t = linspace(start=0, stop=50*pi, endpoint=False, num=size)
x,y = meshgrid(t, t)
img[:,:] = 127 + 127*sin(x)
imshow(img);
In [20]:
F = fft2(img)
# scale image for viewing - do not take log of zero
F_pow = np.abs(F)
F_pow = log10(F_pow.clip(1))
fig, axs = subplots(ncols=2, figsize=(14,5))
plt.setp(axs, xticks=[], yticks=[])
im0 = axs[0].imshow(img, cmap='gray_r')
axs[0].set_title('(a) Image of horizontal sinus wave')
colorbar(im0, ax=axs[0])
im1 = axs[1].imshow(fftshift(F_pow), cmap='gray_r')
axs[1].set_title('(b) Logarithmic Fourier spectrum')
colorbar(im1);
In [10]:
F_pow.max()
e**F_pow.max()
Out[10]:
In [35]:
%pylab inline
from skimage.io import imread
img = imread('two_spots.png')
imshow(img, cmap='gray');
In [36]:
img.shape
Out[36]:
In [37]:
img = img[:,:,0]
In [38]:
img = 255-img # invert
In [39]:
imshow(img, cmap='gray');
In [40]:
from skimage import exposure
In [42]:
eimg = exposure.equalize_adapthist(img)
imshow(eimg);
In [43]:
from skimage import io
In [44]:
io.imsave('spots.png', eimg)
In [398]:
img = io.imread('/Users/arve/Dokumenter/TFY4500/img report/tilescan.png')
In [399]:
img.shape
Out[399]:
In [400]:
timg = transform.resize(img, (img.shape[0]//4, img.shape[1]//4))
for i in range(4):
figure()
imshow(timg[:,:,i])
colorbar()
plt.show()
In [403]:
io.imsave('/Users/arve/Dokumenter/TFY4500/img report/tilescan.png', 1-timg[:,:,0])
In [404]:
img = io.imread('/Users/arve/Dokumenter/git/H14/TFY4500/unstained/view_u0v2ch1z0.png')
io.imsave('/Users/arve/Dokumenter/git/H14/TFY4500/unstained/view_u0v2ch1z0_r.png', 255-img)
In [435]:
img = io.imread('/Users/arve/Dokumenter/TFY4500/unstained/experiment--row 3-5/u10v0ch1z0.tif')
img = transform.resize(img, (img.shape[0]//4, img.shape[1]//4))
img = filters.gaussian_filter(img, 1)
rimg = exposure.equalize_adapthist(img, clip_limit=0.1)
rimg = img_as_ubyte(rimg)
In [436]:
rimg = rimg[:600, :800]
In [437]:
imshow(rimg)
colorbar();
Out[437]:
In [439]:
io.imsave('/Users/arve/Dokumenter/TFY4500/img report/hyd.png', 255-rimg)