Devil With Angel Heart～ MATLABで綺麗な「ボケ」を作ろう～決め手は対数～巡回畳み込み・直線畳み込み、MATLABでブロック処理

・絞り

カメラのレンズには「絞り」という露光量・被写界深度を調整する機構があります。

写真の「ボケ」はこの絞りの形になります。

これは上の7角形絞りレンズで撮った写真です。

こちらは丸絞りレンズで撮った写真。

背景のボケが、絞りの形になっているのがお分かりになるかと思います。

絞りの代わりに、レンズの前に何かの形に切り抜いた黒い紙等を付ければ、その形のボケを出すこともできます。

今回はすでに撮った写真のボケを、信号処理でハート型のボケに変換してみましょう。

被写体の一点から発せられた光は拡散し、レンズと絞りを通ってレンズによって（ピントの合った部分は）再び一点に集約されます。このような現象は、信号処理的には被写体と絞りの「畳み込み演算」で再現できます。

但し、後述しますが、合成マスク（深度マップ相当）の作成にPhotoshopを用いています。

必要なツールは以下の通りです。
MATLAB
Image Processing Toolbox：blockproc()　（FFT版のみ使用）
（Photoshop：マスク作成用）

・対数変換

こういった処理を行う場合、注意点が二つあります。

一つ目は、RGBの振幅にそのまま処理をしても実際と同じ画像は得られないということです。

フィルムの露光量－濃度グラフは、露光量がlogで示されます。
片側対数グラフでほぼ直線となるからです。これは、フィルムに写る「光の強さ」が（人の感覚もそうですが）対数的ということです。

　参考：ネガフィルムの特性曲線（「Logの歩み」タブ）

そのため、実際と同様の効果を得るには、RGB信号の振幅そのままではなく対数変換してから処理を行う必要があります。

0～255レベルのRGB信号をxとすると
　x=0:255;
　y=k.^x;
とします。kはフィルム等の特性によって決まる定数です。

私が以前デジ一で実際のレンズ画像と比較したとき、k=1.035^(R/G/B)でほぼ似たような特性が再現できたので、ここではその数値を使います。

戻すときは
　log_1.035(y)=log(y)/log(1.035);
を取ります。
（log(0)は定義されませんので1でリミットをかけます）

この変換特性は、片側 logスケールで直線となります。

直線の傾きは、kによって変わります。

デジタルカメラでももちろん、このようなlog特性を使ってます。

参考：Canon Log 3について

・演算量削減

もう一つは、処理量削減のためFFTを用いることを検討するということです。
時間軸上でやろうとすると、画像と絞り形状（以下、カーネルと呼びます）を1ピクセルずつズラしながら、処理の重い「畳み込み演算」を行う必要があります。（折り返し雑音とは？～なぜ車のホイールは逆回転するのか～デジタルはアナログの近似ではない　->　数学的補足　->　畳み込み（convolution）参照）

畳み込み演算はFFT（DFT）を用いて周波数軸上に変換すれば、単純な掛け算で同じ処理が行えます。ただし、FFT処理が追加になるため、画像・カーネルサイズや処理フレームワークによってどちらが速いかは異なりますので、その都度確認する必要はあります。

・巡回畳み込みと直線畳み込み

「巡回畳み込み」は、DFTの積のIDFTと等価です。DFTは、フレーム長の信号が無限に繰り返し続いているとして処理を行います。

これは離散化に起因します。これも上に上げた折り返し雑音とは？に少し書いてあるのですが、下の図で示したように時間軸と周波数軸の双対性により、時間軸でアナログ(連続）信号をあるフレーム長の窓で切り取り、サンプリングによって離散化すると周波数軸上で同じ周波数成分が無限に繰り返し、それをさらに周波数軸上でサンプリングすると時間軸上でも窓関数により切り取ったフレーム長で時間信号が無限に繰り返すことになります。これはDFTの導出に他なりません。言い方を変えると、デジタル信号処理においては時間軸・周波数軸双方でデータが離散化されている必要があり、その結果、時間軸・周波数軸とも同じデータが無限に繰り返すことになる、とも言えます。

しかしこれはあくまでもデジタルの世界の話。現実世界はこういう繰り返しにはなっていないので、「直線畳み込み」を求める必要があります。

デジタル信号処理で直線畳み込みを行うためには、以下のようなプロセスを踏みます。
　１．周りを0で埋めて2Nｘ2Nサイズに拡張
　２．FFT
　３．カーネルのFFTと複素掛け算
　４．IFFT＋周波数シフト
　５．振幅に戻す
　６．元のNｘNサイズを切り出す

N×N点の直線畳み込み結果は2N-1点になるのでそれ以上の長さを取っておけば信号が重ならず、巡回畳み込みとなるDFTを使って直線畳み込みと同一のデータを取り出すことができます。FFTしやすいように、2^n点とすることが多いです。（小さな素数の積でも可）

・ブロック処理

MATLABでは、blockproc() というブロック処理用の関数が用意されています。周辺を拡大して0埋め～結果からは拡大分を削除、が簡単にできます。

具体的には以下のようにします。

block_size = [32 32];
border_size = [16 16];

blk_func = @(block_struct) abs(fftshift(ifft2(fft2(block_struct.data) .* fkernel64)));
bR = blockproc(R_log,block_size,blk_func,'BorderSize',border_size, 'PadPartialBlocks',true);

32x32ブロックに分割～その周辺を0埋めして64x64に拡大～ブロック単位処理～処理結果の中心32x32を再切りだし、を blockproc() の一行で行えます。

そのブロック単位処理の中身（上記１～５）は、block_struct を使ってこれも一行で定義してあります。

・マスクの作成

実際のカメラであれば、フォーカス位置と被写界深度によってボケが異なるため自動的に処理がなされるのですが、後から行う場合はその「深度マップ」に相当するマスクが必要となります。本来はボケの大きさによってカーネルサイズも変えた方が良いのですが、今回は一定サイズとしています。
（画面全体を均一にボカして良いのであればマスク自体必要ありません）

iPhoneのポートレートモードだと深度マップも持っています。
Mac版のPhotoshopでは深度マップも読めるので、それを使っても良いですね。

今回は、マスクはPhotoshopの被写体検出＋手作業で作りました。

・ハートを描く

画像として用意しても良いのですが、せっかくのMATLABなので関数で描いてみましょう。

ハート型を描く関数は何種類も提案されているのですが、今回は、【ハートのグラフ選手権】で優秀賞を取られた、CHARTMANさんのを使わせていただきました。

どうやったらこういうのを考えつくんですかね～？(ﾟ～ﾟ)

t=0で大きさ最大としたかったため、分子はcosに変えています。

・MATLABコード

サンプル画像付きMATLABコード

function()にしてあるので以下のようにして呼び出します。

heartBlurer_conv(sourceFile,maskFile,outputFile,auto,we)

heartBlurer_FFT(sourceFile,maskFile,outputFile,auto,we)

sourceFile　　元画像ファイル名
maskFile 　　 マスク画像ファイル名
outputFile　　出力画像ファイル名
auto　　　　  true：ハートの大きさをループで変化させて表示
　　　　　　  false：１枚の静止画のみを表示
we　　　　　 true：静止画とgifファイルを保存
　　　　　　　（ファイル名はコード中で固定）
　　　　　　  false：静止画のみ保存

arguments 宣言しておくとそれだけでデフォルトパラメータとして使えますので、以下のように引数を省略した呼び出しも可能です。

heartBlurer_conv()

heartBlurer_conv("yumemi02a_019.jpg","yumemi02a_019_mask7.jpg")

また、自作関数でもエディター上で引数が表示されるようになります。
So Cool！

畳み込み版ソースコード

function heartBlurer_conv(sourceFile,maskFile,outputFile,auto,we)
arguments
 sourceFile string = "yumemi02a_019.jpg";
 maskFile string = "yumemi02a_019_mask7.jpg";
 outputFile string = "heartBluer.png";
 auto logical = false; % auto play
 we logical = false; % gif file write enable
end

close all

if ~auto
 we = false;
end

step = pi/8 + pi/8*(1-we);
ed = pi + 3*pi*(1-we);

A = double(imread(sourceFile)); % Source image
mask = imread(maskFile); % mask image
mask = double(mask(:,:,1)) / 255;

gain = 2; % kernel gain

filename = "heartsBokeh.gif"; % output gif file name

f1 = figure(1);
f1.Color = [0 0 0]; % background color is black
ax1 = gca;

th = 0:pi/50:2*pi;

t2 = 1.57; % size of a heart

if auto
 t2 = 0:step:ed;
end

for t = t2
 r = (4 + cos(t))./sqrt(1-abs(cos(th)).*sin(th)); % Heart function
 x = r.*cos(th);
 y = r.*sin(th);

 fill(ax1,x,y,'w') % draw heart
 axis(ax1,'equal')
 axis(ax1,'manual')
 xlim(ax1,[-12 12])
 ylim(ax1,[-12 12])
 axis(ax1,'off');

 kernelF = getframe(ax1); % get heart as a image

 kernel64 = double(imresize(kernelF.cdata,[64 64])) / 255;
 s = sum(kernel64(:,:,1),'all');
 kernel64 = kernel64 / s * gain; % normalize and set gain
 kernel32 = kernel64(16:47,16:47,1);

 % transfer to log
 k = 1.035;
 R_log = k.^A(:,:,1);
 G_log = k.^A(:,:,2);
 B_log = k.^A(:,:,3);

 % take 2D convolutions
 bR = conv2(R_log,kernel32);
 bG = conv2(G_log,kernel32);
 bB = conv2(B_log,kernel32);

 % to linear
 R_out = max(bR,1);
 R_lin = (log(R_out) ./ log(k)) / 255;
 G_out = max(bG,1);
 G_lin = (log(G_out) ./ log(k)) / 255;
 B_out = max(bB,1);
 B_lin = (log(B_out) ./ log(k)) / 255;

 bA = cat(3, R_lin, G_lin, B_lin); % to color image
 offset = 16 + 2;
 bA = bA(offset:size(A,1)+offset-1, offset:size(A,2)+offset-1,:); % cut out border redundancy

 blendI = bA; % set size
 % blend source image and bokeh image using mask
 blendI(:,:,1) = bA(:,:,1) .* mask + A(:,:,1)/255 .* (1-mask);
 blendI(:,:,2) = bA(:,:,2) .* mask + A(:,:,2)/255 .* (1-mask);
 blendI(:,:,3) = bA(:,:,3) .* mask + A(:,:,3)/255 .* (1-mask);

 f3 = figure(3);
 f3.Visible = 'on'; % view outside of Live Script window

 if exist('ax2','var')
 image(ax2,'CData',blendI); % change data only
 else % first time
 ax2 = gca;
 imshow(blendI,'Parent',ax2)
 end
 drawnow

 if we
 frame = getframe(f3); % capture with border frame
 % frame = getframe(gca); % take image only
 gr = frame2im(frame);
 [bA,map] = rgb2ind(gr,256);
 if (t==0)
 imwrite(bA,map,filename,'gif','LoopCount',Inf,'DelayTime',1/15);
 else
 imwrite(bA,map,filename,'gif','WriteMode','append','DelayTime',1/15);
 end
 end

end

imwrite(blendI,outputFile)
end

FFT版ソースコード

function heartBlurer_FFT(sourceFile,maskFile,outputFile,auto,we)
arguments
 sourceFile string = "yumemi02a_019.jpg";
 maskFile string = "yumemi02a_019_mask7.jpg";
 outputFile string = "heartBluer_FFT.png";
 auto logical = false; % auto play
 we logical = false; % gif file write enable
end

close all

if ~auto
 we = false;
end

step = pi/8 + pi/8*(1-we);
ed = pi + 3*pi*(1-we);

A = double(imread(sourceFile)); % Source image
mask = imread(maskFile); % mask image
mask = double(mask(:,:,1)) / 255;

gain = 2; % kernel gain

filename = "heartsBokeh.gif"; % output file name

f1 = figure(1);
f1.Color = [0 0 0]; % background color is black
ax1 = gca;

th = 0:pi/50:2*pi;

t2 = 1.57; % size of a heart

if auto
 t2 = 0:step:ed;
end

for t = t2
 r = (4 + cos(t))./sqrt(1-abs(cos(th)).*sin(th)); % Heart function
 x = r.*cos(th);
 y = r.*sin(th);

 fill(ax1,x,y,'w') % draw heart
 axis(ax1,'equal')
 axis(ax1,'manual')
 xlim(ax1,[-12 12])
 ylim(ax1,[-12 12])
 axis(ax1,'off');

 kernelF = getframe(ax1); % get heart as a image

 kernel64 = double(imresize(kernelF.cdata,[64 64])) / 255;
 s = sum(kernel64(:,:,1),'all');
 kernel64 = kernel64 / s * gain; % normalize and set gain

 % frequency domain kernel
 fkernel64 = fft2(kernel64(:,:,1));

 % transfer to log
 k = 1.035;
 R_log = k.^A(:,:,1);
 G_log = k.^A(:,:,2);
 B_log = k.^A(:,:,3);

 % block processing parameters
 block_size = [32 32];
 border_size = [16 16];

 % define core function
 blk_func = @(block_struct) abs(fftshift(ifft2(fft2(block_struct.data) .* fkernel64)));

 % block processing
 bR = blockproc(R_log,block_size,blk_func,'BorderSize',border_size, 'PadPartialBlocks',true);
 bG = blockproc(G_log,block_size,blk_func,'BorderSize',border_size, 'PadPartialBlocks',true);
 bB = blockproc(B_log,block_size,blk_func,'BorderSize',border_size, 'PadPartialBlocks',true);

 % to linear
 R_out = max(bR,1);
 R_lin = (log(R_out) ./ log(k)) / 255;
 G_out = max(bG,1);
 G_lin = (log(G_out) ./ log(k)) / 255;
 B_out = max(bB,1);
 B_lin = (log(B_out) ./ log(k)) / 255;

 bA = cat(3, R_lin, G_lin, B_lin); % to color image
 bA = bA(1:size(A,1), 1:size(A,2),:); % cut out border redundancy

 blendI = bA; % set size
 % blend source image and bokeh image using mask
 blendI(:,:,1) = bA(:,:,1) .* mask + A(:,:,1)/255 .* (1-mask);
 blendI(:,:,2) = bA(:,:,2) .* mask + A(:,:,2)/255 .* (1-mask);
 blendI(:,:,3) = bA(:,:,3) .* mask + A(:,:,3)/255 .* (1-mask);

 f3 = figure(3);
 f3.Visible = 'on'; % view outside of Live Script window

 if exist('ax2','var')
 image(ax2,'CData',blendI); % change data only
 else % first time
 ax2 = gca;
 imshow(blendI,'Parent',ax2)
 end
 drawnow

 if we
 frame = getframe(f3); % capture with border frame
 % frame = getframe(gca); % take image only
 gr = frame2im(frame);
 [bA,map] = rgb2ind(gr,256);
 if (t==0)
 imwrite(bA,map,filename,'gif','LoopCount',Inf,'DelayTime',1/15);
 else
 imwrite(bA,map,filename,'gif','WriteMode','append','DelayTime',1/15);
 end
 end

end

imwrite(blendI,outputFile)
end

処理に多少時間が掛かる（私の環境で1フレーム1秒弱）ので、gifファイル書き込み時とそうでない時（weで制御）でstep, edの設定を変えています。

t2でハートの大きさが制御できます。
ライブスクリプトで、t2を数値スライダーでコントロール可能にしても良いですね。

gain でボケの明るさが制御できます。派手目にしたい場合は大きい値にしてみてください。

カーネル（ハート）は、あらかじめ輝度0のボーダーを含んだ中央部分に書かれるようサイズ調整してあります。

・処理結果

畳み込みだと処理にかなり時間が掛かるだろうから、とりあえず畳み込みで作って確認して、後で分割ブロックのFFTで高速化しようかとも思っていたのですが、なんと、畳み込みの方が速い！

まさに、最初の方に書いた、
「画像・カーネルサイズや処理フレームワークによってどちらが速いかは異なりますので、その都度確認する必要はあります」でした・・。

一応コードは両方載せておきました。カーネルサイズが小さいということもあるのでしょうが、最近のMATLABは余計なコトしない方が速いですね。（；￢＿￢）

カーネルのFFTはあらかじめ一回だけしておけば良いので、カーネルサイズが大きい場合等は、FFTした方が早くなる場合があると思います。

まあ、直線畳み込み・巡回畳み込み、MATLABでのブロック処理の参考にはなるかと。(^-^;

音声のような1次元データであれば、MATLABの DSP System Toolbox に一般的な Overlap-Save/Add 法が関数として用意されています。

Overlap-Save/Addの説明はMATLAB（DSP System Toolbox）のドキュメントの図が一番分かりやすいかと思いますので引用しておきます。

dsp.FrequencyDomainFIRFilter より引用

この関数を使うと自動的に周波数領域に変換して処理を行うため高速化が見込めます。また、普通にやると上記図において FFTLength – length(Numerator) + 1 のレイテンシーが生じますが、低レイテンシーモードも用意されています（同機能のSimulinkブロックもあります）。

冬は空気が澄んで、夜景が綺麗に撮れる季節です。

MATLABでさらにもうひと盛りしてませんか！？（￣ー￣）

編集履歴

2021/12/07
・「被写体の～信号処理的には被写体と絞りの「畳み込み演算」で再現できます。」部分追加
・実際のカメラ画像とシミュレーション画像比較を追加
・畳み込み版ソースコードのカーネルサイズが1ピクセル違っていたのを修正（結果に影響はありません）