Pull to refresh

В поисках изофот

Reading time 5 min
Views 3K
Понадобилось мне однажды вычисление изофот (линий равной интенсивности на изображениях), однако, готовых библиотек я не нашел, а копаться в чужом коде (например, в тех же Octave или Iraf) очень не хотелось. В качестве простейшего алгоритма я нашел метод шагающих квадратов. Однако, этот метод сильно снижает пространственное разрешение при вычислении узлов изофот, поэтому я решил его немного видоизменить и сделать квадраты перекрывающимися.
Первая реализация (кстати, довольно быстрая) была неудачной: т.к. я строил маску, общую для всех уровней, а для конкретного уровня пересчитывал отдельно, в точках, через которые проходит несколько изолиний, получились разрывы. Ковыряния в коде приводили к все большему и большему числу взаимных блокировок и флагов, поэтому я решил пойти в ущерб производительности и вычислять маски отдельно для каждого уровня интенсивности.

Итак, вкратце повторю, в чем заключается метод «шагающих квадратов». Для построения изофот уровня Ii вначале строится специальная маска: для каждого пикселя изображения кодируется относительное положение «соседей» справа, внизу и внизу справа. Кодировать это можно по-разному, я просто представил код в виде битовой маски 000abcd, где a,b,c и d — относительные интенсивности в соответствующих пикселях (a — наш текущий пиксель изображения)
build mask
Если значение в пикселе больше уровня изолинии, соответствующему биту присваивается единица, иначе — ноль. Если w — ширина изображения, а *in — указатель на пиксель a, то значение маски (*out) для этого пикселя вычисляется так:
*out = (unsigned char) ((in[0]>lvl)<<3) | ((in[1]>lvl)<<2) | ((in[w]>lvl)<<1) | (in[w+1]>lvl);

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

Следующим шагом является поочередная проверка всех пикселей маски на неравенство нулю и 15. Как только мы находим такой пиксель, начинаем строить от него изолинию. Для того, чтобы корректно построить изолинии, нам необходимо запоминать направление, с которого мы к данному пикселю подошли. На рисунке выше стрелками обозначены возможные пути «схода» с каждого пикселя маски. Если мы не будем учитывать запрещенное направление (с которого мы пришли), придется либо проверять лишние пиксели маски, либо же (в случае с особыми точками 6 и 9, но о них чуть позже) мы можем вообще неверно выбрать положение следующей точки изофоты.
По мере просмотра пикселей маски (слева направо и сверху вниз, если считать начало координат изображения в левом верхнем углу), пока мы находим лишь пустые ячейки, «запрещенным» направлением является направление «влево». Как только мы находим ячейку, содержащую участок изолинии, мы, в соответствии с заранее подготовленным массивом «направлений» выбираем следующее направление, отсекая «запрещенное». Если таких направлений два, то мы выбираем вправо или вниз (с приоритетом на «вправо») при движении по часовой стрелке и влево или вверх (с приоритетом на «вверх») при движении против часовой стрелки (об этом — позже).
После этого мы вычисляем для конкретного пикселя положение узла изолинии. Это делается при помощи одной из массива функций, в соответствии со значением маски в данном пикселе. Рисунок снизу поясняет, как вычисляется это положение.
point calculation
Итак, пусть изображение имеет интенсивности, изображенные на рисунке (цифры внутри ячеек). Мы проводим изофоту по уровню 30. Конфигурация соседей пикселя (xi, yi) соответствует значению маски, равному семи. Мы вызываем соответствующую функцию, которая методом линейной интерполяции вычисляет координаты опорных точек изолинии на прямых, соединяющих по вертикали правые два пикселя и по горизонтали нижние два. Соответствующие точки отмечены незакрашенными красными кружками. Соединяя их прямой и находя ее центр, мы получим искомый узел изофоты (отмечен красным кружком с зеленой заливкой, «Our point»). Это, конечно, не очень-то корректный метод, но для большинства задач его достаточно. Каждый узел добавляется в «хвост» соответствующего списка.
После того, как соответствующее значение маски использовано, оно обнуляется и мы переходим к следующему соседнему пикселю маски. Если его содержимое соответствует очередному узлу изофоты, мы продолжаем вычисление, иначе — «закрываем» изофоту. Сразу после «закрытия» изофоты мы проверяем, во-первых, не является ли она слишком короткой (скажем, вокруг одного-единственного пикселя), а во-вторых — проверяем ее замкнутость: если начало и конец изофоты лежат не далее, чем, скажем, 2 пикселя друг от друга, мы считаем ее замкнутой.
Процедура «закрывания» изофоты заключается в том, что «на всякий случай» проверяется, не продолжается ли изофота от своей первой точки в противоположном направлении. Для этого мы переходим к нижнему от первого пикселю маски с запретом на движение «вверх». Если мы находим там продолжение изофоты, то проходим вплоть до ее окончания, добавляя узлы к «голове» списка и выбирая приоритетными направления против часовой стрелки.

Все вышесказанное справедливо, пока мы мы не «натыкаемся» на особые пиксели: 6 и 9. Через эти пиксели маски проходит две линии изофоты, поэтому необходимо, на основе «запрещенного» направления выбрать правильное следующее направление. Так, к примеру, если мы «натыкаемся» на значение 6, придя слева, следующим направлением будет «вверх».
Чтобы «не забыть» о еще одной точке изолинии при следующем проходе, мы заменяем особое значение на «не особое», просто «закрасив» уже использованный пиксель (т.е., например, для случая попадания в «шестерку» слева, мы «закрашиваем» верхний левый пиксель, т.е. заменяем 6 на 7). Такой способ, к сожалению, вызывает «промахи» в случае наличия однопиксельных «провалов» в «горбах» интенсивности. Но, к сожалению, исправление этого — очень сложная алгоритмическая задача, которая будет отбирать слишком много машинного времени на построение изофот. А для большинства изображений описанный метод отлично работает.
Таким образом, «особые» пиксели используются дважды, обнуляясь лишь после второго прохода. Вычисление узлов в особых точках выполняется наоборот: для этого берется функция с «закрашиванием» противоположного используемому пикселя.

Проиллюстрирую работу алгоритма на примере обхода «пика» из двух пикселей, изображенного на рисунке ниже.
example
Приведенному участку изображения 4х4 пикселя соответствует маска 3х3 пикселя (справа). Мы начинаем движение с левого верхнего пикселя изображения. Т.к. на маске ему соответствует ноль, мы переходим к соседнему правому пикселю маски (при этом помним, что у нас «запрет» на направление «вправо»). Запрещенные направления изображены на рис. справа красными стрелками, «пустые» переходы — черной стрелкой, переходы с расчетом узлов изолинии — зелеными стрелками.
Итак, мы «натыкаемся» на восьмерку и вычисляем при помощи соответствующей функции положение первого узла изофоты. Далее мы смотрим в заранее подготовленный массив и видим, что направления «схода» с текущего пикселя — либо «вправо», либо «вниз». Мы выбираем направление «вправо» и переходим к следующему пикселю (4), не забывая обнулить уже использованный пиксель маски (обнуление пикселей отмечено красными нулями в правом верхнем углу пикселя).
Далее мы переходим ко значению 1, а после него — к особой точке со значением 6. Т.к. мы пришли справа, то значение узла вычисляем по функции для значения маски, равного семи, а шестерку на маске заменяем на 14. После этого переходим вниз, к единице.
После того, как мы пройдем единицу, двойку и восьмерку, мы вернемся к нашему «особому» пикселю, но там уже будет значение 14, т.е. дальнейшее вычисление будет происходить «как ни в чем ни бывало», после чего нам надо будет переходить вверх. Однако, вверху мы натыкаемся на ноль. В «голову» списка узлов мы тоже ничего не можем добавить, т.к. значение под начальным пикселем изолинии тоже обнулено. Поэтому мы проверяем близость начального и конечного узлов линии. Они довольно близки, и мы можем пометить данную изолинию как замкнутую.

Все исходники к алгоритму можно найти здесь. Я планирую еще добавить распараллеливание вычисления изолиний разных уровней, чтобы запустить хотя бы 4..8 параллельных потоков. Пока, в один поток, для тестового изображения 3000х3000 пикселей построение изофот по 16 уровням занимает около 1.2с.
Tags:
Hubs:
+25
Comments 16
Comments Comments 16

Articles