Peak detection in a 2D array

I detected the peaks using a local maximum filter. Here is the result on your first dataset of 4 paws:

I also ran it on the second dataset of 9 paws and it worked as well.

Here is how you do it:

import numpy as npfrom scipy.ndimage.filters import maximum_filterfrom scipy.ndimage.morphology import generate_binary_structure, binary_erosionimport matplotlib.pyplot as pp#for some reason I had to reshape. Numpy ignored the shape header.paws_data = np.loadtxt("paws.txt").reshape(4,11,14)#getting a list of imagespaws = [p.squeeze() for p in np.vsplit(paws_data,4)]def detect_peaks(image):    """    Takes an image and detect the peaks usingthe local maximum filter.    Returns a boolean mask of the peaks (i.e. 1 when    the pixel's value is the neighborhood maximum, 0 otherwise)    """    # define an 8-connected neighborhood    neighborhood = generate_binary_structure(2,2)    #apply the local maximum filter; all pixel of maximal value     #in their neighborhood are set to 1    local_max = maximum_filter(image, footprint=neighborhood)==image    #local_max is a mask that contains the peaks we are     #looking for, but also the background.    #In order to isolate the peaks we must remove the background from the mask.    #we create the mask of the background    background = (image==0)    #a little technicality: we must erode the background in order to     #successfully subtract it form local_max, otherwise a line will     #appear along the background border (artifact of the local maximum filter)    eroded_background = binary_erosion(background, structure=neighborhood, border_value=1)    #we obtain the final mask, containing only peaks,     #by removing the background from the local_max mask (xor operation)    detected_peaks = local_max ^ eroded_background    return detected_peaks#applying the detection and plotting resultsfor i, paw in enumerate(paws):    detected_peaks = detect_peaks(paw)    pp.subplot(4,2,(2*i+1))    pp.imshow(paw)    pp.subplot(4,2,(2*i+2) )    pp.imshow(detected_peaks)pp.show()

All you need to do after is use scipy.ndimage.measurements.label on the mask to label all distinct objects. Then you'll be able to play with them individually.

Note that the method works well because the background is not noisy. If it were, you would detect a bunch of other unwanted peaks in the background. Another important factor is the size of the neighborhood. You will need to adjust it if the peak size changes (the should remain roughly proportional).

Solution

Data file: paw.txt. Source code:

from scipy import *from operator import itemgettern = 5  # how many fingers are we looking ford = loadtxt("paw.txt")width, height = d.shape# Create an array where every element is a sum of 2x2 squares.fourSums = d[:-1,:-1] + d[1:,:-1] + d[1:,1:] + d[:-1,1:]# Find positions of the fingers.# Pair each sum with its position number (from 0 to width*height-1),pairs = zip(arange(width*height), fourSums.flatten())# Sort by descending sum value, filter overlapping squaresdef drop_overlapping(pairs):    no_overlaps = []    def does_not_overlap(p1, p2):        i1, i2 = p1[0], p2[0]        r1, col1 = i1 / (width-1), i1 % (width-1)        r2, col2 = i2 / (width-1), i2 % (width-1)        return (max(abs(r1-r2),abs(col1-col2)) >= 2)    for p in pairs:        if all(map(lambda prev: does_not_overlap(p,prev), no_overlaps)):            no_overlaps.append(p)    return no_overlapspairs2 = drop_overlapping(sorted(pairs, key=itemgetter(1), reverse=True))# Take the first n with the heighest valuespositions = pairs2[:n]# Print resultsprint d, "\n"for i, val in positions:    row = i / (width-1)    column = i % (width-1)    print "sum = %f @ %d,%d (%d)" % (val, row, column, i)    print d[row:row+2,column:column+2], "\n"

Output without overlapping squares. It seems that the same areas are selected as in your example.

Some comments

The tricky part is to calculate sums of all 2x2 squares. I assumed you need all of them, so there might be some overlapping. I used slices to cut the first/last columns and rows from the original 2D array, and then overlapping them all together and calculating sums.

To understand it better, imaging a 3x3 array:

>>> a = arange(9).reshape(3,3) ; aarray([[0, 1, 2],       [3, 4, 5],       [6, 7, 8]])

Then you can take its slices:

>>> a[:-1,:-1]array([[0, 1],       [3, 4]])>>> a[1:,:-1]array([[3, 4],       [6, 7]])>>> a[:-1,1:]array([[1, 2],       [4, 5]])>>> a[1:,1:]array([[4, 5],       [7, 8]])

Now imagine you stack them one above the other and sum elements at the same positions. These sums will be exactly the same sums over the 2x2 squares with the top-left corner in the same position:

>>> sums = a[:-1,:-1] + a[1:,:-1] + a[:-1,1:] + a[1:,1:]; sumsarray([[ 8, 12],       [20, 24]])

When you have the sums over 2x2 squares, you can use max to find the maximum, or sort, or sorted to find the peaks.

To remember positions of the peaks I couple every value (the sum) with its ordinal position in a flattened array (see zip). Then I calculate row/column position again when I print the results.

Notes

I allowed for the 2x2 squares to overlap. Edited version filters out some of them such that only non-overlapping squares appear in the results.

Choosing fingers (an idea)

Another problem is how to choose what is likely to be fingers out of all the peaks. I have an idea which may or may not work. I don't have time to implement it right now, so just pseudo-code.

I noticed that if the front fingers stay on almost a perfect circle, the rear finger should be inside of that circle. Also, the front fingers are more or less equally spaced. We may try to use these heuristic properties to detect the fingers.

Pseudo code:

select the top N finger candidates (not too many, 10 or 12)consider all possible combinations of 5 out of N (use itertools.combinations)for each combination of 5 fingers:    for each finger out of 5:        fit the best circle to the remaining 4        => position of the center, radius        check if the selected finger is inside of the circle        check if the remaining four are evenly spread        (for example, consider angles from the center of the circle)        assign some cost (penalty) to this selection of 4 peaks + a rear finger        (consider, probably weighted:             circle fitting error,             if the rear finger is inside,             variance in the spreading of the front fingers,             total intensity of 5 peaks)choose a combination of 4 peaks + a rear peak with the lowest penalty

This is a brute-force approach. If N is relatively small, then I think it is doable. For N=12, there are C_12^5 = 792 combinations, times 5 ways to select a rear finger, so 3960 cases to evaluate for every paw.

This is an image registration problem. The general strategy is:

• Have a known example, or some kind of prior on the data.
• Fit your data to the example, or fit the example to your data.
• It helps if your data is roughly aligned in the first place.

Here's a rough and ready approach, "the dumbest thing that could possibly work":

• Start with five toe coordinates in roughly the place you expect.
• With each one, iteratively climb to the top of the hill. i.e. given current position, move to maximum neighbouring pixel, if its value is greater than current pixel. Stop when your toe coordinates have stopped moving.

To counteract the orientation problem, you could have 8 or so initial settings for the basic directions (North, North East, etc). Run each one individually and throw away any results where two or more toes end up at the same pixel. I'll think about this some more, but this kind of thing is still being researched in image processing - there are no right answers!

Slightly more complex idea: (weighted) K-means clustering. It's not that bad.

• Start with five toe coordinates, but now these are "cluster centres".

Then iterate until convergence:

• Assign each pixel to the closest cluster (just make a list for each cluster).
• Calculate the center of mass of each cluster. For each cluster, this is: Sum(coordinate * intensity value)/Sum(coordinate)
• Move each cluster to the new centre of mass.

This method will almost certainly give much better results, and you get the mass of each cluster which may help in identifying the toes.

(Again, you've specified the number of clusters up front. With clustering you have to specify the density one way or another: Either choose the number of clusters, appropriate in this case, or choose a cluster radius and see how many you end up with. An example of the latter is mean-shift.)

Sorry about the lack of implementation details or other specifics. I would code this up but I've got a deadline. If nothing else has worked by next week let me know and I'll give it a shot.