company_banner

Вычисление центра масс за O(1) с помощью интегральных изображений



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

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

    В этой статье я расскажу:

    • Что за задача такая, о которой идет речь;
    • Подробнее об интегральных изображениях;
    • Как использовать интегральные изображения для приближенного решения гравитационной задачи N тел применительно к дискретному полю импульсов (масс-скоростей);
    • Какой недостаток имеет это решение и как его исправить;
    • И, наконец, как за константное время вычислить центр масс для произвольного региона.

    Гравитационная задача N тел


    Строгое решение гравитационного взаимодействия системы из N тел относится к алгоритмам, имеющим квадратичную сложность: на каждое из тел оказывают гравитационное воздействие все остальные тела в системе[1]. Следовательно, силу гравитационного взаимодействия нужно посчитать для каждой из N2/2 пар.

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

    Для вычисления силы тяготения между двумя телами в трехмерном пространстве необходимо выполнить 20 операций с плавающей точкой (FLOP). В Солнечной системе насчитывается около 100 тел размером более 200 километров (включая Солнце, 9 планет, карликовые планеты и спутники)[2]. Приблизительное число астероидов на околоземной орбите составляет около 20 тысяч[3]. В поясе астероидов насчитывается от 1,1 до 1,9 миллионов тел размером больше 1 км[4] и порядка одного миллиона астероидов похожего размера в «троянской»[5] и «греческой» группах Юпитера. В поясе Койпера порядка 100 миллионов объектов размером более 20 км[6] и еще около 2 триллионов ― в облаке Оорта[7].



    Вычислительная мощность, требуемая для строгого решения последней гравитационной задачи, всего на порядок меньше вычислительной мощности, требуемой для симуляции работы человеческого мозга на нейронном уровне (2,5×1026)[8]. Именно поэтому ее строгое решение обычно заменяется приближенным.

    Один из наиболее широко используемых алгоритмов для приближенного решения гравитационной задачи ― Tree Code ― имеет квазилинейную временную сложность O(N*logN)[9]. Tree Code строит пространственное дерево и для каждого узла в этом дереве рассчитывает общую массу и центр масс всех тел, которые в него попадают. В дальнейшем, при расчете гравитационных сил, действующих на каждое тело, Tree Code может учитывать не прямое влияние других тел, а влияние узлов дерева, причем в зависимости от дистанции он выбирает все более крупные узлы, которые содержат параметры все более многочисленного подмножества тел.


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

    Гравитационная задача для поля импульсов


    Поле импульсов ― это дискретное векторное поле mv(p), которое каждой точке рассматриваемого пространства p ставит в соответствие пару масса-скорость. Поле импульсов можно задать для пространства любой размерности. Пара масса-скорость характеризует количество движения и представляет собой вектор размерностью R + 1, где R ― размерность пространства, для которого задано поле импульсов. Например, для R = 2 это может быть вектор {vx, vy, m}.

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


    Небольшой участок векторного поля импульсов: массы обозначены цветом, направления ― стрелками


    Векторное поле импульсов размером 729×729 элементов. Слева ― массы. Справа ― импульсы. Импульсы на этой иллюстрации закодированы цветом в формате HSV (Hue ― направление импульса, а Value линейно отображает порядок его абсолютной величины таким образом, что наименее различимы (Value = 0,05) импульсы порядка 10-7, а самые яркие импульсы (Value = 1,0) имеют порядок 103 и выше).

    Подобные этому ― сеточные ― методы описания физических систем повсеместно используются в астрофизике для моделирования эволюции газовых облаков, формирования звезд и галактик. К ним относятся метод SAMR[10] или сеточная модель гидродинамики[11].

    Как и в случае гравитационной задачи N тел, алгоритму эволюции дискретного поля необходимо учесть гравитационное влияние всех масс, составляющих это самое поле, на каждый его дискретный элемент. При этом необходимо принимать во внимание, что сложность задачи косвенно зависит от размерности пространства R: например, для поля на плоскости, состоящего из 1000×1000 элементов, общее число N составит 106 элементов, а поле такого же размера в трехмерном пространстве будет включать в себя уже 109 элементов.

    Тем не менее, существует ряд приемов, позволяющих приближенно решить гравитационную задачу и для поля импульсов. Метод Multigrid использует иерархическую дискретизацию, поддерживая несколько сеток различного размера[12]. Multipole Expansion трактует группы источников, расположенные близко друг к другу, как один объект[13]. Multigrid имеет квазилинейную вычислительную сложность и логарифмическую сложность по памяти. Один из вариантов Multipole Expansion ― FMM ― демонстрирует линейную вычислительную сложность в обмен на пониженную точность вычислений[14].

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

    Интегральное изображение


    Интегральное изображение дает посчитать за константное время сумму элементов оригинального изображения в произвольной прямоугольной области[15]. В компьютерной графике интегральные изображения первоначально были предложены в качестве альтернативы мипмаппингу и анизотропной фильтрации[16]. Кроме того, они успешно используются для цифровой обработки изображений и в методиках распознавания образов.

    Интегральное изображение ― один из самых наглядных примеров компромисса между вычислительной сложностью и сложностью по памяти. «Наивный» алгоритм подсчета суммы элементов изображения обладает временной сложностью O(M*N) и сложностью по памяти O(1). Интегральное изображение позволяет совершать те же вычисления за время O(1) и обладает сложностью по памяти O(M*N).


    Использование интегрального изображения для вычисления суммы элементов в заданном регионе

    Алгоритм расчета интегрального изображения очень прост, характеризуется линейной вычислительной сложностью и легко параллелизируется под GPGPU[17]. Реализуется он подобно двухпроходному фильтру Гаусса[18]: в первом проходе считаются частные суммы для строк, во втором эти частные суммы аккумулируются по столбцам.


    Вычисление интегрального изображения в два прохода. Слева ― оригинальное изображение. В центре ― частные суммы для строк. Справа ― итоговое интегральное изображение.


    Слева ― изображение, содержащее массы элементов. Справа ― его интегральное изображение. Изображения слева и справа используют разные логарифмические шкалы.

    Приближенное решение гравитационной задачи с помощью интегрального изображения


    Имея в своем распоряжении интегральное изображение масс, несложно адаптировать под нее методику, характерную для Tree Code и Multipole Expansion. Так, для каждого элемента дискретного поля:

    1. Непосредственно учитываем влияние масс только соседних восьми элементов;
    2. Учитываем влияние восьми соседних регионов, состоящих из девяти элементов (3×3), путем вычисления суммы их масс с помощью интегрального изображения;
    3. На каждом последующем шаге увеличиваем размер региона в 3 раза (9×9, 27×27, 81×81, и т.д.) до тех пор, пока этот он не превысит размер всего векторного поля.


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

    Сложность приближенного решения гравитационной задачи для поля импульсов с помощью интегрального изображения растет квазилинейно, как O(N*8*log3(sqrt(N))) для R = 2 и как O(N*26*log3(cbrt(N))) для R = 3.



    Однако в таком виде это решение обладает тем же недостатком, что и методика Multigrid, а именно ― ощутимой дискретной «ригидностью» низкочастотных компонентов функции гравитации. Проще всего продемонстрировать эту проблему наглядно.


    Силы на этой иллюстрации закодированы цветом в формате HSV (Hue ― направление силы, Value линейно отображает порядок его абсолютной величины таким образом, что наименее различимы (Value = 0,05) силы порядка 10-7, а самые яркие силы (Value = 1,0) имеют порядок 103 и выше).

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


    В центре иллюстрации показан экстремум масс при сравнительно равномерном распределении массы в остальном поле

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

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

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


    Интегральное изображение координат

    Возьмем для примера регион в левом нижнем углу с координатами {3; 1} и в верхнем правом с координатами {7; 5}. Сумма координат этого региона {168; 120} ― {18; 45} ― {28; 0} + {3; 0} равна {125; 75}, общее число элементов равно 25, следовательно, координаты его центра будут {5; 3}.

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


    Интегральное изображение взвешенных координат. Увеличенным шрифтом обозначены весовые коэффициенты

    Рассмотрим регион в левом нижнем углу с координатами {2; 1} и в правом верхнем с координатами {5; 4}. Сумма весовых коэффициентов этого региона 4,8 ― 1,0 ― 0,6 + 0,2 равна 3,4, а сумма его взвешенных координат {11,1; 13,2} ― {0,5; 2,0} ― {1,5; 0,0} + {0,1; 0,0} равна {9,2; 11,2}. Таким образом, взвешенные координаты центра этого региона равны {2,7; 3,3}.

    Убедиться в том, что эта схема действительно работает, можно наглядно ― преобразовав интегральное изображение со взвешенными координатами в distance field с помощью distance transform[19].


    Преобразование интегрального изображения со взвешенными координатами в distance field. Слева ― изображение масс. Следующее изображение ― дистанция до центра масс для региона размером 3×3 элемента. Далее ― дистанция до центров масс для регионов размерами 9×9 и 27×27 элементов. Величины дистанций на этой иллюстрации нормализуются к размеру региона, используемого для выборки.


    Преобразование интегрального изображения со взвешенными координатами в directed distance field. Направление и дистанция закодированы цветом в формате HSV, где Hue ― направление, Value ― нормализованная дистанция. Величины дистанций на этой иллюстрации нормализуются к размеру региона, используемого для выборки.

    Суммируя все вышеизложенное:

    1. К исходному изображению масс, используемому для вычисления интегрального изображения, добавляется еще два (или три, если R = 3) компонента ― взвешенные координаты элемента, равные координатам, умноженным на массу.
    2. Координаты центра масс произвольного региона вычисляется путем деления суммы его взвешенных координат на сумму его масс.
    3. Вычисление центра масс региона имеет константную вычислительную сложность O(1).
    4. Сложность по памяти для интегрального изображения со взвешенными координатами остается линейной O(M*N).


    Сравнение методов аппроксимации функции гравитации. Сверху ― для аппроксимации используется интегральное изображение без взвешенных координат. Снизу ― для аппроксимации используется интегральное изображение со взвешенными координатами.

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


    Время обратить внимание на одну из ключевых особенностей интегральных изображений, которая до сих пор ограничивает область их применения, а именно ― ресурсоемкость и численную устойчивость. Так, в зависимости от диапазона величин, которые содержит оригинальное изображение, для его интегрального изображения может потребоваться более «длинный» тип данных. К примеру, при расчете среднеквадратичных отклонений кроме оригинального изображения (диапазон значений 0..255) используется изображение, содержащее квадраты соответствующих величин (диапазон значений 0..65535). В этом случае для вычисления интегральных изображений большого размера 32 бит точности оказывается недостаточно[20].

    Похожая ситуация наблюдается и в случае интегральных изображений со взвешенными координатами. Сама по себе величина суммы координат, которую необходимо хранить в интегральном изображении, растет в зависимости от размера изображения N: квадратично для одномерного случая 0 + 1 + 2 +… + N ― 1 = (N ― 1)*N/2 и кубически для двумерного N*(N ― 1)*N/2. Например, чтобы хранить сумму одной целочисленной координаты для изображения размером 2048×2048 (равную 4292870144), едва хватит 32-битного беззнакового целого (максимальное значение которого 4294967295). При вычислении интегральных изображений большего размера наступит переполнение.

    Использование в интегральных изображениях 32-битных чисел с плавающей точкой отдаляет проблему переполнения на астрономическое расстояние (1038 ― 1010), но при этом точность вычислений со взвешенными координатами значительно снижается из-за накапливаемой при суммировании ошибки усечения[21].


    Абсолютная величина ошибки при вычислении взвешенного центра масс региона размером 4×4 элемента. Интегральное изображение содержит числа с одинарной точностью (32 бита).


    Абсолютная величина ошибки при вычислении взвешенного центра масс региона размером 32×32 элемента. Интегральное изображение содержит числа с одинарной точностью (32 бита).

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


    Увеличение диапазона значений весовых коэффициентов не влияет на абсолютную величину ошибки при расчете взвешенного центра масс. На графике приведены погрешности для «худшего» региона (в верхнем правом углу интегрального изображения). Размер интегрального изображения ― 256×256 элементов.


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

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

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

    Заключение


    Приведенный выше пример использования интегрального изображения, вероятно, можно адаптировать для разработки оптимального алгоритма O(1) выборки по значимости (importance sampling). Существующие ― и крайне ресурсоемкие с точки зрения GPU ― алгоритмы имеют линейную сложность O(N) и находят применение в современных методах глобального освещения[22, 23].

    В целом, на мой взгляд, интегральные изображения ― один из самых недооцененных алгоритмов. Они могут составить великолепную альтернативу одновременно мипмаппингу и анизотропной фильтрации, а с учетом того, как именно реализовано последнее на GPU[24, 25], вполне возможно, что являются уже и более эффективной мерой.

    Ссылки


    Pixonic
    Разрабатываем и издаем игры с 2009 года

    Comments 22

      –4
      Приведенный выше пример использования интегрального изображения, вероятно, можно адаптировать для разработки оптимального алгоритма O(1)
      Так надо было в заголовке так и писать. Нахрена мне читать гору расчётов если это не 0(1). Любая студентота напишет алгоритм с О(n) для центра масс
        +9
        Shiny2, прочитайте еще раз, внимательнее, начиная от заголовка, и заканчивая заключением. Вы удивительно торопливы.
        +4
        Так хотелось увидеть анимацию после всех этих красочных картинок
        +1

        Интересная статья и идея. У вас уже есть реализация? Было бы очень интересно увидеть результы сравнения с другими реализациями.


        В целом, на мой взгляд, интегральные изображения ― один из самых недооцененных алгоритмов.

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


        Ещё интересен вопрос точности. Мне в аналогичной ситуации, при суммировании даже точности double не хватило. Пришлось алгоритм Кохена использовать.

          +1
          Наверное, мне просто нравится интуитивно понятный алгоритм, который превращает сложность O(M*N) в O(1), и вероятно может быть адаптирован для многих вычислительных задач, где нужно складывать и умножать по регионам.
            0
            В любом случае, Fast Multipole Method никакая моя реализация не побьет :) Но в отличие от FFT, которая решает только задачу N тел и более ничего, у интегральных изображений потенциал все же больше.
              +2

              FFT используется в огромном количестве задач. Или вы хотели сказать FMM?


              А вот по скорости надо пробовать, тут может быть очень неоднозначно.

                0
                FMM, конечно, а не FFT :)
            0

            Я несколько озадачен. Кажется, однократный проход по матрице MxN — это О(MxN), а не О(1), как в статье. Или я что-то упустил?

              0

              Для вычисления центра масс по интегральному изображению сложность O(1), но, конечно, нужно ещё это изображение построить.

                +1
                И да и нет. Для двоичного поиска тоже нужно отсортировать массив. А это O(NlogN) в лучшем случае. Но двоичный поиск при этом остается O(logN), не так ли?
                  0

                  Теперь я понял идею, спасибо :)

                +2
                Надеюсь вы знаете о такой фишке как вычитание среднего значения из данных перед интегрированием?
                Потом при получении суммы — просто прибавляете среднее значение необходимое количество раз. Такой трюк особенно полезен для картинок, где значения обычно от 0 до 255 и интегральные значения растут до невообразимых высот :) Если эе сместить range пикселей к интервалу -128...127 интегральное значение будет бегать около нуля, что в случае floating-point значений сильно улучшает точность.

                Но вообще говоря, в image processing community, интегральные изображения давно отошли на задний план. С недавним прогрессом в Convolutional Neural Networks все эти ручные оптимизации и алгоритмические потуги постепенно сходят на нет. Самая тупая нейро-сетка справляется с задачами гораздо лучше чем алгоритмы, разрабатываемые десятками лет. Я думаю и вам нужно смотреть в том-же направлении.
                  +1
                  С недавним прогрессом в Convolutional Neural Networks все эти ручные оптимизации и алгоритмические потуги постепенно сходят на нет.

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

                    0
                    Любопытно. Но у меня как правило получается изображение со сравнительно однородным распределением и несколькими областями-экстремумами, значения которых отличаются на несколько порядков в большую сторону.
                    +1

                    Спасибо за статью, было интересно почитать про подход.


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

                      0
                      Спасибо! Попробую сгенерировать иллюстрацию со строгим решением, хотя не уверен, что GPU это решение осилит. Может быть для поля небольшого размера и получится.
                      +1
                      Кажется в других науках ваши интегральный изображение называются кумулятивной функцией распределения.
                      Я понял, что мы переходом в центр масс мы поднимаем точность до квадратичной, так как убиваем дипольный момент. А не знаете какого-нибудь способа хоть как-нибудь посчитать квадрупольный тензор, чтобы поднять точность до куба? Понятно, что в лоб не вариант, там сложность будет квадрат от числа элементов в секторе. Но может есть какие-нибудь фокусы?
                        0
                        Верно, cumulative distribution function. На второй вопрос я к сожалению не смогу ответить, моя математическая подготовка будет несколько пониже уровнем :(
                        0
                        Но на каждой итерации объекты меняют положения и скорости, а значит и интеральное изображение меняется. Как вы его обновляете, чтобы не тратить O(M*N) на каждом шаге по времени?
                          0
                          Обновлять приходится на каждом шаге, но вычисление интегрального изображения — это самый простой этап по абсолютным вычислительным затратам. Вообще, целиком алгоритм эволюции поля примерно такой: 1) вычисление интегрального изображения масс O(N), 2) вычисление сил — O(NlogN), 3) решение дифференциальных уравнений движения O(N) и 4) перераспределение масс O(N). Здесь N = X*Y — Общее количество дискретных элементов в поле.

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

                          Решение дифференциальных уравнений — это алгоритм RK4, и в абсолютном отношении оно примерно вчетверо сложнее вычисления интегрального изображения.

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