-
Notifications
You must be signed in to change notification settings - Fork 1
/
pchip.js
255 lines (218 loc) · 8.28 KB
/
pchip.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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
// pchip for JavaScript
/*
* Author: Michael Folkerts (http://bit.ly/folkerts)
* Date: July 2012
* Reference: http://www.mathworks.com/moler/interp.pdf
*
* Example Usage:
* var sampling_frequency; //per unit on x-axis
* var pchip_fit_points = pchip.fit(input_points,sampling_frequency); //default is monotone method
* or
* var pchip_fit_points = pchip.monotone(input_points,sampling_frequency); //same as above
* var pchip_fit_points = pchip.fit(input_points,"monotone"); //same as above
* or
* var pchip_fit_points = pchip.shape_preserving(input_points,sampling_frequency);
* var pchip_fit_points = pchip.fit(input_points,sampling_frequency,"shape_preserving");//same as last line
*
*/
var PCHIP = {
version: 0.1,
monotone : function(source,frequency){
/** just wrapping a function **/
return this.fit(source,frequency,"monotone");
},
shape_preserving : function(source,frequency){
/** just wrapping a function **/
return this.fit(source,frequency,"shape_preserving");
},
fit : function(source,frequency,method){
/** this function computes the pchip interpolation method with a monotone constraint or as a "shape preserving" fit
[http://www.mathworks.com/moler/interp.pdf]
source: is assumed to be an array of 2-touples (x,y)
frequency: is the number of samples per unit distance along x-axis
method (optional): defines how to choose slopes,
'monotone' is default
'shape_preserving' is visually appealing and avoids over-shooting data points
return value: a reference to a new 2-touple Array() with length determined by
frequency parameter and spacing of points in source
**/
//set default method:
method = typeof method !== 'undefined' ? method : "monotone";
var m;
switch(method){
case "monotone":
m = this.slopes_monotone(source);
break;
case "shape_preserving":
m = this.slopes_shape_preserving(source);
break;
default:
alert("slope method '"+method+"' is not defined in PCHIP!!");
}
dest = new Array(); //may be slow because we are pushing new elements
offset = 0;
for( var i=0; i<source.length-1; i++){
//fill dest array with sucessive pchip'd points
var thisPoint = source[i];
var nextPoint = source[i+1];
offset += this.chip_v(thisPoint, nextPoint, [ m[i], m[i+1] ], frequency, dest, offset); //faster?
//offset += this.chip_p(thisPoint[1],nextPoint[1],thisPoint[0],nextPoint[0],m[i],m[i+1],frequency,dest,offset);
}
return dest;
},
slopes_monotone: function(source){
var inLen = source.length;
var slope = new Array(inLen);// an array for slopes
var m = new Array(inLen); //an empty array for tangents (new slope parameter)
var a = new Array(inLen); //an empty array for factors
var b = new Array(inLen); //an empty array for factors
var t = new Array(inLen); //an empty array for factors
var to_fix = [];
//compute slopes:
var smallVal =1e-30; //avoid divide by zero
slope[0] = (source[1][1]-source[0][1])/(source[1][0]-source[0][0]+smallVal);
m[0] = slope[0];
if(slope[0] == 0){
to_fix.push(0);
}
for( var i=1; i<inLen-1; i++){
//var y0 = source[i][1];
//var y1 = source[i+1][1];
//var dx = source[i+1][0]-source[i][0]+smallVal;
//slope[i] = (y1-y0)/dx;
slope[i] = (source[i+1][1]-source[i][1])/(source[i+1][0]-source[i][0]+smallVal);
m[i] = (slope[i-1] + slope[i])/2.0;
if(slope[i] == 0.0){
to_fix.push(i);
}
}
m[inLen-1] = slope[inLen-2];
for( var i=0; i<to_fix.length;i++){
m[to_fix[i]] = 0.0;
m[to_fix[i]+1]=0.0;
}
to_fix = [];
for( var i=0; i<inLen-1; i++){
a[i] = m[i]/slope[i];
b[i] = m[i+1]/slope[i];
var dist = a[i]*a[i]+b[i]*b[i];
t[i] = 3.0/Math.sqrt(dist);
if( dist > 9.0 ){
to_fix.push(i);
}
}
for(var i=0; i<to_fix.length;i++){
var j = to_fix[i];
m[j] = t[j]*a[j]*slope[j];
m[j+1] = t[j]*b[j]*slope[j];
}
delete a,b,t,slope; //clean up
return m;
},
slopes_shape_preserving: function (pts){
m = new Array(pts.length); // array for new slopes
//locals
var w1 = 0.0;
var w2 = 0.0;
var h0 = 0.0; //h[k-1]
var h1 = 0.0; //h[k]
var slope0 = 0.0; //d[i-1]
var slope1 = 0.0; //d[i]
slope0 = this.grad(pts[0],pts[1]); //slope at first point
for(i=1; i<pts.length-1; i++){
slope1 = this.grad(pts[i],pts[i+1]); //forward looking slope at point [i]
//if slopes[i-1] or slopes[i] are zero OR opposite signs
if( slope0 == 0.0 || slope1 == 0.0 ){ //don't overshoot
m[i] = 0.0;
}else if( this.sgn(slope0) != this.sgn(slope1) ){ //local min/max
m[i] = 0.0;
}else{
h0 = pts[i][0] - pts[i-1][0];
h1 = pts[i+1][0] - pts[i][0];
w1 = 2.0*h1+h0;
w2 = h1 + 2.0*h0;
m[i] = (w1+w2)/( w1/slope0 + w2/slope1 );
}
slope0 = slope1;//shift slopes
}
//quick fix endpoints... need more research to get this totally right
m[0] = m[1];
m[pts.length-1] = m[pts.length-2];
return m;
},
grad : function (pt1,pt2){
/** the gradient (derivative) between two consecutive points (x1,y1) and (x2,y2)
* with x1<x2 **/
return (pt2[1]-pt1[1])/(pt2[0]-pt1[0]);
},
sgn : function (x){
/** the signum [http://en.wikipedia.org/wiki/Sign_function] function or sign of x **/
if( x < 0.0 ){
return -1;
}
if ( x > 0.0){
return +1;
}
return x; //when x == 0 or not a number
},
chip_p : function (yi,yii,xi,xii,mi,mii,freq,resultArray,offset){
/** This is the Cubic Hermite Interpolating Polynomial (chip) function taking *Point* values as arguments:
* yi, and yii are consecutive y values in a plot (points between which to interpolate
* xi, and xii are consecutive x valuse associated with the above y values
* mi, and mii are the slopes desired at each point
* freq is the number of sample points per unit along the x-axis
* resultArray is a pointer to an exisiting 2-touple (x,y) array to which the sample points are appended
* offset is the starting index at which points are appended
* note: a new 2-touple is pushed for each new point, offset is treated as current length of array
* numPts is returned to report the number of points added **/
var numPts = parseInt((xii-xi)*freq); //this could be zero if spacing is too small
var a = 2.0/Math.pow(xii-xi,3);
var b = (xii-xi)/2.0;
var c = 1.0/Math.pow(xii-xi,2);
for(var i=0; i<numPts; i++){
var X = xi+i/freq;
resultArray.push(Array(2));
resultArray[offset+i][0] = X;
// resultArray[offset+i][1] = yi *( 2.0/Math.pow(xii-xi,3)*Math.pow(X-xii,2)*(X-xi +(xii-xi)/2.0) )
// +yii*(-2.0/Math.pow(xii-xi,3)*Math.pow(X-xi ,2)*(X-xii-(xii-xi)/2.0) )
// + 1.0/Math.pow(xii-xi,2)*( mi*(X-xi)*Math.pow(X-xii,2)+mii*(X-xii)*Math.pow(X-xi,2) );
resultArray[offset+i][1] = yi *( a*Math.pow(X-xii,2)*(X-xi +b) )
+yii*(-a*Math.pow(X-xi ,2)*(X-xii-b) )
+ c *( mi*(X-xi)*Math.pow(X-xii,2)+mii*(X-xii)*Math.pow(X-xi,2) );
}
return numPts;
},
chip_v: function (point1, point2, deriv, freq, resultArray, offset){
/** This is the Cubic Hermite Interpolating Polynomial (chip) function taking *Vectors* as arguments:
* Point1 and Point2 are consecutive 2-touples (x,y) defining the points between which to interpolate
* deriv is a 2-touple (1,2) defining the slope at each point
* freq is the number of sample points per unit along the x-axis
* resultArray is a pointer to an exisiting 2-touple (x,y) array to which the sample points are appended
* offset is the starting index at which points are appended
* note: a new 2-touple is pushed for each new point, offset is treated as current length of array
* numPts is returned to report the number of points added **/
//make code easier to read:
var x1 = point1[0];
var x2 = point2[0];
var y1 = point1[1];
var y2 = point2[1];
var yp1= deriv[0];
var yp2= deriv[1];
///TODO: maybe just wrap the chip_v function
// var a=y1;
// var b=yp1;
var deltax=x2-x1;
var deltay=y2-y1;
var c=(deltay/deltax-yp1)/deltax;
var d=(yp1+yp2-2*deltay/deltax)/(deltax*deltax);
var numPts = parseInt(deltax*freq); //this could be zero if spacing is too small
for(var i=0; i<numPts; i++){
var x = x1+i/freq;
resultArray.push(Array(2));
resultArray[offset+i][0] = x;
// resultArray[offset+i][1] = a+(x-x1)*( b+(x-x1)*( c+d*(x-x2) ) );
resultArray[offset+i][1] = y1+(x-x1)*( yp1+(x-x1)*( c+d*(x-x2) ) );
}
return numPts;
}
}