Détection de crête dans un tableau 2D

Solution

Fichier de données: paw.txt. Code Source:

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"

Sortie sans carrés qui se chevauchent. Il semble que les mêmes zones sont sélectionnées comme dans votre exemple.

Certains commentaires

La partie la plus délicate est de calculer la somme de tous les carrés de 2 x 2. Je suppose que vous avez besoin de tous d'entre eux, alors il pourrait y avoir un certain chevauchement. J'ai utilisé des tranches de couper le premier/dernier les colonnes et les lignes de l'original tableau 2D, puis en les superposant tous ensemble et en calculant les sommes.

Pour mieux la comprendre, de l'imagerie d'une matrice 3x3:

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

Alors vous pouvez prendre ses tranches:

>>> 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]])

Maintenant, imaginez que vous empilez-les les unes sur les autres et la somme des éléments à la même position. Ces sommes seront exactement les mêmes sommes sur la 2x2 places avec le coin en haut à gauche dans la même position:

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

Lorsque vous avez le sommes plus 2x2 places, vous pouvez utiliser max pour trouver le maximum, ou sort, ou sorted pour trouver les pics.

De se rappeler les positions des pics j'chaque couple de valeur (de la somme) avec sa position ordinale dans une stagnation de la matrice (voir zip). Puis-je calculer la ligne/colonne de la position encore quand j'ai imprimer les résultats.

Notes

- Je le droit de la 2x2 places à se chevaucher. Version modifiée d'un filtre certains d'entre eux tels que seuls les non-cumul des carrés apparaissent dans les résultats.

Le choix de doigts (une idée)

Un autre problème est de savoir comment choisir ce qui est susceptible d'être les doigts de tous les sommets. J'ai une idée qui peut ou peut ne pas fonctionner. Je n'ai pas de temps pour la mettre en œuvre dès maintenant, il suffit donc de pseudo-code.

J'ai remarqué que si le front de doigts sur presque un cercle parfait, l'arrière doigt doit être à l'intérieur de ce cercle. Aussi, le front, les doigts sont plus ou moins espacés de façon égale. On peut essayer d'utiliser ces propriétés heuristiques pour détecter les doigts.

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

C'est une force brute approche. Si N est relativement petite, je pense que c'est faisable. Pour N=12, il y a C_12^5 = 792 combinaisons, fois 5 façons de sélectionner un arrière doigt, alors 3960 cas à évaluer pour chaque patte.

Ce problème a été étudié en profondeur par les physiciens. Il y a une bonne mise en œuvre dans RACINE. Regardez les TSpectrum classes (surtout TSpectrum2 pour votre cas) et de la documentation pour eux.

Références:

1. M. Morhac et al.: Arrière-plan de l'élimination méthodes pour multidimensionnelle coïncidence gamma-ray spectres. Nucléaire, Méthodes et Instruments de Recherche en Physique 401 (1997) 113-132.
2. M. Morhac et al.: Efficace à une et deux dimensions de l'Or déconvolution et de son application à la gamma-ray spectres de décomposition. Nucléaire, Méthodes et Instruments de Recherche en Physique 401 (1997) 385-408.
3. M. Morhac et al.: L'Identification des pics dans multidimensionnelle coïncidence gamma-ray spectres. Nucléaire, Méthodes et Instruments de Recherche en Physique A 443(2000), 108-125.

...et pour ceux qui n'ont pas accès à une souscription de NIM:

• Spectrum.doc
• SpectrumDec.ps.gz
• SpectrumSrc.ps.gz
• SpectrumBck.ps.gz

À l'aide de l'homologie persistante à l'analyse de votre ensemble de données, j'obtiens le résultat suivant (cliquez pour agrandir):

C'est de la 2D version de la pointe de détection de la méthode décrite dans cette DONC réponse. La figure ci-dessus montre, tout simplement, 0-dimensions homologie persistante classes triés par la persévérance.

J'ai fait haut de gamme du jeu de données d'origine par un facteur de 2 à l'aide de scipy.misc.imresize(). Toutefois, notez que je n'ai examiner les quatre pattes comme un dataset; de le diviser en quatre, faire plus facilement le problème.

Méthodologie.
L'idée derrière ce assez simple: Considérons la fonction graphique de la fonction qui attribue à chaque pixel de son niveau. Il ressemble à ceci:

Considérons maintenant un niveau d'eau à hauteur de 255 continuellement les descentes vers les niveaux inférieurs. Au maxima locaux les îles pop-up (de la naissance). En selle points deux îles de fusion; nous considérons que le plus faible de l'île de fusion le plus élevé de l'île (la mort). Le soi-disant persistance diagramme (0-th dimensions homologie de classes, nos îles) représente la mort - dessus de la naissance des valeurs de toutes les îles:

La persistance de l'île est alors la différence entre la naissance et de la mort de niveau; la distance verticale d'un point en passant par le gris de la diagonale principale. La figure étiquettes les îles par la diminution de la persistance.

La première photo montre les lieux de naissances les îles. Cette méthode ne donne pas seulement les maxima locaux, mais aussi de quantifier leur "signification" par le ci-dessus mentionné persistance. On aurait alors de filtrer toutes les îles avec une trop faible persistance. Cependant, dans votre exemple, chaque île (c'est à dire, chaque maximum local) est un pic vous recherchez.

Code Python peut être trouvé ici.

Je suis sûr que vous avez assez d'aller sur maintenant, mais je ne peux pas aider, mais vous suggérons d'utiliser le clustering k-means méthode. k-means est un algorithme de clustering non supervisé qui vous permettra de données (dans un nombre quelconque de dimensions - j'arrive de le faire en 3D) et de l'organiser en k clusters avec des frontières distinctes. C'est bien ici parce que vous savez exactement combien d'orteils de ces canidés (dû).

En outre, il est mis en œuvre dans Scipy qui est vraiment sympa (http://docs.scipy.org/doc/scipy/reference/cluster.vq.html).

Voici un exemple de ce qu'il peut faire pour résoudre spatialement 3D clusters:

Ce que vous voulez faire est un peu différent (2D et comprend les valeurs de la pression), mais je pense toujours que vous pourriez donner un coup de feu.

Heres une autre approche que j'ai utilisé quand on fait quelque chose de similaire pour un grand télescope:

1) Recherche de la plus haute pixel.
Une fois que vous avez cela, la recherche autour pour le meilleur ajustement pour les 2x2 (peut-être la maximisation de la 2x2 somme), ou faire un 2d fit gaussien à l'intérieur de la sous-région de dire 4x4 centré sur la plus haute pixel.

Puis la mettre de 2x2 pixels que vous avez trouvé à zéro (ou peut-être 3x3) autour du centre de pointe

revenir à 1) et répétez jusqu'à ce que le plus haut sommet tombe au-dessous d'un seuil de bruit, ou vous disposez de tous les orteils vous avez besoin

une ébauche...

vous voudrez probablement utiliser des composants de l'algorithme d'isoler chaque patte de la région. wiki a une bonne description de ce (avec un peu de code) ici: http://en.wikipedia.org/wiki/Connected_Component_Labeling

vous aurez à prendre une décision sur l'opportunité d'utiliser 4 ou 8 connexité. personnellement, pour la plupart des problèmes, je préfère la 6-connexité. de toute façon, une fois que vous avez séparé chaque "empreinte de patte" comme un connectés à la région, il devrait être assez facile à parcourir la région et de trouver les maxima. une fois que vous avez trouvé la maxima, vous pouvez agrandir de manière itérative la région jusqu'à ce que vous atteignez un seuil prédéterminé, afin de l'identifier comme un "orteil".

un subtil problème, c'est que dès que vous commencez à l'aide de techniques de vision par ordinateur pour identifier quelque chose comme un droit/gauche/avant/à l'arrière de la patte et de commencer à examiner les différents orteils, vous devez commencer à prendre la rotation, l'inclinaison, et les traductions en compte. ceci est accompli grâce à l'analyse de soi-disant "moments". il y a quelques instants à considérer dans les applications de vision:

centrale moments: traduction invariant
normalisé moments: mise à l'échelle et à la traduction de l'invariant de
hu moments: la traduction, l'échelle et la rotation de l'invariant de

plus d'informations à propos de moments peuvent être trouvés par la recherche de l'image "moments" sur le wiki.

Il semble que vous pouvez tricher un peu en utilisant jetxee de l'algorithme. Il est de trouver les trois premiers orteils fine, et vous devriez être en mesure de deviner d'où le quatrième est en fonction de ça.

Bien, voici quelques-uns simple et pas très efficace du code, mais pour cette taille d'un ensemble de données, il est très bien.

import numpy as np
grid = np.array([[0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0.4,0.4,0.4,0,0,0],
[0,0,0,0,0.4,1.4,1.4,1.8,0.7,0,0,0,0,0],
[0,0,0,0,0.4,1.4,4,5.4,2.2,0.4,0,0,0,0],
[0,0,0.7,1.1,0.4,1.1,3.2,3.6,1.1,0,0,0,0,0],
[0,0.4,2.9,3.6,1.1,0.4,0.7,0.7,0.4,0.4,0,0,0,0],
[0,0.4,2.5,3.2,1.8,0.7,0.4,0.4,0.4,1.4,0.7,0,0,0],
[0,0,0.7,3.6,5.8,2.9,1.4,2.2,1.4,1.8,1.1,0,0,0],
[0,0,1.1,5,6.8,3.2,4,6.1,1.8,0.4,0.4,0,0,0],
[0,0,0.4,1.1,1.8,1.8,4.3,3.2,0.7,0,0,0,0,0],
[0,0,0,0,0,0.4,0.7,0.4,0,0,0,0,0,0]])
arr = []
for i in xrange(grid.shape[0] - 1):
for j in xrange(grid.shape[1] - 1):
tot = grid[i][j] + grid[i+1][j] + grid[i][j+1] + grid[i+1][j+1]
arr.append([(i,j),tot])
best = []
arr.sort(key = lambda x: x[1])
for i in xrange(5):
best.append(arr.pop())
badpos = set([(best[-1][0][0]+x,best[-1][0][1]+y)
for x in [-1,0,1] for y in [-1,0,1] if x != 0 or y != 0])
for j in xrange(len(arr)-1,-1,-1):
if arr[j][0] in badpos:
arr.pop(j)
for item in best:
print grid[item[0][0]:item[0][0]+2,item[0][1]:item[0][1]+2]

En gros, je viens de faire un tableau avec la position du coin supérieur gauche et la somme de chaque carré 2x2 et de les trier en fonction de la somme. Je puis prendre la 2x2 carrés avec la plus grande somme de la course, le mettre dans le best tableau, et supprimer tous les autres carrés de 2 x 2 qui a utilisé une partie de cette juste enlevé le carré 2x2.

Il semble bien fonctionner sauf avec la dernière patte (celui avec la plus petite somme à l'extrême droite dans votre première image), il s'avère qu'il y a deux autres admissible 2x2 places avec une somme plus importante (et ils ont une somme égale les unes des autres). L'un d'eux est encore en sélectionne un carré à partir de votre carré 2x2, mais l'autre est sur la gauche. Heureusement, par chance, nous voyons à être le choix de plus de ce que vous voulez, mais cela peut demander quelques autres idées pour être utilisé pour obtenir ce que vous voulez vraiment tout le temps.

Peut-être une approche naïve est suffisant ici: créer une liste de tous 2x2 places sur votre avion, à l'ordre par leur somme (en ordre décroissant).

Tout d'abord, sélectionnez la plus forte valeur carré dans votre "patte liste". Puis, de façon itérative, de choisir les 4 de la meilleure des places qui ne se coupent avec un de la déjà trouvé des places.

Je ne suis pas sûr que cela répond à la question, mais il semble que vous pouvez simplement regarder pour les n sommets les plus élevés qui n'ont pas de voisins.

Voici l'essentiel. Noter que c'est en Ruby, mais l'idée doit être clair.

require 'pp'
NUM_PEAKS = 5
NEIGHBOR_DISTANCE = 1
data = [[1,2,3,4,5],
[2,6,4,4,6],
[3,6,7,4,3],
]
def tuples(matrix)
tuples = []
matrix.each_with_index { |row, ri|
row.each_with_index { |value, ci|
tuples << [value, ri, ci]
}
}
tuples
end
def neighbor?(t1, t2, distance = 1)
[1,2].each { |axis|
return false if (t1[axis] - t2[axis]).abs > distance
}
true
end
# convert the matrix into a sorted list of tuples (value, row, col), highest peaks first
sorted = tuples(data).sort_by { |tuple| tuple.first }.reverse
# the list of peaks that don't have neighbors
non_neighboring_peaks = []
sorted.each { |candidate|
# always take the highest peak
if non_neighboring_peaks.empty?
non_neighboring_peaks << candidate
puts "took the first peak: #{candidate}"
else
# check that this candidate doesn't have any accepted neighbors
is_ok = true
non_neighboring_peaks.each { |accepted|
if neighbor?(candidate, accepted, NEIGHBOR_DISTANCE)
is_ok = false
break
end
}
if is_ok
non_neighboring_peaks << candidate
puts "took #{candidate}"
else
puts "denied #{candidate}"
end
end
}
pp non_neighboring_peaks