## Peak detection in a 2D array

I'm helping a veterinary clinic measuring pressure under a dogs paw. I use Python for my data analysis and now I'm stuck trying to divide the paws into (anatomical) subregions.

I made a 2D array of each paw, that consists of the maximal values for each sensor that has been loaded by the paw over time. Here's an example of one paw, where I used Excel to draw the areas I want to 'detect'. These are 2 by 2 boxes around the sensor with local maxima's, that together have the largest sum.

So I tried some experimenting and decide to simply look for the maximums of each column and row (can't look in one direction due to the shape of the paw). This seems to 'detect' the location of the separate toes fairly well, but it also marks neighboring sensors.

So what would be the best way to tell Python which of these maximums are the ones I want?

**Note: The 2x2 squares can't overlap, since they have to be separate toes!**

Also I took 2x2 as a convenience, any more advanced solution is welcome, but I'm simply a human movement scientist, so I'm neither a real programmer or a mathematician, so please keep it 'simple'.

Here's a version that can be loaded with `np.loadtxt`

##### Results

So I tried @jextee's solution (see the results below). As you can see, it works very on the front paws, but it works less well for the hind legs.

More specifically, it can't recognize the small peak that's the fourth toe. This is obviously inherent to the fact that the loop looks top down towards the lowest value, without taking into account where this is.

Would anyone know how to tweak @jextee's algorithm, so that it might be able to find the 4th toe too?

Since I haven't processed any other trials yet, I can't supply any other samples. But the data I gave before were the averages of each paw. This file is an array with the maximal data of 9 paws in the order they made contact with the plate.

This image shows how they were spatially spread out over the plate.

##### Update:

**I have set up a blog for anyone interested** and I have setup a SkyDrive with all the raw measurements. So to anyone requesting more data: more power to you!

##### New update:

So after the help I got with my questions regarding paw detection and paw sorting, I was finally able to check the toe detection for every paw! Turns out, it doesn't work so well in anything but paws sized like the one in my own example. Off course in hindsight, it's my own fault for choosing the 2x2 so arbitrarily.

Here's a nice example of where it goes wrong: a nail is being recognized as a toe and the 'heel' is so wide, it gets recognized twice!

The paw is too large, so taking a 2x2 size with no overlap, causes some toes to be detected twice. The other way around, in small dogs it often fails to find a 5th toe, which I suspect is being caused by the 2x2 area being too large.

After trying the current solution on all my measurements I came to the staggering conclusion that for nearly all my small dogs it didn't find a 5th toe and that in over 50% of the impacts for the large dogs it would find more!

So clearly I need to change it. My own guess was changing the size of the `neighborhood`

to something smaller for small dogs and larger for large dogs. But `generate_binary_structure`

wouldn't let me change the size of the array.

Therefore, I'm hoping that anyone else has a better suggestion for locating the toes, perhaps having the toe area scale with the paw size?

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 np from scipy.ndimage.filters import maximum_filter from scipy.ndimage.morphology import generate_binary_structure, binary_erosion import 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 images paws = [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 results for 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).

**Find a peak element in a 2D array,** Finding peak element in a 2D Array. #include <bits/stdc++.h>. using namespace std;. const int MAX = 100;. // Function to find the maximum in column 'mid'. Peak Element in 2D array Divide and Conquer Algorithms Data Structure Algorithms An item is said to be a peak element when it is greater than or equal with all four neighbor of that element. The neighbor elements are the top, bottom, left and right elements.

##### Solution

Data file: paw.txt. Source code:

from scipy import * from operator import itemgetter n = 5 # how many fingers are we looking for d = 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 squares def 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_overlaps pairs2 = drop_overlapping(sorted(pairs, key=itemgetter(1), reverse=True)) # Take the first n with the heighest values positions = pairs2[:n] # Print results print 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) ; a array([[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:]; sums array([[ 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.

**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: Peaks detection result. I also ran it on the� As it turns out using the peak detection algorithm from this question Peak detection in a 2D array only complicates matters. After applying the Gaussian filter to the image all that needs to be done is to ask for the maximum bin (as Bi Rico pointed out) and then obtain the maximum in coordinates.

**Find all peaks in 2d array,** I tried peak detector vi's but they all seem to work with 1D signal only (I cannot change my array to 1D, cause the only condition Peak detection in a 2D array. Peak detection in a 2D array I'm helping a veterinary clinic measuring pressure under a dogs paw. I use Python for my data analysis and now I'm stuck trying to divide the paws into (anatomical) subregions. I made a 2D array of each paw, that consists of the maximal values for each sensor that has been loaded by the paw over time.

Using persistent homology to analyze your data set I get the following result (click to enlarge):

This is the 2D-version of the peak detection method described in this SO answer. The above figure simply shows 0-dimensional persistent homology classes sorted by persistence.

I did upscale the original dataset by a factor of 2 using scipy.misc.imresize(). However, note that I did consider the four paws as one dataset; splitting it into four would make the problem easier.

**Methodology.**
The idea behind this quite simple: Consider the function graph of the function that assigns each pixel its level. It looks like this:

Now consider a water level at height 255 that continuously descents to lower levels. At local maxima islands pop up (birth). At saddle points two islands merge; we consider the lower island to be merged to the higher island (death). The so-called persistence diagram (of the 0-th dimensional homology classes, our islands) depicts death- over birth-values of all islands:

The *persistence* of an island is then the difference between the birth- and death-level; the vertical distance of a dot to the grey main diagonal. The figure labels the islands by decreasing persistence.

The very first picture shows the locations of births of the islands. This method not only gives the local maxima but also quantifies their "significance" by the above mentioned persistence. One would then filter out all islands with a too low persistence. However, in your example every island (i.e., every local maximum) is a peak you look for.

Python code can be found here.

**[PDF] 6.006 Introduction to Algorithms Lecture 2: Peak Finding,** Today. • Peak finding (new problem) So can throw away rest of array, 2D Peak Finding. • Given matrix of numbers. • Want an entry not smaller than its (up to). 1) Search for the highest pixel. Once you have that, search around that for the best fit for 2×2 (maybe maximizing the 2×2 sum), or do a 2d gaussian fit inside the sub region of say 4×4 centered on the highest pixel. Then set those 2×2 pixels you have found to zero (or maybe 3×3) around the peak center.

**CodeChef,** python - Peak detection in a 2D array - Stack Overflow. I'm helping a veterinary clinic measuring pressure under a dogs paw. I use Python for my data analysis� I assume a 2D peak is a value that is higher than its neighbors. Also we may define that values at both ends of the array only have one neighbor, and they will be peaks if they are higher than that neighbor. This algorithm will be minimally of O (n) complexity since all elements of the array must be examined in order to find the peaks.

**[PDF] Lecture 20: Peak Finding in 2D - COMS10007,** Peak Finding. Let A = a0,a1,,an-1 be an array of integers of length n. 0. 1. 2. 3. 4 . 5. 6. 7. 8. 9 a0 a1 a2 a3 a4 a5 a6 a7 a8 a9. Definition: (Peak). Hi all, I was looking for a method to find peaks over a 2d array data set (matrix).I found some topics in old threads but noone seems exactly to help me. In fact I was looking for something that returns (i,j) indexes of maximums found over a 2d data array (my seeking area=window).This means that maybe I should consider the derivatives over 2d data (or maybe some 2d surface evaluation for

**Peak detection in a 2D array,** Peak detection in a 2D array. I'm helping a veterinary clinic measuring pressure under a dogs paw. I use Python for my data analysis and now I'm stuck trying to� Basically it searches in the 2D arrays for peaks, and returns back an array of cluster, each cluster contains two integer value (row and column number) where there was a peak detected. Best regards, I.R.

**Finding a peak in 2 Dimensional Array,** A peak in a 2-D array is an element which has left, right, top and bottom elements lower than it. A simple approach will be to extend the 1-D array� The code analyzes noisy 2D images and find peaks using robust local maxima finder (1 pixel resolution) or by weighted centroids (sub-pixel resolution). The code is designed to be as fast as possible, so I kept it pretty basic. It best works when using uint16 \ uint8 images, and assumes that peaks are relatively sparse.