484 lines
		
	
	
		
			13 KiB
		
	
	
	
		
			JavaScript
		
	
	
	
	
	
		
		
			
		
	
	
			484 lines
		
	
	
		
			13 KiB
		
	
	
	
		
			JavaScript
		
	
	
	
	
	
|  | // Extend the Array class
 | ||
|  | Array.prototype.max = function() { | ||
|  |     return Math.max.apply(null, this); | ||
|  | }; | ||
|  | Array.prototype.min = function() { | ||
|  |     return Math.min.apply(null, this); | ||
|  | }; | ||
|  | Array.prototype.mean = function() { | ||
|  |     var i, sum; | ||
|  |     for(i=0,sum=0;i<this.length;i++) | ||
|  | 	sum += this[i]; | ||
|  |     return sum / this.length; | ||
|  | }; | ||
|  | Array.prototype.pip = function(x, y) { | ||
|  |     var i, j, c = false; | ||
|  |     for(i=0,j=this.length-1;i<this.length;j=i++) { | ||
|  | 	if( ((this[i][1]>y) != (this[j][1]>y)) && | ||
|  | 	    (x<(this[j][0]-this[i][0]) * (y-this[i][1]) / (this[j][1]-this[i][1]) + this[i][0]) ) { | ||
|  | 	    c = !c; | ||
|  | 	} | ||
|  |     } | ||
|  |     return c; | ||
|  | } | ||
|  | 
 | ||
|  | var kriging = function() { | ||
|  |     var kriging = {}; | ||
|  | 
 | ||
|  |     var createArrayWithValues = function(value, n) { | ||
|  |         var array = []; | ||
|  |         for ( var i = 0; i < n; i++) { | ||
|  |             array.push(value); | ||
|  |         } | ||
|  |         return array; | ||
|  |     }, | ||
|  | 
 | ||
|  |     // Matrix algebra
 | ||
|  |     kriging_matrix_diag = function(c, n) { | ||
|  |         var Z = createArrayWithValues(0, n * n); | ||
|  |         for(i=0;i<n;i++) Z[i*n+i] = c; | ||
|  |         return Z; | ||
|  |     }, | ||
|  |     kriging_matrix_transpose = function(X, n, m) { | ||
|  | 	var i, j, Z = Array(m*n); | ||
|  | 	for(i=0;i<n;i++) | ||
|  | 	    for(j=0;j<m;j++) | ||
|  | 		Z[j*n+i] = X[i*m+j]; | ||
|  | 	return Z; | ||
|  |     }, | ||
|  |     kriging_matrix_scale = function(X, c, n, m) { | ||
|  | 	var i, j; | ||
|  | 	for(i=0;i<n;i++) | ||
|  | 	    for(j=0;j<m;j++) | ||
|  | 		X[i*m+j] *= c; | ||
|  |     }, | ||
|  |     kriging_matrix_add = function(X, Y, n, m) { | ||
|  | 	var i, j, Z = Array(n*m); | ||
|  | 	for(i=0;i<n;i++) | ||
|  | 	    for(j=0;j<m;j++) | ||
|  | 		Z[i*m+j] = X[i*m+j] + Y[i*m+j]; | ||
|  | 	return Z; | ||
|  |     }, | ||
|  |     // Naive matrix multiplication
 | ||
|  |     kriging_matrix_multiply = function(X, Y, n, m, p) { | ||
|  | 	var i, j, k, Z = Array(n*p); | ||
|  | 	for(i=0;i<n;i++) { | ||
|  | 	    for(j=0;j<p;j++) { | ||
|  | 		Z[i*p+j] = 0; | ||
|  | 		for(k=0;k<m;k++) | ||
|  | 		    Z[i*p+j] += X[i*m+k]*Y[k*p+j]; | ||
|  | 	    } | ||
|  | 	} | ||
|  | 	return Z; | ||
|  |     }, | ||
|  |     // Cholesky decomposition
 | ||
|  |     kriging_matrix_chol = function(X, n) { | ||
|  | 	var i, j, k, sum, p = Array(n); | ||
|  | 	for(i=0;i<n;i++) p[i] = X[i*n+i]; | ||
|  | 	for(i=0;i<n;i++) { | ||
|  | 	    for(j=0;j<i;j++) | ||
|  | 		p[i] -= X[i*n+j]*X[i*n+j]; | ||
|  | 	    if(p[i]<=0) return false; | ||
|  | 	    p[i] = Math.sqrt(p[i]); | ||
|  | 	    for(j=i+1;j<n;j++) { | ||
|  | 		for(k=0;k<i;k++) | ||
|  | 		    X[j*n+i] -= X[j*n+k]*X[i*n+k]; | ||
|  | 		X[j*n+i] /= p[i]; | ||
|  | 	    } | ||
|  | 	} | ||
|  | 	for(i=0;i<n;i++) X[i*n+i] = p[i]; | ||
|  | 	return true; | ||
|  |     }, | ||
|  |     // Inversion of cholesky decomposition
 | ||
|  |     kriging_matrix_chol2inv = function(X, n) { | ||
|  | 	var i, j, k, sum; | ||
|  | 	for(i=0;i<n;i++) { | ||
|  | 	    X[i*n+i] = 1/X[i*n+i]; | ||
|  | 	    for(j=i+1;j<n;j++) { | ||
|  | 		sum = 0; | ||
|  | 		for(k=i;k<j;k++) | ||
|  | 		    sum -= X[j*n+k]*X[k*n+i]; | ||
|  | 		X[j*n+i] = sum/X[j*n+j]; | ||
|  | 	    } | ||
|  | 	} | ||
|  | 	for(i=0;i<n;i++) | ||
|  | 	    for(j=i+1;j<n;j++) | ||
|  | 		X[i*n+j] = 0; | ||
|  | 	for(i=0;i<n;i++) { | ||
|  | 	    X[i*n+i] *= X[i*n+i]; | ||
|  | 	    for(k=i+1;k<n;k++) | ||
|  | 		X[i*n+i] += X[k*n+i]*X[k*n+i]; | ||
|  | 	    for(j=i+1;j<n;j++) | ||
|  | 		for(k=j;k<n;k++) | ||
|  | 		    X[i*n+j] += X[k*n+i]*X[k*n+j]; | ||
|  | 	} | ||
|  | 	for(i=0;i<n;i++) | ||
|  | 	    for(j=0;j<i;j++) | ||
|  | 		X[i*n+j] = X[j*n+i]; | ||
|  | 
 | ||
|  |     }, | ||
|  |     // Inversion via gauss-jordan elimination
 | ||
|  |     kriging_matrix_solve = function(X, n) { | ||
|  | 	var m = n; | ||
|  | 	var b = Array(n*n); | ||
|  | 	var indxc = Array(n); | ||
|  | 	var indxr = Array(n); | ||
|  | 	var ipiv = Array(n); | ||
|  | 	var i, icol, irow, j, k, l, ll; | ||
|  | 	var big, dum, pivinv, temp; | ||
|  | 
 | ||
|  | 	for(i=0;i<n;i++) | ||
|  | 	    for(j=0;j<n;j++) { | ||
|  | 		if(i==j) b[i*n+j] = 1; | ||
|  | 		else b[i*n+j] = 0; | ||
|  | 	    } | ||
|  | 	for(j=0;j<n;j++) ipiv[j] = 0; | ||
|  | 	for(i=0;i<n;i++) { | ||
|  | 	    big = 0; | ||
|  | 	    for(j=0;j<n;j++) { | ||
|  | 		if(ipiv[j]!=1) { | ||
|  | 		    for(k=0;k<n;k++) { | ||
|  | 			if(ipiv[k]==0) { | ||
|  | 			    if(Math.abs(X[j*n+k])>=big) { | ||
|  | 				big = Math.abs(X[j*n+k]); | ||
|  | 				irow = j; | ||
|  | 				icol = k; | ||
|  | 			    } | ||
|  | 			} | ||
|  | 		    } | ||
|  | 		} | ||
|  | 	    } | ||
|  | 	    ++(ipiv[icol]); | ||
|  | 
 | ||
|  | 	    if(irow!=icol) { | ||
|  | 		for(l=0;l<n;l++) { | ||
|  | 		    temp = X[irow*n+l]; | ||
|  | 		    X[irow*n+l] = X[icol*n+l]; | ||
|  | 		    X[icol*n+l] = temp; | ||
|  | 		} | ||
|  | 		for(l=0;l<m;l++) { | ||
|  | 		    temp = b[irow*n+l]; | ||
|  | 		    b[irow*n+l] = b[icol*n+l]; | ||
|  | 		    b[icol*n+l] = temp; | ||
|  | 		} | ||
|  | 	    } | ||
|  | 	    indxr[i] = irow; | ||
|  | 	    indxc[i] = icol; | ||
|  | 
 | ||
|  | 	    if(X[icol*n+icol]==0) return false; // Singular
 | ||
|  | 
 | ||
|  | 	    pivinv = 1 / X[icol*n+icol]; | ||
|  | 	    X[icol*n+icol] = 1; | ||
|  | 	    for(l=0;l<n;l++) X[icol*n+l] *= pivinv; | ||
|  | 	    for(l=0;l<m;l++) b[icol*n+l] *= pivinv; | ||
|  | 
 | ||
|  | 	    for(ll=0;ll<n;ll++) { | ||
|  | 		if(ll!=icol) { | ||
|  | 		    dum = X[ll*n+icol]; | ||
|  | 		    X[ll*n+icol] = 0; | ||
|  | 		    for(l=0;l<n;l++) X[ll*n+l] -= X[icol*n+l]*dum; | ||
|  | 		    for(l=0;l<m;l++) b[ll*n+l] -= b[icol*n+l]*dum; | ||
|  | 		} | ||
|  | 	    } | ||
|  | 	} | ||
|  | 	for(l=(n-1);l>=0;l--) | ||
|  | 	    if(indxr[l]!=indxc[l]) { | ||
|  | 		for(k=0;k<n;k++) { | ||
|  | 		    temp = X[k*n+indxr[l]]; | ||
|  | 		    X[k*n+indxr[l]] = X[k*n+indxc[l]]; | ||
|  | 		    X[k*n+indxc[l]] = temp; | ||
|  | 		} | ||
|  | 	    } | ||
|  | 
 | ||
|  | 	return true; | ||
|  |     }, | ||
|  | 
 | ||
|  |     // Variogram models
 | ||
|  |     kriging_variogram_gaussian = function(h, nugget, range, sill, A) { | ||
|  | 	return nugget + ((sill-nugget)/range)* | ||
|  | 	( 1.0 - Math.exp(-(1.0/A)*Math.pow(h/range, 2)) ); | ||
|  |     }, | ||
|  |     kriging_variogram_exponential = function(h, nugget, range, sill, A) { | ||
|  | 	return nugget + ((sill-nugget)/range)* | ||
|  | 	( 1.0 - Math.exp(-(1.0/A) * (h/range)) ); | ||
|  |     }, | ||
|  |     kriging_variogram_spherical = function(h, nugget, range, sill, A) { | ||
|  | 	if(h>range) return nugget + (sill-nugget)/range; | ||
|  | 	return nugget + ((sill-nugget)/range)* | ||
|  | 	( 1.5*(h/range) - 0.5*Math.pow(h/range, 3) ); | ||
|  |     }; | ||
|  | 
 | ||
|  |     // Train using gaussian processes with bayesian priors
 | ||
|  |     kriging.train = function(t, x, y, model, sigma2, alpha) { | ||
|  | 	var variogram = { | ||
|  | 	    t      : t, | ||
|  | 	    x      : x, | ||
|  | 	    y      : y, | ||
|  | 	    nugget : 0.0, | ||
|  | 	    range  : 0.0, | ||
|  | 	    sill   : 0.0, | ||
|  | 	    A      : 1/3, | ||
|  | 	    n      : 0 | ||
|  | 	}; | ||
|  | 	switch(model) { | ||
|  | 	case "gaussian": | ||
|  | 	    variogram.model = kriging_variogram_gaussian; | ||
|  | 	    break; | ||
|  | 	case "exponential": | ||
|  | 	    variogram.model = kriging_variogram_exponential; | ||
|  | 	    break; | ||
|  | 	case "spherical": | ||
|  | 	    variogram.model = kriging_variogram_spherical; | ||
|  | 	    break; | ||
|  | 	}; | ||
|  | 
 | ||
|  | 	// Lag distance/semivariance
 | ||
|  | 	var i, j, k, l, n = t.length; | ||
|  | 	var distance = Array((n*n-n)/2); | ||
|  | 	for(i=0,k=0;i<n;i++) | ||
|  | 	    for(j=0;j<i;j++,k++) { | ||
|  | 		distance[k] = Array(2); | ||
|  | 		distance[k][0] = Math.pow( | ||
|  | 		    Math.pow(x[i]-x[j], 2)+ | ||
|  | 		    Math.pow(y[i]-y[j], 2), 0.5); | ||
|  | 		distance[k][1] = Math.abs(t[i]-t[j]); | ||
|  | 	    } | ||
|  | 	distance.sort(function(a, b) { return a[0] - b[0]; }); | ||
|  | 	variogram.range = distance[(n*n-n)/2-1][0]; | ||
|  | 
 | ||
|  | 	// Bin lag distance
 | ||
|  | 	var lags = ((n*n-n)/2)>30?30:(n*n-n)/2; | ||
|  | 	var tolerance = variogram.range/lags; | ||
|  | 	var lag = createArrayWithValues(0,lags); | ||
|  | 	var semi = createArrayWithValues(0,lags); | ||
|  | 	if(lags<30) { | ||
|  | 	    for(l=0;l<lags;l++) { | ||
|  | 		lag[l] = distance[l][0]; | ||
|  | 		semi[l] = distance[l][1]; | ||
|  | 	    } | ||
|  | 	} | ||
|  | 	else { | ||
|  | 	    for(i=0,j=0,k=0,l=0;i<lags&&j<((n*n-n)/2);i++,k=0) { | ||
|  | 		while( distance[j][0]<=((i+1)*tolerance) ) { | ||
|  | 		    lag[l] += distance[j][0]; | ||
|  | 		    semi[l] += distance[j][1]; | ||
|  | 		    j++;k++; | ||
|  | 		    if(j>=((n*n-n)/2)) break; | ||
|  | 		} | ||
|  | 		if(k>0) { | ||
|  | 		    lag[l] /= k; | ||
|  | 		    semi[l] /= k; | ||
|  | 		    l++; | ||
|  | 		} | ||
|  | 	    } | ||
|  | 	    if(l<2) return variogram; // Error: Not enough points
 | ||
|  | 	} | ||
|  | 
 | ||
|  | 	// Feature transformation
 | ||
|  | 	n = l; | ||
|  | 	variogram.range = lag[n-1]-lag[0]; | ||
|  | 	 var X = createArrayWithValues(1,2 * n); | ||
|  | 	var Y = Array(n); | ||
|  | 	var A = variogram.A; | ||
|  | 	for(i=0;i<n;i++) { | ||
|  | 	    switch(model) { | ||
|  | 	    case "gaussian": | ||
|  | 		X[i*2+1] = 1.0-Math.exp(-(1.0/A)*Math.pow(lag[i]/variogram.range, 2)); | ||
|  | 		break; | ||
|  | 	    case "exponential": | ||
|  | 		X[i*2+1] = 1.0-Math.exp(-(1.0/A)*lag[i]/variogram.range); | ||
|  | 		break; | ||
|  | 	    case "spherical": | ||
|  | 		X[i*2+1] = 1.5*(lag[i]/variogram.range)- | ||
|  | 		    0.5*Math.pow(lag[i]/variogram.range, 3); | ||
|  | 		break; | ||
|  | 	    }; | ||
|  | 	    Y[i] = semi[i]; | ||
|  | 	} | ||
|  | 
 | ||
|  | 	// Least squares
 | ||
|  | 	var Xt = kriging_matrix_transpose(X, n, 2); | ||
|  | 	var Z = kriging_matrix_multiply(Xt, X, 2, n, 2); | ||
|  | 	Z = kriging_matrix_add(Z, kriging_matrix_diag(1/alpha, 2), 2, 2); | ||
|  | 	var cloneZ = Z.slice(0); | ||
|  | 	if(kriging_matrix_chol(Z, 2)) | ||
|  | 	    kriging_matrix_chol2inv(Z, 2); | ||
|  | 	else { | ||
|  | 	    kriging_matrix_solve(cloneZ, 2); | ||
|  | 	    Z = cloneZ; | ||
|  | 	} | ||
|  | 	var W = kriging_matrix_multiply(kriging_matrix_multiply(Z, Xt, 2, 2, n), Y, 2, n, 1); | ||
|  | 
 | ||
|  | 	// Variogram parameters
 | ||
|  | 	variogram.nugget = W[0]; | ||
|  | 	variogram.sill = W[1]*variogram.range+variogram.nugget; | ||
|  | 	variogram.n = x.length; | ||
|  | 
 | ||
|  | 	// Gram matrix with prior
 | ||
|  | 	n = x.length; | ||
|  | 	var K = Array(n*n); | ||
|  | 	for(i=0;i<n;i++) { | ||
|  | 	    for(j=0;j<i;j++) { | ||
|  | 		K[i*n+j] = variogram.model(Math.pow(Math.pow(x[i]-x[j], 2)+ | ||
|  | 						    Math.pow(y[i]-y[j], 2), 0.5), | ||
|  | 					   variogram.nugget, | ||
|  | 					   variogram.range, | ||
|  | 					   variogram.sill, | ||
|  | 					   variogram.A); | ||
|  | 		K[j*n+i] = K[i*n+j]; | ||
|  | 	    } | ||
|  | 	    K[i*n+i] = variogram.model(0, variogram.nugget, | ||
|  | 				       variogram.range, | ||
|  | 				       variogram.sill, | ||
|  | 				       variogram.A); | ||
|  | 	} | ||
|  | 
 | ||
|  | 	// Inverse penalized Gram matrix projected to target vector
 | ||
|  | 	var C = kriging_matrix_add(K, kriging_matrix_diag(sigma2, n), n, n); | ||
|  | 	var cloneC = C.slice(0); | ||
|  | 	if(kriging_matrix_chol(C, n)) | ||
|  | 	    kriging_matrix_chol2inv(C, n); | ||
|  | 	else { | ||
|  | 	    kriging_matrix_solve(cloneC, n); | ||
|  | 	    C = cloneC; | ||
|  | 	} | ||
|  | 
 | ||
|  | 	// Copy unprojected inverted matrix as K
 | ||
|  | 	var K = C.slice(0); | ||
|  | 	var M = kriging_matrix_multiply(C, t, n, n, 1); | ||
|  | 	variogram.K = K; | ||
|  | 	variogram.M = M; | ||
|  | 
 | ||
|  | 	return variogram; | ||
|  |     }; | ||
|  | 
 | ||
|  |     // Model prediction
 | ||
|  |     kriging.predict = function(x, y, variogram) { | ||
|  | 	var i, k = Array(variogram.n); | ||
|  | 	for(i=0;i<variogram.n;i++) | ||
|  | 	    k[i] = variogram.model(Math.pow(Math.pow(x-variogram.x[i], 2)+ | ||
|  | 					    Math.pow(y-variogram.y[i], 2), 0.5), | ||
|  | 				   variogram.nugget, variogram.range, | ||
|  | 				   variogram.sill, variogram.A); | ||
|  | 	return kriging_matrix_multiply(k, variogram.M, 1, variogram.n, 1)[0]; | ||
|  |     }; | ||
|  |     kriging.variance = function(x, y, variogram) { | ||
|  | 	var i, k = Array(variogram.n); | ||
|  | 	for(i=0;i<variogram.n;i++) | ||
|  | 	    k[i] = variogram.model(Math.pow(Math.pow(x-variogram.x[i], 2)+ | ||
|  | 					    Math.pow(y-variogram.y[i], 2), 0.5), | ||
|  | 				   variogram.nugget, variogram.range, | ||
|  | 				   variogram.sill, variogram.A); | ||
|  | 	return variogram.model(0, variogram.nugget, variogram.range, | ||
|  | 			variogram.sill, variogram.A)+ | ||
|  | 	kriging_matrix_multiply(kriging_matrix_multiply(k, variogram.K, | ||
|  | 							1, variogram.n, variogram.n), | ||
|  | 				k, 1, variogram.n, 1)[0]; | ||
|  |     }; | ||
|  | 
 | ||
|  |     // Gridded matrices or contour paths
 | ||
|  |     kriging.grid = function(polygons, variogram, width) { | ||
|  | 	var i, j, k, n = polygons.length; | ||
|  | 	if(n==0) return; | ||
|  | 
 | ||
|  | 	// Boundaries of polygons space
 | ||
|  | 	var xlim = [polygons[0][0][0], polygons[0][0][0]]; | ||
|  | 	var ylim = [polygons[0][0][1], polygons[0][0][1]]; | ||
|  | 	for(i=0;i<n;i++) // Polygons
 | ||
|  | 	    for(j=0;j<polygons[i].length;j++) { // Vertices
 | ||
|  | 		if(polygons[i][j][0]<xlim[0]) | ||
|  | 		    xlim[0] = polygons[i][j][0]; | ||
|  | 		if(polygons[i][j][0]>xlim[1]) | ||
|  | 		    xlim[1] = polygons[i][j][0]; | ||
|  | 		if(polygons[i][j][1]<ylim[0]) | ||
|  | 		    ylim[0] = polygons[i][j][1]; | ||
|  | 		if(polygons[i][j][1]>ylim[1]) | ||
|  | 		    ylim[1] = polygons[i][j][1]; | ||
|  | 	    } | ||
|  | 
 | ||
|  | 	// Alloc for O(n^2) space
 | ||
|  | 	var xtarget, ytarget; | ||
|  | 	var a = Array(2), b = Array(2); | ||
|  | 	var lxlim = Array(2); // Local dimensions
 | ||
|  | 	var lylim = Array(2); // Local dimensions
 | ||
|  | 	var x = Math.ceil((xlim[1]-xlim[0])/width); | ||
|  | 	var y = Math.ceil((ylim[1]-ylim[0])/width); | ||
|  | 
 | ||
|  | 	var A = Array(x+1); | ||
|  | 	for(i=0;i<=x;i++) A[i] = Array(y+1); | ||
|  | 	for(i=0;i<n;i++) { | ||
|  | 	    // Range for polygons[i]
 | ||
|  | 	    lxlim[0] = polygons[i][0][0]; | ||
|  | 	    lxlim[1] = lxlim[0]; | ||
|  | 	    lylim[0] = polygons[i][0][1]; | ||
|  | 	    lylim[1] = lylim[0]; | ||
|  | 	    for(j=1;j<polygons[i].length;j++) { // Vertices
 | ||
|  | 		if(polygons[i][j][0]<lxlim[0]) | ||
|  | 		    lxlim[0] = polygons[i][j][0]; | ||
|  | 		if(polygons[i][j][0]>lxlim[1]) | ||
|  | 		    lxlim[1] = polygons[i][j][0]; | ||
|  | 		if(polygons[i][j][1]<lylim[0]) | ||
|  | 		    lylim[0] = polygons[i][j][1]; | ||
|  | 		if(polygons[i][j][1]>lylim[1]) | ||
|  | 		    lylim[1] = polygons[i][j][1]; | ||
|  | 	    } | ||
|  | 
 | ||
|  | 	    // Loop through polygon subspace
 | ||
|  | 	    a[0] = Math.floor(((lxlim[0]-((lxlim[0]-xlim[0])%width)) - xlim[0])/width); | ||
|  | 	    a[1] = Math.ceil(((lxlim[1]-((lxlim[1]-xlim[1])%width)) - xlim[0])/width); | ||
|  | 	    b[0] = Math.floor(((lylim[0]-((lylim[0]-ylim[0])%width)) - ylim[0])/width); | ||
|  | 	    b[1] = Math.ceil(((lylim[1]-((lylim[1]-ylim[1])%width)) - ylim[0])/width); | ||
|  | 	    for(j=a[0];j<=a[1];j++) | ||
|  | 		for(k=b[0];k<=b[1];k++) { | ||
|  | 		    xtarget = xlim[0] + j*width; | ||
|  | 		    ytarget = ylim[0] + k*width; | ||
|  | 		    if(polygons[i].pip(xtarget, ytarget)) | ||
|  | 			A[j][k] = kriging.predict(xtarget, | ||
|  | 						  ytarget, | ||
|  | 						  variogram); | ||
|  | 		} | ||
|  | 	} | ||
|  | 	A.xlim = xlim; | ||
|  | 	A.ylim = ylim; | ||
|  | 	A.zlim = [variogram.t.min(), variogram.t.max()]; | ||
|  | 	A.width = width; | ||
|  | 	return A; | ||
|  |     }; | ||
|  |     kriging.contour = function(value, polygons, variogram) { | ||
|  | 
 | ||
|  |     }; | ||
|  | 
 | ||
|  |     // Plotting on the DOM
 | ||
|  |     kriging.plot = function(canvas, grid, xlim, ylim, colors) { | ||
|  | 	// Clear screen
 | ||
|  | 	var ctx = canvas.getContext("2d"); | ||
|  | 	ctx.clearRect(0, 0, canvas.width, canvas.height); | ||
|  | 
 | ||
|  | 	// Starting boundaries
 | ||
|  | 	var range = [xlim[1]-xlim[0], ylim[1]-ylim[0], grid.zlim[1]-grid.zlim[0]]; | ||
|  | 	var i, j, x, y, z; | ||
|  | 	var n = grid.length; | ||
|  | 	var m = grid[0].length; | ||
|  | 	var wx = Math.ceil(grid.width*canvas.width/(xlim[1]-xlim[0])); | ||
|  | 	var wy = Math.ceil(grid.width*canvas.height/(ylim[1]-ylim[0])); | ||
|  | 	for(i=0;i<n;i++) | ||
|  | 	    for(j=0;j<m;j++) { | ||
|  | 		if(grid[i][j]==undefined) continue; | ||
|  | 		x = canvas.width*(i*grid.width+grid.xlim[0]-xlim[0])/range[0]; | ||
|  | 		y = canvas.height*(1-(j*grid.width+grid.ylim[0]-ylim[0])/range[1]); | ||
|  | 		z = (grid[i][j]-grid.zlim[0])/range[2]; | ||
|  | 		if(z<0.0) z = 0.0; | ||
|  | 		if(z>1.0) z = 1.0; | ||
|  | 
 | ||
|  | 		ctx.fillStyle = colors[Math.floor((colors.length-1)*z)]; | ||
|  | 		ctx.fillRect(Math.round(x-wx/2), Math.round(y-wy/2), wx, wy); | ||
|  | 	    } | ||
|  | 
 | ||
|  |     }; | ||
|  | 
 | ||
|  | 
 | ||
|  |     return kriging; | ||
|  | }(); | ||
|  | // if (module && module.exports){
 | ||
|  | //     module.exports = kriging;
 | ||
|  | // }
 |