1. Joe Kington
  2. Paw Analysis

Source

Paw Analysis / analysis.py

import numpy as np

PAW_CODE_LUT = {0:'LF', 1:'RH', 2:'RF', 3:'LH'}
INV_PAW_CODE_LUT = {'LF':0, 'RH':1, 'RF':2, 'LH':3}
STD_NUMX, STD_NUMY = 20, 20

def find_paws(data):
    """Detects and isolates contiguous regions in the input array"""
    from scipy.ndimage import uniform_filter, label, find_objects
    from scipy.ndimage.morphology import binary_fill_holes 
    from scipy.ndimage.morphology import generate_binary_structure

    data = uniform_filter(data, 4)
    # Threshold the blurred data (this needs to be a bit > 0 due to the blur)
    thresh = data > (0.05 * data.max())
    # Fill any interior holes in the paws to get cleaner regions...
    filled = binary_fill_holes(thresh)
    # Label each contiguous paw
    coded_paws, num_paws = label(filled)
    # Isolate the extent of each paw
    data_slices = find_objects(coded_paws)
    return data_slices

def standardize_paw(paw):
    """Standardizes a pawprint onto a STD_NUMYxSTD_NUMX grid. Returns a 1D, 
    flattened version of the paw data resample onto this grid."""
    from scipy.ndimage import map_coordinates
    paw = paw.sum(axis=2)
    ny, nx = paw.shape

    # Trim off any "blank" edges around the paw...
    mask = paw > 0.01 * paw.max()
    y, x = np.mgrid[:ny, :nx]
    ymin, ymax = y[mask].min(), y[mask].max() 
    xmin, xmax = x[mask].min(), x[mask].max() 

    # Make a 20x20 grid to resample the paw pressure values onto
    xi = np.linspace(xmin, xmax, STD_NUMX)
    yi = np.linspace(ymin, ymax, STD_NUMY)
    xi, yi = np.meshgrid(xi, yi)

    # Resample the values onto the 20x20 grid
    coords = np.vstack([yi.flatten(), xi.flatten()])
    zi = map_coordinates(paw, coords)

    # Rescale the pressure values
    zi -= zi.min()
    zi /= zi.max()
    zi -= zi.mean() #<- Helps distinguish front from hind paws...
    return zi.flatten()

def pattern_classification(data_slices):   
    """Attempts to classify a sequence of paw impacts based on their spatial
    and temporal pattern alone. Input "data_slices" must be sorted by initial
    contact time!"""
    #TODO: Add check to see if the paw before the last impact
    #      landed off the sensor...

    # Get the centroid for each paw impact...
    paw_coords = []
    for x,y,z in data_slices:
        paw_coords.append([(item.stop + item.start) / 2.0 for item in (x,y)])
    paw_coords = np.array(paw_coords)

    # Make a vector between each sucessive impact...
    dx, dy = np.diff(paw_coords, axis=0).T

    #-- Group paws -------------------------------------------
    paw_number = np.arange(len(paw_coords))

    # Did we miss the hind paw impact after the first 
    # front paw impact? If so, first dx will be positive...
    if dx[0] > 0: 
        paw_number[1:] += 1

    # Are we starting with the left or right front paw...
    # We assume we're starting with the left, and check dy[0].
    # If dy[0] > 0 (i.e. the next paw impacts to the left), then
    # it's actually the right front paw, instead of the left.
    if dy[0] > 0: # Right front paw impact...
        paw_number += 2
            
    paw_codes = paw_number % 4
    paw_labels = [PAW_CODE_LUT[code] for code in paw_codes]

    problems = paw_pattern_problems(paw_labels, dx, dy)
    return paw_labels, problems

def paw_pattern_problems(paw_labels, dx, dy):
    """Check whether or not the label sequence "paw_labels" conforms to our
    expected spatial pattern of paw impacts."""
    # Check for problems... (This could be written a _lot_ more cleanly...)
    problems = False
    last = paw_labels[0]
    for paw, dy, dx in zip(paw_labels[1:], dy, dx):
        # Going from a left paw to a right, dy should be negative
        if last.startswith('L') and paw.startswith('R') and (dy > 0):
            problems = True
            break
        # Going from a right paw to a left, dy should be positive
        if last.startswith('R') and paw.startswith('L') and (dy < 0):
            problems = True
            break
        # Going from a front paw to a hind paw, dx should be negative
        if last.endswith('F') and paw.endswith('H') and (dx > 0):
            problems = True
            break
        # Going from a hind paw to a front paw, dx should be positive
        if last.endswith('H') and paw.endswith('F') and (dx < 0):
            problems = True
            break
        last = paw
    return problems

class ClassificationSystem(object):
    """Holds the data and methods necessary to project a paw using an
    eigenimage-based method. The basis vectors, etc are stored on disk
    as numpy saved arrays (.npy)."""
    def __init__(self):
        try:
            # Try to load the approriate data files
            self.load_data()
        except IOError:
            # If the data files aren't there, just issue a warning...
            # (presumably this is a first run, and they need to be built)
            import warnings
            warnings.warn('Classification data files not found!', 
                    warnings.RuntimeWarning)

    def __call__(self, paw):
        """Calls self.classify"""
        return self.classify(paw)

    def load_data(self):
        # This is a seperate method so that the data may be reloaded on the
        # fly if it's changed...
        import os
        basedir = os.path.dirname(__file__)
        load = lambda filename: np.load(os.path.join(basedir, filename))
        self.template_paws = load('template_paws.npy')
        self.average_paw = load('average_paw.npy')
        self.basis_stds = load('basis_stds.npy')
        self.basis_vecs = load('basis_vecs.npy')

    def classify(self, paw):
        """Classifies a (standardized) pawprint based on how close its eigenpaw
        score is to known eigenpaw scores of the paws for each leg. Returns a code
        of "LF", "LH", "RF", or "RH" for the left front, left hind, etc paws."""
        # Subtract the average paw (known a-priori from the training dataset)
        paw -= self.average_paw
        # Project the paw into eigenpaw-space
        scores = paw.dot(self.basis_vecs) 
        # "whiten" the score so that all dimensions are equally important
        scores /= self.basis_stds
        # Select which template paw is the closest to the given paw...
        diff = self.template_paws - scores
        diff *= diff
        diff = np.sqrt(diff.sum(axis=1))
        # Return the ascii code for the closest template paw.
        return PAW_CODE_LUT[diff.argmin()]

# Create an instance of ClassificationSystem, so that you can just call
# analysis.classify(paw)
classify = ClassificationSystem()