問題発生
値は殆ど同じではあるけれども、角度が少なめに出る計算式と、多めに出る計算式があるとします。この場合、計算精度を上げるために、これらの２つの値を平均して利用しようと考えるのは自然な流れです。しかし、このような場面で、想定していなかった問題に直面してしまいました。
例えば、それぞれの角度を atan2 関数を使って求めている場合、それらは -$\pi$ ～ $\pi\hspace{0.2em}$ ( -180°～180° ) の範囲の値になります。２つの値が 5°, -5° とか 90°, 100° のときには、これらを単純に平均するだけで、0° とか 95° が求まるので問題はありません。しかし、175°, -175° のような場合には、期待する 180° ではなく、全く方向違いの 0° になってしまいます。これを受けて処理を続ける後続のプログラムは、まともな答を出せるはずがありません。
２つの角度差がいつでも非常に小さいようなケースでは、この不運な条件に出会える確率はかなり低くなります。通り一遍の動作確認だけで、めでたく完成！などと安易な妥協をしていると、あるとき突然、栄えある衆人環視の大舞台でとんでもない答えを曝してしまいます。ということで、プログラム作成者は生き地獄を味わうことになります。
対策（正攻法）
この問題に対処するために、「角度 平均」でググってみると、すぐに次の投稿が見つかります。
角度の平均を求める方法 - Qiita
角度の平均や分散を複素数を用いて求める - Qiita
これらの投稿が紹介している方法の要旨は次のとおりです。
ベクトルや複素数に変換してから平均し、それを角度に戻す。すなわち、平均すべき角度を $\theta_1$, $\theta_2$, $\cdots$, $\theta_n$ 、求める平均値を $\theta$ とすると、
$z_k=\cos\theta_k+i\sin\theta_k$ ($k=1,2,\cdots,n$)
$\displaystyle z=\frac{1}{n}\mathop{\small\sum}_{k=1}^n z_k$
$\displaystyle\theta=\tan^{-1}\frac{\Im(z)}{\Re(z)}$
ただし、$\tan^{-1}$ は atan ではなく、atan2
しかし、殆どのケースでは単純平均でこと足り、アセンブラでも数ステップだけで書けるような処理なのに、特殊ケースの救済だけのために、複雑な関数まで持ち出して対応するしか方法はないのでしょうか？ すこし考えてしまいます。
対策（簡易法-２値）
いま直面しているケースでは、平均化される角度は２つだけで、$n\hspace{-0.2em}=\hspace{-0.2em}2$ です。もっと簡単な方法があるはずです。そこで、問題の本質に立ち返って考え直してみます。
尖端を原点に固定されてはいても、任意に傾くことができるVの字を想定します。問題の根源は、このVの字の内側の中心線の角度を求めたいのに、ある特殊条件下では外側の中心線の角度が求まってしまうということです。この特殊条件というのは、Vの字が、自身の内側に極座標の±180°の線を挟み込んでしまうように傾いた状態です。ですから、この状態のときに限り、単純平均して得られた答とは真逆の角度を選ぶようにすれば、すべては解決します。
すなわち、
$|\theta_1-\theta_2|\leq\pi$（Vの字の内側に±180°の線が無い条件）のとき、
$\theta=(\theta_1+\theta_2)/2$
それ以外のときで、
$\theta$ の範囲をあくまで $-\pi\lt\theta\leq\pi$ に限定したいとき、
$\theta_1+\theta_2\gt 0$ のとき、
$\theta=(\theta_1+\theta_2)/2-\pi$
$\theta_1+\theta_2\leq 0$ のとき、
$\theta=(\theta_1+\theta_2)/2+\pi$
$\theta$ に $\pm 2\pi N$ のバイアスが掛かるのを許容するとき、
$\theta=(\theta_1+\theta_2)/2+\pi$
となります。if 文を使った多行のプログラムを避けたい場合には、可読性は下がりますが下記でもいけます（ただし、th: $\theta$, th1: $\theta_1$, th2: $\theta_2$）。
th=(th1+th2)/2+(abs(th1-th2)>pi)*(((th1+th2)<0)*2-1)*pi;
sign(th1+th2) は、±1 だけでなく 0 を返すこともあるのでここでは使えません。((th1+th2)<0)*2-1 はその代用です。
得られた $\theta$ を後続のプログラムでそのまま使うのではなく、三角関数の引数として使うだけであれば、さらに次のように簡略化できます。（$\theta=\theta\pm2\pi N$ とみなせるゆえ）
th=(th1+th2)/2+(abs(th1-th2)>pi)*pi;
以上では、$\theta_1$ と $\theta_2$ の対が１組しかないことを前提として記述しています。しかし、実際には複数の対について計算しなければならないこともあります。この計算を for ループで処理するのであれば、もうこれ以上何も考える必要はありませんが、MATLABを使うのであれば一気に処理してしまいたいものです。その場合は次のようになります。（ただし、th1(n): $\theta_{1n}$, th2(n): $\theta_{2n}$, th(n): $\theta_{1n}$と$\theta_{2n}$の平均値, n: 対番号）
$\theta$ に $\pm 2\pi N$ のバイアスが掛かるのを許容するとき、
th=(th1+th2)/2;
nn=find(abs(th1-th2)>pi);
th(nn)=th(nn)+pi;
$\theta$ の範囲を $-\pi\lt\theta\leq\pi$ に限定したいときは、上記に続けて次の２行を追加します。
nn=find(th>pi);
th(nn)=th(nn)-2*pi;
図１に無対策の場合の結果、図２に対策済の結果を示します。簡易対策法も正攻法も全く同一の結果を示しており、簡易対策法の正当性が確認できました。処理時間は、なぜか実行回ごとに大幅に異なるので具体値は示しませんが、当然、簡易対策法の方が短くなっています。
対策（簡易法-多値）
以上、平均化される角度が２つだけの場合について、±180°段差問題を回避する簡易方法を考えてみました。さらに発展して、平均化する角度が３つ以上の場合ならどうなのか？という疑問が湧いてきます。
ここで述べる方法を一言でいえば、±180°線を、角度データが集中した領域を避け、データが空いた安全な位置にまで一時的に遠ざけて、平均化処理が済んでから元に戻そうとするものです。具体的な手順はややこしいので割愛しますが、気になる方は添付プログラム中のコメントを参考にしてください。
この方法では、平均化すべき角度データ群がすべて、広がり幅が180deg未満の部分円の範囲に収まっている場合には、唯一の平均値を求めることができます。しかし、これ以外の場合には、データが空いた位置が複数個所に分散されてしまうケースもあるので、どこを選ぶかによって結果が異なってきます。
図３に無対策の場合の結果、図４に対策済の結果を示します。こちらは、簡易対策法と正攻法では結果に多少の差がありますが、無対策のときのような段差問題は排除できています。両対策法による結果の違いはあっても、簡易対策法が間違っているという訳ではなく、平均というものをどのように定義するのか という根本的な問題によるものです。ただし、前記のように、平均化すべきデータが２値だけの場合には、原理的に両対策法の計算結果に差は生じません。
なお、処理時間は、簡易対策法の方が短くなることを期待していましたが、何と、正攻法の方に軍配が上がってしまいました。for ループや分岐は一切使わないように努力したのですが、これでは残念ながら「簡易」の意味が見出せません。
まとめ
２つの角度の平均を求めるだけであれば、正攻法よりも簡易対策法の方が良さそうなことは分かりました。しかし、平均化するデータ数が多い場合には、正攻法を使っておくのが無難なようです。
評価プログラム
% angle_average01.m
%
% 【プログラム概要】
%
% ■■ 角度の平均値の計算（±180度線上の段差障害への対策） ■■
%
% ２つの角度データ（その範囲は-180deg＜θ≦180degとする）があるとき、
% これらの平均値を単純な相加平均として求めるだけでは、特定の条件下で
% 次のような問題が生じる。
%
% -5°と 5°の平均 ⇒ 0°（正解！）
% 90°と 100°の平均 ⇒ 95°（正解！）
% -175°と 175°の平均 ⇒ 0°
% （なんじゃそりゃ！、正解は 180°じゃ！）
%
% この問題を回避する、平均値計算のスマートな正攻法は、
% ２つの角度をth1,th2、平均値をthとしたとき、
% z1=cos(th1)+i*sin(th1); z2=cos(th2)+i*sin(th2); z=(z1+z2)/2;
% th=atan2(imag(z),real(z));
% である。
%
% しかし、これでは複雑すぎる。もっと簡単な計算法を探してみる。
%
% ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
%
% 尖端を原点に拘束されながらも、好き勝手に傾き回るVの字を想定する。
% Vを構成する２本の線分のそれぞれの角度を th1,th2 とする。
% 問題の根源は、このVの字の内側の中心線の角度を求めたいにも拘わらず、
% ある特定の条件下では、外側の中心線の角度が求まってしまうことだ。
% この特定の条件というのは、Vの字が、自身の内側に極座標の±180°の線
% を挟み込むように傾いた状態である。
%
% ゆえに、この状態のときに限り、単純平均して得られた答とは真逆の角度を
% 選ぶようにすれば全ては解決する。
%
% 具体的には、
%
% abs(th1-th2)<=piのとき（Vの字の内側に±180°の線が無いとき）、
% （Vの字には、尖り過ぎたもの や 太り過ぎたもの など色々あるが、ここ
% では、その開き角度は 0°超から 180°未満まで、幅広く受け入れる）
% th=(th1+th2)/2
% それ以外のときで、
% th を、あくまで -pi<th<=pi の範囲内だけで表現したいとき、
% th1+th2>=0 のとき、
% th=(th1+th2)/2-pi
% th1+th2<0$ のとき、
% th=(th1+th2)/2+pi
% th に ±2*pi*N のバイアスがかかっても構わないとき、
% （後続処理で、thを三角関数の引数としてだけ使用するようなとき）
% th=(th1+th2)/2+pi
%
% 場合分けをせず、１行だけで処理したいときには、
%
% th の範囲を -pi<th<=pi に限定したいとき、
% th=(th1+th2)/2+(abs(th1-th2)>pi)*(((th1+th2)<0)*2-1)*pi;
%
% th に ±2*pi*N のバイアスが掛かるのを許容するとき、
% th=(th1+th2)/2+(abs(th1-th2)>pi)*pi;
%
% th1 と th2 をスカラーとしてではなく、複数の角度対を含む配列とみなし、
% これらの角度対を一括処理したいとき、
%
% th に ±2*pi*N のバイアスが掛かるのを許容するとき、
%
% th=(th1+th2)/2;
% nn=find(abs(th1-th2)>pi);
% th(nn)=th(nn)+pi;
%
% th の範囲を -pi<th<=pi に限定したいとき、
% 上記の３行に加えて、その後に次の２行を追加。
%
% nn=find(th>pi);
% th(nn)=th(nn)-2*pi;
%
% これらによる計算結果とスマートな正攻法による計算結果を比較し、両者が
% 等しいことを確認する。
%
% ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
%
% 次に、平均する対象の角度が２つだけではなく、３つ以上ある場合に
% ついても、簡易対策法を考えてみる。
%
% しかし、あまり汎用性に拘り過ぎると難しくなる。ここでは割り切って、
% 平均化すべき角度データ群はすべて、広がり幅が180deg未満の部分円の範
% 囲に収まっているものと仮定する。
% （運さえ良ければ、もっと広がっていても大丈夫）
%
% 以下にその手順を示すが、たぶん、ややこしくて分かり難いと思う。
% 一言でいえば、±180°線を、角度データが集中した領域を避けて、デー
% タが空いた安全な位置にまで一時的に遠ざけ、平均化処理が済んでから元
% に戻そうというもの。
%
% ■ まず、データが空いた位置を探す。
%
% 手順①
% 平均化すべきデータ群で構成された配列を、昇順に並べ変えて、配列 [A]
% とする。
%
% 手順②
% 配列 [A] の後に配列 [A]+2*pi を継ぎ足し、２倍の長さをもつ配列 [B]
% を作る。（最も広い[データの空白部分]が、±180°線が存在する領域に
% あるときには、段差を除くように継ぎ足して考えないと、その空白の広さ
% を認識しにくい）
%
% 手順③
% 配列 [B] を昇順に辿り、要素間の差異が最大になる部分を見つける。
% （最大になる部分が複数あるときには、どちらを採用するかによって平均
% 値の結果が変わってしまうが、ここでは割り切って、最初に見つけた方
% を採用する ）
%
% 手順④
% その部分の直上と直下にある配列要素の角度の中間値を thc とする。
% thc の値が pi を超えるときには、thc-2*pi を新たな thc とする。
%
% ■ ±180°線をデータが空いた所まで回転させ、データ値を読み替える。
%
% 手順⑤
% 今までの±180°線を、この thc の部分に一致させるように、一時的に、
% 極座標を反時計方向に回転させる（回転量は thc+pi になるはず）。
% 結果として、[A] の各要素の値は [A]-(thc+pi) に変換される。これを
% [C] とする。[C] の一部の要素の値は -pi 以下となってしまうので、こ
% れらについてだけ +2*pi の補正を加える。
% 補正済の [C]全体を [D] とする。
%
% ■ この状態で平均化処理をしたあと、±180°線を元の位置に戻す。
%
% 手順⑥
% [D]の全要素の相加平均値を求め th とする。
%
% 手順⑦
% ±180°線を元に戻すと、th は th+(thc+pi) に変わる。これを新たな th
% とする。th が pi を超える場合には、th-2*pi を新たな th とする。
% この th が、求めようとしている最終結果になる。
%
% この結果も、スマートな正攻法による計算結果と比較する。ただし、両者は
% 完全には等しくならない。しかし、簡易対策法が間違っているという訳で
% はなく、平均というものをどのように定義するのか という根本的な問題
% による差なのである。
clear
close all
disp(' ');
disp(' ');
% ■■■ ２値だけの平均の場合
% テスト用のデータの作成。
% [60度の広がり角を持つV字]の中心線の角度th0 を -180度から180度まで、
% 1度ステップで変化させたときの、V字の両線分の角度データを作る。
thd=30*pi/180; % V字の開角の半分。
th0=linspace(-pi,pi,361); % V字の中心線の角度（1°ステップ）。
th1=th0+thd; % V字の左の線(正立の状態で見て)の角度。
th1(th1>pi)=-2*pi+th1(th1>pi); % th1 には180度を超える要素が出てくる
% ので、それらの要素は 2*pi を引い
% てどの th1 も -pi<th1<=pi となる
% ように条件を合わせる。
th2=th0-thd; % V字の右の線(正立の状態で見て)の角度。
th2(th2<-pi)=2*pi+th2(th2<-pi); % th1 に準じて -pi<th2<=pi に抑える。
% ■ 無対策の場合の計算
th=(th1+th2)/2;
th_ng=th; % 間違った答え
% ■ 簡易対策法による計算
tic; % 計時タイマー スタート
th=(th1+th2)/2; % 【プログラム概要】の中の式そのまま。
nn=find(abs(th1-th2)>pi); % 最終の th が -pi<th<=pi となるよ
th(nn)=th(nn)+pi; % うに考慮した汎用性が高い方の式。
nn=find(th>pi);
th(nn)=th(nn)-2*pi;
tt=toc; % 処理時間
th_simp=th; % 簡易対策法による答え
disp(['簡易対策法による２値の平均角度：処理時間 ' num2str(tt) ' 秒']);
% ■ 正攻法による計算
tic;
z1=cos(th1)+i*sin(th1); % 【プログラム概要】の中の式そのまま。
z2=cos(th2)+i*sin(th2);
z=(z1+z2)/2;
th=atan2(imag(z),real(z));
tt=toc;
th_form=th; % 正攻法による答え
disp(['正攻法による２値の平均角度：処理時間 ' num2str(tt) ' 秒']);
% ■ 結果の図示
figure(1)
plot(th0,th_ng,'k','LineWidth',6);
hold on
plot(th0,th1,th0,th2,'LineWidth',2);
xticks(linspace(-pi,pi,13));
yticks(linspace(-pi,pi,13));
xticklabels({'-180','-150','-120','-90','-60','-30','0', ...
'30','60','90','120','150','180',});
yticklabels({'-180','-150','-120','-90','-60','-30','0', ...
'30','60','90','120','150','180',});
xlim([-pi*1.1 pi*1.1]);
ylim([-pi*1.1 pi*1.1]);
grid on
title('２つの角度の平均値（段差障害は無対策）');
xlabel('入力データ群の基準線の向き[deg]')
ylabel('入力データ と それらの平均値[deg]')
legend('平均値(無対策)','入力１','入力２','Location','south');
text(-pi*1.34,pi*1.13,'図１','FontSize',20);
figure(2)
plot(th0,th_simp,'k','LineWidth',6);
hold on
plot(th0,th_form,'r','LineWidth',2);
plot(th0,th1,th0,th2,'LineWidth',2);
xticks(linspace(-pi,pi,13));
yticks(linspace(-pi,pi,13));
xticklabels({'-180','-150','-120','-90','-60','-30','0', ...
'30','60','90','120','150','180',});
yticklabels({'-180','-150','-120','-90','-60','-30','0', ...
'30','60','90','120','150','180',});
xlim([-pi*1.1 pi*1.1]);
ylim([-pi*1.1 pi*1.1]);
grid on
title('２つの角度の平均値（段差障害は対策済）');
xlabel('入力データ群の基準線の向き[deg]')
ylabel('入力データ と それらの平均値[deg]')
legend('平均値(簡易対策法)','平均値(正攻法)', ...
'入力１','入力２','Location','south');
text(-pi*1.34,pi*1.13,'図２','FontSize',20);
% ■■■ ３値以上を平均する場合
% テスト用のデータの作成。
% 平均対象の角度データが多すぎると、問題点があっても顕著な症状が
% 現れないことが考えられるので、３値だけの平均値で確認する。
%
% 実際の用途では、平均化する集団毎に、平均化対象のデータ数が異なるな
% ど、木目細かな工夫が必要になることも多いと思われる。ここでは、そ
% のような配慮は一切していないので、用途に応じた改良が必要である。
thdd=[45 20 -55]*pi/180; % 原点から出る３本の放射線の
% th0 から見た相対角度。
th1=th0+thdd(1);
th1(th1>pi)=-2*pi+th1(th1>pi); % th1 の pi 超の要素だけ -2*pi 補正。
th2=th0+thdd(2);
th2(th2>pi)=-2*pi+th2(th2>pi); % th2 の pi 超の要素だけ -2*pi 補正。
th3=th0+thdd(3);
th3(th3<=-pi)=2*pi+th3(th3<=-pi); % th3 の -pi 以下の要素だけ、
% +2*pi 補正。
% ■ 無対策の場合の計算
th=(th1+th2+th3)/3;
th_ng=th; % 間違った答え
% ■ 簡易対策法による計算
tic;
% 入力データ
Thi=[th1;th2;th3]; % [th1_1 th1_2 th1_3 ... th1_n]
% Thi = [th2_1 th2_2 th2_3 ... th2_n]
% [th3_1 th3_2 th3_3 ... th3_n]
N=size(Thi,2); % 平均化が必要な[３値]の組の数。上式のnに相当。
% 手順①
ThA=sort(Thi); % 各列ごとに、昇順で並べ替え。
% 結果は【プログラム概要】の[A]に対応。
% 明記していなかったが、配列[A]は列ベクトル。
% 行列 ThA の各列がこれに相当。
% 手順②
ThB=[ThA; ThA+2*pi]; % バイアスをかけた配列を継ぎ足して長くする。
% 【プログラム概要】の[B]に対応。
% 手順③の準備
th_bot=-pi*ones(1,N); % 下限角度(-pi)
th_top=3*pi*ones(1,N); % 上限角度(+3*pi)
Th_large=[ThB;th_top]; % Th_large(k)=[th1_k ]
% [th2_k ]
% [th3_k ]
% k:列番号 [th1_k+2*pi]
% [th2_k+2*pi]
% [th3_k+2*pi]
% [3*pi ]
Th_small=[th_bot;ThB]; % Th_small(k)=[-pi ]
% [th1_k ]
% [th2_k ]
% [th3_k ]
% [th1_k+2*pi]
% [th2_k+2*pi]
% [th3_k+2*pi]
% 各角度データ間の間隔
Th_space=Th_large-Th_small; % Th_space(k)=[th1_k - (-pi) ]
% [th2_k - th1_k ]
% [th3_k - th2_k ]
% [th1_k+2*pi - th2_k]
% [th2_k - th1_k ]
% [th3_k - th2_k ]
% [pi - th3_k ]
% 手順③の本番
% 各列内での最大の間隔値 と その行番号
[max_space,id_max]=max(Th_space); % max_space は使しないが・・・。
% 手順④
id_thx=([1:N]-1)*7+id_max;
Thc=(Th_large(id_thx)+Th_small(id_thx))/2;
% 最も広い空白領域の中心の角度。
Thc(Thc>pi)=Thc(Thc>pi)-2*pi; % pi を超えるものは -2*pi を足す。
% 手順⑤
ThC=ThA-(repmat(Thc,3,1)+pi); % 【プログラム概要】の[C]に相当。
% 紛らわしいが ThC≠Thc ゆえ注意。
% 【プログラム概要】との対応は、
% Thc ⇔ thc, ThC ⇔ [C]
ThD=ThC;
ThD(ThD<=-pi)=ThD(ThD<=-pi)+2*pi; % 【プログラム概要】の[D]に相当。
% 手順⑥
Th=mean(ThD); % 平均値 th を並べた行ベクトル。
% 手順⑦
Th=Th+(Thc+pi);
Th(Th>pi)=Th(Th>pi)-2*pi; % 簡易対策法による最終結果。
tt=toc; % 処理時間
th_simp=Th; % 簡易対策法による答え
disp(['簡易対策法による３値の平均角度：処理時間 ' num2str(tt) ' 秒']);
% ■ 正攻法による計算
tic;
z1=cos(th1)+i*sin(th1); % 【プログラム概要】の中の式そのまま。
z2=cos(th2)+i*sin(th2);
z3=cos(th3)+i*sin(th3);
z=(z1+z2+z3)/3;
th=atan2(imag(z),real(z));
tt=toc;
th_form=th; % 正攻法による答え
disp(['正攻法による３値の平均角度：処理時間 ' num2str(tt) ' 秒']);
% ■ 結果の図示
figure(3)
plot(th0,th_ng,'k','LineWidth',6);
hold on
plot(th0,th1,th0,th2,th0,th3,'LineWidth',2);
xticks(linspace(-pi,pi,13));
yticks(linspace(-pi,pi,13));
xticklabels({'-180','-150','-120','-90','-60','-30','0', ...
'30','60','90','120','150','180',});
yticklabels({'-180','-150','-120','-90','-60','-30','0', ...
'30','60','90','120','150','180',});
xlim([-pi*1.1 pi*1.1]);
ylim([-pi*1.1 pi*1.1]);
grid on
title('３つの角度の平均値（段差障害は無対策）');
xlabel('入力データ群の基準線の向き[deg]')
ylabel('入力データ と それらの平均値[deg]')
legend('平均値(無対策)','入力１','入力２','入力３','Location','south');
text(-pi*1.34,pi*1.13,'図３','FontSize',20);
figure(4)
plot(th0,th_simp,'k','LineWidth',6);
hold on
plot(th0,th_form,'r','LineWidth',2);
plot(th0,th1,th0,th2,th0,th3,'LineWidth',2);
xticks(linspace(-pi,pi,13));
yticks(linspace(-pi,pi,13));
xticklabels({'-180','-150','-120','-90','-60','-30','0', ...
'30','60','90','120','150','180',});
yticklabels({'-180','-150','-120','-90','-60','-30','0', ...
'30','60','90','120','150','180',});
xlim([-pi*1.1 pi*1.1]);
ylim([-pi*1.1 pi*1.1]);
grid on
title('３つの角度の平均値（段差障害は対策済）');
xlabel('入力データ群の基準線の向き[deg]')
ylabel('入力データ と それらの平均値[deg]')
legend('平均値(簡易対策法)','平均値(正攻法)', ...
'入力１','入力２','入力３','Location','south');
text(-pi*1.34,pi*1.13,'図４','FontSize',20);
disp(' ');
disp(' ');