# Function Plotting in GLSL

Function plotting is a useful tool to visualize mathematical functions. There are many tools available online. Most of them are written in...

This blog requires WebGL enabled.

Function plotting is a useful tool to visualize mathematical functions. There are many tools available online. Most of them are written in JavaScript with the 2D canvas and line drawing. However, to plot the function in GLSL is not that straight forward. A good anti-aliasing algorithm is the main challenge.
The first approach is to color the pixel if the function passes through the pixel. We sample a few points inside the pixel, $x_{lower} \le x \le x_{upper}$, and if there is a point on the function, the pixel would be colored. Plotting $y=-x^2+2$

For example, to plot $y=-x^2+2$, for each pixel, if the function passes through the pixel, then we color it as red.

// +glslSnippet 300px 300px

uniform vec2 iResolution;

float f(float x) {  return -x*x+2.0;  }

void mainImage(out vec4 fragColor, vec2 fragCoord)
{
vec2 uv = fragCoord / iResolution.xy;
float xc = uv.x * 8.0 - 4.0;
float yc = uv.y * 8.0 - 4.0;
vec2 p = vec2(8.0/iResolution.x, 8.0/iResolution.y);
vec4 x = vec4(xc-p.x/2.0, xc-p.x/4.0, xc+p.x/4.0, xc+p.x/2.0);
vec4 y = vec4(f(x.x), f(x.y), f(x.z), f(x.w));
vec2 center = vec2(xc, f(xc));
vec2 yrange = vec2(yc-p.y/2.0, yc+p.y/2.0);

if((yrange.x <= y.x && y.x <= yrange.y) ||
(yrange.x <= y.y && y.y <= yrange.y) ||
(yrange.x <= center.y && center.y <= yrange.y) ||
(yrange.x <= y.z && y.z <= yrange.y) ||
(yrange.x <= y.w && y.w <= yrange.y)) {
fragColor = vec4(1.0, 0.0, 0.0, 1.0);
return;
}
fragColor = vec4(0.0, 0.0, 0.0, 1.0);
}

The first approach shows artifacts as the derivative of the function increases. Because we only sample 3 points inside the pixel, it also creates breakpoints on the function curve. Increasing the number of samples can reduce the breakpoints but it also significantly slows down the performance. However, it doesn't reduce the artifacts because this method does not differentiate the function passing through the center of the pixel or only a small corner of the pixel.

Let's take a closer look at the artifacts. Ideally, instead of On or Off, we need to weight each pixel, which averages the color and renders a smoother result.

Obviously, for the above image, case 1 and case 2 shall have the highest weight. Case 2 also shows, if the derivative of the function is very large, we are not able to increase the sample get better results, as to have $y$ in the pixel, the $\Delta x$ needs to be very small. Instead of using the samples, the distance to the function will be used as the weight of the given pixel.

// +glslSnippet 300px 300px

uniform vec2 iResolution;

float f(float x) {  return -x*x+2.0;  }

void mainImage(out vec4 fragColor, vec2 fragCoord)
{
vec2 uv = fragCoord / iResolution.xy;
float xc = uv.x * 8.0 - 4.0;
float yc = uv.y * 8.0 - 4.0;
vec2 p = vec2(xc, f(xc));
float dist = length(p - vec2(xc, yc));
float halfWidth = 4.0 / iResolution.x;
float intensity = smoothstep(0., 1., 1.-dist/halfWidth/5.0);
fragColor = vec4(intensity, 0.0, 0.0, 1.0);
}

The second approach is to measure the distance from the center of the pixel to the function value. And we use the distance as a weight to slowly reduce the intensity of the pixel. When the derivative of the function is small, it renders a perfect result, but it fails at the higher derivative part.

The problem is the distance function. In the second approach, we use point to point distance, measuring the center of the pixel to the point on the function. But for the high derivative region, the $y$ value alters very fast. It then results in breakpoints on the graph. Replacing the point distance with tangent distance can easily solve the problem.

// +glslSnippet 300px 300px

uniform vec2 iResolution;

float DistanceToLine(vec2 p0, vec2 p1, vec2 p)
{
// Transform the point to screen space
p0 = (p0 + vec2(4.0, 4.0)) / (vec2(4.0, 4.0) * 2.) * iResolution.xy;
p1 = (p1 + vec2(4.0, 4.0)) / (vec2(4.0, 4.0) * 2.) * iResolution.xy;
p = (p + vec2(4.0, 4.0)) / (vec2(4.0, 4.0) * 2.) * iResolution.xy;
vec2 ld = p0 - p1;
vec2 pd = p - p1;
return length(pd - ld*clamp(dot(pd, ld)/dot(ld, ld), 0.0, 1.0));
}

float f(float x) {  return -x*x+2.0;  }

void mainImage(out vec4 fragColor, vec2 fragCoord)
{
vec2 uv = fragCoord / iResolution.xy;
float width = 8. / iResolution.x;
vec2 p = vec2(uv.x*8.-4., uv.y*8.-4.);
vec2 p0 = vec2(uv.x*8.-4.-width, 0.0);
vec2 p1 = vec2(p0.x+width, 0.0);
vec2 p2 = vec2(p1.x+width, 0.0);
p0 = vec2(p0.x, f(p0.x));
p1 = vec2(p1.x, f(p1.x));
p2 = vec2(p2.x, f(p2.x));
float dist0 = DistanceToLine(p0, p1, p);
float dist1 = DistanceToLine(p1, p2, p);

float intensity = smoothstep(0., 1., 1.-min(dist0, dist1));
fragColor = vec4(intensity, 0.0, 0.0, 1.0);
}

The last approach renders a nice graph. If we take more samples for the tangent, we can get even better results.

Here are more examples:

// +glslSnippet 100% 16:9

uniform vec2 iResolution;

float DistanceToLine(vec2 iPlot, vec2 p0, vec2 p1, vec2 p)
{
p0 = (p0+iPlot/2.) / iPlot * iResolution.xy;
p1 = (p1+iPlot/2.) / iPlot * iResolution.xy;
p = (p+iPlot/2.) / iPlot * iResolution.xy;
vec2 ld = p0 - p1;
vec2 pd = p - p1;
return length(pd - ld*clamp(dot(pd, ld)/dot(ld, ld), 0.0, 1.0));
}

void GetPoints(vec2 iPlot, vec2 fragCoord, out vec4 r0, out vec4 r1)
{
vec2 uv = fragCoord / iResolution.xy;
float width = iPlot.x / iResolution.x;
vec2 p = vec2(uv.x*iPlot.x-iPlot.x/2., uv.y*iPlot.y-iPlot.y/2.);
vec2 p0 = vec2(p.x-width, 0.0);
vec2 p1 = vec2(p.x, 0.0);
vec2 p2 = vec2(p.x+width, 0.0);
r0 = vec4(p, p0);
r1 = vec4(p1, p2);
}

float Attenuate(vec2 iPlot, vec4 r0, vec4 r1)
{
float dist0 = DistanceToLine(iPlot, r0.zw, r1.xy, r0.xy);
float dist1 = DistanceToLine(iPlot, r1.xy, r1.zw, r0.xy);
return smoothstep(0., 1., 1.-min(dist0, dist1));
}

float f(float x) {  return sin(3.14*x)/(3.14*x);  }

void mainImage(out vec4 fragColor, vec2 fragCoord)
{
vec2 iPlot = vec2(20., 20.*iResolution.y/iResolution.x);
vec4 r0, r1;

vec3 f1 = vec3(1., 0., 0.);
GetPoints(iPlot, fragCoord, r0, r1);
r0.w = sin(r0.z); r1.y = sin(r1.x); r1.w = sin(r1.z);
f1 *= Attenuate(iPlot, r0, r1);

vec3 f2 = vec3(0., 1., 0.);
GetPoints(iPlot, fragCoord, r0, r1);
r0.w = cos(r0.z); r1.y = cos(r1.x); r1.w = cos(r1.z);
f2 *= Attenuate(iPlot, r0, r1);

vec3 f3 = vec3(0., 0., 1.);
GetPoints(iPlot, fragCoord, r0, r1);
r0.w = tan(r0.z); r1.y = tan(r1.x); r1.w = tan(r1.z);
f3 *= Attenuate(iPlot, r0, r1);

vec3 f4 = vec3(1., 0., 1.);
GetPoints(iPlot, fragCoord, r0, r1);
r0.w = f(r0.z); r1.y = f(r1.x); r1.w = f(r1.z);
f4 *= Attenuate(iPlot, r0, r1);

fragColor = vec4(f1+f2+f3+f4, 1.0);
}

Have fun :)