Source

Paw Analysis / analysis.py

"""
Paw impact analysis and classification algorithms.
"""

__author__ = 'Joe Kington <jkington@geology.wisc.edu>'
__license__ = """
Copyright (C) 2010 by The Free Software Foundation

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""

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 = np.dot(paw, 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()