第20回 プログラミングについて 『角度はめぐる』
『角度はめぐる』
図形を扱うプログラムでは避けられないものに角度があります。これが実にいやらしいもので一周するとまた０度になってしまうのです。何回まわしても一周すると０度になるのです。
安易にプログラムを書いていると、プログラムの中では３６０度以上の角度になったり、マイナスになったりして思わぬバグになってしまうこともありますが、もっとも厄介なのが、計算した値が０度か否かを判断するときです。
:
:
if(angle>= -0.01 && angle<=0.01){
:
:
}
としたときに、どうしてもこの ｉｆ文に引っかかってくれないことがあります。angleの値を調べてみると 359.99 度だったり、360.01 度だったり、確かに０度の角度なのですが一周以上多かったり、少なかったりするのです。そこで、
:
:
while(angle>=360.0)angle-=360.0;
while(angle< 0.0)angle+=360.0;
if(angle<=0.01 || angle>=360.0-0.01){
:
:
}
:
:
などとするのです。これは、角度を０度以上３６０度未満の範囲に直してその両端であるか否かで判断する方法です。
通常、プログラム内で使用する角度の単位は『度』ではなく、『ラジアン』を使用します。『度』は一周を３６０とするのに対し、『ラジアン』は一周を２πで表現します。『度』という単位は人間の長い歴史のなかで慣習的に使われてきたもので、極論すれば数学的にはあまり意味を持っているものではありません。それに対し『ラジアン』は円と切っても切れない深い縁をもっているπ（円周率）を基準としているのです。もっと角度を扱いやすくしようと『グラジアン』なる単位も考えられました（一周を４００とする）が、結局は普及はしなかったようです。
ちなみに、度とラジアンの間の換算は、度＝３６０／π＊ラジアン です。
さて実際にプログラムのなかでπを使うとしたとき、
angle=π;
などと書いてはコンパイラに怒られてしまいます。
#define PAI 3.14
などとして定数を定義して使いますが、 3.14 では精度が足りないのでプログラムは正常には動作してくれません。実際に使うにはもっと精度が必要です。
#include <math.h>
main()
{
printf("%30.20f\n",acos(-1.0));
}
というプログラムを実行するとπの値が表示されます。
3.14159265358979300000
となりますので、この値を使うのが簡単です。acos(-1.0)をそのまま使っても問題はありませんが、定数を求めるために acos の計算を実行時行なうことになることと、プログラムが理解しにくくなりますので、定数として定義した方がいいでしょう。
それでは実際に直線の角度を考えてみましょう。
直線は、
ＡＹ＝ＢＸ＋Ｃ
と数学的には定義されますが話を簡単にするために、
Ｙ＝ＡＸ＋Ｂ
として話をすすめて行きます。
今、２つの点Ｐ１（ｘ１、ｙ１）とＰ２（ｘ２、ｙ２）が与えられたとき係数ＡとＢ
は、
ｙ１＝Ａｘ１＋Ｂ －－－（１）
ｙ２＝Ａｘ２＋Ｂ －－－（２）
の連立方程式から求められます。
式（２）から式（１）を引くと、
（ｙ２－ｙ１）＝（ｘ２－ｘ１）Ａ
∴Ａ＝（ｙ２－ｙ１）／（ｘ２－ｘ１） －－－（３）
式（３）を（１）に代入すると、
ｙ１＝ｘ１（ｙ２－ｙ１）／（ｘ２－ｘ１）＋Ｂ
∴Ｂ＝ｙ１－ｘ１（ｙ２－ｙ１）／（ｘ２－ｘ１）
＝（ｘ２ｙ１－ｘ１ｙ１－ｘ１ｙ２＋ｘ１ｙ１）／（ｘ２－ｘ１）
＝（ｘ２ｙ１－ｘ１ｙ２）／（ｘ２－ｘ１）
となります。今回は角度が主題ですのでＡの値を考えてみましょう。
実は、
Ａ＝（ｙ２－ｙ１）／（ｘ２－ｘ１）
の式だけではプログラムで使いたい角度を求めるには情報が不足なのです。確かに数学的には直線の角度なのですが、この式では４５度のときも２２５度のときもＡの値が同じになってしまうのです。このＡはあくまでもＸの変化量にたいするＹの変化量の絶対値に過ぎません。
まずＡの式をよく見てみると、Ａは三角関数で言えばタンジェントのことです。いま角度をＴとすると、
Ａ＝ｔａｎ（Ｔ）＝（ｙ２－ｙ１）／（ｘ２－ｘ１）
∴Ｔ＝ａｔａｎ（（ｙ２－ｙ１）／（ｘ２－ｘ１））
となって角度Ｔがとりあえず求まります。ところが先程も言いましたようにこれではプログラムで使えるような角度ではありません。ａｔａｎは－９０度～９０度までの角度しか答えてくれないのです。そこでスパイスが必要になります。
Ｐ２とＰ１の差の成分をそれぞれ次のようにします。
ＤＸ＝Ｘ２－Ｘ１
ＤＹ＝Ｙ２－Ｙ１
この式で，
ＤＸ｜＋｜－｜－｜＋｜
ＤＹ｜＋｜＋｜－｜－｜
－－＋－＋－＋－＋－｜
象限｜１｜２｜３｜４｜
ＤＸとＤＹの符号によってＰ１を原点とするとＰ２の位置が第何象限にあるのかが分かります。この位置関係が分かっているとａｔａｎで求める範囲は０度から９０度までで十分です。
次に１つだけ大きな問題が残っています。ＤＸの値が０のとき（厳密には０に限りなく近いとき）の処理をどうするかです。数学の教科書で見た記憶があると思いますが、角度に対するタンジェント関数は、９０度付近と２７０度付近で正または負の無限大の値になっていて非連続な曲線になっているのです。
ｔａｎ（Ｔ）＝（ｙ２－ｙ１）／（ｘ２－ｘ１）
の式からも分かるようにＤＸの値、すなわち分母の（Ｘ２－Ｘ１）が９０度と２７０度付近でほとんど０に近くなり、９０度と２７０度で０になります。この処理を行わないと直線が垂直のときに浮動少数点数のオーバーフローが発生してしまいます。
これに対処する一般的な方法はＤＸの値がある範囲以内になれば０として特別な処理をすることです。たとえば０．０００１以内を０として特別な処理を行うとすれば、
if(dx>= -0.0001 && dx<=0.0001){
:
:
}
else{
:
:
}
などと記述して浮動少数点数のオーバーフローは避けられることになります。
がしかし、この方法は一見では正しい方法のように見えるのですが実はこの方法では特殊な場合に正常な角度を求められないことがあるのです。その特殊な場合というのは直線の長さが０．００１程度の短いものになると８０度程度でも９０度と判断してしまうのです。当然ながらもっと短い直線になればどんなときでも９０度と判断してしまいます。９０度と判断する範囲を小さくしても結局はイタチごっことなった上、長い直線に対しては精度とか誤差の点で感心できません。
実際、私も上記の方法で角度を求める関数を作って使っていましたが、長い間使っていたプログラムに突然バグが発生しその原因は短い線の角度だったのです。垂直と判断する範囲をいろいろと変えて試してみましたが、この方法では根本的な解決にはならないことが分かりました。
『垂直をどうやって処理するか』。さんざん頭をひねって思い付いたのは、『垂直を処理しないようにすればいい』ということだったのです。このように言ってしまうとこれまでの話とは矛盾するようですが、実に巧妙な抜け道があるのです。ａｔａｎ関数に垂直近辺の計算をさせないようにするのです。先ほども説明しましたが、ａｔａｎは次のように考えられます。
ｔａｎ（Ｔ）＝（ｙ２－ｙ１）／（ｘ２－ｘ１）＝ＤＹ／ＤＸ
ＤＹの絶対値がＤＸの絶対値よりも小さいときには０度以上４５度未満、等しいときには４５度、大きいときには４５度よりも大きく９０度以下となります。この４５度よりも大きな角度になったとき（ＤＹの絶対値がＤＸの絶対値よりも大きくなったとき）に直線を４５度の角度と線対称に反転させるのです。すなわち、
ｔａｎ（Ｔ）＝（ｘ２－ｘ１）／（ｙ２－ｙ１）＝ＤＸ／ＤＹ
として分子と分母を入れ替えるのです。こうすることで精度の低下する垂直近辺の処理はもとより、分子以上の大きさの分母での除算になるので誤差の増大も押さえられるようになり一石二鳥となるのです。結局のところａｔａｎ関数には０度から４５度までの計算しかさせないことになるのです。
このようにして作ったのが次のプログラムです。
#include <math.h>
double angle(x1,y1,x2,y2)
double x1,y1,x2,y2;
{
double dx,dy,t,a;
int area;
dx=x2-x1;
dy=y2-y1;
if(dx>=0.0){
if(dy>=0.0){ area=0; } /* 第１象限 */
else { area=3; t=dx; dx= -dy; dy= t; } /* 第４象限 */
}
else{
if(dy>=0.0){ area=1; t=dx; dx= dy; dy= - t; } /* 第２象限 */
else { area=2; dx= -dx; dy= -dy; } /* 第３象限 */
}
if(dy>dx)a=PAI/2.0-atan(dx/dy); /* ４５度以上 */
else a= atan(dy/dx); /* ４５度以下 */
return a+area*(PAI/2.0);
}
好みによってはいろいろと表現の方法はあるかと思いますが、アルゴリズムはこれまで説明してきた内容を利用しています。この関数の内容からも分かるように、垂直の処理はどこにもなく、ａｔａｎで求める角度は０～４５度（プログラムでの角度の単位はラジアンです）までの範囲でしかありません。
２つの直線の成す角度を求めるとき、それぞれの直線の角度を求めてからその差を計算するはプログラムがめんどうになります。たとえば直線の１つの角度が３５０度で、もう１つのものが５度であったとき人間は直感的に１０度と判断できますが、これをプログラムで行うのは複雑になってしまいます。
以前にも複素数の話のときにも触れましたが、２つの直線をベクトルとして扱いベクトルの割り算をした結果が２つの直線の角度の差の成分となりますので、この成分から単純に角度を求めることができます。
角度はぐるぐる回しても一周すればまた０度から始まります。タンジェント関数は非連続です。単位は度もラジアンもあります。比であらわすこともあります。プログラムのなかで角度をハンドリングしようとするとたくさんの障害を克服しなければなりません。図形を扱うプログラムを初めて作る人たちの最初の壁になるのが以外とこの角度であることが多いようです。でもここでくじけていては、その後にあるたくさんの壁と出会えません。今回は角度についての一面の話でしたが、何かの参考になれば幸いです。
それではまた次回。
