MATLABでPhotoShopの深度ぼかしを再現する。

# はじめに
だいぶ前の話なんですが、大学の授業でMATLABを触ったので、話題のAIと組み合わせてPhotoShopの深度ぼかしフィルターを再現することにしました！今回は画像とその深度マップ、焦点距離を入力とし、その距離に応じてフィルタの強弱を変化させることで、実際の一眼レフカメラで撮影したような画像を生成するフィルタの実装を目指します。

# 目次
1. [そもそも深度ぼかしって？](#Chapter1)
1. [ガウシアンぼかしとレンズのボケの違い](#Chapter2)
1. [ボケの実装](#Chapter3)
1. [焦点距離の実装](#Chapter4)
1. [実装結果](#Chapter5)
1. [参考文献](#reference)

<!-- 各チャプター -->
<a id="#Chapter1"></a>
# そもそも深度ぼかしって？
PhotoShopの「深度ぼかし」機能では、AIによって入力画像の深度を推定し、まるで一眼レフで撮影したような写真に加工することができます。恐らくぼかしの強度を、画素の深度によって変えているんだろうとおもいます。

![deep-blur-s.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3654050/62ebd957-58fd-5e58-9094-c22f09c267db.png)
[Photoshopで立体的なぼかし加工を行う方法](https://321web.link/depth-blur/)

<a id="#Chapter2"></a>
# ガウシアンぼかしとレンズのボケの違い
さて、まずボケの実装をしましょう。MATLABのimgaussfilt関数でボカせばいいと思うかもしれませんがそんなに単純な問題でもありません。下の画像から分かる通り、レンズによってできるボケは明るい部分が暗い部分を侵食するようにボケていることが分かります。

![レンズブラー.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3654050/2ae1e7eb-4f3e-f80e-02bc-c23014fa78a3.png)

つまり、輝度の高い部分から優先的に平均化フィルタやガウシアンフィルタといった畳み込み処理を行う必要があります。これはレンズが対数的な特性を持っていることに由来しています。
レンズ特性については天下のCanon様の記事で詳しいことが説明されているので気になった方は見てみてください。
[【CANON LOG】TECHNICAL SITE](https://canon.jp/business/solution/pro-img-sys/cinema-eos/canon-log/history)

<a id="#Chapter3"></a>
# ボケの実装
上で述べたようなレンズ特性を再現するには、
1. 入力画像に対して指数変換を行う
1. 指数変換した画像に対して畳み込み処理を行う
1. 畳み込み処理後の画像に対して対数変換する

という手順を踏む必要があります。

### 1. 指数変換
X を入力配列とすると指数変換後の配列X ’は、
```math
X' = X^{k}
```
となります。kはフィルムの特性など、撮影画像情報がもつ所定の非線形性で決めなければなりません。k=1.01～1.10 の範囲で調整することで実際のレンズで撮影した画像と似たような結果が得られました。
```MATLAB
I=imread("matlab.png");
Id=double(I);

k = 1.05;
R = k.^Id(:,:,1);
G = k.^Id(:,:,2);
B = k.^Id(:,:,3);
```
### 2. 畳み込み処理
1の画像に対して畳み込み処理をします。今回はより畳み込み処理後に画素付近の情報を残すことのできるガウシアンフィルタを採用しました。Matlabの標準関数を使ってもいいのですが勉強のために自作しました。適当に記述すると黒い枠のようなものが出来上がってしまうため、今回は画像サイズからはみ出さない範囲で処理を行うよう記述します。
```MATLAB
function gauss=gaussian(I,sigma)
sz=3;
[x,y]=meshgrid(-sz:sz,-sz:sz);
X=size(x,1)-1;
Y=size(y,1)-1;

exponential =-(x.^2+y.^2)/(2*sigma*sigma);
Kernel=exp(exponential)/(2*pi*sigma*sigma);

gauss=zeros(size(I));
I=padarray(I,[sz sz]);

%画像サイズからはみ出さない範囲で処理
for i=1:size(I,1)-X
    for j=1:size(I,2)-Y
        Temp=I(i:i+X,j:j+Y).*Kernel;
        gauss(i,j)=sum(Temp(:));
    end
end
end
```
### 3. 対数変換
次に2の画像に対して対数変換をします。指数変換で用いた式に対して逆関数を求め、底の変換公式より
```math
X' = \frac{log{X}}{log{k}}
```
となります。
```MATLAB
R_max = max(gR,1);
R_log = (log(R_max) ./ log(k) /255);
G_max = max(gG,1);
G_log = (log(G_max) ./ log(k) /255);
B_max = max(gB,1);
B_log = (log(B_max) ./ log(k) /255);

%ブラーをかけた画像
bI = cat(3, R_log, G_log, B_log);
```
### 出力結果
ここで、定番のlennaの写真を用いて出力結果を確認してみましょう。

左：入力画像　右：出力画像（σ=3）
![lenna比較.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3654050/53cb9566-33c4-a65c-ec4c-7486e99c5bfe.png)

綺麗にボケていますね、これだけでは本当に輝度の高い部分から畳み込まれているかわからないのでMATLABのimgaussfilt関数と比べてみました。

左：ガウシアンブラー（σ=2）右：提案手法での出力画像（σ=2）
![gauss比較.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3654050/63db2b2b-5d2a-bfcd-d1e6-dcc5854b0d5f.png)

比べると分かりやすいですが、確かに明るい部分が暗い部分を侵食してボケています。

<a id="#Chapter4"></a>
# 焦点距離の実装
さて、焦点距離の実装はどうやるのかというと、まず入力画像の深度マップを用意します。今回は単眼深度推定モデルとして有名なMiDaSを使用しました。以下のサイトを参考に、入力画像の深度マップを出力します。
 [【MiDaS】深度推定モデルをいろんなモノで試してみる](https://farml1.com/midas/)

そして出力した深度マップを深度ごと分解して2値化し、多層マスクを作成します。このマスクに対して異なる強度のボケ画像を適用し、最後に合成するわけです。分かりずらいかもしれないので図を用意しました。

![手法図.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3654050/2001564a-7310-9a68-5059-21dca5368cb4.jpeg)

図のように、例えば深度が奥へ行くほど強い強度のフィルタを適用するようにすれば、奥がボケて手前に焦点が合ったような画像が作れるわけです。（なんで手間が奥になってるんだよ、分かりずらい図だな。）実際のプログラムではユーザーが指定した焦点距離から離れれば離れるほどボケるようなアルゴリズムにしました。

```MATLAB
%多層マスク作成
a=1;
for i = [0:51:204]
depth = im2gray(depth);
mask_list(:,:,a) = i<depth & depth<=i+51;
a=a+1;
end
```
```MATLAB
%指定した焦点距離より手前の処理
for i = [1:r]
a=abs(r-i)+1;
bI=bI_list(:,:,:,a);
layer(:,:,1)=bI(:,:,1).*mask_list(:,:,i);
layer(:,:,2)=bI(:,:,2).*mask_list(:,:,i);
layer(:,:,3)=bI(:,:,3).*mask_list(:,:,i);
layer_list(:,:,:,i)=layer;
end

%指定した焦点距離より奥側の処理
for i = [r:5]
a=abs(r-i)+1;
bI=bI_list(:,:,:,a);
layer(:,:,1)=bI(:,:,1).*mask_list(:,:,i);
layer(:,:,2)=bI(:,:,2).*mask_list(:,:,i);
layer(:,:,3)=bI(:,:,3).*mask_list(:,:,i);
layer_list(:,:,:,i)=layer;
end
```

<a id="#Chapter5"></a>
# 実装結果
rが焦点距離で、ユーザー側から指定できるようになっています。図のように、rを変えることで焦点の合っている物体が、円錐→球→立方体と変わっているのが分かるかと思います。

![焦点距離比較.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3654050/865c92bd-7753-0ee3-1f6e-3c22d5f58b66.png)

今回は簡易的な実装なので5層にしましたが、さらに多層化することでより自然な焦点ぼかしが実装できるでしょう。ひとまず、MATLAB初心者の実装としては満足な出来でした。

<a id="#reference"></a>
# 参考文献
- [【CANON LOG】TECHNICAL SITE](https://canon.jp/business/solution/pro-img-sys/cinema-eos/canon-log/history)
- [デジカメ画像をキレイにボカそうアルゴリズム編](http://www.hirax.net/dekirukana5/bokeboke1/)
- [【MiDaS】深度推定モデルをいろんなモノで試してみる](https://farml1.com/midas/)
- [Photoshopで立体的なぼかし加工を行う方法](https://321web.link/depth-blur/)