Geometria Computacional – Algoritmo de Cyrus-Beck

Neste artigo, vou apresentar uma introdução ao algoritmo de Cyrus-Back para o campo da geometria computacional.

Os cientistas Cyrus e Beck seguiram uma aproximação ao problema do recorte de segmentos de reta no qual conceberam um algoritmo que permite o recorte de um segmento de reta por qualquer polígono convexo, em 2D, ou de um segmento de reta em 3D por qualquer poliedro convexo.

Consideremos a representação paramétrica de um segmento de reta. Como pode observar-se na figura abaixo, existem sempre 4 valores de t que resultam da interseção da reta contendo o segmento do reta a recortar com as 4 arestas do retângulo de recorte.

Baixe o código fonte do projeto aqui: Efetue o Login

Recorte paramétrico: pontos de interseção da reta contendo o segmento de reta a recortar com o retângulo de recorte.

Através da classificação de cada um dos pontos de interseção é possível saber quais os valores de t dos vértices do segmento recortado, caso existam (no exemplo da figura seriam os vértices correspondentes a t2 e t3).

Interseção entre linhas

Considere-se novamente a representação paramétrica de um segmento de reta:

O polígono de recorte é definido pela enumeração dos seus lados (ou vértices) segundo a convenção da circulação direta (sentido contrário ao dos ponteiros do relógio). Para cada aresta define-se a normal, Ni, orientada para o exterior do polígono.

Seja “Pei” um ponto arbitrário pertencente a uma aresta do polígono de recorte. Como pode observar-se na figura abaixo, podemos definir três vectores que partem desse ponto, um deles seguindo a fronteira e os outros dois para o espaço interior e para o espaço exterior ao polígono convexo de recorte.

Algoritmo de Cyrus-Beck: produtos internos dos três vectores com a normal exterior ao polígono de recorte.

O produto interno da normal à aresta do polígono com cada um destes três vectores é, respectivamente, nulo, negativo e positivo. O sinal do produto interno permite determinar se um ponto do segmento de reta coincide com a aresta ou é um ponto interior ou exterior do polígono de recorte.

Para determinar o ponto de interseção do segmento de reta com a aresta que estamos a considerar é necessário obter o valor do parâmetro no ponto de interseção. Como vimos atrás, no ponto de interseção o produto interno da normal à aresta com o vector que une o ponto arbitrário “Pei” com o ponto de interseção é nula, ou seja,

Substituindo P(t) pela equação paramétrica da reta que contem o segmento de reta e explicitando o parâmetro t, obteremos sucessivamente:

e, definindo D=P1-P0,

Para ser possível calcular t é necessário que:

• N_i≠0, o que é sempre verdadeiro;
• P1 e P0 não sejam coincidentes, o que é verdadeiro para qualquer segmento de reta de comprimento não nulo;
• Ni * D ≠ 0, o que aconteceria se o segmento de reta a recortar fosse paralelo à aresta do polígono;

Não existe interseção se o valor de t for inferior a 0 ou superior a 1.

Tendo sido obtido um modo de calcular se um ponto está no exterior ou no interior do polígono convexo de recorte e um modo de calcular o valor de t na interseção está em condições de descrever o algoritmo paramétrico.

Descrição do Algoritmo

O algoritmo de Cyrus-Beck aplica-se sucessivamente a cada uma das arestas do polígono convexo. Para cada uma das arestas, escolhe-se, por exemplo, o primeiro vértice como “Pei”, calcula-se a normal exterior, Ni e o valor de t no ponto de interseção.

É ainda necessário identificar se o primeiro vértice do segmento de reta está no exterior e o 2º no interior, situação em que a interseção é definida como Potencialmente de Entrada, PE, ou vice-versa em que a interseção é definida como Potencialmente de Saída, PS.

Para efetuar essa identificação poderia classificar-se cada vértice do segmento de reta, tal como acontecia no algoritmo da Cohen-Sutherland. No entanto, é possível efetuar a classificação com base no ângulo existente entre o vector P0P1 e o vector Ni. Se o ângulo for inferior a 90º trata-se de uma interseção do tipo PS, se não, será do tipo PE. Esta informação está contida no sinal do produto interno entre aqueles dois vectores e que é necessário para o cálculo do próprio valor de t.

Estando calculados todos os valores t dos pontos de interseção, e devidamente classificados como sendo do tipo PE ou PS, estamos em condições de identificar exatamente os valores de t dos vértices do segmento de reta já recortado.

Identifica-se a interseção do tipo PE que corresponde ao maior t e a interseção do tipo PS que corresponde ao menor t e comparam-se os dois valores de t.

Se o valor de t da interseção do tipo PE for superior ao valor de t da interseção do tipo PS, então todo o segmento de reta pode ser rejeitado.

Se o valor de t da interseção do tipo PE for inferior ao valor de t da interseção do tipo PS, então usam-se os dois valores para obter as coordenadas dos vértices do segmento de reta recortado.

A figura abaixo apresenta um polígono convexo de recorte com 5 lados e um segmento de reta a recortar (P0P1).

Algoritmo de Cyrius-Beck: segmento de reta a recortar com identificação do tipo (PE ou PS) das suas intersecções com o polígono de recorte.

Do cálculo das intersecções resultam 4 valores de t, sendo dois do tipo PE e dois do tipo PS. Como o segmento de reta é paralelo ao lado inferior direito do polígono, o 5º valor de t não foi calculado uma vez que N_i * D = 0.

Escolhendo-se a interseção do tipo PE com maior valor de t e a interseção do tipo PS com o menor valor de t obtêm-se os valores paramétricos de cada um dos vértices do segmento recortado.

Considere-se agora o exemplo apresentado pela figura abaixo. Aplicando o algoritmo, são calculados 5 valores de t dos quais só são classificados os três valores representados na figura uma vez que os restantes têm valores de t inferiores a 0 ou superiores a 1.

Algoritmo de Cyrus-Beck: segmento de reta rejeitado (o maior t das intersecções de tipo PE é superior ao menor t de tipo PS).

Como pode observar-se, a interseção do tipo PS tem um valor de t inferior ao maior valor de t de uma interseção do tipo PE pelo que se conclui que o segmento de reta deve ser rejeitado.

Pseudocódigo (Algoritmo)

Pré-calcular Ni e selecionar “Pei” para cada aresta.

Para cada segmento de reta

Se P1 = P0 então

Segmento de reta degenerado: recortar como ponto

te = 0; tl= 1;

Para cada aresta de recorte

Se Ni . D ¹ 0 então /* ignorar arestas paralelas ao segmento */

calcular t; /* de interseção do segmento com aresta */

usar sinal de Ni . D para caracterizar PE ou PS;

Se PE então tE = max(tE, t);

Se PS então tS = min(tS, t);

FimSe

Se tE > tS então devolve nulo;

Senão devolve [P(tE), P(tS)];

Código em C usando OpenGL

#include <windows.h>   
#include <gl/Gl.h>
#include <gl/glut.h>
#include <iostream>

//********* Estrutura de Dados

struct GLintPoint
{ GLint x,y;
};

struct GLfloatPoint
{ GLfloat x,y;
};

const int MAX = 100;

class GLintPointArray
{ 
public:  
  int num;
  GLintPoint pt[MAX];
};

class GLfloatPointArray
{
public:
  int num;
  GLfloatPoint pt[MAX];
};

//********** Rotinas de suporte
typedef GLfloat colorType[3];

void drawDot (GLint x, GLint y, GLfloat r, GLfloat g, GLfloat b)
{ glColor3f(r,g,b);

  glBegin (GL_POINTS);
  glVertex2i (x,y);
  glEnd();
}

void drawIntPolygon (GLintPointArray P, colorType c)
{ glColor3fv (c);

  glBegin(GL_LINE_LOOP);
    for (int i=0;i<P.num;i++)

	glVertex2i (P.pt[i].x,P.pt[i].y);

  glEnd();
}

void drawFloatPolygon (GLfloatPointArray P, colorType c)
{ glColor3fv (c);

  glBegin(GL_LINE_LOOP);
   for (int i=0; i < P.num; i++)

     glVertex2f (P.pt[i].x,P.pt[i].y);

  glEnd();
}

//**************** Inicialização do programa 
 void myInit(void)
 {
    glClearColor(0.0,0.0,0.0,0.0);  // Coloca fundo branco

	glColor3f (0.0f,0.0f,0.0f);    
    glMatrixMode(GL_PROJECTION); 
    glLoadIdentity();

    gluOrtho2D(0.0, 640.0, 0.0, 480.0);
}

//**************  Subrotinas de desenho.

const int MaxN = 100;
typedef GLfloatPoint Normals[MaxN];

void buildintPolygon (GLintPointArray &P)
{ P.num = 4; 
  P.pt[0].x = 50; P.pt[0].y =  100; 
  P.pt[1].x = 150; P.pt[1].y = 200; 
  P.pt[2].x = 100; P.pt[2].y = 300; 
  P.pt[3].x = 10; P.pt[3].y =  200; 
}

void buildfloatPolygon (GLfloatPointArray &P)
{ P.num = 4;

  P.pt[0].x = 50; P.pt[0].y = 100; 
  P.pt[1].x = 150; P.pt[1].y = 200;

  P.pt[2].x = 100; P.pt[2].y = 300;

  P.pt[3].x = 10; P.pt[3].y = 200; 
}

float DotProduct (GLfloatPoint v1, GLfloatPoint v2)
{ 
   return v1.x*v2.x + v1.y*v2.y;
}

void CalcNormals (GLfloatPointArray p, Normals & n)
{ int i,j,k;

  GLfloatPoint v;

  for (i = 0; i < p.num; i++)
  { j = (i+1)%p.num;

    k = (i+2)%p.num;
    n[i].x = -(p.pt[j].y - p.pt[i].y)/(p.pt[j].x - p.pt[i].x);

    n[i].y = 1.0; 
    v.x = p.pt[k].x - p.pt[i].x;

    v.y = p.pt[k].y - p.pt[i].y;

    if (DotProduct (n[i],v) > 0)    // 

    { n[i].x *= -1;
      n[i].y  = -1;
    }
  }
}

void CBClip (GLfloatPoint p1, GLfloatPoint p2, Normals n, GLfloatPointArray p, 
             bool & visible, GLfloatPoint & rp, GLfloatPoint & q)
{ float t1,t2,t,numer,den;

  GLfloatPoint dirV,F;          // Vetores
  int i;

  t1 = 0.0;

  t2 = 1.0;
    // Calcula a direção do vetor.
  dirV.x = p2.x - p1.x;

  dirV.y = p2.y - p1.y;

  visible = true;
  i = 0;

  while ( (i < p.num) && visible)
  { F.x = p1.x - p.pt[i].x;

    F.y = p1.y - p.pt[i].y;

    numer  = DotProduct (n[i],F);
    den   =  DotProduct (n[i],dirV);

    if (den == 0.0)          
    { if (numer > 0.0)

         visible = false;   
    }
    else 
    { t = -(numer/den);

      if (den < 0.0)        
      { if (t <= 1.0)

          if (t > t1)
             t1 = t;
      }

      else if ( t >= 0.0)    
      { if (t < t2)

          t2 = t;
      }
    }
    i++;
  }
  if ( t1 <= t2)
  { rp.x = p1.x + t1*dirV.x;

    rp.y = p1.y + t1*dirV.y;

    q.x = p1.x + t2*dirV.x;

    q.y = p1.y + t2*dirV.y;
  }

  else
   visible = false;
}

//******************* Exibe na tela **********************
void myDisplay(void)
{     
  GLfloatPointArray Poly;

  colorType C = {1.0f,0.0f,1.0f};
  bool visible;

  GLfloatPoint p1 = {40,20},p2 = {200,400},cp1,cp2;

  Normals Nor;

  glClear(GL_COLOR_BUFFER_BIT);     // Limpa a tela 

  // 1. Cria o polígono
  buildfloatPolygon(Poly);

  drawFloatPolygon (Poly,C);

  // 2. Cria as linhas de interseção
  glColor3f (0.0f,3.0f,0.0f);

  glBegin(GL_LINES);
    glVertex2f (40,20);
    glVertex2f (200,400);

  glEnd();

  // 3. Corta a linha
  CalcNormals (Poly,Nor);
  CBClip (p1,p2,Nor,Poly,visible,cp1,cp2);

  if (visible)
  {	  
	  glColor3f (1.0f,1.0f,0.0f);

	  glBegin(GL_LINES);
	  glVertex2f (cp1.x,cp1.y);

	  glVertex2f (cp2.x,cp2.y);
	  glEnd();
  }

  glFlush();
}

//******************* Principal **********************
// Rotinas do OpenGL
void main(int argc, char** argv)
{

  glutInit(&argc, argv);          
  glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); 
  glutInitWindowSize(640,480);     
  glutInitWindowPosition(100, 150); 
  glutCreateWindow("Cyrus-Beck"); 
  glutDisplayFunc(myDisplay);     
  myInit();          
  glutMainLoop();      
}