Intro
Este documento ilustra la implementación de un algoritmo geométrico sencillo (potencialmente utilizable en la programación de juegos) en los lenguajes C, C++ y Java, efectuando una comparación de la performance relativa. La audiencia potencial son los desarrolladores que utilizan cualquiera de estos lenguajes y desean tener una idea acerca de la performance esperada en las implementaciones que realizan una gran cantidad de cálculos aritméticos.
Revisiones
-
1.0 2007-04-06 Primera versión
|
Note
|
Este documento puede ser leído en compañía del artículo "Ventajas y Desventajas entre C, C++ y Java", disponible en la misma dirección, el cual proporciona diversos aspectos teóricos importantes. El presente material va directamente a la implementación de un caso práctico. |
Planteamiento del Caso
En esta sección describiremos rápidamente lo que deseamos que haga nuestro algoritmo; sin embargo, los detalles matemáticos del mismo serán detallados en un apéndice al final, puesto que no son esenciales para seguir con la discusión de Software.
El problema que deseamos resolver surge con mucha frecuencia en videojuegos: la detección de colisiones entre objetos de dos dimensiones (por ejemplo, un superheroe que es atacado por un mounstro.) Una manera de efectuar este análisis consiste en asumir que los elementos gráficos del juego consisten de un perímetro poligonal; más precísamente, un conjunto de puntos o vértices.
Con estos conceptos, podemos definir una colisión como el caso en que algún vértice de un polígono se encuentre en la región interior de otro polígono.
Evidentemente este problema se puede subdividir en otro un poco más sencillo, a saber, el determinar si un punto está en la región interior de un polígono.
A fin de simplificar la exposición, asumiremos también que los polígonos son convexos (no tienen cavidades o entradas o agujeros.)
El apéndice "A" proporciona una explicación intuitiva pero detallada del algoritmo que será implementado a continuación.
Preparación de datos de Prueba
En tanto deseamos hacer una comparación de implementaciones, es conveniente contar con idénticos datos de prueba. En particular, crearemos 1000 puntos aleatoriamente, los cuales serán analizados contra 10000 polígonos generados aleatoriamente. A fin de simplificar el proceso, utilizaremos una rutina estándar de generación de números pseudoaleatorios de modo tal que en las tres implementaciones se disponga de los mismos datos (y por supuesto, debemos obtener el mismo resultado.)
La siguiente función (en lenguaje C) ha sido tomada de [MN] y se anuncia como un excelente generador de números pseudoaleatorios dada su simplicidad:
#define IA 16807
#define IQ 127773
#define IR 2836
#define AM (1.0/ALEATORIO_MAX)
#define MASK 123459876
#define ALEATORIO_MAX 2147483647L
long semilla=170374;
long aleatorio(void)
{
long k,ans;
semilla^=MASK;
k=semilla/IQ;
semilla=IA*(semilla-k*IQ)-IR*k;
if(semilla<0)
semilla+=ALEATORIO_MAX;
ans=semilla;
semilla^=MASK;
return ans;
}
Nótese que su "traducción" a C++ y Java es trivial.
Esta función será utilizada para "preparar" los puntos de prueba y los polígnos de prueba para el análisis principal. Nótese que no consideraremos la performance de esta etapa de preparación.
Variantes del Algoritmo
Tal como se explica en el apéndice A, el algoritmo básico (al que llamaremos "variante 1") consiste en dividir a los polígonos en un conjunto de triángulos y analizar si el punto está al interior de cualquiera de dichos triángulos.
La "variante 2" consiste en evitar este análisis (que es relativamente pesado) para el caso (afortunado) en que todos los puntos del polígono se ubiquen en ciertos cuadrantes con respecto al punto de prueba. Esto consigue acelerar el proceso en un porcentaje considerable.
La "variante 3" intenta evitar completamente los anteriores análisis para el caso (típicamente frecuente) en que el punto está bastante alejado del polígono.
Con estas variantes obtendremos tres resultados distintos para cada una de los lenguajes de programación. El código que se presenta más adelante implementa la "variante 3"; a partir de ésta, deberá ser muy sencillo para el lector comentar el código necesario a fin de obtener las otras dos variantes.
Implementación de Objetos
Los Puntos y los Polígonos son candidatos naturales para la creación de clases. Para el caso del lenguaje C, crearemos dos estructuras y un conjunto de funciones que hacen un trabajo similar al que harían los métodos correspondientes.
Crearemos también una clase auxiliar utilitaria (llamada Util en C++ y AGUtil en Java) para agregar las rutinas de números pseudoaleatorios y de cálculo de tiempos. Para el C simplemente crearemos un archivo adicional.
Los polígonos de prueba serán heptágonos regulares, lo cual considero que puede ser un buen punto de partida para figuras más complejas [1].
Las implementaciones serán presentadas en el siguiente orden: 1) C++, 2) C ,y 3) Java. El motivo radica en que C++ puede considerarse un "punto medio" a partir de lo cual es fácil comprender las otras dos implementaciones.
Eclipse
Las tres implementaciones han sido elaboradas con ayuda del IDE Eclipse 3.2. Como es sabido, hasta esta versión Eclipse no proporciona automáticamente soporte de C o C++, por lo que le agregamos el CDT (Eclipse C/C++ Development Tools [2].)
Implementación en Lenguaje C++
Rutina principal main
Nótese que proporcionamos el código necesario para crear heptágonos y triángulos (pero sólo el primero se utiliza.)
#include <stdio.h>
#include <math.h>
#include "Point.h"
#include "Polygon.h"
#include "Util.h"
using namespace std;
static void setTriangle(Polygon &t);
static void setHepta(Polygon &t);
int main()
{
int z,j,interior=0;
const int maxpol=10000,maxpnt=1000;
Polygon pol[maxpol];
for(z=0;z<maxpol;z++)
setHepta(pol[z]);
// setTriangle(pol[z]);
Point p[maxpnt];
for(z=0;z<maxpnt;z++)
{
p[z].setX(Util::aleatorio_entero(0,800));
p[z].setY(Util::aleatorio_entero(0,600));
}
double t1,t2;
t1=Util::getTime();
for(j=0;j<maxpnt;j++)
for(z=0;z<maxpol;z++)
{
if(pol[z].isInterior(p[j]))
interior++;
}
t2=Util::getTime();
double ratio=(interior*100.0)/(maxpol*maxpnt);
printf("%d interiores de %d = %g %% (t=%g)\n",
interior,maxpol*maxpnt,ratio,t2-t1);
return 0;
}
static void setTriangle(Polygon &t)
{
Point p1(Util::aleatorio_entero(0,800),
Util::aleatorio_entero(0,600));
Point p2(Util::aleatorio_entero(0,800),
Util::aleatorio_entero(0,600));
Point p3(Util::aleatorio_entero(0,800),
Util::aleatorio_entero(0,600));
t.addPoint(p1);
t.addPoint(p2);
t.addPoint(p3);
}
static void setHepta(Polygon &t)
{
double radius=20+Util::aleatorio_entero(0,800)/8.0;
double h=Util::aleatorio_entero(0,800);
double k=Util::aleatorio_entero(0,600);
int z;
double PI=3.1415926535;
for(z=0;z<7;z++)
t.addPoint(Point(h+radius*cos(2*z*PI/7),
k+radius*sin(2*z*PI/7)));
}Archivos Cabecera
Téngase en cuenta que para los puntos y los polígonos se están proporcionando más operaciones de las estrictamente necesarias, las cuales pueden ser comentadas si se desea.
La clase del Punto
Nótese que las coordenadas corresponden a valores "double", los cuales se han encapsulado en métodos "get/set".
#ifndef POINT_H_
#define POINT_H_
class Point
{
private:
double x;
double y;
public:
Point();
Point(double a,double b);
Point(const Point &a,const Point &b);
Point(const Point &a);
double getX() const;
double getY() const;
void setX(double n);
void setY(double n);
void print() const;
double distance2(const Point &p) const;
double distance(const Point &p) const;
int quadrant(const Point &p) const;
bool insideTriangle(const Point &p1,
const Point &p2, const Point &p3) const;
double prod(const Point &p) const;
Point &add(const Point &p);
Point &sub(const Point &p);
Point &ort();
virtual ~Point();
};
#endif /*POINT_H_*/La clase del Polígono
El polígono hace uso de un vector de la STL:
#ifndef POLYGON_H_
#define POLYGON_H_
#include <vector>
#include "Point.h"
class Polygon
{
private:
double diameter;
std::vector<Point> pt;
void findDiameter();
public:
Polygon();
int getCount() const;
const Point &getPoint(int n) const;
void setPoint(int n,const Point &p);
void addPoint(const Point &p);
void print();
bool isInterior(const Point &p);
virtual ~Polygon();
};
#endif /*POLYGON_H_*/La clase utilitaria
La clase utilitaria implementa las funciones como métodos estáticos.
#ifndef UTIL_H_
#define UTIL_H_
class Util
{
private:
static long aleatorio();
static long semilla;
public:
static long aleatorio_entero(long a,long b);
static double getTime();
};
#endif /*UTIL_H_*/Implementaciones
Implementación del Punto
El único método complicado corresponde a insideTriangle(). Ver el apéndice A para más información.
#include "Point.h"
#include <stdio.h>
#include <math.h>
using namespace std;
Point::Point()
{
x=y=0;
}
Point::Point(double a,double b)
{
x=a; y=b;
}
Point::Point(const Point &a,const Point &b)
{
x=b.getX()-a.getX();
y=b.getY()-a.getY();
}
Point::Point(const Point &a)
{
x=a.getX();
y=a.getY();
}
Point::~Point()
{
}
double Point::getX() const
{
return x;
}
double Point::getY() const
{
return y;
}
void Point::setX(double n)
{
x=n;
}
void Point::setY(double n)
{
y=n;
}
void Point::print() const
{
printf("[%g,%g]",x,y);
}
double Point::distance2(const Point &p) const
{
double dx=x-p.getX();
double dy=y-p.getY();
return dx*dx+dy*dy;
}
double Point::distance(const Point &p) const
{
return sqrt(distance2(p));
}
int Point::quadrant(const Point &p) const
{
if(p.getX()>=x)
{
if(p.getY()>y)
return 1;
else
return 4;
}
else
{
if(p.getY()>y)
return 2;
else
return 3;
}
}
Point &Point::ort()
{
double t=y;
y=x;
x=-t;
return *this;
}
bool Point::insideTriangle(const Point &A,
const Point &B, const Point &C) const
{
Point PA(*this,A);
Point PB(*this,B);
Point PC(*this,C);
Point PAB(PA);PAB.add(PB);
Point PBC(PB);PBC.add(PC);
Point PCA(PC);PCA.add(PA);
Point PAN(B,C);PAN.ort();
Point PBN(C,A);PBN.ort();
Point PCN(A,B);PCN.ort();
Point AB(A,B);
Point BC(B,C);
Point CA(C,A);
if(AB.prod(PAN)*PB.prod(PAN)<0)
return false;
if(BC.prod(PBN)*PC.prod(PBN)<0)
return false;
if(CA.prod(PCN)*PA.prod(PCN)<0)
return false;
return true;
}
Point &Point::add(const Point &p)
{
x+=p.getX();
y+=p.getY();
return *this;
}
Point &Point::sub(const Point &p)
{
x-=p.getX();
y-=p.getY();
return *this;
}
double Point::prod(const Point &p) const
{
return x*p.getX()+y*p.getY();
}Implementación del Polígono y Variantes
El polígono implementa las tres variantes en el método isInterior():
#include "Polygon.h"
#include <stdio.h>
#include <math.h>
#include <exception>
using namespace std;
Polygon::Polygon():diameter(0),pt()
{
}
Polygon::~Polygon()
{
}
int Polygon::getCount() const
{
return pt.size();
}
const Point &Polygon::getPoint(int n) const
{
if(n<0 || n>=(int)pt.size())
throw exception();
return pt[n];
}
void Polygon::setPoint(int n,const Point &p)
{
if(n<0 || n>=(int)pt.size())
return;
pt[n]=p;
findDiameter();
}
void Polygon::addPoint(const Point &p)
{
pt.push_back(p);
findDiameter();
}
void Polygon::print()
{
int z;
for(z=0;z<getCount();z++)
getPoint(z).print();
printf(" D=%g\n",diameter);
}
void Polygon::findDiameter()
{
if(getCount()<3)
return;
int z,j;
double rn;
diameter=0;
for(z=0;z<getCount()-1;z++)
for(j=z+1;j<getCount();j++)
{
rn=getPoint(z).distance2(getPoint(j));
if(rn>diameter)
diameter=rn;
}
diameter=sqrt(diameter);
}
bool Polygon::isInterior(const Point &p)
{
if(getCount()<3)
throw exception();
if(p.distance2(getPoint(0))>diameter*diameter)
return false;
int z,current,Q=0;
for(z=0;z<getCount();z++)
{
current=p.quadrant(getPoint(z));
Q|=(1<<current);
}
if(Q==30) // 4 quadrants
return true;
if(Q!=28 && Q!=26 && Q!=22 && Q!=14 && Q!=10 && Q!=20) // not 3q nor oposed 2q
return false;
for(z=1;z<getCount()-1;z++)
if(p.insideTriangle(getPoint(0),
getPoint(z),getPoint(z+1)))
return true;
return false;
}Las variantes se pueden obtener fácilmente comentando el código que invoca al método distance2() (variante 2) y todo el código relacionado a los cuadrantes "Q" (variante 1.)
Implementación de Clase Utilitaria
A continuación se presenta la implementación de la clase utilitaria. Nótese que las rutinas de tiempo son estándares POSIX.
#include "Util.h"
#include <sys/time.h>
#include <time.h>
#define IA 16807
#define IQ 127773
#define IR 2836
#define AM (1.0/ALEATORIO_MAX)
#define MASK 123459876
#define ALEATORIO_MAX 2147483647L
long Util::semilla=170374;
long Util::aleatorio(void)
{
long k,ans;
semilla^=MASK;
k=semilla/IQ;
semilla=IA*(semilla-k*IQ)-IR*k;
if(semilla<0)
semilla+=ALEATORIO_MAX;
ans=semilla;
semilla^=MASK;
return ans;
}
long Util::aleatorio_entero(long minimo,long maximo)
{
return minimo+(int)((maximo-minimo+1.0)*aleatorio()/(ALEATORIO_MAX+1.0));
}
double Util::getTime()
{
struct timeval tv;
struct timezone tz;
gettimeofday(&tv,&tz);
return tv.tv_sec+tv.tv_usec/1000000.0;
}Implementación en Lenguaje C
Rutina principal main
Análogamente al caso C++, se proporcionan las implementaciones de triángulos y heptágonos. Notar que los métodos de la clase "Point" han sido reemplazados con funciones con prefijo "Point+" a fin de mantener la similitud con el caso C++.
#include <stdio.h>
#include <math.h>
#include "Point.h"
#include "Polygon.h"
#include "Util.h"
static void setTriangle(Polygon *t);
static void setHepta(Polygon *t);
int main()
{
int z,j,interior=0;
const int maxpol=10000,maxpnt=1000;
Polygon pol[maxpol];
for(z=0;z<maxpol;z++)
setHepta(&pol[z]);
// setTriangle(pol[z]);
Point p[maxpnt];
for(z=0;z<maxpnt;z++)
{
p[z].x=aleatorio_entero(0,800);
p[z].y=aleatorio_entero(0,600);
}
double t1,t2;
t1=getTime();
for(j=0;j<maxpnt;j++)
for(z=0;z<maxpol;z++)
{
if(Polygon_isInterior(&pol[z],&p[j]))
interior++;
}
t2=getTime();
double ratio=(interior*100.0)/(maxpol*maxpnt);
printf("%d interiores de %d = %g %% (t=%g)\n",
interior,maxpol*maxpnt,ratio,t2-t1);
return 0;
}
static void setTriangle(Polygon *t)
{
Point p1;Point_set_double(&p1,aleatorio_entero(0,800),
aleatorio_entero(0,600));
Point p2;Point_set_double(&p2,aleatorio_entero(0,800),
aleatorio_entero(0,600));
Point p3;Point_set_double(&p3,aleatorio_entero(0,800),
aleatorio_entero(0,600));
Polygon_set_zero(t);
Polygon_addPoint(t,&p1);
Polygon_addPoint(t,&p2);
Polygon_addPoint(t,&p3);
}
static void setHepta(Polygon *t)
{
double radius=20+aleatorio_entero(0,800)/8.0;
double h=aleatorio_entero(0,800);
double k=aleatorio_entero(0,600);
int z;
double PI=3.1415926535;
Polygon_set_zero(t);
for(z=0;z<7;z++)
{
Point toadd;
Point_set_double(&toadd,
h+radius*cos(2*z*PI/7),
k+radius*sin(2*z*PI/7));
Polygon_addPoint(t,&toadd);
}
}Archivos Cabecera
Declaración del Punto
Esto es muy similar al caso C++, con la diferencia que se han extraído los métodos como funciones globales:
#ifndef POINT_H_
#define POINT_H_
typedef struct
{
double x;
double y;
} Point;
void Point_set_zero(Point *p);
void Point_set_double(Point *p,double a,double b);
void Point_set_2_point(Point *p,
const Point *a,const Point *b);
void Point_set_point(Point *p,const Point *a);
void Point_print(const Point *p);
double Point_distance2(const Point *p1,const Point *p2);
double Point_distance(const Point *p1,const Point *p2);
int Point_quadrant(const Point *me,const Point *p);
void Point_ort(Point *p);
void Point_add(Point *me,const Point *p);
void Point_sub(Point *me,const Point *p);
double Point_prod(const Point *p1,const Point *p2);
int Point_insideTriangle(const Point *p,
const Point *A, const Point *B, const Point *C);
#endif /*POINT_H_*/Declaración del Polígono
La estructura Polygon tiene un detalle importante: almacena un puntero a estructuras Point. En tanto no existen rutinas de vectores en C, asignaremos la memoria dinámicamente (evidentemente, se pudo utilizar una librería auxiliar como "GLIB", pero esto nos alejaría mucho del tema.)
#ifndef POLYGON_H_
#define POLYGON_H_
#include "Point.h"
typedef struct
{
double diameter;
Point *pt;
int count;
} Polygon;
void Polygon_set_zero(Polygon *p);
int Polygon_getCount(Polygon *p);
const Point *Polygon_getPoint(const Polygon *p,int n);
void Polygon_findDiameter(Polygon *p);
void Polygon_setPoint(Polygon *me,int n,const Point *p);
void Polygon_addPoint(Polygon *me,const Point *p);
void Polygon_print(Polygon *p);
int Polygon_isInterior(const Polygon *me,const Point *p);
#endif /*POLYGON_H_*/Declaración de rutinas utilitarias
Finalmente, Util.h sólo contiene la declaración de las funciones a utilizarse externamente:
#ifndef UTIL_H_
#define UTIL_H_
long aleatorio_entero(long a,long b);
double getTime();
#endif /*UTIL_H_*/Implementaciones
Rutinas del Punto
Salvo la apariencia poco estética, esto es casi una reescritura de la versión C++:
#include "Point.h"
#include <stdio.h>
#include <math.h>
void Point_set_zero(Point *p)
{
p->x=p->y=0;
}
void Point_set_double(Point *p,double a,double b)
{
p->x=a; p->y=b;
}
void Point_set_2_point(Point *p,
const Point *a,const Point *b)
{
p->x=b->x-a->x;
p->y=b->y-a->y;
}
void Point_set_point(Point *p,const Point *a)
{
p->x=a->x;
p->y=a->y;
}
void Point_print(const Point *p)
{
printf("[%g,%g]",p->x,p->y);
}
double Point_distance2(const Point *p1,const Point *p2)
{
double dx=p1->x-p2->x;
double dy=p1->y-p2->y;
return dx*dx+dy*dy;
}
double Point_distance(const Point *p1,const Point *p2)
{
return sqrt(Point_distance2(p1,p2));
}
int Point_quadrant(const Point *me,const Point *p)
{
if(p->x>=me->x)
{
if(p->y>me->y)
return 1;
else
return 4;
}
else
{
if(p->y>me->y)
return 2;
else
return 3;
}
}
void Point_ort(Point *p)
{
double t=p->y;
p->y=p->x;
p->x=-t;
}
void Point_add(Point *me,const Point *p)
{
me->x+=p->x;
me->y+=p->y;
}
void Point_sub(Point *me,const Point *p)
{
me->x-=p->x;
me->y-=p->y;
}
double Point_prod(const Point *p1,const Point *p2)
{
return p1->x*p2->x+p1->y*p2->y;
}
int Point_insideTriangle(const Point *p,
const Point *A, const Point *B, const Point *C)
{
Point PA; Point_set_2_point(&PA,p,A);
Point PB; Point_set_2_point(&PB,p,B);
Point PC; Point_set_2_point(&PC,p,C);
Point PAB; Point_set_point(&PAB,&PA); Point_add(&PAB,&PB);
Point PBC; Point_set_point(&PBC,&PB); Point_add(&PBC,&PC);
Point PCA; Point_set_point(&PCA,&PC); Point_add(&PCA,&PA);
Point PAN; Point_set_2_point(&PAN,B,C);Point_ort(&PAN);
Point PBN; Point_set_2_point(&PBN,C,A);Point_ort(&PBN);
Point PCN; Point_set_2_point(&PCN,A,B);Point_ort(&PCN);
Point AB; Point_set_2_point(&AB,A,B);
Point BC; Point_set_2_point(&BC,B,C);
Point CA; Point_set_2_point(&CA,C,A);
if(Point_prod(&AB,&PAN)*Point_prod(&PB,&PAN)<0)
return 0;
if(Point_prod(&BC,&PBN)*Point_prod(&PC,&PBN)<0)
return 0;
if(Point_prod(&CA,&PCN)*Point_prod(&PA,&PCN)<0)
return 0;
return 1;
}Rutinas del Polígono
En la implementación del polígono hay que tener cuidado con el agregado de nuevos puntos: esto requiere asignar memoria dinámicamente, lo cual efectuamos con realloc(). En la primera pasada, el puntero debe contener NULL (de lo contrario el programa podría caer.)
#include "Polygon.h"
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
void Polygon_set_zero(Polygon *p)
{
p->pt=NULL;
p->diameter=0;
p->count=0;
}
int Polygon_getCount(Polygon *p)
{
return p->count;
}
const Point *Polygon_getPoint(const Polygon *p,int n)
{
if(n<0 || n>=(int)p->count)
abort();
return p->pt+n;
}
void Polygon_findDiameter(Polygon *p)
{
if(p->count<3)
return;
int z,j;
double rn;
p->diameter=0;
for(z=0;z < p->count-1;z++)
for(j=z+1;j < p->count;j++)
{
rn=Point_distance2(Polygon_getPoint(p,z),
Polygon_getPoint(p,j));
if(rn>p->diameter)
p->diameter=rn;
}
p->diameter=sqrt(p->diameter);
}
void Polygon_setPoint(Polygon *me,int n,const Point *p)
{
if(n<0 || n>=me->count)
return;
me->pt[n].x=p->x;
me->pt[n].y=p->y;
Polygon_findDiameter(me);
}
void Polygon_addPoint(Polygon *me,const Point *p)
{
me->pt=realloc(me->pt,sizeof(Point)*(me->count+1));
me->count++;
Polygon_setPoint(me,me->count-1,p);
}
void Polygon_print(Polygon *p)
{
int z;
for(z=0;z<p->count;z++)
Point_print(Polygon_getPoint(p,z));
printf(" D=%g\n",p->diameter);
}
int Polygon_isInterior(const Polygon *me,const Point *p)
{
if(me->count<3)
abort();
if(Point_distance2(p,Polygon_getPoint(me,0))
> me->diameter*me->diameter)
return 0;
int z,current,Q=0;
for(z=0;z< me->count;z++)
{
current=Point_quadrant(p,Polygon_getPoint(me,z));
Q|=(1<<current);
}
if(Q==30) // 4 quadrants
return 1;
if(Q!=28 && Q!=26 && Q!=22 && Q!=14 && Q!=10 && Q!=20) // not 3q nor oposed 2q
return 0;
for(z=1;z < me->count-1;z++)
if(Point_insideTriangle(p,Polygon_getPoint(me,0),
Polygon_getPoint(me,z),Polygon_getPoint(me,z+1)))
return 1;
return 0;
}Rutinas Utilitarias
La implementación de las funciones utilitarias es idéntica al caso C++.
#include "Util.h"
#include <sys/time.h>
#include <time.h>
#define IA 16807
#define IQ 127773
#define IR 2836
#define AM (1.0/ALEATORIO_MAX)
#define MASK 123459876
#define ALEATORIO_MAX 2147483647L
static long semilla=170374;
static long aleatorio(void)
{
long k,ans;
semilla^=MASK;
k=semilla/IQ;
semilla=IA*(semilla-k*IQ)-IR*k;
if(semilla<0)
semilla+=ALEATORIO_MAX;
ans=semilla;
semilla^=MASK;
return ans;
}
long aleatorio_entero(long minimo,long maximo)
{
return minimo+(int)((maximo-minimo+1.0)*aleatorio()/(ALEATORIO_MAX+1.0));
}
double getTime()
{
struct timeval tv;
struct timezone tz;
gettimeofday(&tv,&tz);
return tv.tv_sec+tv.tv_usec/1000000.0;
}Implementación en Lenguaje Java
Rutina principal main
Esta clase sólo tiene como finalidad contener el método de arranque "main()" (aunque se pudo agregar a cualquiera de las otras clases, pero considero que así es más ordenado.)
package com.americati.analgeom;
public class Application {
/**
* @param args
*/
public static void main(String[] args) {
// TODO Auto-generated method stub
int z,j,interior=0;
int maxpol=10000,maxpnt=1000;
Polygon[] pol=new Polygon[maxpol];
for(z=0;z<maxpol;z++)
{
pol[z]=new Polygon();
setHepta(pol[z]);
}
// setTriangle(pol[z]);
Point[] p=new Point[maxpnt];
for(z=0;z<maxpnt;z++)
{
p[z]=new Point();
p[z].setX(AGUtil.aleatorio_entero(0,800));
p[z].setY(AGUtil.aleatorio_entero(0,600));
}
double t1,t2;
t1=AGUtil.getTime();
for(j=0;j<maxpnt;j++)
for(z=0;z<maxpol;z++)
{
if(pol[z].isInterior(p[j]))
interior++;
}
t2=AGUtil.getTime();
double ratio=(interior*100.0)/(maxpol*maxpnt);
System.out.println(interior+" interiores de "+
(maxpol*maxpnt)+" = "+ratio+"% (t="+
(t2-t1)+")");
}
static void setTriangle(Polygon t)
{
Point p1=new Point(AGUtil.aleatorio_entero(0,800),
AGUtil.aleatorio_entero(0,600));
Point p2=new Point(AGUtil.aleatorio_entero(0,800),
AGUtil.aleatorio_entero(0,600));
Point p3=new Point(AGUtil.aleatorio_entero(0,800),
AGUtil.aleatorio_entero(0,600));
t.addPoint(p1);
t.addPoint(p2);
t.addPoint(p3);
}
static void setHepta(Polygon t)
{
double radius=20+AGUtil.aleatorio_entero(0,800)/8.0;
double h=AGUtil.aleatorio_entero(0,800);
double k=AGUtil.aleatorio_entero(0,600);
int z;
double PI=3.1415926535;
for(z=0;z<7;z++)
t.addPoint(new Point(h+radius*Math.cos(2*z*PI/7),
k+radius*Math.sin(2*z*PI/7)));
}
}Implementación de Clases
Clase Point
La versión Java de Point es casi idéntica al caso C++:
package com.americati.analgeom;
public class Point {
private double x;
private double y;
public Point()
{
x=y=0;
}
public Point(double a,double b)
{
x=a; y=b;
}
public Point(Point a,Point b)
{
x=b.getX()-a.getX();
y=b.getY()-a.getY();
}
public Point(Point a)
{
x=a.getX();
y=a.getY();
}
public double getX()
{
return x;
}
public double getY()
{
return y;
}
public void setX(double n)
{
x=n;
}
public void setY(double n)
{
y=n;
}
public void print()
{
System.out.println("["+x+","+y+"]");
}
public double distance2(Point p)
{
double dx=x-p.getX();
double dy=y-p.getY();
return dx*dx+dy*dy;
}
public double distance(Point p)
{
return Math.sqrt(distance2(p));
}
public int quadrant(Point p)
{
if(p.getX()>=x)
{
if(p.getY()>y)
return 1;
else
return 4;
}
else
{
if(p.getY()>y)
return 2;
else
return 3;
}
}
public Point ort()
{
double t=y;
y=x;
x=-t;
return this;
}
boolean insideTriangle(Point A,
Point B, Point C)
{
Point PA=new Point(this,A);
Point PB=new Point(this,B);
Point PC=new Point(this,C);
Point PAB=new Point(PA);PAB.add(PB);
Point PBC=new Point(PB);PBC.add(PC);
Point PCA=new Point(PC);PCA.add(PA);
Point PAN=new Point(B,C);PAN.ort();
Point PBN=new Point(C,A);PBN.ort();
Point PCN=new Point(A,B);PCN.ort();
Point AB=new Point(A,B);
Point BC=new Point(B,C);
Point CA=new Point(C,A);
if(AB.prod(PAN)*PB.prod(PAN)<0)
return false;
if(BC.prod(PBN)*PC.prod(PBN)<0)
return false;
if(CA.prod(PCN)*PA.prod(PCN)<0)
return false;
return true;
}
public Point add(Point p)
{
x+=p.getX();
y+=p.getY();
return this;
}
public Point sub(Point p)
{
x-=p.getX();
y-=p.getY();
return this;
}
public double prod(Point p)
{
return x*p.getX()+y*p.getY();
}
}Clase Polygon
El polígono contiene un vector de puntos. Nótese que usamos la sintaxis con "generics" de Java 5 (Vector<Point>.)
package com.americati.analgeom;
import java.util.Vector;
public class Polygon {
private Vector<Point> pt;
private double diameter;
public Polygon()
{
diameter=0;
pt=new Vector<Point>();
}
public int getCount()
{
return pt.size();
}
public Point getPoint(int n)
{
if(n<0 || n>=(int)pt.size())
throw new RuntimeException();
return pt.elementAt(n);
}
public void setPoint(int n,Point p)
{
if(n<0 || n>=pt.size())
return;
pt.setElementAt(p,n);
findDiameter();
}
public void addPoint(Point p)
{
pt.add(p);
findDiameter();
}
public void print()
{
int z;
for(z=0;z<getCount();z++)
getPoint(z).print();
System.out.println(" D="+diameter);
}
private void findDiameter()
{
if(getCount()<3)
return;
int z,j;
double rn;
diameter=0;
for(z=0;z<getCount()-1;z++)
for(j=z+1;j<getCount();j++)
{
rn=getPoint(z).distance2(getPoint(j));
if(rn>diameter)
diameter=rn;
}
diameter=Math.sqrt(diameter);
}
public boolean isInterior(Point p)
{
if(getCount()<3)
throw new RuntimeException();
if(p.distance2(getPoint(0))>diameter*diameter)
return false;
int z,current,Q=0;
for(z=0;z<getCount();z++)
{
current=p.quadrant(getPoint(z));
Q|=(1<<current);
}
if(Q==30) // 4 quadrants
return true;
if(Q!=28 && Q!=26 && Q!=22 && Q!=14 && Q!=10 && Q!=20) // not 3q nor oposed 2q
return false;
for(z=1;z<getCount()-1;z++)
if(p.insideTriangle(getPoint(0),
getPoint(z),getPoint(z+1)))
return true;
return false;
}
}Para no alargar la exposición sólo estamos generando la excepción RuntimeException (para no tener que capturarla.)
Al igual que en el caso C++, el método isInterior() puede comentarse por partes para obtener las variantes 1 y 2.
Nótese que de utilizarse el AWT, sería aconsejable renombrar esta clase pues ya hay una del mismo nombre en dicho paquete (o nos veríamos forzados a usar el nombre calificado.)
Clase AGUtil
La clase AGUtil es también similar a la clase Util de nuestra versión C++. Las constantes del preprocesador son reemplazadas con variables "final". Nótese que hemos evitado el nombre "Util" a fin de no hacer conflicto con java.Util; de lo contrario tendríamos que utilizar el nombre calificado haciendo el código más engorroso.
package com.americati.analgeom;
public class AGUtil {
final static long ALEATORIO_MAX=2147483647L;
final static long IA=16807;
final static long IQ=127773;
final static long IR=2836;
final static long AM=(long)(1.0/ALEATORIO_MAX);
final static long MASK=123459876;
private static long semilla=170374;
private static long aleatorio()
{
long k,ans;
semilla^=MASK;
k=semilla/IQ;
semilla=IA*(semilla-k*IQ)-IR*k;
if(semilla<0)
semilla+=ALEATORIO_MAX;
ans=semilla;
semilla^=MASK;
return ans;
}
public static long aleatorio_entero(long minimo,long maximo)
{
return minimo+(int)((maximo-minimo+1.0)*aleatorio()/(ALEATORIO_MAX+1.0));
}
public static double getTime()
{
return System.currentTimeMillis()/1000.0;
}
}Resultados
Las implementaciones arrojaron 294468 casos en los que los puntos quedaba al interior de los polígonos (de un total de 10000000 tests.) Consideramos de interés para el lector presentar los tiempos obtenidos en estas ejecuciones, sin embargo, debemos advertir que bajo ninguna circunstancia dichos valores representan resultados definitivos o un "benchmark" formal. Los resultados presentados corresponden a los mínimos tiempos obtenidos en un total de cinco ejecuciones por cada caso, y sólo marginalmente hemos intentado evitar que otras aplicaciones distorsionen el consumo de recursos del sistema operativo. En [SPEC] se puede hallar información más completa para hacer un "benchmark" más objetivo; los valores presentados pueden considerarse hasta cierto punto "anecdóticos".
En todos los casos se utilizó Linux Ubuntu Drake con kernel 2.6.15. El compilador de C y C++ fue GCC 4.0.3. Para el caso de Java, los resultados iniciales se obtuvieron con el JDK de Sun 1.5.0 R09; sin embargo, también se proporcionan resultados obtenidos con otras máquinas virtuales.
Los tiempos (en segundos) se resumen en la tabla a continuación.
| Lenguaje | Variante 1 | Variante 2 | Variante 3 |
|---|---|---|---|
C |
6.8 |
4.6 |
3.1 |
C++ |
10.5 |
6.7 |
3.8 |
Java |
90.9 |
11.7 |
7.8 |
Observaciones
Los tiempos obtenidos en relación a las implementaciones de C y C++ son bastante razonables en la medida que el código C++ contiene más validaciones soportadas por el mismo lenguaje (encapsulación) así como controles agregados como parte del código (por ejemplo, índices fuera de rango.) Claramente, el código C++ resulta más comprensible, ordenado y extensible, y de requerirse más velocidad se podría sacrificar gradualmente la orientación a objetos (quizá, hasta convertirlo en la versión en C.)
Dos de los tiempos de Java son razonables si se toma en cuenta su ejecución en "bytecodes", mayor control de excepciones en tiempo de ejecución, y el "garbage collector"; sin embargo, el "90.56" resulta dificil de aceptar.
Analizando la ejecución (esto se puede hacer con herramientas de profiling) hallamos que el método Point.insideTriangle() es el causante principal del retardo. El siguiente código es una reescritura del mismo, la cual ciertamente tiene muy poco de "orientada a objetos":
boolean insideTriangle(Point A, Point B, Point C)
{
double ax=A.getX(),ay=A.getY();
double bx=B.getX(),by=B.getY();
double cx=C.getX(),cy=C.getY();
double PAx=ax-x, PAy=ay-y;
double PBx=bx-x, PBy=by-y;
double PCx=cx-x, PCy=cy-y;
double PANx=-cy+by, PANy=cx-bx;
double PBNx=-ay+cy, PBNy=ax-cx;
double PCNx=-by+ay, PCNy=bx-ax;
double ABx=bx-ax, ABy=by-ay;
double BCx=cx-bx, BCy=cy-by;
double CAx=ax-cx, CAy=ay-cy;
if((ABx*PANx+ABy*PANy)*(PBx*PANx+PBy*PANy)<0)
return false;
if((BCx*PBNx+BCy*PBNy)*(PCx*PBNx+PCy*PBNy)<0)
return false;
if((CAx*PCNx+CAy*PCNy)*(PAx*PCNx+PAy*PCNy)<0)
return false;
return true;
}Como veremos, esta modificación redujo el tiempo de la "variante 1" considerablemente. Esto ciertamente indica que la implementación de las regiones "críticas" deben analizarse con cuidado, y sin temor a seguir enfoques poco ortodoxos.
Nos pareció interesante contrastar este enfoque "no POO" (NPOO) con el original (al que llamaremos POO), utilizando las tres variantes y con diferentes máquinas virtuales (con bytecode generado con sus respectivos compiladores.) Los resultados se muestran a continuación:
Resultados con otras máquinas virtuales
| JVM | Var 1 POO | Var 2 POO | Var 3 POO | Var 1 NPOO | Var 2 NPOO | VAR 3 NPOO |
|---|---|---|---|---|---|---|
Sun 1.5.0 R9 |
90.9 |
11.7 |
7.8 |
21.0 |
10.2 |
6.7 |
Sun 6.0 R1 |
81.1 |
12.0 |
10.6 |
15.1 |
9.3 |
6.4 |
Ibm 5.0-4.0 |
109 |
17.0 |
12.0 |
16.1 |
14.9 |
9.4 |
Jrockit 1.5.0_06 |
71.2 |
22.2 |
16.2 |
18.5 |
15.8 |
12.4 |
Descripción del Algoritmo de Decisión
Primera Variante
Considérese la figura A.1 en la cual tenemos el polígono (P0,P1,P2,P3,P4); si tomamos "P0" como origen, podemos subdividir el polígono en un conjunto de triángulos (recuérdese que sólo consideramos polígonos convexos.) Si podemos determinar si un punto "X" está al interior de cualquiera de estos triángulos, entonces el punto está al interior del polígono y el problema queda resuelto.
Punto al interior del un triángulo
Considérese la figura A.2:
Como se puede apreciar, la recta que contiene al lado del triángulo divide al plano en dos regiones (semiplanos); es claro que en una de ellas se encuentra el tercer vértice. Si ocurre (tal como en la figura) que el punto de análisis está en el "semiplano opuesto", entonces podemos afirmar que dicho punto está al exterior del triángulo.
Por el contrario, si el punto está en el mismo semiplano que el tercer vértice, entonces no podemos afirmar nada al respecto; sin embargo, si hacemos este análisis para los tres lados del triángulo, es claro que si nuestro punto es exterior, entonces deberá necesariamente caer uno de los "semiplanos opuestos"; de lo contrario, es interior.
Con esto hemos reducido el análisis a determinar si el punto está en el semiplano del tercer vértice o en el "semiplano opuesto"; dicho de otro modo, queremos saber si dos puntos están del mismo lado de una recta, o en regiones opuestas.
Puntos con respecto a una recta
La figura A.3 ilustra el caso de dos puntos (A y X) que se encuentran en regiones opuestas a la recta determinada por (B,C). En la figura se ha trazado también un vector ortogonal a la recta (B,C) al que podemos llamar "N".
Ahora podemos trazar dos vectores desde A y X hacia cualquier punto de la recta (por ejemplo, hacia "C"). Estos dos vectores ("v1" y "v2") al proyectarse en "N" van a resultar en sentidos opuestos; sin embargo, si "X" estuviera del mismo lado que "A", entonces los vectores "v1" y "v2" al proyectarse resultarían con el mismo sentido.
Una forma fácil de determinar el sentido de la proyección de un vector sobre otro, consiste en analizar el signo de su producto escalar:
Sean los productos escalares P1=v1.N y P2=v2.n, entonces si tienen signos opuestos tendremos P1.P2<0, en caso contrario, P1.P2>0.
Esto proporciona un criterio de fácil cálculo para determinar la posición relativa, con lo que el problema inicial queda resuelto.
Segunda Variante
Con esta variante estamos intentando aligerar el las operaciones a fin de evitar realizar el análisis antes descrito. Como se verá, no siempre podemos evitarlo, pero los resultados de las pruebas demuestran que vale la pena intentarlo.
Si a partir del punto de análisis trazamos ejes de coordenadas perpendiculares, entonces los puntos del polígono pueden quedar repartidos en diversos cuadrantes. Por ejemplo, si todos los puntos quedan en un único cuadrante (cualquiera de los cuatro) es evidente que el punto es exterior al polígono (fig A.5-a.) De igual modo, si todos los puntos del polígono se ubican en dos cuadrantes 'adyacentes', entonces el punto también es exterior (fig A.5-b.) De igual modo, si los puntos del polígono se ubican en los cuatro cuadrantes, entonces necesariamente el punto será 'interior' (fig A.5-c.)
Esto nos permite determinar rápidamente la ubicación del punto para ciertos casos. Lamentablemente, en ciertas situaciones tendremos que recurrir al análisis de la primera variante. Por ejemplo, en la figura A.6 se muestran dos casos del polígono en tres cuadrantes:
Y la figura A.7 muestra dos casos del polígono en dos cuadrantes no adyacentes:
Son varios los casos a considerar (en total 16-1 pues se descarta el caso en el que no hay puntos del polígono), sin embargo, es extremadamente sencillo y eficiente el analizar los cuadrantes de los puntos (basta ver si sus coordenadas son mayores o menores que las del punto de análisis.)
Tercera Variante
Considérese la distancia máxima entre dos puntos cualesquiera del polígono. A esta distancia le llamaremos "diámetro" del polígono.
Ahora consideremos la distancia que hay entre nuestro punto de análisis y un punto arbitrario del polígono. Si esta distancia es mayor que el diámetro, es claro que el punto debe estar fuera del polígono; lamentablemente, si la distancia es menor entonces el punto puede estar dentro o fuera y debemos proceder a la segunda variante:
A fin de hacer eficiente este análisis, será necesario que cada polígono conserve el valor de su diámetro (lo puede calcular cada vez que se le agrega un nuevo punto) a fin de no tener que calcularlo repetidamente.
Referencias
-
[MN] Press W. H. et. al. "Numerical Recipes in C: The Art of Scientific Computing" Cambridge University Press. 2 edition. 1992.
-
[SPEC] SPEC CPU2006 Run and Reporting Rules http://www.spec.org/cpu2006/Docs/runrules.html