纏を習得する

# はじめに


なんのこっちゃ、と思われた方はこの記事をスルーしていただいてもいいかもしれません。纏（てん）とはハンターハンターにおける念能力の基本技（？）です。作中においては次のように説明されています。




> 体内の精孔からあふれ出ている生命エネルギー（オーラ）を肉体の周りにとどめる！！




漫画の中には描写もあるので、今日はこれを私なりに解釈して実装してみようと思った次第です。アニメ版見れば本当はどんな感じかわかると思うのですが、たぶん私の実装とは違うと思うので見て確認はしていません…。

# 動作環境
MATLABのバージョンについてはR2020b。
また使用しているツールボックスは下記になります。
・Image Processing Toolbox
・Computer Vision Toolbox
・Deep Learning Toolbox

# 成果物


こんな感じになります。




![test.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/b6776fcf-b17b-86ee-59b1-6447a7a72a1b.gif)


# 実装


流れとしては簡単です。



   1.  セマンティックセグメンテーションで人物をセグメンテーション 
   1.  後処理で、人物の周囲の範囲を取得 
   1.  その範囲だけ陽炎のように歪ませる。 



陽炎のように歪ませる、というあたりが私なりの纏の解釈です。オーラが体の周りにとどまるイメージ。




一つずつ見ていきたいと思います。


## セマンティックセグメンテーションで人物をセグメンテーション


セマンティックセグメンテーションモデルを使うだけです。ここはすみません、本当は読者の方々が再現できるようオープンなモデルを使った方がいいのはわかってるのですが、ありもので社内のリソースを使いました。人物のセマンティックセグメンテーションモデル自体はありふれていると思いますので、もしご興味あれば適当にググってインポートもできるかと思います。



```matlab:Code
% モデルの読み込み
load trainedNet.mat
load classNames.mat
img = imread('Shizentai.jpg');
img = imresize(img, 400./min(size(img,1:2))); % モデルが要求する最低サイズ[400 400]に合わせてリサイズ

C = semanticseg(img,net);
msk = C == 'person'; % 人物部分のマスクを作成
imshowpair(img,msk);
```


![figure_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/1e7e9dcc-3cbf-463d-d5df-a6eb6acab16c.png)




少しだけ太目に判定されるので目視でいい感じに絞ります。あと、穴をふさぎます。このあたりはモデルの精度如何かと。ちなみにいつまでたってもdilateとerodeが、どっちが膨張でどっちが収縮か覚えられません。



```matlab:Code
msk2 = imerode(msk,strel('disk',3));
msk2 = imfill(msk2,'holes');
imshowpair(img,msk2)
```


![figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/7b0c04e3-f748-8ce3-7667-0c727f48fa2c.png)


## 後処理で、人物の周囲の範囲を取得


必要なのは人物の周囲の範囲なので、適当に膨張させ、元のマスクを引いて取得します。また、外側にいくにつれて効果が薄くなるように傾斜をつけます。



```matlab:Code
msk3 = imdilate(msk2,strel('disk',61));
msk_around = imgaussfilt(single(msk3),4);
msk_around(msk2) = 0;
imshowpair(img,msk_around)
```


![figure_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/d2c4fd7b-cb2d-cab0-a8f4-1c9adf75edb6.png)



この画像はこの画像ですでにそれっぽいです。


## その範囲だけ陽炎のように歪ませる


こちらは[前記事](https://qiita.com/shun-kusano/items/42c69a753153d0d45018)の成果物を使います。いくつか変更点（バグ修正も含めて）があります。前回はROIの範囲を指定して、そのサイズの変位場行列を作成していました。ですが、この方法だとフレームごとにサイズが変わったときに処理が面倒になります。そこで今回は、画像全体で行列を作成してしまって、必要なROIで切り抜く形にしました。変位場行列作成自体の実装は変わりませんが、前後の処理が少し変わったことになります。あと、事前割当の配列サイズおよび内挿後の配列サイズが間違っていたので直しました。さらに、徐々に変位量を大きくするような仕様にしています。




さらに、陽炎ってぼやけて見えるよな、ってことに思い至ったため、画像のぼやけも追加します。さらに、より分かりやすくするため、色味も加えます。個人的に青色をましましします。「さらに」だらけで申し訳ないです。



```matlab:Code
% 前処理としてぼかしと青味を増す
imEffectedAll = imgaussfilt(img,4); % ぼかす
imEffectedAll(:,:,3) = imEffectedAll(:,:,3)*1.5; % 青味を増す
imEffectedROI = single(img).*(1-msk_around) + ...
    single(imEffectedAll).*(msk_around);
imEffectedROI = uint8(imEffectedROI);
imshow(imEffectedROI)
```


![figure_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/13a1c9ef-af43-5305-6d51-e42aeb858383.png)


```matlab:Code
% 変位場生成オブジェクトを生成
im_Size = size(img,1:2);
defMap = kageroDeformationMap('Sz',im_Size);

fig = figure;
fig.Visible = 'on';

for i=1:100
    [defX,defY] = defMap();
    defX = defX .* msk_around;
    defY = defY .* msk_around;
    warp_img = imwarp(imEffectedROI,cat(3,defX,defY));
    imshow(warp_img);
    drawnow;
end
```


![figure_4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/acb7841a-aff2-e9b4-5ee2-9b58bb38e956.png)




オーラがとどまっているというか、何か妖気のようなものが漏れ出ている感じになっている気もしますが、気にせず進めます。走り切ることが大事。


# 纏の発動条件


これだけだと、常に纏状態となってしまいます。よりリアルにするために発動条件が必要です。作中では自然体が最も適している、みたいなことが書いてあったような気がするので、画像中の人物が自然体になったら、という条件をみたしたら発動、ということにします。




この「自然体になったら」のジャッジはセマンティックセグメンテーションでは難しいので、骨格検出のモデルも使うことにします。骨格検出モデルでは、人体の要所座標をキーポイントとして得ることができますので、姿勢判別には特に有効です。




あまり時間もなかったので、自然体かどうかの判定としては、



   -  肩、肘、手首までがほぼ一直線に真下に下がっている 
   -  腰、膝、足首までがほぼ一直線に真下に下がっている 



こととしました。本当は顔の向きとかも含めるべきなんでしょうが、簡易的に済ませます。モデルには[こちら](https://github.com/matlab-deep-learning/Human-Pose-Estimation-with-Deep-Learning)を使わせていただきました。



```matlab:Code
detector = posenet.PoseEstimator;
Iin = imresize(img,detector.InputSize(1:2));
keypoints = detectPose(detector,Iin);
J = detector.visualizeKeyPoints(Iin,keypoints);
imshow(J);
```


![figure_5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202697/5b8cafbb-e089-1efd-22dc-5572b78d95f7.png)



```matlab:Code
% 自然体の判別
shldrL = keypoints(6,1:2)';
shldrR = keypoints(7,1:2)';
elbwL = keypoints(8,1:2)';
elbwR = keypoints(9,1:2)';
wrstL = keypoints(10,1:2)';
wrstR = keypoints(11,1:2)';
wstL = keypoints(12,1:2)';
wstR = keypoints(13,1:2)';
knL = keypoints(14,1:2)';
knR = keypoints(15,1:2)';
anklL = keypoints(16,1:2)';
anklR = keypoints(17,1:2)';
% パーツごとにまとめる
leftT = [shldrL elbwL wrstL]; % 左腕
leftB = [wstL knL anklL]; % 左脚
rightT = [shldrR elbwR wrstR]; % 右腕
rightB = [wstR knR anklR]; % 右脚
% 差分を求めて角度を出す
dif_leftT = diff(leftT,1,2);
dif_leftB = diff(leftB,1,2);
dif_rightT = diff(rightT,1,2);
dif_rightB = diff(rightB,1,2);

ang_leftT = atan2(dif_leftT(2,:),dif_leftT(1,:));
ang_leftB = atan2(dif_leftB(2,:),dif_leftB(1,:));
ang_rightT = atan2(dif_rightT(2,:),dif_rightT(1,:));
ang_rightB = atan2(dif_rightB(2,:),dif_rightB(1,:));
% 全てsinで変換したときにth以上
th = 0.9;
if all(sin(ang_leftT)>th) ...
        && all(sin(ang_leftB)>th) ...
        && all(sin(ang_rightT)>th) ...
        && all(sin(ang_rightB)>th)
    isTen = true;
else
    isTen = false;
end
```


```text:Output
isTen = 
   1

```

# 最終的な実装


ここまで書いていてなんですが、最終的には青味やボケが徐々に強くなるように、その部分をシステムオブジェクト化したり、纏の判定に姿勢＋その時間安定性を加えたり、といったことをしています。もし興味があればご覧ください。これで冒頭の動画ができました。

<details><summary>コード全文</summary>
<div>

```matlab:Ten_run.mlx
clear, close, clc;
% 骨格推定モデルの読み込み
detector = posenet.PoseEstimator;
% セマンティックセグメンテーションモデルの読み込み
load trainedNet.mat
load classNames.mat

% 1フレーム目を読み込み
vr = vision.VideoFileReader('MOV_0140.mp4');
img = vr();
img = uint8(img*255);
% img = imread('Sizentai.jpg');
img = imresize(img, 400./min(size(img,1:2))); % セマンティックセグメンテーションモデルが要求する最低サイズ[400 400]に合わせてリサイズ

im_Size = size(img,1:2);
% 変位場行列を生成するためのオブジェクト
defMap = kageroDeformationMap('Sz',im_Size);
% キーフレームの平滑化のためのオブジェクト
smoother = keyframeSmoother('wb',9); % 後ろ９フレームで平滑化
% キーフレームが安定化したかを確認するためのオブジェクト
judger = stabilizeJudger('FrameNum',30); % 5フレーム連続で動きがなければ安定とみなす
% 画像ぼかし＋青みがけのためのオブジェクト
effector = imageEffector('TransitionFrameNum',50); % 50フレームかけてエフェクトを徐々につけていく

fig = figure;
fig.Visible = 'on';
ax = gca;

while ~isDone(vr)
    img = vr();
    img = uint8(img*255);
%    img = imread('Sizentai.jpg');
    img = imresize(img, 400./min(size(img,1:2))); % セマンティックセグメンテーションモデルが要求する最低サイズ[400 400]に合わせてリサイズ
    
    % 骨格推定による自然体判定
    Iin = imresize(img,detector.InputSize(1:2));
    keypoints = detectPose(detector,Iin);
    sm_keypoints = smoother(keypoints);
    isStabilized = judger(keypoints);
    J = detector.visualizeKeyPoints(Iin,sm_keypoints);
%    imshow(J,'Parent',ax)
    
    shldrL = sm_keypoints(6,1:2)';
    shldrR = sm_keypoints(7,1:2)';
    elbwL = sm_keypoints(8,1:2)';
    elbwR = sm_keypoints(9,1:2)';
    wrstL = sm_keypoints(10,1:2)';
    wrstR = sm_keypoints(11,1:2)';
    wstL = sm_keypoints(12,1:2)';
    wstR = sm_keypoints(13,1:2)';
    knL = sm_keypoints(14,1:2)';
    knR = sm_keypoints(15,1:2)';
    anklL = sm_keypoints(16,1:2)';
    anklR = sm_keypoints(17,1:2)';
    % パーツごとにまとめる
    leftT = [shldrL elbwL wrstL]; % 左腕
    leftB = [wstL knL anklL]; % 左脚
    rightT = [shldrR elbwR wrstR]; % 右腕
    rightB = [wstR knR anklR]; % 右脚
    % 差分を求めて角度を出す
    dif_leftT = diff(leftT,1,2);
    dif_leftB = diff(leftB,1,2);
    dif_rightT = diff(rightT,1,2);
    dif_rightB = diff(rightB,1,2);
    
    ang_leftT = atan2(dif_leftT(2,:),dif_leftT(1,:));
    ang_leftB = atan2(dif_leftB(2,:),dif_leftB(1,:));
    ang_rightT = atan2(dif_rightT(2,:),dif_rightT(1,:));
    ang_rightB = atan2(dif_rightB(2,:),dif_rightB(1,:));
    % 全てsinで変換したときにth以上
    th = 0.9;
    if all(sin(ang_leftT)>th) ...
            && all(sin(ang_leftB)>th) ...
            && all(sin(ang_rightT)>th) ...
            && all(sin(ang_rightB)>th) ...
            && isStabilized
        isTen = true;
    else
        isTen = false;
    end
    
    if isTen % 纏なら
        C = semanticseg(img,net);
        msk = C == 'person'; % 人物部分のマスクを作成
        msk2 = imerode(msk,strel('disk',3));
        msk2 = imfill(msk2,'holes');
        
        % 人物の周囲を取得
        msk3 = imdilate(msk2,strel('disk',61));
        msk_around = imgaussfilt(single(msk3),5);
        msk_around(msk2) = 0;
        
        % 前処理としてぼかしと青味を増す
        imEffectedAll = effector(img);
        imEffectedROI = single(img).*(1-msk_around) + ...
            single(imEffectedAll).*(msk_around);
        imEffectedROI = uint8(imEffectedROI);
        
        % 変位場を生成し、歪ませる
        [defX,defY] = defMap();
        defX = defX .* msk_around;
        defY = defY .* msk_around;
        warp_img = imwarp(imEffectedROI,cat(3,defX,defY));
        out_img = insertText(warp_img, [10 10], "Ten", "FontSize", 30, "BoxColor", "blue");
    else % そうでなければ
        reset(defMap);
        reset(effector);
        out_img = img;
    end
    
    imshow(out_img,'Parent',ax);
    drawnow;
    

end

delete(smoother);
delete(effector);
delete(defMap);
delete(judger);

```

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

        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.time = 0;
            [obj.ptsX,obj.ptsY] = meshgrid(0:10:ceil(obj.Sz(2)/10)*10,0:10:ceil(obj.Sz(1)/10)*10);
            [obj.ptsX_int,obj.ptsY_int] = meshgrid(linspace(0,obj.Sz(2),obj.Sz(2)),linspace(0,obj.Sz(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)
            obj.time = obj.time + 1;
            
            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 = min(5,5*obj.time/50)*sin(obj.amp_theta_x) .* sin(obj.phase_theta_x);
            y = min(5,5*obj.time/50)*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;
            obj.time = 0;
        end
    end
end
```

```matlab:keyframeSmoother
classdef keyframeSmoother < matlab.System
    properties
        wb
        KeyFrames
    end
    
    methods
        function obj = keyframeSmoother(varargin)
            setProperties(obj,nargin,varargin{:})
        end
    end
    
    methods(Access = protected)
        function setupImpl(obj)
            obj.KeyFrames = nan(17,3,obj.wb+1);
            % Perform one-time calculations, such as computing constants
        end

        function output = stepImpl(obj,keyframe)
            obj.KeyFrames(:,:,1) = keyframe;
            sm_data = smoothdata(obj.KeyFrames,3,'rloess',[obj.wb 0],'omitnan');

            obj.KeyFrames = circshift(obj.KeyFrames,1,3);
            output = sm_data(:,:,1);
        end

        function resetImpl(obj)
            obj.KeyFrames = nan(17,3,obj.wb+1);
        end
    end
end
```

```matlab:imageEffector
classdef imageEffector < matlab.System
    properties
        CurrentTransitionDepth
        TransitionFrameNum
    end
    
    methods
        function obj = imageEffector(varargin)
            setProperties(obj,nargin,varargin{:})
        end
    end
    
    methods(Access = protected)
        function setupImpl(obj)
            obj.CurrentTransitionDepth = 0;
        end

        function imEffectedAll = stepImpl(obj,img)
            obj.CurrentTransitionDepth = obj.CurrentTransitionDepth + 1;
            progress = obj.CurrentTransitionDepth/obj.TransitionFrameNum; 
            imEffectedAll = imgaussfilt(img,4*min(1,progress)); % ぼかす
            imEffectedAll(:,:,3) = imEffectedAll(:,:,3)*(1+0.5*min(1,progress)); % 青味を増す
        end

        function resetImpl(obj)
            obj.CurrentTransitionDepth = 0;
        end
    end
end
```

```matlab:stabilizeJudger
classdef stabilizeJudger < matlab.System
    properties
        FrameNum
        KeyFrames
    end
    
    methods
        function obj = stabilizeJudger(varargin)
            setProperties(obj,nargin,varargin{:})
        end
    end
    
    methods(Access = protected)
        function setupImpl(obj)
            obj.KeyFrames = zeros(17,3,obj.FrameNum);
        end

        function isStabilized = stepImpl(obj,keyframe)
            obj.KeyFrames(:,:,1) = keyframe;
            med_keyframe = median(obj.KeyFrames,3);
            
            diffs = abs(med_keyframe - obj.KeyFrames);
            
            bw = diffs < 20; % less than 20px
            if any(~bw(:))
                isStabilized = false;
            else
                isStabilized = true;
            end
            obj.KeyFrames = circshift(obj.KeyFrames,1,3);
        end

        function resetImpl(obj)
            
            
        end
    end
end
```

</div></details>


# 終わりに


そもそもなんでこんなことしようと思ったかなのですが、時をさかのぼること何年か前、前職でくすぶっていたときに[@takaya901](https://qiita.com/takaya901)さんの[こちらの記事](https://qiita.com/takaya901/items/9f4c87a3d51ea019520f)に妙に感銘を受けてしまいました。見る人が見たらくだらないかもしれませんが、漫画やアニメの世界観を技術で再現するというコンセプトが妙に心を打って「エンジニアリングってこんなに自由なんだな」と思いました。同時に、「そもそも作りたい未来を実現できるからエンジニアに魅かれたのに、日々の業務でそれすら忘れてしまっていたのか俺は」と。その時から、自分もこういうの何か作ってみたいなと思っていたのですが、またまた日々の業務に追われてしまっておりました。この機会にやっと念願果たせてうれしく思っています。




本当は「練」とか「凝」とかも実装したかったのですが、すみません、このタイミングで引っ越ししまして、Qiita書いてる場合でもなく、最低限の纏で妥協しました。

# 謝辞


本記事は [@eigs](https://qiita.com/eigs) さんの[livescript2markdown](https://github.com/minoue-xx/livescript2markdown)を使わせていただいてます。