USACO Training Section 3.4 Closed Fences

[url="http://ace.delos.com/usacoanal2?a=pgBaqh1g9as&S=fence4"]英文原题[/url]

[size=medium]题意[/size]
[size=small]一个闭合的栅栏是平面上的一些不相交的首尾相连的线段形成的多边形，有N个角(顶点) (3 < N < 200)。 顶点不重合，它以逆时针方式以数组{xi, yi}给出(i=1,2,...,N)。每一对相邻的顶点都是一条栅栏。因此共有N条栅栏 (定义xN+1=x1, yN+1=y1)。


编写一个程序实现下面的任务:检查输入的顶点列表{xi,yi}, i=1,2,...,N, 判断它是否为一个合法的闭合栅栏。 找出所有可以被站在点(x,y)处的人所能看到的栅栏(忽略人的高度)，因为有的栅栏会被另外的栅栏所挡住。

只有当存在从(x,y)发射的一条射线第一个穿过栅栏i时，栅栏i是可以被看见的。如果栅栏是平行于目光的，它并不认为是可以看见的。

输出是一个可见栅栏的列表。栅栏用两端的顶点表示，顶点的输出顺序以输入文件中的顺序为准。把栅栏按照最后一个点在输入文件中的顺序排序。如果两条栅栏的最后一个点是一样的，就以它们第一个点的顺序排序。 [/size]

[size=medium]分析和实现[/size]
[size=small]
计算几何题的特点就是边界条件复杂。这道题看上去简单，做起来一旦考虑各种边界条件，非常累人。在网上看到过非常繁琐的题解。标程所使用的二分法看上去简单，但是在极端情况下不正确。

这里提出和实现一个最简单的方式，几乎不用考虑边界条件。

考虑从观察点到边的顶点的射线，这些射线与边一定有交点（最少有两条边：包含该顶点的两条边）。如果边和射线的交点不在端点上，我们也认为在端点上：不妨认为交点也是多边形的一个顶点。这样，只要考虑射线与边交与端点的情况（如下图中的边４和边７）。从而，如图所示：
[img]http://dl.iteye.com/upload/picture/pic/53465/0559c54f-75db-3ba8-88d5-0cb78ba9a1a7.jpg[/img]

射线将平面划分为两部分，从直观上看，为顺时针方向和逆时针方向。从计算角度看，是根据其边向量与射线向量的叉积(cross product)的符号确定的。所以，两半平面的边向量可分别考虑。

对同一半平面上的边向量，有两个参数：
1.与射线的交点的位置，我们在计算时可用射线的轨迹方程x=x_o+(x_p-x_o)*t, y=y_o+(y_p-y_o)*t，其中(x_o,y_o)为观察点，(x_p,y_p)为边顶点。则，交点t>0越大，离观察点越远，必然被更近的向量覆盖。如图中边7,8,4的逆时针半平面部分不可见，因为0,1,2的交点距离观察点更近。
2.相交的角度，这决定了相同的交点下那条边可以被看到。角度可以用向量的点积来表示，也很容易得到。如图中0,1,2中0可见，4,5,6中4可见。（图中45同端点的那个9是个6）。

需要注意两点：
1.将相交的点当作边的顶点来看待不是真正的加入一个顶点，而是在计算时同时更新左右两半平面的值。
2.对一般的栅栏而言，多边形的顶点可以是相交的。从而按照题目要求的输出条件输出是很麻烦的事。标程中有错，无法处理顶点重合的情况。
3.判断多边形边是否相交很平凡，而测试数据中没有这种情况，直接略过了。

[/size]

/*
ID: blackco3
TASK: fence4
LANG: C++
*/
#include <iostream>
#include <math.h>
#include <memory.h>
using namespace std;
const int _max_points_(300) ;
const double _delta = 0.0000001;

inline bool fequal(double x,double y){
	return (x==y) || (x<y &&y-x<_delta) || (x>y && x-y<_delta) ;
}
struct t_point {
	double x,y;
	t_point(double a=0,double b=0):x(a),y(b){};
} observer, poly[_max_points_] ;
bool operator < (t_point & p1 , t_point & p2 ){
	return p1.x==p2.x ? p1.y<p2.y : p1.x<p2.x ;
}
bool operator == (t_point & p1 , t_point & p2 ){
	return p1.x==p2.x && p1.y==p2.y ;
}

/*------------------------------------------------------------------------------------------- 
 *	radial start from p1(time 0), pass p2(time 1), segment <q1(t=0),q2(t=1)> 
 *	p1 is never inside segment <q1,q2>, never equal with p2
 *	if res.type is e_cover, res.t1/t2 means the radio time of segment endpoint 
 *	otherwise if res.type is e_cross, res.t1/t2 means radio/segment time of intersect point
 *------------------------------------------------------------------------------------------*/
typedef enum { e_cover, e_cross, e_passby } t_inter_type ;
struct t_seg_rel {
	t_inter_type type ;
	double t1, t2, cross, dot ; 
};
inline t_seg_rel get_cross(t_point const & p1, t_point const & p2, t_point const & q1, t_point const & q2)
{
	double dx_p , dx_q , dy_p , dy_q , dx_1 , dy_1 ;
	dx_p = p1.x-p2.x , dx_q = q1.x-q2.x ;
	dy_p = p1.y-p2.y , dy_q = q1.y-q2.y ;
	dx_1 = p1.x-q1.x , dy_1 = p1.y-q1.y ;
	double prod = dx_p*dy_q - dx_q*dy_p ;
	double prod_p = dx_1*dy_q - dy_1*dx_q ;

	t_seg_rel res ;
	if( fequal(prod,0) ){
		if( fequal(prod_p,0) ){
			double signal = dx_1*dx_p ;
			if( fequal(signal,0) )
				signal = dy_1*dy_p, res.t1=dy_1/dy_p, res.t2=(q2.y-p1.y)/dy_p ;
			else 
				res.t1 = dx_1/dx_p, res.t2 = (q2.x-p1.x)/dx_p ; 
			res.type = signal > 0 ? e_cover : e_passby ;
		} else 
			res.type=e_passby ;
	} else {
		res.t1 = prod_p / prod ;
		res.t2 = (dx_1*dy_p - dy_1*dx_p)/prod ;
		if( res.t1<=0 || res.t2<0 || res.t2>1 )
			res.type = e_passby ;
		else {
			res.type = e_cross ;
			double len=sqrt( dx_q*dx_q + dy_q*dy_q ) ;
			res.cross = prod / len   ;
			res.dot = (dx_p*dx_q + dy_p*dy_q)/len ;
			if( res.t2==0 )  
				res.cross = -res.cross, res.dot=-res.dot ;
		}
	}
	return res ;
}

int n_pt , visible[_max_points_], pt_no[_max_points_], n_vis, vis_list[_max_points_] ;
int edge_comp( const void *e1, const void *e2 ){
	int e1_v1=pt_no[*(int*)e1], e1_v2=pt_no[(*(int*)e1)+1] , tmp ;
	if( e1_v1 > e1_v2 )
		tmp=e1_v1, e1_v1=e1_v2, e1_v2=tmp ;
	int e2_v1=pt_no[*(int*)e2], e2_v2=pt_no[(*(int*)e2)+1] ;
	if( e2_v1 > e2_v2 )
		tmp=e2_v1, e2_v1=e2_v2, e2_v2=tmp ;
	return e1_v2==e2_v2 ? e1_v1-e2_v1 : e1_v2-e2_v2 ;
}

int main() {
	freopen("fence4.in", "r", stdin) ;
	freopen("fence4.out", "w", stdout) ;
	cin >> n_pt >> observer.x >> observer.y ;
	for( int i=0; i<n_pt; i++ ){
		cin >> poly[i].x >> poly[i].y ;
		for( int j=0; j<=i; j++ )
			if( poly[i]==poly[j] ) {
				pt_no[i]=j;
				break ;
			}
	}
	poly[n_pt] = poly[0] ;

	memset( visible, 0, sizeof(visible) );
	register int radial_pt, seg_pt ;
	for( radial_pt=0; radial_pt<n_pt; radial_pt++ ){
		register double left_t=-1, right_t=-1, left_dot, right_dot ;
		int left_seg=-1,  right_seg = -1 ;
		for( seg_pt=0; seg_pt<n_pt; seg_pt++ ){
			t_seg_rel res=get_cross( observer, poly[radial_pt], poly[seg_pt], poly[seg_pt+1] ) ;
			if( res.type==e_passby ){
				continue ;
			}
			if( res.type==e_cover && (res.t1<1 || res.t2<1) )
				break ;
			else if( res.type==e_cross ) {
				if( res.t1<1 ) 
					break ;
				if( res.t2>0 && res.t2<1 && res.cross>0 )
					res.cross=-res.cross, res.dot=-res.dot ;
				if( res.cross<0 && (left_t==-1 ||  res.t1<left_t || (res.t1==left_t && left_dot<res.dot)) ) 
					left_t = res.t1 , left_dot = res.dot , left_seg = seg_pt ;
				if( res.t2>0 && res.t2<1 )
					res.cross=-res.cross, res.dot=-res.dot ;
				if ( res.cross>0 && (right_t==-1 || res.t1<right_t || ( res.t1==right_t && right_dot<res.dot)) )
					right_t = res.t1 , right_dot = res.dot , right_seg = seg_pt ;
			}
		}
		if( seg_pt!=n_pt )
			continue ;
		if( left_seg!=-1 )
			visible[left_seg]=1 ;
		if( right_seg!=-1 )
			visible[right_seg]=1 ;
	}
	for( int i=0; i<n_pt; i++ )
		vis_list[n_vis]=i, n_vis +=visible[i] ;
	qsort( vis_list, n_vis, sizeof(int), edge_comp );
	cout << n_vis << endl ;
	for( int i=0; i<n_vis; i++ ){
		int s=vis_list[i] , t= vis_list[i]+1;
		if( pt_no[vis_list[i]] > pt_no[vis_list[i]+1] )
			s=vis_list[i]+1 , t= vis_list[i] ;
		cout << poly[s].x << " " << poly[s].y << " " << poly[t].x << " " << poly[t].y << endl;
	}
	return 0;
}