Pregunta Detección de picos en una matriz 2D


Estoy ayudando a una clínica veterinaria a medir la presión debajo de la pata de un perro. Utilizo Python para mi análisis de datos y ahora estoy atascado tratando de dividir las patas en subregiones (anatómicas).

Hice una matriz 2D de cada pata, que consiste en los valores máximos para cada sensor que ha sido cargado por la pata a lo largo del tiempo. Aquí hay un ejemplo de una pata, donde utilicé Excel para dibujar las áreas que quiero 'detectar'. Estos son 2 por 2 cajas alrededor del sensor con máximos locales, que juntos tienen la mayor suma.

alt text

Así que intenté experimentar y decidí simplemente buscar los máximos de cada columna y fila (no puedo mirar en una dirección debido a la forma de la pata). Esto parece 'detectar' bastante bien la ubicación de los dedos separados, pero también marca los sensores vecinos.

alt text

Entonces, ¿cuál sería la mejor manera de decirle a Python cuáles de estos máximos son los que quiero?

Nota: ¡Los cuadrados de 2x2 no se pueden superponer, ya que tienen que estar separados!

También tomé 2x2 como una conveniencia, cualquier solución más avanzada es bienvenida, pero simplemente soy un científico del movimiento humano, así que no soy un programador real o un matemático, así que por favor manténganlo 'simple'.

Aquí está un versión que se puede cargar con np.loadtxt


Resultados

Así que probé la solución de @jextee (ver los resultados a continuación). Como puede ver, funciona muy bien en las patas delanteras, pero funciona menos bien para las patas traseras.

Más específicamente, no puede reconocer el pequeño pico que es el cuarto dedo del pie. Esto es obviamente inherente al hecho de que el ciclo mira hacia arriba hacia el valor más bajo, sin tener en cuenta dónde está.

¿Alguien sabría cómo modificar el algoritmo de @ jextee para que también pueda encontrar el 4º?

alt text

Como todavía no he procesado ningún otro ensayo, no puedo suministrar ninguna otra muestra. Pero los datos que di antes eran los promedios de cada pata. Este archivo es una matriz con los datos máximos de 9 patas en el orden en que hicieron contacto con la placa.

Esta imagen muestra cómo se distribuyeron espacialmente sobre el plato.

alt text

Actualizar:

He creado un blog para cualquier persona interesada y Configuré un SkyDrive con todas las medidas sin procesar. Entonces, para cualquiera que solicite más datos: ¡más poder para usted!


Nueva actualización:

Entonces, después de la ayuda que recibí con mis preguntas sobre detección de pata y clasificación de la pata, ¡Finalmente pude verificar la detección de dedos para cada pata! Resulta que no funciona tan bien en nada más que patas del tamaño de la de mi propio ejemplo. Por supuesto, en retrospectiva, es mi culpa por elegir el 2x2 de forma tan arbitraria.

Aquí hay un buen ejemplo de dónde sale mal: un clavo se reconoce como un dedo del pie y el "talón" es tan ancho que se reconoce dos veces.

alt text

La pata es demasiado grande, por lo que tomar un tamaño de 2x2 sin superposición hace que algunos dedos se detecten dos veces. A la inversa, en perros pequeños a menudo no se encuentra el quinto dedo, que sospecho que es causado por el área 2x2 que es demasiado grande.

Después probando la solución actual en todas mis mediciones Llegué a la sorprendente conclusión de que para casi todos mis perros pequeños no encontró un quinto dedo y que en más del 50% de los impactos para los perros grandes, ¡encontraría más!

Entonces claramente necesito cambiarlo. Mi propia conjetura fue cambiar el tamaño de la neighborhood a algo más pequeño para perros pequeños y más grande para perros grandes. Pero generate_binary_structure no me dejaba cambiar el tamaño de la matriz.

Por lo tanto, espero que alguien más tenga una mejor sugerencia para ubicar los dedos de los pies, quizás teniendo la escala del área del dedo del pie con el tamaño de la pata.


735
2017-09-11 03:38


origen


Respuestas:


Detecté los picos usando un filtro máximo local. Aquí está el resultado en su primer conjunto de datos de 4 patas: Peaks detection result

También lo ejecuté en el segundo conjunto de datos de 9 patas y funcionó tan bien.

Así es como lo haces:

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

Todo lo que necesitas hacer es usar scipy.ndimage.measurements.label en la máscara para etiquetar todos los objetos distintos. Entonces podrás jugar con ellos individualmente.

Nota que el método funciona bien porque el fondo no es ruidoso. Si lo fuera, detectaría un montón de otros picos no deseados en el fondo. Otro factor importante es el tamaño del barrio. Tendrá que ajustarlo si el tamaño máximo cambia (el tamaño debería ser aproximadamente proporcional).


257
2017-09-10 14:09



Solución

Archivo de datos: paw.txt. Código fuente:

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"

Salida sin superponer cuadrados Parece que las mismas áreas se seleccionan como en su ejemplo.

Algunos comentarios

La parte difícil es calcular las sumas de todos los cuadrados de 2x2. Supuse que los necesita a todos, por lo que podría haber cierta superposición. Utilicé cortes para cortar las primeras / últimas columnas y filas de la matriz 2D original, y luego superponerlas todas juntas y calcular las sumas.

Para entenderlo mejor, imaginando una matriz de 3x3:

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

Entonces puedes tomar sus rebanadas:

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

Ahora imagine que los apila uno sobre el otro y suma elementos en las mismas posiciones. Estas sumas serán exactamente las mismas sumas en los cuadrados de 2x2 con la esquina superior izquierda en la misma posición:

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

Cuando tienes las sumas en cuadrados de 2x2, puedes usar max para encontrar el máximo, o sort, o sorted para encontrar los picos

Para recordar las posiciones de los picos, acoplo cada valor (la suma) con su posición ordinal en una matriz aplanada (ver zip) Luego calculo la posición de fila / columna nuevamente cuando imprimo los resultados.

Notas

Permití que los cuadrados 2x2 se superpongan. La versión editada filtra algunos de ellos de forma que solo aparezcan cuadrados no superpuestos en los resultados.

Elegir los dedos (una idea)

Otro problema es cómo elegir lo que es probable que sean dedos de todos los picos. Tengo una idea que puede o no funcionar. No tengo tiempo para implementarlo ahora mismo, así que solo pseudocódigo.

Noté que si los dedos frontales permanecen en un círculo casi perfecto, el dedo posterior debería estar dentro de ese círculo. Además, los dedos frontales están espaciados más o menos equitativamente. Podemos intentar usar estas propiedades heurísticas para detectar los dedos.

Pseudo código:

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

Este es un enfoque de fuerza bruta. Si N es relativamente pequeño, entonces creo que es factible. Para N = 12, hay combinaciones C_12 ^ 5 = 792, multiplicadas por 5 formas de seleccionar un dedo posterior, por lo que se deben evaluar 3960 casos por cada pata.


40
2017-09-10 14:54



Esto es un problema de registro de imagen. La estrategia general es:

  • Tener un ejemplo conocido, o algún tipo de anterior en los datos.
  • Ajuste sus datos al ejemplo o ajuste el ejemplo a sus datos.
  • Ayuda si tus datos son aproximadamente alineado en primer lugar.

Aquí hay un acercamiento rudo y listo, "lo más estúpido que podría funcionar":

  • Comience con cinco coordenadas del dedo del pie en aproximadamente el lugar que espera.
  • Con cada uno, sube iterativamente a la cima de la colina. es decir, dada la posición actual, mueva al píxel máximo vecino, si su valor es mayor que el píxel actual. Deténgase cuando las coordenadas de sus pies hayan dejado de moverse.

Para contrarrestar el problema de orientación, puede tener aproximadamente 8 configuraciones iniciales para las direcciones básicas (Norte, Nordeste, etc.). Ejecute cada uno individualmente y descarte los resultados en los que dos o más dedos terminen en el mismo píxel. Voy a pensar en esto un poco más, pero este tipo de cosas aún se están investigando en el procesamiento de imágenes: ¡no hay respuestas correctas!

Idea un poco más compleja: (ponderado) K-means clustering. No está tan mal.

  • Comience con cinco coordenadas del dedo del pie, pero ahora estos son "centros del grupo".

Luego itera hasta la convergencia:

  • Asigne cada píxel al clúster más cercano (solo haga una lista para cada grupo).
  • Calcule el centro de masa de cada grupo. Para cada grupo, esto es: Suma (valor de intensidad de coordenadas *) / Suma (coordenadas)
  • Mueva cada grupo al nuevo centro de masa.

Es casi seguro que este método dará resultados mucho mejores y obtendrá la masa de cada grupo que puede ayudar a identificar los dedos de los pies.

(De nuevo, ha especificado el número de clústeres por adelantado. Con el clúster debe especificar la densidad de una manera u otra: elija el número de clústeres, apropiado en este caso, o elija un radio de clúster y vea cuántos finaliza con. Un ejemplo de esto último es cambio medio.)

Perdón por la falta de detalles de implementación u otros detalles. Codificaría esto, pero tengo una fecha límite. Si nada más ha funcionado la próxima semana, házmelo saber y lo intentaré.


26
2017-09-10 22:49



Este problema ha sido estudiado en profundidad por los físicos. Hay una buena implementación en RAÍZ. Mira el TSpectrum clases (especialmente TSpectrum2 para su caso) y la documentación para ellos.

Referencias

  1. M.Morhac et al .: Métodos de eliminación de fondo para espectros de rayos gamma de coincidencia multidimensional. Nuclear Instruments and Methods in Physics Research A 401 (1997) 113-132.
  2. M.Morhac et al .: Deconvolución de oro eficiente de una y dos dimensiones y su aplicación a la descomposición de los espectros de rayos gamma. Nuclear Instruments and Methods in Physics Research A 401 (1997) 385-408.
  3. M.Morhac et al .: Identificación de picos en espectros de rayos gamma de coincidencia multidimensional. Instrumentos y métodos nucleares en Física de la investigación A 443 (2000), 108-125.

... y para aquellos que no tienen acceso a una suscripción a NIM:


12
2017-09-10 12:38



Solo un par de ideas fuera de mi cabeza:

  • tomar el gradiente (derivado) del escaneo, ver si eso elimina las llamadas falsas
  • tomar el máximo de los máximos locales

También es posible que desee echar un vistazo a OpenCV, tiene una API de Python bastante decente y podría tener algunas funciones que le resultarían útiles.


9
2017-09-10 13:05



Aquí hay una idea: calcula el Laplaciano (discreto) de la imagen. Esperaría que fuera (negativo y) grande en máxima, de una manera más dramática que en las imágenes originales. Por lo tanto, los máximos podrían ser más fáciles de encontrar.

Aquí hay otra idea: si conoce el tamaño típico de los puntos de alta presión, primero puede suavizar su imagen al convulsionarla con un gaussiano del mismo tamaño. Esto puede darle imágenes más simples para procesar.


9
2017-09-11 01:07



gracias por los datos sin procesar. Estoy en el tren y esto es lo más lejos que he llegado (mi parada se acerca). He masajeado su archivo txt con expresiones regulares y lo he colgado en una página html con algunos javascript para la visualización. Lo estoy compartiendo aquí porque algunos, como yo, pueden encontrarlo más fácilmente pirateable que Python.

Creo que un buen enfoque será la escala y la rotación invariables, y mi próximo paso será investigar las mezclas de gaussianos. (Cada almohadilla de pata es el centro de un gaussiano).

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

alt text


7
2017-09-10 19:24