Photolog

Through the Looking-Glass
2010-10-12: Through the Looking-Glass
My radio speaks is binary!
2010-10-10: My radio speaks is binary!
Gigaminx: (present for my birthday)
2010-09-16: Gigaminx: (present for my birthday)
Trini on bike
2010-09-05: Trini on bike
Valporquero
2010-08-28: Valporquero
My new bike!
2010-08-22: My new bike!
Mario and Ana's wedding
2010-08-13: Mario and Ana's wedding
Canyoning in Guara
2010-08-07: Canyoning in Guara
Trini and Mari in Marbella
2010-08-05: Trini and Mari in Marbella
Trini and Chelo in Tabarca
2010-08-03: Trini and Chelo in Tabarca
Valid XHTML 1.1
Log in
Back to list of problems

Polygons

137.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define EPS 1E-10
#define ZERO(a) (fabs(a) < EPS)
#define SIGN(a) (ZERO(a) ? 0 : (a>0) ? 1 : -1)

typedef struct {
	double x;
	double y;
} punto;

int
minang(int o, int nl, punto *li) {
	int i,m=!o;
	for (i=0; i<nl; i++) {
		if (i==o) continue;
		if ((li[i].x-li[o].x)*(li[m].y-li[o].y)
			> (li[i].y-li[o].y)*(li[m].x-li[o].x)) {
			m=i;
		}
	}
	return m;
}

int
convex_hull(int N, punto *li) {
	int i,m=0;
	punto p;
	if (!N) return 0;
	for (i=1; i<N; i++) {
		double c = li[i].y - li[m].y;
		if (c<0) m=i;
		if (!c && (li[i].x < li[m].x)) m=i;
	}
	p=li[0]; li[0]=li[m]; li[m]=p;
	i=0; m=1;
	while (1) {
		i = minang(m-1, N, li);
		if (i==0) return m;
		p=li[m]; li[m]=li[i]; li[i]=p;
		m++;
	}
}

typedef struct {
	int np;
	punto *p;
} poly;

double
cross(punto p1, punto p2) {
	return p1.x*p2.y - p1.y*p2.x;
}

double
modulo(punto p) {
	return sqrt(p.x*p.x + p.y*p.y);
}

punto
psub(punto p1, punto p2) {
	punto p;
	p.x = p2.x - p1.x;
	p.y = p2.y - p1.y;
	return p;
}

#define dist(p1,p2) modulo(psub(p1,p2))

int
inside(punto pu, poly po) {
	int i;
	int sign=0;

	for (i=0; i<po.np; i++) {
		double d = cross(psub(po.p[i], po.p[(i+1)%po.np]), psub(po.p[i], pu));
		if (ZERO(d)) {
			continue;
		}
		if (!sign) {
			sign = (int)(d / fabs(d));
		} else if (sign != (int)(d / fabs(d))) {
			return 0;
		}
	}
	return 1;
}

double
area(poly po) {
	int i;
	double ar=0.0;

	for (i=1; i<po.np-1; i++) {
		double a,b,c,s;
		a = dist(po.p[0], po.p[i]);
		b = dist(po.p[i], po.p[i+1]);
		c = dist(po.p[0], po.p[i+1]);
		s = (a+b+c)/2.0;
		ar += sqrt(s*(s-a)*(s-b)*(s-c));
	}
	return ar;
}

double
coordy(punto p1, punto p2, double x) {
	return p1.y + (x-p1.x)*(p2.y-p1.y) / (p2.x-p1.x);
}

punto
corte(punto p11, punto p12, punto p21, punto p22) {
	punto p;
	double A,B;

	if (p11.x == p12.x) {
		p.x = p11.x;
		p.y = coordy(p21, p22, p.x);
		return p;
	}
	if (p21.x == p22.x) {
		p.x = p21.x;
		p.y = coordy(p11, p12, p.x);
		return p;
	}
	A = (p12.y - p11.y) / (p12.x - p11.x);
	B = (p22.y - p21.y) / (p22.x - p21.x);
	p.x = (p21.y - p11.y + A*p11.x - B*p21.x) / (A-B);
	p.y = p11.y + A*(p.x-p11.x);
	return p;
}

int
se_cortan(punto p11, punto p12, punto p21, punto p22) {
	if (SIGN(cross(psub(p11,p12), psub(p11,p21))) * SIGN(cross(psub(p11,p12), psub(p11,p22))) == -1 &&
			SIGN(cross(psub(p21,p22), psub(p21,p11))) * SIGN(cross(psub(p21,p22), psub(p21,p12))) == -1) {
		return 1;
	}
	return 0;
}

poly
calc(poly po1, poly po2, int s) {
	int i,j,k;
	poly inter;

	inter.np = 0;
	inter.p = NULL;

	for (i=0; i<po1.np; i++) {
		if (inside(po1.p[i], po2)) {
			int done=0;
#if DEBUG
			printf("point (%d,%d) is inside poly 2\n", (int)po1.p[i].x, (int)po1.p[i].y);
#endif
			for (k=0; k<inter.np; k++) {
				if (ZERO(po1.p[i].x-inter.p[k].x) && ZERO(po1.p[i].y-inter.p[k].y)) {
					done=1;
					break;
				}
			}
			if (!done) {
				inter.p = realloc(inter.p, (inter.np+1)*sizeof(punto));
				if (inter.np>500) abort();
				if (!inter.p) abort();
				inter.p[inter.np] = po1.p[i];
				inter.np++;
			}
		}
	}
	for (i=0; i<po2.np; i++) {
		if (inside(po2.p[i], po1)) {
			int done=0;
#if DEBUG
			printf("point (%d,%d) is inside poly 1\n", (int)po2.p[i].x, (int)po2.p[i].y);
#endif
			for (k=0; k<inter.np; k++) {
				if (ZERO(po2.p[i].x-inter.p[k].x) && ZERO(po2.p[i].y-inter.p[k].y)) {
					done=1;
					break;
				}
			}
			if (!done) {
				inter.p = realloc(inter.p, (inter.np+1)*sizeof(punto));
				if (inter.np>500) abort();
				if (!inter.p) abort();
				inter.p[inter.np] = po2.p[i];
				inter.np++;
			}
		}
	}
	for (i=0; i<po1.np; i++) {
		int i1 = (i+1)%po1.np;
		for (j=0; j<po2.np; j++) {
			int j1 = (j+1)%po2.np;
			if (se_cortan(po1.p[i], po1.p[i1], po2.p[j], po2.p[j1])) {
				int done=0;
				punto p = corte(po1.p[i], po1.p[i1], po2.p[j], po2.p[j1]);
#if DEBUG
				printf("segmentos (%f,%f)-(%f,%f) y (%f,%f)-(%f,%f) se cortan en (%f,%f)\n",
						po1.p[i].x, po1.p[i].y, po1.p[i1].x, po1.p[i1].y,
						po2.p[j].x, po2.p[j].y, po2.p[j1].x, po2.p[j1].y,
						p.x, p.y);
#endif
				for (k=0; k<inter.np; k++) {
					if (ZERO(p.x-inter.p[k].x) && ZERO(p.y-inter.p[k].y)) {
						done=1;
						break;
					}
				}
				if (!done) {
					inter.p = realloc(inter.p, (inter.np+1)*sizeof(punto));
					if (inter.np>500) abort();
					if (!inter.p) abort();
					inter.p[inter.np] = p;
					inter.np++;
				}
			}
		}
	}
#if DEBUG
	printf("=> intersection has %d points:", inter.np);
	for (i=0; i<inter.np; i++) {
		printf(" (%f,%f)", inter.p[i].x, inter.p[i].y);
	}
	printf("\n");
#endif
	convex_hull(inter.np, &inter.p[0]);
#if DEBUG
	printf("=> Convex Hull of intersection:");
	for (i=0; i<inter.np; i++) {
		printf(" (%f,%f)", inter.p[i].x, inter.p[i].y);
	}
	printf("\n");
	printf("area intersection = %f\n", area(inter));
#endif
	return inter;
}

int
main(void) {
	int i;
	int n;
	poly po1, po2, inter;
	double a1, a2, a3;
	int s=0;

	while (1) {
		scanf("%d", &n);
		if (n==0) {
			break;
		}
		po1.np = n;
		po1.p = malloc(n*sizeof(punto));
		if (po1.np>500) abort();
		if (!po1.p) abort();
		for (i=0; i<n; i++) {
			scanf("%lf %lf", &po1.p[i].x, &po1.p[i].y);
		}
		scanf("%d", &n);
		po2.np = n;
		po2.p = malloc(n*sizeof(punto));
		if (po2.np>500) abort();
		if (!po2.p) abort();
		for (i=0; i<n; i++) {
			scanf("%lf %lf", &po2.p[i].x, &po2.p[i].y);
		}
		inter = calc(po1, po2, s);
/*		for (i=0; i<10000000; i++) {n=i;} */
		a1 = area(po1);
		a2 = area(po2);
		a3 = area(inter);
#if DEBUG
		printf("area po1 = %8.2f\n", a1);
		printf("area po2 = %8.2f\n", a2);
		printf("area int = %8.2f\n", a3);
		printf("=>       = %8.2f\n", a1+a2-2.0*a3);
#endif
		printf("%8.2f", a1+a2-2.0*a3);
		free(po1.p);
		free(po2.p);
		free(inter.p);
		s++;
	}
	printf("\n");
	return 0;
}