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

После прочтения статьи про поиск изображения в изображении[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»

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

{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

© Habrahabr.ru