全天球イラスト作画補助ツール(VRビューア)をMATLABで作った話

はじめに

突然ですが、全天球イラストって描くのムズいんですよ。 えぇそれはもう本当に。
だってこんなイラスト、いくら絵心あっても描くのとかムリと思いません？上端(アタマ)のところとか 自分で描いててよくわかんねぇもん。

で、以前はその全天球イラストを描くのに必要な「パースガイド」というものをMATLABで計算するスクリプトについて紹介しました。 もちろんそういったパースガイドを計算することも大事なのですが、それと同じぐらい360度VRビューで表示確認 することが何よりも大事です。
VRブームのおかげで、全天球イラスト(正距円筒図法で描かれた画像)を360度VRビューに変換するサービスやフレームワークなどがちまたにはいっぱいありますが、それらのサービスは当然ですがVRを楽しむためのものであり、 イラスト作画用途としては使いづらい。
というわけで、全天球イラストの作画用途に特化したVRビューアをMATLABで作ってみました。

前にGUIで操作するVRビューアは作ったけど、ペンタブ握りながらGUI操作するのがメンドかったので、VRビュー操作デバイスをつくった(ハコがだせぇ)
シリアル通信(Arduino)でMATLAB Figureを操作しています。 pic.twitter.com/lt6EczcNhZ

— わたなべ (@mana_544) December 6, 2020

VRビューアシステム構成概要

まずは全体構成です。ハード構成はシンプルなんですが、ソフトの構成(MATLAB内部のプログラム)は、ちょっといろんなオブジェクトを操作しているので、ごちゃごちゃしてわかりづらいですね。わかりやすくなるかなと思ってUMLちっくに描きはじめたけど、途中でクラス図描いてるのかオブジェクト図描いてるのかよくわかんなくなりました。でも雰囲気伝わればいいから このままゴリ押しうぇーい(☝ ՞ਊ ՞)☝

今回は全天球イラスト作画用途に特化したVRビューアということで、ペンタブを持ちながらでも操作がしやすいように専用の操作デバイス(スライダーとジョイスティック)を用意しました。
デバイスの通信(シリアル通信)、正距円筒図法イメージからVRビュー(正距方位図法)への変換、結果の表示(描画)はすべてMATLABで行っています。ツールボックスは使っていません。

ここから先は、ハードデバイス〜Arduinoプログラム(シリアル送信)の説明は割愛し(シンプルなので)、MATLAB側のプログラムについて、メインとなるvrviewerオブジェクトを中心にどんな処理をしているかを説明します。(vrviewerオブジェクトはこのVRビューアのために作ったクラスです。後述します)

事前知識

話の内容に入る前に、まずは事前知識として「正距円筒図法とVRビュー」についてお話しします。

正距円筒図法(Equirectangular)

正距円筒図法については、前回記事でも書いていますのでここでは軽くおさらい程度に。
もともとは世界地図の投影法で、360度の視界を円筒状に展開して2次元画像に落とし込んだもので、冒頭のイラストが正距円筒図法で描いたものです。
全天球イラストの作画は、基本的に正距円筒図法で描きます。 描くのツライ。

VRビュー(正距方位図法: Azimuthal Equidistant Projection)

で、実際のVRビューにするためには、正距円筒図法で描かれた画像を変換する必要があります。
画像変換は有名どころでは「心射図法(ちまたのVRビュー再生はだいたいコレ)」「正射影(カーブミラーとかの絵面はコレ)」「ステレオ射影(“リトルプラネット効果”とか、インスタ映え狙うならコレ)」などいろんな手法があるのですが¹ 、今回は 正距方位図法(Azimuthal Equidistant Projection) への変換をします。

下記のイラストは、冒頭の全天球イラスト(正距 円筒 図法の画像)を、正距 方位 図法に画像変換したものになります。(中心方位角(Az)182.9[deg], 中心仰角(Ev)38.3[deg], 拡大率(R)2.6)

オーソドックスな “魚眼レンズ”の見た目 になりますね。変換式については後述します。
では以下から、実装(MATLABスクリプト)の解説になります。

MATLABでシリアル通信

まずはハードウェアとの通信まわりから。これはArduinoをシリアル通信装置として利用し、MATLABのserialportオブジェクトを利用してセンサ値をコールバックで取得しています。

以下に、serialportオブジェクトで呼び出すコールバック関数の中身を示します。

readSerialData.m

function readSerialData(src, vrObj)
% シリアル通信データ(String型)を読み取り、vrviewerオブジェクトに格納するコールバック
% 
% [arguments]
% src : serialportオブジェクト
% vrObj : vrviewerオブジェクト
% vrObj.SerialData : 3x1 double
%     シリアル通信から読み込んだ値
%     [zoom_val; x_val; y_val]
%     データ未着の場合は空[]とかNanとか
%
% [シリアル通信でくるデータフォーマット]
%     "zoom_val,x_val,y_val\r\n"
%
% zoom_val : スライダー可変抵抗値(A0ピン)
% |   0   | 1023  |
% | ----- | ----- |
% | Lower | Upper |
% 
% x_val : ジョイスティック(X)可変抵抗値(A1ピン)
% | 244 |   520   | 778  | 1023 |
% | --- | ------- | ---- | ---- |
% | Up  | Neutral | Down | Push |
% 
% y_val : ジョイスティック(Y)可変抵抗値(A2ピン)
% | 258  |   516   |  763  |
% | ---- | ------- | ----- |
% | Left | Neutral | Right |

% 3x1 String配列に格納{"zoom_val";"x_val";"y_val"}
serial_string = split(readline(src), ",");
% 3x1 数値配列に変換[zoom_val; x_val; y_val]
serial_value = str2double(serial_string);
% vrObj.SerialDataに格納
vrObj.SerialData = serial_value;

シリアル通信でセンサ値をテキスト(String型)としてもらいますので、カンマ分割して数値に変換し、vrviewerオブジェクトのインスタンス変数 vrviewer.SerialData に格納しています。
Arduinoのアナログ入力につないでいるので、[0-5V]が[0-1023]という数値(int)で受け取ります。

正距円筒図法 -> VRビュー(正距方位図法)変換

正距円筒図法からVRビュー(正距方位図法)の変換には、以前にもお世話になったアメリカ地質調査所のガチ作図マニュアルから引用した以下の変換式を利用します。²

\begin{eqnarray*}
\phi &=& \arcsin \left[ \cos c \cdot \sin \phi_1 + \left( y \cdot \sin c \cdot \cos \frac{\phi_1}{\rho} \right) \right] \\
\lambda &=& \lambda_0 + \arctan \left[ \frac{x \cdot \sin c}{\rho \cdot \cos \phi_1 \cdot \cos c - y \cdot \sin \phi_1 \cdot \sin c} \right] (但し、\phi_1 \neq \pm 90)\\
\rho &=& \sqrt{x^2 + y^2} \\
c &=& \frac{\rho}{R}
\end{eqnarray*}

全天球イラストなんて描いてる人いないし、最初どうやって調べればいいかすらわからない中でなんとかこの式見つけ出すまで1ヶ月かかったこととか本当はすっげぇ苦労話したいけど本題と逸脱するのでさらっとスルーして 上記の数式で、ある出力座標の点$(x, y)$は、どこの緯度$\phi$(仰角)と経度$\lambda$(方位角)の点から補間すればよいかがわかるようになりますので、コイツをスクリプトに実装します(下記)。

azeqdist.m

function [inR, inC] = azeqdist(centerEv, outImgSize, inW, inH, R)
% 正距円筒図法 -> 正距方位図法変換インデックス算出
% 
% 正距円筒図法から正距方位図法に変換するためのインデックスを算出します。
% この変換メソッドではcenterAzは180[deg]固定です。
% Az変換はあらかじめ行列の循環シフトを使って変換しておくこと。
% 
% [input]
% centerEv : double
%     画像中心の仰角
%     [-90, 90] [deg]で指定
% outImgSize : [1x2 double]
%     正距方位図法の出力画像サイズ[pixcel]
%     [横ピクセル, 縦ピクセル]で指定。
%     現在は「横ピクセル = 縦ピクセル」のみサポート。
% inW : double
%     入力画像サイズ(横)[pixcel]
% inH : double
%     入力画像サイズ(縦)[pixcel]
% R : double
%     拡大率
%     1 >= R で指定
% [output]
% inR : NxN double
%     正距方位図法変換インデックス(行方向)
%     NはoutImgsizeで指定した出力画像サイズ
% inC : NxN double
%     正距方位図法変換インデックス(列方向)
%     NはoutImgsizeで指定した出力画像サイズ


outW = outImgSize(1);
outH = outImgSize(2);

scale = pi;

[outC, outR] = meshgrid(1:outW, 1:outH);

x = 2 .* scale ./ outW .* outC -scale;
y = -2 .* scale ./ outH .* outR + scale;

rho = sqrt(x.^2 + y.^2);
rho(rho>R*pi) = NaN;
rho(rho==0) = 1e-10;
c = rho ./ R;

phi1 = deg2rad(centerEv);
lmd0 = deg2rad(180);

sinphi1 = sin(phi1);
cosphi1 = cos(phi1);

sinc = sin(c);
cosc = cos(c);

TMP_A = cosc .* sinphi1 + (y .* sinc .* cosphi1) ./ rho;
Ev = rad2deg(asin(TMP_A));
inR = -inH .* Ev ./180 + inH ./ 2;

TMP_Y = x .* sinc;
TMP_X = rho .* cosphi1 .* cosc - y .* sinphi1 .* sinc;
Az = rad2deg(lmd0 + atan2(TMP_Y,TMP_X));
inC = inW ./ 360 .* Az;

% inRとinCの中身(インデックス)を
% inputSize(inW, inH)内になるようにクリッピング
% だけどNaNはそのまま残しておく
inR = min(inH, max(1, inR, 'includenan'), 'includenan');
inC = min(inW, max(1, inC, 'includenan'), 'includenan');

この関数の返り値は行列を2つ返すようにしていますが、( inR, inC )

• 行列のカタチ( size(inR) と size(inC) )は、 “出力座標の点$(x, y)$”
• 行列の中身は、 “どこの緯度$\phi$(仰角)と経度$\lambda$(方位角)の点から補間するか”

表現しています。これはinterp2関数にクエリとして投げるためにこのような表現にしています(次節で説明します)。

interp2関数による画像変換

前述の数式(画像変換)を実装するにあたり、ポイントとなるのは二次元テーブルルックアップinterp2関数の使い方です。

今年のMATLABアドベントカレンダーを見ても、@kotekoさんのアフィン変換の記事とか、@shun-kusanoさんの陽炎の例など、MATLABで画像変換といえば、普通はImage Prossesing Toolboxを使います。

だがしかし、ワタシはMATLABとSimulinkしか持っていない 底辺ユーザー なので(いっぱいツールボックス使えるひとはうらやましいなあちくしょう！！！)、 MATLABの基本機能だけで画像変換を実現したい と考えました。³

そこで活躍するのがinterp2関数です。
interp2関数はいろんな引数スタイル⁴ を受け付けますが、今回は画像の変換で使いたいので、下図のようなイメージで使います。

このような画像変換やFigure再描画などをコールバック関数に書き、Timerオブジェクトで定期実行(ポーリング)しています。以下がTimerコールバックの中身になります。

timerCallbackFcn.m

function timerCallbackFcn(tObj, eObj, vrObj)
% タイマーコールバック
%
% 以下の機能をタイマー実行します。
% 1. Serial値を読み取り
% 2. centerAz, centerEv, Rの計算
% 3. 正距円筒図法 -> 正距方位図法変換(インデックス計算)
% 4. 画像変換(interp2補間)
% 5. image再描画
%
% [arguments]
% tObj    timerオブジェクト
% eObj    (未使用)
% vrObj   vrviewerオブジェクト

%% Serialから値を読み取り
serial_data = vrObj.SerialData;
% serial_dataが未着の場合はreturn
if isempty(serial_data)
    disp('data empty')
    return
elseif length(serial_data)~=3
    disp(serial_data)
    return
end

% zoom_val : スライダー可変抵抗値(A0ピン)
% |   0   | 1023  |
% | ----- | ----- |
% | Lower | Upper |
% 
% x_val : ジョイスティック(X)可変抵抗値(A1ピン)
% | 244 |   520   | 778  | 1023 |
% | --- | ------- | ---- | ---- |
% | Up  | Neutral | Down | Push |
% 
% y_val : ジョイスティック(Y)可変抵抗値(A2ピン)
% | 258  |   516   |  763  |
% | ---- | ------- | ----- |
% | Left | Neutral | Right |
zoom_val = serial_data(1);
x_val = serial_data(2);
y_val = serial_data(3);
% fprintf('zoom_val = %g, x_val = %g, y_val = %g\n', zoom_val, x_val, y_val)

% R(拡大率)算出
% 拡大率Rはスライダーの位置をそのまま変換する
R_CONVERT_X = [0 1023];
R_CONVERT_D = [1 10];
R = intp(R_CONVERT_X, R_CONVERT_D, zoom_val);
R = round(R, 2);

% Az方向変化量算出
% ジョイステック(Y)の傾き = Az変化量
AZDIFF_CONVERT_X = [258,500,521,763];
AZDIFF_CONVERT_D = [-10,  0,  0, 10];
azDiff = intp(AZDIFF_CONVERT_X, AZDIFF_CONVERT_D, y_val);

% Ev方向変化量算出
% ジョイステック(X)の傾き = Ev変化量
EVDIFF_CONVERT_X = [258,500,521,763];
EVDIFF_CONVERT_D = [ 10,  0,  0,-10];
evDiff = intp(EVDIFF_CONVERT_X, EVDIFF_CONVERT_D, x_val);

% fprintf('R = %g, azDiff = %g, evDiff = %g\n', R, azDiff, evDiff)

% 再計算&再描画の判断
% * 前回とR(拡大率)が変化
% * Az方向変化量が0以外
% * Ev方向変化量が0以外
if (vrObj.R ~= R) || ...
    (abs(azDiff) > 0.5) || ...
    (abs(evDiff) > 0.5)

    % 拡大率(R)をオブジェクトに反映
    vrObj.R = R;
    % Azはまず計算
    centerAz = vrObj.centerAz + azDiff;
    % 0 <= centerAz < 360 に収める
    if centerAz < 0
        vrObj.centerAz = 360 + centerAz;
    elseif centerAz >= 360
        vrObj.centerAz = centerAz - 360;
    else
        vrObj.centerAz = centerAz;
    end
    % Evは -90 <= Ev <= 90 に収める
    vrObj.centerEv = min(max(vrObj.centerEv + evDiff, -90), 90);
    
    % 処理負荷を減らすために、方位角方向の変換はイメージ配列を水平方向に回転シフトさせます。
    SHIFT_IMG = azmove(vrObj.IN_IMG, vrObj.centerAz, vrObj.inW);
    % 正距円筒図法 -> 正距方位図法変換 インデックス算出
    % 正距円筒図法から正距方位図法に変換するためのインデックスを算出します。
    [ind_inR, ind_inC] = azeqdist(vrObj.centerEv, vrObj.outImgSize, vrObj.inW, vrObj.inH, vrObj.R);
    % イメージ変換
    % 算出したインデックスを使って、イメージ行列をinterp2関数で補間して変換します。
    SHIFT_IMG = single(SHIFT_IMG);
    OUT_IMG = [];
    for i = 1:3
        IN = SHIFT_IMG(:,:,i);
        OUT = interp2(IN, ind_inC, ind_inR, 'linear');

        OUT_IMG = cat(3,OUT_IMG, OUT);
    end
    OUT_IMG = uint8(OUT_IMG);
    
    % 画像更新
    fprintf('AzEvR=[%.6f,%.6f,%.6f]\n', vrObj.centerAz, vrObj.centerEv, vrObj.R)
    vrObj.imObj.CData = OUT_IMG;
    
end

function vq = intp(x,v,xq)
% interp1関数のラッパー
% 外挿値がNaNになる対策として
% 外挿の場合は端の値を返す
if xq < x(1)
    vq = v(1);
elseif  xq > x(end)
    vq = v(end);
else
    vq = interp1(x, v, xq);
end

vrviewerクラス

このVRビューアのメインとなる「vrviewer」クラスです。⁵

vrviewer.m

classdef vrviewer < handle
    %VRVIEWER VRビューワー
    %   VRビューワーの情報格納用クラス
    %   
    
    properties
        inputFileName    % 入力画像(正距円筒図法)ファイル名
        IN_IMG    % 入力画像(正距円筒図法)[RxCx3]
        inH   % 入力画像サイズ(タテ)[pixcel]
        inW   % 入力画像サイズ(ヨコ)[pixcel]
        outImgSize   % 表示サイズ[outSize outSize][pixcel]
        centerAz   % 表示中心方位角[deg]
        centerEv   % 表示中心仰角[deg]
        R   % 拡大率[-]
        imObj   % イメージオブジェクト
        figObj  % Figureオブジェクト
        SerialData   % Serial通信からの値
        ardObj   % Serial通信(Arduinoデバイス)オブジェクト
        tObj   % Timerオブジェクト
    end
    
    methods
        % コンストラクタ
        function obj = vrviewer(inputFileName, outSize, usbport)
            % vrObj = vrviewer('全天球イラスト.png', 800, '/dev/cu.usbmodem****')
            % 
            % [arguments]
            % inputFileName : char
            %     正距円筒図法画像ファイル名
            % outSize : double
            %     表示サイズ(タテヨコ)[pixcel]
            % usbport : char
            %     USBポート
            
            obj.inputFileName = inputFileName;
            % 画像の読み込み
            % 正距円筒図法(Equirectangular)で描かれた画像を読み込みます。
            obj.IN_IMG = imread(inputFileName);
            [obj.inH, obj.inW, ~] = size(obj.IN_IMG);
            obj.outImgSize = [outSize outSize];
            
            % centerAz, centerEv, R の初期値
            % 中心の方位角(Az)
            obj.centerAz = 180;
            % 中心の仰角(Ev)
            obj.centerEv = 0;
            % 全天球半径・スケール(R)
            obj.R = 1.2;

            fprintf('centerAz: %g, centerEv: %g, R: %g \n', obj.centerAz, obj.centerEv, obj.R)
            
            % Figure生成(イニシャルは黒画面)
            OUT_IMG = zeros(obj.outImgSize(1), obj.outImgSize(2), 3);
            OUT_IMG = uint8(OUT_IMG);
            % イメージオブジェクト
            obj.imObj = imshow(OUT_IMG,'Border','tight');
            axObj = obj.imObj.Parent;
            % Figureオブジェクト
            obj.figObj = axObj.Parent;
            
            % Arduinoつないでるポートを指定してシリアル通信オブジェクトの生成
            obj.ardObj = serialport(usbport, 9600);
            % 終端子CR/LF
            configureTerminator(obj.ardObj,"CR/LF");
            % コールバックの設定
            % Serial通信でArduinoからセンサ値を取得し、vrObj.SerialDataに格納
            configureCallback(obj.ardObj,"terminator",@(src, ~) readSerialData(src, obj));
            
            % タイマーの設定
            % vrObj.SerialDataから[centerAz, centerEv, R]の計算
            % 正距方位図法変換
            % Figure再描画
            obj.tObj = timer('StartDelay', 1, ...
                'Period', 0.2, ...
                'ExecutionMode', 'fixedRate', ...
                'TimerFcn', {@timerCallbackFcn, obj}, ...
                'StopFcn', @timerStopFcn);
            
            start(obj.tObj)
            
            
        end
        % デストラクタ
        function delete(obj)
            disp('Close figure...')
            close(obj.figObj)
            disp('Stop&Delete Timer...')
            stop(obj.tObj)
        end
    end
end

このクラスは、vrviewerオブジェクト(インスタンス)を作ると、そのインスタンス変数に各種のオブジェクト(Serial通信オブジェクト、Figureオブジェクト、Imageオブジェクト、Timerオブジェクト)を紐付けます。また、vrviewerオブジェクトのインスタンス変数を複数のコールバックで(Serialオブジェクトの通信、Timerによるポーリング)非同期で操作・更新しています。

このクラスのポイントは、インスタンス変数をいろんな関数(コールバック)から参照するために、handleクラスを継承して 作っています。

MATLABの既定では、クラスを作ると「値クラス」という種類になります。値クラスは、関数(メソッド)の引数に渡すときにインスタンスのコピーを作るセマンティクスとなります。今回の実装ではコピーオブジェクトを操作しても意味がないので、関数(メソッド)内にオブジェクトを直接渡していじりたいのです。そのためには、 handleクラス を継承して定義します。

またhandleクラスの派生クラスではデストラクタ定義ができるので、vrviewerクラスはデストラクタに「Figureウインドウを閉じる」「Timerオブジェクトの停止と削除」を仕込み、vrviewerオブジェクトの削除と連動して後片付けするようにしました。

これで晴れて、「正距円筒図法のイメージファイル」「Figureのサイズ(pixel)」「ArduinoをつないだUBSポート」を指定したvrviewerオブジェクトのインスタンスを下記のように作れば

vrObj = vrviewer('全天球イラスト.png', 800, '/dev/cu.usbmodem****')

冒頭のTwitter動画のようなVRプレビュー実行ができます。

前にGUIで操作するVRビューアは作ったけど、ペンタブ握りながらGUI操作するのがメンドかったので、VRビュー操作デバイスをつくった(ハコがだせぇ)
シリアル通信(Arduino)でMATLAB Figureを操作しています。 pic.twitter.com/lt6EczcNhZ

— わたなべ (@mana_544) December 6, 2020

終了したい場合は

delete(vrObj)

こんな感じでインスタンスを削除すればFigureやTimerなどバックグラウンドで起動していたオブジェクトも後片付けしてくれるVRビューアの完成です。

おわりに

この記事を通して MATLABは自分のアイデアを素早くカタチにできるよ ということを言いたかったのですが、なんか実装がごちゃついててイマイチ伝わったか微妙な感じになってしまいました・・・。

ですが、たとえば今回のワタシみたいに 「絵を描きたいだけなんだけど！」 って、本来のやりたいことがある人にとっては、煩わしさを排した設計思想であるMATLABはすごくありがたかったりします。
なんたって、ワタシぁリモートリポジトリつくってる最中で git init って打ってたら 何作ろうとしてたか忘れてる ぐらいなんでね。

MATLABは自分のアイデアを素早くカタチにできるよ ということが少しでも伝われば幸いです。 (2回繰り返してごまかす)

1. 各種方位図法(VRビュー)の見た目の変化については、ワタシのFacebookページにイラスト付きで載せましたのでご参照下さい。 https://www.facebook.com/notes/661472258133775/ ↩

2. John, P. S. (1987) "Map projections: A working manual" U.S. Geological Survey Professional Paper 1395, pp. 191–202 ↩

3. そもそもImage Prossesing Toolboxには正距方位図法変換なんてないので、今回の例ではあまり使い道がないと思われます。いや持ってねぇからわからんけど。 ↩

4. interp2関数は、引数の与え方がリリースバージョンによってちょっと実装変わってたりするのでよく確認すること！ ↩

5. 別にわざわざクラス実装でなくても、スクリプトベタ打ちでも、関数スタイルでも何でもよかったのですが、他のオブジェクトと連携して扱いたいとか、App Designer(MATLABのGUI作成環境)へ将来くっつけることとかを考えたら、オブジェクト指向で書いた方がいいかな〜と思ってクラス実装にしてみました。(だけど結局メソッド作らずにコンストラクタに処理をすべて突っ込んでいるので、インスタンス化した途端に実行が走る微妙な仕様・・・) ↩