Geometría

Algoritmos geométricos

El algoritmo de Cyrus-Beck es un algoritmo de recorte de líneas y polígonos convexos. De forma similar al algoritmo de Cohen-Sutherland también utiliza aristas extendidas. Se puede adaptar muy fácilmente entre rayos y polígonos convexos o segmentos y polígonos convexos.

Implementación del algoritmo de Cyrus-Beck en C#

internal sealed class CyrusBeckClipping : IClippingAlgorithm {
    private List<Vector2> _clipArea = new List<Vector2>();
    private List<Vector2> _normals = new List<Vector2>();
 
    public IEnumerable<Vector2> GetBoundingPolygon() {
        return _clipArea;
    }
 
    public void SetBoundingRectangle(Vector2 start, Vector2 end) {
        _clipArea.Clear();
 
        _clipArea.Add(start);
        _clipArea.Add(new Vector2(end.X, start.Y));
        _clipArea.Add(end);
        _clipArea.Add(new Vector2(start.X, end.Y));
 
        computeNormals();
    }
 
    public void SetBoundingPolygon(IEnumerable<Vector2> points) {
        _clipArea.Clear();
        _clipArea.AddRange(points);
        computeNormals();
    }
 
    private void computeNormals() {
        _normals.Clear();
 
        for (int i = 0; i < _clipArea.Count - 1; i++) {
            Vector2 direction = _clipArea[i + 1] - _clipArea[i];
            direction.Normalize();
 
            _normals.Add(new Vector2(-direction.Y, direction.X));
        }
 
        {
            Vector2 direction = _clipArea[0] - _clipArea[_clipArea.Count - 1];
            direction.Normalize();
 
            _normals.Add(new Vector2(-direction.Y, direction.X));
        }
    }
 
    public bool ClipLine(ref Line line) {
        Vector2 P = line.End - line.Start;
        float tMinimum = 0, tMaximum = 1;
        const float epsilon = 0.0001f;
 
        for (int i = 0; i < _clipArea.Count; i++) {
            Vector2 F = _clipArea[i];
            Vector2 N = _normals[i];
            Vector2 Q = line.Start - F;
 
            float Pn = Vector2.DotProduct(P, N);
            float Qn = Vector2.DotProduct(Q, N);
 
            if (Pn < epsilon && Pn > -epsilon) {
                if (Qn < 0) return false;
            }
            else {
                float computedT = -Qn / Pn;
                if (Pn < 0) {
                    if (computedT < tMinimum)
                        return false;
                    if (computedT < tMaximum)
                        tMaximum = computedT;
                }
                else {
                    if (computedT > tMaximum)
                        return false;
                    if (computedT > tMinimum)
                        tMinimum = computedT;
                }
            }
        }
 
        if (tMinimum < tMaximum) {
            if (tMaximum < 1)
                line.End = line.Start + tMaximum * P;
 
            if (tMinimum > 0)
                line.Start = line.Start + tMinimum * P;
        }
        else return false;
 
        return true;
    }
 
    public ClippingCapabilities Capabilities {
        get {
            return ClippingCapabilities.ConvexWindow |
                   ClippingCapabilities.RectangleWindow;
        }
    }
 
    public override string ToString() {
        return "Cyrus-Beck algorithm";
    }
}
// This code was implemented by Grishul Eugeny as part of preparation
// to exam in ITMO university

Ilustración del algoritmo de Cyrus-Beck.







Algoritmo de Cyrus-Beck

^

Este algoritmo se basa en calcular las intersecciones del segmento, P0P1, con cada arista de un polígono convexo de recorte. Se usa la ecuación paramétrica de la línea para determinar los cuatro valores de t. Posteriormente, realizamos varias comparaciones para determinar cuáles de estos cuatro valores de tcorresponden a intersecciones reales. Sólo al tener los valores válidos de t es cuando calculamos los puntos de intersección.


Figura 16 - Clasificación de Puntos

Para encontrar el punto de intersección, el algoritmo de Cyrus-Beck se basa en la ecuación paramétrica de una línea recta. Como tenemos dos segmentos: P0P1 y la arista, que yacen sobre dos líneas rectas, tenemos dos ecuaciones paramétricas que las representan. Al resolver estas dos ecuaciones, averiguaremos un valor de t con el que determinaremos el punto de intersección. La ecuación paramétrica es la siguiente:
P(t) = P0 + (P1-P0) t
donde t queda comprendido en el intervalo [0,1]. Esto implica que P(0) = P0 y P(1) = P1.

Para determinar el valor exacto de t, haremos uso del producto escalar entre dos vectores. El primer vector será el vector normal de la arista que apunta hacia fuera de la región de recorte. El segundo vector será uno que coincide con la arista. Crearemos este segundo vector a partir de un punto cualquiera, Pa, en la arista, y el punto P(t) que es común a la línea del segmento, P0P1, que queremos recortar y a esta misma arista. El cálculo es el siguiente:
N⋅(P(t) - Pa)

Si el producto escalar es positivo, entonces el punto, P(t), yace fuera de la región de recorte y si es negativo, entonces P(t) yace dentro. Como nos interesa averiguar el punto común que yace en la arista y por tanto el borde de la región de recorte, el producto escalar debe ser 0; véase la figura 16. Esto significa que debemos resolver la siguiente ecuación:
N⋅(P(t) - Pa) = 0

Sustituimos P(t) por su definición,
N⋅(P0 + (P1-P0) t - Pa) = 0

Aplicamos la propiedad distributiva,
N⋅(P0 - Pa) + N⋅(P1-P0) t = 0

Dejemos que D = P1-P0, que es el vector de P0 a P1 del segmento. Ahora resolvemos para t:
     N⋅(P0 - Pa) 
t = ------------
        -N⋅D

Este valor de t es válido solamente si el denominador no es cero. Esto significa que las siguientes condiciones deben mantenerse:

N ≠ 0
 El vector normal jamás debería ser nulo. Si ocurriere lo contrario, entonces se trata de un error.
D ≠ 0
 Esto significa que P0 ≠ P1. Si ocurriere lo contrario, entonces se trata de un solo punto y no de un segmento.
N⋅D ≠ 0
 Esto significa que la arista y el segmento no son paralelos. Si fueren paralelos, entonces no existe ninguna intersección.


El algoritmo calcula cada intersección entre el mismo segmento y cada arista de la región de recorte. Para este cálculo, determinamos el vector normal de cada arista y un punto arbitrario que yace en una arista. Lo más sencillo y práctico es elegir un punto conocido como un vértice que forma la arista. Posteriormente, debemos elegir cuáles de los cuatro valores de t corresponden a intersecciones internas a la región de recorte. Primeramente, comprobamos que los valores de t quedan comprendidos en el intervalo [0,1]. Si no se cumple esta condición, entonces significa que el punto de intersección no yace en el segmento y por tanto descartamos tal valor de t. Como cada cálculo que hacemos se basa en encontrar las intersecciones entre dos líneas rectas, esto significa que un valor de t puede indicar una intersección que no yace en una arista.

En algunos casos, no es tan sencillo determinar si una intersección yace dentro o fuera de la región de recorte. Podríamos realizar varias comprobaciones, pero entonces ralentizaríamos el algoritmo y no conseguiríamos una ventaja sobre otros algoritmos, como el de Cohen-Sutherland. Analizando varios casos de intersecciones entre segmentos y aristas, descubrimos un comportamiento que nos servirá para discriminar tal intersección y así determinar si es aceptada o rechazada. Esta técnica se basa en denominar las intersecciones como potencialmente entrantes (PE) y potencialmente salientes (PS) de la región de recorte. Si cruzamos una arista, yendo desde los extremos P0 a P1, que nos hace entrar en la mitad del plano, entonces la intersección lleva la etiqueta de PE. Si por el contrario, al cruzar la arista, nos hace salir del plano interior creado por la arista, entonces tal intersección es PS. Con esta clasificación, notamos que dos intersecciones interiores a la región de recorte tienen etiquetas opuestas.

Para clasificar si una intersección es PE o PS, podemos basarnos en el ángulo formado por los vectores del segmento P0P1 y el normal de una arista. Si el ángulo es mayor de 90°, entonces es PE; y si es menor de 90°, entonces es PS. No requerimos calcular el ángulo, ya que esta información está contenida en el producto escalar del vector normal N y el vector D del segmento P0P1:

N⋅D < 0
 ⇒ ángulo > 90°
 ⇒ PE
N⋅D > 0
 ⇒ ángulo < 90°
 ⇒ PS


Vemos que N⋅D es el denominador en el cálculo de t para determinar el punto de intersección. Por lo tanto, podemos clasificar el valor de t mientras lo calculamos.

Nos interesa encontrar aquellos puntos de intersección etiquetados de PE a PS. El punto de intersección con la etiqueta PE que sea interior a la región de recorte es aquél que tenga el mayor valor de t; lo llamaremos, te. El punto interior que tenga la etiqueta PS es aquél que tenga el menor valor de t, al cual lo llamaremos, ts. El segmento cruzante se define en el intervalo [te,ts]. Como ya mencionamos anteriormente, los valores de t deben estar comprendidos en el intervalo de [0,1]. Por consiguiente, la cota inferior de te = 0 y la cota superior de ts = 1. Si te > ts, entonces el segmento no es interior a la región de recorte y por tanto lo rechazamos. Al final, obtendremos los valores correctos de te y ts para así calcular las coordenadas x e y de los puntos de intersección.


Figura 17 - Ejemplo de Tres Casos

Este algoritmo tiene la ventaja de no estar basado en un bucle comparado con el algoritmo de Cohen-Sutherland.



Ejemplo


Descripción


Figura 18 - Ejemplo

Queremos mostrar un segmento, pero únicamente dentro de nuestra vista representada por un triángulo descrito con los siguientes vértices: \. Tenemos un segmento definido por la siguiente pareja de sus puntos extremos:
A = (-3,-1) y B = (1,1)


Solución


Figura 19 - Calculamos Normales

Calculamos los vectores normales - sin normalizar - de todas las aristas de nuestro triángulo de recorte:
v12 = (  3, -3 ) - (  2,  3 ) = (  1, -6 )  ⇒  N1 = (  6,  1 )
v23 = ( -4, -2 ) - (  3, -3 ) = ( -7,  1 )  ⇒  N2 = ( -1, -7 )
v31 = (  2,  3 ) - ( -4, -2 ) = (  6,  5 )  ⇒  N3 = ( -5,  6 )



Figura 20 - Primera Intersección

Calculamos el valor de t para la intersección de la línea que contiene el segmento AB con la línea generada por v1v2. Primero calculamos y comprobamos el denominador del cálculo de t por si es 0 y por tanto el segmento es paralelo a la arista:
N1⋅D = (  6,  1 ) ⋅ (  4,  2 ) = 6*4 + 1*2 = 26

Como el denominador no es 0, calculamos el valor de t para ts ya que el denominador es positivo por lo que se trata de un punto potencialmente saliente:
     (  6,  1 ) ⋅ (( -3, -1 ) - (  2,  3 ))
t = ---------------------------------------- 
                     -26

     (  6,  1 ) ⋅ ( -5, -4 )
  = -------------------------
               -26

     6*(-5) + 1*(-4)     -34
  = ----------------- = ----- = 1,3077
            -26          -26

Como el valor de t no queda comprendido en [0,1], lo descartamos. El valor de ts sigue siendo 1, en estos momentos, ya que queremos el menor valor para ts. Esto es correcto, ya que 1 < 1,3077.



Figura 21 - Segunda Intersección

Calculamos el valor de t para la intersección de la línea que contiene el segmento AB con la línea generada por v2v3. Primero calculamos y comprobamos el denominador del cálculo de t por si es 0 y por tanto el segmento es paralelo a la arista:
N2⋅D = ( -1, -7 ) ⋅ (  4,  2 ) = (-1)*4 + (-7)*2 = -18

Como el denominador no es 0, calculamos el valor de t parate ya que el denominador es negativo por lo que se trata de un punto potencialmente entrante:
     ( -1, -7 ) ⋅ (( -3, -1 ) - (  3, -3 ))
t = ---------------------------------------- 
                   -(-18)

     ( -1, -7 ) ⋅ ( -6,  2 )
  = -------------------------
               18

     (-1)*(-6) + (-7)*2     -8
  = -------------------- = ---- = -0,4444
             18             18

Como el valor de t no queda comprendido en [0,1], lo descartamos. El valor de te sigue siendo 0, en estos momentos, ya que queremos el mayor valor para te. Esto es correcto, ya que -0,4444 < 0.



Figura 22 - Tercera Intersección

Calculamos el valor de t para la intersección de la línea que contiene el segmento AB con la línea generada por v3v1. Primero calculamos y comprobamos el denominador del cálculo de t por si es 0 y por tanto el segmento es paralelo a la arista:
N3⋅D = ( -5,  6 ) ⋅ (  4,  2 ) = (-5)*4 + 6*2 = -8

Como el denominador no es 0, calculamos el valor de t para te ya que el denominador es negativo por lo que se trata de un punto potencialmente entrante:
     ( -5,  6 ) ⋅ (( -3, -1 ) - ( -4, -2 ))
t = ---------------------------------------- 
                     -(-8)

     ( -5,  6 ) ⋅ (  1,  1 )
  = -------------------------
                8

     (-5)*1 + 6*1     1
  = -------------- = --- = 0,125
           8          8

Como el valor de t queda comprendido en [0,1], actualizamos el valor dete para que sea 0,125, ya que 0,125 > 0, su valor anterior.



Figura 23 - Solución Final

Al terminar de inspeccionar todas las aristas del triángulo de recorte, comprobamos los valores de te y ts, ya que son las cotas inferior y superior, respectivamente, de un intervalo: te ≤ ts. En nuestro caso, obtenemos el siguiente intervalo válido: [0,125, 1].

Pasamos a calcular los nuevos puntos extremos para el segmento recortado:
A′ = ( -3, -1 ) + (( -3, -1 ) - ( 1, 1 )) * 0,125
   = ( -3, -1 ) + (  4,  2 ) * 0,125
   = ( -3 + 0,5, -1 + 0,25 )
   = ( -2,5, -0,75 )

B′ = ( -3, -1 ) + (( -3, -1 ) - ( 1, 1 )) * 1
   = ( -3, -1 ) + ( 4, 2 ) * 1
   = ( -3 + 4, -1 + 2 )
   = (  1,  1 )





Algoritmo

Tenemos varias funciones a destacar para implementar este algoritmo:

La función principal es Recortar(), la cual aceptará tres parámetros: dos puntos extremos, P y Q, donde cada uno contiene las coordenadas (x,y)que representan el segmento a recortar, y un polígono convexo, R, que contiene una lista de vértices \ en el sentido de las agujas del reloj para representar la región de recorte.
La función Recortar_Punto() se invocará para tratar el segmento como un único punto. Básicamente, se realizará el algoritmo explicadoanteriormente.
La función Calcular_Normal() es necesaria si no tenemos una lista de normales de las aristas del polígono de recorte, R.
Las funciones mayor() y menor() determinan cuáles de los dos parámetros dados son el mayor y el menor, respectivamente.
La función Calcular_Punto() sirve para obtener un punto en una línea recta según la ecuación paramétrica. Esta función requiere dos puntos extremos de la línea y el parámetro, t.




Para Recortar(), el algoritmo es:
booleano Recortar( ref Punto P, ref Punto Q, PolígonoConvexo R )

Si P = Q, entonces


resultado ← Recortar_Punto( P, R )
Terminar( resultado )

Si no, entonces


te ← 0
ts ← 1
Para cada vértice, vi, de R, hacer


N ← Calcular_Normal( vi+1 - vi )
numerador ← N⋅(P - vi)
denominador ← -N⋅(Q - P)
Si denominador ≠ 0, entonces


t ← numerador / denominador
Si denominador < 0, entonces


te ← mayor( te, t )

Si no, entonces


ts ← menor( ts, t )


Si no, compruebe que numerador > 0,


Terminar( falso )


Si te > ts, entonces


Terminar( falso )

Si no, entonces


Pe ← Calcular_Punto( P, Q, te )
Ps ← Calcular_Punto( P, Q, ts )
P ← Pe
Q ← Ps


Terminar( verdadero )



Si el segmento es aceptado, entonces el algoritmo terminará con el valor booleano verdadero. Además, los extremos, P y Q, contendrán los valores recortados por este mismo algoritmo. Si el segmento es rechazado, el algoritmo termina con el valor falso, lo cual implica que los parámetros, P y Q, no son alterados. En cualquier caso, se pasan estos dos parámetros por referencia, para así poder modificarlos.

A continuación presentamos el algoritmo para Calcular_Normal(), el cual requiere un vector, v:
Vector Calcular_Normal( Vector v )

vector N ← ( -v.y, v.x )
N ← Normalizar( N )
Terminar( N )



Este cálculo nos sirve si el vector dado sigue el convenio de ir en el sentido de las agujas del reloj. La función Normalizar() se refiere a la operación vectorial para convertir un vector en unitario. Sin embargo, en este algoritmo, no es necesario normalizar estos vectores normales, por lo que podemos conseguir una optimización para el algoritmo de Cyrus-Beck.

Por último, tenemos el algoritmo para Calcular_Punto(), que requiere dos puntos, P y Q, y un valor para t:
Punto Calcular_Punto( Punto P, Punto Q, real t )

Punto A ← P + (Q-P) t
Terminar( A )