点と多角形の内外判定（Winding-Numberアルゴリズム）
概要
点が多角形に含まれているか否かはWinding Numberアルゴリズムを使って判定できる。調べる点から見て多角形の各辺がなす角度を、反時計回りを正・時計回りを負として一周分合計する。合計が内側にある場合は360度になり、外側の場合は0度になる。
コード
調べるポイントごとの処理。input1には多角形のポリラインを入力する。多角形は XZ 平面にあるものとし、上向きを +Y として角度の正負を決めている。ポリラインは閉じていなくてもよい（pts[i-1] は i=0 のとき pts[-1]、つまり末尾の点になるので、コード側で最後の点と最初の点をつないで閉じている）。
//
// Winding Number Algorithm
// Run Over: Points
//
vector up = set(0, 1, 0); // 多角形は XZ 平面にある前提
float result = 0;
int pts[] = primpoints(1, 0);
for(int i = 0; i < len(pts); i++)
{
// i == 0 のとき pts[-1] は末尾の点。これで多角形を閉じる
vector p0 = point(1, "P", pts[i-1]);
vector p1 = point(1, "P", pts[i]);
vector l0 = p0 - @P;
vector l1 = p1 - @P;
// 調べる点から見た辺 p0→p1 の角度。
// 反時計回りなら正、時計回りなら負になる（atan2 なので acos のような NaN は出ない）
float angle = atan2(dot(cross(l0, l1), up), dot(l0, l1));
result += angle;
}
// 合計は 0 か ±2π（自己交差があれば ±4π など）にしかならないので PI で分ける
if(abs(result) > PI)
@Cd = set(1,0.5,0);
else
@Cd = set(0,0.5,1);多角形が自己交差している（8の字のようにねじれている）と、領域によって合計が +360度 になったり −360度 になったりする。そのため abs() で絶対値をとってから判定している。調べる点が頂点や辺の上に乗っている場合は l0 か l1 がゼロベクトルになり、角度が定まらないことに注意。