C++
OpenGL
RBF

Radial Basis Functionによる画像変形

はじめに

Radial Basis Function (RBF:放射基底関数) を使った画像の変形について書きたいと思います。
RBFは非常に便利な補間関数で、主に点群データの滑らかな補間やネットワークのカーネル関数として利用されます。本記事では、RBFを画像や3Dモデルの変形 Radial Basis Function Transformation (RBFT) として適用することを考えます。

rbf.png

example.png

RBFの一般式

control.png

$N$個の位置ベクトル$\boldsymbol{A_i}\in\boldsymbol{R^m}$ (これを制御点と呼びます)が変形して$\boldsymbol{f(A_i)}\in\boldsymbol{R^m}$に移動したとします。このとき任意の位置$\boldsymbol{x}$の変形後の位置は、RBFを用いて以下のように書けます。

\boldsymbol{f(x)} = \sum_{i=1}^N \boldsymbol{c_i} \varphi(||\boldsymbol{x-A_i}||)

ここで、$\varphi(||x||)$は距離のみに依存する実数値関数であり、一般的にはGaussianなどが用いられます。

\varphi(x) = \exp(-x^2/\sigma^2)

各制御点$\boldsymbol{A_i}$による重み$\boldsymbol{c_i}$が未知数なので$\boldsymbol{x=A_i}$と代入して求める必要があります。
したがって、行列の形で書き下すと

\left(
\begin{array}{c}
\boldsymbol{f(A_1)} \\
\vdots \\
\boldsymbol{f(A_i)} \\
\vdots \\
\boldsymbol{f(A_N)}
\end{array}
\right)
=
\left(
\begin{array}{ccccc}
\varphi(||\boldsymbol{A_1-A_1}||) & \cdots & \varphi(||\boldsymbol{A_1-A_i}||) & \cdots & \varphi(||\boldsymbol{A_1-A_N}||)\\
\vdots & \ddots & & & \vdots \\
\varphi(||\boldsymbol{A_i-A_1}||) & & \varphi(||\boldsymbol{A_i-A_i}||) & & \varphi(||\boldsymbol{A_i-A_N}||) \\
\vdots & & & \ddots & \vdots \\
\varphi(||\boldsymbol{A_N-A_1}||) & \cdots & \varphi(||\boldsymbol{A_N-A_i}||) & \cdots & \varphi(||\boldsymbol{A_N-A_N}||)
\end{array}
\right)

\left(
\begin{array}{c}
\boldsymbol{c_1} \\
\vdots \\
\boldsymbol{c_i} \\
\vdots \\
\boldsymbol{c_N} \\
\end{array}
\right)

RBFTの一般式

RBFの式をそのまま適用すると2つの問題があります。

  • $\varphi(||x||)$は距離のみに依存して決まる関数なので、制御点から等距離に存在する点が全て同じ点に収束してしまう
  • $N=3$としたときにアフィン変換に一致しない

そのため、変形に適用する際には以下のように式を変える必要があります。

\boldsymbol{f(x)}_i = \sum_{i=1}^N \boldsymbol{c_i} \varphi(||\boldsymbol{x-A_i}||) + \boldsymbol{\alpha_0} + \sum_{j=1}^m \boldsymbol{x_{ij}\alpha_j}

新しく加えた項はアフィン変換を表します。したがって、先ほどの行列を書き換えると

\left(
\begin{array}{c}
\boldsymbol{f(A_1)} \\
\vdots \\
\boldsymbol{f(A_N)} \\
0 \\
0 \\
\vdots \\
0
\end{array}
\right)
=
\left(
\begin{array}{ccccccc}
\varphi(||\boldsymbol{A_1-A_1}||) & \cdots & \varphi(||\boldsymbol{A_1-A_N}||) & 1 & A_{11} & \cdots & A_{1m} \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
\varphi(||\boldsymbol{A_N-A_1}||) & \cdots & \varphi(||\boldsymbol{A_N-A_N}||) & 1 & A_{N1} & \cdots & A_{Nm} \\
1 & \cdots & 1 & 0 & 0 & \cdots & 0 \\
A_{11} & \cdots & A_{m1} & 0 & 0 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots& \vdots & \ddots & \vdots \\
A_{1m} & \cdots & A_{Nm} & 0 & 0 & \cdots & 0
\end{array}
\right)

\left(
\begin{array}{c}
\boldsymbol{c_1} \\
\vdots \\
\boldsymbol{c_N} \\
\boldsymbol{\alpha_0} \\
\boldsymbol{\alpha_1} \\
\vdots \\
\boldsymbol{\alpha_m}
\end{array}
\right)

さらに、ここでアフィン変換の程度の強さを表す変数$\lambda$を導入し、$(0,0)$ から$(N,N)$ までの部分行列について対角成分にプラスします。$\lambda$がおおきくなることで相対的にRBFによる重みが小さくなっていき、アフィン変換に近づきます。

\left(
\begin{array}{c}
\boldsymbol{f(A_1)} \\
\vdots \\
\boldsymbol{f(A_N)} \\
0 \\
0 \\
\vdots \\
0
\end{array}
\right)
=
\left(
\begin{array}{ccccccc}
\varphi(||\boldsymbol{A_1-A_1}||)+ \lambda & \cdots & \varphi(||\boldsymbol{A_1-A_N}||) & 1 & A_{11} & \cdots & A_{1m} \\
\vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
\varphi(||\boldsymbol{A_N-A_1}||) & \cdots & \varphi(||\boldsymbol{A_N-A_N}||) + \lambda & 1 & A_{N1} & \cdots & A_{Nm} \\
1 & \cdots & 1 & 0 & 0 & \cdots & 0 \\
A_{11} & \cdots & A_{m1} & 0 & 0 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots& \vdots & \ddots & \vdots \\
A_{1m} & \cdots & A_{Nm} & 0 & 0 & \cdots & 0
\end{array}
\right)

\left(
\begin{array}{c}
\boldsymbol{c_1} \\
\vdots \\
\boldsymbol{c_N} \\
\boldsymbol{\alpha_0} \\
\boldsymbol{\alpha_1} \\
\vdots \\
\boldsymbol{\alpha_m}
\end{array}
\right)

実装

今回はC++でOpenCVとEigenを使用してクラスを作成しました。画像の読み込みや座標の扱いにOpenCVを、逆行列の導出にはEigenというライブラリを使っています。
ここでは、上記の行列を
$F=GX$
として$X$ を求めています。

RBFTクラス

ヘッダーファイル

#pragma once

using namespace std;
using namespace cv;

class RBFT {

private:
    int N; // 頂点数
    vector<Point2i> A, B; // 変形前・変形後の座標
    Eigen::MatrixXf F, G, X; // F = GX

    void SetF();
    void SetG();
    double potential(const Point2i Ai, const Point2i Aj);

public:
    RBFT(vector<Point2i> _before, vector<Point2i> _after);
    ~RBFT() {}

    void SetX();
    vector<Point2i> Exe(const vector<Point2i> _inputPoints);
};

ソースファイル

#include "rbft.h"

RBFT::RBFT(std::vector<cv::Point2i> _before, std::vector<cv::Point2i> _after)
{
    if (_before.size() != _after.size()) {
        std::cout << "[warning]: 対応点が一致していません" << std::endl;
    }
    else {
        N = _before.size();
        A = _before;
        B = _after;
    }
}

double RBFT::potential(const Point2i Ai, const Point2i Aj)
{
    double p;
    p = pow(pow(Ai.x - Aj.x, 2) + pow(Ai.y - Aj.y, 2), 0.5);
    //p = exp(-((pow(Ai.x - Aj.x, 2) + pow(Ai.y - Aj.y, 2)) / pow((Width + Height) / 12, 2)));
    return p;
}

void RBFT::SetF()
{
    F.resize(N + 3, 2);
    F.setZero();

    for (int i = 0; i < N; ++i) {
        F(i, 0) = B[i].x;
        F(i, 1) = B[i].y;
    }
}

void RBFT::SetG()
{
    G.resize(N + 3, N + 3);
    G.setZero();

    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            G(i, j) = potential(A[i], A[j]);
        }
        G(i, i) += lamda;
    }

    for (int i = 0; i < N; ++i) {
        G(i, N) = 1;
        G(i, N + 1) = A[i].x;
        G(i, N + 2) = A[i].y;
        G(N, i) = 1;
        G(N + 1, i) = A[i].x;
        G(N + 2, i) = A[i].y;
    }
}

void RBFT::SetX()
{
    /* 行列の宣言 */
    X.resize(N + 3, 2);

    // F,GのSet
    SetF();
    SetG();

    X = G.inverse() * F;
}

std::vector<cv::Point2i> RBFT::Exe(const vector<Point2i> inputPoints)
{
    std::vector<cv::Point2i> outputPoints;
    outputPoints.resize(inputPoints.size());

    for (int i = 0; i < inputPoints.size(); ++i) {

        float x, y;
        x = y = 0;

        for (int j = 0; j < N; ++j) {
            x += potential(inputPoints[i], A[j]) * X(j, 0);
            y += potential(inputPoints[i], A[j]) * X(j, 1);
        }

        x += (X(N, 0) + inputPoints[i].x * X(N + 1, 0) + inputPoints[i].y * X(N + 2, 0));
        y += (X(N, 1) + inputPoints[i].x * X(N + 1, 1) + inputPoints[i].y * X(N + 2, 1));

        outputPoints[i] = Point2i(x, y);
    }
    return outputPoints;
}

使用方法

こんな感じで使えばオッケーです

    RBFT *rbft = new RBFT(controlPointsBefore, controlPointsAfter);
    rbft->SetX();
    allPointsAfter = rbft->Exe(allPointsBefore);
    delete rbft;

実装の詳細はこちら (https://github.com/smygw72/RBFT-Radial-Basis-Function-Transformation-) を見てください。

結果

主にカーネル関数やアフィン変換の強さを表す$\lambda$を変更することで変形を制御することが出来ます。細かい調整は制御点の数を動的に増やすことで解決できます。

カーネル関数による制御

$\varphi(x) = \exp(-x^2/\sigma^2)$ (上段)
$\varphi(x) = |x|$ (下段)

kernel.png

アフィン変換による制御

左から順に、$\lambda=0,0.1,1,10$のときの結果です。$\lambda$ が大きくなるほどアフィン変換による影響が相対的に大きくなっていくのが分かります。

lammda.png

参考文献

"Image Warping by Radial Basis Functions: Application to Facial Expressions" Nur Arad, Nira Dyn, Daniel Reisfeld , Yehezkel Yeshurun