陽炎のような画像処理効果を実装する

はじめに

今考えているおふざけネタを実現するために、陽炎のような画像エフェクトが必要となりました。この記事では、とりあえず陽炎エフェクトのMATLAB実装について書きました。

こんなん実装している人いるだろうと思って本当に軽く検索してみたのですが、アプリとかでは見つかるものの、実装例については見当たりませんでした。もっと深く検索していけば絶対見つかるとは思いますが。そこで自分で作ろうと思ったわけです。なんとなーくそれっぽい結果が出るように作っただけなので、大気中の屈折など物理モデルに基づいたものでは一切ございません。本気の実装例が見たい人はスルーしてくださると助かります。あくまでなんとなくです。

成果物

こんなんができるようになりました。真ん中の部分だけ効果を挿入しています。

そもそも陽炎ってどんなでしたっけ

こちらの動画を勝手に見つけて観察しました。

現象としては、元の像がランダムに上下左右に揺れ動く、みたいな感じでしょうか。ランダムに、と書いてますが物理的に見ていくと決してランダムではなさそうな気もしますが、とりあえずランダムとして考えます。

実装方法を考える

ランダムに上下左右に揺れ動くので、XY方向への変位場を表す配列を作成して値はランダムにセット、imwarpで画像変換したら完成だぜ！くらいに最初は考えました。が、時間方向の連続性を確保しないと1フレーム1フレーム全くランダムな変換となり、とても不自然になることが予想されます。

では、あらかじめ3次元の配列（XYZ）を用意して置き、そこにランダム値をセット、1方向（Z方向）について内挿してあげれば、時間に対して滑らかな変位場の配列が得られるなと思いました。しかし、この方法だと用意した配列分の時間しか変換できません。もちろん、新たな配列を作成して継ぎ足してあげる方法もありますが、継ぎ目の管理が面倒な気がします。

そもそも1フレーム1フレームの管理を簡単にするために、変位の変化分をランダムに作って、どんどん足しこんでいくという方法がよさそうです。しかし、それだと陽炎の滑らかな感じが出なそうですし、そもそもランダムウォークと同じなので、次第に値が発散してしまいそうです。

だらだらと書きましたが、最終的にはランダムな値を足しこんであげる、というところはそのままに、その値を位相として活用して、三角関数で変換する、という方法に落ち着きました。この方法なら、毎フレームごとにやることはランダムに位相を足す、だけで済み、永遠に更新できますし、-1~+1の間をさまようことになるので、発散もしません。とりあえず完璧。ちょとと文章力ないので、何言ってるかわからない人も多いと思いますが、コードで見ていきましょう。

実装（まず1点で考える）

画像の各画素の変位場を作る必要がありますが、まずは超シンプルに1点与えられたときに、その点の位置がゆらゆら揺れ動く、というものを実装します。X方向への変位とY方向への変位を計算して、元の値に足して点を再描画すればいいです。

Code

clear
t = 0; % 使わないけど一応書く。初期時間
dt = 0.03; % 時間の変化分
fp = 1; % 位相の変化周波数
fa = 0.5; % 振幅の変化周波数

% 初期位相
phase_theta_x = rand(1,1)*2*pi;
phase_theta_y = rand(1,1)*2*pi;
% 振幅変化のための初期位相
amp_theta_x = rand(1,1)*2*pi;
amp_theta_y = rand(1,1)*2*pi;

% 描画
p = plot(0,0,'ro'); axis([-1 1 -1 1]);

for i = 1:500
    % 位相の変化分をランダム値で圧縮する
    theta_rnd = rand(1,2); 
    amp_rnd = rand(1,2);
    % 現在位相の計算
    phase_theta_x = phase_theta_x + 2*pi*fp*dt*theta_rnd(1);
    phase_theta_y = phase_theta_y + 2*pi*fp*dt*theta_rnd(2);
    amp_theta_x = amp_theta_x + 2*pi*fa*dt*amp_rnd(1);
    amp_theta_y = amp_theta_y + 2*pi*fa*dt*amp_rnd(2);

    % ポイントの位置
    x = sin(amp_theta_x) * sin(phase_theta_x);
    y = sin(amp_theta_y) * sin(phase_theta_y);
    % 描画に反映
    p.XData = x;
    p.YData = y;
    pause(0.05)
end

結果の絵は省きますが、実装自体はしっかりできていることはわかりますが、陽炎っぽくないな~、と思いつつとりあえず先に進めることにします。1点なんでね、イメージが湧きづらいです。

2次元に拡張する

シンプルに拡張します。コメントなど揃っていませんが…。

Code

clear
t = 0;
dt = 0.01;
fp = 1;
fa = 0.5;
Sz = [30,30]; % 制御する点の数
[ptsX,ptsY] = meshgrid([0:Sz(2)-1],[0:Sz(1)-1]);

% 初期化
phase_theta_x = rand(Sz)*2*pi;
phase_theta_y = rand(Sz)*2*pi;
amp_theta_x = rand(Sz)*2*pi;
amp_theta_y = rand(Sz)*2*pi;

% 描画
fig = figure;
fig.Visible = 'on'
p = plot(ptsX(:),ptsY(:),'r.'); axis([-1 Sz(2)+1 -1 Sz(1)+1]);
for cnt = 1:100
    % 増加分をランダムに生成
    theta_rnd = ones([Sz,2]);%rand(1,2);
    amp_rnd = rand([Sz,2]);
    % 位相
    phase_theta_x = phase_theta_x + 2*pi*fp*dt*theta_rnd(:,:,1); % ランダムにdt分の変化を減らす
    phase_theta_y = phase_theta_y + 2*pi*fp*dt*theta_rnd(:,:,2); % ランダムにdt分の変化を減らす
    amp_theta_x = amp_theta_x + 2*pi*fa*dt*amp_rnd(:,:,1); % ランダムにdt分の変化を減らす
    amp_theta_y = amp_theta_y + 2*pi*fa*dt*amp_rnd(:,:,2); % ランダムにdt分の変化を減らす
    % 元の座標に変化分を足す
    x = ptsX + 0.5*sin(amp_theta_x) .* sin(phase_theta_x);
    y = ptsY + 0.5*sin(amp_theta_y) .* sin(phase_theta_y);
    % 描画に反映
    p.XData = x(:);
    p.YData = y(:);
    drawnow;
end

結果はこんな感じになります。

全然陽炎の動きっぽくないですね。となり合う点同士で動きの相関が無いためだと思いました。

内挿を含める

なのでそのあたりを考慮した実装にしてあげます。具体的にはX,Yの変位量を内挿してあげただけです、はい。

Code

clear
t = 0;
dt = 0.01;
fp = 1;
fa = 0.5;
Sz = [10,10];
% 内挿元の座標
[ptsX,ptsY] = meshgrid([0:Sz(2)-1],[0:Sz(1)-1]);
% 内挿先の座標
[ptsX_int,ptsY_int] = meshgrid([0:1/3:Sz(2)-1],[0:1/3:Sz(1)-1]);

% 初期化
phase_theta_x = rand(Sz)*2*pi;
phase_theta_y = rand(Sz)*2*pi;
amp_theta_x = rand(Sz)*2*pi;
amp_theta_y = rand(Sz)*2*pi;
% 描画
fig = figure;
fig.Visible = 'on'
p = plot(ptsX_int(:),ptsY_int(:),'r.'); axis([-1 Sz(2)+1 -1 Sz(1)+1]);
for cnt = 1:300

    theta_rnd = ones([Sz,2]);%rand(1,2);
    amp_rnd = rand([Sz,2]);
    phase_theta_x = phase_theta_x + 2*pi*fp*dt*theta_rnd(:,:,1);
    phase_theta_y = phase_theta_y + 2*pi*fp*dt*theta_rnd(:,:,2);
    amp_theta_x = amp_theta_x + 2*pi*fa*dt*amp_rnd(:,:,1);
    amp_theta_y = amp_theta_y + 2*pi*fa*dt*amp_rnd(:,:,2);

    x = ptsX + 0.5*sin(amp_theta_x) .* sin(phase_theta_x);
    y = ptsY + 0.5*sin(amp_theta_y) .* sin(phase_theta_y);
    % 内挿
    x_int = interp2(ptsX,ptsY,x,ptsX_int,ptsY_int,'cubic');
    y_int = interp2(ptsX,ptsY,y,ptsX_int,ptsY_int,'cubic');

    p.XData = x_int(:);
    p.YData = y_int(:);
    drawnow;
end

う~ん、どうなんでしょう。陽炎ぽい…かなぁ。まぁなんとなくできてればいいか、ということで先に進めます。

システムオブジェクト化する

1フレームずつ呼び出すと、最新の変位場行列が得られる、みたいにすれば使いやすくなります。こういう場合は、システムオブジェクトが使いやすいです。下記実装例です。私もあまり作り慣れていないので、特に最初の初期設定のところとか、もっとこうした方がいいなどアドバイスいただけるとありがたいです。

Code

classdef kageroDeformationMap < matlab.System
    properties
        t = 0;
        dt = 0.1;
        fp = 1;
        fa = 0.5;
        Sz;
        dataSz;
        ptsX
        ptsY
        ptsX_int
        ptsY_int

        phase_theta_x
        phase_theta_y
        amp_theta_x
        amp_theta_y
    end

    properties(DiscreteState)

    end

    % Pre-computed constants
    properties(Access = private)

    end

    methods
        function obj = kageroDeformationMap(varargin)
            setProperties(obj,nargin,varargin{:})
        end
    end

    methods(Access = protected)
      function setupImpl(obj)
            [obj.ptsX,obj.ptsY] = meshgrid(0:10:obj.Sz(2),0:10:obj.Sz(1));
            [obj.ptsX_int,obj.ptsY_int] = meshgrid(0:1:obj.Sz(2),0:1:obj.Sz(1));

            obj.dataSz = size(obj.ptsX);

            obj.phase_theta_x = rand(obj.dataSz)*2*pi;
            obj.phase_theta_y = rand(obj.dataSz)*2*pi;
            obj.amp_theta_x = rand(obj.dataSz)*2*pi;
            obj.amp_theta_y = rand(obj.dataSz)*2*pi;
        end

        function [x_int,y_int] = stepImpl(obj)
            theta_rnd = ones([obj.dataSz,2]);%rand(1,2);
            amp_rnd = rand([obj.dataSz,2]);
            obj.phase_theta_x = obj.phase_theta_x + 2*pi*obj.fp*obj.dt*theta_rnd(:,:,1); % ランダムにdt分の変化を減らす
            obj.phase_theta_y = obj.phase_theta_y + 2*pi*obj.fp*obj.dt*theta_rnd(:,:,2); % ランダムにdt分の変化を減らす
            obj.amp_theta_x = obj.amp_theta_x + 2*pi*obj.fa*obj.dt*amp_rnd(:,:,1); % ランダムにdt分の変化を減らす
            obj.amp_theta_y = obj.amp_theta_y + 2*pi*obj.fa*obj.dt*amp_rnd(:,:,2); % ランダムにdt分の変化を減らす

%             x = obj.ptsX + 0.5*sin(obj.amp_theta_x) .* sin(obj.phase_theta_x);
%             y = obj.ptsY + 0.5*sin(obj.amp_theta_y) .* sin(obj.phase_theta_y);
            x = 5*sin(obj.amp_theta_x) .* sin(obj.phase_theta_x);
            y = 5*sin(obj.amp_theta_y) .* sin(obj.phase_theta_y);

            x_int = interp2(obj.ptsX,obj.ptsY,x,obj.ptsX_int,obj.ptsY_int,'cubic');
            y_int = interp2(obj.ptsX,obj.ptsY,y,obj.ptsX_int,obj.ptsY_int,'cubic');
        end

        function resetImpl(obj)
            obj.phase_theta_x = rand(obj.dataSz)*2*pi;
            obj.phase_theta_y = rand(obj.dataSz)*2*pi;
            obj.amp_theta_x = rand(obj.dataSz)*2*pi;
            obj.amp_theta_y = rand(obj.dataSz)*2*pi;
        end
    end
end

使い方は読み進めていただけるとわかるかと思います。

画像を準備する

さて、準備ができましたので、適当なテスト画像に対して試してみましょう。

Code

img = imread('peppers.png');
imshow(img);

ROIを設定する

今回は、1枚の画像の中で、所望の箇所にエフェクトを付けてあげるようにします。drawpolygonでだいたい真ん中あたりに適当に四角を作成しました。その範囲内だけに陽炎効果を付加します。

Code

h = drawpolygon; % 
% ポリゴン範囲の座標を取得
minX = floor(min(h.Position(:,1)));
maxX = ceil(max(h.Position(:,1)));
minY = floor(min(h.Position(:,2)));
maxY = ceil(max(h.Position(:,2)));
% createMask関数で領域を含む最小限のマスクを作成するために、オフセットして左上に動かす
h.Position(:,1) = h.Position(:,1)-minX+1;
h.Position(:,2) = h.Position(:,2)-minY+1;
msk = createMask(h,maxY-minY+1,maxX-minX+1); % 最小限のマスク
figure;imshow(msk)

エフェクトを計算する

10pxごとに変位を計算する仕様になっているので、システムオブジェクトに入れるサイズを求めてあげます。

Code

% エフェクトの計算範囲をmskから適当に決める
effectSz = ceil(size(msk)/10) * 10; % 10pxごとに変位を計算（最終的には内挿して1pxごととする）

% 変位場を生成するシステムオブジェクト生成
defmap = kageroDeformationMap('Sz',effectSz); % インスタンス生成
[defX,defY] = defmap(); % 変位場を生成

変位場を配列に

出てきた変位場の配列を元画像と同じサイズの配列の中にはめ込んであげます。ROIの外は変位はゼロです。

Code

defX = defX(1:size(msk,1),1:size(msk,2));
defY = defY(1:size(msk,1),1:size(msk,2));
defX(~msk) = 0;
defY(~msk) = 0;
defXAll = zeros(size(img,[1:2]),'single');
defXAll(minY:minY+size(msk,1)-1,minX:minX+size(msk,2)-1) = defX;
defYAll = zeros(size(img,[1:2]),'single');
defYAll(minY:minY+size(msk,1)-1,minX:minX+size(msk,2)-1) = defY;
imshow(defXAll,[])

画像変換

Code

warpimg = imwarp(img,cat(3,defXAll,defYAll));
imshow(warpimg)

それっぽくできていそうです。それでは以上の内容をループで繰り返します。

Code

fig = figure;
fig.Visible = 'on';
axis tight;
for cnt = 1:30
    [defX,defY] = defmap(); % 変位場を生成
    defX = defX(1:size(msk,1),1:size(msk,2));
    defY = defY(1:size(msk,1),1:size(msk,2));
    defX(~msk) = 0;
    defY(~msk) = 0;
    defXAll = zeros(size(img,[1:2]),'single');
    defXAll(minY:minY+size(msk,1)-1,minX:minX+size(msk,2)-1) = defX;
    defYAll = zeros(size(img,[1:2]),'single');
    defYAll(minY:minY+size(msk,1)-1,minX:minX+size(msk,2)-1) = defY;
    % imshow(defXAll,[])
    warpimg = imwarp(img,cat(3,defXAll,defYAll));
    imshow(warpimg)
    drawnow;
end

出てきた画像が冒頭のものとなります。まぁ陽炎っぽいといえば陽炎っぽいでしょうか。もうちょっとX、Yで変位の最大値などを調整すればさらにそれっぽくできそうですが、今日はこれまで。

終わりに

またアドベントカレンダー後半にもう1記事書きます。そのときにこちらを使います。お楽しみに？

謝辞

本記事は @eigs さんのlivescript2markdownを使わせていただいてます。