はじめに

今回はMATLABを使ってを描く方法を紹介したいと思います。
を描くためには，2つの曲線の間の領域を塗りつぶす必要があります。
簡単そうに思えるのですが，いくつかやり方があるのでそれぞれ紹介します。

完成図

今回は下の図のようにを一つだけ書いてみます。

なお，本文中では，一部処理のみを抜き出したコードを掲載します。
実行可能なスクリプトファイルはこちらのgithubリポジトリにて公開しています。

境界面となる曲線

今回は以下の方程式¹で定義されるを描きます。

x^2 + \left(y - \sqrt[3]{x^2}\right)^2 = 1

ただし，この関数は多価関数なので，2つの境界面の関数に分けてしまいます。
手計算でもよいのですが，せっかくなのでMATLABのSymbolic math toolboxを用いてみます。
方程式を $y$ について解かせるためには，以下のように計算します。

syms x y;
solve(x^2 + (y - (x^2)^(1/3))^2 == 1, y)
% ans =
%  
%  (1 - x)^(1/2)*(x + 1)^(1/2) + (x^2)^(1/3)
%  (x^2)^(1/3) - (1 - x)^(1/2)*(x + 1)^(1/2)

つまり，以下の二つの曲線をプロットすれば，をプロットすることができます。

\begin{align}
y_u &= \sqrt[3]{x^2} + \sqrt{1 - x} \sqrt{x + 1}\\
y_l &= \sqrt[3]{x^2} - \sqrt{1 - x} \sqrt{x + 1}
\end{align}

試しに，Symbolic math toolboxのfplot関数を用いて簡単にプロットしてみると，以下のようになります。

figure;
fplot(solve(x^2 + (y - (x^2)^(1/3))^2 == 1, y), [-1, 1]);

きれいなハートマークになっています。

しかし，これではあまりにシンプルすぎて愛を感じません。
目標とするには程遠いです。
そこで，この2つの曲線の間の領域を塗りつぶしていこうと思います。

領域の塗りつぶし

領域を塗りつぶす方法には，大きく分けて2つの方法があります。

	主な関数	適した用途
1. グリッドデータを用いる方法	contourf	機械学習モデルの判定の可視化
2. 多角形を表示する方法	fill	境界線の内部領域の可視化

1つ目の方法は，グリッドデータを作って等高線を描く方法です。
この方法は，境界線が厳密にわかっていない場合に適しています。
しかし，グリッドデータが膨大になると計算時間がかかります。

2つ目の方法は，内部を塗りつぶした多角形を表示する方法です。
この方法は境界線が既知である必要があります。
しかし，その分計算時間は1つ目の方法よりも短くなります。

今回のの場合は境界線が既知なので，2つ目の方法が適しています。
ですが，紹介のために両方の方法で描いてみます。

1. グリッドデータを用いる方法

今回は簡易的に予測器の関数を作成して，グリッドデータを作成してみます。
適当に構造体とデータを受け取って，logicalベクトルを返す関数を作ります²。

function label = isBetween2Curves(model, data)
    validateattributes(model, {'struct'}, {'scalar'});
    validateattributes(model.yu, {'function_handle'}, {'scalar'});
    validateattributes(model.yl, {'function_handle'}, {'scalar'});
    validateattributes(data, {'numeric'}, {'2d'});

    % 変数の置きなおし
    xi = data(:, 1);
    yi = data(:, 2);

    % falseですべて埋めたベクトルを用意
    label = false(size(xi));

    % データのxに対する上界と下界を計算
    yu = model.yu(xi);
    yl = model.yl(xi);

    % 内部領域にある場合のみ，trueを代入
    label(yl < yi & yi < yu) = true;
end

次に，meshgrid関数を用いて，平面を間隔 $dx$ ごとに細かく区切ります。
そして，区切ってサンプリングした全ての点において判別を行わせます。

% 分類学習器のようなモデルを作成
model = struct(...
    'yu', @(x) (1 - x).^(1/2).*(x + 1).^(1/2) + (x.^2).^(1/3),...
    'yl', @(x) (x.^2).^(1/3) - (1 - x).^(1/2).*(x + 1).^(1/2));

% グリッドデータを作成
[X, Y] = meshgrid(-1.2:dx:1.2, -1.2:dx:2);

% 全ての組み合わせについて判別させる
Z = isBetween2Curves(model, [X(:), Y(:)]);

% 判別結果をXやYと同じサイズの行列に整形する
Z = reshape(Z, size(X));

これにより，平面上の点X,Yに対応する判定結果をZに格納することができました。
最後に，これを用いて，contourf関数を用いてプロットします。
また，colormap関数とpbaspect関数を用いて，色と縦横比を整えておきます。

fig_grid = figure('Name', 'グリッドデータを用いる方法');
contourf(X, Y, Z);          % 塗りつぶし等高線をプロット
pbaspect([1, 1, 1]);        % それぞれの軸の比を1:1:1に固定
colormap([1,1,1; 1,0,0]);   % false(0)を白([1,1,1])，true(1)を赤([1,0,0])へ

これにより，燃えるように真っ赤なをプロットすることができました。

ただし，この方法は重い処理となります。
サンプル点を増やした場合にはメモリ不足になる可能性もあります。
ですので，可能であれば次に紹介する方法を使います。

2. 多角形を表示する方法

ここでは，多角形を表示することで先ほどと同じような見た目を実現します。

多角形を表示するためには，多角形の頂点の座標が必要となります。
そこで，頂点のx座標とy座標をそれぞれ作成します。
ここで，fliplr関数³を用いることに注意です。
これは，fill関数が最初の頂点と最後の頂点を自動的につないで，多角形を閉じるため必要となります。

vX = [x, fliplr(x)];    % 頂点のx座標
vY = [yu, fliplr(yl)];  % 頂点のy座標

頂点座標を作成できれば，後は塗りつぶし色を指定してfill関数を実行するだけです。
また，グリッドを用いた方法と同様に，縦横比と表示範囲を整えておきます。

fig_fill = figure('Name', '多角形を表示する方法');
fill(vX, vY, 'r');      % 多角形の表示
pbaspect([1, 1, 1]);    % それぞれの軸の比を1:1:1に固定
xlim([-1.2, 1.2]);      % x方向の表示範囲
ylim([-1.2, 2.0]);      % y方向の表示範囲

これにより，重くない愛を作成することができました。

おわりに

読んでくださった皆様に私のMATLAB は伝わったでしょうか？

今回は冗談っぽくハートの方程式を塗りつぶしましたが，実際の研究でも以下のような場面で塗りつぶしをしたいときがありました。

SVM の識別範囲の可視化
ロバスト性能を保証可能な領域を表示⁴

他にもきっと役立つ場面があると思うので，その際に皆様の参考になれば幸いです。

参考文献

美しすぎるハートの数式（itchynyさんのブログ）
Heart Curve (英語サイト)
グラフ間の影の領域(shade area between graphs)（MATLAB Answers:英語）

おまけ

補足1：多角形を表示する方法でfliplrを忘れたらどうなる？

実際にやってみると，こちらのようになりました。
始点と終点の対応関係が逆になるため，塗られる領域も逆になってしまっていますね。

補足2：area関数について

塗りつぶしができそうな関数として，今回紹介した以外にarea関数があります。
area関数では，任意のBaseValueから1つの曲線までの領域を塗りつぶすことができます。
しかし，これを2つの曲線の間の領域を塗りつぶすために使うのは得策ではありません。

確かに今回紹介したハートの方程式程度であれば，area関数でもプロットすることができます。

fig_area = figure('Name', 'Area');
area(x, yu, min(yu), 'FaceColor', 'r', 'ShowBaseLine', 'off', 'LineStyle', 'none');
hold on;
area(x, yl, max(yl), 'FaceColor', 'r', 'ShowBaseLine', 'off', 'LineStyle', 'none');
hold on;
plot(x, yu, 'k', x, yl, 'k');
pbaspect([1, 1, 1]);    % それぞれの軸の比を1:1:1に固定
xlim([-1.2, 1.2]);      % x方向の表示範囲
ylim([-1.2, 2.0]);      % y方向の表示範囲

しかし，area関数にはBaseValueが柔軟に変更できない問題があります。
どうやら，BaseValueはAreaクラスの静的メンバとなっているようです。
そのため，$\min y_u < \max y_l$となる場合には正しくプロットができません。

いろんな数式が考案されているようですが，一番有名っぽいやつを選びました。 ↩
実際の用途では，学習済みのモデルを受け取ってくれるpredict関数などを利用します。 ↩
配列の左右を反転する関数（flip Left to Right）です。xやyが縦ベクトルなら，flipup関数(flip Up to Down)を利用すれば良いです。 ↩
千葉工業大学の熱海先生が研究されているRBode plotのトレース ↩

We love MATLAB（領域の塗りつぶし）