shepseg

Python implementation of the image segmentation algorithm described by Shepherd et al [1]

Implemented using scikit-learn’s K-Means algorithm [2], and using numba [3] compiled code for the other main components.

Main entry point is the doShepherdSegmentation() function.

Examples

Read in a multi-band image as a single array, img, of shape (nBands, nRows, nCols). Ensure that any null pixels are all set to a known null value in all bands. Failure to correctly identify null pixels can result in a poorer quality segmentation.

>>> from pyshepseg import shepseg
>>> segRes = shepseg.doShepherdSegmentation(img, imgNullVal=nullVal)

The segimg attribute of the segRes object is an array of segment ID numbers, of shape (nRows, nCols).

Resulting segment ID numbers start from 1, and null pixels are set to zero.

Efficient segment location

After segmentation, the location of individual segments can be found efficiently using the object returned by makeSegmentLocations().

>>> segSize = shepseg.makeSegSize(segimg)
>>> segLoc = shepseg.makeSegmentLocations(segimg, segSize)

This segLoc object is indexed by segment ID number (must be of type shepseg.SegIdType), and each element contains information about the pixels which are in that segment. This information can be returned as a slicing object suitable to index the image array

>>> segNdx = segLoc[segId].getSegmentIndices()
>>> vals = img[0][segNdx]

This example would give an array of the pixel values from the first band of the original image, for the given segment ID.

This can be a very efficient way to calculate per-segment quantities. It can be used in pure Python code, or it can be used inside numba jit functions, for even greater efficiency.

References

1

Shepherd, J., Bunting, P. and Dymond, J. (2019). Operational Large-Scale Segmentation of Imagery Based on Iterative Elimination. Remote Sensing 11(6). https://www.mdpi.com/2072-4292/11/6/658

2

https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html

3

https://numba.pydata.org/

class pyshepseg.shepseg.RowColArray(length)

This data structure is used to store the locations of every pixel in a given segment. It will be used for entries in a jit dictionary. This means we can quickly find all the pixels belonging to a particular segment.

Attributes
idxint

Index of most recently added pixel

rowcolsuint32 ndarray (length, 2)

Row and col numbers of pixels in the segment

append(row, col)

Add the coordinates of a new pixel in the segment

Parameters
rowint

Row number of pixel

colint

Column number of pixel

getSegmentIndices()

Return the row and column numbers of the segment pixels as a tuple, suitable for indexing the image array. This supports selection of all pixels for a given segment.

class pyshepseg.shepseg.SegmentationResult

Results of the segmentation process

Attributes
segimgnumpy array (nRows, nCols)

Elements are segment ID numbers (starting from 1)

kmeanssklearn.cluster.KMeans

Fitted KMeans object

maxSpectralDifffloat

The value used to limit segment merging

singlePixelsEliminatedint

Number of single pixels merged to adjacent segments

smallSegmentsEliminatedint

Number of small segments merged into adjacent segments

pyshepseg.shepseg.applySpectralClusters(kmeansObj, img, imgNullVal)

Use the given KMeans object to predict spectral clusters on a whole image array.

Parameters
kmeansObjsklearn.cluster.KMeans

A fitted instance, as returned by fitSpectralClusters().

imgint ndarray (nBands, nRows, nCols)

The image to predict on

imgNullValint

Any pixels in img which have value imgNullVal will be set to SEGNULLVAL (i.e. zero) in the output cluster image.

Returns
segimgint ndarray (nRows, nCols)

The initial segment image, each element being the segment ID value for that pixel

pyshepseg.shepseg.autoMaxSpectralDiff(km, maxSpectralDiff, distPcntile)

Work out what to use as the maxSpectralDiff.

If current value is ‘auto’, then return the median spectral distance between cluster centres from the KMeans clustering object km.

If current value is None, return 10 times the largest distance between cluster centres (i.e. too large ever to make a difference)

Otherwise, return the given current value.

Parameters
kmsklearn.cluster.KMeans

KMeans clustering object

maxSpectralDiffstr or float

It is given in the units of the spectral space of img. If maxSpectralDiff is ‘auto’, a default value will be calculated from the spectral distances between cluster centres, as a percentile of the distribution of these (distPcntile). The value of distPcntile should be lowered when segementing an image with a larger range of spectral distances.

distPcntileint

See maxSpectralDiff

Returns
maxSpectralDiffint

The value to use as maxSpectralDiff.

pyshepseg.shepseg.buildSegmentSpectra(seg, img, maxSegId)

Build an array of the spectral statistics for each segment. Return an array of the sums of the spectral values for each segment, for each band

Parameters
segSegIdType ndarray (nRows, nCols)

Segmentation image

imgInteger ndarray (nBands, nRows, nCols)

Input multi-band image

maxSegIdint

Largest segment ID number in seg

Returns
spectSumfloat32 ndarray (numSegments+1, nBands)

Sums of all pixel values. Element [i, j] is the sum of all values in img for the j-th band, which have segment ID i. The row for i==0 is unused, as zero is not a valid segment ID.

pyshepseg.shepseg.clump(img, ignoreVal, fourConnected=True, clumpId=1)

Implementation of clumping using Numba.

Parameters
imgint ndarray (nRows, nCols)

Image array containing the data to be clumped.

ignoreValint

should be the “no data” value for the input

fourConnectedbool

If True, use 4-way connected, otherwise 8-way

clumpIdint

The start clump id to use

Returns
clumpimgSegIdType ndarray (nRows, nCols)

Image array containing the clump IDs for each pixel

clumpIdint

The highest clumpid used + 1

pyshepseg.shepseg.diagonalClusterCentres(xSample, numClusters)

Calculate an array of initial guesses at cluster centres. This will be given to the KMeans constructor as the init parameter.

The centres are evenly spaced along the diagonal of the bounding box of the data. The end points are placed 1 step in from the corners.

Parameters
xSampleint ndarray (numPoints, numBands)

A sample of data to be used for fitting

numClustersint

Number of cluster centres to be calculated

Returns
centresint ndarray (numPoints, numBands)

Initial cluster centres in spectral space

pyshepseg.shepseg.doMerge(segId, nbrSegId, seg, segSize, segLoc, spectSum)

Carry out a single segment merge. The segId segment is merged to the neighbouring nbrSegId. Modifies seg, segSize, segLoc and spectSum in place.

Parameters
segIdSegIdType

Segment ID of the segment to be merged. Modified in place

nbrSegIdSegIdType

Segment ID of the segment into which segId will be merged. Modified in place

segSegIdType ndarray (nRows, nCols)

Segment ID image array

segSizeint ndarray (numSegments+1, )

Counts of pixels in each segment, indexed by segment ID. Modified in place with new counts for both segments

segLocnumba.typed.Dict

Dictionary of per-segment pixel coordinates. As computed by makeSegmentLocations()

spectSumfloat32 ndarray (numSegments+1, nBands)

Sums of all pixel values. As computed by buildSegmentSpectra(). Updated in place with new sums for both segments

pyshepseg.shepseg.doShepherdSegmentation(img, numClusters=60, clusterSubsamplePcnt=1, minSegmentSize=50, maxSpectralDiff='auto', imgNullVal=None, fourConnected=True, verbose=False, fixedKMeansInit=False, kmeansObj=None, spectDistPcntile=50)

Perform Shepherd segmentation in memory, on the given multi-band img array.

Parameters
imginteger ndarray of shape (nBands, nRows, nCols)
numClustersint

Number of clusters to create with k-means clustering

clusterSubsamplePcntint

Passed to fitSpectralClusters(). See there for details

minSegmentSizeint

the minimum segment size (in pixels) which will be left after eliminating small segments (except for segments which cannot be eliminated).

maxSpectralDiffstr or float

sets a limit on how different segments can be and still be merged. It is given in the units of the spectral space of img. If maxSpectralDiff is ‘auto’, a default value will be calculated from the spectral distances between cluster centres, as a percentile of the distribution of these (spectDistPcntile). The value of spectDistPcntile should be lowered when segementing an image with a larger range of spectral distances.

spectDistPcntileint

See maxSpectralDiff

fourConnectedbool

If True, use 4-way connectedness when clumping, otherwise use 8-way

imgNullValint or None

If not None, then pixels with this value in any band are set to zero (SEGNULLVAL) in the output segmentation. If there are null values in the image array, it is important to give this null value, as it can stringly affect the initial spectral clustering, which in turn strongly affects the final segmenation.

fixedKMeansInitbool

If fixedKMeansInit is True, then choose a fixed set of cluster centres to initialize the KMeans algorithm. This can be useful to provide strict determinacy of the results by avoiding sklearn’s multiple random initial guesses. The default is to allow sklearn to guess, which is good for avoiding local minima.

kmeansObjsklearn.cluster.KMeans object

By default, the spectral clustering step will be fitted using the given img. However, if kmeansObj is not None, it is taken to be a fitted instance of sklearn.cluster.KMeans, and will be used instead. This is useful when enforcing a consistent clustering across multiple tiles (see the pyshepseg.tiling module for details).

Returns
segResultSegmentationResult object

Notes

Default values are mostly as suggested by Shepherd et al.

Segment ID numbers start from 1. Zero is a NULL segment ID.

The return value is an instance of SegmentationResult class.

pyshepseg.shepseg.eliminateSinglePixels(img, seg, segSize, minSegId, maxSegId, fourConnected)

Approximate elimination of single pixels, as suggested by Shepherd et al (section 2.3, page 6). This step suggested as an efficient way of removing a large number of segments which are single pixels, by approximating the spectrally-nearest neighbouring segment with the spectrally-nearest neighouring pixel.

Parameters
imgint ndarray (nBands, nRows, nCols)

The original spectral image

segSegIdType ndarray (nRows, nCols)

The image of segment IDs

segSizeint array (numSeg+1, )

Array of pixel counts for every segment

minSegIdSegIdType

Smallest segment ID

maxSegIdSegIdType

Largest segment ID

fourConnectedbool

If True use 4-way connectedness, otherwise 8-way

Notes

Segment ID numbers start at 1 (i.e. 0 is not valid)

Modifies seg array in place.

pyshepseg.shepseg.eliminateSmallSegments(seg, img, maxSegId, minSegSize, maxSpectralDiff, fourConnected, minSegId)

Eliminate small segments. Start with smallest, and merge them into spectrally most similar neighbour. Repeat for larger segments.

Modifies seg array in place.

Parameters
segSegIdType ndarray (nRows, nCols)

Segment ID image array. Modified in place as segments are merged.

imgInteger ndarray (nBands, nRows, nCols)

Input multi-band image

maxSegIdSegIdType

Largest segment ID number in seg

minSegSizeint

Size (in pixels) of the smallest segment which will NOT be eliminated

maxSpectralDifffloat

Limit on how different segments can be and still be merged. It is given in the units of the spectral space of img.

fourConnectedbool

If True, use four-way connectedness to judge neighbours, otherwise use eight-way.

minSegIdSegIdType

Minimum valid segment ID number

Returns
numEliminatedint

Number of segments eliminated

pyshepseg.shepseg.findMergeSegment(segId, segLoc, seg, segSize, spectSum, maxSpectralDiff, fourConnected)

For the given segId, find which neighboring segment it should be merged with. The chosen merge segment is the one which is spectrally most similar to the given one, as measured by minimum Euclidean distance in spectral space.

Called by eliminateSmallSegments().

Parameters
segIdSegIdType

Segment ID number of segment to merge

segLocnumba.typed.Dict

Dictionary of per-segment pixel coordinates. As computed by makeSegmentLocations()

segSegIdType ndarray (nRows, nCols)

Segment ID image array

segSizeint ndarray (numSegments+1, )

Counts of pixels in each segment, indexed by segment ID

spectSumfloat32 ndarray (numSegments+1, nBands)

Sums of all pixel values. As computed by buildSegmentSpectra()

maxSpectralDifffloat

Limit on how different segments can be and still be merged. It is given in the units of the spectral space of img

fourConnectedbool

If True, use four-way connectedness to judge neighbours, otherwise use eight-way

pyshepseg.shepseg.findNearestNeighbourPixel(img, seg, i, j, segSize, fourConnected)

For the (i, j) pixel, choose which of the neighbouring pixels is the most similar, spectrally.

Returns the row and column of the most spectrally similar neighbour, which is also in a clump of size > 1. If none is found, return (-1, -1)

Parameters
imgint ndarray (nBands, nRows, nCols)

Input multi-band image

segSegIdType ndarray (nRows, nCols)

Partially completed segmentation image (values are segment ID numbers)

iint

Row number of target pixel

jint

Column number of target pixel

segSizeint ndarray (numSegments+1, )

Pixel counts, indexed by segment ID number (i.e. a histogram of the seg array)

fourConnectedbool

If True, use four-way connectedness to judge neighbours, otherwise use eight-way.

Returns
iiint

Row number of the selected neighbouring pixel (-1 if not found)

jjint

Column number of the selected neighbouring pixel (-1 if not found)

pyshepseg.shepseg.fitSpectralClusters(img, numClusters, subsamplePcnt, imgNullVal, fixedKMeansInit)

First step of Shepherd segmentation. Use K-means clustering to create a set of “seed” segments, labelled only with their spectral cluster number.

Parameters
imgint ndarray (nBands, nRows, nCols).
numClustersint

The number of clusters for the KMeans algorithm to find (i.e. it is ‘k’)

subsamplePcntint

The percentage of the pixels to actually use for KMeans clustering. Shepherd et al find that only a very small percentage is required.

imgNullValint or None

If imgNullVal is not None, then pixels in img with this value in any band are set to segNullVal in the output.

fixedKMeansInitbool

If True, then use a simple algorithm to determine the fixed set of initial cluster centres. Otherwise allow the sklearn routine to choose its own initial guesses.

Returns
kmeansObjsklearn.cluster.KMeans

A fitted object of class sklearn.cluster.KMeans. This is suitable to use with the applySpectralClusters() function.

pyshepseg.shepseg.makeSegSize(seg)

Return an array of segment sizes, essentially a histogram for the segment ID values.

Parameters
segSegIdType ndarray (nRows, nCols)

Image array of segment ID values

Returns
segSizeint ndarray (numSegments+1, )

Array is indexed by segment ID. Each element is the number of pixels in that segment.

pyshepseg.shepseg.makeSegmentLocations(seg, segSize)

Create a data structure to hold the locations of all pixels in all segments.

Parameters
segSegIdType ndarray (nRows, nCols)

Segment ID image array

segSizeint ndarray (numSegments+1, )

Counts of pixels in each segment, indexed by segment ID

Returns
segLocnumba.typed.Dict

Indexed by segment ID number, each entry is a RowColArray object, giving the pixel coordinates of all pixels for that segment

pyshepseg.shepseg.mergeSinglePixels(img, seg, segSize, segToElim, fourConnected)

Search for single-pixel segments, and decide which neighbouring segment they should be merged with. Finds all to eliminate, then performs merge on all selected. Modifies seg and segSize arrays in place, and returns the number of segments eliminated.

Parameters
imgint ndarray (nBands, nRows, nCols)

the original spectral image

segint ndarray (nRows, nCols)

the image of segments

segSizeint array (numSeg+1, )

Array of pixel counts for every segment

segToElimint ndarray (3, maxSegId)

Temporary storage for segments to be eliminated

fourConnectedbool

If True use 4-way connectedness, otherwise 8-way

Returns
numEliminatedint

Number of segments eliminated

pyshepseg.shepseg.relabelSegments(seg, segSize, minSegId)

The given seg array is an image of segment labels, with some numbers unused, due to elimination of small segments. Go through and find the unused numbers, and re-label segments above these so that segment labels are contiguous.

Modifies the seg array in place. The segSize array is not updated, and should be recomputed.

Parameters
segSegIdType ndarray (nRows, nCols)

Segmentation image. Updated in place with new segment ID values

segSizeint array (numSeg+1, )

Array of pixel counts for every segment

minSegIdint

Smallest valid segment ID number