/*
this program for calculating area of closed polygon
*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include<ran3.c>
double roundintagral(double *x, double *y, int line);
double montecalro(double *x, double *y, int line);
int pnpoly(int npol, float *xp, float *yp, float x, float y);
void swap(float *a, float *b);
int main(int argc, char **argv)
{
FILE *ofp;
char str[20];
ofp=fopen(argv[1], "r");
double *x, *y;
int line;
int i;
line=0;
while(fgets(str, sizeof(str), ofp) && ftell(ofp) != EOF)
line++;
rewind(ofp);
x=calloc(line, sizeof(double));
y=calloc(line, sizeof(double));
for(i=0; i<line; i++)
fscanf(ofp, "%lf %lf", &x[i], &y[i]);
printf("%s\t%lf\t%lf\n", argv[1],roundintagral(x, y, line), montecalro(x, y, line));
free(x);
free(y);
fclose(ofp);
return 0;
}
double roundintagral(double *x, double *y, int line)
{
double centerX, centerY;
double area;
int i;
area=0;
centerX=0; centerY=0;
for(i=0; i<line; i++)
{
centerX+=x[i];
centerY+=y[i];
}
centerX=centerX/line;
centerY=centerY/line;
for(i=0; i<line-1; i++)
area += (0.5 * ( (x[i]-centerX)*(y[i+1]-centerY) - (x[i+1]-centerX)*(y[i]-centerY) ) );
return fabs(area);
}
double montecalro(double *x, double *y, int line)
{
double area;
float *xsort, *ysort;
float *xunsort, *yunsort;
double fraction;
float xrand, yrand;
int i,j;
long seed = time(0);
fraction=0;
xsort=calloc(line, sizeof(float));
ysort=calloc(line, sizeof(float));
xunsort=calloc(line, sizeof(float));
yunsort=calloc(line, sizeof(float));
for(i=0; i<line; i++)
{
xsort[i]=x[i];
ysort[i]=y[i];
xunsort[i]=x[i];
yunsort[i]=y[i];
}
for(i=0; i<line; i++)
for(j=i+1; j<line; j++)
{
if(xsort[i] > xsort[j]) swap(&xsort[i], &xsort[j]);
if(ysort[i] > ysort[j]) swap(&ysort[i], &ysort[j]);
}
area=(xsort[line-1]-xsort[0])*(ysort[line-1]-ysort[0]);
for(i=0; i<10000; i++)
{
xrand=xsort[0]+ran3(&seed)*(xsort[line-1]-xsort[0]);
yrand=ysort[0]+ran3(&seed)*(ysort[line-1]-ysort[0]);
if(pnpoly(line, xunsort, yunsort, xrand, yrand))
fraction++;
}
area=area*(fraction/10000);
free(xsort);
free(ysort);
free(xunsort);
free(yunsort);
return area;
}
int pnpoly(int npol, float *xp, float *yp, float x, float y)
{
int i, j, c = 0;
for (i = 0, j = npol-1; i < npol; j = i++)
{
if ((((yp[i] <= y) && (y < yp[j])) || ((yp[j] <= y) && (y < yp[i]))) && (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
c = !c;
}
return c;
}
void swap(float *a, float *b)
{
float temp;
temp = *a;
*a = *b;
*b = temp;
}
여기선 두가지 방법으로 넓이를 구합니다.
구면 적분을 한것과 몬테카를로방법을 사용한것입니다.
사용방법은 컴파일을 하신다음에 명령행(리눅스) 기준으로 쓰면
$ ./area 읽어들일데이터파일
이렇게 간단하게 쓰시면 되구요. 데이터파일의 모양은
0.0 0.0
0.0 1.0
1.0 1.0
1.0 0.0
0.0 0.0
이런 식이 되면 좌표계에서 사각형의 넓이를 구하게 되는 것이죠.
값은 1이 나오게 됩니다.
ps. ran3.c라는 파일은 numerical recipe에 있는 것인데요. 그냥 지우고 rand함수를 써도 무관합니다. 그래도 ran3가 더 정확하겠죠?
