-
Notifications
You must be signed in to change notification settings - Fork 0
/
rbf.js
129 lines (102 loc) · 3.02 KB
/
rbf.js
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
var distanceFunctions = {
'linear': distanceLinear,
'cubic': distanceCubic,
'quintic': distanceQuintic,
'thin-plate': distanceThinPlate,
'gaussian': distanceGaussian,
'inverse-multiquadric': distanceInverseMultiquadric,
'multiquadric': distanceMultiquadric
};
export default function RBF(points, values, distanceFunction, epsilon) {
var distance = distanceFunctions.linear;
if(distanceFunction) {
distance = typeof distanceFunction !== 'string'
? distanceFunction
: distanceFunctions[distanceFunction];
}
// `identity` here serves as a shorthand to allocate
// the matrix nested array.
var M = numeric.identity(points.length);
// First compute the point to point distance matrix
// to allow computing epsilon if it's not provided
for(var j=0; j<points.length; j++) {
for(var i=0; i<points.length; i++) {
M[j][i] = norm(points[i], points[j]);
}
}
// if needed, compute espilon as the average distance between points
if(epsilon === undefined) {
epsilon = numeric.sum(M) / (Math.pow(points.length, 2) - points.length);
}
// update the matrix to reflect the requested distance function
for(var j=0; j<points.length; j++) {
for(var i=0; i<points.length; i++) {
M[j][i] = distance(M[j][i], epsilon);
}
}
// determine dimensionality of target values
var sample = values[0];
var D = typeof sample === 'number'
? 1
: sample.length;
// generalize to vector values
if(D === 1) {
values = values.map(function(value) {
return [value];
});
}
// reshape values by component
var tmp = new Array(D);
for(var i=0; i<D; i++) {
tmp[i] = values.map(function(value) {
return value[i];
});
}
values = tmp;
// Compute basis functions weights by solving
// the linear system of equations for each target component
var w = new Array(D);
var LU = numeric.LU(M);
for(var i=0; i<D; i++) {
w[i] = numeric.LUsolve(LU, values[i]);
}
// The returned interpolant will compute the value at any point
// by summing the weighted contributions of the input points
function interpolant(p) {
var distances = new Array(points.length);
for(var i=0; i<points.length; i++) {
distances[i] = distance(norm(p, points[i]), epsilon);
}
var sums = new Array(D);
for(var i=0; i<D; i++) {
sums[i] = numeric.sum(numeric.mul(distances, w[i]));
}
return sums;
}
return interpolant;
}
function norm(pa, pb) {
return numeric.norm2(numeric.sub(pb, pa));
}
function distanceLinear(r) {
return r;
}
function distanceCubic(r) {
return Math.pow(r, 3);
}
function distanceQuintic(r) {
return Math.pow(r, 5);
}
function distanceThinPlate(r) {
if(r === 0) return 0;
return Math.pow(r, 2) * Math.log(r);
}
function distanceGaussian(r, epsilon) {
return Math.exp(- Math.pow(r / epsilon, 2));
}
function distanceInverseMultiquadric(r, epsilon) {
return 1.0 / Math.sqrt(Math.pow(r / epsilon, 2) + 1);
}
function distanceMultiquadric(r, epsilon) {
return Math.sqrt(Math.pow(r / epsilon, 2) + 1);
}