[ZOJ][FFT快速傅立葉] 1637 - Fast Image Match
For example, we have
A (3 * 3): | a1 | a2 | a3 | B (2 * 2): | b1 | b2 |
a4 | a5 | a6 | b4 | b5 | ||
a7 | a8 | a9 |
When B is placed on position a5, the difference of them is ((b1-a5)^2 + (b2-a6)^2 + (b4-a8)^2 + (b5-a9)^2). Now we hope to have the position of the top left corner of B that gives the minimum difference. (B must completely reside on A)
It is clear that a simple solution will appear with very low efficiency when A and B have too many elements. But we can use 1-dimensional repeat convolution, which can be computed by Fast Fourier Transform (FFT), to improve the performance.
A program with explanation of FFT is given below:
/** * Given two sequences {a1, a2, a3.. an} and {b1, b2, b3... bn}, * their repeat convolution means: * r1 = a1*b1 + a2*b2 + a3*b3 + ... + an*bn * r2 = a1*bn + a2*b1 + a3*b2 + ... + an*bn-1 * r3 = a1*bn-1 + a2*bn + a3*b1 + ... + an*bn-2 * ... * rn = a1*b2 + a2*b3 + a3*b4 + ... + an-1*bn + an*b1 * Notice n >= 2 and n must be power of 2. */ #include <vector> #include <complex> #include <cmath> #define for if (0); else for using namespace std; const int MaxFastBits = 16; int **gFFTBitTable = 0; int NumberOfBitsNeeded(int PowerOfTwo) { for (int i = 0;; ++i) { if (PowerOfTwo & (1 << i)) { return i; } } } int ReverseBits(int index, int NumBits) { int ret = 0; for (int i = 0; i < NumBits; ++i, index >>= 1) { ret = (ret << 1) | (index & 1); } return ret; } void InitFFT() { gFFTBitTable = new int *[MaxFastBits]; for (int i = 1, length = 2; i <= MaxFastBits; ++i, length <<= 1) { gFFTBitTable[i - 1] = new int[length]; for (int j = 0; j < length; ++j) { gFFTBitTable[i - 1][j] = ReverseBits(j, i); } } } inline int FastReverseBits(int i, int NumBits) { return NumBits <= MaxFastBits ? gFFTBitTable[NumBits - 1][i] : ReverseBits(i, NumBits); } void FFT(bool InverseTransform, vector<complex<double> >& In, vector<complex<double> >& Out) { if (!gFFTBitTable) { InitFFT(); } // simultaneous data copy and bit-reversal ordering into outputs int NumSamples = In.size(); int NumBits = NumberOfBitsNeeded(NumSamples); for (int i = 0; i < NumSamples; ++i) { Out[FastReverseBits(i, NumBits)] = In[i]; } // the FFT process double angle_numerator = acos(-1.) * (InverseTransform ? -2 : 2); for (int BlockEnd = 1, BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) { double delta_angle = angle_numerator / BlockSize; double sin1 = sin(-delta_angle); double cos1 = cos(-delta_angle); double sin2 = sin(-delta_angle * 2); double cos2 = cos(-delta_angle * 2); for (int i = 0; i < NumSamples; i += BlockSize) { complex<double> a1(cos1, sin1), a2(cos2, sin2); for (int j = i, n = 0; n < BlockEnd; ++j, ++n) { complex<double> a0(2 * cos1 * a1.real() - a2.real(), 2 * cos1 * a1.imag() - a2.imag()); a2 = a1; a1 = a0; complex<double> a = a0 * Out[j + BlockEnd]; Out[j + BlockEnd] = Out[j] - a; Out[j] += a; } } BlockEnd = BlockSize; } // normalize if inverse transform if (InverseTransform) { for (int i = 0; i < NumSamples; ++i) { Out[i] /= NumSamples; } } } vector<double> convolution(vector<double> a, vector<double> b) { int n = a.size(); vector<complex<double> > s(n), d1(n), d2(n), y(n); for (int i = 0; i < n; ++i) { s[i] = complex<double>(a[i], 0); } FFT(false, s, d1); s[0] = complex<double>(b[0], 0); for (int i = 1; i < n; ++i) { s[i] = complex<double>(b[n - i], 0); } FFT(false, s, d2); for (int i = 0; i < n; ++i) { y[i] = d1[i] * d2[i]; } FFT(true, y, s); vector<double> ret(n); for (int i = 0; i < n; ++i) { ret[i] = s[i].real(); } return ret; } int main() { double a[4] = {1, 2, 3, 4}, b[4] = {1, 2, 3, 4}; vector<double> r = convolution(vector<double>(a, a + 4), vector<double>(b, b + 4)); // r[0] = 30 (1*1 + 2*2 + 3*3 + 4*4) // r[1] = 24 (1*4 + 2*1 + 3*2 + 4*3) // r[2] = 22 (1*3 + 2*4 + 3*1 + 4*2) // r[3] = 24 (1*2 + 2*3 + 3*4 + 4*1) return 0; }
Input
The first line contains n (1 <= n <= 10), the number of test cases.
For each test case, the first line contains four integers m, n, p and q, where A is a matrix of m * n, B is a matrix of p * q (2 <= m, n, p, q <= 500, m >= p, n >= q). The following m lines are the elements of A and p lines are the elements of B.
Output
For each case, print the position that gives the minimum difference (the top left corner of A is (1, 1)). You can assume that each test case has a unique solution.
Sample Input
2 2 2 2 2 1 2 3 4 2 3 1 4 3 3 2 2 0 5 5 0 5 5 0 0 0 5 5 5 5
Sample Output
1 1 1 2
兩張圖片的快速匹配,找差平方和最小的匹配位置。
// sigma((x[i][j]-y[p][q])^2)
brute force => O(n*m*p*q) 相當於兩個總像素個數積。
// FFT O(N *log N) // N = n*m
使用 convolution 旋積計算。
-----
不知道為什麼 AC 而且還是最快的,已經兩年沒人寫這一題。
(*´∀`)っ゛FFT 是什麼,問神去吧。
(*´∀`)っ゛題目提供的 code 加些計算就提交了。
#include <cstdio>
#include <cstring>
#include <vector>
#include <complex>
#include <cmath>
#define for if (0); else for
using namespace std;
const int MaxFastBits = 16;
int **gFFTBitTable = 0;
int NumberOfBitsNeeded(int PowerOfTwo) {
for (int i = 0;; ++i) {
if (PowerOfTwo & (1 << i)) {
return i;
}
}
}
int ReverseBits(int index, int NumBits) {
int ret = 0;
for (int i = 0; i < NumBits; ++i, index >>= 1) {
ret = (ret << 1) | (index & 1);
}
return ret;
}
void InitFFT() {
gFFTBitTable = new int *[MaxFastBits];
for (int i = 1, length = 2; i <= MaxFastBits; ++i, length <<= 1) {
gFFTBitTable[i - 1] = new int[length];
for (int j = 0; j < length; ++j) {
gFFTBitTable[i - 1][j] = ReverseBits(j, i);
}
}
}
inline int FastReverseBits(int i, int NumBits) {
return NumBits <= MaxFastBits ? gFFTBitTable[NumBits - 1][i] : ReverseBits(i, NumBits);
}
void FFT(bool InverseTransform, vector<complex<double> >& In, vector<complex<double> >& Out) {
if (!gFFTBitTable) { InitFFT(); }
// simultaneous data copy and bit-reversal ordering into outputs
int NumSamples = In.size();
int NumBits = NumberOfBitsNeeded(NumSamples);
for (int i = 0; i < NumSamples; ++i) {
Out[FastReverseBits(i, NumBits)] = In[i];
}
// the FFT process
double angle_numerator = acos(-1.) * (InverseTransform ? -2 : 2);
for (int BlockEnd = 1, BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) {
double delta_angle = angle_numerator / BlockSize;
double sin1 = sin(-delta_angle);
double cos1 = cos(-delta_angle);
double sin2 = sin(-delta_angle * 2);
double cos2 = cos(-delta_angle * 2);
for (int i = 0; i < NumSamples; i += BlockSize) {
complex<double> a1(cos1, sin1), a2(cos2, sin2);
for (int j = i, n = 0; n < BlockEnd; ++j, ++n) {
complex<double> a0(2 * cos1 * a1.real() - a2.real(), 2 * cos1 * a1.imag() - a2.imag());
a2 = a1;
a1 = a0;
complex<double> a = a0 * Out[j + BlockEnd];
Out[j + BlockEnd] = Out[j] - a;
Out[j] += a;
}
}
BlockEnd = BlockSize;
}
// normalize if inverse transform
if (InverseTransform) {
for (int i = 0; i < NumSamples; ++i) {
Out[i] /= NumSamples;
}
}
}
vector<double> convolution(vector<double> a, vector<double> b) {
int n = a.size();
vector<complex<double> > s(n), d1(n), d2(n), y(n);
for (int i = 0; i < n; ++i) {
s[i] = complex<double>(a[i], 0);
}
FFT(false, s, d1);
s[0] = complex<double>(b[0], 0);
for (int i = 1; i < n; ++i) {
s[i] = complex<double>(b[n - i], 0);
}
FFT(false, s, d2);
for (int i = 0; i < n; ++i) {
y[i] = d1[i] * d2[i];
}
FFT(true, y, s);
vector<double> ret(n);
for (int i = 0; i < n; ++i) {
ret[i] = s[i].real();
}
return ret;
}
double a[262144 + 10], b[262144 + 10];
double grid_a[512][512];
int binaryResize(int n) {
int v = 2;
while(v < n)
v <<= 1;
return v;
}
double getArea(int lx, int ly, int rx, int ry) {
double ret = 0;
ret += grid_a[rx][ry];
if(lx - 1 >= 0)
ret -= grid_a[lx - 1][ry];
if(ly - 1 >= 0)
ret -= grid_a[rx][ly - 1];
if(lx - 1 >= 0 && ly - 1 >= 0)
ret += grid_a[lx - 1][ly - 1];
return ret;
}
int main() {
int testcase;
scanf("%d", &testcase);
while(testcase--) {
int m, n, p, q;
scanf("%d %d %d %d", &m, &n, &p, &q);
memset(a, 0, sizeof(a));
memset(b, 0, sizeof(b));
memset(grid_a, 0, sizeof(grid_a));
int N;
N = binaryResize(m);
N = max(N, binaryResize(n));
N = max(N, binaryResize(p));
N = max(N, binaryResize(q));
for(int i = 0; i < m; i++)
for(int j = 0; j < n; j++)
scanf("%lf", &a[i*N + j]);
double sqsum = 0;
for(int i = 0; i < p; i++) {
for(int j = 0; j < q; j++) {
scanf("%lf", &b[i*N + j]);
sqsum += b[i*N + j]*b[i*N + j];
}
}
double column_sum[512] = {};
for(int i = 0; i < N; i++) {
for(int j = 0; j < N; j++) {
column_sum[j] += a[i*N + j]*a[i*N + j];
grid_a[i][j] += column_sum[j];
if(j)
grid_a[i][j] += grid_a[i][j-1];
}
}
vector<double> r = convolution(vector<double>(a, a + N*N), vector<double>(b, b + N*N));
int qx = m - p;
int qy = n - q;
double diff = 1e+30;
int bestX = 0, bestY = 0;
for(int i = 0; i <= qx; i++) {
for(int j = 0; j <= qy; j++) {
double c = r[i*N + j];
double sq_a = getArea(i, j, i + p - 1, j + q - 1);
double sq_b = sqsum;
if(sq_a - 2*c + sq_b < diff) {
diff = sq_a - 2*c + sq_b;
bestX = i, bestY = j;
}
//printf("%lf %lf c = %lf\n", sq_a, sq_a - 2*c + sq_b, c);
}
}
printf("%d %d\n", bestX + 1, bestY + 1);
}
return 0;
}
上一篇:[高等演算法] 雜II
下一篇:[高演][攤銷] 資料結構