TopCoderOpen 2017 Round1A Hard PolygonRotation

あれ、今回Parallelないの？
https://community.topcoder.com/stat?c=problem_statement&pm=10691

問題

2次元の多角形が与えられる。
この図形をY軸を中心に回転させた図形の体積を求めよ。
なお、多角形は以下の条件を満たす。

• 凸である。
• 各頂点は格子点上にある。
• 図形のY座標の最大値をYmax、最小値をYminとすると、(0,Ymax)、(0,Ymin)には頂点がある。

解法

この図形のうち、X座標が負の部分をY軸に線対称に反転させた図形を考えよう。
元々X座標が正の部分との論理和を取った図形を回転させることになる。

そのような図形を直接求めてもよいだろうが、自分は別のアプローチで解いた。
点が格子点であることから、図形をY座標1毎に分割しよう。
すると、個々の分割領域は、1つの台形か、2つの台形を合わせた形状となる。
それぞれの形状を求め、回転させればよい。

回転後の図形の体積は、パップス=ギュルダンで求めてもよいが、重心を求めるのが手間だったので普通に数値積分した。

class PolygonRotation {
	public:
	double LX[202];
	double RX[202];
	
	double integral(double a,double b,double t) {
		double pi=atan(1)*4;
		
		return pi*(a*a*t*t*t/3.0+a*b*t*t+b*b*t);
	}
	
	double getVolume(vector <int> x, vector <int> y) {
		int N=x.size();
		int yma,ymi;
		int i;
		FOR(i,N) y[i]+=100;
		FOR(i,N) {
			if(i==0) yma=y[i];
			else if(x[i]>0) {
				if(y[i]==y[i-1]) RX[y[i]]=x[i];
				else {
					for(int ty=y[i-1]-1;ty>=y[i];ty--) RX[ty] = (x[i-1]*(ty-y[i])+x[i]*(y[i-1]-ty))*1.0/(y[i-1]-y[i]);
				}
			}
			else if(x[i]==0) {
				ymi=y[i];
				if(y[i-1]!=y[i]) {
					for(int ty=y[i-1]-1;ty>=y[i];ty--) RX[ty] = (x[i-1]*(ty-y[i])+x[i]*(y[i-1]-ty))*1.0/(y[i-1]-y[i]);
				}
			}
			else if(x[i]<0) {
				if(y[i]==y[i-1]) LX[y[i]]=-x[i];
				else {
					for(int ty=y[i-1]+1;ty<=y[i];ty++) LX[ty] = -(x[i-1]*(y[i]-ty)+x[i]*(ty-y[i-1]))*1.0/(y[i]-y[i-1]);
				}
			}
		}
		if(y[0]!=y[N-1]) {
			for(int ty=y[N-1]+1;ty<=y[0];ty++) LX[ty] = -(x[N-1]*(y[0]-ty)+x[0]*(ty-y[N-1]))*1.0/(y[0]-y[N-1]);
		}
		
		double ret=0;
		for(int y=ymi;y<yma;y++) {
			if(LX[y]>=RX[y] && LX[y+1]>=RX[y+1]) {
				ret += integral(LX[y+1]-LX[y],LX[y],1);
			}
			else if(RX[y]>=LX[y] && RX[y+1]>=LX[y+1]) {
				ret += integral(RX[y+1]-RX[y],RX[y],1);
			}
			else if(LX[y]>RX[y] && LX[y+1]<RX[y+1]) {
				double a=LX[y]-RX[y];
				double b=RX[y+1]-LX[y+1];
				double p=a/(a+b);
				ret += integral(LX[y+1]-LX[y],LX[y],p);
				ret += integral(RX[y]-RX[y+1],RX[y+1],1-p);
			}
			else if(LX[y]<RX[y] && LX[y+1]>RX[y+1]) {
				double a=RX[y]-LX[y];
				double b=LX[y+1]-RX[y+1];
				double p=a/(a+b);
				ret += integral(RX[y+1]-RX[y],RX[y],p);
				ret += integral(LX[y]-LX[y+1],LX[y+1],1-p);
			}
		}
		
		return ret;
		
	}
}

まとめ

今回参加者の大半が正の点を取ってRound1通過してしまったけど、次回以降参加者確保できるのかな…。