The example will demonstrate where pyparty
fits into the image processing workflow, without altering core algorithms. This example is adapted directly from scikit-image:
http://scikit-image.org/docs/dev/auto_examples/plot_watershed.html
Configure notebook style (see NBCONFIG.ipynb), add imports and paths. The %run magic used below requires IPython 2.0 or higher.
from __future__ import division
%run NBCONFIG.ipynb
Populating the interactive namespace from numpy and matplotlib
Before getting starting, let me refactor the original example a bit. I want to extract the portion of the code related to imports, watershedding and plotting so that I can reuse it.
from __future__ import division
from scipy import ndimage
from skimage.morphology import watershed
from skimage.segmentation import random_walker
from skimage.feature import peak_local_max
from pyparty import showim, splot
def watershed_and_plot(binaryimage, figsize=(8,2.7) ):
""" multiplot of original image, distance map, watershed results."""
dist = ndimage.distance_transform_edt(binaryimage)
local_maxi = peak_local_max(dist,
indices=False,
footprint=np.ones((3, 3)),
labels=binaryimage)
markers = ndimage.label(local_maxi)[0]
labels = watershed(-dist, markers, mask=binaryimage)
UNIQUELABELS = len(np.unique(labels)) - 1 #background is 1
ax0, ax1, ax2 = splot(1,3, figsize=figsize)
#Note I'm using pyparty's showim() plotting utility
showim(binaryimage, ax0, 'gray', interpolation='nearest', title='Overlapping objects')
showim(-dist, ax1, 'jet', interpolation='nearest', title='Distances')
showim(labels, ax2, 'spectral', interpolation='nearest',
title='%s Watershedded objects' % UNIQUELABELS)
for ax in (ax0, ax1, ax2):
ax.axis('off')
plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1,
bottom=0, left=0,right=1)
plt.show();
Generating particles as masks
First, let's generate two overlapping circles as masks using pure numpy arrays
# Generate an initial image with two overlapping circles
x, y = np.indices((80, 80))
x1, y1, x2, y2 = 28, 28, 44, 52
r1, r2 = 16, 20
mask_circle1 = (x - x1)**2 + (y - y1)**2 < r1**2
mask_circle2 = (x - x2)**2 + (y - y2)**2 < r2**2
image = np.logical_or(mask_circle1, mask_circle2)
watershed_and_plot(image);
Regenerate the overlapping circles in pyparty
While we could generate separate circles, pyparty
basic n-circle shapes builtin. (See canvas.available()) For the task at hand, let's use the "dimer" particle and set its orientation and overlap. I'm going to estimate that the orientation is 30 degrees and the overlap is 18%; the slight difference in results is due to this estimation
from pyparty import Canvas
c=Canvas(background='black', rez=(80,80))
c.add('dimer', center=(35,35), radius_1=20, radius_2=16,
overlap=.18, phi=-30.0)
watershed_and_plot(c.pbinary)
In this case, pyparty
doesn't save a whole lot of work. This is because we were only generating one particle, and once masked, our problem was done.
This is a bit of a digression to show some basic features of pyparty particles. Skip to the next section to continue the watershedding
Let's look at some of the functionality that pyparty
adds. First, it is easier to access the particles independently, as well as basic descriptors (area, perimieter etc...). Admittedly, this is not too difficult to do using pure scikit image and a masked array.
print c.perimeter, c.area
c.particles
[ 189.43860018] [ 1962.]
NAME PTYPE CENTER PHI 0 dimer_0 dimer (35, 35) -30.0
What if we needed to change some aspect of the particle, or the image. For example, let's change orientation and radius of one of the particles, as well as the background and particle color:
# # Access individual particles via "plist"
p1 = c.plist[0]
#Image center is half of the canvas resolution
img_center = c.rx / 2, c.ry / 2
ax = c.show('gray')
ax.scatter(*img_center, color='r')
ax.scatter(*p1.center, color='g')
ax.scatter(*p1.centroid, color='b');
Notice that the particle center and centroid are not the asme. The Centroid is the central point of the two radii in the absence of overlap. The image center is marked in red.
Let's set the center of the particle to the center of the image. In addition, let's rotate the particle, and play with the colors.
p1.phi = -90
p1.color = 'orange'
c.background= 'purple'
p1.center = img_center[0], img_center[1]
ax = c.show(title='Rotated and recolored')
ax.scatter(*p1.center, color='red');
The dimer center is not the center of mass; it is the central point of the two radii in the absence of overlap. The center of mass is available as a descriptor, however
(Allow me to deviate from the task at hand to embrace my inner-dork)
c.add('circle', name='lefteye', radius=2, center=(35,20), color='black')
c.add('circle', name='righteye', radius=2, center=(45,20),color='black')
c.add('square', name='hat', center=(40,3), length=5, color='red')
c.add('rectangle', name='hat2', center=(40,7), length=15, width=7, color='black')
c.add('line', center=(57, 21), length=30, width=1, phi=70, color='black')
c.add('line', center=(22, 21), length=30, width=1, phi=110, color='black')
from skimage.data import moon
c.set_bg(moon(), keepres=True, inplace=True)
c.patchshow(title='greetings from the moon');
01-30 01:21:34 WARNING pyparty.utils: background color has been converted (1-channel to 3-channel RGB)
Ok, enough playing with the peanut. Let's see how watershedding varies as a function of partcle size
Let's set 6 x 6 grid and apply 36 dimers
c = Canvas()
c.grid.xdiv = c.grid.ydiv = 10
c.patchshow(gcolor='b', gstyle='--')
c.grid
CartesianGrid (512 X 512) at 0x3d97b90: X --> 10 divisions (51.2 pixels / div) Y --> 10 divisions (51.2 pixels / div)
Notice that there's just over 50 pixels per grid tile. Thus, we will have our dimers be radii of 12.5 with random orientations and overlap varying from 1-100. Let's also shade their color to reflect overlap, with black being no overlap and green being full overlap.
from random import randint
def randphi(): return randint(0,360)
numtiles = len(c.gpairs('centers'))
for idx, (cx, cy) in enumerate(c.gpairs('centers')):
# 0-.99
normval = idx / numtiles
c.add('dimer', r1=12.5, r2=12.5, center=(cx, cy),
phi=randphi(), color=(0, normval ,0), overlap=normval )
c.show(gcolor='red');
c
Canvas at 0x2c9f050: background --> (512 X 512) : default particles --> 100 particles : 1 types xygrid[100] --> (10p X 10p) : (51.2 X 51.2) [pix/tile]
watershed_and_plot(c.pbinary, figsize=(10,5))
We can see that this particular watershedding call led to several spurious particles, up to 3 per dimer.