時速3kmと時速5kmの平均
僕が小学校のころに塾の先生に出題された問題で、
「A君は家から学校まで行きは時速3kmで行き、帰りは時速5kmで帰りました。A君は平均時速何kmで行ったでしょう」
というのがあった。もちろん引っ掛け問題ではないので、行きと全く同じ道を辿って帰ったとしてもらいたい。
初めのうちは「え?家から学校までの距離は?」と思ったが、一向にその情報が与えられない。手がかりがつかめない。どこから手をつけていいかわからない。しばし途方にくれたあげく、最終的に「そうか、情報が与えられないと言うことは、別に距離は何だって同じだということだな」と悟り、勝手に15kmと仮定して答えを出した覚えがある。今から考えると危険である。
それでも「時速4km」じゃないことに納得いかない人のために、最近になって、たぶん一番わかりやすい説明を考えた。
ちなみに正しくは3と5の「調和平均」と言って、「逆数を(相加)平均したものの逆数」、つまり、 である。
まず確認。平均の速さは「全体の距離」÷「全体の時間」である。20kmの道を5時間かけて行くなら、途中どれだけ走ろうが休もうが、平均時速は4kmである。
時速3kmと時速5kmの平均が時速4kmになるのは、例えば以下の場合である。
「A君は時速3kmで1時間行き、そのあと時速5kmで1時間行きました。平均時速は何kmでしょう」
この場合は、つまりは2時間で8km行くわけで、平均時速は kmとなる。時間が一緒ならば普通の(相加)平均でいいわけだ。
これが時間が違ってくると、例えば問題が
「A君は時速3kmで2時間行き、そのあと時速5kmで1時間行きました。平均時速は何kmでしょう」
となると、つまりは3時間で3*2+5=11km行くわけで、平均時速は kmとなる。距離の平均は時間の比率によって加重がかかるのである。
冒頭の問題では、「距離が等しい」ということは、時間の比は速さの比の逆比である。
「A君は時速3kmで5時間行き、そのあと時速5kmで3時間行きました。平均時速は何kmでしょう」
と同じなので、平均時速は kmとなる。
一般に速さと速さで、それぞれ同じ距離を行く場合、その時間の比はなので、その平均の速さは である。
はじめて手をつけるVisualC#
簡易更新。
今日(昨日)はじめて Visual C# に手をつける。VisualStudio2008Expressをネットダウンロード。仕事で使わない勉強なので、延び延びになっていた。
少しやった感想は、Windowsプログラミングにおいて、C#はめちゃ便利だな、と。Ruby/SDLで何かと大変だった直線描画・多角形描画もとても簡単。Graphicsクラスが既にあるから。
参考:
http://ufcpp.net/study/csharp/lib_drawing.html
http://dobon.net/vb/dotnet/graphics/drawline.html
そうとは知らず、同じ計算幾何をやらせている後輩(仕事でC#もしている)に「描画できるだろうか」と勝手に心配して、適当に調べた知識をメールで送った自分のピエロっぷりに泣いた。恥ずかしい。
これからの自分の計算幾何の勉強をRubyで進めるかC#で進めるかは、決めかねています。次の更新で明らかになるとおもいます。
それにしても、C#クラスの見やすい網羅型リファレンスが欲しいなあ・・・
Ruby/SDLで計算幾何(5) 点が多角形の中に内部にあるか判定する2
その前に前の日記のPolygonクラス絡みで拡張したいところやらバグやらがあったので修正します。
class GeoBase class Polygon attr_accessor :vertices, :size def initialize(points) @size = points.size # 頂点数を追加 @vertices = points.push points[0], points[1] end end def polygon(x1, *arg) if arg if arg.size % 2 == 0 raise ArgumentError, 'num of coodinates is odd.' else arg.unshift x1 vertices = [] (0..(arg.size/2-1)).each do |i| # このあたり修正 vertices << p(arg[2*i], arg[2*i+1]) end end else vertices = x1 end Polygon.new vertices end alias poly polygon end
そしてPointクラスのメソッドinsideです。もっとスマートに書けそうなものですが・・・
class GeoBase class Point def inside(poly) x_inf = Point.new(65535,y) #とりあえず lt= Line.new(self, x_inf) p = poly.vertices j = 0 count = 0 (1..poly.size).each do |i| lp=Line.new(p[i],p[i]) unless (GeoUtil.ccw(lt.p1, lt.p2, p[i]) == 0) if (i == j + 1) lp.p2 = p[j] count += 1 if GeoUtil.intersect(lp, lt) else count += 1 if (GeoUtil.ccw(lt.p1, lt.p2, p[i]) * GeoUtil.ccw(lt.p1, lt.p2, p[j])) < 0 end j = i end end (count&1 == 0)? false : true end end end
では、テストをしてみます。テストデータとして下図のような、凹凸のある七角形を用意します。
class GeoBaseTest < GeoBase def start poly = poly(50,50,80,100,50,150,100,120,150,150,120,100,150,50) draw poly def set_p(x, y) p1 = p(x,y) draw p1 p1 end puts set_p(100,100).inside(poly) #=> true puts set_p(80,120).inside(poly) #=> true puts set_p(50,100).inside(poly) #=> false end end GeoBaseTest.start
正常に動いているようです。
ところで、点が辺上にある場合、頂点である場合はどういう答えが返ってくるかを見ていませんでした。七頂点を設定してみます。
puts set_p(50,50).inside(poly) #=> true puts set_p(80,100).inside(poly) #=> false puts set_p(50,150).inside(poly) #=> false puts set_p(100,120).inside(poly) #=> true puts set_p(150,150).inside(poly) #=> false puts set_p(120,100).inside(poly) #=> true puts set_p(150,50).inside(poly) #=> true
黒点は多角形の内部、白点は多角形の外部を表しています。
実は、一点の例外を除き、これらの頂点の判定は、各頂点を少しx軸負方向(図でいう左側方向)にずらした点の判定と同一です。これは、ccw,intersectによる判定が、線分の端点を含むことに由来しています。
では、例外の一点(50,50)は何かと言いますと、、、配列の先頭の処理に由来すると思うのですが、、、疲れましたのでまた今度です。
Ruby/SDLで計算幾何(4) 点が多角形の中に内部にあるか判定する1
さて、『アルゴリズムC第二巻』ISBN:4764902567 p.185にある、点が多角形に包含されているか判定するコードを見ていきましょう。
判定したい点から無限遠(扱う領域の端)まで、一方向に線分を延ばします。ここではx軸正方向とします。
そして、その線分が多角形の辺、もしくは頂点と交差する回数を数え、その数が偶数なら多角形の外部、奇数なら多角形の内部、としたいのですが、これだけでは問題があるわけです。
上図の2点について、ともに延ばした線分が1辺と1頂点と交差しているので条件は同じはずなのですが、多角形の内部か外部かは異なっています。
辺については問題ないのですが、頂点については、その交わり方にはカウントしていいものといけないものがあります。ここでは上図のように、「横断的」「非横断的」と名付けておきます。*1「非横断的」な交差は交点としてカウントする必要はありませんので、その検出も必要です。
以上のことが意識された関数 inside (p.185) は以下です。ただ、自分なりにコーディングスタイルを変えたり、ちょこちょこコメントを付したりしています。*2
int inside ( struct point t, /* 点 */ struct point p[], /* 多角形の頂点集合 */ int N /* 多角形の頂点数 */ ) { int i; int count = 0; int j = 0; struct line lt, lp; p[0] = p[N]; p[N+1] = p[1]; /* ltの設定 */ lt.p1 = t; lt.p2 = t; lt.p2.x = INT_MAX; for (i = 1; i <= N; i++) { lp.p1 = p[i]; lp.p2 = p[i]; if (ccw(lt.p1, lt.p2, p[i]) != 0) /* p[i]がlt上に乗ってないとき */ { if (i == j+1) { lp.p2 = p[j]; /* p[i]p[j]は多角形の一辺 */ if (intersect(lp, lt)) count++; } else { /* p[i]とp[j]の間の頂点が、すべてlt上に乗っている */ /* 横断的か非横断的かを判定 */ if (ccw(lt.p1, lt.p2, p[i]) * ccw(lt.p1, lt.p2, p[j]) < 0) count++; } j = i; } } return count & 1; /* 偶奇を返す */ }
現在これをRubyに書き換えてる最中なのですが・・・さらっと出来ないのが当方のオツムの至らないところでございまして・・・。
Array#injectが使えるのかな・・・。
続きはまた今度です。
*1:ちなみにこの用語は、自分が学生時代にかじっていた微分トポロジーの用語「横断的に交わる」intersect transversally からの連想で名付けました。http://en.wikipedia.org/wiki/Transversal_intersection
*2:また、多角形にもちょっと条件が必要になってくるらしいのですが、そこの分析も今度です。
Ruby/SDLで計算幾何(3) 線分交差判定のディテール
intersectの謎の分析をしようと思いましたが、ccwの分析ができてればそう悩むこともなくてですね。
再掲。
module GeoUtil module_function # 線分p0p1を基準にして、線分p0p2が # 反時計回りの向きにあれば+1 # 時計回りの向きにあれば-1 def ccw(p0, p1, p2) dx1 = p1.x - p0.x; dy1 = p1.y - p0.y dx2 = p2.x - p0.x; dy2 = p2.y - p0.y return +1 if dx1*dy2 > dy1*dx2 return -1 if dx1*dy2 < dy1*dx2 return -1 if (dx1*dx2 < 0) || (dy1*dy2 < 0) #(1) return +1 if dx1*dx1+dy1*dy1 < dx2*dx2+dy2*dy2 #(2) return 0 #(3) end # 2線分が交差しているかを判定する def intersect(l1, l2) return ((ccw(l1.p1, l1.p2, l2.p1) * ccw(l1.p1, l1.p2, l2.p2)) <= 0 ) && ((ccw(l2.p1, l2.p2, l1.p1) * ccw(l2.p1, l2.p2, l1.p2)) <= 0 ) end end
「一方の線分から見て、他方の線分の両端が左右に分かれている」&&「他方から見ても同様」なら、交差しているとみなすのがintersectメソッド。
そして、他方の端点が一方の線分上にある場合は、これは交差しているとみなされます。
線分が同一直線上にある場合、端点位置のパターンとしては上図の3つしかないのですが、その場合の交差・非交差の状況は上の通り。
うん、直感には反してないと思います。
ただ、「一方の線分の端点が同一(線分が一点)の場合」は、個人的には直感に反する結果が出てきます。
class GeoBaseTest < GeoBase def start puts GeoUtil.intersect(l(30,30,30,30), l(20,20,40,40)) end def self.start geo = self.new geo.start end end GeoBaseTest.start #=> false
Ruby/SDLで計算幾何(2)線分が交差しているか判定する
とりあえずテキスト『アルゴリズムC第二巻』amazon:4764902567を写しました。
module GeoUtil module_function # 線分p0p1を基準にして、線分p0p2が # 反時計回りの向きにあれば+1 # 時計回りの向きにあれば-1 def ccw(p0, p1, p2) dx1 = p1.x - p0.x; dy1 = p1.y - p0.y dx2 = p2.x - p0.x; dy2 = p2.y - p0.y return +1 if dx1*dy2 > dy1*dx2 return -1 if dx1*dy2 < dy1*dx2 return -1 if (dx1*dx2 < 0) || (dy1*dy2 < 0) #(1) return +1 if dx1*dx1+dy1*dy1 < dx2*dx2+dy2*dy2 #(2) return 0 #(3) end # 2線分が交差しているかを判定する def intersect(l1, l2) return ((ccw(l1.p1, l1.p2, l2.p1) * ccw(l1.p1, l1.p2, l2.p2)) <= 0 ) && ((ccw(l2.p1, l2.p2, l1.p1) * ccw(l2.p1, l2.p2, l1.p2)) <= 0 ) end end
このccwについて、なぜこういう関数にしてあるのか分析します。まだ途中です。
(1)-(3)はp0,p1,p2の3点が同一直線上にある場合(もしくは3点が同一の場合)の処理をしている箇所です。
以下、p0,p1,p2の3点が同一直線上にあり、p0とp1が異なる点であるとします。(つまり、線分p0p1が存在するとします)
- (1)により値が-1となるのは、p0から見てp2がp1と異なる側にある場合、つまりp2-p0-p1の順に並んでいる場合です。
- xについての条件とyについての条件がorで結びついているのは、直線がx軸に平行、またはy軸に平行な場合を考慮してのことです。
- (2)により値が+1となるのは、p0から見てp2がp1の向こう側にある場合、つまりp0-p1-p2の順に並んでいる場合です。
- (3)により値が0となるのは、p2が線分p0p1上にある場合、つまりp0-p2-p1の順に並んでいる場合です。これは両端点も含みます。
p0とp1が同じ点ならどうでしょうか。関数に dx1=dy1=0 を代入して追ってみますと、上の結果から類推できるような結果が出ます。実行した方が早いかもしれません。
このようなccwの仕様を、intersectにどう活かしているかは、時間がありませんので明日。
Ruby/SDLで計算幾何(1)線分が交差しているか判定する
とりあえずテキスト『アルゴリズムC第二巻』amazon:4764902567を写しました。
module GeoUtil module_function # 線分p0p1を基準にして、線分p0p2が # 反時計回りの向きにあれば+1 # 時計回りの向きにあれば-1 def ccw(p0, p1, p2) dx1 = p1.x - p0.x; dy1 = p1.y - p0.y dx2 = p2.x - p0.x; dy2 = p2.y - p0.y return +1 if dx1*dy2 > dy1*dx2 return -1 if dx1*dy2 < dy1*dx2 return -1 if (dx1*dx2 < 0) || (dy1*dy2 < 0) #(1) return +1 if dx1*dx1+dy1*dy1 < dx2*dx2+dy2*dy2 #(2) return 0 #(3) end # 2線分が交差しているかを判定する def intersect(l1, l2) return ((ccw(l1.p1, l1.p2, l2.p1) * ccw(l1.p1, l1.p2, l2.p2)) <= 0 ) && ((ccw(l2.p1, l2.p2, l1.p1) * ccw(l2.p1, l2.p2, l1.p2)) <= 0 ) end end
このccwについて、なぜこういう関数にしてあるのか分析します。まだ途中です。
(1)-(3)はp0,p1,p2の3点が同一直線上にある場合(もしくは3点が同一の場合)の処理をしている箇所です。
以下、p0,p1,p2の3点が同一直線上にあり、p0とp1が異なる点であるとします。(つまり、線分p0p1が存在するとします)
- (1)により値が-1となるのは、p0から見てp2がp1と異なる側にある場合、つまりp2-p0-p1の順に並んでいる場合です。
- xについての条件とyについての条件がorで結びついているのは、直線がx軸に平行、またはy軸に平行な場合を考慮してのことです。
- (2)により値が+1となるのは、p0から見てp2がp1の向こう側にある場合、つまりp0-p1-p2の順に並んでいる場合です。
- (3)により値が0となるのは、p2が線分p0p1上にある場合、つまりp0-p2-p1の順に並んでいる場合です。これは両端点も含みます。
p0とp1が同じ点ならどうでしょうか。関数に dx1=dy1=0 を代入して追ってみますと、上の結果から類推できるような結果が出ます。実行した方が早いかもしれません。
このようなccwの仕様を、intersectにどう活かしているかは、時間がありませんので明日。