31:Tc72,md`# Displaying 3D RBFs with ray marching A [Radial basis function (RBF)](https://en.wikipedia.org/wiki/Radial_basis_function) provides a neat way to solve the so-called *scattered data interpolation* problem. In a nutshell, suppose you want a function ${tex`f(X)`} about which you nothing other than its value at a few points of the the domain, i.e., you know that ${tex`f(X_i) = Y_i`} for sets of points ${tex`\{X_0, X_1 \cdots X_n \}`} and values ${tex`\{Y_0, Y_1 \cdots Y_n\}`}. The "radial" in RBF is because ${tex`f(X)`} depends on the distance from ${tex`X`} to points ${tex`X_i`}. The *kernel* of a RBF controls how distances are measured. For instance, a \`linear\` kernel means that the Euclidean norm ${tex`d (X) = \| X - X_i \|`} is used, whereas a \`quintic\` kernel employs the Euclidean norm to the fifth power, and so forth. A surface can be modeled using RBFs by adding points ${tex`X_i`} that we wish to interpolate and associating them with ${tex`Y_i = 0`}. In the implementation below, these are generated at random within a cube ${tex`-1 \leq x,y,z \leq 1`}. We must also add at least one point not associated with zero so that the equation system is fully determined. Here, we attempt to use [Ray Marching](http://jamie-wong.com/2016/07/15/ray-marching-signed-distance-functions/) to produce a rendering of such surfaces. The technique is used to produce renderings of *signed distance functions* (SDFs), i.e., scalar fields that map points in space to values that represent the distance to an object bounded by a closed surface, using the convention that points inside the surface are mapped to negative values, whereas points outside the object are mapped to positive values. In order to make our RBFs resemble an SDF, we add 8 points at the corners of a bigger cube and map them to the distance to the closest point that is mapped to zero. Unfortunately, this scheme alone does not quite work, simply because RBF surfaces *are not distance functions*. We only get a (sort of) decent picture for linear kernels and only for some point/value sets. In a last ditch attempt to make ray marching work with RBF surfaces using non-linear kernels, we introduce a *correction function* ${tex`g(f(X)) = a f(X)^b`} which is used to estimate the ray marching steps. In other words, for a ray with origin at ${tex`X`} such that ${tex`f(X) > 0`}, the origin for the next step is not set to ${tex`X^\prime = X + f(X)`} as with normal ray marching, but rather to ${tex`X^{\prime\prime} = X + g(f(X))`}. This sort of works in the sense that we can at least distinguish the shape of the surface, but the lighting still comes off wrong since the finite differences approach used to compute the gradient of *f* returns very imprecise values. In the interface below, you can select the desired kernel and the linear (*a*) and power (*b*) coefficients of *g* as well as the number of points used to generate the surface. In order to visualize the disparity between ${tex`f(X)`} and the actual distance between ${tex`X`} and the surface, we show a scatterplot of these for a regular sampling of the cube space. We also show in red the behavior of the correction function. `32:T12af,fragShader = ` precision highp float; const int MAX_MARCHING_STEPS = 255; const float MIN_DIST = 0.0; const float MAX_DIST = 100.0; const float EPSILON = 1.0e-3; const float fieldOfView = radians(45.0); // Vertical camera angle const vec3 lightPos = vec3(1.0, 1.0, 1.0); // Position of the light source const int maxn=${npoints + 8}; // Max number of points uniform vec2 u_resolution; // Size of the canvas uniform mat4 transform; // Rotation matrix //uniform int u_n; // Number of points / weights uniform float correctionPower; // power factor to correct rbf estimates of distance uniform float correctionLinear; // linear factor to correct rbf estimates of distance uniform vec3 u_points[maxn]; // Array of points uniform float u_weights[maxn]; // Array of point weights vec3 eye = 10.0*(transform*vec4(0.0, 0.,1.0,0.0)).xyz; // Position of the eye bool pointsFunc = false; // When true, ray marches the points sdf rather than the RBF sdf float kernel(float r) { return ${glslKernels[kernel]}; } float rbfSDF (vec3 samplePoint) { float total = 0.; for (int i = 0; i < maxn; i++) { //if (i>=u_n) return total; total += kernel(length(samplePoint-u_points[i]))*u_weights[i]; } return total; } float pointsSDF (vec3 p) { float d = MAX_DIST; for (int i = 8; i < maxn; i++) { //if (i>=u_n) return d; d = min(d,length(p-u_points[i])-0.05); } return d; } float cubeSDF(vec3 p) { vec3 d = abs(p) - vec3(1.0, 1.0, 1.0); float insideDistance = min(max(d.x, max(d.y, d.z)), 0.0); float outsideDistance = length(max(d, 0.0)); return insideDistance + outsideDistance; } float sphereSDF(vec3 samplePoint) { return length(samplePoint) - 1.0; } float sceneSDF(vec3 samplePoint) { return pointsFunc ? pointsSDF(samplePoint) : rbfSDF(samplePoint); } float shortestDistanceToSurface(vec3 marchingDirection, float start, float end) { float depth = start; for (int i = 0; i < MAX_MARCHING_STEPS; i++) { float dist = sceneSDF(eye + depth * marchingDirection); if (dist < EPSILON) { return depth; } if (!pointsFunc) dist = correctionLinear * pow(dist,correctionPower); depth += dist; if (depth >= end) { return end; } } return depth; } vec3 rayDirection(vec2 fragCoord) { vec2 xy = fragCoord - u_resolution / 2.0; float z = u_resolution.y / tan(fieldOfView / 2.0); return normalize(vec3(xy, -z)); } vec3 normal(vec3 p) { return normalize (vec3 (sceneSDF(p+vec3(EPSILON,0.0,0.0))-sceneSDF(p-vec3(EPSILON,0.0,0.0)), sceneSDF(p+vec3(0.0,EPSILON,0.0))-sceneSDF(p-vec3(0.0,EPSILON,0.0)), sceneSDF(p+vec3(0.0,0.0,EPSILON))-sceneSDF(p-vec3(0.0,0.0,EPSILON)))); } /** * Returns the lighting contribution of a single point light source via Phong illumination. * * k_a: Ambient color * k_d: Diffuse color * k_s: Specular color * alpha: Shininess coefficient * p: position of point being lit * N: normal at p * lightPos: the position of the light * lightIntensity: color/intensity of the light */ vec3 phong (vec3 k_a, vec3 k_d, vec3 k_s, float alpha, vec3 p, vec3 N, vec3 lightPos, vec3 lightIntensity) { vec3 L = normalize(lightPos - p); vec3 V = normalize(eye - p); vec3 R = normalize(reflect(-L, N)); float dotLN = abs(dot(L, N)); float dotRV = dot(R, V); if (dotLN < 0.0) { // Light not visible from this point on the surface return k_a * lightIntensity; } if (dotRV < 0.0) { // Light reflection in opposite direction as viewer, apply only diffuse // component return lightIntensity * (k_a + k_d * dotLN); } return lightIntensity * (k_a + k_d * dotLN + k_s * pow(dotRV, alpha)); } void main( ) { vec3 dir = (transform*(vec4(rayDirection(gl_FragCoord.xy),1.0))).xyz; pointsFunc = false; float distRBF = shortestDistanceToSurface(dir, MIN_DIST, MAX_DIST); pointsFunc = true; float distPoints = shortestDistanceToSurface(dir, MIN_DIST, MAX_DIST); float dist = min(distPoints, distRBF); if (dist > MAX_DIST - EPSILON) { discard; } else { vec3 p = eye + dist * dir; pointsFunc = dist == distPoints; vec3 N = normal(p); vec3 color = pointsFunc ? vec3(1.,1.,0.) : vec3(1.,0.,0.); gl_FragColor = vec4(phong (0.2*color, // ambient 0.5*color, // diffuse vec3 (0.5, 0.5, 0.5), // specular 50.0, // shininess p, // point N, // normal at point 10.0*lightPos, // light position vec3(1.0,1.0,1.0)), // light color/intensity 1.0); } }`33:Te9f,/** * This is essentially the rbf module by Thibaut Séguy, but modified to return the RBF weights, rather than a * function **/ RBF = { var distanceFunctions = { linear: distanceLinear, cubic: distanceCubic, quintic: distanceQuintic, "thin-plate": distanceThinPlate, gaussian: distanceGaussian, "inverse-multiquadric": distanceInverseMultiquadric, multiquadric: distanceMultiquadric }; var normFunctions = { euclidean: euclideanNorm }; var distance = distanceFunctions.linear; var norm = normFunctions.euclidean; function RBF(points, values, distanceFunction, normFunction, epsilon) { if (distanceFunction) { distance = typeof distanceFunction !== "string" ? distanceFunction : distanceFunctions[distanceFunction]; } if (normFunction) { norm = typeof normFunction !== "string" ? normFunction : normFunctions[normFunction]; } // `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 { w, interpolant }; } function euclideanNorm(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); } return RBF; }34:T475,class Arcball { constructor(canvas) { this.canvas = canvas; this.center = [canvas.width / 2, canvas.height / 2]; this.radius = (Math.min(canvas.width, canvas.height) / 2) * 0.9; this.prev = []; this.transform = mat4.create(); this.mouseIsDown = false; } ballPos(p) { let u = vec2.sub([], p, this.center); if (vec2.len(u) >= this.radius) return vec3.normalize([], [...u, 0]); let z = Math.sqrt(this.radius ** 2 - u[0] ** 2 - u[1] ** 2); return vec3.normalize([], [...u, z]); } mousedown(x, y) { this.mouseIsDown = true; this.prev = [x, y]; } mousemove(x, y) { if (!this.mouseIsDown) return this.transform; let mouse = [x, y]; let u = this.ballPos(this.prev); let v = this.ballPos(mouse); this.prev = mouse; let axis = vec3.cross([], v, u); let ang = Math.asin(vec3.len(axis)); let rot = mat4.fromRotation([], ang, vec3.normalize([], axis)); [this.u, this.v, this.ang, this.axis] = [u, v, ang, axis]; if (rot) mat4.mul(this.transform, this.transform, rot); return this.transform; } mouseup(x, y) { this.mouseIsDown = false; } }35:[] 2f:{"title":"Displaying 3D RBFs with ray marching","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":307,"value":"viewof kernel = Inputs.select([\"linear\", \"cubic\", \"quintic\", \"thin-plate\"], {\n label: \"RBF Kernel\",\n value: \"linear\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":623,"value":"viewof npoints = Inputs.select([4, 8, 16, 32], { label: \"points\", value: 8 })","pinned":false,"hidden":false,"mode":"ojs"},{"id":602,"value":"viewof correctionPower = Inputs.range([0.15, 1], {\n label: \"distance correction power\",\n value: correctionDefaults[kernel][0],\n step: 0.05\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":609,"value":"viewof correctionLinear = Inputs.range([0.1, 5], {\n label: \"distance correction linear coefficient\",\n value: correctionDefaults[kernel][1],\n step: 0.05\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":581,"value":"viewof randomize = html`