Help us understand the problem. What is going on with this article?

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

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

判定アルゴリズムは、その点から正の方向に伸ばした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) {
  var windingNumber=0,x=point[0],y=point[1];
  var i=poly.length-1; // indexLast(poly)
  for (var 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))
            ==(yj*(poly[i][0]-x)>yi*(poly[j][0]-x)))
      windingNumber+=counterClockwise?1:-1;
  }
  return windingNumber%2==1; // ray casting algorithm
}
//  return windingNumber!==0;  // winding number algorithm

var points = [], 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 (var j=0; j<=7; j++) {
  for (var i=0; i<=7; i++) {
    var p = [i, j];
    if (pointInPolygon(p, poly)) {points.push(p)}
  }
}

function mytransform(p) {
    var scale = 32, offset = 4;
    return [p[0]*scale+offset, (7-p[1])*scale+offset]
}

var 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 の利用へ置き換え可能です。

// var 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) {
  var 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. 複数の多角形領域で平面を余さず分割している場合は、与えられた点の一意の帰属判定にも使えます。 

kkdd
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
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  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
ユーザーは見つかりませんでした