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

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

# 完成図
今回は下の図のように :heart: を一つだけ書いてみます。
![heart_fill.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/248335/aa42f9d4-f41c-b710-fc35-f8586e6f61a3.png)
なお，本文中では，一部処理のみを抜き出したコードを掲載します。
実行可能なスクリプトファイルは[こちらのgithubリポジトリ](https://github.com/larking95/heartCurveOnMatlab)にて公開しています。

# 境界面となる曲線
今回は以下の方程式[^eqn]で定義される :heart: を描きます。

```math
x^2 + \left(y - \sqrt[3]{x^2}\right)^2 = 1
```
ただし，この関数は[多価関数](https://ja.wikipedia.org/wiki/%E5%A4%9A%E4%BE%A1%E9%96%A2%E6%95%B0)なので，2つの境界面の関数に分けてしまいます。
手計算でもよいのですが，せっかくなのでMATLABのSymbolic math toolboxを用いてみます。
方程式を $y$ について解かせるためには，以下のように計算します。

```matlab
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)
```

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

```math
\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関数](https://jp.mathworks.com/help/symbolic/fplot.html)を用いて簡単にプロットしてみると，以下のようになります。

```matlab
figure;
fplot(solve(x^2 + (y - (x^2)^(1/3))^2 == 1, y), [-1, 1]);
```
![fplot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/248335/d4257484-36d7-f865-6be9-2b1c415b440e.png)

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

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

[^eqn]: いろんな数式が考案されているようですが，一番有名っぽいやつを選びました。

# 領域の塗りつぶし

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

|  | 主な関数 | 適した用途 |
|:-:|:-:|:-:|
| 1. [グリッドデータを用いる方法](https://jp.mathworks.com/help/matlab/surface-and-mesh-plots-1.html) | [contourf](https://jp.mathworks.com/help/matlab/ref/contourf.html) | [機械学習モデルの判定の可視化](https://jp.mathworks.com/help/stats/support-vector-machines-for-binary-classification.html) |
| 2. [多角形を表示する方法](https://jp.mathworks.com/help/matlab/geometry.html) | [fill](https://jp.mathworks.com/help/matlab/ref/fill.html) | [境界線の内部領域の可視化](https://jp.mathworks.com/matlabcentral/answers/180829-shade-area-between-graphs) |

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

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

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

## 1. グリッドデータを用いる方法
今回は簡易的に予測器の関数を作成して，グリッドデータを作成してみます。
適当に構造体とデータを受け取って，logicalベクトルを返す関数を作ります[^predict]。

```matlab
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$ ごとに細かく区切ります。
そして，区切ってサンプリングした全ての点において判別を行わせます。

```matlab 
% 分類学習器のようなモデルを作成
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関数](https://jp.mathworks.com/help/matlab/ref/colormap.html)と[pbaspect関数](https://jp.mathworks.com/help/matlab/ref/pbaspect.html)を用いて，色と縦横比を整えておきます。

```matlab
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])へ
```

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

![heart_grid.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/248335/e7cfb77f-c73d-3fb3-38d1-a0a83c201d16.png)


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

[^predict]: 実際の用途では，学習済みのモデルを受け取ってくれる[predict](https://jp.mathworks.com/help/stats/compactclassificationensemble.predict.html)関数などを利用します。

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

多角形を表示するためには，多角形の頂点の座標が必要となります。
そこで，頂点のx座標とy座標をそれぞれ作成します。
ここで，[fliplr関数](https://jp.mathworks.com/help/matlab/ref/fliplr.html)[^fliplr]を用いることに注意です。
これは，fill関数が最初の頂点と最後の頂点を自動的につないで，多角形を閉じるため必要となります。

```matlab
vX = [x, fliplr(x)];    % 頂点のx座標
vY = [yu, fliplr(yl)];  % 頂点のy座標
```
頂点座標を作成できれば，後は塗りつぶし色を指定してfill関数を実行するだけです。
また，グリッドを用いた方法と同様に，縦横比と表示範囲を整えておきます。

```matlab
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方向の表示範囲
```

これにより，重くない愛 :heart: を作成することができました。
![heart_fill.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/248335/aa42f9d4-f41c-b710-fc35-f8586e6f61a3.png)


[^fliplr]: 配列の左右を反転する関数（flip Left to Right）です。xやyが縦ベクトルなら，[flipup関数](https://jp.mathworks.com/help/matlab/ref/flipud.html)(flip Up to Down)を利用すれば良いです。

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

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

- SVM の識別範囲の可視化
- ロバスト性能を保証可能な領域を表示[^Rbode]

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

# 参考文献
1. [美しすぎるハートの数式](https://itchyny.hatenablog.com/entry/20130214/1360847516)（itchynyさんのブログ）
2. [Heart Curve](http://www.mathematische-basteleien.de/heart.htm) (英語サイト)
2. [グラフ間の影の領域(shade area between graphs)](https://jp.mathworks.com/matlabcentral/answers/180829-shade-area-between-graphs)（MATLAB Answers:英語）


#おまけ
## 補足1：多角形を表示する方法でfliplrを忘れたらどうなる？
実際にやってみると，こちらのようになりました。
始点と終点の対応関係が逆になるため，塗られる領域も逆になってしまっていますね。

![heart_fill.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/248335/2b65841a-bd4d-2f0d-2900-1842c858014e.png)


## 補足2：area関数について
塗りつぶしができそうな関数として，今回紹介した以外に[area関数](https://jp.mathworks.com/help/matlab/ref/area.html)があります。
area関数では，任意のBaseValueから1つの曲線までの領域を塗りつぶすことができます。
しかし，これを2つの曲線の間の領域を塗りつぶすために使うのは得策ではありません。

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

```matlab
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$となる場合には正しくプロットができません。

[^Rbode]: 千葉工業大学の熱海先生が研究されている[RBode plot](https://sites.google.com/p.chibakoudai.jp/mechatronics-lab/%E7%A0%94%E7%A9%B6%E5%86%85%E5%AE%B9/%E5%8F%AF%E8%A6%96%E5%8C%96%E3%81%95%E3%82%8C%E3%81%9F%E5%88%B6%E5%BE%A1%E7%B3%BB%E8%A8%AD%E8%A8%88?authuser=0)のトレース