python image algorithm matlab image-processing

python - Algoritmo para agrupar conjuntos de puntos que siguen una dirección



image algorithm (3)

Estoy usando una versión recortada de la imagen publicada como entrada. Aquí solo abordo el caso en el que se puede pensar que la orientación de la cuadrícula es casi horizontal / vertical. Es posible que esto no aborde por completo su alcance, pero creo que puede brindarle algunos consejos.

Binarize la imagen para que los cuadrados distorsionados se llenen. Aquí uso un umbral de Otsu simple. Luego tome la transformación de distancia de esta imagen binaria.

En la imagen transformada de distancia, vemos las brechas entre los cuadrados como picos.

Para obtener líneas orientadas horizontalmente, tome los máximos locales de cada una de las columnas de la imagen de distancia y luego encuentre los componentes conectados.

Para obtener líneas orientadas verticalmente, tome los máximos locales de cada una de las filas de la imagen de distancia y luego encuentre los componentes conectados.

Las imágenes a continuación muestran las líneas horizontales y verticales así encontradas con los puntos de las esquinas como círculos.

Para componentes conectados razonablemente largos, puede ajustar una curva (una línea o un polinomio) y luego clasificar los puntos de las esquinas, por ejemplo, según la distancia a la curva, en qué lado de la curva se encuentra el punto, etc.

Hice esto en Matlab. No probé las piezas de ajuste y clasificación de curvas.

clear all; close all; im = imread(''0RqUd-1.jpg''); gr = rgb2gray(im); % any preprocessing to get a binary image that fills the distorted squares bw = ~im2bw(gr, graythresh(gr)); di = bwdist(bw); % distance transform di2 = imdilate(di, ones(3)); % propagate max corners = corner(gr); % simple corners % find regional max for each column of dist image regmxh = zeros(size(di2)); for c = 1:size(di2, 2) regmxh(:, c) = imregionalmax(di2(:, c)); end % label connected components ccomph = bwlabel(regmxh, 8); % find regional max for each row of dist image regmxv = zeros(size(di2)); for r = 1:size(di2, 1) regmxv(r, :) = imregionalmax(di2(r, :)); end % label connected components ccompv = bwlabel(regmxv, 8); figure, imshow(gr, []) hold on plot(corners(:, 1), corners(:, 2), ''ro'') figure, imshow(di, []) figure, imshow(label2rgb(ccomph), []) hold on plot(corners(:, 1), corners(:, 2), ''ro'') figure, imshow(label2rgb(ccompv), []) hold on plot(corners(:, 1), corners(:, 2), ''ro'')

Para obtener estas líneas para cuadrículas orientadas arbitrariamente, puede pensar en la imagen de distancia como un gráfico y encontrar caminos óptimos. Vea here un buen enfoque basado en gráficos.

Nota: Estoy colocando esta pregunta en las etiquetas MATLAB y Python ya que soy el más competente en estos idiomas. Sin embargo, doy la bienvenida a las soluciones en cualquier idioma.

Pregunta Preámbulo

He tomado una imagen con una lente ojo de pez. Esta imagen consiste en un patrón con un montón de objetos cuadrados. Lo que quiero hacer con esta imagen es detectar el centroide de cada uno de estos cuadrados, luego usar estos puntos para realizar una distorsión de la imagen, específicamente, estoy buscando los parámetros correctos del modelo de distorsión. Cabe señalar que no todos los cuadrados deben ser detectados. Mientras una buena mayoría de ellos lo estén, está totalmente bien ... pero ese no es el objetivo de esta publicación. El algoritmo de estimación de parámetros que ya he escrito, pero el problema es que requiere puntos que aparecen colineales en la imagen.

La pregunta base que quiero formular es la siguiente: ¿cuál es la mejor manera de agruparlos para que cada grupo consista en una línea horizontal o una línea vertical?

Antecedentes de mi problema

Esto no es realmente importante con respecto a la pregunta que estoy formulando, pero si desea saber de dónde obtuve mis datos y para comprender mejor la pregunta que estoy formulando, lea. Si no está interesado, puede pasar directamente a la sección Configuración del problema a continuación.

Un ejemplo de una imagen que estoy tratando se muestra a continuación:

Es una imagen de 960 x 960. La imagen era originalmente de mayor resolución, pero submuestreo de la imagen para facilitar un procesamiento más rápido. Como puede ver, hay un montón de patrones cuadrados que están dispersos en la imagen. Además, los centroides que he calculado están con respecto a la imagen submuestreada anterior.

La tubería que configuré para recuperar los centroides es la siguiente:

  1. Realizar una detección de bordes Canny
  2. Concéntrese en una región de interés que minimice los falsos positivos. Esta región de interés son básicamente los cuadrados sin ninguna cinta negra que cubra uno de sus lados.
  3. Encuentra todos los distintos contornos cerrados
  4. Para cada contorno cerrado distinto ...

    a. Realizar una detección de esquina de Harris

    segundo. Determine si el resultado tiene 4 puntos de esquina

    do. Si esto sucede, entonces este contorno pertenece a un cuadrado y encuentra el centroide de esta forma

    re. Si no es así, omita esta forma

  5. Coloque todos los centroides detectados del paso n. ° 4 en una matriz para un examen más detallado.

Aquí hay un resultado de ejemplo con la imagen de arriba. Cada cuadrado detectado tiene los cuatro puntos codificados por color de acuerdo con la ubicación de donde está con respecto al cuadrado en sí. Para cada centroide que he detectado, escribo una identificación justo donde ese centroide está en la imagen misma.

Con la imagen de arriba, hay 37 cuadrados detectados.

Configuración de problemas

Supongamos que tengo algunos puntos de píxeles de imagen almacenados en una matriz N x 3 . Las dos primeras columnas son las coordenadas x (horizontal) e y (vertical) donde, en el espacio de coordenadas de la imagen, la coordenada y está invertida , lo que significa que el positivo y mueve hacia abajo. La tercera columna es una ID asociada con el punto.

Aquí hay un código escrito en MATLAB que toma estos puntos, los traza en una cuadrícula 2D y etiqueta cada punto con la tercera columna de la matriz. Si lees el fondo anterior, estos son los puntos que fueron detectados por mi algoritmo descrito anteriormente.

data = [ 475. , 605.75, 1.; 571. , 586.5 , 2.; 233. , 558.5 , 3.; 669.5 , 562.75, 4.; 291.25, 546.25, 5.; 759. , 536.25, 6.; 362.5 , 531.5 , 7.; 448. , 513.5 , 8.; 834.5 , 510. , 9.; 897.25, 486. , 10.; 545.5 , 491.25, 11.; 214.5 , 481.25, 12.; 271.25, 463. , 13.; 646.5 , 466.75, 14.; 739. , 442.75, 15.; 340.5 , 441.5 , 16.; 817.75, 421.5 , 17.; 423.75, 417.75, 18.; 202.5 , 406. , 19.; 519.25, 392.25, 20.; 257.5 , 382. , 21.; 619.25, 368.5 , 22.; 148. , 359.75, 23.; 324.5 , 356. , 24.; 713. , 347.75, 25.; 195. , 335. , 26.; 793.5 , 332.5 , 27.; 403.75, 328. , 28.; 249.25, 308. , 29.; 495.5 , 300.75, 30.; 314. , 279. , 31.; 764.25, 249.5 , 32.; 389.5 , 249.5 , 33.; 475. , 221.5 , 34.; 565.75, 199. , 35.; 802.75, 173.75, 36.; 733. , 176.25, 37.]; figure; hold on; axis ij; scatter(data(:,1), data(:,2),40, ''r.''); text(data(:,1)+10, data(:,2)+10, num2str(data(:,3)));

De manera similar en Python, usando numpy y matplotlib , tenemos:

import numpy as np import matplotlib.pyplot as plt data = np.array([[ 475. , 605.75, 1. ], [ 571. , 586.5 , 2. ], [ 233. , 558.5 , 3. ], [ 669.5 , 562.75, 4. ], [ 291.25, 546.25, 5. ], [ 759. , 536.25, 6. ], [ 362.5 , 531.5 , 7. ], [ 448. , 513.5 , 8. ], [ 834.5 , 510. , 9. ], [ 897.25, 486. , 10. ], [ 545.5 , 491.25, 11. ], [ 214.5 , 481.25, 12. ], [ 271.25, 463. , 13. ], [ 646.5 , 466.75, 14. ], [ 739. , 442.75, 15. ], [ 340.5 , 441.5 , 16. ], [ 817.75, 421.5 , 17. ], [ 423.75, 417.75, 18. ], [ 202.5 , 406. , 19. ], [ 519.25, 392.25, 20. ], [ 257.5 , 382. , 21. ], [ 619.25, 368.5 , 22. ], [ 148. , 359.75, 23. ], [ 324.5 , 356. , 24. ], [ 713. , 347.75, 25. ], [ 195. , 335. , 26. ], [ 793.5 , 332.5 , 27. ], [ 403.75, 328. , 28. ], [ 249.25, 308. , 29. ], [ 495.5 , 300.75, 30. ], [ 314. , 279. , 31. ], [ 764.25, 249.5 , 32. ], [ 389.5 , 249.5 , 33. ], [ 475. , 221.5 , 34. ], [ 565.75, 199. , 35. ], [ 802.75, 173.75, 36. ], [ 733. , 176.25, 37. ]]) plt.figure() plt.gca().invert_yaxis() plt.plot(data[:,0], data[:,1], ''r.'', markersize=14) for idx in np.arange(data.shape[0]): plt.text(data[idx,0]+10, data[idx,1]+10, str(int(data[idx,2])), size=''large'') plt.show()

Obtenemos:

Volver a la pregunta

Como puede ver, estos puntos están más o menos en un patrón de cuadrícula y puede ver que podemos formar líneas entre los puntos. Específicamente, puede ver que hay líneas que se pueden formar horizontal y verticalmente.

Por ejemplo, si hace referencia a la imagen en la sección de fondo de mi problema, podemos ver que hay 5 grupos de puntos que se pueden agrupar de forma horizontal. Por ejemplo, los puntos 23, 26, 29, 31, 33, 34, 35, 37 y 36 forman un grupo. Los puntos 19, 21, 24, 28, 30 y 32 forman otro grupo, y así sucesivamente. De manera similar, en un sentido vertical, podemos ver que los puntos 26, 19, 12 y 3 forman un grupo, los puntos 29, 21, 13 y 5 forman otro grupo y así sucesivamente.

Pregunta para preguntar

Mi pregunta es esta: ¿Cuál es un método que puede agrupar con éxito puntos en agrupaciones horizontales y agrupaciones verticales por separado, dado que los puntos podrían estar en cualquier orientación?

Condiciones

  1. Debe haber al menos tres puntos por línea. Si hay algo menos que eso, entonces esto no califica como un segmento. Por lo tanto, los puntos 36 y 10 no califican como una línea vertical, y de manera similar, el punto aislado 23 no debe ser de calidad como una línea vertical, sino que es parte de la primera agrupación horizontal.

  2. El patrón de calibración anterior puede estar en cualquier orientación. Sin embargo, por lo que estoy tratando, el peor tipo de orientación que puede obtener es lo que ve arriba en la sección de antecedentes.

Rendimiento esperado

La salida sería un par de listas donde la primera lista tiene elementos donde cada elemento le da una secuencia de identificadores de punto que forman una línea horizontal. Del mismo modo, la segunda lista tiene elementos en los que cada elemento proporciona una secuencia de ID de punto que forman una línea vertical.

Por lo tanto, el resultado esperado para las secuencias horizontales sería algo como esto:

MATLAB

horiz_list = {[23, 26, 29, 31, 33, 34, 35, 37, 36], [19, 21, 24, 28, 30, 32], ...}; vert_list = {[26, 19, 12, 3], [29, 21, 13, 5], ....};

Pitón

horiz_list = [[23, 26, 29, 31, 33, 34, 35, 37, 36], [19, 21, 24, 28, 30, 32], ....] vert_list = [[26, 19, 12, 3], [29, 21, 13, 5], ...]

Lo que he intentado

Algorítmicamente, lo que he intentado es deshacer la rotación que se experimenta en estos puntos. Realicé Análisis de Componentes Principales e intenté proyectar los puntos con respecto a los vectores de base ortogonal calculados para que los puntos estuvieran más o menos en una cuadrícula rectangular recta.

Una vez que tengo eso, es solo una cuestión de hacer un poco de procesamiento de escaneo donde podría agrupar puntos en función de un cambio diferencial en las coordenadas horizontales o verticales. Debería ordenar las coordenadas por los valores y , luego examinar estas coordenadas ordenadas y buscar un gran cambio. Una vez que encuentre este cambio, puede agrupar puntos entre los cambios para formar sus líneas. Hacer esto con respecto a cada dimensión le daría las agrupaciones horizontales o verticales.

Con respecto a PCA, esto es lo que hice en MATLAB y Python:

MATLAB

%# Step #1 - Get just the data - no IDs data_raw = data(:,1:2); %# Decentralize mean data_nomean = bsxfun(@minus, data_raw, mean(data_raw,1)); %# Step #2 - Determine covariance matrix %# This already decentralizes the mean cov_data = cov(data_raw); %# Step #3 - Determine right singular vectors [~,~,V] = svd(cov_data); %# Step #4 - Transform data with respect to basis F = V.''*data_nomean.''; %# Visualize both the original data points and transformed data figure; plot(F(1,:), F(2,:), ''b.'', ''MarkerSize'', 14); axis ij; hold on; plot(data(:,1), data(:,2), ''r.'', ''MarkerSize'', 14);

Pitón

import numpy as np import numpy.linalg as la # Step #1 and Step #2 - Decentralize mean centroids_raw = data[:,:2] mean_data = np.mean(centroids_raw, axis=0) # Transpose for covariance calculation data_nomean = (centroids_raw - mean_data).T # Step #3 - Determine covariance matrix # Doesn''t matter if you do this on the decentralized result # or the normal result - cov subtracts the mean off anyway cov_data = np.cov(data_nomean) # Step #4 - Determine right singular vectors via SVD # Note - This is already V^T, so there''s no need to transpose _,_,V = la.svd(cov_data) # Step #5 - Transform data with respect to basis data_transform = np.dot(V, data_nomean).T plt.figure() plt.gca().invert_yaxis() plt.plot(data[:,0], data[:,1], ''b.'', markersize=14) plt.plot(data_transform[:,0], data_transform[:,1], ''r.'', markersize=14) plt.show()

El código anterior no solo reproyecta los datos, sino que también traza los puntos originales y los puntos proyectados en una sola figura. Sin embargo, cuando intenté reproyectar mis datos, esta es la trama que obtengo:

Los puntos en rojo son las coordenadas de la imagen original, mientras que los puntos en azul se reproyectan en los vectores de base para intentar eliminar la rotación. Todavía no acaba de hacer el trabajo. Todavía hay cierta orientación con respecto a los puntos, por lo que si traté de hacer mi algoritmo de escaneo, los puntos de las líneas a continuación para el trazado horizontal o al costado para el trazado vertical se agruparían inadvertidamente y esto no es correcto.

Tal vez estoy pensando demasiado en el problema, pero cualquier idea que tenga con respecto a esto sería muy apreciada. Si la respuesta es realmente excelente, me inclinaría a otorgar una gran recompensa, ya que he estado atascado en este problema durante bastante tiempo.

Espero que esta pregunta no sea larga. Si no tienes una idea de cómo resolver esto, te agradezco tu tiempo para leer mi pregunta independientemente.

Esperando cualquier idea que pueda tener. ¡Muchas gracias!


Nota 1: Tiene una serie de configuraciones -> que para otras imágenes pueden necesitar modificarse para obtener el resultado que desea ver% Configuración - jugar con estos valores

Nota 2: No encuentra todas las líneas que desea -> pero es un punto de partida ....

Para llamar a esta función, invoque esto en el símbolo del sistema:

>> [h, v] = testLines;

Obtenemos:

>> celldisp(h) h{1} = 1 2 4 6 9 10 h{2} = 3 5 7 8 11 14 15 17 h{3} = 1 2 4 6 9 10 h{4} = 3 5 7 8 11 14 15 17 h{5} = 1 2 4 6 9 10 h{6} = 3 5 7 8 11 14 15 17 h{7} = 3 5 7 8 11 14 15 17 h{8} = 1 2 4 6 9 10 h{9} = 1 2 4 6 9 10 h{10} = 12 13 16 18 20 22 25 27 h{11} = 13 16 18 20 22 25 27 h{12} = 3 5 7 8 11 14 15 17 h{13} = 3 5 7 8 11 14 15 h{14} = 12 13 16 18 20 22 25 27 h{15} = 3 5 7 8 11 14 15 17 h{16} = 12 13 16 18 20 22 25 27 h{17} = 19 21 24 28 30 h{18} = 21 24 28 30 h{19} = 12 13 16 18 20 22 25 27 h{20} = 19 21 24 28 30 h{21} = 12 13 16 18 20 22 24 25 h{22} = 12 13 16 18 20 22 24 25 27 h{23} = 23 26 29 31 33 34 35 h{24} = 23 26 29 31 33 34 35 37 h{25} = 23 26 29 31 33 34 35 36 37 h{26} = 33 34 35 37 36 h{27} = 31 33 34 35 37 >> celldisp(v) v{1} = 33 28 18 8 1 v{2} = 34 30 20 11 2 v{3} = 26 19 12 3 v{4} = 35 22 14 4 v{5} = 29 21 13 5 v{6} = 25 15 6 v{7} = 31 24 16 7 v{8} = 37 32 27 17 9

También se genera una figura que dibuja las líneas a través de cada conjunto adecuado de puntos:

function [horiz_list, vert_list] = testLines global counter; global colours; close all; data = [ 475. , 605.75, 1.; 571. , 586.5 , 2.; 233. , 558.5 , 3.; 669.5 , 562.75, 4.; 291.25, 546.25, 5.; 759. , 536.25, 6.; 362.5 , 531.5 , 7.; 448. , 513.5 , 8.; 834.5 , 510. , 9.; 897.25, 486. , 10.; 545.5 , 491.25, 11.; 214.5 , 481.25, 12.; 271.25, 463. , 13.; 646.5 , 466.75, 14.; 739. , 442.75, 15.; 340.5 , 441.5 , 16.; 817.75, 421.5 , 17.; 423.75, 417.75, 18.; 202.5 , 406. , 19.; 519.25, 392.25, 20.; 257.5 , 382. , 21.; 619.25, 368.5 , 22.; 148. , 359.75, 23.; 324.5 , 356. , 24.; 713. , 347.75, 25.; 195. , 335. , 26.; 793.5 , 332.5 , 27.; 403.75, 328. , 28.; 249.25, 308. , 29.; 495.5 , 300.75, 30.; 314. , 279. , 31.; 764.25, 249.5 , 32.; 389.5 , 249.5 , 33.; 475. , 221.5 , 34.; 565.75, 199. , 35.; 802.75, 173.75, 36.; 733. , 176.25, 37.]; figure; hold on; axis ij; % Change due to Benoit_11 scatter(data(:,1), data(:,2),40, ''r.''); text(data(:,1)+10, data(:,2)+10, num2str(data(:,3))); text(data(:,1)+10, data(:,2)+10, num2str(data(:,3))); % Process your data as above then run the function below(note it has sub functions) counter = 0; colours = ''bgrcmy''; [horiz_list, vert_list] = findClosestPoints ( data(:,1), data(:,2) ); function [horiz_list, vert_list] = findClosestPoints ( x, y ) % calc length of points nX = length(x); % set up place holder flags modelledH = false(nX,1); modelledV = false(nX,1); horiz_list = {}; vert_list = {}; % loop for all points for p=1:nX % have we already modelled a horizontal line through these? % second last param - true - horizontal, false - vertical if modelledH(p)==false [modelledH, index] = ModelPoints ( p, x, y, modelledH, true, true ); horiz_list = [horiz_list index]; else [~, index] = ModelPoints ( p, x, y, modelledH, true, false ); horiz_list = [horiz_list index]; end % make a temp copy of the x and y and remove any of the points modelled % from the horizontal -> this is to avoid them being found in the % second call. tempX = x; tempY = y; tempX(index) = NaN; tempY(index) = NaN; tempX(p) = x(p); tempY(p) = y(p); % Have we found a vertial line? if modelledV(p)==false [modelledV, index] = ModelPoints ( p, tempX, tempY, modelledV, false, true ); vert_list = [vert_list index]; end end end function [modelled, index] = ModelPoints ( p, x, y, modelled, method, fullRun ) % p - row in your original data matrix % x - data(:,1) % y - data(:,2) % modelled - array of flags to whether rows have been modelled % method - horizontal or vertical (used to calc graadients) % fullRun - full calc or just to get indexes % this could be made better by storing the indexes of each horizontal in the method above % Settings - play around with these values gradDelta = 0.2; % find points where gradient is less than this value gradLimit = 0.45; % if mean gradient of line is above this ignore numberOfPointsToCheck = 7; % number of points to check when look along the line % to find other points (this reduces chance of it % finding other points far away % I optimised this for your example to be 7 % Try varying it and you will see how it effect the result. % Find the index of points which are inline. [index, grad] = CalcIndex ( x, y, p, gradDelta, method ); % check gradient of line if abs(mean(grad))>gradLimit index = []; return end % add point of interest to index index = [p index]; % loop through all points found above to find any other points which are in % line with these points (this allows for slight curvature combineIndex = []; for ii=2:length(index) % Find inedex of the points found above (find points on curve) [index2] = CalcIndex ( x, y, index(ii), gradDelta, method, numberOfPointsToCheck, grad(ii-1) ); % Check that the point on this line are on the original (i.e. inline -> not at large angle if any(ismember(index,index2)) % store points found combineIndex = unique([index2 combineIndex]); end end % copy to index index = combineIndex; if fullRun % do some plotting % TODO: here you would need to calculate your arrays to output. xx = x(index); [sX,sOrder] = sort(xx); % Check its found at least 3 points if length ( index(sOrder) ) > 2 % flag the modelled on the points found modelled(index(sOrder)) = true; % plot the data plot ( x(index(sOrder)), y(index(sOrder)), colours(mod(counter,numel(colours)) + 1)); counter = counter + 1; end index = index(sOrder); end end function [index, gradCheck] = CalcIndex ( x, y, p, gradLimit, method, nPoints2Consider, refGrad ) % x - data(:,1) % y - data(:,2) % p - point of interest % method (x/y) or (y/x) % nPoints2Consider - only look at N points (options) % refgrad - rather than looking for gradient of closest point -> use this % - reference gradient to find similar points (finds points on curve) nX = length(x); % calculate gradient for g=1:nX if method grad(g) = (x(g)-x(p))/(y(g)-y(p)); else grad(g) = (y(g)-y(p))/(x(g)-x(p)); end end % find distance to all other points delta = sqrt ( (x-x(p)).^2 + (y-y(p)).^2 ); % set its self = NaN delta(delta==min(delta)) = NaN; % find the closest points [m,order] = sort(delta); if nargin == 7 % for finding along curve % set any far away points to be NaN grad(order(nPoints2Consider+1:end)) = NaN; % find the closest points to the reference gradient within the allowable limit index = find(abs(grad-refGrad)<gradLimit==1); % store output gradCheck = grad(index); else % find the points which are closes to the gradient of the closest point index = find(abs(grad-grad(order(1)))<gradLimit==1); % store gradients to output gradCheck = grad(index); end end end


Si bien no puedo sugerir un mejor enfoque para agrupar cualquier lista dada de puntos centroides que la que ya probaste, espero que la siguiente idea te pueda ayudar:

Dado que usted es muy específico sobre el contenido de su imagen (que contiene un campo de cuadrados), me preguntaba si de hecho necesita agrupar los puntos del centroide a partir de los datos proporcionados en la problem setup , o si puede usar los datos descritos en Background to the problem también. Como ya determinaste las esquinas de cada cuadrado detectado así como su posición en ese cuadrado dado, me parece que sería muy preciso determinar un vecino de un cuadrado dado comparando las coordenadas de las esquinas.

Entonces, para encontrar cualquier candidato para un vecino de la derecha de cualquier cuadrado, sugiero que compare la esquina superior derecha e inferior derecha de ese cuadrado con la esquina superior izquierda e inferior de cualquier otro cuadrado (o cualquier cuadrado dentro de una cierta distancia). Teniendo en cuenta solo pequeñas diferencias verticales y diferencias horizontales ligeramente mayores, puedes "unir" dos cuadrados, si ambos puntos de esquina correspondientes están lo suficientemente cerca.

Al usar un límite superior a la diferencia vertical / horizontal permitida entre esquinas, incluso podría asignar el mejor cuadrado coincidente dentro de estos límites como vecino.

Un problema podría ser que no se detectan todos los cuadrados, por lo que hay un espacio bastante grande entre el cuadrado 30 y 32. Como dijiste que necesitas ''al menos'' 3 cuadrados por fila, puede ser viable que simplemente ignores cuadrado 32 en esa línea horizontal. Si esa no es una opción para ti, podrías intentar hacer coincidir tantos cuadrados como sea posible y luego asignar los cuadrados "faltantes" a un punto en tu cuadrícula usando los datos calculados previamente:

En el ejemplo sobre el cuadrado 32, habría detectado que tiene vecinos superiores e inferiores 27 y 37. También debería haber podido determinar que el cuadrado 27 se encuentra dentro de la fila 1 y 37 se encuentra dentro de la fila 3, por lo que puede asignar un cuadrado 32 a la fila "mejor coincidencia" en el medio, que obviamente es 2 en este caso.

Este enfoque general es básicamente el enfoque que ya ha intentado, pero debería ser mucho más preciso ya que ahora compara la orientación y la distancia de dos líneas en lugar de simplemente comparar la ubicación de dos puntos en una cuadrícula.

También como un sidenode en sus intentos anteriores, ¿puede usar las líneas negras del borde para corregir un poco la rotación inicial de la imagen? Esto podría hacer que otros algoritmos de distorsión (como los que discutiste con knedlsepp en los comentarios) sean mucho más precisos. (EDIT: ya leí los comentarios de Parag, comparar los puntos por el ángulo de las líneas es, por supuesto, básicamente lo mismo que rotar la imagen de antemano)