Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationEventAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
19
Help us understand the problem. What are the problem?
@kkdd

二次元多角形領域内部に与点が含まれるかどうかを判定

こんにちは。
二次元平面で、与えられた点が、多角形領域の内部に含まれるかどうかを判定しました(内外判定)。下記の pointInPolygon() 関数です。

判定アルゴリズムは、その点から正の方向に伸ばしたx軸(半直線)が多角形境界線と交差する回数を符号付きで数えます1 2 3 4。下記の二種類の条件を与えました。

実行例で、内部と判定された点をプロットしました(下図内の赤い点)。
ray_casting.png
winding_number.png

pointInPolygon.html
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8">
<title>pointInPolygon</title>
<script src="http://d3js.org/d3.v3.min.js"></script>
</head>
<body>
ray casting algorithm
<div id="result"></div>
<script type="text/javascript">
function pointInPolygon(point, poly) {
  let windingNumber=0,x=point[0],y=point[1];
  let i=poly.length-1; // indexLast(poly)
  for (let j=0,yi,yj,counterClockwise; i>=0; j=i--) {
    yi=poly[i][1]-y; yj=poly[j][1]-y;
    if ((yi>0)!=(yj>0) &&
      (counterClockwise=(yj>yi)?1:-1)*(yj*(poly[i][0]-x)-yi*(poly[j][0]-x))>0)
      windingNumber+=counterClockwise;
  }
  return windingNumber%2==1; // ray casting algorithm
}
//  return windingNumber!==0;  // winding number algorithm

const points = [], max = 7, poly = [[0,0], [7,0], [7,2], [7,3], [5,5], [2,5], [2,3], [6,3], [6,1], [4,1], [4,7], [2,7], [0,5]];
for (let j=0; j<=max; j++) {
  for (let i=0; i<=max; i++) {
    const p = [i, j];
    if (pointInPolygon(p, poly)) {points.push(p)}
  }
}

function mytransform(p) {
    const scale = 32;
    return [(p[0]+1)*scale, (max-p[1]+1)*scale]
}

const svg = d3.select("#result")
    .append("svg")
    .attr("width",300)
    .attr("height",300);

svg.append("polygon")
    .data([poly.map(mytransform)])
    .attr("points", function(d){return d.join(" ")})
    .attr("stroke", "gray")
    .attr("fill", "#fff0e0")
    .attr("stroke-width", 1.5);

svg.selectAll("points")
    .data(points.map(mytransform))
    .enter().append("circle")
    .attr("fill", "#ff4040")
    .attr("r", 2.5)
    .attr("transform", function(d){return `translate(${d})`});
</script>
</body>
</html>

補足

なお ray casting algorithm の方ならば、下記のように偶奇を判定する Boolean 型変数 inside の利用へ置き換え可能です。

// let inside = false; // initialization
if ((yi>0)!=(yj>0) && (yj>yi)==(yj*(poly[i][0]-x)>yi*(poly[j][0]-x)))
      inside = !inside;
// return inside; // ray casting algorithm

また今回の用途には、sign(a, b)は下記の定義を使わず、上記のように (a > b) の値を使えば十分です。

function sign(a, b) {
    return (a > b) ? 1 : (a < b) ? -1 : 0;
}

また、indexLast(poly)は下記の定義を使わず、上記のように poly.length-1 の値を使えば十分です。

function indexLast(arr) {
  let last = arr.length - 1;
  if (arr[0][0]==arr[last][0]&&arr[0][1]==arr[last][1]) last -= 1;
  return last;
}

  1. ソースコードから分かるように、winding number を求める計算も、ray casting algorithm とほぼ同じアルゴリズムに帰着でき、角度計算も不要で極めて効率が良いです。  

  2. 参考記事として、"Inclusion of a point in a polygon" (Dan Sunday)、"point-in-polygon" (GitHub)。 

  3. ただし多角形の全辺を走査することになります。辺数が多い多角形に対して、多数の与点を判定したい場合は、計算量低減(毎回全辺を走査しない)を考える必要があります。例えば、"Expedicious and Exact Extracts with Osmium" (Jochen Topf's Blog)。 

  4. 左側の辺、および下側の水平な辺上も内部と判定されます。 

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
19
Help us understand the problem. What are the problem?