【MATLAB】曲面上を指し示す「理想的なマーカー」を求めて試行錯誤

１．はじめに

　三次元のグラフにおいて、時々、曲面上の特定の点に注意を向けさせたいことがあります。そのようなときに最も効果的な方法は、その部分をマーカーで指し示すことです。しかし、二次元グラフの場合とは異なり、plot3(x,y,z,'o') のような単純な方法では、期待どおりの効果が得られない場合もあります。

　例えば、図１は、最もシンプルな方法で、球面上の多数の点にマーカーを重ね描きしたものです。見てのとおり、全身が正しく表示できているのは、真正面に描かれているマーカーだけです。これではなかなか満足はできません。

　図２は、最も簡単に思いつく対策を講じて描いた図です。視線に沿ってマーカーを手前の方に引き寄せたものです。視線方向に移動しているだけなので、このグラフ上だけで見れば、位置ズレを感じさせずにマーカーの全身を表示できるようになりました。

　しかし、実際にはその操作だけでは不十分です。引き寄せ操作によって、最適サイズだった座標圏内からマーカーがはみ出てしまうため、それをカバーしようとして座標範囲が勝手に拡大します。その結果、主役の球体の表示が相対的に小さくなり、訴求力が失われてしまいます。これを防ぐために xlim, ylim, zlim コマンドによって座標範囲の勝手な拡大を禁止する必要があります。しかし、こうすると、座標圏内から飛び出したマーカーは画面から消えてしまいます。この対策として、マーカーの 'Clipping' 属性を 'off' に変更してやる必要があります。また、手前に引き寄せる量が大きすぎると、本来は見えないはずの球体の裏側にあるマーカーまでシャシャリ出てきて大混乱になります。これを防ぐための裏側マーカーの削除処理も必要になります。ここまで手を加えても、マウスでグリグリとドラッグしながら見る角度を変えると、マーカーはとんでもない位置に移動し、舞台裏が丸見えになってしまいます。どう見ても「理想的なマーカー」とは言い難いものです。

　「見る角度を変える操作」をしたときに割り込み（コールバック）を発生させて、マーカーを正規の位置に描き変えてやれば舞台裏は見られなくても済むはずです。しかし、いくら捜し回ってみても、「見る角度を変える操作」によるコールバック機能はないようです。

２．特製マーカー

　上記のような簡単な方法でも、問題なく利用できる場面も多いとは思われます。しかし、見る角度を変えても正規の位置を保てるようなマーカーが欲しいとなると、ある程度の工夫が必要になります。

　この工夫のために、試行錯誤を繰り返して時間の浪費が続きました。本来であれば、この間の苦労話を披露して共感を求めたいところです。しかし、終わってみれば「コロンブスの卵」的なことが多く、読み手にとって面白い話にはなりそうもありません。残念ではありますが、現段階での最終結果を示すだけで済ますことにします。

２.１　[真球面] & [全軸 等スケール] 専用マーカー

　図３は、特製のマーカーによるマーキングの結果です。指定された位置にシールのように貼り付くので、どこから見ても位置ずれはありません。なお、この図のように、マーカーの形が真円の場合には気になりませんが、上下左右が非対称のマーカーの場合、その向きが位置によってバラバラというのでは美しくありません。そのため、図４に示すように、マーカーの頂部が常に最大の登り傾斜の方を向いて配置されるようにも考慮しています。

　これを実現するのは次の２つのローカル関数です。詳細は本記事の末尾に添付している「プログラム」の項を参照ください。２つの関数はプログラムリストの最後の部分に記載しています。本項の【 [真球面] & [全軸 等スケール] 専用マーカー 】としてだけではなく、次項の【 [任意曲面] & [各軸 異スケール] 対応マーカー 】にも共通して使えるように記述されています。

• function marker
  　「曲面上の貼付け位置、そこの法線ベクトル、マーカーのサイズ、その形」を入力すると、それに見合ったマーカー画像を形作る三次元折線データを出力します。この中では、下記のローカル関数 rot_points_3d を呼び出して利用しています。
  　この他に、次項での利用のために、「各軸のスケール幅、マーカーサイズの基準軸名」の入力もできるようにしています。
• function rot_points_3d
  　「回転軸のベクトルの始点、そのベクトルの向き、回転される点群の座標、回転角度」を入力すると、回転された点群の座標を出力します。この中では「ロドリゲスの三次元回転行列」を使っています。

２.２　[任意曲面] & [各軸 異スケール] 対応マーカー

　図３,４は真球面のため、都合の良いことが多く、特製マーカーのプログラムもそれほど複雑にはなりません。例えば、マーカーの位置を指定する数値が、そのまま法線ベクトルとして利用できるという大きな利点があります。また、全面が凸面であることも好都合です。あとから分かるように、凹面は貼付けマーカーにとっては鬼門です。

　その一方、球面は「有効なMESHGRID」（碁盤目状の行儀の良いXY格子）による曲面ではないため、曲面の補間（内挿）コマンド interp2 が使えないという欠点があります。この欠点は、真球に限らず、媒介変数を使って描いた複雑な曲面に共通の問題です。我がアイコンの「クラインの壺」などは、真球のような長所は一つも無いのに、この欠点だけはあるため、マーカーにとってはかなり難しい貼付け対象になります。まだ実現できてはいません。

　ここで言う「任意曲面」とは、真の「任意」ではなく、「有効なMESHGRID」の X,Y を元に描かれた Z の値が一価だけの曲面のことを指します。そのため、interp2 コマンドが利用でき、複雑なアルゴリズムでも楽にプログラミングできるようになります。また、ここでは、「axis equal」ではない、各軸がバラバラのスケールの図も扱えるようにするという課題も追加しました。このような環境でもマーカーの形が崩れないようにするため、より難しい対応が必要になります。

　この方針に沿って、MATLAB でお馴染みの peaks 曲面をモデルにして、山や谷の色々な場所に貼り付けられたマーカーを図５,６に示しています。これらの図では、ローカル関数 marker から受け取ったデータにさらに手を加え、補間処理によって、マーカーを曲面に圧着しています。この処理をしないと、起伏が激しい曲面の凹んだ部分では、マーカーの一部が曲面に隠れてしまうことがあるからです。ただし、貼付け位置によっては、この処理をしない図７,８の方が自然に見えることもあります。

３．おわりに

　この最終結果も、未だ「理想的なマーカー」とは言えません。圧着されたマーカーの不自然な変形は見るに堪えません。これからも試行錯誤は必要なようです。

　しかし、もう少々飽きてきました(^-^;。　axis equal の真球面に〇印を貼付けられるようになった段階で、既に私の喫緊の要求は満たされていたのです。それなのに、余計な領域まで深入りし過ぎてしまったようです。

４．プログラム

marker_on_surf01.m

% marker_on_surf01.m

% 「曲面に貼り付けるマーカー」を生成するローカル関数の動作確認。

clear
close all

% ■■■■■■■■■■　地球の描画。
figure(1);
Re=1;                                      % 地球の半径。
[X,Y,Z] = sphere(50);                      % 50×50メッシュの地球
surf(X*Re,Y*Re,Z*Re,'EdgeColor','#aaa');   % 描画。線色は淡く。
colormap(gca,parula*0.3+0.7);              % 表面色も淡くする。
axis equal;
xlabel('X-axis'); ylabel('Y-axis'); zlabel('Z-axis');
hold on;

% 図の複製。
h=get(gcf,'Children');
figure(2);
copyobj(h,gcf);
figure(3);
copyobj(h,gcf);
figure(4);
copyobj(h,gcf);

% 座標軸（右手系直交デカルト座標）をx,y,zとする。
% 　原点は地球の中心。
% 　地軸をz軸（正の向きは北極側）。
% 　z軸に直交する経度0°の向きをx軸。
% 　x軸をz軸周りに右ネジ方向に90°回してy軸。

% マーカーを貼り付ける位置の座標群。
% 緯度（北緯が正）: [-90:30:90]
% 経度（東経が正）: 
% 　　　緯度 = 0 のとき　→　[-180+30:30:180]+16
% 　　　緯度 = ±30 のとき　→　[-180+36:36:180]+16
% 　　　緯度 = ±60 のとき　→　[-180+60:60:180]+16
% 　　　緯度 = ±90 のとき　→　[0]+16

a=[-90:30:90];                % 緯度
b{1}=[-180+30:30:180]+16;     % 経度
b{2}=[-180+36:36:180]+16;
b{3}=[-180+60:60:180]+16;
b{4}=[0]+16;

points=[];                    % 緯度・経度のリスト

% マーカーを貼り付ける点の緯度・経度のリスト(points)を作成（多行２列）
for n=1:length(a)
  lat=a(n);                   % 緯度
  if lat==0
    nb=1;
  elseif lat==30 || lat==-30
    nb=2;
  elseif lat==60 || lat==-60
    nb=3;
  elseif lat==90 || lat==-90
    nb=4;
  end
  n_lon=length(b{nb});        % 経度の個数
  for m=1:n_lon
    lon=b{nb}(m);             % 経度
    points=[points;lat lon];  % [緯度 経度]の組を縦に並べる。
  end
end
points=points*pi/180;         % deg → rad

% マーカー貼付け点の緯度・経度のリストを直交座標データに変換。
px=Re*cos(points(:,1)).*cos(points(:,2));
py=Re*cos(points(:,1)).*sin(points(:,2));
pz=Re*sin(points(:,1));
pp=[px py pz];                % マーカー貼付け点のx,y,z座標（多行３列）

% ==========
figure(1);  % MATLAB標準のマーカー（小細工なし）

for p=1:size(pp,1)
  plot3(pp(p,1),pp(p,2),pp(p,3),'sb','MarkerSize',18,'LineWidth',2.5)
end
title('図１　　標準のマーカー''s''（小細工なし）','FontSize',14)

% ==========
figure(2);  % MATLAB標準のマーカー（視線方向に引き寄せ）

% 禁じ手を使うため、マーカーが座標圏外にはみ出し、
% 　軸スケールが勝手に広がるので、それを抑制する。
xlim([-1 1]);
ylim([-1 1]);
zlim([-1 1]);

% 視線ベクトル[nx ny nz]の計算（自分に向かう方向が正）
[AZ,EL] = view;                  % figureの視点の取得
AZ=AZ*pi/180; EL=EL*pi/180;
nx=cos(EL)*cos(AZ-pi/2);
ny=cos(EL)*sin(AZ-pi/2);
nz=sin(EL);
K=0.5;                           % 適度な引き寄せ量

% マーカーを視線方向に適度に引き寄せ、球面への食い込みを外す。
for p=1:size(pp,1)
  if dot(pp(p,:),[nx ny nz])>0   % 球の表側のマーカーだけ表示する。
    hp=plot3(pp(p,1)+K*nx, ...
             pp(p,2)+K*ny, ...
             pp(p,3)+K*nz,'sb','MarkerSize',18,'LineWidth',2.5);
    set(hp,'Clipping','off');    % 手前に飛び出した座標圏外のマーカーも
                                 % 　クリップせずに表示させる。
  end
end
title('図２　　標準のマーカー''s''（視線方向に引き寄せ）','FontSize',14)

% 完成図をグリグリと回して、
% 　横から覗かれるとボロが出るので、卑怯にもこれを禁止！
ax=gca;
ax.Interactions = [rotateInteraction];
disableDefaultInteractivity(ax)

% ==========
figure(3);  % 特製〇印マーカー。

for p=1:size(pp,1)
  % マーカーを形成する三次元折線データを取得。
  % 　このプログラムの後の方にあるローカル関数の呼び出し。
  A=[pp(p,1) pp(p,2) pp(p,3)];    % 貼り付ける位置の座標。
  [xyz]=marker(A,A,0.1,1);
      % 第1引数: マーカー貼付け点の座標
      % 第2引数: マーカー貼付け点での曲面の法線ベクトル
      % 第3引数: マーカーの大きさ
      % 第4引数: マーカーの形（ 1: 〇、2: △ ）

  % マーカーの描画。
  plot3(xyz(:,1),xyz(:,2),xyz(:,3),'b','LineWidth',2)
end
title('図３　　特製マーカー〇','FontSize',14)

% ==========
figure(4);  % 特製△印マーカー。

for p=1:size(pp,1)
  A=[pp(p,1) pp(p,2) pp(p,3)];
  [xyz]=marker(A,A,0.12,2);
  plot3(xyz(:,1),xyz(:,2),xyz(:,3),'b','LineWidth',2)
end
title('図４　　特製マーカー△','FontSize',14)

% ■■■■■■■■■■　変形peaks曲面の描画。
% peaks曲面データを取得してサイズを変更
% （意地悪テストのため、axis equal には不向きなサイズにする）
figure(5);
[X,Y,Z] = peaks(50);              % 50×50メッシュの峰々のデータ
X=X/3;                            % x軸のサイズを変更。
surf(X,Y,Z,'EdgeColor','#aaa');   % 描画。線色は淡く。
colormap(gca,parula*0.3+0.7);     % 表面色も淡くする。
xlabel('X-axis'); ylabel('Y-axis'); zlabel('Z-axis');
hold on;

% 図の複製。
h=get(gcf,'Children');
figure(6);
copyobj(h,gcf);
figure(7);
copyobj(h,gcf);
figure(8);
copyobj(h,gcf);

% 曲面の法線ベクトルのデータを取得。
[Nx,Ny,Nz]=surfnorm(X,Y,Z);

% マーカーを貼り付ける場所の設定と、その場所の法線ベクトルを取得。
pp=[];     % マーカーの貼付け場所のx,y,z座標の組み合わせ。
nn=[];     % その場所の曲面の法線ベクトル。
for n=[-0.75:0.25:0.9]                % 貼付け場所のx座標群
  for m=[-0.75:0.25:0.9]*3            % 貼付け場所のy座標群
    pp=[pp; n m interp2(X,Y,Z,n,m)];  % 補間してz座標値を加える
    nn=[nn; interp2(X,Y,Nx,n,m) ...
            interp2(X,Y,Ny,n,m) ...
            interp2(X,Y,Nz,n,m)];     % 法線ベクトルの各成分も補間。
                                      % 　（単位化は少し崩れる）
  end
end

% 各軸の定義範囲を取得
a=gca;
scale=[a.XLim(2)-a.XLim(1) ...
       a.YLim(2)-a.YLim(1) ...
       a.ZLim(2)-a.ZLim(1)];

% ==========
figure(5)       % 球面以外の曲面（peaks曲面）への貼付け（△圧着済み）。

for p=1:size(pp,1)
  A=[pp(p,1) pp(p,2) pp(p,3)];     % マーカー位置。
  B=[nn(p,1) nn(p,2) nn(p,3)];     % 法線。まだ単位化されていないが、
                                   % 　次のmarker関数内で単位化される。

  [xyz]=marker(A,B,0.3,2,scale,2);               % マーカーの取得
      % 第1引数: マーカー貼付け点の座標
      % 第2引数: マーカー貼付け点での曲面の法線ベクトル
      % 第3引数: マーカーの大きさ
      % 第4引数: マーカーの形（ 1: 〇、2: △ ）
      % 第5引数: 各軸の定義範囲 [ xの範囲  yの範囲  zの範囲 ]
      % 第6引数: マーカーの大きさの基準軸（ 1: x軸、2: y軸、3: z軸 ）

  xyz(:,3)=interp2(X,Y,Z,xyz(:,1),xyz(:,2));     % マーカーの圧着処理
      % 曲面の起伏が激しいと、特に凹面部でマーカーが隠れるので、
      % 　マーカー折れ線の各点のz方向位置を曲面上に引き込む。
      % 　しかし、マーカー位置によっては形が大きく歪むことあり。

  plot3(xyz(:,1),xyz(:,2),xyz(:,3),'b','LineWidth',1)
end
title('図５　　球面以外の曲面への貼付け（△）圧着済み','FontSize',14)

% ==========
figure(6);  % 球面以外の曲面（peaks曲面）への貼付け（〇圧着済み）。

for p=1:size(pp,1)
  A=[pp(p,1) pp(p,2) pp(p,3)];                   % マーカー位置。
  B=[nn(p,1) nn(p,2) nn(p,3)];                   % 法線ベクトル。
  [xyz]=marker(A,B,0.23,1,scale,2);              % マーカーの取得
  xyz(:,3)=interp2(X,Y,Z,xyz(:,1),xyz(:,2));     % マーカーの圧着処理
  plot3(xyz(:,1),xyz(:,2),xyz(:,3),'b','LineWidth',1)
end
title('図６　　球面以外の曲面への貼付け（〇）圧着済み','FontSize',14)

% ==========
figure(7);  % 球面以外の曲面（peaks曲面）への貼付け（△未圧着）。

for p=1:size(pp,1)
  A=[pp(p,1) pp(p,2) pp(p,3)];                   % マーカー位置。
  B=[nn(p,1) nn(p,2) nn(p,3)];                   % 法線ベクトル。
  [xyz]=marker(A,B,0.3,2,scale,2);               % マーカーの取得
  plot3(xyz(:,1),xyz(:,2),xyz(:,3),'b','LineWidth',1)
end
title('図７　　球面以外の曲面への貼付け（△）未圧着','FontSize',14)

% ==========
figure(8);  % 球面以外の曲面（peaks曲面）への貼付け（〇未圧着）。

for p=1:size(pp,1)
  A=[pp(p,1) pp(p,2) pp(p,3)];                   % マーカー位置。
  B=[nn(p,1) nn(p,2) nn(p,3)];                   % 法線ベクトル。
  [xyz]=marker(A,B,0.23,1,scale,2);              % マーカーの取得
  plot3(xyz(:,1),xyz(:,2),xyz(:,3),'b','LineWidth',1)
end
title('図８　　球面以外の曲面への貼付け（〇）未圧着','FontSize',14)

% ==========================================
% ■■■ ローカル関数 marker　■■■
% ==========================================
% 曲面に貼り付けるマーカー（三次元折線データ）を生成する関数

function [xyz]=marker(p0,n0,rc,mark,scale,axis0)

  % 例えば、MATLAB標準マーカーの「plot3(...'o'...)」では、
  % 　〇の一部が曲面の背後に隠れて半円にしか見えないので、その対策。

  %【入力（必須）】
  % p0:    曲面上の貼り付け位置（x,y,z座標値の行ベクトル）
  % n0:    貼付け面の外向き法線ベクトル（x,y,z座標値の行ベクトル）
  %        　（単位化されていなくても可）
  % rc:    マーカーの大きさ
  % mark:  マーカーの形（　1: 〇印、　2: 変形三角形　）
  % 
  %【入力（任意）】
  % scale: 各軸のスケールの図上の全幅 [x_span y_spax z_span]
  %          例: x_span = XLim(2) - XLim(1)
  %          ただし、数値をそのまま使う訳ではない。
  %        　利用するのは、軸間のスケールの比率だけ。
  % axis0: マーカの大きさのスケールの基準として採用する軸の指定。
  %        　（　1: x軸、　2: y軸、　3: z軸　）
  % 
  %【出力】
  % xyz: マーカーを形成する折線の座標値（多行３列の行列）

  % [異スケール軸系]の軸間の換算係数を求める。
  if nargin>=5                  % 引数に scale と axis0 を含むとき
    scale=scale/scale(axis0);   % scaleの規格化
    Kx=scale(1); Ky=scale(2); Kz=scale(3); 
  else
    Kx=1; Ky=1; Kz=1; 
  end

  % [異スケール軸系]による法線を [等スケール軸系]に換算する。
  % 　マーカーの原型は等スケール軸系上で準備する。
  nor_p=[n0(1)*Kx  n0(2)*Ky  n0(3)*Kz];

  nor_p=nor_p/norm(nor_p);      % 貼付け面の法線ベクトルの単位化。
  nor_z=[0 0 1];                % xy面の法線ベクトル。
                                % マーカーの原型はxy面上に作る。
  % マーカーの原型を作る。
  % 　（曲面形状にフィットできるように折れ点数は多めにする）
  if mark==1
    % 〇印マーカーの座標値を生成
    a=linspace(0, 2*pi, 33);
    x=rc*sin(a);                % xy平面上に、原点を中心にした、
    y=rc*cos(a);                % 　指定半径の円形のデータを作る。
    z=0*a;
  else
    % 変形△印マーカーの座標値を生成（ x軸:右、y軸:上 ）
    x=(rc/2)*[linspace(       0,  sqrt(3), 11) ...
              linspace( sqrt(3), -sqrt(3), 11) ...
              linspace(-sqrt(3),        0, 11) ...
              linspace(       0,        0,  6) ...
              linspace(       0,      0.5,  3)];
    y=(rc/2)*[linspace(       2,       -1, 11) ...
              linspace(      -1,       -1, 11) ...
              linspace(      -1,        2, 11) ...
              linspace(       2,        0,  6) ...
              linspace(       0,        0,  3)];
    z=y*0;                      % x,y,z: すべて行ベクトル
  end

  % 貼付け後のマーカーの上下の向きを揃えるために、xy面内で事前回転。
  th=atan2(nor_p(2),nor_p(1))+pi/2;            % 回転させるべき角度。
  A=[cos(th) -sin(th);sin(th) cos(th)]*[x;y];  % 二次元回転行列で回転。
  x=A(1,:);
  y=A(2,:);                                    % zデータは変更なし。

  % 法線nor_zを法線nor_pに一致させるための回転中心軸と回転角度を求める。
  % 　スカラー積(sp)の正負で処理を分けるところが最大のキーポイント。
  sp=dot(nor_z,nor_p);          % 両ベクトルの角度差が90°以下で正。
  if sp>=0     % 90°以下
    axrot=cross(nor_z,nor_p);   % 両者の外積をとり、回転中心軸を得る。
    ang=asin(norm(axrot));      % 一致させるために回転すべき角度。
  else         % 90°超
    axrot=cross(nor_z,-nor_p);  % 回転中心軸をsp>=0のときとは逆向きに。
    ang=asin(norm(axrot))+pi;   % 一致させるために回転すべき角度。
  end

  % マーカーの面を、曲面の貼付け点の接平面と平行になるように傾ける。
  if norm(axrot)~=0    % [nor_z]≠[nor_p] のとき
    [xyz]=rot_points_3d([0 0 0],axrot,[x' y' z'],ang);
                       % このプログラムのすぐ後にある、
                       % 　多点群の三次元回転用のローカル関数の呼び出し。
                       % 　各引数の意味については関数内に詳述あり。
  else                 % [nor_z]＝[nor_p] のとき。
    xyz=[x' y' z'];    % 　回転は不要
  end

  % マーカーの座標を[等スケール軸系]から[異スケール軸系]に変換
  xyz(:,1) = xyz(:,1)*Kx;
  xyz(:,2) = xyz(:,2)*Ky;
  xyz(:,3) = xyz(:,3)*Kz;

  xyz=xyz+p0;          % 回転後のマーカーを貼付け点まで平行移動。
end

% ==========================================
% ■■■　ローカル関数 rot_points_3d　■■■
% ==========================================
% 多点群を一括して三次元回転させる関数

function [p2]=rot_points_3d(v00,v0,p1,ang)

  %【入力】
  % v00:  回転軸ベクトルの始点の座標（１行３列の行ベクトル）
  % v0:   回転軸ベクトルの向き
  %       　（１行３列の行ベクトル。単位化されていなくても可）
  % p1:   回転される前の点群の座標（多行３列の行列）
  % ang:  回転角[rad]（ v0の右ねじ方向が正。スカラー）
  % 
  %【出力】
  % p2:   回転された後の点群の座標（多行３列の行列）

  p1t=p1.';        % 点群の座標を転置して３行多列化。
  nn=v0/norm(v0);  % 回転軸ベクトルの単位化
  n1=nn(1); n2=nn(2); n3=nn(3);  % 回転軸ベクトルのx,y,z成分
  sa=sin(ang);
  ca=cos(ang);
  ca1=1-ca;        % 1-cos(ang)
  % ロドリゲスの三次元回転行列
  R=[n1^2*ca1+ca       n1*n2*ca1-n3*sa  n3*n1*ca1+n2*sa
     n1*n2*ca1+n3*sa   n2^2*ca1+ca      n2*n3*ca1-n1*sa
     n3*n1*ca1-n2*sa   n2*n3*ca1+n1*sa  n3^2*ca1+ca    ];
  p1t=p1t-v00.';   % 回転ベクトルの始点が原点に来るように、
                   % 　座標全体を平行移動する。
  p2t=R*p1t;       % v0を中心軸にして点群を回転させる。
  p2t=p2t+v00.';   % 回転前に行った平行移動を元に戻す。
  p2=p2t.';        % 転置して多行３列化したうえ、最終出力とする。
end