1065 - Raising the Roof

You have agreed to help a roofer estimate the amount of material required to protect the surface of a roof from sun damage by computing the surface area of the roof that can be seen when viewed from directly above.

The geometry of a roof is defined by the vertices of the triangles that enclose the plane surfaces of the roof. The x and y coordinates of a vertex are determined from an aerial view of the roof, and the zcoordinate is the height of the vertex above the ground, assuming that the ground is level (which we do). For example, a line segment with vertices (10,10,10) and (10,20,10) is parallel to the ground, and is 10 units long. A line segment with vertices (10,10,10) and (10,18,16) slopes from 10 units to 16 units above the ground; in the aerial view of the roof it is 8 units long, but the actual length is 10 units.

In an aerial view of a roof, one region of the roof can sometimes obscure a lower region. Only visible portions of the roof should be included in the total. For example, suppose you view a two-story building with a second floor that is smaller than the first floor. Obviously the area of the second floor's roof must be included in the total, but the visible surface area of the first floor's roof will not include the area directly underneath the second floor's roof.

Input 

The input consists of several test cases corresponding to roof descriptions. A test case begins with a line containing integers V (1V300) , which is the number of vertices, and T (1T1000) , which is the number of triangles. Each of the next V lines contains the coordinates x , y , z of a vertex. The test case concludes with T triangle descriptions, each consisting of the three vertex indices describing the vertices of a roof triangle. Vertices are numbered sequentially in the order they appear in the input, starting with 1.

All coordinates are positive integers no greater than 100. No roof triangles are degenerate, and all pairs of roof triangles have disjoint interiors -- that is, they may touch but they will not overlap or intersect.

The last test case is followed by a line containing two zeroes.

Output 

For each test case, print a line containing the test case number (beginning with 1) followed by the total visible surface area of the roof, accurate to two fractional digits. Use the sample output format and print a blank line after each case.

Sample Input 

6 2 
10 10 10 
10 20 10 
20 10 10 
10 10 20 
10 20 20
20 20 20
1 2 3 
4 5 6 
3 1 
10 10 
10 18 16 
20 10 10 
1 2 3 
0 0

Sample Output 

Case 1: 75.00 

Case 2: 50.00

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
int N,T;
long double pos[305][3];
int Graph[305][305];
int Tri[1005][3];
long double Area[1005];
long double allX[1000005];
long double pool[2][1005];
long double line[1005][2][2];
int Order[2005][2];

int Heap[1005],Size;
int pointer[1005];
long double tmp[2];
long double *cnt,*last;
int s,t;
long double Ans;
long double m1,m2,m3,m4;
int Test=0;
#define Eps 1e-13
#define Equal(a,b) (fabsl((a)-(b))<Eps)
#define Max(a,b) ((a)>(b)?(a):(b))
#define Min(a,b) ((a)<(b)?(a):(b))
#define sqr(a) ((a)*(a))
#define Luck -112901989

long double Cross(long double *pa,long double *pb,long double *pc)
{
	return (pb[0]-pa[0])*(pc[1]-pa[1]) - (pb[1]-pa[1])*(pc[0]-pa[0]);
}

int cmp(const void *a,const void *b)
{
	return (*(long double *)a > *(long double *)b ?1 :-1);
}

int cmp2(const void *a,const void *b)
{
	if(!Equal(line[((int *)a)[0]][((int *)a)[1]][0], line[((int *)b)[0]][((int *)b)[1]][0]))
		return (line[((int *)a)[0]][((int *)a)[1]][0] - line[((int *)b)[0]][((int *)b)[1]][0]?1:-1);
	return (((int *)a)[1] < ((int *)b)[1] ? 1:-1);
}

int Check()
{
	if(!Equal(tmp[0],Luck))
	{
		if(m1>tmp[0])
		{
			m1=tmp[0];
			m3=tmp[1];
		}
		if(m2<tmp[0])
		{
			m2=tmp[0];
			m4=tmp[1];
		}
	}
	return 0;
}

int Inner(long double *pa,long double *pb, long double t)
{
	long double p1,p2;
	p1=pb[0]-t;
	p2=t-pa[0];
	if(p1<-Eps && p2> Eps || p1>Eps && p2<-Eps)
		tmp[0]=Luck;
	else
	{
		tmp[0]=(p1*pa[1] + p2*pb[1])/(p1+p2);
		tmp[1]=(p1*pa[2] + p2*pb[2])/(p1+p2);
	}
	return 0;
}

int Swap(int x)
{
	int i=Heap[x];
	Heap[x]=Heap[(x-1)/2];
	Heap[(x-1)/2]=i;

	pointer[Heap[x]]=x;
	pointer[Heap[(x-1)/2]]=(x-1)/2;
	return 0;
}

int Greater(int a,int b)
{
	long double p1,p2,std;
	std=(Max(line[a][0][0],line[b][0][0]) + Min(line[a][1][0],line[b][1][0]))/2;

	Inner(line[a][0],line[a][1],std);
	p1=tmp[0];
	Inner(line[b][0],line[b][1],std);
	p2=tmp[0];
	if(Equal(p1,Luck) || Equal(p2,Luck))
		while(1);
	return (p1-p2);
}

int modifyUp(int x)
{
	while(x&&Greater(Heap[x],Heap[(x-1)/2]))
	{
		Swap(x);
		x=(x-1)/2;
	}
	return 0;
}

int modifyDown(int x)
{
	int i;
	while(1)
	{
		i=x;
		if(2*x+1<Size && Greater(Heap[2*x+1],Heap[i]))
			i=2*x+1;
		if(2*x+2<Size && Greater(Heap[2*x+2],Heap[i]))
			i=2*x+2;
		if(i==x)
			break;
		Swap(i);
		x=i;
	}
	return 0;
}

int Insert(int x)
{
	Heap[pointer[x]=Size++]=x;
	modifyUp(Size-1);
	return 0;
}

int Delete(int x)
{
	Heap[x]=Heap[--Size];
	pointer[Heap[x]]=x;
	if(x<Size)
	{
		modifyUp(x);
		modifyDown(x);
	}
	return 0;
}

int main()
{
	int i,j,k,l;
	long double *tmp2;
	long double u[3],v[3],vec[3];
	long double a;
	long double t1,t2,b1,b2,c1,c2;
	while(scanf("%d%d",&N,&T)&&(N||T))
	{
		for(i=0;i<N;i++)
			for(j=0;j<3;j++)
			{
				double p;
				scanf("%lf",&p);
				pos[i][j]=p;
			}
		for(i=0;i<N;i++)
			for(j=0;j<N;j++)
				Graph[i][j]=0;
		for(i=0;i<T;i++)
		{
			scanf("%d%d%d",&j,&k,&l);
			j--;
			k--;
			l--;
			if(!Equal(Cross(pos[j],pos[k],pos[l]),0))
			{
				Tri[i][0]=j;
				Tri[i][1]=k;
				Tri[i][2]=l;
				Graph[j][k]=Graph[k][j]=Graph[j][l]=Graph[l][j]=Graph[k][l]=Graph[l][k]=1;
			}
			else
			{
				T--;
				i--;
			}
		}
		s=0;
		for(i=0;i<N;i++)
			for(j=i+1;j<N;j++) if(Graph[i][j])
				for(k=i+1;k<N;k++)
					for(l=k+1;l<N;l++)
						if(Graph[k][l]&&Cross(pos[i],pos[j],pos[k])*Cross(pos[i],pos[j],pos[l])<0
							&&Cross(pos[k],pos[l],pos[i])*Cross(pos[k],pos[l],pos[j])<0)
						{
							t1=pos[l][0]-pos[k][0];
							t2=pos[l][1]-pos[k][1];
							b1=pos[i][0]-pos[j][0];
							b2=pos[j][0]-pos[k][0];
							c1=pos[i][1]-pos[j][1];
							c2=pos[j][1]-pos[k][1];
							a=(t2*b2-t1*c2)/(t1*c1-t2*b1);
							allX[s++]=a*(pos[i][0]-pos[j][0]) + pos[j][0];
						}
		for(i=0;i<N;i++)
		{
			allX[s++]=pos[i][0]-1.2*Eps;
			allX[s++]=pos[i][0]+1.2*Eps;
		}
		qsort(allX,s,sizeof(long double),cmp);
		memset(pool,0,sizeof(pool));
		memset(Area,0,sizeof(Area));
		last=pool[0];
		cnt=pool[1];
		for(i=j=1;i<s;i++)
			if(!Equal(allX[j-1],allX[i]))
				allX[j++]=allX[i];
		s=j;
		for(i=0;i<s;i++)
		{
			a=allX[i];
			k=0;
			l=0;
			for(j=0;j<T;j++)
			{
				m1=1E20;
				m2=-1E20;
				Inner(pos[Tri[j][0]],pos[Tri[j][1]],a);

				Check();
				Inner(pos[Tri[j][1]],pos[Tri[j][2]],a);

				Check();
				Inner(pos[Tri[j][0]],pos[Tri[j][2]],a);

				Check();
				if(m1<m2-Eps)
				{
					line[j][0][0]=m1;
					line[j][0][1]=m3;
					line[j][1][0]=m2;
					line[j][1][1]=m4;
					Order[l][0]=Order[l+1][0]=j;
					Order[l][1]=0;
					Order[l+1][1]=1;
					l+=2;
				}
			}
			qsort(Order,l,sizeof(int)*2,cmp2);
			t=l;
			memset(cnt,0,sizeof(pool[0]));
			Size=0;
			for(j=0;j<t;j++)
			{
				if(Size)
					cnt[Heap[0]]+=line[Order[j][0]][Order[j][1]][0]-line[Order[j-1][0]][Order[j-1][1]][0];
				k=Order[j][0];
				l=Order[j][1];
				if(l==0)
					Insert(k);
				else
					Delete(pointer[k]);
			}
			if(i)
				for(j=0;j<T;j++)
					Area[j]+=(cnt[j]+ last[j])/2*(allX[i]-allX[i-1]);
			tmp2=cnt;
			cnt=last;
			last=tmp2;
		}
		Ans=0;
		for(i=0;i<T;i++)
		{
			for(j=0;j<3;j++)
			{
				u[j]=pos[Tri[i][1]][j]-pos[Tri[i][0]][j];
				v[j]=pos[Tri[i][2]][j]-pos[Tri[i][0]][j];
			}
			vec[2]=1;
			vec[0]=(u[1]*v[2] - v[1]*u[2])/(v[1]*u[0] - u[1]*v[0]);
			vec[1]=(u[0]*v[2] - v[0]*u[2])/(v[0]*u[1] - u[0]*v[1]);
			Ans+=Area[i] * sqrtl(sqr(vec[0]) + sqr(vec[1]) + 1);
		}
		printf("Case %d: %.2llf\n\n",++Test,Ans);
	}
	return 0;
}