第20回 プログラミングについて 『角度はめぐる』

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