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

    Понадобилось мне однажды вычисление изофот (линий равной интенсивности на изображениях), однако, готовых библиотек я не нашел, а копаться в чужом коде (например, в тех же 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с.
    Ads
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More

    Comments 16

      +2
      Поделитесь, пожалуйста, результатом работы алгоритма. И, если не секрет, для чего вы искали эти изофоты?
        +6
        Вот — как выглядит в 2D: 2D а вот так — в 3D: 3D

        Этот алгоритм нужен был мне для улучшения метода определения центра тяжести «точечных» объектов, а также для кое-каких других операций (например, вычисления энергии в кружке) — в рамках постепенно разрабатываемой мной программки для простого просмотра/редактирования FITS-файлов и (главное, для чего я ей занялся) восстановления поверхности зеркала по гартманограммам.
          +2
          Простите, не нашел определения «гартманограммы». Это способ контроля качества оптики?
            +4
            Да, контроль качества оптики методом Гартманна (с использованием специальной диафрагмы-маски). Кроме того, если известно состояние оптики и есть быстрая камера, можно анализировать текущее состояние атмосферы и при помощи блока адаптивной оптики (например, активного вторичного зеркала) частично компенсировать атмосферные искажения. Но это еще в будущем для БТА.
        0
        А можно вопрос?
        А чем изофоты отличаются вот от этого?
        image
          +1
          У вас — обычные границы, найденные дифференциальным фильтром (изображение), а мне необходимо получить структуру с координатами каждой изолинии.
          А методов выделения границ, конечно, полным-полно. Если изофоты просто надо нарисовать, не вычисляя узлы, то достаточно квантовать изображение на нужное количество уровней, а затем вычислить модуль градиента, который и отобразит нам все изолинии.
            0
            а такой вариант?
            image
            изображение бинарное 0,1

            и найти замкнутые области

            почему спрашиваю.
            Эта штука невероятно быстрая.
            3000x3000 бдет штробить в риал тайме.
              +1
              Вы можете пройтись по ссылочке на мой ЖЖ, где описывается первый, неудачный подход.
              Считалось очень быстро (т.к. маска была всего одна), но только для изображений, не содержащих проблемных участков (несколько изолиний через 1 пиксель; сложный профиль изолиний и т.п.).
              Можно бинаризовать изображение по каждому искомому уровню и найти затем контуры, но все равно ведь придется пробегаться по каждому контуру и вычислять более-менее точные координаты.
              Так что, основная проблема со временем здесь не в поиске контуров изолиний, а в вычислении их узлов.
                0
                Можно уточнение?
                «придется пробегаться по каждому контуру и вычислять более-менее точные координаты.»

                получив такое изображение как приведено выше, можно говорить о точных координатах белых точек/линий.
                как мы можем видеть, ширина этих линий точно в 1 пиксель.

                А что есть «узлы» в вашем понимании?
                  +1
                  получив такое изображение как приведено выше, можно говорить о точных координатах белых точек/линий.
                  как мы можем видеть, ширина этих линий точно в 1 пиксель.

                  При том, что координата изолинии может проходить, скажем, через точку (1.001, 100.097). А округление до пикселя слишком угрубит вычисления. Так можно делать только в очень редких случаях.
                  Ну, а если вам, например, надо с точностью до одной тысячной пикселя посчитать центр тяжести изображения, придется и контуры вычислять более-менее надежно.

                  А что есть «узлы» в вашем понимании?

                  Это просто некоторые точки с высокой надежностью близкие к реально проходящей изолинии.
                    +1
                    понятно
                    я упустил тот факт что у вас «изображение» это уже «производная» от некоторых данных.

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

                    Но метод о котором я говорю один из самых точных, что я видел.
                    img-service.com/overview/cellular_neural_boolean_filtering.html
                    habrahabr.ru/blogs/image_processing/128803/

                    Там же ссылки на научную основу.
                      +1
                      я упустил тот факт что у вас «изображение» это уже «производная» от некоторых данных

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

                      Сегодня не могу: с работы уже ушел, а дома интернет очень медленный. Изображения же — порядка 18МБ в несжатом виде, в сжатом — магабайт по 10 на изображение будет.
                      Но метод о котором я говорю один из самых точных, что я видел

                      Не забывайте, что у меня — не поиск границ и/или резких перепадов интенсивности, а поиск линий одинакового уровня. Это несколько иная задача.
                        0
                        Да, добавлю еще, что зачастую после отыскания узлов каждой изофоты применяют еще и сглаживание — B-сплайнами, например.
            0
            С того момента, как вы упомянули зеркала, я был уверен, что вы связаны с астрофизикой и ваш жж это только подтверждает. Я учусь на физика-теоретика и лично мне было бы очень интересно знать, какие задачи может найти себе программист в космосе.
              +2
              Задач достаточно много: и в области теоретической физики, и в прикладной области. Обработка изображений до сих пор активно развивается, нужны новые методы, позволяющие упростить рутинную работу и повысить точность получаемых данных.
              Всякое машинное зрение, теории управления и т.д., и т.п. тоже нужны.

              В общем, много чего надо сделать, и как всегда времени на все не хватает, вот так и пользуемся полусобранными приборами, полуфункционирующим софтом и т.п.
              +1
              Попробовал вчера распараллелить (вычислять в четырех отдельных потоках изолинии для каждого уровня) — выигрыш довольно небольшой. Если без распараллеливания изолинии тестового изображения 3Кх3К по 16 подуровням строились за 1.2с, то с распараллеливанием — за 0.75с.
              Похоже, все-таки довольно много времени отбирает построение маски (которое и так распараллелено).
              Буду думать на досуге, как видоизменить алгоритм, чтобы маски строились быстрее (например, маску уровня i+1 можно строить на основе маски уровня i, проверяя лишь ненулевые пиксели предыдущей маски — тогда вычислений и разыменования указателей будет меньше).

              Only users with full accounts can post comments. Log in, please.