Одномерный поиск образца с использованием дискретного преобразования Фурье

После прочтения статьи про поиск изображения в изображении[1], осталось множество вопросов к формулам, да и к самому коду где преобразование массивов мне показалось не прозрачным из-за использования множества сторонних библиотечных функций.

Потому занялся дополнительным поиском готовых реализаций, но к сожалению не смотря на обилие упоминаний идеи 1974 года[2], реализаций алгоритма, даже на законодателе моды в вычислительной математике Фортране я не обнаружил. В семинарах и лекциях да и в диссертациях описание не блистало целостностью, потому собрав с десяток статей и обсуждений в кучу появилось желание написать статью для тех кто простейшую реализацию поиска подстроки хочет "подержать в руках".

И так, написание алгоритмических программ обычно произвожу сначала в математических пакетах, и только после выяснения численной устойчивости и корректности работы алгоритма перевожу в код на компилируемые или байт ориентированные языки. Таков уж мой опыт, - или считать медленно но точно, или быстро но то что уже практически известно. Потому для отладочного иллюстративного кода использовал Wolfram Language и набор функций пакета Mathematica V 12.

Собственно в чем ценность идеи: использование дискретного Фурье преобразования (DFT) сокращает сложность вычисления от "наивного" O(n*m) до O(n Log(n)), где n - размер текста а m - размер искомого образца. Более того расширения метода позволяют производить поиск с "Джокером", - символом обозначающим любой другой символ в алфавите, в то время как суффиксные реализации это сделать не позволяют.

Описание "наивного" подхода:

Для массива T размером n и образца P размером m, вычислим функцию квадратов разницы значений элементов. Нумерация в массиве начинается с нуля.

S_i^F = \sum\nolimits_{j = 0}^{m - 1} {({t_{i + j}}}  - {p_j}{)^2} = \sum\nolimits_{j = 0}^{m - 1} {t_{i + j}^2}  - 2\sum\nolimits_{j = 0}^{m - 1} {{t_{i + j}}} {p_j} + \sum\nolimits_{j = 0}^{m - 1} {p_j^2}  = T{T_i} - 2P{T_i} +P{P_i}

Очевидно что в точке соответствия искомая сумма показывает минимум, если точнее обнуляется. После раскрытия квадрата под суммой получаются три слагаемых, последнее из которых имеет постоянное значение. Соответственно для поиска минимума необходимо вычислить только первые две суммы. Можно увидеть что прямое вычисление всех элементов S требует O((n-m+1)*m) операций, или оценочно O(n*m).

В качестве массива используем строчку из картинки на половине ее высоты, а образец который будем искать пусть будет отрезком в первой половине этой строчки:

"Test.png"
"Test.png"

Произведем прямое вычисление искомой функции:

{S_i} = \sum\nolimits_{j = 0}^{m - 1} {t_{i + j}^2}  - 2\sum\nolimits_{j = 0}^{m - 1} {{t_{i + j}}} {p_j} = T{T_i} - 2P{T_i}
Img = ColorConvert[Import["Test.png"], "Grayscale"];
{W, H} = ImageDimensions[Img];   

y = IntegerPart[H/2];                                (*Pattern width and coordinates*)
x = IntegerPart[W/4]; 
w = IntegerPart[W/8];            

L = Part[ImageData[ImageTake[Img, {y, y}]],1];       (*Stripe Array*)

T[i_] := Table[Part[L[[i + j - 1]], 1], {j, 1, w}] ;      (*Sample data*)
P[i_] := Table[Part[L[[j + x - 1]], 1], {j, 1, w}] ;      (*Pattern data*)

TT = Table[Sum[(T[i][[j]]* T[i][[j]]), {j, 1, w}], {i, 1, W - w + 1}];
PT = Table[Sum[(P[i][[j]]* T[i][[j]]), {j, 1, w}], {i, 1, W - w + 1}];

ListPlot[TT - 2 PT, Joined -> True, PlotRange -> All]
Результат вычисления квадрата разницы без постоянного слагаемого
Результат вычисления квадрата разницы без постоянного слагаемого

Как видно в точке (x=175), где был взят образец, функция показала минимальное значение и повторила его значение ведь изображение почти дублируется.

Свойства дискретного Фурье преобразования.

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

Вычисление PT

PolyT = {1, 2, 3, 4, 5};               LT = Length[PolyT];
PolyP = {1, 2, 3};                     LP = Length[PolyP];
PolyR = {};                            LR = LT + LP - 1;

eT = Table[If[i > LT, 0, PolyT[[i]]], {i, 1, LR}]
eP = Table[If[i > LP, 0, PolyP[[i]]], {i, 1, LR}]

PolyR = InverseFourier[
   Fourier[eT, FourierParameters -> {1, 1}]*
   Fourier[eP, FourierParameters -> {1, 1}]
  , FourierParameters -> {1, 1}]
PolyR[[LP ;; LT]]

результат действия такого кода:

{1, 2, 3, 4, 5, 0, 0} 						(* PolyT *)
{1, 2, 3, 0, 0, 0, 0} 						(* PolyP *)
{1., 4., 10., 16., 22., 22., 15.} (* PolyR = PolyT * PolyP *)
{10., 16., 22.}                   (* PolyR starting at LP to LT*)	

Итак, если массив значений представить полиномами, то получаемое значение тоже полином размером n+m-1.

\left( {1 + 2x + 3{x^2} + 4{x^3} + 5{x^4}} \right)\left( {1 + 2x + 3{x^2}} \right) = 1 + 4x + 10{x^2} + 16{x^3} + 22{x^4} + 22{x^5} + 15{x^6}

Более того, начиная с позиции m (если нумерация начинается с единицы) мы получаем сумму перекрестного произведения элементов для окна длинны m:

10 = 1*3+2*2+3*1
16 = 2*3+3*2+4*1
...

Потому для вычисления элементов PT искомый образец P разворачивается. Всего получается n-m+1 значений.

Вычисление TT

PolyT = {1, 2, 3, 4, 5};            LT = Length[PolyT];
PolyP = {1, 1, 1};                  LP = Length[PolyP];
PolyR = {};                         LR = LT + LP - 1;

eT = Table[If[i > LT, 0, PolyT[[i]]], {i, 1, LR}]
eP = Table[If[i > LP, 0, PolyP[[i]]], {i, 1, LR}]

PolyR = InverseFourier[
   Fourier[eT, FourierParameters -> {1, 1}]*
   Fourier[eP, FourierParameters -> {1, 1}]
  , FourierParameters -> {1, 1}]
PolyR[[LP ;; LT]]

результат действия кода:

{1, 2, 3, 4, 5, 0, 0}                 (* PolyT *)
{1, 1, 1, 0, 0, 0, 0}                 (* PolyP *)
{1., 3., 6., 9., 12., 9., 5.}         (* PolyR = PolyT * PolyP *)
{6., 9., 12.}                         (* PolyR starting at LP to LT*)	
6 = 1*1+2*1+3*1
9 = 2*1+3*1+4*1
...

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

Сборка и сравнение

Вычисление PP и TT с использованием DFT:

Tt = Table[If[1 <= i <= W, Part[L[[i]], 1], 0], {i, 1, Wt}] ;    (*Sample data*)
Ft = Fourier[Tt, FourierParameters -> {1, 1}];

Ts = Table[If[1 <= i <= W, (Part[L[[i]], 1])^2, 0], {i, 1, Wt}]; (*Sample squared data*)
Fs = Fourier[Ts, FourierParameters -> {1, 1}];

Es = Table[If[1 <= i <= w, 1, 0], {i, 1, Wt}] ;                  (*m size unit array*)
Fe = Fourier[Es, FourierParameters -> {1, 1}];

Pa = Table[If[1 <= i <= w, Part[L[[x + w - i]], 1], 0], {i, 1, Wt}];                                                             \
Fp = Fourier[Pa, FourierParameters -> {1, 1}];                   (*Inverse pattern data*)

TTf = InverseFourier[Fs*Fe, FourierParameters -> {1, 1}][[w ;; W]];
PTf = InverseFourier[Ft*Fp, FourierParameters -> {1, 1}][[w ;; W]];

Сравниваем полученные значения:

ListPlot[{TT - 2 PT, TTf - 2 PTf, TT - 2 PT - TTf + 2 PTf}, Joined -> True, PlotRange -> All]
Три графика, два совпадающих и один показывающий разницу между ними, совпадает с осью.
Три графика, два совпадающих и один показывающий разницу между ними, совпадает с осью.

Выводы

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

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

  1. https://habr.com/ru/post/266129/

  2. https://eduardgorbunov.github.io/assets/files/amc_778_seminar_08.pdf

Комментарии 24

    +1
    Верно ли, что ваш алгоритм вычисляет взаимнокорреляционную последовательность? И поэтому, по теореме о свёртке, такую взаимную корреляцию можно вычислять с помощью преобразования Фурье?
      0
      Да алгоритм как промежуточный результат для вычисления PT это использует, вычисление TT тоже имеет свое математическое название хоть и является частным случаем первого.

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

        0
        Если вы обратите внимание, то в определении свертки и взаимной корреляционной последовательности есть различие, в интегральной форме смещение незаметно так как пределы интегрирования бесконечны, но для понимания алгоритма нахождения искомой суммы определение свертки менее подходящее. Потому я тоже избегал упоминания теоремы о свертке так как это требует пусть не значительной но замены переменных и дополнительных преобразований.
          +1
          Если вы обратите внимание, то в определении свертки и взаимной корреляционной последовательности есть различие
          Единственное различие между свёрткой и корреляцией — это знак функции времени второго аргумента, а так они вполне однозначно выражаются через друг друга.

          в интегральной форме смещение незаметно так как пределы интегрирования бесконечны
          Это очень спорное утверждение, сможете его подтвердить реальным примером, или я что-то не так понял? Потому что смещение во времени — это просто линейный наклон фаз вне зависимости от того, дискретное или непрерывное преобразование используется.
            0
            сможете его подтвердить реальным примером,

            разумеется:
            — использование инверсии и взаимнокорреляционной функции приводят началу последовательности PT с позиции m
            использование свертки и сопряжения, приводят началу последовательности PT с позиции 1

            Такое смещение в точности соответствует преобразованию для изменения «знака функции времени второго аргумента».

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

                Теперь про интегрирование(приведу без комплексного варианта),
                — в сверточном определении f(x)g(y-x)dx
                — для взаимнокорреляционной функции f(x+y)g(x)dx
                переход от одного к другому это замена переменных
                x=x+y — вот откуда смещение.

                Если вы будете искать t[i+j]p[j] в произведении полиномов
                A×B=InvDFT(DFT(A)×DFT(B))
                то получите начальную позицию m.

                Я использую именно свойство DFT так как оно определено аналитически точно, но в пределах машинной точности вычислений.

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

                Дополнительные действия есть как при использовании DFT, так и при FT, но для первого их меньше. В названии статьи приведено DFT и именно его свойствами я и пользовался в создании алгоритма, хотя можно использовать и свойства FT

                  +2
                  Что-то я не пойму, о чём спор. Есть же дискретный вариант теоремы о свёртке, там всё прекрасно выражается через ДПФ. А знак времени для одной из функций (то, чем различается свёртка от автокорреляции) — это уже мелочи.
                    0
                    Спор о выборе теории для иллюстрации действия алгоритма. Автокорреляция это еще один термин который от темы уводит дальше чем свертка.
                      0
                      Дело в разных размерах массивов, на практике можно получить нужную последовательность несколькими способами, но они отличаются разными способами подготовки массивов и позициями в них. Что говорить на остаток и саму позицию здесь и разбираем.

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

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

          Этот метод легко обобщается до 2+ измерений, просто решили не тратить место в статье?
          Существует еще такая вариация, как БПФ на кольцах по модулю простых чисел, она точная. (Хотя, вариант с плавающей точкой зачастую все равно выигрывает по времени, проигрывая по используемой памяти).

            0
            Относительно утверждения про легко, я сомневаюсь, и именно поиску изображения в изображении планирую создать еще одну публикацию. Насколько видно из [1] там автор с легкостью то и не справился, забыв указать про инверсию массива образца, оставив только не существенное комплексное сопряжение которое к стати в приведенном одномерном алгоритме исключено. Как и используется другая функция для определения позиции образца.
              0

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


              Для целочисленного вычисления сверки по вместо длинной арифметики есть еще такой вариант, как сделать операцию по модулю нескольких различных простых, а потом использовать как-то вернуть обратно (https://scask.ru/g_book_s_alg.php?id=20)

                0
                Да актуально, но приведенная статья посвящена больше внутренней работе функций преобразования, я сразу перешел к следующей главе, но и там не заметил упоминания или построения функции сравнения.
              0
              «Существует еще такая вариация, как БПФ на кольцах по модулю простых чисел, она точная.» аналитически может и точная, но при использовании численного расчета на действительный числах обойти минимум машинную точность вычисления не получится. Поэтому пишу про приближение и возможный диапазон яркостей.
                +1
                БПФ на кольцах вычисляется с применением исключительно целочисленных операций. Результат тоже целочисленный. Поэтому оно точное. В английской терминологии обозначается как «Number-theoretic transform».
                  0
                  Мне хватило точности и без использования целочисленного словаря. Но спасибо за точное название метода, может кому то пригодится.
                0
                Если использовать сопряжение и прямой порядок, тогда нумерация PT начинается с первой позиции, длинна же как и в публикации будет n-m+1
                (*PT*)
                PolyT = {1, 2, 3, 4, 5}; LT = Length[PolyT];
                PolyP = {3, 2, 1}; LP = Length[PolyP];
                PolyR = {}; LR = LT + LP - 1;

                eT = Table[If[i > LT, 0, PolyT[[i]]], {i, 1, LR}]
                eP = Table[If[i > LP, 0, PolyP[[i]]], {i, 1, LR}]

                PolyR = InverseFourier[
                Fourier[eT, FourierParameters -> {1, 1}]*
                Conjugate[Fourier[eP, FourierParameters -> {1, 1}]]
                , FourierParameters -> {1, 1}]


                Результат:
                {1, 2, 3, 4, 5, 0, 0}
                {3, 2, 1, 0, 0, 0, 0}
                {10., 16., 22., 22., 15., 1., 4.}
                  0
                  Причина в разнице начальной позиции в использовании различных функций в определении. При теории основанной на свойствах полиномов необходима инверсия, тогда как использование теоремы о свертке подразумевает дополнительную замену переменных, которая и приводит к величине изменения позиции.

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

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

                  При использовании дискретного преобразования Фурье, матрица M содержит также элементы от циклического сдвига изображений между собой. Поэтому, если не требуется анализировать циклический сдвиг кадров, то поиск максимального элемента в матрице M нужно ограничить областью (0,0)-(N1-M1, N2-M2).


                  Именно из-за таких мелочей которые будучи смешаны в кучу дают «абракадабру» я и написал статью.
                    0
                    Если у вас под рукой есть ссылка на теорию по построению функции сравнения для варианта 2D напишите ее пожалуйста. Вариант TT-2PT очевиден именно для одномерного случая, тогда как для двумерного случая уже не понятно, ведется суммирование по строке или столбцу, а если и то и другое то как это соотносится с точечной разницей яркости.
                      0
                      В своей практике я в основном работал с рядами и интегралами Фурье, они удобны при подстановке в дифференциальные уравнения тем что превращают их в систему линейных уравнений, потому ожидал что обрезка ряда по частотам приведет к некоторым ошибкам которые можно будет уточнить другим способом.

                      Однако как оказалось несмотря на близость названий, DFT имеет другой набор свойств хотя и совместим со свойствами обычного преобразования. Само определение DFT для полинома это его значение в точке x=Exp(2*i*Pi/n), — точное, как и формула преобразования, которая вроде похожа на сверточную версию но все же без сопряжения. Соответственно не имеет ошибки связанной с конечностью ряда, а содержит только ошибку связанную точностью машинного вычисления. Которую, как вы упомянули, возможно устранить при сведении вычислений к целочисленным преобразованиям.

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

                      В итоге как оказалось вроде бы очевидное и легкое дополнение для каждой лекции и статьи которые просматривал, в реализации дюжины строчек кода, показали страшную ретивость и потребовали уйму времени для их отладки и согласования.
                      0
                      Не понял что с чем сравниваете, точнее что ищите?
                        0
                        Спасибо за замечание, тоже думаю внести правку в статью, так как это написано только в коде. Ищется отрезок длинной W/8 с началом в точке W/4 на линии рисунка с вертикальной координатой H/2.
                        где W-ширина (width), H-высота(height) рисунка в пикселях.

                        ...
                        y = IntegerPart[H/2]; (*Pattern width and coordinates*)
                        x = IntegerPart[W/4];
                        w = IntegerPart[W/8];
                        0
                        Задача которая привела к написанию статьи была связана с нахождением величины смещения двух изображений, меня не устроила скорость «наивной» реализации. Потому такой развлекательный вариант исходного массива.

                        Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                        Самое читаемое