public class RBF {
double K[][] = new double[100][100];
double ax[] = new double[100];
double ay[] = new double[100];
double tx[] = new double[100];
double ty[] = new double[100];
int np;
double r0, r2;
public RBF() {
r0 = 10000;
r2 = r0 * r0;
}
// in this function, we solve K*a = pt_out
public void createFransform(double[] x_in, double[] y_in, int[] x_out,
int[] y_out, int np) {
// initialize K
double d, e, f;
for (int i = 0; i < np; i++) {
K[i][i] = 1.0f;
for (int j = 0; j < i; j++) {
d = x_in[i] - x_in[j];
e = y_in[i] - y_in[j];
f = d * d + e * e;
if (f < r2) {
f = Math.sqrt(f);
e = f / r0;
d = 1 - e;
K[i][j] = d * d * d * d * (1 + 4 * e);
} else
K[i][j] = 0;
}
}
// cholesky decomposition (in place)
// after which, we have an upper triangular matrix K
// satisfying K(old) = K^T * K
for (int j = 0; j < np; j++) {
d = Math.sqrt(K[j][j]);
K[j][j] = d;
for (int i = j + 1; i < np; i++) {
e = K[i][j] / d;
K[j][i] = e;
for (int k = j + 1; k <= i; k++)
K[i][k] -= e * K[j][k];
}
}
// solve K^T * t = pt_out;
for (int i = 0; i < np; i++) {
d = 0;
e = 0;
for (int j = 0; j < i; j++) {
d += K[j][i] * tx[j];
e += K[j][i] * ty[j];
}
tx[i] = (x_out[i] - d) / K[i][i];
ty[i] = (y_out[i] - e) / K[i][i];
}
// then, solve K * a = t
for (int i = np - 1; i >= 0; i--) {
d = 0;
e = 0;
for (int j = np - 1; j > i; j--) {
d += K[i][j] * ax[j];
e += K[i][j] * ay[j];
}
ax[i] = (tx[i] - d) / K[i][i];
ay[i] = (ty[i] - e) / K[i][i];
tx[i] = x_in[i];
ty[i] = y_in[i];
}
this.np = np;
}
public void transformPoints(double[] x_in, double[] y_in, int n,
int[] x_out, int[] y_out) {
double d, e, f = 0, g, h;
for (int i = 0; i < n; i++) {
g = 0;
h = 0;
for (int j = 0; j < np; j++) {
d = x_in[i] - tx[j];
e = y_in[i] - ty[j];
f = d * d + e * e;
if (f < r2) {
f = Math.sqrt(f);
e = f / r0;
d = 1 - e;
d = d * d * d * d * (1 + 4 * e);
g += ax[j] * d;
h += ay[j] * d;
}
}
x_out[i] = (int) g;
y_out[i] = (int) h;
}
}
int x_out = 0;
int y_out = 0;
public void transformPoint(double x_in, double y_in) {
double d = 0, e = 0, f = 0, g = 0, h = 0;
for (int j = 0; j < np; j++) {
d = x_in - tx[j];
e = y_in - ty[j];
f = d * d + e * e;
if (f < r2) {
f = Math.sqrt(f);
e = f / r0;
d = 1 - e;
d = d * d * d * d * (1 + 4 * e);
g += ax[j] * d;
h += ay[j] * d;
}
}
x_out = (int) g;
y_out = (int) h;
}
}
