Radial Basis Functionによる画像変形

# はじめに
本記事ではカーネル回帰による画像の変形について書きたいと思います。カーネル回帰を行うことで点群データを滑らかに通るような関数を求めることが出来ます（カーネル回帰そのものの説明については省略します）。本記事ではカーネルとしてRadial Basis Function（RBF）を用いて、画像や3Dモデルの変形に適用することを考えます(Radial Basis Function Transformation, RBFT) 。

![rbf.png](https://qiita-image-store.s3.amazonaws.com/0/176122/1a7dceb2-af9b-6ef0-7bd6-1d91c9caa630.png)


![example.png](https://qiita-image-store.s3.amazonaws.com/0/176122/38903b82-fccf-a3ef-4392-522d736b3cc1.png)

# RBFの一般式
![control.png](https://qiita-image-store.s3.amazonaws.com/0/176122/f97a4a3e-74cc-d4a0-20f4-ef276dbc1831.png)

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

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

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

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

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

```math
\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$としたときにアフィン変換に一致しない

そのため、変形に適用するためには以下のように式を変える必要があります。RBFの項とバイアスの項と線形項を組み合わせたカーネル関数とみなせますね。

```math
\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}
```
新しく加えた項はアフィン変換を表します。したがって、先ほどの行列を書き換えると

```math
\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による重みが小さくなっていき、アフィン変換に近づくことがわかります。

```math
\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クラス
ヘッダーファイル

```c++

#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);
};
```
ソースファイル

```c++

#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;
}

```

## 使用方法
こんな感じで使えばオッケーです

```c++

	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](https://qiita-image-store.s3.amazonaws.com/0/176122/bd5b0a3a-309b-0a64-4d4f-ab0f1cad8e6e.png)


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

![lammda.png](https://qiita-image-store.s3.amazonaws.com/0/176122/5c849f0d-0b56-a582-fc52-0a8c2349d772.png)

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