ストロークを平均化～ MATLAB の figure に callback を設定してお絵かき

leftbank

2022年8月27日 09:33

数式を平均化

フーリエ級数展開

任意の曲線は、（一定の条件下で）フーリエ級数展開によって sin/cos の和で表せます。

それを応用した、「手書きひらがなの平均を求める」、というのがありました。

みんなの手書きの「あ」を集めて平均するとどうなる？

ひらがなの平均手書き文字は綺麗

面白いな、と思ってMATLABでちょっと試してみたので解説を。

論文では、ストロークデータを折り返して閉曲線にしたり、元との差が十分小さくなる次数まで展開したりしています。

それならばDCTでやってしまった方が処理が楽そうなのでそうします。
（DCTは、入力を偶関数化してDFTしたのと等価ですし、可逆変換です）

メインの処理は数分で書けましたが、GUIのコールバックをどうすればシンプルでうまく動作するかはかなり試行錯誤しました。(^^;)

一画バージョン

まずは簡単化のために、「一画」だけに限定してみます。

マウスインタラクティブ

figureには、マウス／キーボード入力に関するコールバック関数がいくつかあります。それぞれどの動作で呼び出されるか決まっており、そこにその時の動作内容を記述することで、場合に応じた処理を自動的に行うことができます。

ライブエディター内の figure では使えないので、ライブスクリプトでは以下のようにして外に出す必要があります。

f=figure;
f.Visible = 'on';

（今回のスクリプトでは 'waitforbuttonpress' で外に出しています）

マウスに関する callback は以下の3つで、今回は全て使います。

WindowButtonDownFcn：マウスクリック時に呼び出される
WindowButtonMotionFcn：マウスを移動したときに呼び出される
WindowButtonUpFcn：クリックを離したときに呼び出される

コールバック関数の定義方法としては、以下のように直接書く方法、

f = figure(1);
f.WindowButtonDownFcn = 'pos = ax.CurrentPoint(1,1:2);';
f.WindowButtonMotionFcn = 'movFlag = true;';
f.WindowButtonUpFcn = 'start = false;';

今回は処理的に不要ですが、無名関数を使って直接書く方法、

ユーザー定義関数を呼び出す方法があります。

ユーザー定義関数を呼び出すには、以下のようにします。

f.WindowButtonDownFcn = @userFunc;

function userFunc(src,event)
    % 実際の処理
end

ユーザー定義関数を使った方が処理としてはスッキリしますが function化が必要で、function化してしまうとライブスクリプトでは動かないようなので、通常のMATLABスクリプトとして動かすことになります。

ライブスクリプトの方が、ファイル保存しなくても動くし、ワークスペースに変数が入るし、とにかく色々と便利なので、スクリプトじゃないとダメなとき以外、最近はほぼライブスクリプトを使っています。

キーボードスキャン

今回は使っていませんが、f.KeyPressFcn でキースキャンもできます。

公式ドキュメントには「KeyPressFcn コールバックはライブエディターではサポートされていません。」と書かれていますが、先ほど書いたように、figure を外に出せば使えます。

以下をコマンドラインで打ち込んで、適当にキーを押してみてください。
（無名関数で定義しているので、関数名は不要です）

f = figure;
f.KeyPressFcn = @(src,kdata) disp(kdata.Key);

実行結果例

t
q
w
space
escape
shift
capslock
leftarrow
rightarrow
uparrow
downarrow

このようにキー情報が取れます。'Key' は最後に押されたキー情報のみで、大文字小文字は区別されませんが装飾キーも取れます。普通にキャラクターとして取りたい場合は、'Character' を使います。

f.KeyPressFcn = @(src,kdata) disp(kdata.Character);

e
E

"Modifier" では、押されている複数の装飾キーが cell 配列で同時に取れます。

f.KeyPressFcn = @(src,kdata) disp(kdata.Modifier);

{'shift'}    {'control'}

スクリプト

ライブスクリプト・スクリプト共通版

ライブスクリプトを新規作成して、そこにペーストしてください。
通常のスクリプトとしても動きます。

clear
close all

f = figure(1);
f.Pointer = 'hand';  % change cursor shape
f.Color = [0.75 0.75 0.75];  % background color
ax = axes;  % create axes
axis equal  % square aspect ratio

DCTn = 32;  % DCT length

for k=1:5  % do 5 times drawing
    pos = ax.CurrentPoint(1,1:2);  % get cursor position
    p(k) = plot(ax, pos(:,1),pos(:,2), 'b', Linewidth=2);  % make strok object
    xlim(ax, [0 600]);  ylim(ax, [0 600]);
    set(ax,'xtick',[],'ytick',[])
    box(ax,'on')

    % set each callback
    f.WindowButtonDownFcn = 'pos = ax.CurrentPoint(1,1:2);';  % get 1st position
    f.WindowButtonUpFcn = 'start = false;';  % end of one strok
    f.WindowButtonMotionFcn = 'movFlag = true;';  % draw strok

    waitforbuttonpress  % wait for the each 1st click

    start = true;  movFlag = false;
    while start
        if movFlag
            pos = [pos; ax.CurrentPoint(1,1:2,1)];  % add current position to strok
            p(k).XData = pos(:,1);  % update the strok object
            p(k).YData = pos(:,2);  % update the strok object
            movFlag = false;
        end
        drawnow  % interrupt for event
    end

    % clear callbacks
    f.WindowButtonDownFcn = '';
    f.WindowButtonMotionFcn = '';
    f.WindowButtonUpFcn = '';

    n = length(pos);
    t0 = (1:n)';  % index of orginal strok data
    t = (0:n-1)'/(n-1);  % for interpolation input
    tt = 0:1/DCTn:1-1/DCTn;  % for interpolation output
    xx = spline(t,pos(t0,1),tt);  % transform data to suitable for DCT
    yy = spline(t,pos(t0,2),tt);  % transform data to suitable for DCT
    fx = dct(xx);  fy = dct(yy);  % get DCT coefficients

    if k == 1
        % set as mean coefficients
        fxm = fx;
        fym = fy;
    else
        % update running mean
        fxm = (fxm * (k-1) + fx) / k;
        fym = (fym * (k-1) + fy) / k;
    end

    hold(ax,'on')

    for pk = k-1:-1:1
        p(pk).Color = [p(pk).Color pk/k/4];  % decay transparancy of old strok
        pm(pk).Color = [pm(pk).Color 0];  % erase old mean strok
    end

    pm(k) = plot(ax, idct(fxm), idct(fym), 'r', Linewidth=2);  % draw mean strok
end

hold(ax,'off')
p(k).Color = [p(k).Color k/k/4];  % decay transparancy of the last strok
f.Pointer = 'arrow';  % restore cursor shape
disp('finished.')

MATLABはメモリの動的確保ができます。
一応このような警告が出ますが、ここはお行儀悪く初期化せずに使います。

平均のアップデート

平均のアップデートは、Simulink であれば Mean ブロックの "Running mean" をチェックするだけですが、関数にはその機能はないので自前で計算します。

とは言っても、平均の定義式を

$$
\mu_N=\displaystyle\sum_{i=1}^n x_i / N
$$

とすると、一つ前の平均 $${\mu_{N-1}}$$ と新たなデータ $${x_N}$$ から、

$$
\mu_N=((N-1) \mu_{N-1} + x_N) / N
$$

で求まりますね。

（Simulink の "Running mean" オプションは削除予定で、新たに機能拡張された "Moving Average" ブロックが用意されています）

標準偏差とかの場合はもう少し複雑になります。以下をご参照ください。

MATLAB Functionブロックを使おう！～ Running σ

透明度変更

1回目を描き終えると、過去のストロークの透明度を変えて薄くしています。

これはなぜか公式ドキュメントには書かれていないのですが、α チャネルを直接指定することで透明度が変えられます。

詳しくはこちらをご覧ください。

MATLAB で透明度付きプロット～円は rectangle で

動作例

動作例（一画バージョン）

クリックしながら一画だけ描いてクリックを離すと、スプライン補間を行い、さらに平均化されたストロークが赤で描かれます。

動画は最初に録ったので過去の平均ストロークも残っていますが、上記スクリプトでは消しています。

スプライン補間は、ある程度滑らかにするのと、その度違うポイント数になるストロークデータを、DCT長に合わせる目的の両方を兼ねています。

DCTした後はDCT係数をそれまでのストロークと平均化して、逆DCTして平均化ストロークとして描いています。

ローカル動作では動画のようにスムーズに動きますが、MATLAB Online の場合は早く動かすとデータが取れないようなので、軌跡がマウスに追いつくよう、ゆ～～～っくり動かしてみてください。

一応、スクリプト版も載せておきます。

スクリプト版

function strok_DCT
clear
close all

f = figure(1);
f.Pointer = 'hand';  % change cursor shape
f.Color = [0.75 0.75 0.75];  % background color
f.WindowButtonDownFcn = '';
f.WindowButtonMotionFcn = '';
f.WindowButtonUpFcn = '';

ax = axes;  % create axes
axis equal  % square aspect ratio

DCTn = 32;  % DCT length

for k=1:5  % do 5 times drawing
    pos = ax.CurrentPoint(1,1:2);  % get cursor position
    p(k) = plot(ax, pos(:,1),pos(:,2), 'b', Linewidth=2);  % make strok object
    xlim(ax, [0 600]);  ylim(ax, [0 600]);
    set(ax,'xtick',[],'ytick',[])
    box(ax,'on')
    hold(ax,'on')

    % set each callback
    f.WindowButtonDownFcn = @startStrok;
    f.WindowButtonUpFcn = @endStrok;

    waitforbuttonpress  % wait for the each 1st click

    start = true;
    while start
        drawnow  % interrupt for event
    end

    % clear callbacks
    f.WindowButtonDownFcn = '';
    f.WindowButtonMotionFcn = '';
    f.WindowButtonUpFcn = '';

    n = length(pos);
    t0 = (1:n)';  % index of orginal strok data
    t = (0:n-1)'/(n-1);  % for interpolation input
    tt = 0:1/DCTn:1-1/DCTn;  % for interpolation output
    xx = spline(t,pos(t0,1),tt);  % transform data to suitable for DCT
    yy = spline(t,pos(t0,2),tt);  % transform data to suitable for DCT
    fx = dct(xx);  fy = dct(yy);  % get DCT coefficients

    if k == 1
        % set as mean coefficients
        fxm = fx;
        fym = fy;
    else
        % update running mean
        fxm = (fxm * (k-1) + fx) / k;
        fym = (fym * (k-1) + fy) / k;
    end

    hold(ax,'on')

    for pk = k-1:-1:1
        p(pk).Color = [p(pk).Color pk/k/4];  % decay transparancy of old strok
        pm(pk).Color = [pm(pk).Color 0];  % erase old mean strok
    end

    pm(k) = plot(ax, idct(fxm), idct(fym), 'r', Linewidth=2);  % draw mean strok
end

hold(ax,'off')
p(k).Color = [p(k).Color k/k/4];  % decay transparancy of the last strok
f.Pointer = 'arrow';  % restore cursor shape
disp('finished.')


    function startStrok(src,event)
        pos = ax.CurrentPoint(1,1:2);
        f.WindowButtonMotionFcn = @drawStrok;
    end

    function drawStrok(src,event)
        pos = [pos; ax.CurrentPoint(1,1:2)];
        % add current position to strok
        p(k).XData = pos(:,1);  % update the strok object
        p(k).YData = pos(:,2);  % update the strok object
    end

    function endStrok(src,event)
        start = false;
        f.WindowButtonMotionFcn = '';
    end
end

strok_DCT.m として保存して実行してください。

複数画（かく）対応

次に、ライブスクリプト版を、複数画（かく）対応にしてみます。

ストロークのポイント数は当然毎回違うので、大きさの違う行列をまとめて扱える cell 配列を使うと便利です。

このように、違うサイズの行列をまとめて扱えます。
変数の型が違ってもOKです。
cell 内部にアクセスするには {} を使います。

以下のようにすれば、毎回サイズの分からない配列を追加収納できます。

pos{:,cN} = [pos{:,cN}; ax.CurrentPoint(1,1:2)];

中身にアクセスするには例えばこのようにします。

pos{cN}(:,1)  % x data

複数画（かく）対応ライブスクリプト

clear
close all

f = figure(1);
f.ToolBar = 'none';
f.Pointer = 'hand';  % change cursor shape
f.Color = [0.75 0.75 0.75];  % background color
ax = axes;  % create axes
axis equal  % square aspect ratio

DCTn = 32;  % DCT length

minSt = 100;
for k=1:5  % do 5 times drawing
    nSt = 0;  % number of strokes
    while true   
        nSt = nSt + 1;
        pos{:,nSt} = ax.CurrentPoint(1,1:2);  % get cursor position
        p(k,nSt) = plot(ax, pos{nSt}(:,1),pos{nSt}(:,2), 'b', Linewidth=2);  % make strok object
        xlim(ax, [0 600]);  ylim(ax, [0 600]);
        set(ax,'xtick',[],'ytick',[])
        box(ax,'on')
        hold(ax,'on')

        % set each callback
        f.WindowButtonDownFcn = 'pos{:,nSt} = ax.CurrentPoint(1,1:2);';  % get 1st position
        f.WindowButtonUpFcn = 'start = false;';  % end of one strok
        f.WindowButtonMotionFcn = 'movFlag = true;';  % draw strok

        w = waitforbuttonpress;  % wait for the each 1st click
        if w; break; end  % key pressed -> next try

        start = true;  movFlag = false;
        while start
            if movFlag
                pos{:,nSt} = [pos{:,nSt}; ax.CurrentPoint(1,1:2)];  % add current position to strok
                p(k,nSt).XData = pos{nSt}(:,1);  % update the strok x data
                p(k,nSt).YData = pos{nSt}(:,2);  % update the strok y data
                movFlag = false;
            end
            drawnow  % interrupt for event
        end

        % clear callbacks
        f.WindowButtonDownFcn = '';
        f.WindowButtonMotionFcn = '';
        f.WindowButtonUpFcn = '';
    end

    % delete too short strok
    for st = 1:nSt
        n = size(pos{st},1)
        if (n < 4)
            pos(st) = [];  % delete too short strok
        end
    end
    nSt = size(pos,2);
    minSt = min(minSt, nSt);  % set to minimum strok

    % calculate mean stroks
    if k == 1
        fxm = zeros(minSt,DCTn);  fym = zeros(minSt,DCTn);
    end
    for st = 1:minSt
        n = size(pos{st},1)
            t0 = (1:n)';  % index of orginal strok data
            t = (0:n-1)'/(n-1);  % for interpolation input
            tt = 0:1/DCTn:1-1/DCTn;  % for interpolation output
            xx = spline(t,pos{st}(t0,1),tt);  % transform data to suitable for DCT
            yy = spline(t,pos{st}(t0,2),tt);  % transform data to suitable for DCT
            fx = dct(xx);  fy = dct(yy);  % get DCT coefficients

            % update running mean
            fxm(st,:) = (fxm(st,:) * (k-1) + fx) / k;
            fym(st,:) = (fym(st,:) * (k-1) + fy) / k;

            for pk = k-1:-1:1
                p(pk,st).Color = [p(pk,st).Color pk/k/4];  % decay transparancy of old stroks
                pm(pk,st).Color = [pm(pk,st).Color 0];  % erase old mean stroks
            end
            p(k,st).Color = [p(k,st).Color k/k/4];  % decay transparancy of old strok
            pm(k,st) = plot(ax, idct(fxm(st,:)), idct(fym(st,:)), 'r', Linewidth=2);  % draw mean strok
    end
end

hold(ax,'off')
for st = 1:nSt
    p(k,st).Color = [p(k,st).Color k/k/4];  % decay transparancy of the last strok
end

f.Pointer = 'arrow';  % restore cursor shape
disp('finished.')

これも、ライブスクリプトを新規作成してそこにペーストしてください。

クリック＆ドラッグして一画を描き、クリックを離すと二画目に移行します。あまり短いとエラーが出るかもしれません。

一文字描き終わったら、任意のキーを押してください。二回目に移行します。スペースキーでもカーソルキーでも大丈夫です。

一画バージョンと同じく、5回で終わるようになっています。
別の操作で終わらすようにしてもよいですね。

動作例

動作例（任意画数バージョン）

タッチスクリーンPC であれば、そのままペンや指が使えます。

動作例（タッチスクリーンPCでの例）

ただし、書き順や書く方向が同じでないとダメです！

任意書き順・任意書き方向対応版

ということで、一応、任意書き順・任意書き方向対応版も作ってみました。

任意書き順・任意書き方向対応ライブスクリプト

clear
close all

f = figure(1);
f.ToolBar = 'none';
f.Pointer = 'hand';  % change cursor shape
f.Color = [0.75 0.75 0.75];  % background color
ax = axes;  % create axes
axis equal  % square aspect ratio

DCTn = 32;  % DCT length

minSt = 100;
for k=1:5  % do 5 times drawing
    nSt = 0;  % number of strokes
    while true   
        nSt = nSt + 1;
        pos{:,nSt} = ax.CurrentPoint(1,1:2);  % get cursor position
        p(k,nSt) = plot(ax, pos{nSt}(:,1),pos{nSt}(:,2), 'b', Linewidth=2);  % make strok object
        xlim(ax, [0 600]);  ylim(ax, [0 600]);
        % axis(ax, 'off')  % erase axes
        set(ax,'xtick',[],'ytick',[])
        box(ax,'on')
        hold(ax,'on')

        % set each callback
        f.WindowButtonDownFcn = 'pos{:,nSt} = ax.CurrentPoint(1,1:2);';  % get 1st position
        f.WindowButtonUpFcn = 'start = false;';  % end of one strok
        f.WindowButtonMotionFcn = 'movFlag = true;';  % draw strok

        w = waitforbuttonpress;  % wait for the each 1st click
        if w; break; end  % key pressed -> next try

        start = true;  movFlag = false;
        while start
            if movFlag
                pos{:,nSt} = [pos{:,nSt}; ax.CurrentPoint(1,1:2)];  % add current position to strok
                p(k,nSt).XData = pos{nSt}(:,1);  % update the strok x data
                p(k,nSt).YData = pos{nSt}(:,2);  % update the strok y data
                movFlag = false;
            end
            drawnow  % interrupt for event
        end

        % clear callbacks
        f.WindowButtonDownFcn = '';
        f.WindowButtonMotionFcn = '';
        f.WindowButtonUpFcn = '';
    end

    % delete too short strok
    for st = 1:size(pos,2) %nSt
        n = size(pos{st},1);
        if (n < 4)
            pos(st) = [];  % delete too short strok
        end
    end
    nSt = size(pos,2);
    minSt = min(minSt, nSt);  % set to minimum strok

    % calculate mean stroks
    if k == 1
        fxm = zeros(minSt,DCTn);  fym = zeros(minSt,DCTn);
    end
    for st = 1:minSt
        n = size(pos{st},1);
            t0 = (1:n)';  % index of orginal strok data
            t = (0:n-1)'/(n-1);  % for interpolation input
            tt = 0:1/DCTn:1-1/DCTn;  % for interpolation output
            xx = spline(t,pos{st}(t0,1),tt);  % transform data to suitable for DCT
            yy = spline(t,pos{st}(t0,2),tt);  % transform data to suitable for DCT
            fx = dct(xx);  fy = dct(yy);  % get DCT coefficients
            fxi = dct(flip(xx));  fyi = dct(flip(yy));  % filped strok

            % find nearest strok
            if k==1
                idx = st;
            else
                fxmidct = idct(fxm,DCTn,2);  fymidct = idct(fym,DCTn,2);
                dist = mean(sqrt((fxmidct - xx ) .^ 2 + (fymidct - yy ) .^ 2), 2);
                dsiv = mean(sqrt((fxmidct - flip(xx)) .^ 2 + (fymidct - flip(yy)) .^ 2), 2);
                dss =[dist; dsiv];
                [~,idx] = min(dss,[],'all');

                if idx > minSt  % flip strok
                    fx = fxi;
                    fy = fyi;
                    idx = idx - minSt
                end
            end

            % update running mean
            fxm(idx,:) = (fxm(idx,:) * (k-1) + fx) / k;
            fym(idx,:) = (fym(idx,:) * (k-1) + fy) / k;

            for pk = k:-1:1
                p(pk,st).Color = [p(pk,st).Color pk/k/4];  % decay transparancy of old stroks
            end
            if k ~= 1
                pm(st).Color = [pm(pk,st).Color 0];  % erase old mean stroks
            end
            pm(st) = plot(ax, idct(fxm(idx,:)), idct(fym(idx,:)), 'r', Linewidth=2);  % draw mean strok
    end
end

hold(ax,'off')
for st = 1:nSt
    p(k,st).Color = [p(k,st).Color k/k/4];  % decay transparancy of the last strok
end

f.Pointer = 'arrow';  % restore cursor shape
disp('finished.')

dct 前の今書いたストロークと、idct してストロークデータに戻した平均ストロークのノルム（ユークリッド距離）を計算して、一番近いストロークと平均を取るようにしています。このとき、flipで書き方向が逆のデータも作っておき、それも含めて距離計算を行っています。位置が近いと間違える可能性はあります。

たまにエラーが出るかもしれませんが、バグを見つけたら直して教えてください。（￣ー￣；

動作例

動作例（任意書き順・任意書き方向対応版）

まとめ

確かに、「平均化すると綺麗になる」感じはありますね。

何か他にも面白い使い道がありそうななさそうな・・。（￢＿￢）

何か思い付いたら教えてください。(^ｰ^)ﾉ

タイトル画像モデル：海老澤一恵

この記事が気に入ったらサポートをしてみませんか？