Matchings

Explicación

Queremos colocar la imagen más pequeña (imagen del robot) sobre la imagen más grande (imagen del suelo) de modo que se maximice el número de píxeles coincidentes.

Sea la imagen del robot (máscara) una grilla $n \times m$ , y la imagen del suelo una grilla $n_1 \times m_1$ , donde $n \le n_1$ y $m \le m_1$ .

Para una colocación fija $(x, y)$ , contamos:

\sum_{i=0}^{n-1} \sum_{j=0}^{m-1} [\text{mask}_{i,j} = \text{floor}_{i+y, j+x}]

Como los valores son solo $0$ o $1$ , separamos esto en dos contribuciones independientes:

posiciones donde ambos son $1$
posiciones donde ambos son $0$

Definimos arreglos indicadores:

$A_{i,j} = \text{mask}_{i,j}$
$B_{i,j} = \text{floor}_{i,j}$

Entonces las coincidencias de $1$ son:

\sum A_{i,j} \cdot B_{i+y, j+x}

Para contar las coincidencias de $0$ , observamos:

(1 - A_{i,j}) \cdot (1 - B_{i+y, j+x})

Así las coincidencias totales quedan:

\sum A_{i,j} \cdot B_{i+y, j+x} + \sum (1 - A_{i,j}) \cdot (1 - B_{i+y, j+x})

Computamos ambos términos por separado y los sumamos.

¿Se puede resolver con una sola FFT?

En lugar de computar dos convoluciones separadas (una para coincidencias de $1$ s y otra para coincidencias de $0$ s), podemos reducir el problema a una sola convolución.

Transformamos los arreglos así:

reemplazar $0 \rightarrow -1$
dejar $1 \rightarrow 1$

Ahora consideremos la contribución de una sola posición:

si los píxeles coinciden, obtenemos $1 \cdot 1 = 1$ o $(-1) \cdot (-1) = 1$
si no coinciden, obtenemos $1 \cdot (-1) = -1$

Así la convolución computa directamente:

\text{sum} = \text{matches} - \text{mismatches}

También sabemos:

\text{matches} + \text{mismatches} = n \cdot m

Combinando ambos,

\text{matches} = \frac{\text{sum} + n \cdot m}{2}

Así una sola convolución basta para recuperar el número de píxeles coincidentes.

Convertir a multiplicación de polinomios

Reducimos el problema 2D a 1D usando aplanado por filas (row-major).

Mapeamos la posición $(i, j)$ al índice:

\text{index} = i \cdot m_1 + j

Construimos arreglos:

Máscara aplanada → a
Suelo aplanado → b

Para simular correctamente el deslizamiento, revertimos el arreglo a.

Ahora la contribución:

\sum A_{i,j} \cdot B_{i+y, j+x}

se convierte en una convolución estándar.

Tratar los bordes 2D en una convolución 1D

Cuando aplanamos los arreglos 2D a 1D, perdemos la noción de fronteras de fila. Una convolución 1D estándar computa valores para todos los desplazamientos posibles, incluidos los inválidos en los que la máscara cruza fronteras de fila o se sale del suelo.

Estos desplazamientos inválidos pueden producir valores sin sentido, pero eso no afecta la corrección. Esto es porque la convolución computa cada desplazamiento de forma independiente, así que los alineamientos inválidos no interfieren con los válidos.

Al restringirnos a los índices que corresponden a colocaciones válidas de la esquina superior izquierda $(0 \le x \le m_1 - m,\ 0 \le y \le n_1 - n)$ , solo leemos entradas que representan superposiciones 2D correctas. Todos los casos inválidos de wrap-around quedan en índices distintos y se ignoran con seguridad.

Después de la multiplicación, cada índice del arreglo resultante corresponde a un desplazamiento específico. Para una colocación $(x, y)$ , el índice correspondiente es:

i = (y \cdot m_1 + x) + (|a| - 1)

Este valor da el número de píxeles $1$ coincidentes para esta colocación.

Repetimos el mismo proceso para $(1 - A)$ y $(1 - B)$ para contar los píxeles $0$ coincidentes, y sumamos los resultados:

\text{matches} = c_1[i] + c_2[i]

Iteramos sobre todas las colocaciones válidas $(x, y)$ y nos quedamos con las que alcanzan el valor máximo.

Optimizar con FFT

Hacerlo por fuerza bruta corre en tiempo $\mathcal{O}(n \cdot m \cdot n_1 \cdot m_1)$ .

Con la Transformada Rápida de Fourier (FFT), la convolución corre en tiempo $\mathcal{O}(N \log N)$ , donde $N = n_1 \cdot m_1$ .

¿Se puede reducir el tamaño de la FFT?

En la multiplicación de polinomios estándar, elegimos el tamaño de la FFT de al menos $|a| + |b|$ para evitar el solapamiento de la convolución circular.

Sin embargo, en este problema solo nos importan los desplazamientos en los que la máscara queda completamente dentro de la imagen del suelo.

El desplazamiento válido máximo se queda dentro del rango de índices que consultamos, que está acotado por $|b|$ , así que no necesitamos el rango completo de la convolución.

Así, basta elegir el tamaño de la FFT como la menor potencia de $2$ tal que:

n \ge |b|

en lugar de $|a| + |b|$ .

Esto reduce tanto el uso de memoria como el tiempo de ejecución aproximadamente a la mitad.

Implementación

Complejidad temporal: $\mathcal{O}(N \log N)$ donde $N = n_1 \cdot m_1$

C++


 
#include <bits/stdc++.h>
using namespace std;
 
using cd = complex<double>;
const double PI = acos(-1);
 
// BeginCodeSnip{FFT Template}
void fft(vector<cd> &a, bool inv) {
	int n = a.size();
	if (n == 1) return;
 
	vector<cd> a0(n / 2), a1(n / 2);
	for (int i = 0; 2 * i < n; i++) a0[i] = a[2 * i], a1[i] = a[2 * i + 1];
 
	fft(a0, inv);
	fft(a1, inv);
 
	double ang = 2 * PI / n * (inv ? -1 : 1);
	cd w(1), wn(cos(ang), sin(ang));
 
	for (int i = 0; 2 * i < n; i++) {
		a[i] = a0[i] + w * a1[i];
		a[i + n / 2] = a0[i] - w * a1[i];
		if (inv) a[i] /= 2, a[i + n / 2] /= 2;
		w *= wn;
	}
}
 
vector<int> mul(vector<int> a, vector<int> b) {
	int n = 1;
	while (n < (int)a.size() + (int)b.size()) n <<= 1;
 
	vector<cd> fa(a.begin(), a.end()), fb(b.begin(), b.end());
	fa.resize(n);
	fb.resize(n);
 
	fft(fa, 0);
	fft(fb, 0);
 
	for (int i = 0; i < n; i++) fa[i] *= fb[i];
 
	fft(fa, 1);
 
	vector<int> r(n);
	for (int i = 0; i < n; i++) r[i] = round(fa[i].real());
 
	return r;
}
// EndCodeSnip
 
int main() {
	ios::sync_with_stdio(0);
	cin.tie(0);
 
	int m, n;
	cin >> m >> n;
 
	vector<vector<int>> mask(n, vector<int>(m));
	for (int i = 0; i < n; ++i)
		for (int j = 0; j < m; ++j) cin >> mask[i][j];
 
	int m1, n1;
	cin >> m1 >> n1;
 
	vector<int> a((n - 1) * m1 + m, 0), b(n1 * m1);
 
	for (int i = 0; i < n1; i++)
		for (int j = 0; j < m1; j++) cin >> b[i * m1 + j];
 
	for (int i = 0; i < n; i++)
		for (int j = 0; j < m; j++) a[i * m1 + j] = mask[i][j];
 
	reverse(a.begin(), a.end());
	auto c1 = mul(a, b);
 
	for (int i = 0; i < n; i++)
		for (int j = 0; j < m; j++) a[i * m1 + j] = 1 - mask[i][j];
 
	reverse(a.begin(), a.end());
	for (int &x : b) x = 1 - x;
 
	auto c2 = mul(a, b);
 
	vector<int> c(c1.size());
	for (int i = 0; i < (int)c.size(); i++) c[i] = c1[i] + c2[i];
 
	int best = 0;
 
	for (int y = 0; y <= n1 - n; y++) {
		for (int x = 0; x <= m1 - m; x++) {
			int ind = y * m1 + x;
			int i = ind + (int)a.size() - 1;
			best = max(best, c[i]);
		}
	}
 
	vector<pair<int, int>> ans;
 
	for (int y = 0; y <= n1 - n; y++) {
		for (int x = 0; x <= m1 - m; x++) {
			int ind = y * m1 + x;
			int i = ind + (int)a.size() - 1;
			if (c[i] == best) ans.push_back({x, y});
		}
	}
 
	sort(ans.begin(), ans.end());
	for (auto &p : ans) cout << p.first << " " << p.second << '\n';
}