1.1.1.1 凸包の計算と多角形の面積計算
問題 5‑18 猫の行動範囲問題
猫の出没する範囲を調べました。猫が出没する場所を最小の面積で囲む境界を考え、その面積を求めてください。
入力として、最初に頂点の数が与えられます。頂点の数は1から100までです。次に、頂点の座標が与えられます。頂点の座標は空白と改行で区切られます。
出力として、猫の行動範囲の面積を小数点以下2桁までの精度で示してください。
サンプル入力
14
0 0
3 0
4 0
5 1
4 2
5 3
6 6
4 4
0 4
0 6
6 0
2 2
0 7
1 6
サンプル出力
39.00
まず、猫が出没する場所を最小の面積で囲む境界を求めなくてはなりません。このような境界のことを凸包(convex hull)といいます。
なお、凸包は面積が最小であるだけでなく、面積を囲うための最小の長さの線分となります。凸包は、すべての点を囲む最小の凸多角形です。
凸包を求める包装アルゴリズム(package wrapping algorithm)を次に示します。
１．y座標が最小の点を求めます。y座標が同じ点があったときは、x座標が大きいほうを選びます。（直感的には右下の点を選ぶということです。）
２．nowを0度にします。求めた点をp1とします。
３．すべての点とp1について、水平線と成す角を求めます。その中で、now以上で最小の角度を持つ点を選びます。同じ角度の点があったときは、なるべく遠くの点を選びます。
４．点が見つからないか、最初の点に戻ってきたら終了します。
５．３で選んだ点を凸包に加えます。そして、nowをその点と水平線が成す角度とします。その後、選んだ点を３で調べる点の候補からはずし、点を新たなp1とします。
６．３へ戻ります。
このアルゴリズムで、y座標の最小の点を求めるのは、水平線と成す角度が最小のものをみつけるためです。
イメージを説明します。最初に水平で右に向かっているベクトルがあるとします。一番下にある点を基準として、そのベクトルと水平線との成す角を徐々に増やしていきます。増やしていくと、やがて別の点とベクトルがぶつかります。その点は、凸包を構成する点です。次に、その点を中心として、また角度を増やしていきます。また点とぶつかります。その点も凸包を構成する点です。これを最初の点に戻ってくるまで繰り返します。
リスト 5‑37 包装アルゴリズムのプログラム
class TPoint {
public:
int x,y;
TPoint( int x1=0, int y1=0 ) : x(x1), y(y1) {};
friend bool operator == ( const TPoint &p1, const TPoint &p2 ) {
return p1.x == p2.x && p1.y == p2.y;
}
friend bool operator != ( const TPoint &p1, const TPoint &p2 ) {
return !(p1 == p2);
}
};
//------------------------------------------------------
double getTheta(const TPoint &p1, const TPoint &p2) {
double dx, dy, ax, ay;
double t;
dx = p2.x - p1.x; ax = dx >= 0 ? dx : -dx;
dy = p2.y - p1.y; ay = dy >= 0 ? dy : -dy;
t = (ax+ay == 0) ? 0 : (double) dy/(ax+ay);
if ( dx < 0) t = 2-t; else if (dy < 0) t = 4+t;
return t*90.0;
}
//------------------------------------------------------
double distance( const TPoint &p1, const TPoint &p2 ) {
double dx = (p2.x-p1.x);
double dy = (p2.y-p1.y);
return sqrt( dx*dx+ dy*dy );
}
//------------------------------------------------------
void makeConvexHull( TPoint poi[], int &poiNum ) {
int minIndex;
int minY, maxX;
int i,j;
minY = poi[0].y; maxX = poi[0].x;
minIndex = 0;
// x座標最大、y座標最小の点を探す
for( i=1; i<poiNum; i++ ) {
if ( minY > poi[i].y || ( minY == poi[i].y && maxX < poi[i].x )) {
maxX = poi[i].x; minY = poi[i].y;
minIndex = i;
}
}
double rad, find;
double now = 0;
double dist = 0;
// 番兵のセット
poi[ poiNum ] = poi[ minIndex ];
for( j=0; j<poiNum; j++ ) {
// 先頭に見つけた要素を持ってくる
TPoint tmp = poi[ j ];
poi[ j ] = poi[ minIndex ];
poi[ minIndex ] = tmp;
find = 360; // 絶対にない値
for( i=j+1; i<=poiNum; i++ ) {
// 前の点とでできる線分の成す角を求める
rad = getTheta( poi[j], poi[i] );
if ( poi[j] != poi[i] && find >= rad && rad >= now ) {
assert( rad >= now );
// 同じ角度だったときは、最も遠くの点を選ぶ
if ( find != rad || find == rad && dist < distance( poi[j], poi[i]) ) {
dist = distance( poi[j], poi[i] );
find = rad;
minIndex = i;
}
}
}
// 一つも見つからないか、番兵に行ったら終了
if ( find == 360 || minIndex == poiNum )
break;
now = find;
}
j++;
// 末尾と先頭が同じ点だったら、削除する。
while( poi[0] == poi[j-1] )
j--;
// 最終的な点の数をセットする。
poiNum = j;
}
番兵というものを用います。番兵として先頭と同じ要素を末尾にも置いておくことで、配列を超えて探索していないかのチェックを省くことができ、高速化できます。
なお、配列poiは、番兵の分だけ１つ余分に領域を確保しておかねばなりません。また求めた凸包は、引数の配列を上書きして保存されます。
getTheta関数は、擬似角度を求める関数です。この関数を使うと、正確な角度は求まりませんが、どちらの角の方が大きいかなどを判定することができます。角度の大小関係はちゃんと成り立っています。
次に凸包が求まったとして、その凸包の面積を求めるプログラムを次に示します。
リスト 5‑38 多角形の面積を求めるプログラム
double sumPolygonArea( const TPoint poi[], const int poiNum ) {
int i;
double sum;
sum = 0;
for( i=0; i< poiNum; i++ ) {
sum += poi[i].x*poi[i+1].y - poi[i+1].x*poi[i].y;
}
sum += poi[i].x*poi[0].y - poi[0].x*poi[i].y;
return sum / 2;
}
引数は頂点の座標と頂点の数です。
なお、多角形は凸多角形である必要はなく、任意の多角形の面積を求めることができます。多角形の頂点が左回りに配置されているときは、計算結果はプラスの値になります。右回りに配置されているときは、計算結果はマイナスの値となってしまいますが、絶対値をとればそれが求める多角形の面積です。
このプログラムがなぜうまくいくのか、私にはまだ理解できていません。しかし、簡単なプログラムのため、暗記しています。
リスト 5‑37とリスト 5‑38を用いて問題を解いたプログラムを示します。
リスト 5‑39 猫の行動問題の解プログラム
#include <iostream.h>
#include <assert.h>
#include <iomanip.h>
// 点を表す構造体
class TPoint {
public:
int x,y;
TPoint( int x1=0, int y1=0 ) : x(x1), y(y1) {};
friend bool operator == ( const TPoint &p1, const TPoint &p2 ) {
return p1.x == p2.x && p1.y == p2.y;
}
friend bool operator != ( const TPoint &p1, const TPoint &p2 ) {
return !(p1 == p2);
}
};
//------------------------------------------------------
class TCalc {
TPoint point[101]; // 頂点の座標
int pointNum; // 頂点の数
public:
void calc();
void input( istream &is );
};
//------------------------------------------------------
// 擬似角度を求める
double getTheta(const TPoint &p1, const TPoint &p2) {
int dx, dy, ax, ay;
double t;
dx = p2.x - p1.x; ax = dx >= 0 ? dx : -dx;
dy = p2.y - p1.y; ay = dy >= 0 ? dy : -dy;
t = (ax+ay == 0) ? 0 : (double) dy/(ax+ay);
if ( dx < 0) t = 2-t; else if (dy < 0) t = 4+t;
return t*90.0;
}
//------------------------------------------------------
// 点間の距離の計算
double distance( const TPoint &p1, const TPoint &p2 ) {
double dx = (p2.x-p1.x);
double dy = (p2.y-p1.y);
return sqrt( dx*dx+ dy*dy );
}
//------------------------------------------------------
// 凸包を求める関数
void makeConvexHull( TPoint poi[], int &poiNum ) {
int minIndex;
int minY, maxX;
int i,j;
minY = poi[0].y; maxX = poi[0].x;
minIndex = 0;
// x座標最大、y座標最小の点を探す
for( i=1; i<poiNum; i++ ) {
if ( minY > poi[i].y || ( minY == poi[i].y && maxX < poi[i].x )) {
maxX = poi[i].x; minY = poi[i].y;
minIndex = i;
}
}
double rad, find;
double now = 0;
double dist = 0;
// 番兵のセット
poi[ poiNum ] = poi[ minIndex ];
for( j=0; j<poiNum; j++ ) {
// 先頭に見つけた要素を持ってくる
TPoint tmp = poi[ j ];
poi[ j ] = poi[ minIndex ];
poi[ minIndex ] = tmp;
find = 360; // 絶対にない値
for( i=j+1; i<=poiNum; i++ ) {
// 前の点とでできる線分の成す角を求める
rad = getTheta( poi[j], poi[i] );
if ( poi[j] != poi[i] && find >= rad && rad >= now ) {
assert( rad >= now );
// 同じ角度だったときは、最も遠くの点を選ぶ
if ( find != rad || find == rad && dist < distance( poi[j], poi[i]) ) {
dist = distance( poi[j], poi[i] );
find = rad;
minIndex = i;
}
}
}
// 一つも見つからないか、番兵に行ったら終了
if ( find == 360 || minIndex == poiNum )
break;
now = find;
}
j++;
// 末尾と先頭が同じ点だったら、削除する。
while( poi[0] == poi[j-1] )
j--;
// 最終的な点の数をセットする。
poiNum = j;
}
//------------------------------------------------------
// 多角形の面積を求める
double sumPolygonArea( const TPoint poi[], const int poiNum ) {
int i;
double sum;
sum = 0;
for( i=0; i< poiNum-1; i++ ) {
sum += poi[i].x*poi[i+1].y - poi[i+1].x*poi[i].y;
}
sum += poi[i].x*poi[0].y - poi[0].x*poi[i].y;
return sum / 2;
}
//------------------------------------------------------
void TCalc::calc() {
makeConvexHull( point, pointNum );
cout << setprecision(2) << setiosflags(ios::fixed);
cout << sumPolygonArea( point, pointNum ) << endl;
}
//------------------------------------------------------
void TCalc::input( istream &is ) {
int i;
is >> pointNum;
for( i=0; i<pointNum; i++ )
is >> point[i].x >> point[i].y;
}
//------------------------------------------------------
int main() {
TCalc calc1;
calc1.input( cin );
calc1.calc();
return 0;
}
#include <iostream.h>
#include <assert.h>
#include <iomanip.h>
// 点を表す構造体
struct TPoint {
int x,y;
TPoint( int x1=0, int y1=0 ) : x(x1), y(y1) {};
friend bool operator == ( const TPoint &p1, const TPoint &p2 ) {
return p1.x == p2.x && p1.y == p2.y;
}
friend bool operator != ( const TPoint &p1, const TPoint &p2 ) {
return !(p1 == p2);
}
};
//------------------------------------------------------
TPoint g_point[101]; // 頂点の座標
int g_pointNum; // 頂点の数
//------------------------------------------------------
// 擬似角度を求める
double getTheta(const TPoint &p1, const TPoint &p2) {
int dx, dy, ax, ay;
double t;
dx = p2.x - p1.x; ax = dx >= 0 ? dx : -dx;
dy = p2.y - p1.y; ay = dy >= 0 ? dy : -dy;
t = (ax+ay == 0) ? 0 : (double) dy/(ax+ay);
if ( dx < 0) t = 2-t; else if (dy < 0) t = 4+t;
return t*90.0;
}
//------------------------------------------------------
// 点間の距離の計算
double distance( const TPoint &p1, const TPoint &p2 ) {
double dx = (p2.x-p1.x);
double dy = (p2.y-p1.y);
return sqrt( dx*dx+ dy*dy );
}
//------------------------------------------------------
// 凸包を求める関数
void makeConvexHull( TPoint poi[], int &poiNum ) {
int minIndex;
int minY, maxX;
int i,j;
minY = poi[0].y; maxX = poi[0].x;
minIndex = 0;
// x座標最大、y座標最小の点を探す
for( i=1; i<poiNum; i++ ) {
if ( minY > poi[i].y || ( minY == poi[i].y && maxX < poi[i].x )) {
maxX = poi[i].x; minY = poi[i].y;
minIndex = i;
}
}
double rad, find;
double now = 0;
double dist = 0;
// 番兵のセット
poi[ poiNum ] = poi[ minIndex ];
for( j=0; j<poiNum; j++ ) {
// 先頭に見つけた要素を持ってくる
TPoint tmp = poi[ j ];
poi[ j ] = poi[ minIndex ];
poi[ minIndex ] = tmp;
find = 360; // 絶対にない値
for( i=j+1; i<=poiNum; i++ ) {
// 前の点とでできる線分の成す角を求める
rad = getTheta( poi[j], poi[i] );
if ( poi[j] != poi[i] && find >= rad && rad >= now ) {
assert( rad >= now );
// 同じ角度だったときは、最も遠くの点を選ぶ
if ( find != rad || find == rad && dist < distance( poi[j], poi[i]) ) {
dist = distance( poi[j], poi[i] );
find = rad;
minIndex = i;
}
}
}
// 一つも見つからないか、番兵に行ったら終了
if ( find == 360 || minIndex == poiNum )
break;
now = find;
}
j++;
// 末尾と先頭が同じ点だったら、削除する。
while( poi[0] == poi[j-1] )
j--;
// 最終的な点の数をセットする。
poiNum = j;
}
//------------------------------------------------------
// 多角形の面積を求める
double sumPolygonArea( const TPoint poi[], const int poiNum ) {
int i;
double sum;
sum = 0;
for( i=0; i< poiNum-1; i++ ) {
sum += poi[i].x*poi[i+1].y - poi[i+1].x*poi[i].y;
}
sum += poi[i].x*poi[0].y - poi[0].x*poi[i].y;
return sum / 2;
}
//------------------------------------------------------
void calc() {
makeConvexHull( g_point, g_pointNum );
cout << setprecision(2) << setiosflags(ios::fixed);
cout << sumPolygonArea( g_point, g_pointNum ) << endl;
}
//------------------------------------------------------
void input() {
int i;
cin >> g_pointNum;
for( i=0; i<g_pointNum; i++ )
cin >> g_point[i].x >> g_point[i].y;
}
//------------------------------------------------------
int main() {
input();
calc();
return 0;
}
凸包を求める方法、多角形の面積を求める方法の２つの方法を示しました。
なお、点を扱う関数を書くときは、点の座標をdoubleに入れるか、intに値を入れるかを注意してください。