algorithm - Cómo dibujar una cuadrícula Perspective-Correct en 2D

graphics geometry (11)

Tengo una aplicación que define un rectángulo del mundo real en la parte superior de una imagen / fotografía, por supuesto en 2D, puede que no sea un rectángulo porque lo estás mirando desde un ángulo.

El problema es decir que el rectángulo debe tener líneas de cuadrícula dibujadas, por ejemplo si es 3x5, así que necesito dibujar 2 líneas del lado 1 al lado 3 y 4 líneas del lado 2 al lado 4.

A partir de ahora estoy dividiendo cada línea en partes equidistantes, para obtener el punto inicial y final de todas las líneas de la cuadrícula. Sin embargo, cuanto más anguloso es el rectángulo, más "incorrectas" se vuelven estas líneas, ya que las líneas horizontales más alejadas de usted deben estar más juntas.

¿Alguien sabe el nombre del algoritmo que debería estar buscando?

Sí, sé que puedes hacer esto en 3D, sin embargo, estoy limitado a 2D para esta aplicación en particular.

El problema es que es la transformación de 3D a 2D lo que te atrapa.

Aquí hay un tutorial sobre cómo se hace.

Lo que necesita hacer es representarlo en 3D (mundo) y luego proyectarlo a 2D (pantalla).

Esto requerirá que uses una matriz de transformación 4D que hace la proyección en un 4D homogéneo hasta un vector homogéneo 3D, que luego puedes convertir a un vector de espacio de pantalla 2D.

Tampoco pude encontrarlo en Google, pero un buen libro de gráficos por computadora tendrá los detalles.

Las palabras clave son matriz de proyección, transformación de proyección, transformación afín, vector homogéneo, espacio mundial, espacio de pantalla, transformación de perspectiva, transformación 3D

Y, por cierto, esto generalmente requiere algunas conferencias para explicar todo eso. Buena suerte.

Si bien mi Google-Fu no ha podido producir ninguna solución matemática sólida como la roca, tal vez este dibujo que encontré podría ayudar un poco.

http://studiochalkboard.evansville.edu/lp-diminish.html

Creo que en realidad podría ser bastante difícil encontrar la matemática correcta por tu cuenta, es probable que sea una especie de expresión logarítmica o de suma. Es de esperar que el dibujo y los términos en ese enlace puedan proporcionarle algo un poco más buscable.

Al usar el método de subdivisión de Breton (que está relacionado con el método de extensión de Mongo), obtendrás divisiones precisas arbitrarias de poder de dos. Para dividirse en divisiones que no sean de potencia de dos usando esos métodos, deberá subdividir el espacio entre subpíxeles, que puede ser computacionalmente costoso.

Sin embargo, creo que puedes aplicar una variación del Teorema de Haga (que se usa en origami para dividir un lado en Nths dado un lado dividido en (N-1) ths) a las subdivisiones de perspectiva-cuadrado para producir divisiones arbitrarias desde la potencia más cercana a 2 sin tener que continuar subdividiendo.

Aquí está la solución: http://freespace.virgin.net/hugo.elias/graphics/x_persp.htm

La idea básica es que puede encontrar el "centro" correcto de la perspectiva de su rectángulo conectando las esquinas diagonalmente. La intersección de las dos líneas resultantes es su centro correcto de perspectiva. Desde allí, subdivide su rectángulo en cuatro rectángulos más pequeños, y repite el proceso. El número de veces depende de qué tan preciso lo desee. Puede subdividir justo por debajo del tamaño de un píxel para obtener una perspectiva efectivamente perfecta.

Luego, en sus subrectangles, simplemente aplica sus triángulos, rectángulos o lo que sea, "rectángulos" sin corregir estándar.

Puede realizar este algoritmo sin tener que enfrentarse al complejo problema de construir un mundo 3D "real". también es bueno si tienes un mundo real en 3D modelado, pero tus textriangles no tienen perspectiva corregida en el hardware, o necesitas una forma eficiente de obtener planos correctos de perspectiva sin trucos de representación por píxel.

En el caso especial cuando miras perpendicular a los lados 1 y 3, puedes dividir esos lados en partes iguales. Luego dibuja una diagonal y dibuja paralelos al lado 1 a través de cada intersección de la diagonal y las líneas divisorias dibujadas anteriormente.

La solución más elegante y rápida sería encontrar la matriz de homografía, que asigna las coordenadas del rectángulo a las coordenadas de la foto.

Con una biblioteca matricial decente no debería ser una tarea difícil, siempre y cuando conozca sus cálculos matemáticos.

Palabras clave: Collineación, Homografía, Transformación Lineal Directa

Sin embargo, el algoritmo recursivo anterior debería funcionar, pero probablemente si sus recursos son limitados, la geometría proyectiva es la única manera de hacerlo.

Imagen: Ejemplo de transformación bilineal y perspectiva (Nota: la altura de las líneas horizontales superior e inferior de la cuadrícula es en realidad la mitad de la altura de las líneas de descanso, en ambos dibujos)

========================================

Sé que esta es una vieja pregunta, pero tengo una solución genérica, así que decidí publicarla y creo que será útil para los futuros lectores. El código siguiente puede dibujar una cuadrícula de perspectiva arbitraria sin la necesidad de cálculos repetitivos.

Comienzo en realidad con un problema similar: dibujar una cuadrícula de perspectiva 2D y luego transformar la imagen de subrayado para restaurar la perspectiva.

Empecé a leer aquí: http://www.imagemagick.org/Usage/distorts/#bilinear_forward

y luego aquí (la Biblioteca de Leptonica): http://www.leptonica.com/affine.html

donde encontré esto

Cuando miras un objeto en un plano desde una dirección arbitraria a una distancia finita, obtienes una distorsión "trapezoidal" adicional en la imagen. Esta es una transformada proyectiva, que mantiene rectas las líneas rectas, pero no conserva los ángulos entre las líneas. Esta deformación no puede describirse mediante una transformación afín lineal, y de hecho difiere según los términos dependientes de xey en el denominador.

La transformación no es lineal, como mucha gente ya señaló en este hilo. Implica resolver un sistema lineal de 8 ecuaciones (una vez) para calcular los 8 coeficientes requeridos y luego puede usarlos para transformar tantos puntos como desee.

Para evitar incluir toda la biblioteca de Leptonica en mi proyecto, extraje algunos fragmentos de código, eliminé todos los tipos de datos y macros especiales de Leptonica, arreglé algunas pérdidas de memoria y las convertí a una clase C ++ (principalmente por razones de encapsulación) que solo hace una cosa: asigna una coordenada flotante Q (Qt) QPointF (x, y) a la Coordenada Perspectiva correspondiente.

Si desea adaptar el código a otra biblioteca C ++, lo único que debe redefinir / sustituir es la clase de coordenadas QPointF.

Espero que algunos lectores futuros lo encuentren útil. El código a continuación se divide en 3 partes:

A. Un ejemplo sobre cómo usar la clase genImageProjective C ++ para dibujar una cuadrícula 2D en perspectiva

B. archivo genImageProjective.h

Archivo C. genImageProjective.cpp

//============================================================ // C++ Code Example on how to use the // genImageProjective class to draw a perspective 2D Grid //============================================================ #include "genImageProjective.h" // Input: 4 Perspective-Tranformed points: // perspPoints[0] = top-left // perspPoints[1] = top-right // perspPoints[2] = bottom-right // perspPoints[3] = bottom-left void drawGrid(QPointF *perspPoints) { (...) // Setup a non-transformed area rectangle // I use a simple square rectangle here because in this case we are not interested in the source-rectangle, // (we want to just draw a grid on the perspPoints[] area) // but you can use any arbitrary rectangle to perform a real mapping to the perspPoints[] area QPointF topLeft = QPointF(0,0); QPointF topRight = QPointF(1000,0); QPointF bottomRight = QPointF(1000,1000); QPointF bottomLeft = QPointF(0,1000); float width = topRight.x() - topLeft.x(); float height = bottomLeft.y() - topLeft.y(); // Setup Projective trasform object genImageProjective imageProjective; imageProjective.sourceArea[0] = topLeft; imageProjective.sourceArea[1] = topRight; imageProjective.sourceArea[2] = bottomRight; imageProjective.sourceArea[3] = bottomLeft; imageProjective.destArea[0] = perspPoints[0]; imageProjective.destArea[1] = perspPoints[1]; imageProjective.destArea[2] = perspPoints[2]; imageProjective.destArea[3] = perspPoints[3]; // Compute projective transform coefficients if (imageProjective.computeCoeefficients() != 0) return; // This can actually fail if any 3 points of Source or Dest are colinear // Initialize Grid parameters (without transform) float gridFirstLine = 0.1f; // The normalized position of first Grid Line (0.0 to 1.0) float gridStep = 0.1f; // The normalized Grd size (=distance between grid lines: 0.0 to 1.0) // Draw Horizonal Grid lines QPointF lineStart, lineEnd, tempPnt; for (float pos = gridFirstLine; pos <= 1.0f; pos += gridStep) { // Compute Grid Line Start tempPnt = QPointF(topLeft.x(), topLeft.y() + pos*width); imageProjective.mapSourceToDestPoint(tempPnt, lineStart); // Compute Grid Line End tempPnt = QPointF(topRight.x(), topLeft.y() + pos*width); imageProjective.mapSourceToDestPoint(tempPnt, lineEnd); // Draw Horizontal Line (use your prefered method to draw the line) (...) } // Draw Vertical Grid lines for (float pos = gridFirstLine; pos <= 1.0f; pos += gridStep) { // Compute Grid Line Start tempPnt = QPointF(topLeft.x() + pos*height, topLeft.y()); imageProjective.mapSourceToDestPoint(tempPnt, lineStart); // Compute Grid Line End tempPnt = QPointF(topLeft.x() + pos*height, bottomLeft.y()); imageProjective.mapSourceToDestPoint(tempPnt, lineEnd); // Draw Vertical Line (use your prefered method to draw the line) (...) } (...) } ========================================== //======================================== //C++ Header File: genImageProjective.h //======================================== #ifndef GENIMAGE_H #define GENIMAGE_H #include <QPointF> // Class to transform an Image Point using Perspective transformation class genImageProjective { public: genImageProjective(); int computeCoeefficients(void); int mapSourceToDestPoint(QPointF& sourcePoint, QPointF& destPoint); public: QPointF sourceArea[4]; // Source Image area limits (Rectangular) QPointF destArea[4]; // Destination Image area limits (Perspectivelly Transformed) private: static int gaussjordan(float **a, float *b, int n); bool coefficientsComputed; float vc[8]; // Vector of Transform Coefficients }; #endif // GENIMAGE_H //======================================== //======================================== //C++ CPP File: genImageProjective.cpp //======================================== #include <math.h> #include "genImageProjective.h" // ---------------------------------------------------- // class genImageProjective // ---------------------------------------------------- genImageProjective::genImageProjective() { sourceArea[0] = sourceArea[1] = sourceArea[2] = sourceArea[3] = QPointF(0,0); destArea[0] = destArea[1] = destArea[2] = destArea[3] = QPointF(0,0); coefficientsComputed = false; } // -------------------------------------------------------------- // Compute projective transform coeeeficients // RetValue: 0: Success, !=0: Error /*-------------------------------------------------------------* * Projective coordinate transformation * *-------------------------------------------------------------*/ /*! * computeCoeefficients() * * Input: this->sourceArea[4]: (source 4 points; unprimed) * this->destArea[4]: (transformed 4 points; primed) * this->vc (computed vector of transform coefficients) * Return: 0 if OK; <0 on error * * We have a set of 8 equations, describing the projective * transformation that takes 4 points (sourceArea) into 4 other * points (destArea). These equations are: * * x1'' = (c[0]*x1 + c[1]*y1 + c[2]) / (c[6]*x1 + c[7]*y1 + 1) * y1'' = (c[3]*x1 + c[4]*y1 + c[5]) / (c[6]*x1 + c[7]*y1 + 1) * x2'' = (c[0]*x2 + c[1]*y2 + c[2]) / (c[6]*x2 + c[7]*y2 + 1) * y2'' = (c[3]*x2 + c[4]*y2 + c[5]) / (c[6]*x2 + c[7]*y2 + 1) * x3'' = (c[0]*x3 + c[1]*y3 + c[2]) / (c[6]*x3 + c[7]*y3 + 1) * y3'' = (c[3]*x3 + c[4]*y3 + c[5]) / (c[6]*x3 + c[7]*y3 + 1) * x4'' = (c[0]*x4 + c[1]*y4 + c[2]) / (c[6]*x4 + c[7]*y4 + 1) * y4'' = (c[3]*x4 + c[4]*y4 + c[5]) / (c[6]*x4 + c[7]*y4 + 1) * * Multiplying both sides of each eqn by the denominator, we get * * AC = B * * where B and C are column vectors * * B = [ x1'' y1'' x2'' y2'' x3'' y3'' x4'' y4'' ] * C = [ c[0] c[1] c[2] c[3] c[4] c[5] c[6] c[7] ] * * and A is the 8x8 matrix * * x1 y1 1 0 0 0 -x1*x1'' -y1*x1'' * 0 0 0 x1 y1 1 -x1*y1'' -y1*y1'' * x2 y2 1 0 0 0 -x2*x2'' -y2*x2'' * 0 0 0 x2 y2 1 -x2*y2'' -y2*y2'' * x3 y3 1 0 0 0 -x3*x3'' -y3*x3'' * 0 0 0 x3 y3 1 -x3*y3'' -y3*y3'' * x4 y4 1 0 0 0 -x4*x4'' -y4*x4'' * 0 0 0 x4 y4 1 -x4*y4'' -y4*y4'' * * These eight equations are solved here for the coefficients C. * * These eight coefficients can then be used to find the mapping * (x,y) --> (x'',y''): * * x'' = (c[0]x + c[1]y + c[2]) / (c[6]x + c[7]y + 1) * y'' = (c[3]x + c[4]y + c[5]) / (c[6]x + c[7]y + 1) * */ int genImageProjective::computeCoeefficients(void) { int retValue = 0; int i; float *a[8]; /* 8x8 matrix A */ float *b = this->vc; /* rhs vector of primed coords X''; coeffs returned in vc[] */ b[0] = destArea[0].x(); b[1] = destArea[0].y(); b[2] = destArea[1].x(); b[3] = destArea[1].y(); b[4] = destArea[2].x(); b[5] = destArea[2].y(); b[6] = destArea[3].x(); b[7] = destArea[3].y(); for (i = 0; i < 8; i++) a[i] = NULL; for (i = 0; i < 8; i++) { if ((a[i] = (float *)calloc(8, sizeof(float))) == NULL) { retValue = -100; // ERROR_INT("a[i] not made", procName, 1); goto Terminate; } } a[0][0] = sourceArea[0].x(); a[0][1] = sourceArea[0].y(); a[0][2] = 1.; a[0][6] = -sourceArea[0].x() * b[0]; a[0][7] = -sourceArea[0].y() * b[0]; a[1][3] = sourceArea[0].x(); a[1][4] = sourceArea[0].y(); a[1][5] = 1; a[1][6] = -sourceArea[0].x() * b[1]; a[1][7] = -sourceArea[0].y() * b[1]; a[2][0] = sourceArea[1].x(); a[2][1] = sourceArea[1].y(); a[2][2] = 1.; a[2][6] = -sourceArea[1].x() * b[2]; a[2][7] = -sourceArea[1].y() * b[2]; a[3][3] = sourceArea[1].x(); a[3][4] = sourceArea[1].y(); a[3][5] = 1; a[3][6] = -sourceArea[1].x() * b[3]; a[3][7] = -sourceArea[1].y() * b[3]; a[4][0] = sourceArea[2].x(); a[4][1] = sourceArea[2].y(); a[4][2] = 1.; a[4][6] = -sourceArea[2].x() * b[4]; a[4][7] = -sourceArea[2].y() * b[4]; a[5][3] = sourceArea[2].x(); a[5][4] = sourceArea[2].y(); a[5][5] = 1; a[5][6] = -sourceArea[2].x() * b[5]; a[5][7] = -sourceArea[2].y() * b[5]; a[6][0] = sourceArea[3].x(); a[6][1] = sourceArea[3].y(); a[6][2] = 1.; a[6][6] = -sourceArea[3].x() * b[6]; a[6][7] = -sourceArea[3].y() * b[6]; a[7][3] = sourceArea[3].x(); a[7][4] = sourceArea[3].y(); a[7][5] = 1; a[7][6] = -sourceArea[3].x() * b[7]; a[7][7] = -sourceArea[3].y() * b[7]; retValue = gaussjordan(a, b, 8); Terminate: // Clean up for (i = 0; i < 8; i++) { if (a[i]) free(a[i]); } this->coefficientsComputed = (retValue == 0); return retValue; } /*-------------------------------------------------------------* * Gauss-jordan linear equation solver * *-------------------------------------------------------------*/ /* * gaussjordan() * * Input: a (n x n matrix) * b (rhs column vector) * n (dimension) * Return: 0 if ok, 1 on error * * Note side effects: * (1) the matrix a is transformed to its inverse * (2) the vector b is transformed to the solution X to the * linear equation AX = B * * Adapted from "Numerical Recipes in C, Second Edition", 1992 * pp. 36-41 (gauss-jordan elimination) */ #define SWAP(a,b) {temp = (a); (a) = (b); (b) = temp;} int genImageProjective::gaussjordan(float **a, float *b, int n) { int retValue = 0; int i, icol=0, irow=0, j, k, l, ll; int *indexc = NULL, *indexr = NULL, *ipiv = NULL; float big, dum, pivinv, temp; if (!a) { retValue = -1; // ERROR_INT("a not defined", procName, 1); goto Terminate; } if (!b) { retValue = -2; // ERROR_INT("b not defined", procName, 1); goto Terminate; } if ((indexc = (int *)calloc(n, sizeof(int))) == NULL) { retValue = -3; // ERROR_INT("indexc not made", procName, 1); goto Terminate; } if ((indexr = (int *)calloc(n, sizeof(int))) == NULL) { retValue = -4; // ERROR_INT("indexr not made", procName, 1); goto Terminate; } if ((ipiv = (int *)calloc(n, sizeof(int))) == NULL) { retValue = -5; // ERROR_INT("ipiv not made", procName, 1); goto Terminate; } for (i = 0; i < n; i++) { big = 0.0; for (j = 0; j < n; j++) { if (ipiv[j] != 1) { for (k = 0; k < n; k++) { if (ipiv[k] == 0) { if (fabs(a[j][k]) >= big) { big = fabs(a[j][k]); irow = j; icol = k; } } else if (ipiv[k] > 1) { retValue = -6; // ERROR_INT("singular matrix", procName, 1); goto Terminate; } } } } ++(ipiv[icol]); if (irow != icol) { for (l = 0; l < n; l++) SWAP(a[irow][l], a[icol][l]); SWAP(b[irow], b[icol]); } indexr[i] = irow; indexc[i] = icol; if (a[icol][icol] == 0.0) { retValue = -7; // ERROR_INT("singular matrix", procName, 1); goto Terminate; } pivinv = 1.0 / a[icol][icol]; a[icol][icol] = 1.0; for (l = 0; l < n; l++) a[icol][l] *= pivinv; b[icol] *= pivinv; for (ll = 0; ll < n; ll++) { if (ll != icol) { dum = a[ll][icol]; a[ll][icol] = 0.0; for (l = 0; l < n; l++) a[ll][l] -= a[icol][l] * dum; b[ll] -= b[icol] * dum; } } } for (l = n - 1; l >= 0; l--) { if (indexr[l] != indexc[l]) { for (k = 0; k < n; k++) SWAP(a[k][indexr[l]], a[k][indexc[l]]); } } Terminate: if (indexr) free(indexr); if (indexc) free(indexc); if (ipiv) free(ipiv); return retValue; } // -------------------------------------------------------------- // Map a source point to destination using projective transform // -------------------------------------------------------------- // Params: // sourcePoint: initial point // destPoint: transformed point // RetValue: 0: Success, !=0: Error // -------------------------------------------------------------- // Notes: // 1. You must call once computeCoeefficients() to compute // the this->vc[] vector of 8 coefficients, before you call // mapSourceToDestPoint(). // 2. If there was an error or the 8 coefficients were not computed, // a -1 is returned and destPoint is just set to sourcePoint value. // -------------------------------------------------------------- int genImageProjective::mapSourceToDestPoint(QPointF& sourcePoint, QPointF& destPoint) { if (coefficientsComputed) { float factor = 1.0f / (vc[6] * sourcePoint.x() + vc[7] * sourcePoint.y() + 1.); destPoint.setX( factor * (vc[0] * sourcePoint.x() + vc[1] * sourcePoint.y() + vc[2]) ); destPoint.setY( factor * (vc[3] * sourcePoint.x() + vc[4] * sourcePoint.y() + vc[5]) ); return 0; } else // There was an error while computing coefficients { destPoint = sourcePoint; // just copy the source to destination... return -1; // ...and return an error } } //========================================

Esta es una solución geométrica que pensé. No sé si el ''algoritmo'' tiene un nombre.

Digamos que quiere comenzar dividiendo el ''rectángulo'' en n piezas con líneas verticales primero.

El objetivo es colocar los puntos P1..Pn-1 en la línea superior, que podemos utilizar para trazar líneas a través de ellos hasta los puntos donde la línea izquierda y derecha se encuentran o son paralelos a ellos cuando dicho punto no existe.

Si las líneas superior e inferior son paralelas entre sí, solo coloque esos puntos para dividir la línea superior entre las esquinas de forma equidistante.

De lo contrario, coloque n puntos Q1..Qn en la línea izquierda de modo que estas y la esquina superior izquierda sean equidistantes e i <j => Qi esté más cerca de la esquina superior izquierda que Qj. Para mapear los puntos Q a la línea superior, encuentre la intersección S de la línea desde Qn a través de la esquina superior derecha y la línea paralela a la izquierda a través de la intersección de las líneas superior e inferior. Ahora conecta S con Q1..Qn-1. La intersección de las nuevas líneas con la línea superior son los P-puntos deseados.

Haz este análogo para las líneas horizontales.

Dada una rotación alrededor del eje y, especialmente si las superficies de rotación son planas, la perspectiva se genera mediante gradientes verticales. Estos se acercan progresivamente en perspectiva. En lugar de usar diagonales para definir cuatro rectángulos, que pueden funcionar dados potencias de dos ... defina dos rectángulos, a la izquierda y a la derecha. Serán más altos que anchos, eventualmente, si uno continúa dividiendo la superficie en segmentos verticales más estrechos. Esto puede acomodar superficies que no son cuadradas. Si una rotación está alrededor del eje x, entonces se necesitan gradientes horizontales.

Creo que la respuesta seleccionada no es la mejor solución disponible. Una mejor solución es aplicar la transformación en perspectiva (proyectiva) de un rectángulo a una cuadrícula simple siguiendo el guión de Matlab y la presentación de imágenes. Puede implementar este algoritmo con C ++ y OpenCV también.

function drawpersgrid sz = [ 24, 16 ]; % [x y] srcpt = [ 0 0; sz(1) 0; 0 sz(2); sz(1) sz(2)]; destpt = [ 20 50; 100 60; 0 150; 200 200;]; % make rectangular grid [X,Y] = meshgrid(0:sz(1),0:sz(2)); % find projective transform matching corner points tform = maketform(''projective'',srcpt,destpt); % apply the projective transform to the grid [X1,Y1] = tformfwd(tform,X,Y); hold on; %% find grid for i=1:sz(2) for j=1:sz(1) x = [ X1(i,j);X1(i,j+1);X1(i+1,j+1);X1(i+1,j);X1(i,j)]; y = [ Y1(i,j);Y1(i,j+1);Y1(i+1,j+1);Y1(i+1,j);Y1(i,j)]; plot(x,y,''b''); end end hold off;