点と多角形の内外判定（Winding-Numberアルゴリズム）

Winding Numberアルゴリズムを使う。多角形の各頂点を周回して得られる角度の合計が、もし内側にある場合は360度になり、外側の場合は0度になる。これを利用して判定する。

//
// Winding Number Algorithm
// 
// Run Over: Points
//
float result = 0;

int pts[] = primpoints(1, 0);
for(int i = 0; i < len(pts); i++)
{
    vector p0 = point(1, "P", pts[i-1]);
    vector p1 = point(1, "P", pts[i]);
    
    vector l0 = p0 - @P;
    vector l1 = p1 - @P;
    
    float angle = acos(dot(normalize(l1), normalize(l0)));
    
    // 進行方向が逆転しているか外積でチェックする
    vector cross = cross(l0, l1);
    if(dot(cross, set(0,1,0)) < 0)
        angle *= -1;
    
    result += angle;
}

if(abs(result) >= 0.01)
    @Cd = set(1,0.5,0);
else
    @Cd = set(0,0.5,1);

ポリラインにねじれがあると進行方向が逆転して正負が変わる。

国境だけのポリラインと、国名を含む首都のポイントをお互いに突き合わせて、国境に国名を設定しているコード。

//
// Winding Number Algorithm
// 
// Run Over: Primitive
//

int pts[] = primpoints(0, @primnum);

for(int i = 0; i < npoints(1); i++)
{
    float result = 0;

    vector center = point(1, "P", i);

    for(int j = 0; j < len(pts); j++)
    {
        vector p1 = point(0, "P", pts[j-1]);
        vector p2 = point(0, "P", pts[j]);
        
        vector l1 = p1 - center;
        vector l2 = p2 - center;
        
        float angle = acos(dot(normalize(l2), normalize(l1)));
        
        vector cross = cross(l1, l2);
        if(dot(cross, set(0,1,0)) < 0)
            angle *= -1;
        
        result += angle;
    }
    
    if(abs(result) >= 0.01)
    {
        s@name = point(1, "Country", i);
        break;
    }
}