Détection des pics dans un réseau 2D

J'aide une clinique vétérinaire à mesurer la pression sous la patte d'un chien. J'utilise Python pour l'analyse de mes données et maintenant je suis coincé à essayer de diviser les pattes en sous-régions (anatomiques).

J'ai fait un tableau 2D de chaque patte, qui se compose des valeurs maximales pour chaque capteur qui a été chargé par la patte au fil du temps. Voici un exemple d'une patte, où j'ai utilisé Excel pour dessiner les zones que je veux «détecter». Ce sont des cases 2 par 2 autour du capteur avec des maxima locaux, qui ont ensemble la plus grande somme.

J'ai donc essayé quelques expériences et j'ai décidé de simplement chercher les maximums de chaque colonne et ligne (je ne peux pas regarder dans une direction à cause de la forme de la patte). Cela semble assez bien «détecter» l'emplacement des orteils séparés, mais cela marque également les capteurs voisins.

Alors, quelle serait la meilleure façon de dire à Python quels sont ces maximums que je veux?

Remarque: Les carrés 2x2 ne peuvent pas se chevaucher, car ils doivent être des orteils séparés!

J'ai également pris 2x2 pour plus de commodité, toute solution plus avancée est la bienvenue, mais je suis simplement un scientifique du mouvement humain, donc je ne suis ni un vrai programmeur ni un mathématicien, alors s'il vous plaît restez simple.

Voici une version qui peut être chargée avecnp.loadtxt

Résultats

J'ai donc essayé la solution de @ jextee (voir les résultats ci-dessous). Comme vous pouvez le voir, cela fonctionne très bien sur les pattes avant, mais cela fonctionne moins bien pour les pattes arrière.

Plus précisément, il ne peut pas reconnaître le petit pic qui est le quatrième orteil. Ceci est évidemment inhérent au fait que la boucle regarde de haut en bas vers la valeur la plus basse, sans tenir compte de son emplacement.

Quelqu'un sait-il comment modifier l'algorithme de @ jextee, afin qu'il puisse également trouver le 4e orteil?

Comme je n'ai pas encore traité d'autres essais, je ne peux pas fournir d'autres échantillons. Mais les données que j'ai données auparavant étaient les moyennes de chaque patte. Ce fichier est un tableau avec les données maximales de 9 pattes dans l'ordre où elles ont pris contact avec la plaque.

Cette image montre comment ils ont été répartis dans l'espace sur la plaque.

Mise à jour:

J'ai créé un blog pour toute personne intéressée et j'ai configuré un SkyDrive avec toutes les mesures brutes. Donc à tous ceux qui demandent plus de données: plus de pouvoir pour vous!

Nouvelle mise à jour:

Donc, après l'aide que j'ai reçue avec mes questions concernant la détection et le tri des pattes , j'ai finalement pu vérifier la détection des orteils pour chaque patte! Il s'avère que cela ne fonctionne pas si bien que les pattes de la taille de celle de mon propre exemple. Bien sûr, avec le recul, c'est ma faute si j'ai choisi arbitrairement le 2x2.

Voici un bel exemple de problème: un ongle est reconnu comme un orteil et le «talon» est si large qu'il est reconnu deux fois!

La patte est trop grande, donc en prenant une taille 2x2 sans chevauchement, certains orteils sont détectés deux fois. Dans l'autre sens, chez les petits chiens, il ne parvient souvent pas à trouver un 5e orteil, ce qui, je le soupçonne, est dû à une zone 2x2 trop grande.

Après avoir essayé la solution actuelle sur toutes mes mesures, je suis arrivé à la conclusion stupéfiante que pour presque tous mes petits chiens, il n'a pas trouvé de 5e orteil et que dans plus de 50% des impacts pour les grands chiens, il en trouverait plus!

Donc, clairement, je dois le changer. Ma propre supposition changeait la taille de la neighborhoodpour quelque chose de plus petit pour les petits chiens et plus grand pour les gros chiens. Mais generate_binary_structurene me laisse pas changer la taille du tableau.

Par conséquent, j'espère que quelqu'un d'autre aura une meilleure suggestion pour localiser les orteils, peut-être avoir l'échelle de la zone des orteils avec la taille de la patte?

J'ai détecté les pics en utilisant un filtre maximum local . Voici le résultat sur votre premier jeu de données de 4 pattes:

Je l'ai également exécuté sur le deuxième ensemble de données de 9 pattes et cela a également fonctionné .

Voici comment procéder:

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

Il vous suffit ensuite d'utiliser scipy.ndimage.measurements.labelle masque pour étiqueter tous les objets distincts. Ensuite, vous pourrez jouer avec eux individuellement.

Notez que la méthode fonctionne bien car le fond n'est pas bruyant. Si c'était le cas, vous détecteriez un tas d'autres pics indésirables en arrière-plan. Un autre facteur important est la taille du quartier . Vous devrez l'ajuster si la taille du pic change (la devrait rester à peu près proportionnelle).

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 soient sélectionnées comme dans votre exemple.

Certains commentaires

La partie délicate consiste à calculer les sommes de tous les carrés 2x2. J'ai supposé que vous en aviez besoin, donc il pourrait y avoir des chevauchements. J'ai utilisé des tranches pour couper les premières / dernières colonnes et lignes du tableau 2D d'origine, puis les recouvrir toutes ensemble et calculer des sommes.

Pour mieux le comprendre, imaginez une matrice 3x3:

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

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

Imaginez maintenant que vous les empilez les uns au-dessus des autres et que vous additionnez les éléments aux mêmes positions. Ces sommes seront exactement les mêmes sur les carrés 2x2 avec le coin supérieur 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 les sommes sur 2x2 carrés, vous pouvez utiliser maxpour trouver le maximum, ou sort, ou sortedpour trouver les pics.

Pour me souvenir des positions des pics, je couple chaque valeur (la somme) avec sa position ordinale dans un tableau aplati (voir zip). Ensuite, je calcule à nouveau la position de la ligne / colonne lorsque j'imprime les résultats.

Remarques

J'ai permis aux carrés 2x2 de se chevaucher. La version modifiée en filtre certains, de sorte que seuls les carrés qui ne se chevauchent pas apparaissent dans les résultats.

Choisir les doigts (une idée)

Un autre problème est de savoir comment choisir ce qui est susceptible d'être des doigts parmi tous les pics. J'ai une idée qui peut ou non fonctionner. Je n'ai pas le temps de l'implémenter en ce moment, donc juste un pseudo-code.

J'ai remarqué que si les doigts avant restent sur un cercle presque parfait, le doigt arrière devrait être à l'intérieur de ce cercle. De plus, les doigts avant sont plus ou moins également espacés. Nous pouvons 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

Il s'agit d'une approche par force brute. Si N est relativement petit, je pense que c'est faisable. Pour N = 12, il existe C_12 ^ 5 = 792 combinaisons, fois 5 façons de sélectionner un doigt arrière, donc 3960 cas à évaluer pour chaque patte.

Il s'agit d'un problème d'enregistrement d'image . La stratégie générale est la suivante:

• Ayez un exemple connu, ou une sorte de préalable sur les données.
• Ajustez vos données à l'exemple ou ajustez l'exemple à vos données.
• Cela aide si vos données sont à peu près alignées en premier lieu.

Voici une approche approximative et prête , "la chose la plus stupide qui pourrait éventuellement fonctionner":

• Commencez avec cinq coordonnées d'orteil à peu près à l'endroit que vous attendez.
• Avec chacun, montez itérativement au sommet de la colline. c'est-à-dire une position actuelle donnée, déplacez-vous vers le pixel voisin maximum, si sa valeur est supérieure au pixel actuel. Arrêtez-vous lorsque vos coordonnées d'orteil ont cessé de bouger.

Pour contrer le problème d'orientation, vous pouvez avoir environ 8 réglages initiaux pour les directions de base (Nord, Nord-Est, etc.). Exécutez chacun individuellement et jetez tous les résultats où deux ou plusieurs orteils se retrouvent au même pixel. J'y penserai un peu plus, mais ce genre de chose est toujours à l'étude dans le traitement d'image - il n'y a pas de bonnes réponses!

Idée un peu plus complexe: (pondéré) clustering K-means. C'est pas si mal.

• Commencez avec cinq coordonnées d'orteil, mais maintenant ce sont des "centres de cluster".

Itérer ensuite jusqu'à convergence:

• Attribuez chaque pixel au cluster le plus proche (faites simplement une liste pour chaque cluster).
• Calculez le centre de masse de chaque groupe. Pour chaque cluster, c'est: Sum (coordonnée * valeur d'intensité) / Sum (coordonnée)
• Déplacez chaque cluster vers le nouveau centre de masse.

Cette méthode donnera presque certainement de bien meilleurs résultats et vous obtiendrez la masse de chaque grappe, ce qui peut aider à identifier les orteils.

(Encore une fois, vous avez spécifié le nombre de clusters à l'avance. Avec le clustering, vous devez spécifier la densité d'une manière ou d'une autre: choisissez le nombre de clusters, approprié dans ce cas, ou choisissez un rayon de cluster et voyez combien vous finissez Un exemple de ce dernier est le décalage moyen .)

Désolé pour le manque de détails d'implémentation ou d'autres détails. Je coderais cela, mais j'ai un délai. Si rien d'autre n'a fonctionné la semaine prochaine, faites le moi savoir et je vais essayer.

En utilisant l'homologie persistante pour analyser votre ensemble de données, j'obtiens le résultat suivant (cliquez pour agrandir):

Il s'agit de la version 2D de la méthode de détection des pics décrite dans cette réponse SO . La figure ci-dessus montre simplement les classes d'homologie persistante à 0 dimension triées par persistance.

J'ai mis à l'échelle l'ensemble de données d'origine par un facteur de 2 en utilisant scipy.misc.imresize (). Cependant, notez que j'ai considéré les quatre pattes comme un ensemble de données; le diviser en quatre faciliterait le problème.

Méthodologie. L'idée derrière cela est assez simple: considérons le graphique de fonction de la fonction qui attribue à chaque pixel son niveau. Cela ressemble à ceci:

Considérons maintenant un niveau d'eau à la hauteur 255 qui descend en continu vers des niveaux inférieurs. Aux îles maxima locales apparaissent (naissance). Aux points de selle, deux îles fusionnent; nous considérons l'île inférieure comme fusionnant avec l'île supérieure (mort). Le soi-disant diagramme de persistance (des classes d'homologie de dimension 0, nos îles) représente les valeurs de mortalité par rapport à la naissance de toutes les îles:

La persistance d'une île est alors la différence entre les niveaux de naissance et de mort; la distance verticale d'un point à la diagonale grise principale. La figure marque les îles en diminuant la persistance.

La toute première photo montre les lieux de naissance des îles. Cette méthode donne non seulement les maxima locaux mais quantifie également leur «signification» par la persistance mentionnée ci-dessus. On filtrerait alors toutes les îles avec une persistance trop faible. Cependant, dans votre exemple, chaque île (c'est-à-dire chaque maximum local) est un pic que vous recherchez.

Le code Python peut être trouvé ici .

Ce problème a été étudié en profondeur par les physiciens. Il y a une bonne implémentation dans ROOT . Regardez les classes TSpectrum (en particulier TSpectrum2 pour votre cas) et leur documentation.

Références:

1. M.Morhac et al.: Méthodes d'élimination de fond pour les spectres de rayons gamma de coïncidence multidimensionnelle. Instruments et méthodes nucléaires en recherche physique A 401 (1997) 113-132.
2. M.Morhac et al.: Déconvolution d'or unidimensionnelle et bidimensionnelle efficace et son application à la décomposition des spectres de rayons gamma. Instruments et méthodes nucléaires en recherche physique A 401 (1997) 385-408.
3. M.Morhac et al.: Identification des pics dans les spectres de rayons gamma à coïncidence multidimensionnelle. Instruments et méthodes nucléaires en physique de la recherche A 443 (2000), 108-125.

... et pour ceux qui n'ont pas accès à un abonnement à NIM:

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

Voici une idée: vous calculez le laplacien (discret) de l'image. Je m'attendrais à ce qu'il soit (négatif et) grand au maximum, d'une manière plus dramatique que dans les images originales. Ainsi, les maxima pourraient être plus faciles à trouver.

Voici une autre idée: si vous connaissez la taille typique des taches à haute pression, vous pouvez d'abord lisser votre image en la convoluant avec un gaussien de la même taille. Cela peut vous donner des images plus simples à traiter.

Juste quelques idées du haut de ma tête:

• prendre le gradient (dérivé) du scan, voir si cela élimine les faux appels
• prendre le maximum des maxima locaux

Vous voudrez peut-être également jeter un œil à OpenCV , il dispose d'une API Python assez décente et pourrait avoir certaines fonctions que vous trouveriez utiles.

Je suis sûr que vous en avez assez pour continuer, mais je ne peux m'empêcher de suggérer d'utiliser la méthode de clustering k-means. k-means est un algorithme de clustering non supervisé qui vous prendra des données (dans n'importe quel nombre de dimensions - il se trouve que je le fais en 3D) et les organisera en k clusters avec des limites distinctes. C'est bien ici parce que vous savez exactement combien d'orteils ces canines (devraient) avoir.

De plus, il est implémenté dans Scipy, ce 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 les clusters 3D:

Ce que vous voulez faire est un peu différent (2D et inclut des valeurs de pression), mais je pense toujours que vous pourriez essayer.

merci pour les données brutes. Je suis dans le train et c'est aussi loin que je suis arrivé (mon arrêt arrive). J'ai massé votre fichier txt avec des expressions rationnelles et l'ai déposé dans une page html avec du javascript pour la visualisation. Je le partage ici parce que certains, comme moi, pourraient le trouver plus facilement piratable que python.

Je pense qu'une bonne approche sera l'échelle et la rotation invariantes, et ma prochaine étape sera d'étudier les mélanges de gaussiens. (chaque patte étant le centre d'un gaussien).

    <html>
<head>
    <script type="text/javascript" src="http://vis.stanford.edu/protovis/protovis-r3.2.js"></script> 
    <script type="text/javascript">
    var heatmap = [[[0,0,0,0,0,0,0,4,4,0,0,0,0],
[0,0,0,0,0,7,14,22,18,7,0,0,0],
[0,0,0,0,11,40,65,43,18,7,0,0,0],
[0,0,0,0,14,61,72,32,7,4,11,14,4],
[0,7,14,11,7,22,25,11,4,14,65,72,14],
[4,29,79,54,14,7,4,11,18,29,79,83,18],
[0,18,54,32,18,43,36,29,61,76,25,18,4],
[0,4,7,7,25,90,79,36,79,90,22,0,0],
[0,0,0,0,11,47,40,14,29,36,7,0,0],
[0,0,0,0,4,7,7,4,4,4,0,0,0]
],[
[0,0,0,4,4,0,0,0,0,0,0,0,0],
[0,0,11,18,18,7,0,0,0,0,0,0,0],
[0,4,29,47,29,7,0,4,4,0,0,0,0],
[0,0,11,29,29,7,7,22,25,7,0,0,0],
[0,0,0,4,4,4,14,61,83,22,0,0,0],
[4,7,4,4,4,4,14,32,25,7,0,0,0],
[4,11,7,14,25,25,47,79,32,4,0,0,0],
[0,4,4,22,58,40,29,86,36,4,0,0,0],
[0,0,0,7,18,14,7,18,7,0,0,0,0],
[0,0,0,0,4,4,0,0,0,0,0,0,0],
],[
[0,0,0,4,11,11,7,4,0,0,0,0,0],
[0,0,0,4,22,36,32,22,11,4,0,0,0],
[4,11,7,4,11,29,54,50,22,4,0,0,0],
[11,58,43,11,4,11,25,22,11,11,18,7,0],
[11,50,43,18,11,4,4,7,18,61,86,29,4],
[0,11,18,54,58,25,32,50,32,47,54,14,0],
[0,0,14,72,76,40,86,101,32,11,7,4,0],
[0,0,4,22,22,18,47,65,18,0,0,0,0],
[0,0,0,0,4,4,7,11,4,0,0,0,0],
],[
[0,0,0,0,4,4,4,0,0,0,0,0,0],
[0,0,0,4,14,14,18,7,0,0,0,0,0],
[0,0,0,4,14,40,54,22,4,0,0,0,0],
[0,7,11,4,11,32,36,11,0,0,0,0,0],
[4,29,36,11,4,7,7,4,4,0,0,0,0],
[4,25,32,18,7,4,4,4,14,7,0,0,0],
[0,7,36,58,29,14,22,14,18,11,0,0,0],
[0,11,50,68,32,40,61,18,4,4,0,0,0],
[0,4,11,18,18,43,32,7,0,0,0,0,0],
[0,0,0,0,4,7,4,0,0,0,0,0,0],
],[
[0,0,0,0,0,0,4,7,4,0,0,0,0],
[0,0,0,0,4,18,25,32,25,7,0,0,0],
[0,0,0,4,18,65,68,29,11,0,0,0,0],
[0,4,4,4,18,65,54,18,4,7,14,11,0],
[4,22,36,14,4,14,11,7,7,29,79,47,7],
[7,54,76,36,18,14,11,36,40,32,72,36,4],
[4,11,18,18,61,79,36,54,97,40,14,7,0],
[0,0,0,11,58,101,40,47,108,50,7,0,0],
[0,0,0,4,11,25,7,11,22,11,0,0,0],
[0,0,0,0,0,4,0,0,0,0,0,0,0],
],[
[0,0,4,7,4,0,0,0,0,0,0,0,0],
[0,0,11,22,14,4,0,4,0,0,0,0,0],
[0,0,7,18,14,4,4,14,18,4,0,0,0],
[0,4,0,4,4,0,4,32,54,18,0,0,0],
[4,11,7,4,7,7,18,29,22,4,0,0,0],
[7,18,7,22,40,25,50,76,25,4,0,0,0],
[0,4,4,22,61,32,25,54,18,0,0,0,0],
[0,0,0,4,11,7,4,11,4,0,0,0,0],
],[
[0,0,0,0,7,14,11,4,0,0,0,0,0],
[0,0,0,4,18,43,50,32,14,4,0,0,0],
[0,4,11,4,7,29,61,65,43,11,0,0,0],
[4,18,54,25,7,11,32,40,25,7,11,4,0],
[4,36,86,40,11,7,7,7,7,25,58,25,4],
[0,7,18,25,65,40,18,25,22,22,47,18,0],
[0,0,4,32,79,47,43,86,54,11,7,4,0],
[0,0,0,14,32,14,25,61,40,7,0,0,0],
[0,0,0,0,4,4,4,11,7,0,0,0,0],
],[
[0,0,0,0,4,7,11,4,0,0,0,0,0],
[0,4,4,0,4,11,18,11,0,0,0,0,0],
[4,11,11,4,0,4,4,4,0,0,0,0,0],
[4,18,14,7,4,0,0,4,7,7,0,0,0],
[0,7,18,29,14,11,11,7,18,18,4,0,0],
[0,11,43,50,29,43,40,11,4,4,0,0,0],
[0,4,18,25,22,54,40,7,0,0,0,0,0],
[0,0,4,4,4,11,7,0,0,0,0,0,0],
],[
[0,0,0,0,0,7,7,7,7,0,0,0,0],
[0,0,0,0,7,32,32,18,4,0,0,0,0],
[0,0,0,0,11,54,40,14,4,4,22,11,0],
[0,7,14,11,4,14,11,4,4,25,94,50,7],
[4,25,65,43,11,7,4,7,22,25,54,36,7],
[0,7,25,22,29,58,32,25,72,61,14,7,0],
[0,0,4,4,40,115,68,29,83,72,11,0,0],
[0,0,0,0,11,29,18,7,18,14,4,0,0],
[0,0,0,0,0,4,0,0,0,0,0,0,0],
]
];
</script>
</head>
<body>
    <script type="text/javascript+protovis">    
    for (var a=0; a < heatmap.length; a++) {
    var w = heatmap[a][0].length,
    h = heatmap[a].length;
var vis = new pv.Panel()
    .width(w * 6)
    .height(h * 6)
    .strokeStyle("#aaa")
    .lineWidth(4)
    .antialias(true);
vis.add(pv.Image)
    .imageWidth(w)
    .imageHeight(h)
    .image(pv.Scale.linear()
        .domain(0, 99, 100)
        .range("#000", "#fff", '#ff0a0a')
        .by(function(i, j) heatmap[a][j][i]));
vis.render();
}
</script>
  </body>
</html>

Solution du physicien:
Définissez 5 marqueurs de patte identifiés par leur position X_iet initiez-les avec des positions aléatoires. Définir une fonction énergétique combinant une récompense pour l'emplacement des marqueurs dans la position des pattes avec une punition pour le chevauchement des marqueurs; Disons:

E(X_i;S)=-Sum_i(S(X_i))+alfa*Sum_ij (|X_i-Xj|<=2*sqrt(2)?1:0)

( S(X_i)est la force moyenne en carré de 2x2 environ X_i, alfaest un paramètre à atteindre expérimentalement)

Il est maintenant temps de faire de la magie Metropolis-Hastings:
1. Sélectionnez un marqueur aléatoire et déplacez-le d'un pixel dans une direction aléatoire.
2. Calculez dE, la différence d'énergie provoquée par ce mouvement.
3. Obtenez un nombre aléatoire uniforme de 0 à 1 et appelez-le r.
4. Si dE<0ou exp(-beta*dE)>r, acceptez le mouvement et passez à 1; sinon, annulez le mouvement et passez à 1.
Cela devrait être répété jusqu'à ce que les marqueurs convergent vers les pattes. La bêta contrôle l'analyse pour optimiser le compromis, elle doit donc également être optimisée expérimentalement; il peut également être constamment augmenté avec le temps de simulation (recuit simulé).

Voici une autre approche que j'ai utilisée en faisant quelque chose de similaire pour un grand télescope:

1) Recherchez le pixel le plus élevé. Une fois que vous avez cela, recherchez autour de lui le meilleur ajustement pour 2x2 (peut-être en maximisant la somme 2x2), ou faites un ajustement gaussien 2d à l'intérieur de la sous-région de 4x4 centrée sur le pixel le plus élevé.

Ensuite, définissez ces pixels 2x2 que vous avez trouvés à zéro (ou peut-être 3x3) autour du centre du pic

revenir à 1) et répéter jusqu'à ce que le pic le plus élevé tombe en dessous d'un seuil de bruit, ou que vous ayez tous les orteils dont vous avez besoin

Cela vaut probablement la peine d'essayer avec les réseaux de neurones si vous êtes capable de créer des données d'entraînement ... mais cela nécessite de nombreux échantillons annotés à la main.

un aperçu approximatif ...

vous voudrez probablement utiliser un algorithme de composants connectés pour isoler chaque région de patte. wiki a une description décente de cela (avec du code) ici: http://en.wikipedia.org/wiki/Connected_Component_Labeling

vous devrez décider si vous souhaitez utiliser la connectivité 4 ou 8. personnellement, pour la plupart des problèmes, je préfère la connectivité 6. de toute façon, une fois que vous avez séparé chaque "empreinte de patte" en tant que région connectée, il devrait être assez facile de parcourir la région et de trouver les maxima. une fois que vous avez trouvé les maxima, vous pouvez agrandir la région de manière itérative jusqu'à ce que vous atteigniez un seuil prédéterminé afin de l'identifier comme un "orteil" donné.

un problème subtil ici est que dès que vous commencez à utiliser des techniques de vision par ordinateur pour identifier quelque chose comme une patte droite / gauche / avant / arrière et que vous commencez à regarder les orteils individuels, vous devez commencer à prendre en compte les rotations, les asymétries et les traductions. ceci est accompli grâce à l'analyse des soi-disant «moments». il y a quelques moments différents à considérer dans les applications de vision:

moments centraux: invariant de translation moments normalisés: invariant de mise à l'échelle et de translation moments hu: invariant de translation, d'échelle et de rotation

plus d'informations sur les moments peuvent être trouvées en recherchant "moments image" sur wiki.

Vous pouvez peut-être utiliser quelque chose comme des modèles de mélange gaussiens. Voici un package Python pour faire des GMM (vient de faire une recherche Google) http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/em/

Il semble que vous puissiez tricher un peu en utilisant l'algorithme de jetxee. Il trouve que les trois premiers orteils vont bien, et vous devriez pouvoir deviner où le quatrième est basé.

Problème intéressant. La solution que j'essaierais est la suivante.

1. Appliquez un filtre passe-bas, tel que la convolution avec un masque gaussien 2D. Cela vous donnera un tas de valeurs (probablement, mais pas nécessairement en virgule flottante).

2. Effectuez une suppression 2D non maximale en utilisant le rayon approximatif connu de chaque patte (ou orteil).

Cela devrait vous donner les positions maximales sans avoir plusieurs candidats qui sont proches les uns des autres. Juste pour clarifier, le rayon du masque à l'étape 1 devrait également être similaire au rayon utilisé à l'étape 2. Ce rayon pourrait être sélectionnable, ou le vétérinaire pourrait le mesurer explicitement à l'avance (cela variera avec l'âge / la race / etc.).

Certaines des solutions suggérées (décalage moyen, réseaux neuronaux, etc.) fonctionneront probablement dans une certaine mesure, mais sont trop compliquées et probablement pas idéales.

Eh bien, voici un code simple et pas terriblement efficace, mais pour cette taille d'un ensemble de données, c'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]

Je fais simplement un tableau avec la position du coin supérieur gauche et la somme de chaque carré 2x2 et je le trie par la somme. Je prends ensuite le carré 2x2 avec la somme la plus élevée hors de contention, le place dans le besttableau et supprime tous les autres carrés 2x2 qui ont utilisé n'importe quelle partie de ce carré 2x2 juste supprimé.

Cela semble bien fonctionner sauf avec la dernière patte (celle avec la plus petite somme à l'extrême droite dans votre première photo), il s'avère qu'il y a deux autres carrés 2x2 éligibles avec une somme plus grande (et ils ont une somme égale à L'une et l'autre). L'un d'eux sélectionne toujours un carré de votre carré 2x2, mais l'autre est à gauche. Heureusement, par chance, nous pensons choisir plus de celui que vous voudriez, mais cela peut nécessiter d'autres idées pour obtenir tout le temps ce que vous voulez réellement.

Je veux juste vous dire qu'il y a une bonne option pour trouver du local maximadans des images avec python:

from skimage.feature import peak_local_max

ou pour skimage 0.8.0:

from skimage.feature.peak import peak_local_max

http://scikit-image.org/docs/0.8.0/api/skimage.feature.peak.html

Peut-être qu'une approche naïve est suffisante ici: Construisez une liste de tous les carrés 2x2 sur votre avion, triez-les par leur somme (en ordre décroissant).

Tout d'abord, sélectionnez le carré le plus apprécié dans votre "liste de pattes". Ensuite, choisissez de manière itérative 4 des meilleurs carrés suivants qui ne se croisent avec aucun des carrés trouvés précédemment.

Il existe plusieurs et vastes logiciels disponibles auprès de la communauté de l'astronomie et de la cosmologie - il s'agit d'un domaine de recherche important à la fois historique et actuel.

Ne vous inquiétez pas si vous n'êtes pas astronome - certains sont faciles à utiliser en dehors du terrain. Par exemple, vous pouvez utiliser l'astropie / photutils:

https://photutils.readthedocs.io/en/stable/detection.html#local-peak-detection

[Il semble un peu grossier de répéter leur court exemple de code ici.]

Une liste incomplète et légèrement biaisée des techniques / packages / liens qui pourraient être intéressants est donnée ci-dessous - ajoutez-en plus dans les commentaires et je mettrai à jour cette réponse si nécessaire. Bien sûr, il existe un compromis entre la précision et les ressources de calcul. [Honnêtement, il y en a trop pour donner des exemples de code dans une seule réponse comme celle-ci, donc je ne sais pas si cette réponse volera ou non.]

Extracteur de source https://www.astromatic.net/software/sextractor

MultiNest https://github.com/farhanferoz/MultiNest [+ pyMultiNest]

Défi de recherche de source ASKAP / EMU: https://arxiv.org/abs/1509.03931

Vous pouvez également rechercher des défis d'extraction de source Planck et / ou WMAP.

...

Que faire si vous procédez étape par étape: vous localisez d'abord le maximum global, traitez si nécessaire les points environnants en fonction de leur valeur, puis définissez la région trouvée sur zéro et répétez pour la suivante.

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

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

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