One-Dimensional Sample Search Using Discrete Fourier Transform

After reading the article about searching for an image in an image [1], there are many questions left to the formulas, and to the code itself, where the transformation of arrays seemed to me not transparent due to the use of many third-party library functions.





Therefore, I took up an additional search for ready-made implementations, but unfortunately, despite the abundance of references to the idea of ​​1974 [2], I did not find implementations of the algorithm even on the trendsetter in computational mathematics Fortran. In seminars and lectures, and even in dissertations, the description did not shine with integrity, therefore, having collected a dozen articles and discussions in a heap, there was a desire to write an article for those who want to "hold in their hands" the simplest implementation of substring search.





And so, I usually write algorithmic programs first in mathematical packages, and only after finding out the numerical stability and correctness of the algorithm's operation I translate it into code into compiled or byte-oriented languages. This is my experience - either count slowly but accurately, or quickly, but what is already practically known. Therefore, for the debugging illustrative code I used the Wolfram Language and the set of functions of the Mathematica V 12 package.





Actually, what is the value of the idea: the use of discrete Fourier transform (DFT) reduces the complexity of the calculation from "naive" O (n * m) to O (n Log (n)), where n is the size of the text and m is the size of the desired sample. Moreover, method extensions allow you to search with "Joker" - a symbol denoting any other character in the alphabet, while suffix implementations do not allow this.





Description of the "naive" approach:





For an array T of size n and a sample P of size m, we calculate the function of the squares of the difference between the values ​​of the elements. The array is numbered starting at zero.





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]
      
      



The result of calculating the squared difference without a constant term

(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.









Calculating PP and TT using 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]];
      
      



Let's compare the obtained values:





ListPlot[{TT - 2 PT, TTf - 2 PTf, TT - 2 PT - TTf + 2 PTf}, Joined -> True, PlotRange -> All]
      
      



Three graphs, two coinciding and one showing the difference between them, coincides with the axis.
Three graphs, two coinciding and one showing the difference between them, coincides with the axis.

conclusions





Despite the fact that the method is approximate, its accuracy is more than enough for working with text and most ordinary images where the brightness values ​​differ significantly.





The given code does not pretend to be optimized for performance, but is intended more for the convenience of understanding and evaluating the accuracy of the algorithm.





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





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








All Articles