Обнаружение пиков в двумерном массиве

874

Я помогаю ветеринарной клинике измерять давление под лапой собаки. Я использую Python для анализа данных, и теперь я застрял, пытаясь разделить лапы на (анатомические) субрегионы.

Я сделал двумерный массив каждой лапы, который состоит из максимальных значений для каждого датчика, который был загружен лапой с течением времени. Вот пример одной лапы, где я использовал Excel, чтобы нарисовать области, которые я хочу «обнаружить». Это 2 на 2 поля вокруг датчика с локальными максимумами, которые вместе имеют наибольшую сумму.

альтернативный текст

Поэтому я попытался поэкспериментировать и решил просто искать максимумы каждого столбца и строки (не могу смотреть в одном направлении из-за формы лапы). Похоже, это «хорошо» определяет местоположение отдельных пальцев, но также отмечает соседние датчики.

альтернативный текст

Итак, что будет лучшим способом сказать Python, какие из этих максимумов те, которые я хочу?

Примечание: квадраты 2х2 не могут перекрываться, так как они должны быть отдельными пальцами!

Кроме того, я взял 2x2 для удобства, любое более продвинутое решение приветствуется, но я просто ученый по человеческому движению, поэтому я не являюсь настоящим программистом или математиком, поэтому, пожалуйста, держите его «простым».

Вот версия, которая может быть загружена сnp.loadtxt


Результаты

Поэтому я попробовал решение @ jextee (см. Результаты ниже). Как вы можете видеть, он очень хорошо работает на передних лапах, но хуже работает на задних лапах.

Более конкретно, он не может распознать маленький пик, который является четвертым пальцем ноги. Это, очевидно, присуще тому факту, что цикл смотрит сверху вниз в направлении наименьшего значения, не принимая во внимание, где это находится.

Кто-нибудь знает, как настроить алгоритм @ jextee, чтобы он тоже мог найти 4-й палец?

альтернативный текст

Поскольку я еще не обработал другие испытания, я не могу предоставить другие образцы. Но данные, которые я давал раньше, были средними для каждой лапы. Этот файл представляет собой массив с максимальными данными 9 лап в порядке их контакта с пластиной.

Это изображение показывает, как они были пространственно распределены по пластине.

альтернативный текст

Обновить:

Я создал блог для всех, кто интересуется, и я установил SkyDrive со всеми необработанными измерениями. Так что для любого, кто запрашивает больше данных: больше власти для вас!


Новое обновление:

Таким образом , после помощи я получил мои вопросы , касающиеся обнаружений лапы и лапа сортировка , я, наконец , смог проверить обнаружение схождения для каждой лапы! Оказывается, это ни на что не годится, кроме лап размером с ту, что была в моем собственном примере. Конечно, задним числом, я сам виноват в том, что выбрал 2х2 произвольно.

Вот хороший пример того, как все идет не так, как надо: гвоздь распознается как носок, а пятка такая широкая, что его распознают дважды!

альтернативный текст

Лапа слишком большая, поэтому размер 2х2 без перекрытия приводит к тому, что некоторые пальцы ног обнаруживаются дважды. С другой стороны, у маленьких собак часто не удается найти пятый палец, что, как я подозреваю, вызвано слишком большой площадью 2х2.

После пробного текущего решения на всех моих измерениях я пришел к ошеломляющему выводу , что для почти всех моих маленьких собак он не нашел 5 - й палец и что в более чем 50% от последствий для больших собак были бы найти больше!

Ясно, что мне нужно это изменить. Мое собственное предположение было изменение размера neighborhoodна что-то меньшее для маленьких собак и больше для больших собак. Но generate_binary_structureне позволил бы мне изменить размер массива.

Поэтому я надеюсь, что у кого-то еще есть лучшее предложение для расположения пальцев ног, возможно, с масштабом области пальцев ног в зависимости от размера лап?

Ivo Flipse
источник
Я так понимаю, что запятые - это десятичные разряды, а не разделители значений?
MattH
Да, они запятые. И @Christian, я пытаюсь вставить его в легко читаемый файл, но даже это не помогает мне :(
Ivo Flipse
3
Поскольку я делаю технико-экономическое обоснование, все идет действительно. Поэтому я ищу как можно больше способов определения давления, включая субрегионы. Также мне нужно уметь различать стороны «большого пальца» и «маленького пальца», чтобы оценить ориентацию. Но так как это не было сделано раньше, мы не можем сказать, что мы можем найти :-)
Ivo Flipse
2
@ Рон: одна из целей этого исследования - выяснить, для какого размера / веса собак подходит система, так что да, пока эта собака весила около 20 кг. У меня есть некоторые, которые значительно меньше (и больше) и ожидаю, что я не смогу сделать то же самое для настоящих маленьких.
Иво Флипс
2
@ frank лапы измеряются во времени, отсюда и 3-е измерение. Тем не менее, они не двигаются со своего места (условно говоря), поэтому меня больше всего интересует, где пальцы находятся в 2D. После этого 3D-аспект становится бесплатным
Ivo Flipse,

Ответы:

332

Я обнаружил пики, используя локальный фильтр максимума . Вот результат для вашего первого набора данных из 4 лап: Результат обнаружения пиков

Я также управлял им на втором наборе данных из 9 лап, и это работало также .

Вот как вы это делаете:

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

Все, что вам нужно сделать, это использовать scipy.ndimage.measurements.labelмаску для маркировки всех отдельных объектов. Тогда вы сможете играть с ними индивидуально.

Обратите внимание, что метод работает хорошо, потому что фон не шумный. Если бы это было так, вы бы обнаружили кучу других нежелательных пиков на заднем плане. Другим важным фактором является размер окрестности . Вам нужно будет отрегулировать его, если размер пика изменяется (он должен оставаться примерно пропорциональным).

Иван
источник
1
Существует более простое решение, чем (eroded_background ^ local_peaks). Просто сделай (передний план и местные пики)
Райан Сокласки
53

Решение

Файл данных: paw.txt . Исходный код:

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"

Вывод без перекрывающихся квадратов. Кажется, что выделены те же области, что и в вашем примере.

Некоторые комментарии

Сложная задача - вычислить суммы всех квадратов 2х2. Я предположил, что вам нужны все из них, так что могут быть некоторые совпадения. Я использовал срезы, чтобы вырезать первые / последние столбцы и строки из исходного 2D-массива, а затем накладывать их все вместе и вычислять суммы.

Чтобы понять это лучше, представим массив 3х3:

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

Тогда вы можете взять его ломтиками:

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

Теперь представьте, что вы сложили их один над другим и суммировали элементы в одинаковых позициях. Эти суммы будут точно такими же суммами для квадратов 2x2 с верхним левым углом в той же позиции:

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

Когда у вас есть суммы более 2х2 квадратов, вы можете использовать, maxчтобы найти максимум или sort, или, sortedчтобы найти пики.

Чтобы запомнить положения вершин, я соединяю каждое значение (сумму) с его порядковым положением в уплощенном массиве (см. zip). Затем я снова вычисляю положение строки / столбца при печати результатов.

Ноты

Я учел перекрывать квадраты 2х2. Отредактированная версия отфильтровывает некоторые из них так, что в результатах отображаются только непересекающиеся квадраты.

Выбор пальцев (идея)

Другая проблема состоит в том, как выбрать то, что, вероятно, будет пальцами из всех пиков. У меня есть идея, которая может или не может работать. У меня нет времени, чтобы реализовать это прямо сейчас, так что просто псевдокод.

Я заметил, что если передние пальцы остаются на почти идеальном круге, задний палец должен быть внутри этого круга. Кроме того, передние пальцы расположены более или менее равномерно. Мы можем попытаться использовать эти эвристические свойства для обнаружения пальцев.

Псевдокод:

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

Это подход грубой силы. Если N относительно мало, то я думаю, что это выполнимо. Для N = 12, есть C_12 ^ 5 = 792 комбинации, умноженные на 5 способов выбора заднего пальца, поэтому 3960 случаев нужно оценить для каждой лапы.

sastanin
источник
Ему придется отфильтровать лапы вручную, учитывая ваш список результатов ... выбор четырех лучших результатов даст ему четыре возможности построить квадрат 2х2, содержащий максимальное значение 6,8
Йоханнес Чарра
Ящики 2x2 не могут перекрываться, так как, если я хочу делать статистику, я не хочу использовать один и тот же регион, я хочу сравнивать регионы :-)
Ivo Flipse
Я отредактировал ответ. Теперь в результатах нет перекрывающихся квадратов.
Састанин
1
Я попробовал это, и это, кажется, работает для передних лап, но меньше для задних. Думаю, нам придется попробовать что-то, что знает, где искать
Иво Флипс
1
Я объяснил свою идею, как можно обнаружить пальцы в псевдокоде. Если вам это нравится, я могу попытаться реализовать это завтра вечером.
Састанин
34

Это проблема с регистрацией изображений . Общая стратегия:

  • Имейте известный пример, или какой-то предварительный на данных.
  • Подгоните ваши данные к примеру или подгоните пример к вашим данным.
  • Это помогает, если ваши данные примерно выровнены в первую очередь.

Вот грубый и готовый подход , «самая глупая вещь, которая могла бы работать»:

  • Начните с пяти пальцевых координат примерно в том месте, где вы ожидаете.
  • С каждым итеративно поднимайтесь на вершину холма. то есть, учитывая текущую позицию, перейдите к максимальному соседнему пикселю, если его значение больше текущего пикселя. Остановитесь, когда ваши координаты пальцев перестанут двигаться.

Чтобы противодействовать проблеме ориентации, у вас может быть около 8 начальных настроек для основных направлений (Север, Северо-Восток и т. Д.). Запустите каждый из них по отдельности и отбросьте любые результаты, когда два или более пальцев оказываются в одном пикселе. Я подумаю об этом еще немного, но такие вещи все еще исследуются в обработке изображений - нет правильных ответов!

Немного более сложная идея: (взвешенная) K-означает кластеризацию. Это не так уж плохо.

  • Начните с пяти пальцевых координат, но теперь это «кластерные центры».

Затем итерации до сходимости:

  • Назначьте каждый пиксель ближайшему кластеру (просто составьте список для каждого кластера).
  • Рассчитайте центр масс каждого кластера. Для каждого кластера это: сумма (координата * значение интенсивности) / сумма (координата)
  • Переместите каждый кластер в новый центр масс.

Этот метод почти наверняка даст гораздо лучшие результаты, и вы получите массу каждого кластера, которая может помочь в определении пальцев.

(Опять же, вы указали количество кластеров заранее. С кластеризацией вы должны указать плотность так или иначе: либо выберите количество кластеров, подходящее в этом случае, либо выберите радиус кластера и посмотрите, сколько вы заканчиваете вверх с. Примером последнего является среднее смещение .)

Извините за отсутствие деталей реализации или других особенностей. Я бы написал это, но у меня есть крайний срок. Если к следующей неделе ничего не получится, дайте мне знать, и я сделаю это.

CakeMaster
источник
1
Проблема в том, что лапы меняют свою ориентацию, и у меня нет калибровки / базовой линии правильной лапы для начала. Кроме того, я боюсь, что многие алгоритмы распознавания изображений немного не в моей лиге.
Иво Флипс
«Грубый и готовый» подход довольно прост - возможно, я не очень хорошо понял эту идею. Я добавлю псевдокод для иллюстрации.
CakeMaster
У меня есть ощущение, что ваше предложение поможет исправить распознавание задних лап, я просто не знаю, «как»
Иво Флипс
Я добавил еще одну идею. Кстати, если у вас есть куча хороших данных, было бы здорово выложить их где-нибудь в сети. Это может быть полезно для людей, изучающих обработку изображений / машинное обучение, и вы можете получить еще немного кода из этого ...
CakeMaster
1
Я просто думал о том, чтобы записать свою обработку данных в простой блог Wordpress, просто чтобы быть полезным для других, и я все равно должен записать это. Мне нравятся все ваши предложения, но я боюсь, что мне придется ждать кого-то без крайнего срока ;-)
Ivo Flipse
18

Используя постоянную гомологию для анализа вашего набора данных, я получаю следующий результат (нажмите, чтобы увеличить):

Результат

Это 2D-версия метода обнаружения пиков, описанного в этом ответе SO . На приведенном выше рисунке просто показаны 0-мерные классы постоянных гомологий, отсортированные по постоянству.

Я увеличил исходный набор данных в 2 раза, используя scipy.misc.imresize (). Тем не менее, обратите внимание, что я рассматривал четыре лапы как один набор данных; разделение его на четыре облегчит проблему.

Методология. Идея этого довольно проста: рассмотрим функциональный график функции, которая назначает каждому пикселю свой уровень. Это выглядит так:

График функции 3D

Теперь рассмотрим уровень воды на высоте 255, который постоянно опускается до более низких уровней. На местных максимумах всплывают острова (роды). В седловых точках два острова сливаются; мы считаем, что нижний остров сливается с высшим островом (смерть). Так называемая диаграмма персистентности (классов гомологии 0-го уровня, наши острова) отображает значения смертности по рождению всех островов:

Диаграмма постоянства

настойчивость островка тогда разница между birth- и смерти уровня; вертикальное расстояние от точки до серой главной диагонали. Фигура обозначает острова, уменьшая постоянство.

Самая первая картинка показывает места рождения островов. Этот метод не только дает локальные максимумы, но также количественно определяет их «значимость» по вышеупомянутой стойкости. Затем можно было бы отфильтровать все острова со слишком низкой устойчивостью. Однако в вашем примере каждый остров (т. Е. Каждый локальный максимум) - это вершина, которую вы ищете.

Код Python можно найти здесь .

С. Хубер
источник
16

Эта проблема была подробно изучена физиками. В ROOT есть хорошая реализация . Посмотрите на классы TSpectrum (особенно TSpectrum2 для вашего случая) и документацию для них.

Ссылки:

  1. М. Морхач и др .: Методы устранения фона для многомерного совпадения гамма-спектров. Ядерные приборы и методы в физических исследованиях A 401 (1997) 113-132.
  2. М. Морхак и др .: Эффективная одно- и двумерная деконволюция золота и ее применение к разложению гамма-спектров. Ядерные приборы и методы в физических исследованиях А 401 (1997) 385-408.
  3. М. Морхач и др .: Идентификация пиков в многомерных спектрах гамма-излучения совпадений. Ядерные приборы и методы в физике исследований А 443 (2000), 108-125.

... и для тех, у кого нет доступа к подписке на NIM:

dmckee --- котенок экс-модератора
источник
Если взглянуть на статью, то здесь описана та же обработка данных, что и здесь, но я боюсь, что она значительно превзошла мои навыки программирования :(
Ivo Flipse
@ Ivo: я никогда не пытался реализовать это сам. Я просто использую ROOT. Тем не менее, существуют привязки Python, но имейте в виду, что ROOT - довольно тяжелый пакет.
dmckee --- котенок экс-модератора
@ Ivo Flipse: я согласен с dmckee. У вас много перспективных ответов в других ответах. Если все они терпят неудачу и вам хочется потратить некоторое время, вы можете вникнуть в ROOT, и он (вероятно) сделает то, что вам нужно. Я никогда не знал никого, кто пытался изучать ROOT через привязки Python (а не на естественном C ++), поэтому я желаю вам удачи.
физика Майкл
13

Вот идея: вы вычисляете (дискретный) лапласиан изображения. Я ожидал бы, что он будет (отрицательным и) большим в максимумах, таким способом, который будет более драматичным, чем на исходных изображениях. Таким образом, максимумы легче найти.

Вот еще одна идея: если вы знаете типичный размер пятен высокого давления, вы можете сначала сгладить ваше изображение, свернув его с гауссианом того же размера. Это может дать вам более простые изображения для обработки.

Эрик О Лебиго
источник
11

Просто пара идей с моей головы:

  • взять градиент (производную) сканирования, посмотреть, устраняет ли это ложные вызовы
  • взять максимум локальных максимумов

Возможно, вы также захотите взглянуть на OpenCV , у него довольно приличный Python API и могут быть некоторые функции, которые вы найдете полезными.

ChrisC
источник
С градиентом, вы имеете в виду, что я должен рассчитывать крутизну склонов, как только эта величина превысит определенное значение, я знаю, что есть «пик»? Я попробовал это, но некоторые пальцы имеют очень низкие пики (1,2 Н / см) по сравнению с некоторыми другими (8 Н / см). Так, как я должен обращаться с пиками с очень низким градиентом?
Иво Флипс
2
Что работало для меня в прошлом, если я не мог использовать градиент напрямую, так это смотреть на градиент и максимумы, например, если градиент является локальным экстремумом, а я на локальных максимумах, то я нахожусь в точке интерес.
ChrisC
11

Я уверен, что у вас есть достаточно, чтобы продолжить, но я не могу не предложить использовать метод кластеризации k-средних. k-means - это алгоритм неконтролируемой кластеризации, который будет собирать ваши данные (в любом количестве измерений - мне случается делать это в 3D) и распределять их по k кластерам с различными границами. Здесь хорошо, потому что вы точно знаете, сколько пальцев должно быть у этих клыков.

Кроме того, это реализовано в Scipy, что очень приятно ( http://docs.scipy.org/doc/scipy/reference/cluster.vq.html ).

Вот пример того, что он может сделать для пространственного разрешения 3D-кластеров: введите описание изображения здесь

То, что вы хотите сделать, это немного по-другому (2D и включает в себя значения давления), но я все еще думаю, что вы могли бы попробовать.

astromax
источник
10

спасибо за необработанные данные. Я нахожусь в поезде, и это так далеко, как я получил (моя остановка подходит). Я помассировал ваш txt-файл с помощью регулярных выражений и поместил его на html-страницу с некоторым javascript для визуализации. Я делюсь этим здесь, потому что некоторые, как и я, могут найти его более легко взломать, чем Python.

Я думаю, что хорошим подходом будет инвариант масштаба и ротации, и мой следующий шаг - исследовать смеси гауссиан. (каждая подушечка лап является центром гауссиана).

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

альтернативный текст

Рон
источник
1
Я считаю, что это является доказательством того, что рекомендуемые методики Гаусса могут работать, теперь, если бы кто-то мог доказать это с помощью Python ;-)
Ivo Flipse
8

Решение физика:
определите 5 маркеров лап, идентифицированных по их позициям, X_iи начните их со случайных позиций. Определите некоторую энергетическую функцию, сочетающую некоторую награду за расположение маркеров в положениях лап и некоторое наказание за наложение маркеров; скажем так:

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

( S(X_i)средняя сила в квадрате 2x2 X_i, alfaявляется параметром, который должен быть экспериментально достигнут)

Теперь пора заняться магией Метрополис-Гастингс:
1. Выберите случайный маркер и переместите его на один пиксель в случайном направлении.
2. Рассчитайте dE, разницу энергии, вызванную этим движением.
3. Получите равномерное случайное число от 0 до 1 и назовите его r.
4. Если dE<0или exp(-beta*dE)>r, примите движение и перейдите к 1; если нет, отмените движение и перейдите к 1.
Это должно повторяться до тех пор, пока маркеры не сойдутся в лапы. Бета контролирует сканирование для оптимизации компромисса, поэтому его также следует оптимизировать экспериментально; оно также может постоянно увеличиваться со временем моделирования (имитация отжига).

МБк
источник
Хотите показать, как это будет работать на моем примере? Поскольку я действительно не в математике высокого уровня, мне уже трудно разобраться с предложенной вами формулой :(
Ivo Flipse
1
Это математика в старшей школе, вероятно, моя запись просто запутана. У меня есть план, чтобы проверить это, так что следите за обновлениями.
МБк
4
Я физик элементарных частиц. В течение долгого времени инструмент программного обеспечения в нашей дисциплине назывался PAW, и у него была сущность, связанная с графами, называемая «маркер». Вы можете себе представить, как я запутался в этом ответе в первые пару раз ...
dmckee --- котенок экс-модератора
6

Вот еще один подход, который я использовал, когда делал нечто подобное для большого телескопа:

1) Поиск самого высокого пикселя. Как только у вас это получится, найдите это для наилучшего соответствия для 2x2 (возможно, максимизируя сумму 2x2), или сделайте 2d гауссовое соответствие внутри подобласти, скажем, 4x4 с центром в верхнем пикселе.

Затем установите те 2х2 пикселя, которые вы нашли на ноль (или, может быть, 3х3) вокруг центра пика

вернитесь к пункту 1) и повторяйте до тех пор, пока самый высокий пик не опустится ниже порога шума, или пока у вас не появятся все необходимые пальцы

Паулюс
источник
Хотите поделиться примером кода, который делает это? Я могу следить за тем, что ты пытаешься сделать, но понятия не имею, как самому это кодировать
Ivo Flipse
Я на самом деле пришел работать с Matlab, так что да, это уже поможет. Но если вы используете действительно посторонние функции, мне может быть сложно воспроизвести его на Python
Ivo Flipse
6

Вероятно, стоит попробовать с нейронными сетями, если вы можете создать некоторые обучающие данные ... но для этого нужно много образцов, аннотированных вручную.

antirez
источник
Если это того стоит, я бы не стал комментировать большой образец вручную. Моя проблема была бы: как мне это реализовать, поскольку я ничего не знаю о программировании нейронных сетей
Ivo Flipse
6

грубая схема ...

Вы, вероятно, захотите использовать алгоритм связанных компонентов, чтобы изолировать каждую область лапы. В вики есть достойное описание этого (с некоторым кодом) здесь: http://en.wikipedia.org/wiki/Connected_Component_Labeling

вам придется принять решение о том, использовать 4 или 8 подключений. лично для большинства проблем я предпочитаю 6-связность. в любом случае, после того, как вы выделите каждый «отпечаток лапы» в качестве связанной области, должно быть достаточно легко выполнить итерацию по области и найти максимумы. Найдя максимумы, вы можете многократно увеличивать регион, пока не достигнете заранее определенного порога, чтобы идентифицировать его как заданный «носок».

одна тонкая проблема здесь заключается в том, что как только вы начнете использовать методы компьютерного зрения для определения чего-либо как правой / левой / передней / задней лапы и начнете смотреть на отдельные пальцы ног, вы должны будете принимать во внимание повороты, перекосы и переводы. это достигается путем анализа так называемых «моментов». Есть несколько различных моментов, которые следует учитывать в приложениях для зрения:

центральные моменты: трансляционно-инвариантные нормализованные моменты: скейлинг и трансляционно-инвариантные моменты h: трансляция, масштаб и инвариант вращения

Более подробную информацию о моментах можно найти, выполнив поиск «моменты изображения» в вики.

Джошуа
источник
4

Кажется, вы можете немного обмануть, используя алгоритм Jetxee. Он находит первые три пальца в порядке, и вы должны быть в состоянии угадать, где четвертый основан на этом.

Geoff
источник
4

Интересная проблема. Решение, которое я бы попробовал, заключается в следующем.

  1. Примените фильтр нижних частот, такой как свертка, с двумерной гауссовой маской. Это даст вам кучу (возможно, но не обязательно с плавающей запятой) значений.

  2. Выполните двумерное немаксимальное подавление, используя известный приблизительный радиус каждой подушечки лап (или носка).

Это должно дать вам максимальные позиции без нескольких кандидатов, которые находятся близко друг к другу. Просто чтобы уточнить, радиус маски на шаге 1 должен быть аналогичен радиусу, используемому на шаге 2. Этот радиус может быть выбран, или ветеринар может явно измерить его заранее (он будет зависеть от возраста / породы / и т. Д.).

Некоторые из предложенных решений (среднее смещение, нейронные сети и т. Д.), Вероятно, будут работать в некоторой степени, но слишком сложны и, вероятно, не идеальны.

Боб Моттрам
источник
У меня 0 опыта работы с матрицами свертки и гауссовыми фильтрами, так что вы хотели бы показать, как это будет работать на моем примере?
Иво Флипс
3

Ну, вот несколько простых и не очень эффективных кодов, но для такого размера набора данных это хорошо.

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]

Я просто делаю массив с положением верхнего левого угла и суммой каждого квадрата 2х2 и сортирую его по сумме. Затем я беру квадрат 2х2 с наибольшей суммой из спора, кладу его вbest массив и удаляю все остальные квадраты 2x2, которые использовали любую часть этого только что удаленного квадрата 2x2.

Кажется, что он работает нормально, за исключением последней лапы (с самой маленькой суммой в правом нижнем углу на вашей первой картинке), оказывается, что есть два других приемлемых квадрата 2x2 с большей суммой (и они имеют равную сумму друг с другом). Один из них по-прежнему выбирает один квадрат из вашего квадрата 2x2, а другой - слева. К счастью, по счастливой случайности мы видим, что выбираем больше того, что вам хотелось бы, но для этого может потребоваться использовать другие идеи, чтобы получить то, что вы на самом деле хотите, все время.

Джастин Пил
источник
Я считаю, что ваши результаты такие же, как в ответе @ Jextee. Или, по крайней мере, мне кажется, что я это тестирую.
Иво Флипс
3

Просто хочу сказать вам, ребята, что есть хорошая опция для поиска локальных maximaизображений в python:

from skimage.feature import peak_local_max

или для лыжного мага 0.8.0:

from skimage.feature.peak import peak_local_max

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

Thomio
источник
1

Возможно, здесь достаточно наивного подхода: составьте список всех квадратов 2х2 на плоскости, упорядочите их по сумме (в порядке убывания).

Сначала выберите самый ценный квадрат в вашем «списке лап». Затем, итеративно выберите 4 из следующих лучших квадратов, которые не пересекаются ни с одним из ранее найденных квадратов.

Йоханнес Чарра
источник
Я на самом деле составил список со всеми суммами 2х2, но когда я их заказал, я не знал, как их итеративно сравнить. Моя проблема заключалась в том, что, когда я сортировал это, я потерял трек координат. Возможно, я мог бы вставить их в словарь с координатами в качестве ключа.
Иво Флипс
Да, какой-то словарь был бы необходим. Я бы предположил, что ваше представление сетки уже является своего рода словарем.
Йоханнес Чарра
Ну, изображение, которое вы видите выше, является массивом. Остальное в настоящее время хранится в многомерных списках. Вероятно, было бы лучше прекратить это делать, хотя я не так хорошо разбираюсь в словарях
Иво Флипс
1

Существует несколько обширных программных продуктов, доступных от астрономического и космологического сообщества - это значительная область исследований как исторически, так и в настоящее время.

Не пугайтесь, если вы не астроном - некоторые из них легко использовать за пределами поля. Например, вы можете использовать astropy / photutils:

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

[Кажется немного грубым повторять их короткий пример кода здесь.]

Неполный и слегка предвзятый список техник / пакетов / ссылок, которые могут представлять интерес, приведен ниже - добавьте больше комментариев, и я буду обновлять этот ответ по мере необходимости. Конечно, есть компромисс между точностью и вычислительными ресурсами. [Честно говоря, их слишком много, чтобы привести примеры кода в одном ответе, таком как этот, поэтому я не уверен, будет ли этот ответ вылетать или нет.]

Source Extractor https://www.astromatic.net/software/sextractor

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

Задача ASKAP / EMU по поиску источника: https://arxiv.org/abs/1509.03931

Вы также можете искать проблемы с извлечением источника в Планке и / или WMAP.

...

jtlz2
источник
0

Что делать, если вы продолжите шаг за шагом: сначала найдите глобальный максимум, обработайте, если необходимо, окружающие точки, учитывая их значение, затем установите найденный регион на ноль и повторите для следующего.

Седрик Х.
источник
Хм, что установка на ноль, по крайней мере, удалит его из любых дальнейших расчетов, что будет полезно.
Иво Флипс
Вместо установки на ноль, вы можете вычислить гауссову функцию с выбранными вручную параметрами и вычесть найденные значения из первоначальных значений давления. Таким образом, если палец сжимает ваши датчики, то, найдя самую высокую точку нажатия, вы используете его, чтобы уменьшить воздействие этого пальца на датчики, тем самым устраняя соседние ячейки с высокими значениями давления. en.wikipedia.org/wiki/File:Gaussian_2d.png
Данияр,
Хотите показать пример, основанный на моих данных образца @Daniyar? Поскольку я действительно не знаком с такой обработкой данных
Ivo Flipse
0

Я не уверен, что это отвечает на вопрос, но кажется, что вы можете просто искать n самых высоких вершин, у которых нет соседей.

Вот суть. Обратите внимание, что это в Ruby, но идея должна быть ясной.

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
Рик Халл
источник
Я собираюсь попробовать и посмотреть, смогу ли я преобразовать его в код Python :-)
Ivo Flipse
Пожалуйста, включите код в сам пост, а не ссылку на суть, если это разумная длина.
agf