Wykrywanie pików w układzie 2D

Pomagam klinice weterynaryjnej mierzącej ciśnienie pod łapą psa. Używam Pythona do analizy danych i teraz utknąłem próbując podzielić łapy na (anatomiczne) podregiony.

Zrobiłem tablicę 2D każdej łapy, która składa się z maksymalnych wartości dla każdego czujnika, który został obciążony przez łapę w czasie. Oto przykład jednej łapy, w której użyłem programu Excel do narysowania obszarów, które chcę „wykryć”. Są to 2 na 2 pola wokół czujnika z lokalnymi maksimami, które razem mają największą sumę.

Więc spróbowałem trochę eksperymentować i postanowiłem po prostu poszukać maksimów każdej kolumny i rzędu (nie mogę patrzeć w jednym kierunku z powodu kształtu łapy). Wydaje się, że dość dobrze „wykrywa” to położenie oddzielnych palców, ale oznacza także sąsiednie czujniki.

Więc jaki byłby najlepszy sposób, aby powiedzieć Pythonowi, które z tych maksimów są tymi, których chcę?

Uwaga: kwadraty 2x2 nie mogą się pokrywać, ponieważ muszą być oddzielnymi palcami!

Również wziąłem 2x2 jako wygodę, każde bardziej zaawansowane rozwiązanie jest mile widziane, ale jestem po prostu naukowcem zajmującym się ludzkim ruchem, więc nie jestem ani prawdziwym programistą, ani matematykiem, więc proszę, zachowaj prostotę.

Oto wersja, którą można załadowaćnp.loadtxt

Wyniki

Wypróbowałem więc rozwiązanie @ jextee (zobacz wyniki poniżej). Jak widać, działa bardzo na przednich łapach, ale działa gorzej na tylne nogi.

Mówiąc dokładniej, nie rozpoznaje małego piku, który jest czwartym palcem. Jest to oczywiście nieodłączne od faktu, że pętla wygląda od góry w dół w kierunku najniższej wartości, bez względu na to, gdzie to jest.

Czy ktokolwiek wiedziałby, jak dostosować algorytm @ jextee, aby mógł znaleźć również czwarty palec u nogi?

Ponieważ nie przetworzyłem jeszcze żadnych prób, nie mogę dostarczyć żadnych innych próbek. Ale dane, które podałem wcześniej, były średnimi dla każdej łapy. Ten plik to tablica z maksymalnymi danymi 9 łap w kolejności, w której zetknęły się z płytką.

Ten obraz pokazuje, jak zostały one rozmieszczone przestrzennie na talerzu.

Aktualizacja:

Założyłem blog dla wszystkich zainteresowanych i skonfigurowałem SkyDrive ze wszystkimi nieprzetworzonymi pomiarami. Tak więc każdemu, kto prosi o więcej danych: więcej mocy dla Ciebie!

Nowa aktualizacja:

Więc po otrzymaniu pomocy z pytaniami dotyczącymi wykrywania i sortowania łap , w końcu mogłem sprawdzić wykrywanie palców u każdej łapy! Okazuje się, że nie działa tak dobrze w niczym innym, jak łapy wielkości podobne do tej z mojego przykładu. Oczywiście z perspektywy czasu to moja wina, że ​​wybrałem 2x2 tak arbitralnie.

Oto dobry przykład, w którym idzie źle: paznokieć jest rozpoznawany jako palec u nogi, a „pięta” jest tak szeroka, że ​​zostaje rozpoznana dwukrotnie!

Łapa jest zbyt duża, więc przyjęcie rozmiaru 2x2 bez zachodzenia na siebie powoduje dwukrotne wykrycie niektórych palców. Odwrotnie, u małych psów często nie znajduje się 5. palec u nogi, co, jak podejrzewam, jest spowodowane zbyt dużym obszarem 2x2.

Po wypróbowaniu obecnego rozwiązania we wszystkich moich pomiarach doszedłem do oszałamiającego wniosku, że dla prawie wszystkich moich małych psów nie znalazł piątego palca u nogi i że w ponad 50% uderzeń dużych psów znalazłby więcej!

Tak wyraźnie muszę to zmienić. Domyślam się, że zmieniłem rozmiar neighborhoodna coś mniejszego dla małych psów i większego dla dużych psów. Ale generate_binary_structurenie pozwoliłbym zmienić rozmiaru tablicy.

Dlatego mam nadzieję, że ktokolwiek ma lepszą propozycję lokalizacji palców, być może mając skalę obszaru palców wraz z rozmiarem łapy?

Wykryłem piki za pomocą lokalnego filtru maksymalnego . Oto wynik pierwszego zestawu danych 4 łap:

Uruchomiłem go również na drugim zestawie danych 9 łap i również zadziałało .

Oto jak to zrobić:

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

Wszystko, co musisz zrobić po, to użyć scipy.ndimage.measurements.labelmaski do oznaczenia wszystkich różnych obiektów. Wtedy będziesz mógł grać z nimi indywidualnie.

Należy pamiętać, że metoda działa dobrze, ponieważ tło nie jest głośne. Gdyby tak było, wykryłbyś w tle kilka innych niepożądanych szczytów. Innym ważnym czynnikiem jest wielkość dzielnicy . Będziesz musiał go dostosować, jeśli zmieni się rozmiar piku (powinien pozostać w przybliżeniu proporcjonalny).

Rozwiązanie

Plik danych: paw.txt . Kod źródłowy:

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"

Wyjście bez nakładających się kwadratów. Wygląda na to, że wybrane są te same obszary, co w twoim przykładzie.

Niektóre komentarze

Trudną częścią jest obliczenie sumy wszystkich kwadratów 2x2. Zakładałem, że potrzebujesz ich wszystkich, więc mogą się nakładać. Użyłem wycinków, aby wyciąć pierwsze / ostatnie kolumny i wiersze z oryginalnej tablicy 2D, a następnie nałożyć je wszystkie na siebie i obliczyć sumy.

Aby to lepiej zrozumieć, zobrazuj tablicę 3x3:

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

Następnie możesz wziąć jego plasterki:

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

Teraz wyobraź sobie, że układasz je jeden na drugim i sumujesz elementy w tych samych pozycjach. Te sumy będą dokładnie takie same sumy na kwadratach 2x2 z lewym górnym rogiem w tej samej pozycji:

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

Kiedy masz sumy powyżej 2x2 kwadratów, możesz użyć, maxaby znaleźć maksimum lub sort, lub sortedznaleźć szczyty.

Aby zapamiętać pozycje pików, łączę każdą wartość (sumę) z jej porządkową pozycją w spłaszczonej tablicy (patrz zip). Następnie ponownie drukuję pozycję wiersza / kolumny podczas drukowania wyników.

Notatki

Pozwoliłem, aby kwadraty 2x2 zachodziły na siebie. Edytowana wersja odfiltrowuje niektóre z nich, tak że w wynikach pojawiają się tylko nie nakładające się kwadraty.

Wybór palców (pomysł)

Kolejnym problemem jest to, jak wybrać palce, które prawdopodobnie będą palcami ze wszystkich szczytów. Mam pomysł, który może, ale nie musi, działać. Nie mam teraz czasu na jego wdrożenie, więc po prostu pseudo-kod.

Zauważyłem, że jeśli przednie palce pozostaną na prawie idealnym okręgu, tylny palec powinien znajdować się wewnątrz tego koła. Ponadto przednie palce są mniej więcej równomiernie rozmieszczone. Możemy próbować wykorzystać te właściwości heurystyczne do wykrycia palców.

Pseudo kod:

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

To podejście oparte na brutalnej sile. Jeśli N jest stosunkowo mały, to myślę, że jest to wykonalne. Dla N = 12 istnieją kombinacje C_12 ^ 5 = 792 razy 5 razy, aby wybrać tylny palec, więc 3960 przypadków do oceny dla każdej łapy.

Jest to problem z rejestracją obrazu . Ogólna strategia to:

• Miej znany przykład lub jakiś rodzaj wcześniejszych danych.
• Dopasuj dane do przykładu lub dopasuj przykład do danych.
• Pomaga, jeśli Twoje dane są z grubsza wyrównane.

Oto szorstkie i gotowe podejście , „najgłupsza rzecz, która mogłaby ewentualnie zadziałać”:

• Zacznij od pięciu współrzędnych palca mniej więcej w oczekiwanym miejscu.
• Z każdym z nich iteracyjnie wspinaj się na szczyt wzgórza. tzn. biorąc pod uwagę bieżącą pozycję, przejdź do maksymalnego sąsiedniego piksela, jeśli jego wartość jest większa niż bieżący piksel. Zatrzymaj się, gdy przestaną się poruszać Twoje współrzędne palca.

Aby przeciwdziałać problemowi z orientacją, możesz mieć około 8 ustawień początkowych dla podstawowych kierunków (północ, północ wschód itp.). Uruchom każdy z osobna i wyrzuć wszystkie wyniki, w których dwa lub więcej palców kończy się na tym samym pikselu. Zastanowię się nad tym jeszcze trochę, ale tego rodzaju rzeczy są wciąż badane w przetwarzaniu obrazu - nie ma prawidłowych odpowiedzi!

Nieco bardziej złożony pomysł: (ważony) K-oznacza grupowanie. Nie jest tak źle.

• Zacznij od pięciu współrzędnych palców, ale teraz są to „centra skupień”.

Następnie iteruj aż do konwergencji:

• Przypisz każdy piksel do najbliższego klastra (po prostu utwórz listę dla każdego klastra).
• Oblicz środek masy każdej gromady. Dla każdego klastra jest to: Suma (współrzędna * wartość intensywności) / Suma (współrzędna)
• Przenieś każdą gromadę do nowego środka masy.

Ta metoda prawie na pewno da znacznie lepsze wyniki, a otrzymasz masę każdego skupienia, co może pomóc w identyfikacji palców.

(Ponownie podałeś liczbę klastrów z góry. W przypadku klastrowania musisz określić gęstość w ten czy inny sposób: albo wybierz liczbę klastrów, odpowiednią w tym przypadku, albo wybierz promień klastra i zobacz, ile zakończymy w górę. Przykładem tego drugiego jest zmiana średniej ).

Przepraszamy za brak szczegółów implementacyjnych lub innych szczegółów. Kodowałbym to, ale mam termin. Jeśli nic nie zadziała do następnego tygodnia, daj mi znać, a dam ci szansę.

Używając trwałej homologii do analizy zestawu danych, otrzymuję następujący wynik (kliknij, aby powiększyć):

Jest to wersja 2D metody wykrywania pików opisanej w tej odpowiedzi SO . Powyższy rysunek pokazuje po prostu 0-wymiarowe klasy homologii trwałej posortowane według trwałości.

Zrobiłem przeskalować oryginalny zestaw danych o współczynnik 2 używając scipy.misc.imresize (). Zauważ jednak, że cztery łapy uznałem za jeden zestaw danych; podzielenie go na cztery ułatwiłoby problem.

Metodologia. Idea tego dość prosta: rozważ wykres funkcji, która przypisuje poziom każdemu pikselowi. To wygląda tak:

Teraz rozważ poziom wody na wysokości 255, który stale schodzi do niższych poziomów. Na lokalnych wyspach maksymalnych pojawiają się (narodziny). W punktach siodłowych łączą się dwie wyspy; uważamy dolną wyspę za połączoną z wyższą wyspą (śmierć). Tak zwany diagram trwałości (klas homologii 0-wymiarowej, nasze wyspy) przedstawia wartości śmiertelne z powodu narodzin wszystkich wysp:

Trwałość z wysp jest to różnica pomiędzy narodzenie- i śmierci poziomie; pionowa odległość kropki od szarej głównej przekątnej. Liczba ta określa wyspy, zmniejszając ich wytrzymałość.

Pierwsze zdjęcie pokazuje miejsca narodzin wysp. Ta metoda nie tylko podaje lokalne maksima, ale także określa ich „znaczenie” na podstawie wyżej wspomnianej trwałości. Następnie odfiltrowano by wszystkie wyspy o zbyt niskiej trwałości. Jednak w twoim przykładzie każda wyspa (tj. Każde lokalne maksimum) jest szczytem, ​​którego szukasz.

Kod Pythona można znaleźć tutaj .

Problem ten został dogłębnie zbadany przez fizyków. Istnieje dobra implementacja w ROOT . Spójrz na klasy TSpectrum (szczególnie TSpectrum2 dla twojego przypadku) i dokumentację dla nich.

Bibliografia:

1. M. Morhac i wsp .: Metody eliminacji tła dla wielowymiarowych widm koincydencji gamma. Nuclear Instruments and Methods in Physics Research A 401 (1997) 113-132.
2. M. Morhac i in .: Skuteczna jedno- i dwuwymiarowa dekonwolucja złota i jej zastosowanie do rozkładu widm promieniowania gamma. Nuclear Instruments and Methods in Physics Research A 401 (1997) 385-408.
3. M. Morhac i wsp .: Identyfikacja pików w widmach gamma promieniowania wielowymiarowego. Nuclear Instruments and Methods in Research Physics A 443 (2000), 108-125.

... a dla tych, którzy nie mają dostępu do subskrypcji NIM:

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

Oto pomysł: obliczasz (dyskretny) Laplaciana obrazu. Spodziewałbym się, że będzie on (ujemny i) duży przy maksimach, w sposób bardziej dramatyczny niż na oryginalnych obrazach. Tak więc maksima mogłyby być łatwiejsze do znalezienia.

Oto kolejny pomysł: jeśli znasz typowy rozmiar plam wysokociśnieniowych, możesz najpierw wygładzić obraz, zwijając go gaussianem tego samego rozmiaru. Może to ułatwić przetwarzanie obrazów.

Kilka pomysłów z mojej głowy:

• weź gradient (pochodną) skanu, sprawdź, czy to eliminuje fałszywe wywołania
• weź maksimum lokalnych maksimów

Możesz także rzucić okiem na OpenCV , ma dość przyzwoity interfejs API Pythona i może mieć pewne funkcje, które mogą ci się przydać.

Jestem pewien, że masz już wystarczająco dużo czasu, ale nie mogę nie poradzić, ale sugeruję użycie metody klastrowania k-średnich. K-średnich jest nie nadzorowanym algorytmem grupowania, który zabierze cię danych (w dowolnej liczbie wymiarów - robię to w 3D) i uporządkuje je w klastry z wyraźnymi granicami. Jest tu miło, ponieważ dokładnie wiesz, ile palców mają te kły (powinny) mieć.

Dodatkowo jest zaimplementowany w Scipy, co jest naprawdę miłe ( http://docs.scipy.org/doc/scipy/reference/cluster.vq.html ).

Oto przykład tego, co może zrobić, aby przestrzennie rozwiązać klastry 3D:

To, co chcesz zrobić, jest nieco inne (2D i zawiera wartości ciśnienia), ale nadal myślę, że możesz spróbować.

dzięki za surowe dane. Jestem w pociągu i to jest tak daleko, jak to możliwe (zbliża się mój przystanek). Zmasowałem twój plik txt za pomocą regexps i umieściłem go na stronie HTML z javascript do wizualizacji. Udostępniam go tutaj, ponieważ niektórzy, podobnie jak ja, mogą łatwiej go zhakować niż python.

Myślę, że dobrym podejściem będzie niezmienność skali i rotacji, a moim następnym krokiem będzie zbadanie mieszanin gaussów. (każda podkładka łapy jest środkiem gaussa).

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

Rozwiązanie fizyka:
Zdefiniuj 5 znaczników łap zidentyfikowanych według ich pozycji X_ii zainicjuj je losowymi pozycjami. Zdefiniuj jakąś funkcję energetyczną, łącząc pewną nagrodę za lokalizację znaczników w pozycjach łap z pewną karą za nakładanie się znaczników; powiedzmy:

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

( S(X_i)to średnia siła w kwadracie 2x2 wokół X_i, alfajest parametrem, który należy eksperymentalnie zwiększyć)

Teraz czas na magię Metropolis-Hastings:
1. Wybierz losowy znacznik i przesuń go o jeden piksel w losowym kierunku.
2. Oblicz dE, różnicę energii spowodowaną tym ruchem.
3. Uzyskaj jednolity losowy numer od 0-1 i nazwij go r.
4. Jeśli dE<0lub exp(-beta*dE)>r, zaakceptuj ruch i przejdź do 1; jeśli nie, cofnij ruch i przejdź do 1. Powtórz
tę czynność, aż znaczniki zbiegną się w łapy. Beta kontroluje skanowanie w celu optymalizacji kompromisu, dlatego należy go również zoptymalizować eksperymentalnie; można go również stale zwiększać wraz z czasem symulacji (symulowane wyżarzanie).

Oto inne podejście, które zastosowałem, robiąc coś podobnego dla dużego teleskopu:

1) Wyszukaj najwyższy piksel. Gdy już to zrobisz, poszukaj tego, aby uzyskać najlepsze dopasowanie dla 2x2 (być może maksymalizując sumę 2x2), lub wykonaj dopasowanie 2d gaussa w podregionie powiedzmy 4x4 wyśrodkowanym na najwyższym pikselu.

Następnie ustaw znalezione 2x2 piksele na zero (a może 3x3) wokół środka piku

wróć do 1) i powtarzaj, aż najwyższy szczyt spadnie poniżej progu hałasu lub będziesz mieć wszystkie palce, których potrzebujesz

Prawdopodobnie warto wypróbować sieci neuronowe, jeśli jesteś w stanie stworzyć dane treningowe ... ale to wymaga wielu próbek opatrzonych adnotacjami ręcznie.

przybliżony zarys ...

prawdopodobnie chcesz użyć algorytmu połączonych komponentów do izolacji każdego regionu łapy. wiki ma porządny opis tego (z pewnym kodem) tutaj: http://en.wikipedia.org/wiki/Connected_Component_Labeling

musisz podjąć decyzję, czy użyć połączenia 4 czy 8. osobiście, dla większości problemów wolę połączenie 6. tak czy inaczej, po oddzieleniu każdego „odcisku łapy” jako połączonego regionu, powinno być łatwo przejść przez ten region i znaleźć maksima. po znalezieniu maksimów możesz iteracyjnie powiększać region, aż osiągniesz z góry określony próg, aby zidentyfikować go jako dany „palec u nogi”.

jednym subtelnym problemem jest to, że jak tylko zaczniesz korzystać z komputerowych technik wizyjnych, aby zidentyfikować coś jako prawą / lewą / przednią / tylną łapę i zaczniesz patrzeć na poszczególne palce u stóp, musisz zacząć brać pod uwagę obroty, pochylenia i tłumaczenia. osiąga się to poprzez analizę tak zwanych „momentów”. w aplikacjach wizyjnych należy rozważyć kilka różnych momentów:

momenty centralne: niezmienne tłumaczenie znormalizowane momenty: skalowanie i tłumaczenie niezmienne momenty hu: niezmienne tłumaczenie, skala i rotacja

więcej informacji o chwilach można znaleźć, przeszukując „momenty obrazów” na wiki.

Być może możesz użyć czegoś takiego jak Gaussian Mixture Models. Oto pakiet Pythona do wykonywania GMM (właśnie przeprowadziłem wyszukiwanie w Google) http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/em/

Wygląda na to, że możesz trochę oszukiwać za pomocą algorytmu jetxee. Uważa, że ​​pierwsze trzy palce są w porządku i powinieneś być w stanie zgadnąć, gdzie opiera się czwarty.

Ciekawy problem. Rozwiązanie, które chciałbym wypróbować, jest następujące.

1. Zastosuj filtr dolnoprzepustowy, taki jak splot z maską gaussowską 2D. To da ci garść (prawdopodobnie, ale niekoniecznie zmiennoprzecinkowych) wartości.

2. Wykonaj nie-maksymalne tłumienie 2D, używając znanego przybliżonego promienia każdej poduszki (lub palca).

To powinno dać ci maksymalną liczbę pozycji bez posiadania wielu kandydatów, którzy są blisko siebie. Aby to wyjaśnić, promień maski w kroku 1 powinien być również podobny do promienia zastosowanego w kroku 2. Promień ten można wybrać, lub weterynarz może go wyraźnie zmierzyć wcześniej (zmienia się w zależności od wieku / rasy / itp.).

Niektóre z sugerowanych rozwiązań (średnia zmiana, sieci neuronowe itp.) Prawdopodobnie będą do pewnego stopnia działały, ale są zbyt skomplikowane i prawdopodobnie nie są idealne.

Oto prosty i niezbyt wydajny kod, ale dla tego rozmiaru zestawu danych jest w porządku.

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]

Zasadniczo po prostu tworzę tablicę z pozycją lewego górnego rogu i sumą każdego kwadratu 2x2 i sortuję według sumy. Następnie biorę kwadrat 2x2 z najwyższą sumą ze sporu, umieszczam go w besttablicy i usuwam wszystkie pozostałe kwadraty 2x2, które wykorzystały dowolną część tego właśnie usuniętego kwadratu 2x2.

Wygląda na to, że działa dobrze, z wyjątkiem ostatniej łapy (tej z najmniejszą sumą po prawej stronie na pierwszym zdjęciu), okazuje się, że istnieją dwa inne kwalifikujące się kwadraty 2x2 z większą sumą (i mają one taką samą kwotę jak wzajemnie). Jeden z nich nadal wybiera jeden kwadrat z twojego kwadratu 2x2, ale drugi znajduje się po lewej stronie. Na szczęście na szczęście widzimy, że wybieramy więcej tego, czego chcesz, ale może to wymagać innych pomysłów, aby uzyskać to, czego naprawdę chcesz przez cały czas.

Chcę wam tylko powiedzieć, że istnieje fajna opcja znalezienia lokalnego maximaobrazu na Pythonie:

from skimage.feature import peak_local_max

lub do 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

Może wystarczające jest naiwne podejście: Zbuduj listę wszystkich kwadratów 2x2 na swoim samolocie, uporządkuj je według ich sumy (w kolejności malejącej).

Najpierw wybierz kwadrat o najwyższej wartości z „listy łap”. Następnie, iteracyjnie wybierz 4 z następnych najlepszych kwadratów, które nie przecinają się z żadnym z wcześniej znalezionych kwadratów.

Istnieje wiele obszernych programów dostępnych w społeczności astronomii i kosmologii - jest to znaczący obszar badań zarówno historycznych, jak i obecnie.

Nie przejmuj się, jeśli nie jesteś astronomem - niektóre są łatwe w użyciu poza polem. Na przykład możesz użyć astropy / photutils:

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

[Powtarzanie krótkiego przykładowego kodu tutaj jest trochę niegrzeczne.]

Niekompletna i nieco stronnicza lista technik / pakietów / linków, które mogą być interesujące, znajduje się poniżej - dodaj więcej w komentarzach, a ja w razie potrzeby zaktualizuję tę odpowiedź. Oczywiście istnieje kompromis między dokładnością a zasobami obliczeniowymi. [Szczerze mówiąc, jest zbyt wielu, aby podać przykłady kodu w jednej odpowiedzi, takiej jak ta, więc nie jestem pewien, czy ta odpowiedź będzie latać, czy nie.]

Ekstraktor źródłowy https://www.astromatic.net/software/sextractor

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

Wyzwanie w poszukiwaniu źródła ASKAP / EMU: https://arxiv.org/abs/1509.03931

Możesz także wyszukać wyzwania związane z wydobyciem źródła przez Planck i / lub WMAP.

...

Co się stanie, jeśli wykonasz krok po kroku: najpierw zlokalizujesz globalne maksimum, w razie potrzeby przetworzysz otaczające punkty na podstawie ich wartości, następnie ustaw znaleziony region na zero i powtórz dla następnego.

Nie jestem pewien, czy to odpowiada na pytanie, ale wygląda na to, że możesz po prostu wyszukać n najwyższych szczytów, które nie mają sąsiadów.

Oto sedno. Zauważ, że jest w Ruby, ale pomysł powinien być jasny.

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