Generate the actual test data to be used in our papers. These images are very precisely tuned to be representative of our real samples as described in the Considerations/Background
cell.
Configure notebook style (see NBCONFIG.ipynb), add imports and paths. The %run magic used below requires IPython 2.0 or higher!
% run NBCONFIG.ipynb
from pyparty import Grid
from pyparty.utils import grayhist, zoomshow, any2rgb, rgb2uint
Populating the interactive namespace from numpy and matplotlib
There are many variables that were carefully tuned based on measurements.
Gray Background
Particles
Smoothing
Noise
First, let's explicitly define all of the tunable parameters ahead of time. This will ensure reusability if the test images need altered later:
######################
## PARAMETERS #######
#####################
# BACKGROUND
#-------------------
BGRES = (3072, 2304) #Image resolution
BGRANGE = (80,170) #BG Gray range
# PARTICLES
#------------------
DIAM_SINGLE = 30 #Mean pixel diameter single particles
STD_DEV = 3 #Std. dev in diameter norm. distribution
CRANGE = (190,230) #Color range single particles
PERC_SINGLES = .50 #Fraction singles
PERC_DOUBLES = .25 #Frac. dimers
PERC_TRIMERS = .125 #Frac. trimers
PERC_TETRA = .125 #Frac. tetra
MAX_OVERLAP = .2 #Maximum particle % overlap
_GSPACE = (1 + np.sqrt(2)) + .12 #Grid spacing based on max particle size
_CLUSBIAS = 8 #Cluster bias; added N times to average color
_SAMPLES = 10000 #Random samples when generating distribution
# SMOOTHING
#------------------
SIGMA_SMOOTH = 0.1 * DIAM_SINGLE #std dev. gaussian smoothing (pixels)
# NOISE RELATED
NOISE_FRAC = 0.25 #percent coverage
NOISE_STDEV = .25 # % std deviation in noise color
# PLOTTING RELATED
ZOOM_MED = (0,0, 1000,1300) # Mid-scale zoom bounds
ZOOM_NEAR = (200,200, 650,650) # Small-scale zoom bounds
# SAVE PATHS
from os.path import join as opj
SAVE = False # To save or not to save
OUTDIR = '/home/glue/Desktop/testimages/' # Folder to save in
SEM_BG = opj(OUTDIR, 'SEM_test_bg.png') #Just background
SEM_NOSMOOTH = opj(OUTDIR, 'SEM_test_nosmooth.png') #Rasterized image
SEM_SMOOTH = opj(OUTDIR, 'SEM_test_smooth.png') #Smoothed image
SEM_NOISE = opj(OUTDIR, 'SEM_test_noise.png') #Smoothed + noise
SEM_BINARY = opj(OUTDIR, 'SEM_test_binary.png') #Non-smooth; binary
SEM_LABELS = opj(OUTDIR, 'SEM_test_labels.png') #Non-smooth; labels
TEM_BG = opj(OUTDIR, 'TEM_test_bg.png')
TEM_NOSMOOTH = opj(OUTDIR, 'TEM_test_nosmooth.png') #Non-smooth; TEM
TEM_SMOOTH = opj(OUTDIR, 'TEM_test_smooth.png') #Smoothed; TEM
TEM_NOISE = opj(OUTDIR, 'TEM_test_noise.png') #Smoothed + noise; TEM
It is also advantageous to define some plotting utilties that will be used throughout this exercise.
Let's first define a histogram wrapper, with regions of interest labeled by vertical lines, for example the background color upper and lower limits (BGRANGE). This also plots the cumulative distribution function in purple. Also be aware that the histogram limits are automatically adjusted to show the non-zero bins; eg if the lowest pixel has an intenisty of 80, the histograms lower limit will be just below 80.
tripshow() returns two zoomed plots so that we can look at changes in the image at various scales, and the gray histogram defined above.
def labeled_hist(grayimg, *histargs, **histkwargs):
""" Gray histogram; cdf; BGRANGE and CRANGE labeled w/ vlines"""
annotate = histkwargs.pop('annotate', True)
histkwargs.setdefault('cdf', 'purple')
histkwargs.setdefault('ls', ':')
histkwargs.setdefault('lw', 3)
histkwargs.setdefault('xlim', 'auto')
histkwargs.setdefault('alpha', .7)
histkwargs.setdefault('color', 'black')
grayhist(grayimg, *histargs, **histkwargs)
if annotate:
plt.vlines(BGRANGE, 0, 1, colors='orange', linestyles='--')
plt.vlines(CRANGE, 0, 1, colors='blue', linestyles='--')
def tripshow(img, grayimg=None, **histkwargs):
""" In place plot of two zooms and histogram;
title is added to histogram."""
if grayimg is None:
grayimg = img
zoomshow(img, ZOOM_MED, color='pink')
zoomshow(img, ZOOM_NEAR, color='pink')
labeled_hist(grayimg, **histkwargs)
Grid() provides utilties for making a 2d mesh. We merely provide a simple trig function to pattern an uneven contrast along the image. The trig function will yield values ranging from 0-1, so Evelyn came up with a nice scaling that will scale the 2d values up to our desired range (BGRANGE).
power
option to warped().def warped(xx, yy, power=1.0):
""" 2d function to generate uneven contrasts """
return abs(np.cos(xx**power) + np.sin(yy**power))
c = Canvas(rez=BGRES)
g2 = Grid(rez=BGRES)
g2.set_zfcn(warped)
warpedbg = g2.blend_bg(c.image, weight=.9969)
# Scale data from 0-1 to BGCOLOR
k = (BGRANGE[1]-BGRANGE[0])/(warpedbg.max()-warpedbg.min())
warpedbg = (warpedbg-0.75) *k+BGRANGE[0]
graymin, graymax = warpedbg.min(), warpedbg.max()
# Gray image vmin,vmax doesn't work, so plot rgb but show gray limits
ax1, ax2 = splot(1,2)
showim(any2rgb(warpedbg), ax1, title='Showing rgb (%.1f-%.1f)'%(graymin,graymax))
labeled_hist(warpedbg, ax2, title='Gray Histogram');
02-26 19:08:29 WARNING pyparty.tools.grids: 3-channel image convert to 8-bit uint 02-26 19:08:34 WARNING pyparty.utils: color has been converted (1-channel to 3-channel RGB)
First, we add the pre-generated background to a canvas and adjust the grid spacing and preview the result. (patchshow() is called to give a better grid render)
c=Canvas(rez=BGRES, background=warpedbg)
c.grid.xspacing = _GSPACE * DIAM_SINGLE
c.grid.yspacing = _GSPACE * DIAM_SINGLE
CENTERS = c.gpairs('centers')
COUNT = len(CENTERS)
c.patchshow(title='%s available sites' % COUNT, gcolor='r');
02-26 19:08:39 WARNING pyparty.utils: background color has been converted (1-channel to 3-channel RGB)
Next, we define several functions used to generate particles with all of the random distributions and properties espoused earlier. This gets tedious, keeping track of all the random numbers and distributions, so feel free to skip to results.
from random import randint, random
def rint(x): return int(round(x, 0))
# Normal distribution of diameters and colors
RAND_DIAM = np.random.normal(loc=DIAM_SINGLE,
scale=STD_DEV,
size=_SAMPLES+1)
# Scaled beta distribution of brighntess
RAND_COLOR = np.random.beta(2, 2, size=_SAMPLES+1)
RAND_COLOR = CRANGE[0] + RAND_COLOR*(CRANGE[1]-CRANGE[0])
def shuffled_indicies():
""" Shuffled list of indicies 1-4 of same size as COUNTS.
Randomizes order of 1,2,3,4 while locking proportions."""
if PERC_SINGLES + PERC_DOUBLES + PERC_TRIMERS + PERC_TETRA != 1.0:
raise Exception("SINGLES, DOUBLES, TRIMERS... must sum to 1.0")
csize = COUNT+1 # Seems to fix everything
ones = np.empty(rint( csize*PERC_SINGLES ), dtype=int)
twos = np.empty(rint( csize*PERC_DOUBLES ), dtype=int)
threes = np.empty(rint( csize*PERC_TRIMERS ), dtype=int)
fours = np.empty(rint( csize*PERC_TETRA ), dtype=int)
ones.fill(1); twos.fill(2); threes.fill(3); fours.fill(4)
arrays = np.concatenate((ones,twos,threes,fours))
np.random.shuffle(arrays)
# Sometimes under counting
diff = len(arrays) - csize
if diff > 0:
np.append(arrays, arrays[0:diff])
return arrays
def rran(samples=1):
""" Random radius and corresponding color from RAND_DIAM/RAND_COLOR """
radii = [rint(RAND_DIAM[randint(0,_SAMPLES)]/2.0) for i in range(samples)]
colors = [RAND_COLOR[randint(0,_SAMPLES)] for i in range(samples)]
# Average color, bias for clusters
color = np.array(colors).mean()
if color < CRANGE[0]:
color = CRANGE[0]
elif color > CRANGE[1]:
color = CRANGE[1]
if samples > 1:
color = color + samples*_CLUSBIAS
if color > 255:
color = 255
return radii, color
def rphi():
""" Random angle 0-360"""
return rint(random() * 360.0)
def roverlap():
""" Random intra-particle overlap 0-.2 """
return 0 + (random() * MAX_OVERLAP )
Let's look at the diameter and brightness distribution of the single particles. The diameters are distributed normally; however, the brightness is distributed according to a beta distribution (see RELATED) to explore the color range more thoroughly.
ax1, ax2 = splot(1,2, figsize=(8,5))
ax1.hist(RAND_DIAM, bins=25, color='r', alpha=.5);
ax1.set_title("Single Diameter Dist. (Normal)")
ax1.set_xlabel("Diameter (pixels)")
ax2.hist(RAND_COLOR, bins=25, color='b', alpha=.5);
ax2.set_title("Single Color Dist. (Beta)")
ax2.set_xlabel("Pixel Brightness");
With all of the necessary functionality specified above, we now add the particles to the grid centers and call the tripshow() plot to examine it at length scales.
c.clear_particles()
ps = ('circle', 'dimer', 'trimer', 'tetramer')
RAND_IDXS = shuffled_indicies()
for i, center in enumerate(CENTERS):
idx = RAND_IDXS[i]
pname = ps[idx-1]
radii, pcolor = rran(idx)
shared = {'center':center, 'phi':rphi(), 'overlap':roverlap()}
if idx == 1:
c.add(pname, radius=radii[0], color=pcolor, **shared)
else:
c.add(pname, *radii, color=pcolor, **shared)
tripshow(c.image);
02-26 19:08:58 WARNING pyparty.utils: 3-Channel converted to 1-channel (gray).
We nicely see separate bins for our different particle species. Notice that the regions separating the particle classes have large grayscale gaps. Smoothing should remove these gaps by blurring particle boundaries.
from skimage.filter import gaussian_filter
smoothed = gaussian_filter(c.image, SIGMA_SMOOTH, mode='nearest')
tripshow(smoothed)
02-26 19:09:11 WARNING pyparty.utils: 3-Channel converted to 1-channel (gray).
Indeed the smoothing has filled in the color gaps without transforming the background much.
We can see exactly where smoothing acted by taking the absolute difference of the smoothed image from the original. I'll use the jet colormap to really make the edges pop; since this is a grayscale plot, the color is exagerrated!
grayhalos = rgb2uint(np.abs(c.image - smoothed))
zoomshow(grayhalos, ZOOM_NEAR, color='pink', cmap=plt.cm.jet);
Next, we add the guassian noise as explained above.
import pyparty.noise as pnoise
MIDBG = (BGRANGE[0] + BGRANGE[1]) / 2.0
graysmooth = rgb2uint(smoothed)
graynoise = pnoise.normal(graysmooth, NOISE_FRAC,
stddev=NOISE_STDEV, mean=MIDBG )
noisyimg = any2rgb(graynoise)
tripshow(noisyimg)
02-26 19:09:29 WARNING pyparty.utils: color has been converted (1-channel to 3-channel RGB) 02-26 19:09:33 WARNING pyparty.utils: 3-Channel converted to 1-channel (gray).
Just an aside, we can easily invert the image to get TEM Test data. - Literally just doing [(1,1,1) - image] or [255 - image].
from pyparty.utils import invert
tripshow(invert(noisyimg), annotate=False);
02-26 19:09:43 WARNING pyparty.utils: 3-Channel converted to 1-channel (gray).
The histogram is spred now because the noise is widely distributed over the grayscale.
We can now use pyparty's builtin constructs to get a quantitative overview of the particles:
print c.ptype_count
print c
{'tetramer': 146, 'circle': 585, 'dimer': 293, 'trimer': 146} Canvas at 0x7f1b741a2650: background --> (3072 X 2304) : ndarray particles --> 1170 particles : 4 types xygrid[1200] --> (30p X 40p) : (76.8 X 76.8) [pix/tile]
c.show('pbinary', title='%.2f perc. coverage' % (100* c.pixarea) );
We can indirectly check for the particle overlap by comparing the sum of the individual particle areas to the total pixel area in the binary image.
particle_area = c.area.sum()
binary_area = c.pixarea
print 'Overlap error: %.6f percent of total picture' % \
(100.0* ((particle_area / c.pixcount) - binary_area))
Overlap error: 0.000000 percent of total picture
It is very straightforward to label by particle type. Let's use the following color conventions throughout:
Until a canonical solution for storing multiple canvases in employed, the order here WILL AFFECT the order of plotting in the deconstructed histogram
from collections import OrderedDict
COLORS = OrderedDict([ ('noise', 'teal'), ('single','r'), ('dimer','g'), ('trimer','y'), ('tetramer','magenta'),
('background', 'gray')])
clab = Canvas.copy(c)
def label(p):
ptype = p.ptype
if ptype == 'circle':
ptype = 'single'
if ptype not in COLORS:
raise Exception("Invalid ptype %s" % p.ptype)
p.color = COLORS[ptype]
return p
clab.pmap(label, inplace=True)
zoomshow(clab.image, (0,0,500,500), color='white' );
The of_ptypes methods makes it easy to grab groups by name from the canvas. Let's separate the species
singles, dimers = c.of_ptypes('circle'), c.of_ptypes('dimer')
trimers, tetramers = c.of_ptypes('trimer'), c.of_ptypes('tetramer')
GROUPS = [singles, dimers, trimers, tetramers]
ax1, ax2 = splot(1,2, figsize=(8,4))
chartkwds = {'shadow':True,'autopct':'%1.1f%%','labels':COLORS.keys()[1:5],
'colors':COLORS.values()[1:5]}
netcount = [len(x) for x in GROUPS]
netarea = [x.area.sum() for x in GROUPS]
ax1.pie(netcount, **chartkwds)
ax2.pie(netarea, explode=(0,0,0,0.1), **chartkwds)
ax1.set_title("Particle distribution", size=15)
ax2.set_title("Area distribution", size=15);
A very informative illustration is the contribution to the grayscale histogram by the various components; that is, being able to separate the noise, particle contributions, background etc... this is only possible with test data. One must be careful in the decomposition. For example, since noise and particles are being plotted independently, the background must then be plotted sans particles and sans noise. This is done in the next cell with boolean:
Gray background:
Smoothed Singles, Dimers etc.. cast on a black background:
Noise sans gray background:
These set operations are done below and plotted for clarity:
#from pyparty.tools.arraytools import intersect, differ
from pyparty.utils import crop
def differ(img1, img2):
return (img1 != img2) * img1
def intersect(img1, img2):
return (img1 == img2) * img1
# Noise as difference between noisy smoothed and noise-free smooth
just_noise = differ(graynoise, graysmooth)
# Harder is to separate background from noise AND particles
# Step 1: Get background with noise removed
nonoise = intersect(graysmooth, graynoise)
# Step 2: Get just smoothed particles on a gray background
cblank = Canvas.copy(c, background='black')
just_particles = rgb2uint(gaussian_filter(cblank.image, SIGMA_SMOOTH, mode='nearest'))
# Step 3: Mask the smoothed particles to make "holes"
# - Cutoff for masking is mid-point of BG
holes = np.invert(just_particles > MIDBG)
# Step 4: Keep pixels not appearing over particles
just_bg = nonoise * holes
axes = splot(2,2, figsize=(8,8))
SETARRAYS = [('just noise', just_noise), ('sans noise',nonoise),
('holes', holes), ('just bg', just_bg)]
for idx, entry in enumerate(SETARRAYS):
name, img = entry
ax = axes[idx]
showim(crop(any2rgb(img), ZOOM_NEAR), ax, title=name);
02-26 19:10:29 WARNING pyparty.utils: color has been converted (1-channel to 3-channel RGB) 02-26 19:10:29 WARNING pyparty.utils: color has been converted (1-channel to 3-channel RGB) 02-26 19:10:30 WARNING pyparty.utils: color has been converted (1-channel to 3-channel RGB) 02-26 19:10:30 WARNING pyparty.utils: color has been converted (1-channel to 3-channel RGB)
Now, we will store the individual components (noise, background, particles) in a container of the same order as we saw in the Ordered dict above.
ALL = [just_noise]
for cvs in GROUPS:
cvs.background ='black'
out = rgb2uint(gaussian_filter(cvs.image, SIGMA_SMOOTH, mode='nearest'))
out = (out > MIDBG) * out
ALL.append(out)
ALL.append(just_bg)
You may have noticed the following operation was used several times so far:
This is the correction to a subtle artifact. In the decomposition, we project the particles onto a black background and then smooth them. Smoothing against a gray background results in a smaller halo since smoothing can't lower the contrast below the local background. Therefore, we cut the blur off at the middle background value, MIDBG. In the histograms below, you'll see a sharp cutoff of particle signals around this value.
AXES = splot(2,2, figsize=(12,18))
BINS=256
LOWSTART = 30
HIGHSTART = 170
END= 250
STYLE1 = {'histtype':'stepfilled', 'alpha':.9}
STYLE2 = {'histtype':'step', 'lw':2}
for idx, ax in enumerate(AXES):
# Plot full histogram; store its ylim; cdf for top row
if idx < 2:
labeled_hist(noisyimg, ax, annotate=False)
else:
labeled_hist(noisyimg, ax, annotate=False, cdf=False)
ymin, ymax = ax.get_ylim()
if idx % 2:
style = STYLE1
else:
style = STYLE2
ax.hist(ALL, color=COLORS.values(), label=COLORS.keys(),
bins=BINS, **style)
ax.set_ylim(ymin, ymax)
ax.set_xlim(LOWSTART,END)
ax.minorticks_on()
ax.set_title('Histogram Smoothed Components')
l = ax.legend(loc=2, prop={'size':12})
l.set_zorder(20) # Ensurse legend on top
# Halve y-axis; change xaxis
for ax in AXES[2:4]:
ax.set_ylim(ymin, ymax/4) #Just a guess at optimal value
ax.set_xlim(HIGHSTART, END);
02-26 19:19:36 WARNING pyparty.utils: 3-Channel converted to 1-channel (gray). 02-26 19:19:41 WARNING pyparty.utils: 3-Channel converted to 1-channel (gray). 02-26 19:19:45 WARNING pyparty.utils: 3-Channel converted to 1-channel (gray). 02-26 19:19:50 WARNING pyparty.utils: 3-Channel converted to 1-channel (gray).
We see that the effects of smoothing can be seen throughout the image histogram, even at low intensities!
from skimage.io import imsave
clab.background='black' #for saving
if SAVE:
imsave(SEM_BG, c.background)
imsave(SEM_NOSMOOTH, c.image)
imsave(SEM_SMOOTH, smoothed)
imsave(SEM_NOISE, noisyimg)
imsave(SEM_BINARY, c.pbinary)
imsave(SEM_LABELS, clab.image)
imsave(TEM_BG, invert(c.background))
imsave(TEM_NOSMOOTH, invert(c.image))
imsave(TEM_SMOOTH, invert(smoothed))
imsave(TEM_NOISE, invert(noisyimg))
print 'SAVING COMPLETE'
else:
print 'SKIPPING SAVE STEP'
SAVING COMPLETE