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

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

from skimage.feature import peak_local_max

или для 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

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

Не пугайтесь, если вы не астроном - некоторые из них легко использовать вне поля зрения. Например, вы можете использовать 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

Вы также можете искать задачи извлечения источников Planck и / или WMAP.

...

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

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

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

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

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

настойчивость острова - это тогда разница между уровнем рождения и смерти; расстояние от точки до серой главной диагонали по вертикали. На рисунке островки помечаются уменьшением стойкости.

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

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

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

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

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

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

спасибо за необработанные данные. Я в поезде, и это все, что у меня есть (приближается моя остановка). Я обработал ваш 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>

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

Ссылки:

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

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

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

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

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

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

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

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

Решение физика:
Определите 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.
Так повторять до тех пор, пока маркеры не сойдутся в лапы. Бета-версия управляет сканированием для оптимизации компромисса, поэтому его также следует оптимизировать экспериментально; его также можно постоянно увеличивать со временем моделирования (имитация отжига).

Возможно, вы можете использовать что-то вроде Gaussian Mixture Models. Вот пакет Python для выполнения GMM (только что поискал в Google) http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/em/

приблизительный набросок ...

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

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

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

центральные моменты: инвариант переноса нормализованные моменты: инвариант масштабирования и трансляции hu моменты: перенос, масштаб и инвариант вращения

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

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

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

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

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

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

Я не уверен, что это отвечает на вопрос, но похоже, что вы можете просто найти 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

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

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

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

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

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

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

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

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

Затем повторяйте до сходимости:

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

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

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

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

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

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]

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

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

Решение

Файл данных: 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"

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

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

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

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

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

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

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

Примечания

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

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

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

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

Псевдокод:

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 случаев для оценки для каждой лапы.

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

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

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

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

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

Пара идей, которые мне не приходят в голову:

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

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