[Scala] 直行座標と極座標の相互変換の罠

直行座標上の点を回転させたいときに、極座標を使ってハマってしまいました。どんなとこにハマっちゃったかというと、θ方向に少しづつ加算していって点を回転させていった時のことです。θが180°に近づいたあたりで、そこからθの値がほとんど変わらなくなっちゃったんです。

なぜ！？

ここで少し情報を整理しましょう。

まず、これは３次元空間の座標系でのお話です。下の図のような感じです。原点からの距離をr、y軸とのなす角をθ、x軸とのなす角をφと定義しています。


もろもろの都合上、基本的には、座標系は直行座標系で表現することにしていました。Scalaでは、タプルで。

1	type Point = (Double, Double, Double)

なので、ある点を１回だけ回転させるためには、
直行座標 ⇒ 極座標 ⇒ 回転操作 ⇒ 直行座標
という操作を行っていて、複数回の回転はこの操作の連続適用していました。

次の関数は、直行座標系上の点をθ方向にd°だけ回転して移動して直行座標上の点を返す関数です。

1  2  3  4

val addTheta = ( p: Point, d: Double ) => {       val p1 = toPolar( p )       fromPolar( (p1._1, p1._2 + d.toRadians, p1._3 ) )     }

この関数をθが180°付近で適用すると、180°付近を行き来するようになるんです。

ちなみに、addTheta関数で使われているtoPolarとfromPolarは、それぞれの座標系へ変換する関数で次のように定義しました。

直行座標から極座標への変換

1  2  3  4  5  6

val toPolar: Point => Point = { case (x, y, z) =>        val r = sqrt( x*x + y*y + z*z )       val theta = acos( y/r )       val phi = if( z > 0 ) atan2( z, x ) else Pi*2 + atan2( z, x )       (r, theta, phi)     }

極座標系から直行座標への変換

1  2  3

val fromPolar: Point => Point = { case (r, theta, phi) =>         ( r*sin( theta )*cos( phi ), r*cos( theta ), r*sin( theta )*sin( phi ) )     }

さて、上記のaddTheta関数を使って、点 (1.0, -10.0, 1.0)(極座標表示では、(10.1, 172.0, 45))をθ方向に20°づつ4回だけ回転してみます。

1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23

// 度表示に変換する。     val toDegrees: Point => Point = { case (r, theta, phi) =>         ( r, theta.toDegrees, phi.toDegrees )     }     val p0 = (1.0, -10.0, 1.0)     println( toDegrees( toPolar( p0 ) ) )     // １度目の回転     val p1 = addTheta( p0, 20 )     println( toDegrees( toPolar( p1 ) ) )     // 2度目の回転     val p2 = addTheta( p1, 20 )     println( toDegrees( toPolar( p2 ) ) )     // 3度目の回転     val p3 = addTheta( p2, 20 )     println( toDegrees( toPolar( p3 ) ) )     // 4度目の回転     val p4 = addTheta( p3, 20 )     println( toDegrees( toPolar( p4 ) ) )

20°の回転を４回適用したので、172.0°から80°回転されて、θは252.0°(つまり、0 < θ < 180なので108°)になっているはず！　と思っていたんですが、結果は、

(10.099504938362077,171.95053302447164,45.0) (10.099504938362077,168.04946755734255,225.0) (10.099504938362077,171.95053302447164,44.99999999999999) (10.099504938362077,168.04946755734255,225.0) (10.099504938362077,171.95053302447164,44.99999999999999)

となりました。うーむ、171°と168°を行き来してますね。なんでこんなことが起きたんだー！

addTheta関数をもう１度見てみましょう。
さっきも言ったように、「直行座標 ⇒ 極座標 ⇒ 回転操作 ⇒ 直行座標」って操作をしているんですね。具体的にどのように値が変化しているのかを見るために、前述の点 (1.0, -10.0, 1.0)の値の変化を追ってみます。まず、点 (1.0, -10.0, 1.0)を極座標に変換するとだいたい(10.1, 172.0, 45)になります。これをθ方向に20°回転させると、172+20=192になり、この点は極座標上で(10.1, 192.0, 45)という位置に移動しました。そして、この点を直行座標に戻すと、(-1.478,-9.880,-1.478)になりました。ここまでは、問題ないでしょう。予想通り。では、この点(-1.478,-9.880,-1.478)を極座標表示で確認してみましょう。toPolar関数を適用しますね。きっと(10.1, 192.0, 45)に変換されることでしょう。ところが期待に反して、(10.099,168.049,225.0)という計算結果になっちゃってます。toPolarにバグがあるのでしょうか？　でも、ちょっと待ってください。もうちょっと結果をよく見てみましょう。rが10.1なのは予想通り。θとφは、予想と違いますが、θは、192も168も180°を挟んでともに12°だけ、φは、45も225もちょうど180°だけずれているではありませんか！つまり、これはどういうことかというと、(10.1, 192.0, 45)も(10.099,168.049,225.0)も数字こそ違いますが示している場所は同じだったんですね。紛らわしいよ～。じゃあ、どうしてこんな紛らわしいことになっちゃったの？？という疑問に答えるために、toPolar関数の定義を見てみましょう。その中で、直行座標からθの値を得るときに、acosを使ってますね。acosは、0°～180°までしか返しません。この定義のためにθは180°より大きい値は扱えないんです。その代わり、φが0～360°まで取れるので、θが180°を超えた分は、φを180°回転させて表現していたんですね。ということで、ここまでの計算過程は、問題なさそうですね。次のステップに進みましょう。さらにこの点(-1.478,-9.880,-1.478)をθ方向に20°回転させましょう。先ほどと同じようにaddTheta関数をこの点に適用します。同じことの繰り返しですね。さっさと極座標に変換してθに20°だけ加算しちゃいましょう。この点(-1.478,-9.880,-1.478)の極座標のθは、168°だから、これに20°足して178°にして・・・ん？なんかおかしいぞ。ちょっと整理しよう。最初に20°回転させて、さらに20°回転させたんだから、計40°足されて212°になり、θは180°までで表現しなければいけないので、142°になってないとおかしいはずですけど～！なるほど。いくら回転させても180°付近を往来する原因は、ここだったのか。θは180°以内で表さなければいけないというtoPolarの定義によるものだったのですね。そして、いちいち極座標に変換して、また直行座標に変換するというaddThetaがいけなかったんですね。172°+20°=168°　⇒　168°+20°=172°　⇒　172°+20°=168°・・・となっちゃうんですか。。。ということで、θを少しづつ増加させて360°回転するのは、このやり方では、不可能！

では、どうすればいいのか？

いちいち直行座標に変換なんてさせずに、ずっと極座標の値を保持していくのが一番簡単なんでしょうね。直行座標が必要になったときだけ、極座標から変換して使うようにすればいいんでしょう。
または、一貫して直行座標にこだわるのであれば、回転行列を使って直行座標を変換するのもいいのかもしれません。