MATLABのmesh図は、特別な操作をしなくても、陰線は見えないようにうまく処理されています。それなのに、「この人、何を勘違いしているの？」と疑問に思われることでしょう。
この投稿は、昔のMATLABで、印刷品質を上げるために必要に迫られて作ったユーザーfunctionの紹介です。
古いバージョンでは、ある程度の密度以上のmesh図は、粗いビットマップ図としてしか印刷できないという難点がありました。何とか印刷品質を上げようと、eps出力から抽出した線情報をCADに読み込んで、苦労して手作業で陰線消去し、印刷用のファイルを作ったという苦い思い出があります。
これに懲りて、次の機会にはもっと作業効率を上げようと作ったのがこのfunctionです。自己流の力わざでの労作です。しかし、次の機会は一度も訪れませんでした。さらに、最近のバージョンでは、密度の高いmesh図でもギザギザのない綺麗な線画として印刷できるようになり、もう永久に出番はなくなってしまいました。
そのまま削除するのも悔しいので、ここで最後の紹介をさせていただく次第です。mesh用データをplot3用データに変換するfunctionです。どういう訳か、plot3で描いた図はビットマップ化されずに綺麗に印刷できたことを利用したものです。 今では何の役にも立たないものですが、出来上がった図を別の角度から見ると、破れ傘状態になるのが面白いところです。
一見すると全く同じに見えます（右の図のZ=0のグリッド線の一部が消えているのは、今になって気づいたバグ）。
しかし、視点を変えるとこんなことになってしまいます。
mesh2plot3.m
% mesh2plot3.m
% mesh画面の高品質印刷のためのfunction。mesh用データをplot3用データに変換。
% 使用方法については mファイル内のコメントを参照（末尾に、呼出し元のプログラム例もあり）。
function[XX,YY,ZZ,Xg,Yg,Zg]=mesh2plot3(X,Y,Z,V_view,V_grid,data_aspect)
% 含まれているSubfunction一覧：
% line_reform : マスク更新用データの準備（断面折れ線間の橋渡し）（呼出し元 : mesh2plot3）
% mask_renew : マスク更新 （呼出し元 : mesh2plot3, line_reform, quasi_grid）
% cross_point : 交点計算 （呼出し元 : mask_renew, erase_hidden）
% erase_hidden : 陰線消去 （呼出し元 : mesh2plot3, quasi_grid）
% line_combine : マスク上下の折れ線の結合 （呼出し元 : mesh2plot3）
% restore_3D : 二次元投影されたmesh面を三次元に復元 （呼出し元 : mesh2plot3）
% quasi_grid : 疑似grid線の作成 （呼出し元 : mesh2plot3）
% grid_2D_to_3D : 二次元疑似grid線を三次元に復元 （呼出し元 : quasi_grid）
% 目的：
% MATLABで、meshコマンドによる図をそのまま印刷しても、残念なことに下記の不満がある。
% （ただし、この現象については、バージョンR2007bにおいて確認しているだけであり、他のバー
% ジョンでの動作は不明である）
%
% meshgrid行列のサイズが19行×19列程度以上になると、印刷品質がアウトラインベクトル
% 仕様ではなく、勝手にビットマップ仕様に落とされる。これでは、プレゼンテーション資料
% としての価値がかなり損なわれる。
% mesh図の品質が低下するだけならともかく、同一のfigure内に他のaxes（描画領域）も含ま
% れている場合には、それらの内容についても、文字フォントも含め全てビットマップ品質に
% 落とされてしまう。
%
% ところが、plot3で描いた図には、印刷品質に関する問題はみられない。これを利用し、上記問題
% を解決するために、mesh用のデータをplot3コマンドで描画できる形式に変換するのが本function
% の目的である。手順としては、先ずはmeshコマンドによりお手本となるmesh図を描かせ、この図
% から取得した各種propertyに従って、X,Y,Zデータをplot3用に変換する。この方式で描いた図を、
% 以降、疑似mesh図と呼ぶことにする。処理内容の中心は陰線消去操作となる。
%
% しかし、描画後にマウスで視点変更するなどの便利なGUI操作は無意味なものになる。（操作自体
% は可能だが、あまり役には立たない舞台裏が見えるだけである。）
%
% せっかく陰線消去をしても、plot3画面でgridをonにすると、図形の裏に隠れて見えないはず
% のgrid線まで表示されてしまうため、本来のmesh図とはイメージが異なったものになる。
% これを救済するため、陰線処理済の疑似grid線も出力できるようにした。
% 入力：
% meshgridコマンドで作成された、x-y平面の網目座標を表す行列 ：X, Y
% 網目座標を覆うmesh曲面を表す関数値の行列 ：Z
% 視点位置を表す方位角az[deg]と伏角el[deg]から成るベクトル ：V_view （=[az el]）
% 方位角は、-y軸を0[deg]の基準線とし、z軸の周りに観測者の視点を右へ移動する方向を正と
% する。伏角は、z=0の地平面を基準とし、高所から見下ろす角度を正、地下から見上げる角度
% を負とする。自動設定による az, el の値が不明のときには、お手本となるmesh図がアクテ
% ィブな状態で、[az,el]=view とすることにより取得できる。
% 陰線処理されたgrid線を作成するための情報 ：V_grid。
% この値は、お手本となるmesh図を描かせたあと、下記のコマンドを実行することによって取
% 得することができる。
% V_grid=get(gca,{'XLim' 'XTick' 'YLim' 'YTick' 'ZLim' 'ZTick'});
% x,y,z各軸への割り付け尺度の比のベクトル ：data_aspect
% 視点の修正のために必要となるデータである。お手本となるmesh図を描かせたあと、
% data_aspect=daspect とすることで取得できる。
% ただし、x,y,zの各要素をdaspect(1)などとして直接取得しようとしてもエラーとなる。
% いったん、data_aspectなどのような変数に取り込む必要がある。
%
% 出力：
% mesh図の可視線分だけを集めたベクトル ：XX,YY,ZZ
% これらの出力を使って、plot3(XX,YY,ZZ)などとすれば、meshコマンドによるものと同一形
% 状の曲面が描ける。
% grid線の可視線分だけを集めたベクトル ：Xg,Yg,Zg
% これらの出力を使って、plot3(Xg,Yg,Zg,'-','Color',[0.85 0.85 0.85])などとすれば、
% meshコマンドによるものと類似仕様のgrid線が描ける。
%
% ただし、plot3でお手本どおりの疑似mesh図を得るには、上記で描画した後の処理も重要である。
% function呼出し時に使用した V_view, V_grid の値を利用して下記のように設定すること。
% なお、data_aspectについては設定不要である。
% view(V_view);
% set(gca,'XLim',V_grid{1},'XTick',V_grid{2},'YLim',V_grid{3},'YTick',V_grid{4}, ...
% 'ZLim',V_grid{5},'ZTick',V_grid{6});
% また、お手本を描かせたときに 'axis equal' 等を指示した場合には、plot3上でも同様に
% これを指示すること。
% アルゴリズムを徒に複雑にしないよう、次の条件を設定している。
% mesh面を表すZは、X,Yの一価関数とする。
% mesh面には穴が開いていないこと（Zデータの途中にNaNが入っていないということ）。
% mesh面の色は単色とする。z座標の大きさに応じて線の色を変えるようなことはできない。
% 苦手で特異な視点を指定されたときには、気付かれない程度の微小量（0.01[deg]）だけ移動させ
% て逃げることとする。
% 疑似grid線の線種は、正規のgrid線の仕様とは完全には合わせられないが、これを容認する。
% mesh曲面の縁の一部が疑似grid線と重なって、画面上で正しく表示されないことがある。しかし、
% 印刷は正常に行われるので問題なしとする。
% 伏角が-90[deg]～90[deg]の範囲を超えたときに軸ラベルが誤表示される現象は、MATLABの仕様の
% ようである。これを修正したい場合には、'xlabel'コマンド等の代りに'text'コマンドを使用
% し、文字を本来の正しい位置に強制配置してやれば良い。
% 残る課題：
% az=90*Nのときに疑似grid線を描くと、mesh曲面の左右の縁が表示も印刷もされなくなる。
% 原因の追及ができないまま終了するが、function 'quasi_grid' 内に改良途中の状態をコメント
% として残してある。
% ＝＝＝＝＝＝＝＝＝＝＝＝
% 作成者 ：tsubolabo
% 作成期間 ：2016-5/1～7/20
% 使用Version ：MATLAB R2007b
% 文字code ：Shift-JIS
% 作成環境 ：Windows 7 Professional 32bit, 秀丸
% 利用制限 ：なし
% 賠償責任 ：負いません
% ＝＝＝＝＝＝＝＝＝＝＝＝
az=V_view(1); % 方位角az[deg], azimath
el=V_view(2); % 伏角 el[deg], elevation
% 視点の方位角azと伏角elを、-180<=az<180度,-180<=el<180度の範囲に制限
az = mod(az-180,360) - 180;
el = mod(el-180,360) - 180;
% 視点を真正面や真横（az=90*N）にすると、meshの網目の線が垂直になり、interp1コマンドによる
% 横方向の分割補間などができなくなるため。この条件から僅かにずらす。さらに、az=45+90*Nの場合
% にも一部の変数の'=='判定が微妙になることがあるため、同様に対処する（しかし、後者は他の原因
% による不具合だったのかもしれず、最終的な要否については再確認の必要はある）。
% また、視点が真上や真下になるときも、tan(el)やcot(el)がinfや0になるので、この条件からも
% 僅かにずらす。
da=0.01; % 僅かにずらす角度[deg]
% azをずらす。
spc_ang=[-180 -135 -90 -45 0 45 90 135 180]; % azの特異点
n1ang=find( (spc_ang-da) <= az & az < spc_ang ); % azは上記特異点の近傍の下半分にある？
% ある場合はどの特異点？
n2ang=find( spc_ang <= az & az < (spc_ang+da) ); % azは上記特異点の近傍の上半分にある？
% ある場合はどの特異点？
if ~isempty(n1ang) % azが特異点の近傍の下半分にある。
az=spc_ang(n1ang)-da; % azを近傍の下端に追いやる。
end
if ~isempty(n2ang) % azが特異点の近傍の上半分にある。
az=spc_ang(n2ang)+da; % azを近傍の上端に追いやる。
end
% elをずらす。（azと同様）
spc_ang=[-180 -90 0 90 180]; % elの特異点
n1ang=find( (spc_ang-da) <= el & el < spc_ang );
n2ang=find( spc_ang <= el & el < (spc_ang+da) );
if ~isempty(n1ang)
el=spc_ang(n1ang)-da;
end
if ~isempty(n2ang)
el=spc_ang(n2ang)+da;
end
% 3D図の軸の位置を判定するための、伏角elの象限分け。（実際には、特異点の近傍の微妙なデータは
% 上記で排除されてしまったあとなので、'<=' と '<', '>=' と '>' の厳密な使い分けは必要ない）
if 0<=el && el<90
q_el=1;
elseif 90<=el && el<180
q_el=2;
elseif -180<=el && el<-90
q_el=3;
else % -90<=el && el<=0
q_el=4;
end
% 視点の伏角が-90～90度の範囲を超えたとき、方位角を反対側に移動し、伏角を-90～90度内に収める。
% これはMATLABの動作仕様に合わせるための処理。
if el > 90
az = mod(az,360) - 180;
el = 180 - el;
elseif el < -90
az = mod(az,360) - 180;
el = -180 - el;
end
% 以上で述べた視点位置制限の要点を再掲すると次のようになる。
%
% ① 視点の方位角azと伏角elを-180～180度の範囲に制限。
% ② 視点の伏角が-90～90度の範囲を超えたとき、方位角を反対側に移動し、伏角を-90～90度内に収める。
%
% これらの制限は、本functionのアルゴリズムを分かりやすくするための処置だが、呼出し元のazや
% elの値にまで影響を与えるものではない。呼出し元では、難しいことは考えずに、オリジナルのaz
% とelの値をそのまま使用するだけで良い。
% 以降、方位角と伏角の単位は、[度]ではなく[rad]とする
el = el*pi/180;
az = az*pi/180;
% x,y,z軸の設定をオートスケールにしているときには、DataAspectRatioの影響を受けて、実際の視
% 点の方向は上記のaz,elとは異なるものとなる。よって、下記でその換算を行う。
azz=az;
ell=el;
K1=data_aspect(1)/data_aspect(2);
K2=data_aspect(3)/sqrt(data_aspect(2)^2*cos(azz)^2+data_aspect(1)^2*sin(azz)^2);
az=atan2(sin(azz)*K1,cos(azz));
el=atan2(sin(ell)*K2,cos(ell));
% 3D図を考える時、頭が混乱しやすいことがあるので。後日の見直しに備えて記録しておく。
%
% 例えば、x=[0:4]; y=[2:6]; [X,Y]=meshgrid(x,y); で作成された行列X,Yのデータや、
% Z=X+Y で作成された行列Zのデータは次のような要素配置になる。
%
% 行列X 行列Y 行列Z(=X+Y)
%
% [ 0 1 2 3 4 [ 2 2 2 2 2 [ 2 3 4 5 6
% 0 1 2 3 4 3 3 3 3 3 3 4 5 6 7
% 0 1 2 3 4 4 4 4 4 4 4 5 6 7 8
% 0 1 2 3 4 5 5 5 5 5 5 6 7 8 9
% 0 1 2 3 4 ] 6 6 6 6 6 ] 6 7 8 9 10 ]
%
% これがmeshコマンドでグラフ化されるときには、下図のように、x-y平面上で画面の手前側に高さZの
% データが積み上がる。（下図は、視点を az=0[deg], el=90[deg] として見た状態を示す）
% グラフ化で表示されるデータの並びは、行列Zのデータの並びをそのままx-y平面に貼り付けた形には
% ならず、上下が反転していることに注意！。
%
% Yの列 y
% ベクトル |
% 6- 6 7 8 9 10
% 5- 5 6 7 8 9
% 4- 4 5 6 7 8 ← meshコマンドで描かれる図の例（az=0[deg], el=90[deg]）
% 3- 3 4 5 6 7
% 2- 2 3 4 5 6
% 1- |
% 0- ---+--+--+--+-- x
% 0 1 2 3 4 Xの行ベクトル
%
% 蛇足：
% y=[2:6]; の代りに y=[6:-1:2]; と降順にした場合、グラフ表示は上図のままだが、
% 行列X,Y,Zは下記のようになり、行列の要素配置とグラフのデータ配置が一致する。
%
% 行列X 行列Y 行列Z(=X.*Y)
%
% [ 0 1 2 3 4 [ 6 6 6 6 6 [ 6 7 8 9 10
% 0 1 2 3 4 5 5 5 5 5 5 6 7 8 9
% 0 1 2 3 4 4 4 4 4 4 4 5 6 7 8
% 0 1 2 3 4 3 3 3 3 3 3 4 5 6 7
% 0 1 2 3 4 ] 2 2 2 2 2 ] 2 3 4 5 6 ]
%
% 行列の要素配置とグラフのデータ配置の違いが気になって困るという場合には、この手で逃げる
% ことはできるが、yだけを降順にするというのは少々不自然な感じがする。
% 行列表示とグラフ表示とでは、上下が反転するということを素直に受け入れて、x,yとも昇順に
% 設定しておくのが無難と考えられる。
% 本functionのアルゴリズムの検討も、「x,yとも昇順」を前提にして行っている。
% なお、fliplr, flipud, rot90 などのコマンドにより行列の要素配置を変更する場合、行列X,Y,Zの全
% てに、同時に同一の種類の配置変更を施すのであれば、どのように変更しても、グラフ表示の形は上
% 図から変化はしない。
% なぜなら、全行列に同一の並べ替えをしている限り、X,Y,Zの各要素間の対応関係は崩れることがな
% いので、指定されたx,y座標におけるデータzの値は変わりようがないからである。
% 3D表示の座標軸は、視点がどの象限にあるかによって下表のように変化する。
% ただし、「右軸」・「左軸」は次のように定義する。
% 右軸：指定のaz,elで3D表示されたとき、グラフの右下に右上がり線として示される軸を指す。
% 左軸： 〃 〃 左下に左上がり線 〃
% これに応じて、行列X,Y,Zの要素の並び方を変えて、3D表示したときに最も手前のコーナーにくる網目
% が、行列の第１行・第１列目になるようにする（上図のグラフ化の例では、Z=2の要素が視点に最も近
% くなるように並べ変える）。
% こうすることで、後述の行列A1,B1などの処理を、行番や列番の順に行ってさえいけば、手前の曲線か
% ら順番に陰線消去していくことができて都合が良い。
% （既述のように、この並び替えを行っても、グラフの表示状態が改変される心配はない）
%
% ＝＝＝視点＝＝＝＝ ＝＝＝対応するx,y軸＝＝＝
% 存在 方位角az 左軸 右軸 行列の配置変更操作
% 象限 [deg] 目盛増方向 目盛増方向
%
% 第１ 90 ～ 180 y x X,Y,Zを180度回転
% → ←
% 第２ -180 ～ -90 x y X,Y,Zを上下反転
% ← ←
% 第３ -90 ～ 0 y x そのまま
% ← →
% 第４ 0 ～ 90 x y X,Y,Zを左右反転
% → →
% その前に、陰線消去用のマスクとして利用するため、既存のX,Yよりもやや手前に新たなmeshgridを
% 作成し、その点におけるZの値を一次補間で求めておく。その後、これらも含めて、上表のような
% 行列の配置変更操作を施す。
if az >= pi/2 && az < pi % 第1象限
quad=1; % 元の座標に戻すときのために象限を記録する
Xm=X+(X(1,2)-X(1,1))/10000; % X網目を手前方向（x増加方向）へ微小量だけ移動させXmとする。
Xm(:,size(Xm,2))=X(:,size(X,2)); % 領域はみ出しを制限
Ym=Y+(Y(2,1)-Y(1,1))/10000; % Y網目を手前方向（y増加方向）へ微小量だけ移動させYmとする。
Ym(size(Ym,1),:)=Y(size(Y,1),:); % 領域はみ出しを制限
Zxm=interp2(X,Y,Z,X,Ym,'linear'); % x平行断面処理マスク用
Zym=interp2(X,Y,Z,Xm,Y,'linear'); % y 〃
Xm=rot90(Xm,2); % X,Y,Zを180度回転
Ym=rot90(Ym,2);
Zxm=rot90(Zxm,2);
Zym=rot90(Zym,2);
X=rot90(X,2);
Y=rot90(Y,2);
Z=rot90(Z,2);
elseif az >= -pi && az <= -pi/2 % 第2象限
quad=2;
Xm=X-(X(1,2)-X(1,1))/10000; % X網目を手前方向（x減少方向）へ微小量だけ移動させXmとする。
Xm(:,1)=X(:,1); % 領域はみ出しを制限
Ym=Y+(Y(2,1)-Y(1,1))/10000; % Y網目を手前方向（y増加方向）へ微小量だけ移動させYmとする。
Ym(size(Ym,1),:)=Y(size(Y,1),:); % 領域はみ出しを制限
Zxm=interp2(X,Y,Z,X,Ym,'linear');
Zym=interp2(X,Y,Z,Xm,Y,'linear');
Xm=flipud(Xm); % X,Y,Zを上下反転
Ym=flipud(Ym);
Zxm=flipud(Zxm);
Zym=flipud(Zym);
X=flipud(X);
Y=flipud(Y);
Z=flipud(Z);
elseif az <= 0 && az > -pi/2 % 第3象限
quad=3;
Xm=X-(X(1,2)-X(1,1))/10000; % X網目を手前方向（x減少方向）へ微小量だけ移動させXmとする。
Xm(:,1)=X(:,1); % 領域はみ出しを制限
Ym=Y-(Y(2,1)-Y(1,1))/10000; % Y網目を手前方向（y減少方向）へ微小量だけ移動させYmとする。
Ym(1,:)=Y(1,:); % 領域はみ出しを制限
Zxm=interp2(X,Y,Z,X,Ym,'linear');
Zym=interp2(X,Y,Z,Xm,Y,'linear');
% 第3象限では、配置変更は不要
else % 第4象限（az > 0 & az < pi/2）
quad=4;
Xm=X+(X(1,2)-X(1,1))/10000; % X網目を手前方向（x増加方向）へ微小量だけ移動させXmとする。
Xm(:,size(Xm,2))=X(:,size(X,2)); % 領域はみ出しを制限
Ym=Y-(Y(2,1)-Y(1,1))/10000; % Y網目を手前方向（y減少方向）へ微小量だけ移動させYmとする。
Ym(1,:)=Y(1,:); % 領域はみ出しを制限
Zxm=interp2(X,Y,Z,X,Ym,'linear');
Zym=interp2(X,Y,Z,Xm,Y,'linear');
Xm=fliplr(Xm); % X,Y,Zを左右反転
Ym=fliplr(Ym);
Zxm=fliplr(Zxm);
Zym=fliplr(Zym);
X=fliplr(X);
Y=fliplr(Y);
Z=fliplr(Z);
end
% 視線と垂直な二次平面（直交a,b座標）に、三次元mesh曲面データを投影する。
% （a軸はx-y平面と平行で、視点から見て右向きの軸。b軸は、a軸を反時計方向に90度回転させた
% 向きの軸）
% オリジナルのmesh折れ線を投影して二次元化
A1 = X*cos(az) + Y*sin(az); % mesh曲面上の各網目点の、横軸aの座標
B1 = sin(el)*( Y*cos(az) - X*sin(az) ) + Z*cos(el); % 〃 縦軸bの座標
B_range=[min(B1) max(B1)]; % 縦軸（b座標）のデータ存在領域の最下端と最上端
% ここで、今後多用する用語の定義をしておく：
% x軸平行断面折れ線：
% mesh曲面を、z-x平面と平行な面でスライスしたときの断面形状を表す折れ線
% y軸平行断面折れ線：
% mesh曲面を、y-z平面と平行な面でスライスしたときの断面形状を表す折れ線
% x軸平行断面のマスク用mesh折れ線を投影して二次元化
A1xm = X*cos(az) + Ym*sin(az);
B1xm = sin(el)*( Ym*cos(az) - X*sin(az) ) + Zxm*cos(el);
% y軸平行断面のマスク用mesh折れ線を投影して二次元化
A1ym = Xm*cos(az) + Y*sin(az);
B1ym = sin(el)*( Y*cos(az) - Xm*sin(az) ) + Zym*cos(el);
% この変換により、x軸平行断面折れ線は、a-b面に投影後に、行列A1,B1等の各行のベクトルとして記
% 録されたことになる。A1,B1等の1行目は、3D図上で最も手前に描かれる折れ線に相当し、行番号が
% 大きくなるほど奥の折れ線を表す。
% 同様に、y軸平行断面折れ線は、行列A1,B1等の各列のベクトルとして記録され、1列目が最も手前の
% 折れ線を表す。
% なお、3D図のX,Y,Z行列のn行m列の点に対応する点は、2D図のA1,B1行列等においても、同位置のn行
% m列となっていることに注意。
% まずは、行列A1,B1等の各行の行頭からの処理（x軸平行断面折れ線の処理）を行うことになるが、
% 現在のA1,B1行列等のままだと、これらの処理は次の操作を行うことに相当する。
%
% 視点が第1象限のとき： 右軸に平行な断面折れ線を左から右へ処理して行くことに相当
% 〃 第2象限 〃 ： 左軸 〃 右から左 〃
% 〃 第3象限 〃 ： 右軸 〃 左から右 〃
% 〃 第4象限 〃 ： 左軸 〃 右から左 〃
%
% 次に、行列A1,B1等の各列の列頭からの処理（y軸平行断面折れ線の処理）を行うが、上記と同様に、
% 現在の行列配置のままだと、次の操作に相当していることになる。
%
% 視点が第1象限のとき： 左軸に平行な断面折れ線を右から左へ処理して行くことに相当
% 〃 第2象限 〃 ： 右軸 〃 左から右 〃
% 〃 第3象限 〃 ： 左軸 〃 右から左 〃
% 〃 第4象限 〃 ： 右軸 〃 左から右 〃
%
% しかし、どのケースにおいても、処理が進むにつれてa軸の値が増えて行くように揃えた方が、デー
% タを定形的に扱いやすくて都合が良い。
% そこで、処理開始前に、行列A1,B1等にさらに次のような並べ替えを施す。なお、「y軸平行断面折れ
% 線の処理」における行列A1,B1等の「各列の列頭からの処理」は、転置行列A1',B1'等の「各行の行頭
% からの処理」と等価であることも考慮して変形すると、列に沿っての処理は不要となり、全ての処理
% は行に沿って行うだけでよいことになる。
%
% ＝＝＝＝＝＝A1,B1等の要素の並べ替え＝＝＝＝＝
% 視点の位置 X軸平行断面折れ線 Y軸平行断面折れ線
% を処理するとき を処理するとき
% （A1ym,B1ymは不使用） （A1xm,B1xmは不使用）
% ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
% 第1象限： そのまま 転置後、左右反転
% 第2象限： 左右反転 転置
% 第3象限： そのまま 転置後、左右反転
% 第4象限： 左右反転 転置
for para=[1 2] % 1 or 2: x or y軸に平行な断面折れ線の処理
if para==1
switch quad % 行列A1,B1等を上表のように並べ替えてA2,B2等とする
case 1
A2=A1; B2=B1; Am=A1xm; Bm=B1xm;
case 2
A2=fliplr(A1); B2=fliplr(B1); Am=fliplr(A1xm); Bm=fliplr(B1xm);
case 3
A2=A1; B2=B1; Am=A1xm; Bm=B1xm;
case 4
A2=fliplr(A1); B2=fliplr(B1); Am=fliplr(A1xm); Bm=fliplr(B1xm);
end
else
switch quad % 行列A1,B1等を上表のように並べ替えてA2,B2等とする
case 1
A2=fliplr(A1'); B2=fliplr(B1'); Am=fliplr(A1ym'); Bm=fliplr(B1ym');
case 2
A2=A1'; B2=B1'; Am=A1ym'; Bm=B1ym';
case 3
A2=fliplr(A1'); B2=fliplr(B1'); Am=fliplr(A1ym'); Bm=fliplr(B1ym');
case 4
A2=A1'; B2=B1'; Am=A1ym'; Bm=B1ym';
end
end
% 陰線処理済の折れ線のベクトルを収納するためのセルAout,Boutを準備する（行毎にベクトルの
% 要素数が異なってくるので、行列として統合することはできない）。
% 最も手前の断面折れ線には、視線を遮られるものがないので、A2,B2の1行目の値をそのまま処理
% 済みベクトルとして登録する。
Aout={}; % Aout,Boutのサイズは、paraごとに異なるので、毎回リセットが必要。
Bout={};
Aout{1}=A2(1,:);
Bout{1}=B2(1,:);
% 隠線消去のためのマスクの設定（初期値として、最も手前の折れ線を代入しておく）
% これも、上縁と下縁でベクトルの要素数が異なるので、行列でなくセルで扱う。
% {1}は上縁、{2}は下縁を示す。
Ma{1} = A2(1,:); % マスクの上縁の折れ線のa座標のベクトル
Mb{1} = B2(1,:); % 〃 b座標 〃
Ma{2} = A2(1,:); % マスクの下縁の折れ線のa座標 〃
Mb{2} = B2(1,:); % 〃 b座標 〃
% 行列A2,B2などの2行目以降を処理していく
for line=2:size(A2,1) % 処理すべき行番号（ = 処理すべき断面番号）
for ulx=[1 2] % マスクの上縁(1) or 下縁(2)による陰線処理とマスク更新
% 今回処理すべき折れ線よりも僅か手前にある仮想折れ線Am,Bmでマスクを更新するが、
% 事前に、その極大/極小点で折れ線を補正しておく（網目線が垂直近くになったときに、
% スリットの隙間に隠れている線の消去が不完全になるのを防止するため）。
% 関数呼び出し ================================================
[Amc,Bmc]=line_reform(Am(line,:),Bm(line,:),Am(line-1,:),Bm(line-1,:),ulx,B_range);
% =============================================================
% 補正された折れ線Amc,Bmcで、マスクMa,Mbを更新する。
% 関数呼び出し ================================================
[Ma{ulx},Mb{ulx}]=mask_renew(Ma{ulx},Mb{ulx},Amc,Bmc,ulx,B_range);
% =============================================================
% マスクの背後に隠される断面折れ線を削除する。
% 関数呼び出し ================================================
[A3{ulx},B3{ulx}]=erase_hidden(Ma{ulx},Mb{ulx},A2(line,:),B2(line,:),ulx);
% =============================================================
end % end of 'for ulx'
% マスクより上の折れ線A3{1},B3{1}と、マスクより下の折れ線A3{2},B3{2}の結合
% 関数呼び出し ================================================
[A4,B4]=line_combine(A3{1},B3{1},A3{2},B3{2});
% =============================================================
Aout{line}=A4;
Bout{line}=B4;
if line==size(A2,1) % x軸 or y軸に平行な全ての断面折れ線の処理を完了
if para==1
Ma1=Ma; % x軸平行断面折れ線処理完了時の最大マスクの
Mb1=Mb; % 上縁（Ma1{1},Mb1{1}）と下縁（Ma1{2},Mb1{2}）
else
Ma2=Ma; % y軸平行断面折れ線処理完了時の最大マスクの
Mb2=Mb; % 上縁（Ma2{1},Mb2{1}）と下縁（Ma2{2},Mb2{2}）
end
end
end % end of 'for line'
if para==1
Aout1=Aout; % x軸平行断面折れ線のa座標値
Bout1=Bout; % 〃 b座標値
else
Aout2=Aout; % y軸平行断面折れ線のa座標値
Bout2=Bout; % 〃 b座標値
end
end % end of 'for para'
% a,b座標の２次元データから、plot3用のx,y,z座標の３次元データを復元する。
% 関数呼び出し ================================================
[XX,YY,ZZ]=restore_3D(Aout1, Bout1, Aout2, Bout2, X(1,:), Y(:,1), az, el, quad, V_grid);
% =============================================================
% 疑似grid線の生成
% 関数呼び出し ================================================
[Xg,Yg,Zg]=quasi_grid(V_grid, az, el, quad, q_el, Ma1, Mb1, Ma2, Mb2);
% =============================================================
end % end of function 'mesh2plot3'（functionの'end'の後ろに';'をつけるとエラー）
% ■■■===============================================================
function[A3,B3]=line_reform(A1,B1,A2,B2,ulx,B_range)
% 視点azが90[deg]の整数倍近くになると、網目の線のうち垂直に近くなる方の陰線消去が不完全になる。
% これを防止するために、マスク修正用に使用する断面折れ線に、事前に下記のように手を加える。
% 断面折れ線A1,B1の極大点/極小点（ulx=1/ulx=2のとき）を検出し、そこから、一つ手前の断面折れ
% 線A2,B2上の隣接する点に線を引き、A_bridge,B_bridgeとする。（ここでは、新旧断面折れ線間の
% 架け橋に見立てて、「ブリッジ線」と呼ぶことにする）
% 断面折れ線A1,B1とブリッジ線A_bridge,B_bridgeとの大きい方/小さい方（ulx=1/ulx=2のとき）の
% 包絡線をA3,B3として出力する。
% 入力：
% 現在のマスク修正用断面折れ線のa,b座標ベクトル ：A1, B1
% 一つ手前のマスク修正用断面折れ線のa,b座標ベクトル ：A2, B2
% 包絡線の作り方 ：ulx = 1 上側優先
% = 2 下側優先
% 二次元化mesh曲面の最下端と最上端 ：B_range
%
% 出力：
% リフォーム後のマスク修正用折れ線のa,b座標ベクトル ：A3, B3
% 極大点（ulx=1の場合）/極小点（ulx=2の場合）の検出
if ulx==1
sgn=1; % 包絡線の拡張方向に応じて、判定条件の不等式を反転させるために利用
else
sgn=-1;
end
grad=B1(2:length(B1))-B1(1:length(B1)-1); % B1ベクトルの各要素間の差
F_peak=find(sgn*grad(2:length(grad))<=0 & sgn*grad(1:length(grad)-1)>=0);
% ただし、F_peak中の連続する要素番号（多区間に渡る水平線がある場合に生じる）は先頭だけ残し
% て削除
F_cont=find(F_peak-[-inf F_peak(1:length(F_peak)-1)]==1);
F_peak(F_cont)=[]; % A1,B1ベクトルのF_peak+1番目の要素が極大/極小点である。
% 断面折れ線の両端も、極大/極小点に準じて扱う
F_peak=[0 F_peak length(B1)-1];
if ulx==1 % 折れ線領域より十分に下の方に、ブリッジ線の橋脚の基礎を設定。
base=B_range(1)-100*(B_range(2)-B_range(1)+1); % 「+1」は、狭Rangeの場合のエラー防止。
else % 折れ線領域より十分に上の方に、ブリッジ線の吊り下げ点を設定。
base=B_range(2)+100*(B_range(2)-B_range(1)+1);
end
tilt=(A1(2)-A1(1))/1000; % 橋脚/吊り縄の広がり代（垂直線が存在すると具合が悪いため）
A_left=A1(1); % 断面折れ線領域の左端
A_right=A1(length(A1)); % 同じく、右端
A3x=A1;
B3x=B1;
for n=F_peak+1 % 各極大/極小点について処理を繰り返す。
% ブリッジ線を作る。
A_bridge=[A1(n) A2(n)];
B_bridge=[B1(n) B2(n)];
[A_bridge,ID]=sort(A_bridge); % 視点によつてA_bridgeの並びの方向が変わるので、昇順に統一。
B_bridge=B_bridge(ID);
A_bsl=A_bridge(1)-tilt; % bridge_sustain_left
A_bbl=A_bridge(1); % bridge_body_left
A_bbr=A_bridge(length(A_bridge)); % bridge_body_right
A_bsr=A_bridge(length(A_bridge))+tilt; % bridge_sustain_right
B_bbl=B_bridge(1);
B_bbr=B_bridge(2);
% ブリッジを橋脚/吊り縄で支える。
if ( A_left<A_bsl && A_bsr<A_right )
A_bridge=[A_bsl A_bbl A_bbr A_bsr];
B_bridge=[base B_bbl B_bbr base];
elseif ( A_bsl<=A_left && A_bsr<A_right )
A_bridge=[A_bbl A_bbr A_bsr];
B_bridge=[B_bbl B_bbr base];
elseif ( A_left<A_bsl && A_right<=A_bsr )
A_bridge=[A_bsl A_bbl A_bbr];
B_bridge=[base B_bbl B_bbr];
else
A_bridge=[A_bbl A_bbr];
B_bridge=[B_bbl B_bbr];
end
% 断面折れ線とブリッジ線の、大きい方/小さい方 の包絡線を求める
% 関数呼び出し ================================================
[A3x,B3x]=mask_renew(A3x,B3x,A_bridge,B_bridge,ulx,B_range);
% =============================================================
end
A3=A3x;
B3=B3x;
end % end of function 'line_reform'（functionの'end'の後ろに';'をつけるとエラー）
% ■■■===============================================================
function[A3,B3]=mask_renew(A1,B1,A2,B2,ulx,B_range)
% 既存マスクの折れ線データと、マスク修正用折れ線データを入力すると、それらの包絡線のデータ
% を出力する。上の包絡線を採るか、下の包絡線を採るかは、ulxで指定する。
% 2本の折れ線の両端のa座標が揃っていない場合には、まず、下記のようにして揃える。
% 端までの長さが不足している折れ線に線分を一本追加する。それによりできる新たな端点のa座
% 標は、相手の折れ線の端点と同一とし、b座標は条件に応じて次のように設定する。
% 上の包絡線を出力に指定した場合 → 非常に小さい値（気持ち的には-inf）
% 下の 〃 → 非常に大きい値（ 〃 +inf）
% （要するに、追加した線分が包絡線出力に影響を与えないように選ぶ）
% 入力：
% 既存マスクの折れ線のa,b座標ベクトル ：A1, B1
% マスク修正用折れ線のa,b座標ベクトル ：A2, B2
% A1,A2の要素は昇順に並んでいること。また、A1,A2,B1,B2にNaNは含まれないこと。
% 両折れ線とも、垂直線を含んでいないこと。
% マスク更新の方法 ：ulx = 1 上に拡張する方向に更新
% = 2 下に拡張する方向に更新
% 二次元化mesh曲面の最下端と最上端 ：B_range
%
% 出力：
% 包絡線のa,b座標ベクトル ：A3, B3
if ulx==1
sgn=1; % 包絡線の拡張方向に応じて、判定条件の不等式を反転させるために利用。
big=B_range(1)-100*(B_range(2)-B_range(1)+1); % 領域外のダミー巨大データ
% 「+1」は、狭Rangeの場合のエラー防止。
else
sgn=-1;
big=B_range(2)+100*(B_range(2)-B_range(1)+1);
end
if A1(1)~=A2(1) % 左端が揃っていない
if A1(1)>A2(1) % 折れ線A1,B1の左側のデータが不足
A1=[A2(1) A1]; % 不足データを追加
B1=[big B1];
else % 折れ線A2,B2の左側のデータが不足
A2=[A1(1) A2]; % 不足データを追加
B2=[big B2];
end
end
if A1(length(A1))~=A2(length(A2)) % 右端が揃っていない
if A1(length(A1))<A2(length(A2)) % 折れ線A1,B1の右側のデータが不足
A1=[A1 A2(length(A2))]; % 不足データを追加
B1=[B1 big];
else % 折れ線A2,B2の右側のデータが不足
A2=[A2 A1(length(A1))]; % 不足データを追加
B2=[B2 big];
end
end
% 両折れ線のa座標を混ぜ合わせて昇順に並べ替え、重複するものは一つにまとめる。
A_merge=sort([A1 A2]); % 混合してソート
n=find( ([A_merge inf]-[-inf A_merge] )==0); % 重複する要素番号を拾い出す
A_merge(n)=[]; % 重複要素を削除
% 折れ線座標B1,B2に共通のA座標を持たせてB1x,B2xとする。
% （両線の形状をそのままに、必要最少数の共通折れ点による表現に変更する）
B1x=interp1(A1,B1,A_merge);
B2x=interp1(A2,B2,A_merge);
Ad_point=[]; % 追加する折れ点
for n=1:length(A_merge)-1 % A_mergeの全区間をスキャン
if (B2x(n)-B1x(n))*(B2x(n+1)-B1x(n+1))<0 % ２線が交叉する区間だけを拾い出す
% 関数呼び出し ================================================
[x,y,pa,pb]=cross_point(A_merge(n),A_merge(n+1),B1x(n),B1x(n+1),...
A_merge(n),A_merge(n+1),B2x(n),B2x(n+1)); % 交点の計算
% =============================================================
if pa==1 % 計算誤差が疑われるような微妙な交点ではない。
Ad_point=[Ad_point x]; % 交点のa座標を追加折れ点として登録していく。
end
end
end
% 交点座標を追加して、A_merge,B1x,B2x を再構成
A_merge=sort([A_merge Ad_point]); % 混合してソート
n=find([A_merge inf]-[-inf A_merge]==0); % 重複する要素番号を拾い出す
A_merge(n)=[]; % 重複要素を削除
B1x=interp1(A1,B1,A_merge);
B2x=interp1(A2,B2,A_merge);
% 共通a座標A_mergeの各点は、元々のA1やA2に含まれていたのと同一点かどうか？の記録。
% 後続処理で、直線の過剰な１本化を防止するために使用。
% Na1=1 → A1に含まれていた点
% =0 → A1には含まれていなかった点
% Na2=1 → A2に含まれていた点
% =0 → A2には含まれていなかった点
Na1=[];
n=1; % A_mergeの要素番号ポインタ
for m=1:length(A1) % A1の全要素をスキャン
while A_merge(n)<=A1(m)
if A_merge(n)==A1(m)
Na1=[Na1 1];
else
Na1=[Na1 0];
end
n=n+1;
if n>length(A_merge)
break;
end
end
end
Na1([1 length(Na1)])=[]; % 後続のN_delのfind対象とのベクトルサイズ合わせ
Na2=[];
n=1; % A_mergeの要素番号ポインタ
for m=1:length(A2) % A2の全要素をスキャン
while A_merge(n)<=A2(m)
if A_merge(n)==A2(m)
Na2=[Na2 1];
else
Na2=[Na2 0];
end
n=n+1;
if n>length(A_merge)
break;
end
end
end
Na2([1 length(Na2)])=[]; % 後続のN_delのfind対象とのベクトルサイズ合わせ
% 以上の処理により、A_mergeで区切られた区間内では、２本の線分の中途半端な関係は存在できなく
% なり、相手より完全に上か？下か？を明確に判別できる。（完全に同一という関係は残るが）
% しかし、計算誤差により、想定外の問題が生じる可能性はある。例えば、区間の片方の端点だけの
% b座標の大小関係だけでは、上下の判定が狂うことがあり得る。そこで、両端点のb座標の差の合計
% で大小関係を判定することにする。
B3x=[]; % 暫定的な更新マスクのb座標ベクトル
chosen=[]; % 既存マスク(1)・修正折れ線(2)のどちらが採用されたか？を記録
for n=1:length(A_merge)-1 % A_mergeの全区間をスキャン
if sgn*( (B2x(n)-B1x(n)) + (B2x(n+1)-B1x(n+1)) ) >= 0 % 修正折れ線を新規マスクに採用
if n==1
B3x(n)=B2x(n);
chosen(n)=2;
end
B3x(n+1)=B2x(n+1);
chosen(n+1)=2;
else % 既存マスクをそのまま維持
if n==1
B3x(n)=B1x(n);
chosen(n)=1;
end
B3x(n+1)=B1x(n+1);
chosen(n+1)=1;
end
end
chosen([1 length(chosen)])=[]; % 後続のN_delのfind対象とのベクトルサイズ合わせ
% 以上の処理によって、本来は１本の直線だったものが分断されてしまうことがある。
% これによって、メモリ消費も増えるし、後処理段階において、何らかの想定外の事態を生じる可能
% 性もある。これを避けるために、本来の１本の直線に戻す。
% 具体的には、隣接区間同士で線分の傾きが同一の場合には、両区間の境界にある折れ点を削除する。
% しかし、削除が行き過ぎると、長区間で滑らかな折れ線は、強制的に直線化されて歪んでしまう可
% 能性がある。そこで、オリジナルのA1,A2の線分単位を超えての１本化は行わないことにする。
delta_A = A_merge(2:length(A_merge)) - A_merge(1:length(A_merge)-1);
delta_B = B3x(2:length(A_merge)) - B3x(1:length(A_merge)-1);
grad=atan(delta_B./delta_A);
N_del=find(...
( abs(grad(2:length(grad))-grad(1:length(grad)-1)) < 0.001*pi/180 )...
& ~( (chosen==1 & Na1==1) | (chosen==2 & Na2==1) )...
);
A_merge(N_del+1)=[]; % 折れ点のa座標を消去
B3x(N_del+1)=[]; % 折れ点のb座標を消去
A3=A_merge;
B3=B3x;
end % end of function 'mask_renew'（functionの'end'の後ろに';'をつけるとエラー）
% ■■■===============================================================
function[x,y,pa,pb]=cross_point(x1a,x2a,y1a,y2a,x1b,x2b,y1b,y2b)
% 2本の線分の両端点の座標を入力すると、交点の座標を出力する。
% （注：他のfunctionとは異なり、ここでのa,bは、平面座標a,bを示す文字ではなく、二本の各線分名
% を示すものである。座標はx,yで表している。）
%
% 入力：
% 線分aの左側の端点の座標 ：x1a, y1a
% 〃 右側 〃 ：x2a, y2a
% 線分bの左側の端点の座標 ：x1b, y1b
% 〃 右側 〃 ：x2b, y2b
%
% 出力：
% 交点の座標 ：x, y
% a線、b線に対する交点の位置 ：pa, pb
% =0 線の左端よりも左にある
% =1 線上にある
% =2 線の右端よりも右にある
% =3 交点は無い
% 傾斜の計算
Ga=(y2a-y1a)/(x2a-x1a);
Gb=(y2b-y1b)/(x2b-x1b);
% 交点座標の計算
if abs(Ga)==inf && abs(Gb)==inf % 両線とも垂直線
x=NaN; % 交点は無い
y=NaN;
elseif abs(Ga)==inf % a線だけが垂直線
x=x2a;
y=Gb*(x-x1b)+y1b;
elseif abs(Gb)==inf % b線だけが垂直線
x=x2b;
y=Ga*(x-x1a)+y1a;
else % どちらの線も垂直線ではない
if Ga==Gb % 両線は同一傾斜
x=NaN; % 交点は無い
y=NaN;
else % 両線とも特殊性は無く、交点は一般式で計算できる。
x=(Ga*x1a-Gb*x1b-y1a+y1b)/(Ga-Gb);
y=Ga*(x-x1a)+y1a;
end
end
% 交点の位置を調べる
if isnan(x) && isnan(y)
pa=3;
pb=3;
else
if x<x1a
pa=0;
elseif x>x2a
pa=2;
else
pa=1;
end
if x<x1b
pb=0;
elseif x>x2b
pb=2;
else
pb=1;
end
end
end % end of function 'cross_point'（functionの'end'の後ろに';'をつけるとエラー）
% ■■■===============================================================
function[A3,B3]=erase_hidden(A1,B1,A2,B2,ulx)
% マスクの縁の折れ線と陰線処理対象の断面折れ線を入力し、見せたいのはマスクの上側か？下側か？
% をulxで指定することによって、陰線処理された断面折れ線のデータが返される（隠された部分のb座
% 標はNaNで置き換えられる）。
% 2本の折れ線の両端のa座標が揃っていない場合は、次のように加工してから処理する。
% 両折れ線の左端同士、右端同士を直線で結び、追加した各直線は、端のデータが不足していた
% 方の折れ線に属していたものとする。
%
% 入力：
% マスクの縁の折れ線のa,b座標ベクトル ：A1, B1
% 処理対象折れ線のa,b座標ベクトル ：A2, B2
% A1,A2の要素は昇順に並んでいること。また、A1,A2,B1,B2にNaNは含まれないこと。
% 両折れ線とも、垂直線を含んでいないこと。
% 見せたい線は上か下かの指定 ：ulx
% = 1 マスク線の上側を可視線として残す。
% = 2 〃 下側 〃
%
% 出力：
% 陰線処理済の折れ線のa,b座標ベクトル ：A3, B3
% function 'mask_renew'と類似の処理が多く含まれるが、一つに纏めることに拘り過ぎると想定外の
% バグが出そうなので、別functionとして独立したままとする。
minA2=[]; % 処理対象折れ線の左端が不足していた場合、追加前の左端のa座標
maxA2=[]; % 処理対象折れ線の右端が不足していた場合、追加前の右端のa座標
if A1(1)~=A2(1) % 左端が揃っていない
if A1(1)>A2(1) % 折れ線A1,B1の左側のデータが不足
A1=[A2(1) A1]; % 不足データを追加
B1=[B2(1) B1];
else % 折れ線A2,B2の左側のデータが不足
minA2=A2(1);
A2=[A1(1) A2]; % 不足データを追加
B2=[B1(1) B2];
end
end
if A1(length(A1))~=A2(length(A2)) % 右端が揃っていない
if A1(length(A1))<A2(length(A2)) % 折れ線A1,B1の右側のデータが不足
A1=[A1 A2(length(A2))]; % 不足データを追加
B1=[B1 B2(length(B2))];
else % 折れ線A2,B2の右側のデータが不足
maxA2=A2(length(A2));
A2=[A2 A1(length(A1))]; % 不足データを追加
B2=[B2 B1(length(B1))];
end
end
% 両折れ線のa座標を混ぜ合わせて昇順に並べ替え、重複するものは一つにまとめる。
A_merge=sort([A1 A2]); % 混合してソート
n=find([A_merge inf]-[-inf A_merge]==0); % 重複する要素番号を拾い出す
A_merge(n)=[]; % 重複要素を削除
% 折れ線座標B1,B2に共通のA座標を持たせてB1x,B2xとする。
% （両線の形状をそのままに、必要最少数の共通折れ点による表現に変更する）
B1x=interp1(A1,B1,A_merge);
B2x=interp1(A2,B2,A_merge);
Ad_point=[]; % 追加する折れ点
for n=1:length(A_merge)-1 % A_mergeの全区間をスキャン
if (B2x(n)-B1x(n))*(B2x(n+1)-B1x(n+1))<0 % ２線が交叉する区間だけを拾い出す
% 関数呼び出し ================================================
[x,y,pa,pb]=cross_point(A_merge(n),A_merge(n+1),B1x(n),B1x(n+1),...
A_merge(n),A_merge(n+1),B2x(n),B2x(n+1)); % 交点の計算
% =============================================================
if pa==1 % 計算誤差が疑われるような微妙な交点ではない。
Ad_point=[Ad_point x]; % 交点のa座標を追加折れ点として登録していく。
end
end
end
% 交点座標を追加して、A_merge,B1x,B2x を再構成
A_merge=sort([A_merge Ad_point]); % 混合してソート
n=find([A_merge inf]-[-inf A_merge]==0); % 重複する要素番号を拾い出す
A_merge(n)=[]; % 重複要素を削除
B1x=interp1(A1,B1,A_merge);
B2x=interp1(A2,B2,A_merge);
% 共通a座標A_mergeの各点は、元々のA1やA2に含まれていたのと同一点かどうか？の記録。
% 後続処理で、直線の過剰な１本化を防止するために使用。
% Na1=1 → A1に含まれていた点
% =0 → A1には含まれていなかった点
% Na2=1 → A2に含まれていた点
% =0 → A2には含まれていなかった点
Na1=[];
n=1; % A_mergeの要素番号ポインタ
for m=1:length(A1) % A1の全要素をスキャン
while A_merge(n)<=A1(m)
if A_merge(n)==A1(m)
Na1=[Na1 1];
else
Na1=[Na1 0];
end
n=n+1;
if n>length(A_merge)
break;
end
end
end
Na1([1 length(Na1)])=[]; % 後続のN_delのfind対象とのベクトルサイズ合わせ
Na2=[];
n=1; % A_mergeの要素番号ポインタ
for m=1:length(A2) % A2の全要素をスキャン
while A_merge(n)<=A2(m)
if A_merge(n)==A2(m)
Na2=[Na2 1];
else
Na2=[Na2 0];
end
n=n+1;
if n>length(A_merge)
break;
end
end
end
Na2([1 length(Na2)])=[]; % 後続のN_delのfind対象とのベクトルサイズ合わせ
% 以上の処理により、A_mergeで区切られた区間内では、２本の線分の中途半端な関係は存在できなく
% なり、相手より完全に上か？下か？を明確に判別できる。（完全に同一という関係は残るが）
% しかし、計算誤差により、想定外の問題が生じる可能性はある。例えば、区間の片方の端点だけの
% b座標の大小関係だけでは、上下の判定が狂うことがあり得る。そこで、両端点のb座標の差の合計
% で大小関係を判定することにする。
if ulx==1
sgn=1; % 見せたい線が上下のどちらかに応じて、判定条件の不等式を反転させるために利用
else
sgn=-1;
end
B3x=[]; % 暫定的な陰線処理結果のb座標ベクトル
chosen=[]; % 陰線(1)・見える線(2)のどちらに割り付けられたか？を記録
for n=1:length(A_merge)-1 % A_mergeの全区間をスキャン
if sgn*( (B2x(n)-B1x(n)) + (B2x(n+1)-B1x(n+1)) ) >= 0 % この線分は見える
if n==1
B3x(n)=B2x(n);
chosen(n)=2;
end
B3x(n+1)=B2x(n+1);
chosen(n+1)=2;
else % この線分は見えない（陰線）
if n==1
B3x(n)=NaN;
chosen(n)=1;
end
if sum(Ad_point==A_merge(n+1))==0
B3x(n+1)=NaN;
else
B3x(n+1)=B1x(n+1); % そこが交点なら、次の区間では可視線分ゆえ、折れ点として有効。
end
chosen(n+1)=1;
end
end
chosen([1 length(chosen)])=[]; % 後続のN_delのfind対象とのベクトルサイズ合わせ
% 以上の処理によって、本来は１本の直線だったものが分断されてしまうことがある。
% これによって、メモリ消費も増えるし、後処理段階において、何らかの想定外の事態を生じる可能
% 性もある。これを避けるために、本来の１本の直線に戻す。
% 具体的には、隣接区間同士で線分の傾きが同一の場合には、両区間の境界にある折れ点を削除する。
% しかし、削除が行き過ぎると、長区間で滑らかな折れ線は、強制的に直線化されて歪んでしまう可
% 能性がある。そこで、オリジナルのA1,A2の線分単位を超えての１本化は行わないことにする。
delta_A = A_merge(2:length(A_merge)) - A_merge(1:length(A_merge)-1);
delta_B = B3x(2:length(A_merge)) - B3x(1:length(A_merge)-1);
grad=atan(delta_B./delta_A);
N_del=find(...
( abs(grad(2:length(grad))-grad(1:length(grad)-1)) < 0.01*pi/180 )...
& ~( (chosen==1 & Na1==1) | (chosen==2 & Na2==1) )...
);
A_merge(N_del+1)=[]; % 折れ点のa座標を消去
B3x(N_del+1)=[]; % 折れ点のb座標を消去
% さらに、B3xに連続した複数のNaNがある場合には、一つのNaN要素に集約する。
C=isnan(B3x);
Ca=find(C==1);
Cb=[Ca inf]-[inf Ca];
Cb(length(Cb))=[];
Ca(Cb~=1)=[];
A_merge(Ca)=[];
B3x(Ca)=[];
% さらに、B3xの最初と最後にNaN要素がある場合（あるとしても、それぞれに１つずつのはず）には、
% それらの点は削除する。
if isnan(B3x(1))
A_merge(1)=[];
B3x(1)=[];
end
if ~isempty(B3x) && isnan(B3x(length(B3x)))
A_merge(length(A_merge))=[];
B3x(length(B3x))=[];
end
% 最初の処理の段階で両端揃えのためにA2,B2に追加した線分に相当する部分は削除。
if ~isempty(minA2)
N=find(A_merge<minA2);
A_merge(N)=[];
B3x(N)=[];
end
if ~isempty(maxA2)
N=find(A_merge>maxA2);
A_merge(N)=[];
B3x(N)=[];
end
A3=A_merge;
B3=B3x;
end % end of function 'erase_hidden'（functionの'end'の後ろに';'をつけるとエラー）
% ■■■===============================================================
function[A3,B3]=line_combine(A1,B1,A2,B2)
% 陰線処理済の、マスクより上にある折れ線と下にある折れ線を結合して、一つのベクトルに集約す
% る。陰線処理用function 'erase_hidden'の出力の仕様を満たしていることを前提にしたアルゴリズ
% ムなので、汎用性は小さい。
%
% 入力：
% マスクより上にある折れ線のa,b座標ベクトル ：A1, B1
% 〃 下 〃 ：A2, B2
%
% 出力：
% 結合された折れ線のa,b座標ベクトル ：A3, B3
m=0;
AG={};
BG={};
for ulx=[1 2]
if ulx==1
A=A1;
B=B1;
else
A=A2;
B=B2;
end
C=isnan(B); % B1やB2について、NaNを'1'、その他を'0'として置き換えたベクトル
D=ones(1,length(C)); % 全要素が'1'のベクトル
D=conv(C,D); % Cを積分
D=D(1:length(C)); % ここまでの C→D変換例
% C=[0 0 0 0 1 0 0 1 0 0 0 1 0 0]
% D=[0 0 0 0 1 1 1 2 2 2 2 3 3 3]
F=find(C==1); % F=[5 8 12]
A(F)=[]; % B1,B2にNaNを含む部分のa座標を除去
B(F)=[]; % B1,B2のNaN要素を除去
D(F)=[]; % 例：D=[0 0 0 0 1 1 2 2 2 3 3]
for n=min(D):max(D) % NaNを含んだ一つの分断折れ線ベクトルを、
% 複数の連続折れ線ベクトルに分離してセルに収納。
% マスクより上か下かにはお構いなく順不同で一つのセルにぶち
% 込む。
m=m+1;
H=find(D==n);
AG{m}=A(H);
BG{m}=B(H);
end
end
% この段階では、どの折れ線ベクトルにもNaN要素は含まれなくなっている。
% 複数の連続折れ線をa座標の昇順に結合し、一つのベクトルに結合する。
% そのとき、各連続線の境界にはNaN要素を挿入する。
for n=1:m
J(n)=min(AG{n}); % 各連続折れ線の最小のa座標値（nはセル番号に対応）
end
[K,L]=sort(J); % L: a座標が小さい順に並べるためのセル番号のインデックス。
% K: 各折れ線のa座標の最小値を小さい順に並べたもの（使用せず）。
A3=[]; % 出力用のベクトルを準備。
B3=[];
n=0;
for m=L % a座標が小さい順にセルの内容を繋ぎ合わせる。
n=n+1;
if m~=L(1) && AG{n-1}(length(AG{n-1}))~=AG{n}(1)
A3=[A3 ( AG{n-1}(length(AG{n-1}))+AG{n}(1) )/2]; % 連続線の境界にNaN要素を挿入
B3=[B3 NaN];
end
A3=[A3 AG{m}];
B3=[B3 BG{m}];
end
end % end of function 'line_combine'（functionの'end'の後ろに';'をつけるとエラー）
% ■■■===============================================================
function[XX,YY,ZZ]=restore_3D(Aout1, Bout1, Aout2, Bout2, X_v, Y_v, az, el, quad, V_grid)
% ２次元データから、plot3用の３次元データを復元する。
% 入力：
% x軸平行断面折れ線のa,b座標値 ：Aout1, Bout1
% y軸平行断面折れ線のa,b座標値 ：Aout2, Bout2
% 網目点のx座標,y座標のベクトル ：X_v, Y_v
% 視点の方位角と伏角 ：az, el
% 視点azによる象限分け ：quad
% grid線用の座標情報 ：V_grid
%
% 出力：
% 三次元に変換・合成された折れ線 ：XX, YY, ZZ
% アルゴリズムは下記のとおり。
% ただし、A1,B1は、
% x軸平行断面折れ線の処理の場合は、セルAout1,Bout1の各要素のベクトル、
% y軸平行断面折れ線の処理の場合は、セルAout2,Bout2の各要素のベクトル。
% B0は、Z行列の全要素が0の場合のB1の値を表すベクトル。
% Kyは、行列Yの各行の要素の値（全て同一値）を示すスカラー。
% Kxは、行列Xの各列の 〃
% az,elは、視点の位置（既述）
%
% x軸平行断面折れ線の処理
% Zの値を0とすれば、
% A1=X*cos(az)+Ky*sin(az) .....(1)
% B0=(Ky*cos(az)-X*sin(az))*sin(el) .....(2)
% (1)式より、
% X=(A1-Ky*sin(az))/cos(az) .....(3)
% (3)式を(2)式に代入すると、
% B0=(Ky*cos(az)-(A1-Ky*sin(az))*sin(az)/cos(az))*sin(el)
% =-A1*tan(az)*sin(el) + Ky*(cos(az)+sin(az)*tan(az))*sin(el) .....(4)
% (4)式にA1を代入すればB0が求まる（tan(az)=infを避けるように、最初の方でazを操作済）。
% A1,B0の点での本来のZの値は
% Z=(B1-B0)/cos(el) .....(5)
% x,y平面上の各点の座標は、
% |a1| | cos(az) sin(az) | |x| |x|
% | | = | |*| | = Mc*| |
% |b0| |-sin(az)*sin(el) cos(az)*sin(el)| |y| |y|
% より、
% |x| |a1|
% | | = Mc^(-1)*| | .....(6)
% |y| |b0|
%
% y軸平行断面折れ線の処理
% Zの値を0とすれば、
% A1=Kx*cos(az)+Y*sin(az) .....(7)
% B0=(Y*cos(az)-Kx*sin(az))*sin(el) .....(8)
% (7)式より、
% Y=(A1-Kx*cos(az))/sin(az) .....(9)
% (9)式を(8)式に代入すると、
% B0=((A1-Kx*cos(az))*cos(az)/sin(az)-Kx*sin(az))*sin(el)
% =A1*cot(az)*sin(el) - Kx*(sin(az)+cos(az)*cot(az))*sin(el) .....(10)
% (5),(6)式についてはx軸平行断面折れ線の場合と同じ。
%
% sin,cos,cotが微小になったり、tanが極端に大きくなることもあり、計算精度が心配だったが、
% MATLAB演算の有効桁数が十分なためか、今のところ問題なし。
C3=cos(el);
Mc=[cos(az) sin(az) ; -sin(az)*sin(el) cos(az)*sin(el)];
Mcc=inv(Mc);
% x軸平行断面折れ線の処理
C1=-tan(az)*sin(el);
C2=(cos(az)+sin(az)*tan(az))*sin(el);
X1={};
Y1={};
Z1={};
for line=1:size(Aout1,2)
Ky=Y_v(line);
B0=C1*Aout1{line}+C2*Ky;
zz=(Bout1{line}-B0)/C3;
N=length(Aout1{line});
XY=[];
for n=1:N
xy=Mcc*[Aout1{line}(n) ; B0(n)];
XY=[XY xy];
end
X1{line}=XY(1,:);
Y1{line}=XY(2,:);
Z1{line}=zz;
end
% y軸平行断面折れ線の処理
C1=cot(az)*sin(el);
C2=-(sin(az)+cos(az)*cot(az))*sin(el);
X2={};
Y2={};
Z2={};
for line=1:size(Aout2,2)
Kx=X_v(line);
B0=C1*Aout2{line}+C2*Kx;
zz=(Bout2{line}-B0)/C3;
N=length(Aout2{line});
XY=[];
for n=1:N
xy=Mcc*[Aout2{line}(n) ; B0(n)];
XY=[XY xy];
end
X2{line}=XY(1,:);
Y2{line}=XY(2,:);
Z2{line}=zz;
end
XX=[];
YY=[];
ZZ=[];
for line=1:size(Aout1,2)
XX=[XX NaN X1{line}];
YY=[YY NaN Y1{line}];
ZZ=[ZZ NaN Z1{line}];
end
for line=1:size(Aout2,2)
XX=[XX NaN X2{line}];
YY=[YY NaN Y2{line}];
ZZ=[ZZ NaN Z2{line}];
end
% 峰の頂上や谷底が陰線となって消えたとき、Z軸のスケールが勝手に変わるのを防止するため、
% 本来のZ軸の上下端にダミーplotデータを加えておく。
% （この問題は、以下の処理をしなくても、疑似grid線も同時に描いてやれば解決できる。
% しかし、grid線を描かない応用もあることを考えれば、この処理はしておいた方が良い。）
Xmin=V_grid{1}(1);
Xmax=V_grid{1}(2);
Ymin=V_grid{3}(1);
Ymax=V_grid{3}(2);
Zmin=V_grid{5}(1);
Zmax=V_grid{5}(2);
switch quad
case 1
XX=[XX NaN Xmax NaN Xmax];
YY=[YY NaN Ymin NaN Ymin];
ZZ=[ZZ NaN Zmin NaN Zmax];
case 2
XX=[XX NaN Xmax NaN Xmax];
YY=[YY NaN Ymax NaN Ymax];
ZZ=[ZZ NaN Zmin NaN Zmax];
case 3
XX=[XX NaN Xmin NaN Xmin];
YY=[YY NaN Ymax NaN Ymax];
ZZ=[ZZ NaN Zmin NaN Zmax];
case 4
XX=[XX NaN Xmin NaN Xmin];
YY=[YY NaN Ymin NaN Ymin];
ZZ=[ZZ NaN Zmin NaN Zmax];
end
end % end of function 'restore_3D'（functionの'end'の後ろに';'をつけるとエラー）
% ■■■===============================================================
function[Xg,Yg,Zg]=quasi_grid(V_grid, az, el, quad, q_el, Ma1, Mb1, Ma2, Mb2)
% grid線の描画はユーザが制御できる自由度が小さく、plot3で描いた図上でgrid onすると、陰線処
% 理した図形の背後にあって本来は表示されるべきではないgrid線もそのまま表示されてしまう。
% この対策には、次の２通りが考えられる。
%
% ① 最外殻のマスク領域を、視点からx-y,y-z,z-x平面に投影して３つの多角形を作る。それらに囲
% まれた領域を、fill3コマンドによって不透明な白で塗り潰し、grid線の不要な部分を覆い隠す。
% ② 本来のgrid線は使用せず、陰線処理された疑似grid線を作り、plot3コマンドで描く。
%
% しかし、①の場合には、本来は覆い隠される筈のgrid線の一部が見えたままになり実用に耐えない。
% これは、MATLABのバグによるものと思われる。
% また、②の場合には、plot3による疑似grid線のスタイルが、正規のgrid線のような極細線や微小ピ
% ッチの破線にはならず、目立ち過ぎてやや不快。
%
% どちらも、十分に満足できる解決策ではないが、無難な方の②を採用する。線の色を'k'ではなく、
% 灰色の[0.85 0.85 0.85]程度とすれば、目立ち方は抑えられる。
% 入力：
% grid線用の座標情報 ：V_grid{1} （x軸の[下端 上端]）
% {2} （x軸の目盛を示すベクトル）
% {3},{4} （y軸について同上）
% {5},{6} （z軸について同上）
% 視点の方位角と伏角 ：az, el
% 視点azによる象限分け ：quad
% 視点elによる象限分け ：q_el
% x軸平行断面折れ線処理を
% 完了時の最広マスク ：Ma1, Mb1 （ {1}が上縁、{2}が下縁 ）
% y軸平行断面折れ線処理を
% 完了時の最広マスク ：Ma2, Mb2 （ 〃 ）
%
% 出力：
% 三次元に変換・合成されたgrid線 ：Xg, Yg, Zg
% 疑似grid線の処理
% x目盛
x=sort([V_grid{1} V_grid{2}]); % 目盛部だけでなく、上下端にもgrid線をひく
n=find( ([x inf]-[-inf x] )==0); % 重複する要素番号を拾い出す
x(n)=[]; % 重複要素を削除
% y目盛
y=sort([V_grid{3} V_grid{4}]);
n=find( ([y inf]-[-inf y] )==0); % 重複する要素番号を拾い出す
y(n)=[]; % 重複要素を削除
% z目盛
z=sort([V_grid{5} V_grid{6}]);
n=find( ([z inf]-[-inf z] )==0); % 重複する要素番号を拾い出す
z(n)=[]; % 重複要素を削除
switch quad % 象限に応じた座標面や軸位置の設定
case 1
p_yz=min(x); p_zx=min(y); % p_yz: 背景y_z面のx座標, p_zx: 背景z_x面のy座標
ax_x=max(y); ax_y=max(x); % ax_x: x軸位置のy座標, ax_y: y軸位置のx座標
% （axの'x'はx軸やx座標とは無関係。axはaxisの略）
case 2
p_yz=max(x); p_zx=min(y);
ax_x=max(y); ax_y=min(x);
case 3
p_yz=max(x); p_zx=max(y);
ax_x=min(y); ax_y=min(x);
case 4
p_yz=min(x); p_zx=max(y);
ax_x=min(y); ax_y=max(x);
end
if q_el==1 || q_el==2
p_xy=min(z); % p_xy: 背景x_y面のz座標
else % q_el==3 | q_el==4
p_xy=max(z);
end
% x-y平面上のオリジナル疑似grid線（x軸に平行なgrid線）
G{1,1}=repmat([min(x) max(x)],length(y),1); % x成分
G{1,2}=[y' y']; % y成分
G{1,3}=ones(length(y),2)*p_xy; % z成分
% x-y平面上のオリジナル疑似grid線（y軸に平行なgrid線）
G{2,1}=[x' x']; % x成分
G{2,2}=repmat([min(y) max(y)],length(x),1); % y成分
G{2,3}=ones(length(x),2)*p_xy; % z成分
% y-z平面上のオリジナル疑似grid線（y軸に平行なgrid線）
G{3,1}=ones(length(z),2)*p_yz; % x成分
G{3,2}=repmat([min(y) max(y)],length(z),1); % y成分
G{3,3}=[z' z']; % z成分
% y-z平面上のオリジナル疑似grid線（z軸に平行なgrid線）
G{4,1}=ones(length(y),2)*p_yz; % x成分
G{4,2}=[y' y']; % y成分
G{4,3}=repmat([min(z) max(z)],length(y),1); % z成分
% z-x平面上のオリジナル疑似grid線（z軸に平行なgrid線）
G{5,1}=[x' x']; % x成分
G{5,2}=ones(length(x),2)*p_zx; % y成分
G{5,3}=repmat([min(z) max(z)],length(x),1); % z成分
% z-x平面上のオリジナル疑似grid線（x軸に平行なgrid線）
G{6,1}=repmat([min(x) max(x)],length(z),1); % x成分
G{6,2}=ones(length(z),2)*p_zx; % y成分
G{6,3}=[z' z']; % z成分
% x,y,zの各軸の目盛線と一致する疑似grid線は描かないようにする（重なると見苦しいので）。
% そのために、まずは、MATLAB標準の軸の目盛線の位置を象限別に把握しておく。
% なお、伏角elの象限q_elが2と3の場合には、方位角azは実際に指定した方向ではなく、180度反転し
% た方向（quad=1,2,3,4 → quad=3,4,1,2）に切り替わっていることに留意する必要がある。
% （蛇足： 規則性を十分に探れば、もっとスマートなケース分けができるのかもしれないが、面倒
% なので力技での遣っ付け仕事。まことに冗長。）
switch quad
case 1
switch q_el
case 1
ax_x=[max(y) min(z)];
ax_y=[max(x) min(z)];
ax_z=[max(x) min(y)];
case 2
ax_x=[min(y) max(z)];
ax_y=[min(x) max(z)];
ax_z=[min(x) max(y)];
case 3
ax_x=[max(y) max(z)];
ax_y=[max(x) max(z)];
ax_z=[min(x) max(y)];
case 4
ax_x=[min(y) min(z)];
ax_y=[min(x) min(z)];
ax_z=[max(x) min(y)];
end
case 2
switch q_el
case 1
ax_x=[max(y) min(z)];
ax_y=[min(x) min(z)];
ax_z=[max(x) max(y)];
case 2
ax_x=[min(y) max(z)];
ax_y=[max(x) max(z)];
ax_z=[min(x) min(y)];
case 3
ax_x=[max(y) max(z)];
ax_y=[min(x) max(z)];
ax_z=[min(x) min(y)];
case 4
ax_x=[min(y) min(z)];
ax_y=[max(x) min(z)];
ax_z=[max(x) max(y)];
end
case 3
switch q_el
case 1
ax_x=[min(y) min(z)];
ax_y=[min(x) min(z)];
ax_z=[min(x) max(y)];
case 2
ax_x=[max(y) max(z)];
ax_y=[max(x) max(z)];
ax_z=[max(x) min(y)];
case 3
ax_x=[min(y) max(z)];
ax_y=[min(x) max(z)];
ax_z=[max(x) min(y)];
case 4
ax_x=[max(y) min(z)];
ax_y=[max(x) min(z)];
ax_z=[min(x) max(y)];
end
case 4
switch q_el
case 1
ax_x=[min(y) min(z)];
ax_y=[max(x) min(z)];
ax_z=[min(x) min(y)];
case 2
ax_x=[max(y) max(z)];
ax_y=[min(x) max(z)];
ax_z=[max(x) max(y)];
case 3
ax_x=[min(y) max(z)];
ax_y=[max(x) max(z)];
ax_z=[max(x) max(y)];
case 4
ax_x=[max(y) min(z)];
ax_y=[min(x) min(z)];
ax_z=[min(x) min(y)];
end
end
% x軸の目盛線に一致する疑似grid線の消去
for n=[1 6]
m=find( G{n,2}(:,1)==ax_x(1) & G{n,3}(:,1)==ax_x(2) );
G{n,1}(m,:)=[];
G{n,2}(m,:)=[];
G{n,3}(m,:)=[];
end
% y軸の目盛線に一致する疑似grid線の消去
for n=[2 3]
m=find( G{n,1}(:,1)==ax_y(1) & G{n,3}(:,1)==ax_y(2) );
G{n,1}(m,:)=[];
G{n,2}(m,:)=[];
G{n,3}(m,:)=[];
end
% z軸の目盛線に一致する疑似grid線の消去
for n=[4 5]
m=find( G{n,1}(:,1)==ax_z(1) & G{n,2}(:,1)==ax_z(2) );
G{n,1}(m,:)=[];
G{n,2}(m,:)=[];
G{n,3}(m,:)=[];
end
B_range=[min([Mb1{1} Mb1{2}]) max([Mb1{1} Mb1{2}])];
% 最広マスク領域の上縁の曲線を求める
% 関数呼び出し ================================================
[A_ceil,B_ceil]=mask_renew(Ma1{1},Mb1{1},Ma2{1},Mb2{1},1,B_range);
% =============================================================
% 最広マスク領域の下縁の曲線を求める
% 関数呼び出し ================================================
[A_floor,B_floor]=mask_renew(Ma1{2},Mb1{2},Ma2{2},Mb2{2},2,B_range);
% =============================================================
% 計算誤差によって、mesh曲面の縁の一部が、疑似grid線で上書きされてしまうことがある。
% この場合、疑似grid線の色を灰色（目立ち過ぎずに丁度良い色）に設定していると、mesh曲面の一部
% が消えたように見える。これは、mesh曲面と疑似grid線を描く順序を入れ替えても改善されない。
% ただし、zoom inすると見えるようになることもあるし、さらに殆どの場合、印刷すれば正常に描画
% される。
% しかし、印刷されるとおりに表示されないのでは困るので、この状況をできるだけ避けるために、
% マスク領域を上下に僅かだけ広げてやる（これでも、稀には好ましくない表示がみられるが、印刷
% 結果が正常なら良しとしたい）。
Bd0=(B_range(2)-B_range(1))/1000; % マスクを上下に広げる量
B_ceil =B_ceil +Bd0;
B_floor=B_floor-Bd0;
% また、az=90*Nのときには、mesh曲面の左右の縁が表示も印刷もされなくなることがある。
% これを避けるつもりで、マスクの左右を延長してみた。しかし、残念ながらこれでは全く解決しな
% かった。しかし、後日の参考用として消去しないでコメント化したまま残す。
%
% Ad0=(A_ceil(length(A_ceil))-A_ceil(1))/1000
% A_ceil=[A_ceil(1)-Ad0 A_ceil A_ceil(length(A_ceil))+Ad0];
% B_ceil=[B_ceil(1) B_ceil B_ceil(length(B_ceil)) ];
% A_floor=[A_floor(1)-Ad0 A_floor A_floor(length(A_floor))+Ad0];
% B_floor=[B_floor(1) B_floor B_floor(length(B_floor)) ];
% 疑似grid線を視線と垂直な面（直交a,b座標）に投影し、陰線消去。
for n=1:6
X=G{n,1};
Y=G{n,2};
Z=G{n,3};
% 投影
Ba{n} = X*cos(az) + Y*sin(az);
Bb{n} = sin(el)*( Y*cos(az) - X*sin(az) ) + Z*cos(el);
if Ba{n}(1,1)>Ba{n}(1,2) % 視点によつてBa{n}の並びの方向が変わるので、昇順に統一。
Ba{n}=fliplr(Ba{n});
Bb{n}=fliplr(Bb{n});
end
% 陰線消去
if ( n==1 || n==2 || n==3 || n==6 ) % 疑似gridがz軸に平行でない場合
for m=1:size(Ba{n},1)
AB{n,m}=[];
% 関数呼び出し ================================================
[A1,B1]=erase_hidden(A_ceil,B_ceil,Ba{n}(m,:),Bb{n}(m,:),1);
[A2,B2]=erase_hidden(A_floor,B_floor,Ba{n}(m,:),Bb{n}(m,:),2);
% =============================================================
AB{n,m}=[AB{n,m};[A1' B1';NaN NaN;A2' B2']];
end
else % 疑似gridがz軸に平行な場合（ n==4 || n==5 ）
for m=1:size(Ba{n},1)
AB{n,m}=[];
A_cross=Ba{n}(m,1);
A1=[A_cross A_cross];
B_upper=interp1(A_ceil,B_ceil,A_cross); % マスク上縁との交点
B_lower=interp1(A_floor,B_floor,A_cross); %
E_upper=max(Bb{n}(m,:)); % upper-end of original grid-line
E_lower=min(Bb{n}(m,:)); % lower-end of original grid-line
if (E_lower<B_lower && B_upper<E_upper)
B1=[B_upper E_upper];
B2=[E_lower B_lower];
AB{n,m}=[AB{n,m};[A1' B1';NaN NaN;A1' B2']];
elseif B_upper<E_upper
B1=[B_upper E_upper];
AB{n,m}=[AB{n,m};[A1' B1']];
elseif E_lower<B_lower
B2=[E_lower B_lower];
AB{n,m}=[AB{n,m};[A1' B2']];
else
% 見える疑似grid線は無い
end
end
end
end
% 二次元の疑似grid線を三次元データに復元
Xg=[];
Yg=[];
Zg=[];
for n=1:6
if n==1 || n==2 % x-y平面上の疑似grid線
C0=p_xy;
elseif n==3 || n==4 % y-z平面上の疑似grid線
C0=p_yz;
else % n==5 || n==6 z-x平面上の疑似grid線
C0=p_zx;
end
for m=1:size(AB,2)
if ~isempty(AB{n,m})
AB_nm=AB{n,m};
% 関数呼び出し ================================================
[X1,Y1,Z1]=grid_2D_to_3D(n, AB_nm, C0, az, el);
% =============================================================
if ~isempty(X1) || ~isempty(Y1) || ~isempty(Z1)
Xg=[Xg NaN X1 NaN];
Yg=[Yg NaN Y1 NaN];
Zg=[Zg NaN Z1 NaN];
end
end
end
end
end % end of function 'quasi_grid'（functionの'end'の後ろに';'をつけるとエラー）
% ■■■===============================================================
function[X1,Y1,Z1]=grid_2D_to_3D(n, AB_nm, C0, az, el)
% 二次元の疑似grid線を三次元データに復元
% 入力：
% 疑似grid線が描かれる平面と向き
% による分類番号 ：n n=1: x-y平面上のx軸に平行な疑似grid線
% n=2: 〃 y軸 〃
% n=3: y-z平面上のy軸に平行な疑似grid線
% n=4: 〃 z軸 〃
% n=5: z-x平面上のz軸に平行な疑似grid線
% n=6: 〃 x軸 〃
% 二次元疑似grid線のa,b座標の行列 ：AB_nm
% 疑似grid線が描かれる平面の位置
% （x,y,z座標いずれかの一定値） ：C0
% 視点の方位角と伏角 ：az, el
%
% 出力：
% 三次元に変換・合成された疑似grid線 ：X1, Y1, Z1
% function 'restore_3D' と類似の処理だが、各grid線は各座標平面に張り付いているので、もっと
% 簡単な処理が可能。
% 下記は既出の基本式を変形したものだが、このままでは3変数X,Y,Zが2変数A,Bに置き換えられてし
% まっているので、A,Bを単純にX,Y,Zに逆変換することはできない。しかし、grid線の場合は、X,Y,Z
% のいずれかが既知の定数なので、簡単に逆変換できる。
%
% A = X*cos(az) + Y*sin(az)
% B = -X*sin(az)*sin(el) + Y*cos(az)*sin(el) + Z*cos(el)
%
% x-y平面のgrid線（n=1,2）は、Zが既知の定数C0だから、
%
% |X| | cos(az) sin(az) | | A |
% | |=inv | | * | |
% |Y| |-sin(az)*sin(el) cos(az)*sin(el)| |B - C0*cos(el)|
%
% y-z平面のgrid線（n=3,4）は、Xが既知の定数C0だから、
%
% |Y| | sin(az) 0 | | A - C0*cos(az) |
% | |=inv | | * | |
% |Z| |cos(az)*sin(el) cos(el)| |B + C0*sin(az)*sin(el)|
%
% z-x平面のgrid線（n=5,6）は、Yが既知の定数C0だから、
%
% |Z| | 0 cos(az) | | A - C0*sin(az) |
% | |=inv | | * | |
% |X| |cos(el) -sin(az)*sin(el)| |B - C0*cos(az)*sin(el)|
%
X1=[];
Y1=[];
Z1=[];
if n==1 || n==2
Cinv=inv([cos(az) sin(az);-sin(az)*sin(el) cos(az)*sin(el)]);
AB_line=[AB_nm(:,1) AB_nm(:,2)-C0*cos(el)]';
XY=Cinv*AB_line;
X1=XY(1,:); Y1=XY(2,:); Z1=C0*ones(1,length(X1));
j=find(isnan(X1) | isnan(Y1));
X1(j)=NaN; Y1(j)=NaN; Z1(j)=NaN;
elseif n==3 || n==4
Cinv=inv([sin(az) 0;cos(az)*sin(el) cos(el)]);
AB_line=[AB_nm(:,1)-C0*cos(az) AB_nm(:,2)+C0*sin(az)*sin(el)]';
YZ=Cinv*AB_line;
Y1=YZ(1,:); Z1=YZ(2,:); X1=C0*ones(1,length(Y1));
j=find(isnan(Y1) | isnan(Z1));
X1(j)=NaN; Y1(j)=NaN; Z1(j)=NaN;
elseif n==5 || n==6
Cinv=inv([0 cos(az);cos(el) -sin(az)*sin(el)]);
AB_line=[AB_nm(:,1)-C0*sin(az) AB_nm(:,2)-C0*cos(az)*sin(el)]';
ZX=Cinv*AB_line;
Z1=ZX(1,:); X1=ZX(2,:); Y1=C0*ones(1,length(Z1));
j=find(isnan(Z1) | isnan(X1));
X1(j)=NaN; Y1(j)=NaN; Z1(j)=NaN;
end
end % end of function 'grid_2D_to_3D'（functionの'end'の後ろに';'をつけるとエラー）
% ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
% 以下には、本functionの使用例を示すscriptをコメント化して掲載しています。
% '% % test_for_mesh2plot3.m' 以下の全行をコピーし、全行頭の'% 'の2文字を削除してから、
% 適当なmファイル名で保存して実行させてみてください（Editorの文字のEncodeの種類は、MATLABの
% コマンドラインから'slCharacterEncoding()'として調べ、それに合わせます）。
% ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
% % test_for_mesh2plot3.m
%
% % mesh図と疑似mesh図を左右で対比。視点を変更しながら繰り返す。
%
% clear
% close all
%
% % === モデル選択のための指示待ち ===
% % ここでinputコマンドを使用しておかないと、何故か、初回だけ、figure(1)画面がMATLABウィンド
% % ウの裏に隠れてしまうので、2枚のfigureを並べて表示できない。
%
% a=input('モデルの選択：通常のpeak曲面は"1"、変形peak曲面は"2" を入力してください。\n','s');
%
% [X,Y,Z]=peaks; % モデル曲面
% if a=='1'
% % そのまま
% else
% % バグ出しのため、素直なデータを避ける。
% X(:,[1:12 38:49])=[]; % データ定義領域を正方形でなく長方形にする。
% Y(:,[1:12 38:49])=[];
% Z(:,[1:12 38:49])=[];
% X=X-1; % X軸をずらせてみる。
% end
%
% az=-37.5; el=30; % デフォルトの視点
% %az=0; el=0; % 必要に応じて'%'を外し、視点をデフォルト以外の値に変える
%
% while true % 連続繰り返し。中断は'ctrl-c'による。
%
% figure(1);
% mesh(X,Y,Z,'EdgeColor','blue'); % お手本となるmesh図を描かせる。
% %set(gcf,'Position',[200 350 800 600]); % 画面サイズ1920*1080用の位置決め（要・適宜変更）
% set(gcf,'Position',[200 350 800 600]/1.5); % 画面サイズ1920*1080用の位置決め（要・適宜変更）
% xlabel('X');
% ylabel('Y');
% zlabel('Z');
% title('meshコマンドによる作図','FontSize',18);
% %axis equal; % 各軸長を比例尺関係にするときは'%'を外す。
% view(az,el); % 視点の設定
% % これ以外にも、好みに応じて各種propertyの設定を行う。
%
% % お気に入りのmesh図が完成したら、function呼び出し用のargumentを準備する。
%
% V_view=[az,el]; % 視点
% V_grid=get(gca,{'XLim' 'XTick' 'YLim' 'YTick' 'ZLim' 'ZTick'});
% % grid線用の情報
% data_aspect=daspect; % x,y,z軸のアスペクト比
%
% % 視点により軸のスケールが変化しないように固定（疑似mesh図と条件を合わせて比較するため）
% set(gca,'XLim',V_grid{1},'XTick',V_grid{2},'YLim',V_grid{3},'YTick',V_grid{4},'ZLim',...
% V_grid{5},'ZTick',V_grid{6}); % 各軸の尺度の設定
%
% % 関数呼び出し ================================================
% [X1,Y1,Z1,Xg,Yg,Zg]=mesh2plot3(X,Y,Z,V_view,V_grid,data_aspect);
% % =============================================================
%
% figure(2);
% %set(gcf,'Position',[1050 350 800 600]); % 画面サイズ1920*1080用の位置決め（要・適宜変更）
% set(gcf,'Position',[1050 350 800 600]/1.5); % 画面サイズ1920*1080用の位置決め（要・適宜変更）
% plot3(X1,Y1,Z1,'b'); % 疑似mesh曲面の描画
% hold on
% plot3(Xg,Yg,Zg,'-','Color',[0.85 0.85 0.85]); % 疑似グリッド線のプロツト
% xlabel('X');
% ylabel('Y');
% zlabel('Z');
% title('plot3コマンドによる疑似mesh','FontSize',18);
% % 以下、お手本のmesh図との様式合わせ
% %axis equal; % 各軸長を比例尺関係にするときは'%'を外す。
% view(V_view); % 視点の設定
% set(gca,'XLim',V_grid{1},'XTick',V_grid{2},'YLim',V_grid{3},'YTick',V_grid{4},'ZLim',...
% V_grid{5},'ZTick',V_grid{6}); % 各軸の尺度の設定
%
% % === 方位角azと伏角elを増減するための指示待ち ===
%
% disp(['現在の視点は az=' num2str(az) '[deg], el=' num2str(el) '[deg] です。' ...
% ' 実行の中断は"ctrl-c"、'])
% a=input('方位角azの増/減には "s/a"、 伏角elの増/減には "w/z" を入力してください。\n','s');
% if a=='a' % Decrease az
% az=az-10;
% elseif a=='s' % Increase az
% az=az+10;
% elseif a=='z' % Decrease el
% el=el-10;
% elseif a=='w' % Increase el
% el=el+10;
% end
%
% close all
%
% end % End of "while true"