１．概要
MATLABには標記の機能の組み込みコマンドは無いようなので、必要に迫られてユーザー関数を作ってみました。実行結果は図１のとおりです。青い曲線群と赤い曲線群の交点が緑色と黒色の〇印です。各曲線群の実体は、複数の子曲線を、NaN折点([x y]=[NaN NaN])の継手で繋げて１本にまとめたものです。なお、「曲線の多数の交点」と書いてはおりますが、もっと簡単な２本の直線どうしの一つの交点の検出にももちろん使えます。
私自身の当面の用途にはほとんど問題なく使えていますので、お役に立てればと思いプログラムを添付しておきました。
図１
２．使い方：
MATLABのコマンドラインから、「help intersection_of_two_curves」と打ち込めば、次のように表示されます。
manual.txt
intersection_of_two_curves.m
平面上の２本の曲線(折線)の交点座標を返すユーザー関数。
（動作確認バージョン：MATLAB Home R2019a）
それぞれの曲線は、複数の子曲線を１本にまとめたものでもよい。
その場合、子曲線間には、区切り折点としてNaNを挿入すること。
出力されるのは、第１曲線と第２曲線どうしの交点のみ。同一曲線
内の子曲線どうしの交叉は無視される。
[x,y,uc_code]=intersection_of_two_curves(x1,y1,x2,y2)
【入力】
x1: 第１曲線の点列のx座標（行ベクトル）
y1: 同じく、y座標（行ベクトル）
（普通は、折点数の少ない方を第１曲線にすると処理が速い）
x2,y2: 第２曲線の、上記同様のx.y座標（行ベクトル）
【出力】
x: 交点のx座標（交点の個数に相当する長さの行ベクトル）
y: 同じく、y座標
uc_code: 交点の判定が不確実だった場合、その情報（5行×[件数]列）。
1行目： 不確実だった第１曲線の線分の若番側の折点番号
2行目： 同じく、 第２曲線の線分の若番側の折点番号
3,4行目： それらの各線分のどちら側の端が問題だったのか？
0: 両端とも問題なし、
1: 若番側の端、 2: 後番側の端、
3: 端の問題ではなく、両線分の重なり合い。
5行目： 疑いながらも、交点座標を出力済のときはその交点番号。
交点座標の出力を見合わせたときは0。
３．残された課題
そこそこには動作しますが、まだ改善の余地はあります。
３－１ 検出漏れや多重検出の可能性
双方の線分が次の関係になったとき、現状の処理だけでは、計算誤差によって、交点の検出漏れや多重検出を起こす心配があります（図２参照）。
・片方の曲線の折点の直上を、他方の曲線の線分が交叉
・双方の曲線の折点どうしが一致して交叉
・双方の曲線の線分どうしが重なりながら交叉 等々
図２
図１のサンプル曲線でも多重検出が発生しています。緑色の〇印は普通に検出された交点ですが、10個の黒色の〇印は重複して２個検出されたものです。また、図の左側の二重線部分には二ヶ所の検出漏れがあります。しかし、これらの問題について考え始めると、色々と複雑なケースが想定されるため、自動的に最適解を出力するようにプログラミングするのは大変そうです。正確な交点数が必要なときには、「使い方」に書いた uc_code 出力を利用して、人工ではなく生身の知能で検討・修正いただくことになります。図３には、uc_codeを利用して不確実部分を抽出して描いたグラフを示しています。
しかし、考えようによっては、元々は滑らかだった解析的曲線を折線で代用したときから、誤差の影響は覚悟していたはずです。折線にしてからの交点の問題を、そこまで厳格に詮索するのも無益なことかもしれません。
図３
３－２ 処理時間には短縮の余地あり
曲線群の相互間で、含まれる線分どうしの総当たり検査をするような単純な方法では、図１の結果を出すのに約55分もかかります。現在は、ある程度の工夫をして余計な処理を回避するようにしてはいますが、それでも約18秒かかります（テスト条件....各群内の折点数：22933個 + 24383個、実在するはずの交点数：347個、CPU：Core i5-8250U,1.60GHz）。
例えば、第１曲線内の長めの線分の周辺を、多数の短かい線分で構成された第２曲線が延々と取り巻いているような場合、多量のチェックが必要になって無駄な時間がかかっています。長い方の線分が、水平や垂直から大きく傾くとさらに厳しくなります。折点数が少ない曲線を扱うときにはそれほど気にはなりませんが、多いときにはもう一工夫が必要です。
３－３ その他
不確定情報 uc_code の内容を吟味すると、横長の青色線分に、赤色の縦に短い線分が絡んでいるときに多重検出が起こりやすくなっていることが分かります。しかし、テスト用の曲線モデルは、y軸を等間隔に刻んで計算したx座標の値を使ってプロットしたものなので、水平に近い線分は長く、垂直に近い線分は短くなる傾向にあります。それゆえ、多重検出は線分の縦横には関係なく、２本の長さが極端に違うときに現れやすい現象と考えた方が良さそうです。この原因については引き続き考えてみたいと思います。
４．プログラム
プログラムの最後の方には、コメント化した評価用のプログラムも掲載しています。この評価用プログラムを実行すると、上に示した図１や図３のグラフが表示されます。もし、実行されるのであれば、intersection_of_two_curves.m の 219～222行目あたりにある途中経過表示用のdispコマンド（コメント化されている）の「%」は外してください。結果の表示までに時間がかかるため、そのままではフリーズしたかと勘違いします。
intersection_of_two_curves.m
% intersection_of_two_curves.m
%
% 平面上の２本の曲線(折線)の交点座標を返すユーザー関数。
%
% （動作確認バージョン：MATLAB Home R2019a）
%
% それぞれの曲線は、複数の子曲線を１本にまとめたものでもよい。
% その場合、子曲線間には、区切り折点としてNaNを挿入すること。
% 出力されるのは、第１曲線と第２曲線どうしの交点のみ。同一曲線
% 内の子曲線どうしの交叉は無視される。
%
% [x,y,uc_code]=intersection_of_two_curves(x1,y1,x2,y2)
%
% 【入力】
% x1: 第１曲線の点列のx座標（行ベクトル）
% y1: 同じく、y座標（行ベクトル）
% （普通は、折点数の少ない方を第１曲線にすると処理が速い）
% x2,y2: 第２曲線の、上記同様のx.y座標（行ベクトル）
%
% 【出力】
% x: 交点のx座標（交点の個数に相当する長さの行ベクトル）
% y: 同じく、y座標
% uc_code: 交点の判定が不確実だった場合、その情報（5行×[件数]列）。
% 1行目： 不確実だった第１曲線の線分の若番側の折点番号
% 2行目： 同じく、 第２曲線の線分の若番側の折点番号
% 3,4行目： それらの各線分のどちら側の端が問題だったのか？
% 0: 両端とも問題なし、
% 1: 若番側の端、 2: 後番側の端、
% 3: 端の問題ではなく、両線分の重なり合い。
% 5行目： 疑いながらも、交点座標を出力済のときはその交点番号。
% 交点座標の出力を見合わせたときは0。
function[x,y,uc_code]=intersection_of_two_curves(x1,y1,x2,y2)
% MATLABには、この機能を持つコマンドが無いようなので自作した。
%
% ■残る問題点（検出漏れ、多重検出、解釈の多様性）■
% 一方の曲線のちょうど折点上を、他方の曲線が正確に横切るとき、
% 交点が検出できなかったり、二重検出するような不都合が考えられる。
% ２本の曲線の折点同士が重なったときは、最悪四重検出もあり得る。
% ２本の曲線のそれぞれの線分が完全に重なってしまったときも、
% 「交点なし」ではなく、「無限個の交点」と解釈したい人もいる。
%
% 現状では、uc_code の内容により、利用者の最終判断で修正する必要あり。
%
% ■残る問題点（処理時間の問題）■
% 第１曲線の長めの線分（特に45度程度の傾斜）の縄張り領域に、
% 多数の短小な線分で構成された第２曲線が延々と続いているとき、
% 交点の検出にかなりの時間がかかる。
% 折線数が少ない曲線を扱うときには、あまり気にはならないが、
% 多いときには、もう一工夫が必要。
uc_code=[]; % 交点の判定が不確実だった場合、その情報を収容する変数。
if (size(x1,1)~=1 | size(y1,1)~=1 | size(y2,1)~=1 | size(y2,1)~=1)
disp('エラー： 入力データは行ベクトル形式にしてください')
x=NaN;
y=NaN;
return
end
N1=length(x1); % 入力曲線データのベクトルの長さを取得
N2=length(x2);
if N1~=length(y1)||N2~=length(y2)
disp('エラー： xとyのデータ長が一致しません')
x=NaN;
y=NaN;
return
end
if ~isreal([x1 y1 x2 y2]); % NaNも実数扱いゆえ、処理に問題なし。
disp('エラー： 複素数のデータは使えません')
x=NaN;
y=NaN;
return
end
ind2 = [1:N2]; % 入力された第２曲線の折点要素に番号づけ。処理の途中で
% 不要要素を削除しても、元の折点番号が辿れるように。
% 座標値に対する相対精度の基準（ eps(mar)を折点位置の判定精度にする ）
mar=max(abs([x1 y1 x2 y2]));
% 入力された全曲線の折点の最大の座標値（絶対値）
% NaN要素は無視した数値要素だけの最大値。
x=[]; % 交点座標の出力データ格納用ベクトルを準備。
y=[];
% 交点の判定が不確実だった場合、その該当線分の情報を記録する変数
%（UnCertain）
uc_n1=[]; % 不確実だった第１曲線の線分の若番側の折点番号
uc_n2=[]; % 同じく、第２曲線の線分の若番側の折点番号
uc_x1=[]; % それらの各線分のどちら側の端が問題だったのか？
uc_x2=[]; % どちらの端も問題なし:0、
% 若番側の端:1、後番側の端:2、両線分の重畳:3。
uc_ot=[]; % 不確実ながら交点座標を出力済のときは、その交点
% 番号。交点座標の出力を見合わせたときは0。
for n1=[1:N1-1]; % 第１曲線の折線の線分を、一本ずつ順番に処理。
% 第１曲線の現在の処理線分の、両端点の情報を拾い出す。
% （p,s: primary，secondary）。
% （ここで、p，s各点の座標の大小関係は不定であることに留意しておく
% 必要がある）
x1p=x1(n1); % 現処理線分の、若番側の端点のx座標
y1p=y1(n1); % 同じく、 y座標
x1s=x1(n1+1); % 同じく、 後番側の端点のx座標
y1s=y1(n1+1); % 同じく、 y座標
if isnan(x1p+y1p+x1s+y1s)
continue; % 座標値に一つでもNaNがあれば、この線分の処理はスキッ
% プする。（NaNを持つ線分はグラフ上には
% 描画されない＝存在しない）
end
x2t=x2; % 第２曲線を編集するため、
y2t=y2; % オリジナルからコピーを作る。
ind2t=ind2;
% x2,y2座標値のどちらかにでもNaNがある折点は、
% その両座標値ともNaNにする。
idnan = ~( ~isnan(x2t) .* ~isnan(y2t) ); % NaNがある要素だけ1,他は0。
x2t(idnan==1) = NaN;
y2t(idnan==1) = NaN;
% 処理時間の短縮のため、第２曲線から、判定が必要な線分だけを選び出す。
% 第２曲線の各折線を、交点有無の判定対象とすべきか否かの識別用フラグ
% のベクトル（1:要、判定対象とする、 0:不要、判定の対象外）。
need2 = ones(1,N2); % まずは、折線の全線分を「要」としておく。
% 第１曲線の現在の処理線分の縄張り領域（その線分を包含する傾き0の
% 長方形。周囲に若干のマージンをとる）
st = max([y1p y1s]) + eps(mar); % top の辺のy座標
sb = min([y1p y1s]) - eps(mar); % bottom の辺のy座標
sl = min([x1p x1s]) - eps(mar); % left の辺のy座標
sr = max([x1p x1s]) + eps(mar); % right の辺のy座標
% 第２曲線の線分のうち、この縄張り領域に掛かる可能性がないものを排除
% していく。具体的には、上下左右の４本の境界線（st,sb,sl,srを通る
% 直線）それぞれについて、第２曲線の各折点が境界のどちら側にあるの
% かを判定する。縄張りの反対側にある折点はほぼ不要なので、とりあえ
% ず、need2の該当要素のフラグを1⇒0に書き換える。
% （該当要素だけ0で他は1のベクトルを論理積する）。
% しかし、フラグが0の折点でも、隣の折点のフラグが1であれば、その両
% 折点に挟まれた線分上には交点が存在する可能性がある。この場合、こ
% のフラグ0の折点を無視すると、交点を見落とすことになるので、特例
% の救済措置として、1に書き換える。４本の境界線すべてについて、こ
% の処理を順次繰り返し、不要な要素番号のフラグを0に確定していく。
need2a = (y2t<st) & need2; % 縄張り外の上側にある折点を削除候補化。
need2a = expand_range(need2a); % 上記の救済措置（ローカル関数呼出）
need2b = (y2t>sb) & need2; % 縄張り外の下側にある折点を削除候補化。
need2b = expand_range(need2b);
need2c = (x2t>sl) & need2; % 〃 左側 〃
need2c = expand_range(need2c);
need2d = (x2t<sr) & need2; % 〃 右側 〃
need2d = expand_range(need2d);
% 上の4つの条件を全て満たす第２曲線の折線の線分を求め、
% 交点探索の対象にする。
need2 = need2a & need2b & need2c & need2d;
% しかし、MATLABの文法上、このままでは縄張り領域内にあるNaNの折点
% (消去不可)は、フラグが0に書き換えられて消去対象になってしまって
% いる。NaNの折点が消えると、意識して分断しようとしていた曲線が一
% 本に繋がったように認識される。すると、そこにできた余計な線分によ
% り、誤った位置に架空の交点が検知されるという不都合が起こる。
% これを防止するために、NaNの折点のフラグをすべて1に書き替える。
need2(isnan(x2t)==1) = 1;
% このままでは、縄張り外に、不要なNaNを多く残すことになり、余計な処
% 理が増える。そこで、複数のNaN折点が連続している部分は一つのNaN折
% 点に集約する。
nan1=isnan(x2t);
nan2=circshift(nan1,1);
nan2(1)=0;
need2((nan1 & nan2)==1) = 0;
% 以上により、第２曲線の消去対象折点が確定したので、いよいよ実行。
x2t(need2==0) = []; % need2のフラグが0になっている全折点要素につ
y2t(need2==0) = []; % いて、x,y座標値と折点番号を強制削除する。
ind2t(need2==0) = []; % これで、チェックの対象線分数が大幅に削減
% され、処理が高速化できる。
N2x=length(x2t); % 抽出後の第２曲線の折点の数。
for n2=[1:N2x-1]; % 第２曲線からの抽出線分だけを対象にして、こ
% れらを順番に処理していく。
% 第２曲線の処理線分の両端点の情報を拾い出す
x2p=x2t(n2);
y2p=y2t(n2);
x2s=x2t(n2+1);
y2s=y2t(n2+1);
if isnan(x2p+y2p+x2s+y2s)
continue; % 座標値に一つでもNaNがあれば、この線分はスキップ。
end
% 第１曲線と第２曲線の、各１本の線分同士の交点をチェック。
[xc,yc,pos_1,pos_2,uc_1,uc_2]= ...
line_cross(x1p,y1p,x1s,y1s,x2p,y2p,x2s,y2s,mar);
% ■注意■
% このプログラムと上のローカル関数内では、線の識別記号が異なる。
% このプログラムでの線番号1,2 ⇒ ローカル関数では、a,b。
% このプログラムでの線分の端番号p,s ⇒ ローカル関数では、1,2。
cpout = 0; % 交点座標の出力済みフラグ。
% とりあえずは、未出力とおく（=0）。
if pos_1==0 && pos_2==0; % 真の交点があったときだけ、座標を出力。
% 延長線が作る架空の交点は出力対象外。
x=[x xc];
y=[y yc];
cpout = 1; % 交点の座標を出力済み
% disp(['発見交点数：' num2str(length(x)) '個']); % コマンドラ
% disp([' 線分位置：' num2str(n1) '/' num2str(N1)]); % インに進
% lap=toc; % 行状況を
% disp([' 経過時間：' num2str(lap) '秒']); % 表示。
end
% 交点の有無の判定に自信が持てないとき、それらの情報を記録する。
if ((uc_1~=0 | uc_2~=0) & ~isnan(uc_1));
uc_n1=[uc_n1 n1]; % 第１曲線の折点番号
uc_n2=[uc_n2 ind2t(n2)]; % 第２曲線 〃
uc_x1=[uc_x1 uc_1]; % 第１曲線の折線のどちらの端か？
uc_x2=[uc_x2 uc_2]; % 第２曲線 〃（問題なし:0, 始端:1,
% 終端:2, 両線分の重なり:3）
if cpout==1; % 交点座標を出力済みのとき、
uc_ot=[uc_ot length(x)]; % 交点の要素番号を記録する。
else
uc_ot=[uc_ot 0]; % 交点座標の出力を見合わせたとき0。
end
end
end; % 第２曲線線分によるスキャンの終了
end; % 第１曲線線分によるスキャンの終了
% 出力形式を整えるため、不確実情報を一つの行列にまとめる。
uc_code=[uc_n1;uc_n2;uc_x1;uc_x2;uc_ot];
% =======================================================
% line_cross
%
% ■■ ２本の直線の交点を求めるローカル関数
%
% [x,y,pos_a,pos_b,uc_a,uc_b]= ...
% line_cross(xa1,ya1,xa2,ya2,xb1,yb1,xb2,yb2,xymax)
%
% 【入力】
% xa1,ya1,xa2,ya2 : 線分aの両端の座標。
% xb1,yb1,xb2,yb2 : 線分bの両端の座標。
% xymax : 最大座標値（線分a,bの重なり度合の判定基準の参考値にする）
%
% 【出力】
% x,y : 交点の座標。
% pos_a,pos_b : 線分a,bに対し、交点がどの位置にあるかを示す値。
% （0:線分上、1:端点1側の延長線上、2:端点2側の延長線上）
% uc_a,uc_b : 各線分について、交点有無の判定が確実か？微妙か？
% （0:確実、 3:線分a,b重なり ）
% （1:端点1側で微妙、2:端点2側で微妙）
%
% 平行線同士や、少なくともどちらか一方の線長が0（=点）のときは、
% x, y, pos_a, pos_b は全てNaN。
function [x,y,pos_a,pos_b,uc_a,uc_b] = ...
line_cross(xa1,ya1,xa2,ya2,xb1,yb1,xb2,yb2,xymax)
% 線分a上の任意の点、 線分b上の任意の点を
% x = xa1+ka*(xa2-xa1) x = xb1+kb*(xb2-xb1)
% y = ya1+ka*(ya2-ya1) y = yb1+kb*(yb2-yb1)
% 0 <= ka <= 1 0 <= kb <= 1
% とする。
%
% 両線のx,yを等しいと置くと、
% xa1 + ka*(xa2-xa1) = xb1 + kb*(xb2-xb1)
% ya1 + ka*(ya2-ya1) = yb1 + kb*(yb2-yb1)
% ⇒
% ka*(xa2-xa1) - kb*(xb2-xb1) = xb1-xa1
% ka*(ya2-ya1) - kb*(yb2-yb1) = yb1-ya1
% ⇒
% [xa2-xa1 xb1-xb2] [ka] [xb1-xa1]
% [ ]・[ ] = [ ]
% [ya2-ya1 yb1-yb2] [kb] [yb1-ya1]
%
% このka,kbを計算し、上のx,yの式に代入すれば、交点座標が求まる。
% ただし、ka,kbが上記の範囲から外れているときには、実際の交点は存在
% せず、線分の延長線上でしか交わらない名目だけの交点になる。
A=[xa2-xa1 xb1-xb2 ; ya2-ya1 yb1-yb2];
B=[xb1-xa1;yb1-ya1];
overlap = 0; % a,b両線は重なり合ってはいない（とりあえずの設定）
if rcond(A)<eps(1); % 両線分はほとんど平行で、交点は無さそうなとき。
% rcond(A): Aから精度の良い逆行列が作れない
% とき、微少量を返すコマンド。
% det(A)よりも判定が公正。
% しかし、ほとんど重なり合ってしまっている場合には、判定が微妙なこ
% とを宣言する必要があるので、まずは線分の間隔を計算してみる。
aax=xa2-xa1; % 線分aベクトルのx成分
aay=ya2-ya1; % 〃 y成分
% 線分aの始点から線分bの始点へのベクトルをベクトルcとする。
ccx=xb1-xa1; % cベクトルのx成分
ccy=yb1-ya1; % y成分
% a,bを平行とすれば、
%
% abs(| a × c |) | aax*ccy - aay*ccx |
% 線の間隔 h = ----------------- = ---------------------
% | a | sqrt( aax^2 + aay^2 )
%
h = abs(aax*ccy - aay*ccx) / sqrt(aax^2 + aay^2 ); % 二線の間隔
x=NaN; % まずは、交点が無いことを前提に諸量を設定。
y=NaN; % 線の間隔に関わらず、平行な場合には、策を
pos_a=NaN; % 弄せず、素直に「交点はない」と判定するこ
pos_b=NaN; % とを基本方針とする。
if h<2*eps(xymax); % しかし、線分a,bが殆ど重なっているときには
overlap = 1; % 「交点なし」の一方的な強制には問題もある。
% 交点座標の出力はできないにしても、人為的
% な再判定を促すための不確定警告を出す。
else % 線分a,bが平行で、十分に離れているときは、これ以上
uc_a=NaN; % の吟味は不要ゆえ、呼び出し元にreturnする。
uc_b=NaN;
return;
end
end
if overlap==1; % 線分a,bが殆ど重なっているとき
uc_a=3;
uc_b=3;
else % 交点が存在するはずのとき、その座標を計算。
kk=inv(A)*B; % 逆行列が作れないケースは排除済みゆえ、
ka=kk(1); % ここでエラーは出ないはず。
kb=kk(2);
x = xa1+ka*(xa2-xa1); % 交点の座標（ここでの最重要情報）。
y = ya1+ka*(ya2-ya1);
% 交点の位置が線分の端のギリギリにあって、判定に自信がないとき。
% 不確定警告のためのデータを作る。
% 線分aに関して：
uc_a=0; % まずは0とおいて、判定に自信を表明。
if abs(ka)<eps(1); % しかし、端点1は気になる。
uc_a=1; % 「端点1に自信なし」
elseif abs(1-ka)<eps(1) % あるいは、端点2が気になる。
uc_a=2; % 「端点2に自信なし」
end
% 線分bに関して： 線分aと同様に処理。
uc_b=0;
if abs(kb)<eps(1)
uc_b=1;
elseif abs(1-kb)<eps(1)
uc_b=2;
end
% 自信はあろうがなかろうが、どのような交叉状態だったのかを明示。
if ka>=0 && ka<=1
pos_a=0; % 交点は線分a上にあり。
elseif ka<0
pos_a=1; % 交点は線分aの端点1側の延長線上にあり。
elseif ka>1
pos_a=2; % 交点は線分aの端点2側の延長線上にあり。
else
pos_a=NaN;
end
if kb>=0 && kb<=1
pos_b=0; % 交点は線分b上にあり。
elseif kb<0
pos_b=1; % 交点は線分bの端点1側の延長線上にあり。
elseif kb>1
pos_b=2; % 交点は線分bの端点2側の延長線上にあり。
else
pos_b=NaN;
end
end
end
% =======================================================
% =======================================================
% expand_range
%
% ■■ 処理対象領域を広げるためにベクトルを加工するローカル関数
%
% [flg_new]=expand_range(flg1)
%
% 【入力】
% flg1 : 複数のフラグ要素(0か1)からなる行ベクトル
%
% 【出力】
% flg_new : 加工後のフラグを格納した行ベクトル
% 加工内容： 1に隣り合った0要素を強制的に1に書き換える。
function [flg_new]=expand_range(flg1)
flg2 = circshift(flg1,-1); % フラグベクトルを左シフト。
flg2(length(flg2)) = 0; % 左端は捨てて、右端の空きには0を注入。
flg2 = flg1 | flg2; % オリジナルベクトルとの論理和をとり、
% 1と隣り合った若番側の0を1に書き換え
flg_new = circshift(flg2,1); % 次にこれを右シフト。
flg_new(1) = 0; % 右端は捨てて、左端の空きには0を注入。
flg_new = flg2 | flg_new; % 右シフト前のベクトルとの論理和をとり、
% 1と隣り合った後番側の0を1に書き換え
end
% =======================================================
end % End of "intersection_of_two_curves.m"
% ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
% 以下には、本functionの評価用プログラムをコメント化して掲載している。
% '% % xxtest_intersection_of_two_curves.m' 以下の全行をコピーし、
% 全行頭の'% 'の2文字を削除してから、適当なmファイル名で保存して実行
% します。（Editorの文字のEncodeの種類は、MATLABのコマンドラインから
% 'slCharacterEncoding()'として調べ、それに合わせます。
% ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
% % xxtest_intersection_of_two_curves.m
%
% % function「intersection_of_two_curves」のテスト
%
% close all
% clear
%
% % =========================
% % 評価用のモデル曲線群の作成（ここから） ↓↓↓↓↓
% %（いろいろな離心率の二次曲線）
%
% xd=-10; % 準線の座標
% xf=10; % 焦点の座標
% yy=[0:0.1:60]; % y軸の等間隔サンプリング点
% xyall=[]; % １つの曲線群の原型を収容する変数
%
% % 各種の離心率でスキャン（それぞれが子曲線となる）
% for e=[0.2 0.38 0.52 0.62 0.71 0.78 0.88 1 ...
% 1.2 1.4 1.65 2 2.5 3.2 4.6 8 30]
% if e==1; % 放物線
% x = yy.^2/(2*(xf-xd)) + (xd+xf)/2;
% xy = [x;yy];
% xy2 = xy;
% xy2(:,1) = [];
% xy2 = fliplr(xy2);
% xy2(2,:) = -xy2(2,:);
% xy = [xy2 xy];
% else
% wx = e*(xf-xd)/(1-e^2);
% wy = e*(xf-xd)/sqrt(abs(1-e^2));
% bx = (xd*e^2-xf)/(1-e^2);
% if e<1; % 楕円
% % 上半の曲線
% x1 = -sqrt( wx^2*( 1 - (yy/wy).^2 )) - bx; % 上半左側の曲線
% x2 = sqrt( wx^2*( 1 - (yy/wy).^2 )) - bx; % 上半右側の曲線
% xy = [ [x1 x2];[yy yy] ]; % 左右折点の混合（順不同）
% nf = find( imag(xy(1,:)) ~= 0 ); % 線が存在しない
% xy( :,nf ) = []; % 複素解の折点を削除
% xy = sortrows(xy'); % 左右を混合した折点を
% xy = xy'; % x座標の昇順に並べ替え
% % 下半の曲線
% xy2 = xy; % 上半の曲線をコピー
% xy2(:,1) = []; % 重複点を削除したうえ
% xy2 = fliplr(xy2); % x座標の降順に並べ替え
% xy2(2,:) = -xy2(2,:); % y座標の符号を反転
% xy = [xy2 xy]; % 上半の曲線と下半の曲線を接続
% elseif e>1; % 双曲線
% % 左側
% x1 = -sqrt( wx^2*( 1 + (yy/wy).^2 )) - bx; % 左側上半の曲線
% xy = [x1;yy];
% xy2 = xy; % これをコピーして
% xy2(:,1) = []; % 楕円と同様に下半の曲線
% xy2 = fliplr(xy2); % を作り、上下の曲線を
% xy2(2,:) = -xy2(2,:); % 接続する（左側完成）
% xy = [xy2 xy];
% xym = xy;
% % 右側
% x2 = sqrt( wx^2*( 1 + (yy/wy).^2 )) - bx;
% xy = [x2;yy];
% xy2 = xy;
% xy2(:,1) = [];
% xy2 = fliplr(xy2);
% xy2(2,:) =-xy2(2,:); % 右側も完成
% xy = [xy2 xy [NaN;NaN] xym]; % 左右は分離した別曲線なので
% % NaN折点を介して接続
% end
% end
% xyall = [xyall [NaN;NaN] xy]; % 子曲線を、NaN折点を介して
% % 次々に接続。
% end; % 一つの曲線群の原型が完成
%
% % 曲線群の原型をコピーし、位置を移動した第１曲線と第２曲線を作る
% x1 = xyall(1,:) - 15; % 第１曲線
% y1 = xyall(2,:) + 5;
% x2 = xyall(1,:) + 15; % 第２曲線
% y2 = xyall(2,:) - 5;
%
% CC1 = [x1;y1]; % 第１曲線の折点のx,y座標の行列
% CC2 = [x2;y2]; % 第２曲線 〃
% % グラフから無駄な範囲を省くために、外周を適当にトリミング
% CC1(:,(x1>80|x1<-80|y1>60|y1<-60)==1) = [];
% CC2(:,(x2>80|x2<-80|y2>60|y2<-60)==1) = [];
%
% % 別の曲線を追加するためのスペースを空ける
% % （第１曲線の一部を削除）
% CC1(:,(CC1(1,:)<-48 & CC1(2,:)>-30 & CC1(2,:)<25)==1) = NaN;
%
% % これにより、多くの無駄なNaN折点ができるので、
% % 複数のNaN折点が連続している部分は一つのNaN折点に集約する。
% nan1 = isnan(CC1(1,:));
% nan2 = circshift(nan1,1);
% nan2(1) = 0;
% CC1(:,(nan1 & nan2)==1) = [];
%
% % =============
% % 別種のテスト曲線を追加（二重線の間隔、二重線の平行度 色々）
% % グラフ上では１本に重なっていて見分けがつかないほどの違い
%
% % 第１曲線（傾き45度の線、20本）
% addo = [-75 -71;-28 -24]; % 線のひな型。
% add1 = addo; % これをコピーして上に平行移動し、
% Add1 = []; % 10本の線を配置する。
% for n=1:10
% add1(2,:) = addo(2,:)+4*(n-1);
% Add1 = [Add1 add1 [NaN;NaN]];
% end
% add1 = addo; % 右隣にも同様の10本の線を配置する。
% add1(1,:) = addo(1,:) + 12;
% for n=1:10
% add1(2,:) = addo(2,:) + 4*(n-1);
% Add1 = [Add1 add1 [NaN;NaN]];
% end
% CC1 = [CC1 [NaN;NaN] Add1]; % 第１曲線群に以上の評価曲線を追加
%
% % 第２曲線（平行線の間隔いろいろ）
% add2 = addo; % 左側に作った10本の第１曲線と平行に、
% Add2 = []; % それぞれ第２曲線を引く。
% % 最下の線は第１曲線と完全に重なった線。
% % 上に行くほど少しずつ間隔が大きくなる。
% for n=1:10
% add2(2,:) = addo(2,:) + 4*(n-1) - eps(80)*(n-1);
% Add2 = [Add2 add2 [NaN;NaN]];
% end
%
% % 第２曲線（交叉線の交叉角、いろいろな微少値）
% add2 = addo; % 右側に作った10本の第１曲線の各中央部分に、
% % 非常に狭い角度で交叉する第２曲線を引く。
% % 最下の線は第１曲線と完全に重なった線。
% % 上に行くほど交叉角度は少しずつ大きくなる。
% add2(1,:) = addo(1,:) + 12;
% for n=1:10
% add2(2,1) = addo(2,1) + 4*(n-1) + eps(80)*(n-1)/45;
% add2(2,2) = addo(2,2) + 4*(n-1) - eps(80)*(n-1)/45;
% Add2 = [Add2 add2 [NaN;NaN]];
% end
% CC2 = [CC2 [NaN;NaN] Add2]; % 第２曲線群に以上の評価曲線を追加。
% % =============
%
% % テスト関数への入力用に、曲線座標ベクトルとして分離
% x1 = CC1(1,:);
% y1 = CC1(2,:);
% x2 = CC2(1,:);
% y2 = CC2(2,:);
%
% % 評価用のモデル曲線群の作成（ここまで） ↑↑↑↑↑
% % =========================
%
% figure(1); % 「２つの曲線群どうしの多数の交点の検出」図
% plot(x1,y1,'b',x2,y2,'r'); % 評価用の曲線群を描画
% axis equal
% hold on
% plot(x1,y1,'.b',x2,y2,'.r','MarkerSize',4); % 折点の追加描画
%
% tic; % 処理時間計測用タイマーをスタート
%
% % ■■■ いよいよ、評価対象のユーザー関数を呼び出し
% [x,y,uc_code] = intersection_of_two_curves(x1,y1,x2,y2);
%
% delay=toc; % 処理時間を測定
%
% disp(' ')
% disp(' ')
% disp(' ')
% disp('【不確定情報】')
% disp(' 第１曲線の線分番号')
% disp(' 第２曲線 〃')
% disp(' 第１線分の問題がある端 （0:問題なし、1:始端、2:終端')
% disp(' 第２線分 〃 （3:線の重なり')
% disp(' 問題の交点番号。交点情報なしのときは0')
% % 不確定コードをコマンドラインに書き出す。
% uc_code
% disp(['第１曲線折点数：' num2str(length(x1)) '個'])
% disp(['第２曲線折点数：' num2str(length(x2)) '個'])
% disp(['見つけた交点数：' num2str(length(x)) '個'])
% disp(['要した処理時間：' num2str(delay) '秒'])
%
% % figure(1)に交点やコメントなどを追加描画
% plot(x,y,'o','MarkerEdgeColor',[0 0.6 0],'LineWidth',1)
% nuc = [2 37 39 51 68 104 106 108 110 234]; % 二重検出交点番号
% plot(x(nuc),y(nuc),'o','MarkerEdgeColor','k','LineWidth',1)
% text(-78,20,'各種二重線')
% text(-77,15,'間隔')
% text(-66,15,'交叉角')
% title('２つの曲線群どうしの多数の交点の検出')
% legend('第１曲線','第２曲線','折点１','折点２', ...
% '交点','重複交点','Location','SouthEast')
% ax_fig1 = axis; % figure(2)とのスケール合わせのためのデータを取得
%
% % 不確定部分の周辺要素だけを抜き出して、グラフ上に表示。
% figure(2); % 「交点の判定が不確実だった部分だけを拾い出し」図
% nnn = size(uc_code,2);
% for n=1:nnn
% n1 = uc_code(1,n);
% n2 = uc_code(2,n);;
% h1 = plot([x1(n1) x1(n1+1)],[y1(n1) y1(n1+1)],'b');
% hold on
% h2 = plot([x2(n2) x2(n2+1)],[y2(n2) y2(n2+1)],'r');
% if uc_code(5,n)~=0
% cp = uc_code(5,n);
% h3 = plot(x(cp),y(cp),'o','MarkerEdgeColor',[0 0.6 0], ...
% 'LineWidth',1);
% text(x(cp),y(cp)+n*0.2,num2str(cp))
% else
% h4 = plot((x1(n1)+x1(n1+1))/2,(y1(n1)+y1(n1+1))/2,'xk', ...
% 'LineWidth',1);
% text((x1(n1)+x1(n1+1))/2,(y1(n1)+y1(n1+1))/2+n*0.2, ...
% ['line(' num2str(n1) ')'])
% end
% end
% axis equal;
% axis([ax_fig1(1) ax_fig1(2) ax_fig1(3) ax_fig1(4)]);
% title('交点の判定が不確実だった部分だけを拾い出し')
% legend([h1,h2,h3,h4],'第１曲線','第２曲線','検出した交点', ...
% '不確実第１線分','Location','NorthWest')
% text(-75,-55,['単独の数字は交点番号、' ...
% 'lineに続く数字は第１曲線の不確実線分番号'])
%
% figure(1); % メインのグラフを最上層レイヤーに移動。
% % ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■