segmentation examples

binary_segmentation_example.py

#!/usr/bin/env python

import numpy
from scipy import ndimage
from image_processing.segmentation.segmentation import SegmentVisitor

"""
This example looks at segmentation from the results of a binary threshold.
The pixel quality code in the ga-neo-landsat-processor repo is a production
example of this methodolgy.

The terminolgy, labels, segments, objects, blobs, used by various people and
packages are basically one and the same.  So if you're familiar with any of
those, then this will be nothing new.
"""

# 1st create a 2D array of random data and set the first pixel to a NaN
img = numpy.random.randint(0, 256, (10, 10)).astype('float32')
img[0, 0] = numpy.nan
print "Random image"
print img

# Simiulate a binary segmented array such as a result from a threshold
b = numpy.zeros(img.shape, dtype='bool')
b[0:2, 0:2] = True
b[-2:, -2:] = True
b[0:2, -2:] = True
b[-2:, 0:2] = True
b[4:6, 4:6] = True
print "Binary segments"
print b

# Now uniquely label each segment/blob/object with an ID
# Some circles also call this an enumerated array
labelled_array, n_labels = ndimage.label(b)
print "Segmented/Labelled array"
print labelled_array

# Define an instance of a SegmentVisitor which we'll use to quickly traverse
# the segments
segments = SegmentVisitor(labelled_array)

# An example of calculating the mean; Note the appearance of a NaN
xbar = segments.segment_mean(img)
print "Segment means not accounting for NaN's"
print xbar

# Now account for the presence of any NaN values
xbar_nan = segments.segment_mean(img, nan=True)
print "Segment means accounting for NaN's"
print xbar_nan
print "Is the mean value for segment 1 finite?"
print numpy.isfinite(xbar[1])
print "Is the mean value for segment 1 finite?"
print numpy.isfinite(xbar_nan[1])

# Get the minimum bounding box of each segment
# With no args this will return the bbox for all segments
bbox = segments.segment_bounding_box()
print "Bounding box of each segment"
print bbox

# 1st segment; Note the shape will be 2D
yse, xse = bbox[1]
subset = img[yse[0]:yse[1], xse[0]:xse[1]]
print "2D subset for segment 1"
print subset
print "2D subset shape for segment 1"
print subset.shape

# 3D example; contains 3 bands; Note the shape will be 3D
img2 = numpy.random.randint(0, 256, (3, 10, 10))
subset2 = img2[:, yse[0]:yse[1], xse[0]:xse[1]]
print "3D subset for segment 1"
print subset2
print "3D subset shape for segment 1"
print subset2.shape

# 3rd segment
yse, xse = bbox[3]
subset = img[yse[0]:yse[1], xse[0]:xse[1]]
print "2D subset for segment 3"
print subset
subset2 = img2[:, yse[0]:yse[1], xse[0]:xse[1]]
print "3D subset for segment 3"
print subset2

"""
I've only included some basic functions, but essentially using the
get_segment_locations and get_segment_data one can do and calculate almost
anything.
I have code elsewhere to calculate compactness, centroid, rectangularity and
roundness of each segment/object/blob, but I'll expose them through the
SegmentVisitor soon enough.

Check the docstrings for get_segment_data, get_segment_locations,
segment_total, segment_mean, segment_max, segment_min, segment_stddev,
segment_area, reset_segment_ids, sieve_segments and segment_bounding_box
for more information.
"""

rasterisation_segmentation_example.py

#!/usr/bin/env python

import matplotlib.cm as cm
import matplotlib.pyplot as plt
import numpy
import pandas

from image_processing.segmentation.rasterise import Rasterise
from image_processing.segmentation.segmentation import SegmentVisitor

"""
This example will take a vector file and a base image file and rasterise the
polygons into an array based on the dimensions, projection and geolocation
info from the base image.

The end result is a segmented/labelled/enumerated array upon which the examples
given in binary_segmentation_example.py can used.
"""

img_fname = 'test_image'
vec_fname = 'rois.shp'

# Initialise the class with a base image file and the vector file
segments_ds = Rasterise(raster_fname=img_fname, vector_fname=vec_fname)

# Execute the rasterisation process
# The FID's of the vector file are used to define the segment/object ID
# but with a twist!
# FID's are positive, so to keep the same dynamic range of an array, without
# higher datatypes, use unsigned integers.
# As such 0 is assigned a fill value, and the segment ID's are the FID's + 1
# i.e. a FID of 10 has a segment ID of 11.
segments_ds.rasterise()

# Retrieve the segmented array
segmented_array = segments_ds.segmented_array

# Define an instance of a SegmentVisitor which we'll use to quickly traverse
# the segments
segments = SegmentVisitor(segmented_array)

# Create a random image with 10 time slices
rows, cols = segmented_array.shape
img = numpy.random.randint(0, 256, (10, rows, cols))

# Retrieve the bounding boxes; Note despite the shape of the segment/object
# a rectangular array is returned
bboxes = segments.segment_bounding_box()
print "Bounding box for segment 2"
print bboxes[2]
print "Subset of segmented_array for segment 2"
yse, xse = bboxes[2]
subset = segmented_array[yse[0]:yse[1], xse[0]:xse[1]]
plt.imshow(subset)
plt.title("Subset of segmented_array for segment 2")
plt.show()
plt.ion()

# An example of summarising by time (or band) then space
xbar_zaxis = numpy.mean(img, axis=0)
segment_stddev = segments.segment_stddev(xbar_zaxis)
print "Segment standard deviations of the mean time."
print segment_stddev

# Perform a sanity check that our method is correct
wh = segmented_array == 1
print "Sanity check for standard deviation of mean time for polygon 1:"
print numpy.std(xbar_zaxis[wh], ddof=1)

# An alternative method that will access only the pixels for a given segment
# and compute across time
idx = segments.get_segment_locations(1)
wh = numpy.zeros((rows, cols), dtype='bool')
wh[idx] = True
xbar_zaxis_id1 = numpy.mean(img[:, wh], axis=0)
segment_stddev2 = numpy.std(xbar_zaxis_id1, ddof=1)
print "Alternative method for standard deviation of the mean time:"
print segment_stddev2
print "Are methods 1 and 2 evaluate to the same value?"
print segment_stddev2 == segment_stddev[1]

# An example of summarising by space then time
# We'll get the available keys from the previous example and initialise a
# list to hold the mean value across time for each segment
xbar = {}
for key in segment_stddev:
    xbar[key] = []

# Calculate the mean value for each segment across time
for band in img:
    xbar_segments = segments.segment_mean(band)
    for key in xbar_segments:
        xbar[key].append(xbar_segments[key])

print "Mean value for each segment across time"
for key in xbar:
    print "Polgon {}:".format(key)
    print xbar[key]

# We'll use a pandas DataFrame to hold and display the data
df = pandas.DataFrame()
for key in xbar:
    k_str = 'Polygon {}'.format(key)
    df[k_str] = xbar[key]

# Now plot the mean value for each segment across time
# Note; If this doesn't plot launch ipython --pylab
df.plot()
plt.ion()
raw_input('Press Enter to exit')