EXPLAIN EXTENDED

In my New Year posts I usually try to recap and summarize the past year. It won't take long this time:

Fuck you, coronavirus!

Now that I've gotten that off my chest, I have to think of something to write about in this New Year's post.

So I was thinking, why not put a face to the name we all hate so much?

Let's use SQL to do some ray tracing and draw a 3D picture of the dreaded virus.

By now, I believe we are all familiar with the picture of the virus. It looks like a ball covered with spikes. The spikes look something like the solar corona, which is what gave the virus its name. They have this distinct triangular shape.

We'll create a sphere covered with several dozens of spikes.

Every spike will be a small pyramid, with an equilateral triangle as a base and isosceles triangles as lateral faces. This means it will be a right pyramid.

The pyramids will be "standing" on their apexes, upside-down. The height of every pyramid will be perpendicular to the sphere surface and continue the sphere's radius at the apex.

Then we will implement the pinhole camera model and use ray tracing algorithms to calculate the lighting of the sphere and the spikes.

Types

3D modeling heavily uses vector algebra. Of course pure SQL offers enough math functions to get around. But functions and routines are not first class citizens in SQL, which means we would have to copy-paste the bulky vector manipulation formulas every time we will need them, which would make our query unwieldy really fast.

This is a good chance to get familiar with PostgreSQL's rich system of custom types and custom operators. It lets users define their own types, create functions to work with them and even overload the operators.

Before we start working with our model, we will create a type for 3D vectors and several functions and operators to work with that type.

CREATE TYPE VEC3 AS (x DOUBLE PRECISION, y DOUBLE PRECISION, z DOUBLE PRECISION);



CREATE OR REPLACE FUNCTION VEC3_ADD (a VEC3, b VEC3)
RETURNS vec3
AS
$$
SELECT (a.x + b.x, a.y + b.y, a.z + b.z)::VEC3;
$$
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

CREATE OPERATOR + (FUNCTION = VEC3_ADD, LEFTARG = VEC3, RIGHTARG = VEC3, COMMUTATOR = +);



CREATE OR REPLACE FUNCTION VEC3_NEG(a VEC3)
RETURNS VEC3
AS
$$
SELECT  (-a.x, -a.y, -a.z)::VEC3;
$$
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

CREATE OPERATOR - (FUNCTION = VEC3_NEG, RIGHTARG = VEC3);



CREATE OR REPLACE FUNCTION VEC3_SUB(a VEC3, b VEC3)
RETURNS VEC3
AS
$$
SELECT (a.x - b.x, a.y - b.y, a.z - b.z)::VEC3;
$$
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

CREATE OPERATOR - (FUNCTION = VEC3_SUB, RIGHTARG = VEC3, LEFTARG = VEC3);



CREATE OR REPLACE FUNCTION VEC3_MUL(a VEC3, k DOUBLE PRECISION)
RETURNS VEC3
AS
$$
SELECT  (a.x * k, a.y * k, a.z * k)::VEC3;
$$
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

CREATE OR REPLACE FUNCTION VEC3_MUL(k DOUBLE PRECISION, a VEC3)
RETURNS VEC3
AS
$$
SELECT  VEC3_MUL(a, k);
$$
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

CREATE OPERATOR * (FUNCTION = VEC3_MUL, LEFTARG = VEC3, RIGHTARG = DOUBLE PRECISION, COMMUTATOR = *);
CREATE OPERATOR * (FUNCTION = VEC3_MUL, LEFTARG = DOUBLE PRECISION, RIGHTARG = VEC3, COMMUTATOR = *);



CREATE OR REPLACE FUNCTION VEC3_DIV(a VEC3, k DOUBLE PRECISION)
RETURNS VEC3
AS
$$
SELECT  (a * (1 / k))::VEC3
$$
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

CREATE OPERATOR / (FUNCTION = VEC3_DIV, LEFTARG = VEC3, RIGHTARG = DOUBLE PRECISION);



CREATE OR REPLACE FUNCTION VEC3_DOT(a VEC3, b VEC3)
RETURNS DOUBLE PRECISION
AS
$$
SELECT  a.x * b.x + a.y * b.y + a.z * b.z;
$$
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

CREATE OPERATOR * (FUNCTION = VEC3_DOT, LEFTARG = VEC3, RIGHTARG = VEC3, COMMUTATOR = *);



CREATE OR REPLACE FUNCTION VEC3_MODULE(a VEC3)
RETURNS DOUBLE PRECISION
AS
$$
SELECT  SQRT(a.x^2 + a.y^2 + a.z^2)
$$
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

CREATE OPERATOR | (FUNCTION = VEC3_MODULE, RIGHTARG = VEC3);



CREATE OR REPLACE FUNCTION VEC3_NORM(a VEC3)
RETURNS VEC3
AS
$$
SELECT  a / |a
$$
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

CREATE OPERATOR || (FUNCTION = VEC3_NORM, RIGHTARG = VEC3);



CREATE OR REPLACE FUNCTION VEC3_CROSS(a VEC3, b VEC3)
RETURNS VEC3
AS
$$
SELECT  (a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x)::VEC3
$$
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

CREATE OPERATOR ** (FUNCTION = VEC3_CROSS, LEFTARG = VEC3, RIGHTARG = VEC3);



CREATE TYPE SPHERICAL AS (radius DOUBLE PRECISION, theta DOUBLE PRECISION, phi DOUBLE PRECISION);



CREATE OR REPLACE FUNCTION VEC3_TO_SPHERICAL (a VEC3)
RETURNS SPHERICAL
AS
$$
SELECT  (|a, CASE WHEN |a = 0 THEN 0 ELSE ACOS(a.z / |a) END, ATAN2(a.y, a.x))::SPHERICAL
$$
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

CREATE CAST (VEC3 AS SPHERICAL) WITH FUNCTION VEC3_TO_SPHERICAL(VEC3) AS ASSIGNMENT;



CREATE OR REPLACE FUNCTION SPHERICAL_TO_VEC3 (a SPHERICAL)
RETURNS VEC3
AS
$$
SELECT  a.radius * (SIN(a.theta) * COS(a.phi), SIN(a.theta) * SIN(a.phi), COS(a.theta))::VEC3
$$
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT;

CREATE CAST (SPHERICAL AS VEC3) WITH FUNCTION SPHERICAL_TO_VEC3(SPHERICAL) AS ASSIGNMENT;

The type VEC3 is just a fancy name for a tuple of x, y, and z.

The operators +, -, *, / do scalar arithmetic operations on the vector, they accept a vector and a scalar.

The operators * and ** are dot product and cross product, respectively. They both accept two vectors on the either side, but the first one returns a number, and the second one returns another vector.

The unary operator | calculates the vector module (length), and || normalizes the vector (returns a vector with the same direction and the module of 1).

Finally, we'll create another type, SPHERICAL, which is a representation of a vector in spherical coordinates. SPHERICAL and VEC3 can be cast into each other, as they are just different representations of the same thing. We will need both of them, because some operations are easier in cartesian representation, and some are in spherical.

The apexes

Modeling a sphere is quite easy: all we need is its radius and the vector for its center.

For the spikes, we need to come up with a list of points on the sphere where the apexes will be. Then we'll need to create the base at a certain distance from the apex.

For our model to look close enough to the real thing, we need the apexes to be distributed more or less evenly around the sphere. I'm saying "more or less" because this problem, known as Thomson problem, does not have an exact solution except for some particular cases.

We don't need an exact solution, though, because we're modeling nature and it's a long way from being exact. All we need is a good enough solution. Fortunately, good enough is something nature itself has long since figured out. Turns out, there's a really simple algorithm, called the Fibonacci lattice, which allows to place points on a disc or a sphere so that their average density remains about the same.

The trick is to place the points on a Fermat's spiral at a certain fixed angle between every two adjacent points. This angle, called the golden angle, divides the circle in the golden ratio. So the points having the polar coordinates (k * √i, i * ϑ), where ϑ is the golden angle, will having density quite close to uniform. This is how the sunflower grows its florets. It makes sure that all the seeds will have equal space.

The same principle works for the sphere, except that the flat spiral has to be mapped onto a 3-dimensional surface. To preserve the density, we need to map the spiral in such a way that the flat area becomes proportional to the sphere area. We will need to use the cylindrical equal-area projection for that.

The polar coordinates on the plane will map to these spherical coordinates: (1, ACOS(2 * i / (N + 1) - 1), i * ϑ), where N is the number of points.

Let's build the set of 40 apexes for our future spikes on the unit sphere:

WITH    spikes AS
        (
        SELECT  i, spherical.*
        FROM    (
                SELECT  40 AS spikes
                ) constants
        CROSS JOIN LATERAL
                GENERATE_SERIES(0, spikes - 1) AS i
        CROSS JOIN LATERAL
                (
                SELECT  ACOS(2 * i / spikes::DOUBLE PRECISION - 1) AS theta,
                        PI() * (3 - SQRT(5)) * i AS phi
                ) spherical
        )
SELECT  *
FROM    spikes;

Let's see how good our distribution is. To do this, we'll take a hemisphere with the center in each apex, and see how many spikes are located on this hemisphere. If the distribution is uniform, every hemisphere should contain half of all the spikes.

All the points on the hemisphere are within the straight line distance of √2 from its center. So we'll just have to convert the spherical coordinates into cartesian (using our cool types and conversion operators) and select all of them which are within the distance. The distance between two apexes is the module of their vector difference.

Let's make it into a query:

WITH    spikes AS
        (
        SELECT  i, spherical.*
        FROM    (
                SELECT  40 AS spikes
                ) constants
        CROSS JOIN LATERAL
                GENERATE_SERIES(0, spikes - 1) AS i
        CROSS JOIN LATERAL
                (
                SELECT  ACOS(2 * i / spikes::DOUBLE PRECISION - 1) AS theta,
                        PI() * (3 - SQRT(5)) * i AS phi
                ) spherical
        )
SELECT  MIN(hemisphere_neighbors), MAX(hemisphere_neighbors)
FROM    spikes s1
CROSS JOIN LATERAL
        (
        SELECT  COUNT(*) hemisphere_neighbors
        FROM    spikes s2
        WHERE   |((1, s1.theta, s1.phi)::SPHERICAL::VEC3 - (1, s2.theta, s2.phi)::SPHERICAL::VEC3) <= SQRT(2)
        ) s2

It looks like every hemisphere contains from 18 to 21 apexes (out of 40 total) which is just about all right for our purpose.

The spikes

We have our apexes and we need to build the spike models from them. Every spike will be defined by 4 triangles: the base and the three faces of the pyramid. Every triangle will be defined as a set of three vectors, its vertices. So we'll have 4 records, three fields in each, for every spike.

The easiest way to build 40 spikes is to build one "prototype" spike and then clone it.

For the reasons that will become obvious later, it's best to build the prototype spike on the apex with both inclination and azimuth being zero: (1, 0, 0) in spherical coordinates. Note that this apex is not on our list. It is just a prototype, which will not make it to the final model.

We'll make our spikes 30% of the sphere's radius long and 10% wide. This way, the will be big enough to look cool on an ASCII art page, while still preserving that distinctive triangular shape.

The model apex has cartesian coordinates (0, 0, 1) which means that its radius is the Z-axis. The altitude (or height) dropped from the apex to the base should continue the radius, and its length is 0.3, so the center of the base will have the coordinates (0, 0, 1.3). The vertices of the base have to be at the distance of 0.1 from the center of the base and should be 120 degrees apart from each other. The base is perpendicular to the height (which is the Z-axis), so it's parallel to the XY plane.

It means that the Z component of the all the base vertices will be the same (i.e. 1.3), and the X and Y components will be points on the 0.1 circle around the origin, at 0, 120 and 240 degrees.

Let's create the prototype spike:

WITH    spike_parameters AS
        (
        SELECT  0.3 AS height, 0.1 AS width
        ),
        spike_origin_vertices AS
        (
        SELECT  1 AS j, (0, 0, 1)::VEC3 AS origin_vertex
        UNION ALL
        SELECT  i + 2 AS j, (0, 0, 1 + height)::VEC3 + (SIN(2 * PI() * i / 3), COS(2 * PI() * i / 3), 0)::VEC3 * width AS origin_vertex
        FROM    spike_parameters
        CROSS JOIN
                GENERATE_SERIES(0, 2) i
        )
SELECT  *
FROM    spike_origin_vertices

Barring the rounding errors, the vertices look about right: the apex is on the sphere, the base vertices are 0.1 units away from the Z-axis, all having the same 1.3 units Z component.

Now we have to take this prototype and graft it on every apex in the Fibonacci lattice we had created earlier. We also need to make sure the base is continuing the radius at the apex, and adjust all the angles accordingly.

One way to do that is to rotate the prototype about the center of the sphere, so that prototype apex becomes the new apex. It's doable, but the formula for calculating this rotation (the Rodrigues' rotation formula) is very bulky, and we would also need an extra step to calculate the axis of the rotation.

However, we can make the rotation much easier, if we do it in two steps. Every apex is a point on the unit sphere, so it will have the radius 1, and some inclination and azimuth: (1, θ, φ). Now, remember how we located our prototype apex at (1, 0, 0)? This means that we can rotate it in two easy, successive steps: first, about the Y-axis by θ (the inclination), then about the Z-axis by φ (the azimuth): (1, 0, 0) → (1, θ, 0) → (1, θ, φ). The rest of the vertices will tag along.

A rotation about a basis axis is much easier (since it's effectively 2-d), and we already know the angles (they come straight from the Fibonacci lattice formula).

Here's the query to clone all the vertices:

WITH    spikes AS
        (
        SELECT  i, spherical.*
        FROM    (
                SELECT  40 AS spikes
                ) constants
        CROSS JOIN LATERAL
                GENERATE_SERIES(0, spikes - 1) AS i
        CROSS JOIN LATERAL
                (
                SELECT  ACOS(2 * i / spikes::DOUBLE PRECISION - 1) AS theta,
                        PI() * (3 - SQRT(5)) * i AS phi
                ) spherical
        ),
        spike_parameters AS
        (
        SELECT  0.3 AS height, 0.1 AS width
        ),
        spike_origin_vertices AS
        (
        SELECT  1 AS j, (0, 0, 1)::VEC3 AS origin_vertex
        UNION ALL
        SELECT  i + 2 AS j, (0, 0, 1 + height)::VEC3 + (SIN(2 * PI() * i / 3), COS(2 * PI() * i / 3), 0)::VEC3 * width AS origin_vertex
        FROM    spike_parameters
        CROSS JOIN
                GENERATE_SERIES(0, 2) i
        ),
        spike_vertices AS
        (
        SELECT  i, rotated
        FROM    spikes
        CROSS JOIN
                spike_origin_vertices
        CROSS JOIN LATERAL
                (
                SELECT  ((origin_vertex).x * COS(theta) + (origin_vertex).z * SIN(theta), (origin_vertex).y, (origin_vertex).z * COS(theta) - (origin_vertex).x * SIN(theta))::VEC3 AS rotated_y
                ) y
        CROSS JOIN LATERAL
                (
                SELECT  ((rotated_y).x * COS(phi) - (rotated_y).y * SIN(phi), (rotated_y).y * COS(phi) + (rotated_y).x * SIN(phi), (rotated_y).z)::VEC3 AS rotated
                ) rotated
        )
SELECT  *
FROM    spike_vertices

For every spike, we have a set of 4 vertices that define it.

Let's combine these vertices into triangles. A triangular pyramid has this nice property that every three of its vertices make a face, so we just need to take all 3-combinations of 4 vertices.

Here's the query:

WITH    spikes AS
        (
        SELECT  i, spherical.*
        FROM    (
                SELECT  40 AS spikes
                ) constants
        CROSS JOIN LATERAL
                GENERATE_SERIES(0, spikes - 1) AS i
        CROSS JOIN LATERAL
                (
                SELECT  ACOS(2 * i / spikes::DOUBLE PRECISION - 1) AS theta,
                        PI() * (3 - SQRT(5)) * i AS phi
                ) spherical
        ),
        spike_parameters AS
        (
        SELECT  0.3 AS height, 0.1 AS width
        ),
        spike_origin_vertices AS
        (
        SELECT  1 AS j, (0, 0, 1)::VEC3 AS origin_vertex
        UNION ALL
        SELECT  i + 2 AS j, (0, 0, 1 + height)::VEC3 + (SIN(2 * PI() * i / 3), COS(2 * PI() * i / 3), 0)::VEC3 * width AS origin_vertex
        FROM    spike_parameters
        CROSS JOIN
                GENERATE_SERIES(0, 2) i
        ),
        spike_vertices AS
        (
        SELECT  i, j, rotated
        FROM    spikes
        CROSS JOIN
                spike_origin_vertices
        CROSS JOIN LATERAL
                (
                SELECT  ((origin_vertex).x * COS(theta) + (origin_vertex).z * SIN(theta), (origin_vertex).y, (origin_vertex).z * COS(theta) - (origin_vertex).x * SIN(theta))::VEC3 AS rotated_y
                ) y
        CROSS JOIN LATERAL
                (
                SELECT  ((rotated_y).x * COS(phi) - (rotated_y).y * SIN(phi), (rotated_y).y * COS(phi) + (rotated_y).x * SIN(phi), (rotated_y).z)::VEC3 AS rotated
                ) rotated
        ),
        spike_triangles AS
        (
        SELECT  i, sv1.rotated AS one, sv2.rotated AS two, sv3.rotated AS three
        FROM    spike_vertices sv1
        JOIN    spike_vertices sv2
        USING   (i)
        JOIN    spike_vertices sv3
        USING   (i)
        WHERE   sv2.j > sv1.j
                AND sv3.j > sv2.j
        )
SELECT  *
FROM    spike_triangles

Ray tracing

For the ray tracing, we will use the most simple mapping: the pinhole camera model. Every pixel on the resulting image would be lit by exactly one point on the surface of the model. There are no lenses, which means there are no geometric distortions on the image: all straight lines map to straight lines.

In a real pinhole camera, the image plane is located behind the pinhole, and the image is reversed. However, we're dealing with a mathematical model, so we are free to locate the image in front of the pinhole (or the viewpoint). All the mathematics stay the same, but this way, the image is not reversed.

The brightness of every pixel on our image is defined by intensity of light which the corresponding point on the model surface reflects towards this pixel.

For our model, we will be using a single source of light, located quite far away from the scene. We will also assume some amount of ambient light, less bright than the main source, so that we could see the parts of our model which are not being directly illuminated by the source.

Here's what we need to do to calculate the brightness of every pixel:

1. Define if the pixel ray intersects any part of the model. If it doesn't, there is nothing to reflect the light, and the pixel is completely black
2. If the pixel ray intersects several parts of the model, find out which one is the closest. Since it's blocking the light from the other surfaces, the closest surface will be the only one which defines the pixel intensity
3. Find out the amount of light reflected by surface at the intersection point in the direction of the pixel ray.

Intersections

For every pixel, we know the parameters of the pixel ray: it can be defined by a unit vector from the origin (the viewpoint) towards the pixel. We have two types of objects in our model: spheres (one sphere, actually) and triangles. For every type of object, we know its parameters: the radius and the center for the sphere, and the coordinates of all the vertices for the triangles.

So we have everything we need to implement two different algorithms: intersection of a ray with a sphere and intersection of a ray with a triangle.

Intersections with the sphere

To figure out whether or not the ray intersects with the sphere, we need to calculate the distances from the center of the sphere to the ray. This distance is equal to the length of the perpendicular dropped from the center of the sphere. We need to see if its length is more, equal or less than the sphere's radius. If it's more, then the ray does not intersect the sphere. If it's equal, the ray touches the sphere (intersects it in exactly one point). Finally, if it's less, then the ray intersects the sphere in two points.

The intersection point can be presented in the vector form as the origin point (V) plus pixel unit vector (P) times the distance from the viewpoint (d): (V̅ + d * P̅). The difference between this vector and the sphere center vector is another vector, whose module should be equal to the radius of the sphere.

This way, we get a quadratic equation in regard to d: |V̅ + d * P̅ - C̅|2 = R2, or d2 + 2 * d * P̅(V̅ - C̅) + (|V̅ - C̅|2 - R2) = 0

We know that if the discriminant of this equation is less than zero, there is no intersection. We can use this fact to make our query a little bit more efficient: we will calculate the determinant first (we need to do that anyway to solve the equation), and just skip the rows where it's less than zero.

For the pixels whose rays do intersect the sphere, we need to find the closest intersection, because this is the only one which we will need. The farthest one, first, will face away from the pixel, and, second, will be blocked by the closest one. This way we can optimize our query yet a little bit.

Let's find all the intersection points for every pixel.

We'll generate a grid of pixels first, at the distance of 100 units from the origin. For each pixel, we'll calculate a unit vector from the viewpoint.

Then we will find the discriminant of the quadratic equation above and throw away the pixels with the negative discriminant.

We then will solve the quadratic equation above, throwing away all the negative results (not that there's gonna be any). A negative result would mean an intersection behind the pixel grid, and all of the sphere is all in front of the grid. But we will want to keep our query tidy. Every solution would represent the distance to the intersection point.

Finally, if there's more that one solution for a pixel, we select the smallest one, because this one would be the closest intersection.

Here's the query (limiting the results for demonstration purposes):

WITH    origin AS
        (
        SELECT  (0, 0, 0)::VEC3 AS origin
        ),
        sphere AS
        (
        SELECT 120 AS radius, (0, 0, 300)::VEC3 AS center
        ),
        pixels AS
        (
        SELECT  *
        FROM    GENERATE_SERIES(-70, 70) x
        CROSS JOIN
                GENERATE_SERIES(-70, 70) y
        CROSS JOIN LATERAL
                (
                SELECT  ||((x, y, 100)::VEC3) AS pixel_unit
                ) pixel_unit
        ),
        sphere_intersection_coefficients AS
        (
        SELECT  *
        FROM    pixels
        CROSS JOIN
                sphere
        CROSS JOIN
                origin
        CROSS JOIN LATERAL
                (
                SELECT  1 AS a, 2 * pixel_unit * (origin - center) AS b, (|(origin - center))^2 - radius^2 AS c
                ) q3
        CROSS JOIN LATERAL
                (
                SELECT b ^ 2 - 4 * a * c AS discriminant
                ) q4
        WHERE   discriminant > 0
        ),
        sphere_intersections AS
        (
        SELECT  x, y, pixel_unit, intersection_distance
        FROM    sphere_intersection_coefficients
        CROSS JOIN LATERAL
                (
                SELECT  t AS intersection_distance
                FROM    (VALUES (1), (-1)) q (sign)
                CROSS JOIN LATERAL
                        (
                        SELECT  (-b + SQRT(discriminant) * sign) / (2 * a) AS t
                        WHERE   discriminant > 0
                        ) q2
                WHERE   t > 0
                ORDER BY
                        t
                LIMIT   1
                ) q
        )
SELECT  *
FROM    sphere_intersections
LIMIT 10;

We got the pixel unit and distance to every intersection with the sphere, which is enough to reconstruct the cartesian coordinates for the intersection points. We will need them in the future.

Intersections with the triangles

Remember the spikes? We have 40 spikes, each one made of 4 triangles. This means we have 160 triangles. Same as with the sphere, we will need to find the intersections of pixel rays with every one of these triangles.

There is quite an efficient algorithm for finding the intersections of a ray with a triangle. It is called Möller–Trumbore intersection algorithm.

Those interested in the math and principles behind it can read its detailed description on Scratchapixel: Möller-Trumbore algorithm. Scratchapixel is an awesome online book about computer graphics in general and ray tracing in particular.

I will just copy the implementation from Wikipedia and adapt it for SQL. This is quite a straightforward algorithm, with no recursion and almost no branching. It translates to SQL very easily. It's a little bit heavy on math, though.

When we first constructed our spikes, we had modeled them around a unit sphere. This time, we will have to model them around the actual sphere. This is very easy using the vector algebra. We will multiply every spike vertex vector by the radius (enlargement) and add the sphere center vector to it (translation). This is doable just with our basic vector operators, without any fancy matrix calculations. I'll add one more field to the CTE called spike_vertices to do that.

Now that we know the pixel ray parameters and spike parameters, we can just feed them to the algorithm and get the list of intersections back. As in the previous case, it will return us the list of distances to the intersection for every pixel.

Note that unlike the sphere, a ray can have zero, one or infinity intersections with a triangle. Infinity means that the ray lies in the same plane as the triangle and the intersection is not a point, but a line segment. Such triangles would be invisible, as they are facing the viewer with their infinitely thin, non-existing third dimension. We will have to throw away cases like that, which is something that the algorithm takes care of.

For the demonstration purposes, we'll take just one triangle out of just one spike.

Here's the query:

WITH    origin AS
        (
        SELECT  (0, 0, 0)::VEC3 AS origin
        ),
        pixels AS
        (
        SELECT  *
        FROM    GENERATE_SERIES(-70, 70) x
        CROSS JOIN
                GENERATE_SERIES(-70, 70) y
        CROSS JOIN LATERAL
                (
                SELECT  ||((x, y, 100)::VEC3) AS pixel_unit
                ) pixel_unit
        ),
        sphere AS
        (
        SELECT 120 AS radius, (0, 0, 300)::VEC3 AS center
        ),
        spikes AS
        (
        SELECT  i, spherical.*
        FROM    (
                SELECT  40 AS spikes
                ) constants
        CROSS JOIN LATERAL
                GENERATE_SERIES(0, spikes - 1) AS i
        CROSS JOIN LATERAL
                (
                SELECT  ACOS(2 * i / spikes::DOUBLE PRECISION - 1) AS theta,
                        PI() * (3 - SQRT(5)) * i AS phi
                ) spherical
        ),
        spike_parameters AS
        (
        SELECT  0.3 AS height, 0.1 AS width
        ),
        spike_origin_vertices AS
        (
        SELECT  1 AS j, (0, 0, 1)::VEC3 AS origin_vertex
        UNION ALL
        SELECT  i + 2 AS j, (0, 0, 1 + height)::VEC3 + (SIN(2 * PI() * i / 3), COS(2 * PI() * i / 3), 0)::VEC3 * width AS origin_vertex
        FROM    spike_parameters
        CROSS JOIN
                GENERATE_SERIES(0, 2) i
        ),
        spike_vertices AS
        (
        SELECT  i, j, vertex
        FROM    spikes
        CROSS JOIN
                spike_origin_vertices
        CROSS JOIN LATERAL
                (
                SELECT  ((origin_vertex).x * COS(theta) + (origin_vertex).z * SIN(theta), (origin_vertex).y, (origin_vertex).z * COS(theta) - (origin_vertex).x * SIN(theta))::VEC3 AS rotated_y
                ) y
        CROSS JOIN LATERAL
                (
                SELECT  ((rotated_y).x * COS(phi) - (rotated_y).y * SIN(phi), (rotated_y).y * COS(phi) + (rotated_y).x * SIN(phi), (rotated_y).z)::VEC3 AS rotated
                ) rotated
        CROSS JOIN LATERAL
                (
                SELECT  (sphere.radius * rotated) + sphere.center AS vertex
                FROM    sphere
                ) mapped
        ),
        spike_triangles AS
        (
        SELECT  i,
                sv1.j + sv2.j + sv3.j - 5 AS j,
                sv1.vertex AS one, sv2.vertex AS two, sv3.vertex AS three
        FROM    spike_vertices sv1
        JOIN    spike_vertices sv2
        USING   (i)
        JOIN    spike_vertices sv3
        USING   (i)
        WHERE   sv2.j > sv1.j
                AND sv3.j > sv2.j
        ),
        spike_intersections AS
        (
        SELECT  i, j, one, two, three, x, y, pixel_unit, intersection_distance
        FROM    pixels
        CROSS JOIN
                origin
        CROSS JOIN
                spike_triangles
        CROSS JOIN LATERAL
                (
                SELECT  two - one AS edge1, three - one AS edge2
                ) edges
        CROSS JOIN LATERAL
                (
                SELECT  pixel_unit ** edge2 AS p_vector
                ) normal
        CROSS JOIN LATERAL
                (
                SELECT  edge1 * p_vector AS determinant
                ) det
        CROSS JOIN LATERAL
                (
                SELECT  origin - one AS t_vector
                ) t
        CROSS JOIN LATERAL
                (
                SELECT  (t_vector * p_vector) / determinant AS u
                WHERE   determinant NOT BETWEEN -1e-8 AND 1e-8
                ) u
        CROSS JOIN LATERAL
                (
                SELECT  t_vector ** edge1 AS q_vector
                WHERE   u BETWEEN 0 AND 1
                ) q
        CROSS JOIN LATERAL
                (
                SELECT  (pixel_unit * q_vector) / determinant AS v
                ) v
        CROSS JOIN LATERAL
                (
                SELECT  (edge2 * q_vector) / determinant AS intersection_distance
                WHERE   v >= 0 AND u + v <= 1
                ) intersection
        )
SELECT  *
FROM    spike_intersections
WHERE   i = 39
        AND j = 1;

This is just one triangle from the farthest spike from the point of view, and we only got seven intersections (which means that if this triangle is visible at all and not blocked by something else, it will occupy only 7 pixels on the screen). But if we look at the values of x and y of the intersections, we can see that it will have a distinctive triangular shape.

Combining the intersections

For each pixel, we have a list of its ray intersections with different objects (the sphere and the triangles). It is possible that one ray intersects with different objects. In this case, we should get rid of the duplicates and just leave one intersection, the closest one to the viewpoint. This is a usual max-row-per-group task. We can do that by combining the queries with UNION ALL, assigning every row a ROW_NUMBER() with partition by the pixel and ordering by the distance, and getting rid of all records with the row numbers higher than 1.

For demonstration purposes, we will just list the intersection for a narrow 10x2 rectangle in the center of the screen. We will also output the origin of each intersection (sphere or spike triangle), and, if the latter, which spike and which triangle.

WITH    origin AS
        (
        SELECT  (0, 0, 0)::VEC3 AS origin
        ),
        pixels AS
        (
        SELECT  *
        FROM    GENERATE_SERIES(-5, 4) x
        CROSS JOIN
                GENERATE_SERIES(-5, -4) y
        CROSS JOIN LATERAL
                (
                SELECT  ||((x, y, 100)::VEC3) AS pixel_unit
                ) pixel_unit
        ),
        sphere AS
        (
        SELECT 120 AS radius, (0, 0, 300)::VEC3 AS center
        ),
        spikes AS
        (
        SELECT  i, spherical.*
        FROM    (
                SELECT  40 AS spikes
                ) constants
        CROSS JOIN LATERAL
                GENERATE_SERIES(0, spikes - 1) AS i
        CROSS JOIN LATERAL
                (
                SELECT  ACOS(2 * i / spikes::DOUBLE PRECISION - 1) AS theta,
                        PI() * (3 - SQRT(5)) * i AS phi
                ) spherical
        ),
        spike_parameters AS
        (
        SELECT  0.3 AS height, 0.1 AS width
        ),
        spike_origin_vertices AS
        (
        SELECT  1 AS j, (0, 0, 1)::VEC3 AS origin_vertex
        UNION ALL
        SELECT  i + 2 AS j, (0, 0, 1 + height)::VEC3 + (SIN(2 * PI() * i / 3), COS(2 * PI() * i / 3), 0)::VEC3 * width AS origin_vertex
        FROM    spike_parameters
        CROSS JOIN
                GENERATE_SERIES(0, 2) i
        ),
        spike_vertices AS
        (
        SELECT  i, j, vertex
        FROM    spikes
        CROSS JOIN
                spike_origin_vertices
        CROSS JOIN LATERAL
                (
                SELECT  ((origin_vertex).x * COS(theta) + (origin_vertex).z * SIN(theta), (origin_vertex).y, (origin_vertex).z * COS(theta) - (origin_vertex).x * SIN(theta))::VEC3 AS rotated_y
                ) y
        CROSS JOIN LATERAL
                (
                SELECT  ((rotated_y).x * COS(phi) - (rotated_y).y * SIN(phi), (rotated_y).y * COS(phi) + (rotated_y).x * SIN(phi), (rotated_y).z)::VEC3 AS rotated
                ) rotated
        CROSS JOIN LATERAL
                (
                SELECT  (sphere.radius * rotated) + sphere.center AS vertex
                FROM    sphere
                ) mapped
        ),
        spike_triangles AS
        (
        SELECT  i,
                sv1.j + sv2.j + sv3.j - 5 AS j,
                sv1.vertex AS one, sv2.vertex AS two, sv3.vertex AS three
        FROM    spike_vertices sv1
        JOIN    spike_vertices sv2
        USING   (i)
        JOIN    spike_vertices sv3
        USING   (i)
        WHERE   sv2.j > sv1.j
                AND sv3.j > sv2.j
        ),
        spike_intersections AS
        (
        SELECT  i, j, one, two, three, x, y, pixel_unit, intersection_distance
        FROM    pixels
        CROSS JOIN
                origin
        CROSS JOIN
                spike_triangles
        CROSS JOIN LATERAL
                (
                SELECT  two - one AS edge1, three - one AS edge2
                ) edges
        CROSS JOIN LATERAL
                (
                SELECT  pixel_unit ** edge2 AS p_vector
                ) normal
        CROSS JOIN LATERAL
                (
                SELECT  edge1 * p_vector AS determinant
                ) det
        CROSS JOIN LATERAL
                (
                SELECT  origin - one AS t_vector
                ) t
        CROSS JOIN LATERAL
                (
                SELECT  (t_vector * p_vector) / determinant AS u
                WHERE   determinant NOT BETWEEN -1e-8 AND 1e-8
                ) u
        CROSS JOIN LATERAL
                (
                SELECT  t_vector ** edge1 AS q_vector
                WHERE   u BETWEEN 0 AND 1
                ) q
        CROSS JOIN LATERAL
                (
                SELECT  (pixel_unit * q_vector) / determinant AS v
                ) v
        CROSS JOIN LATERAL
                (
                SELECT  (edge2 * q_vector) / determinant AS intersection_distance
                WHERE   v >= 0 AND u + v <= 1
                ) intersection
        ),
        sphere_intersection_coefficients AS
        (
        SELECT  *
        FROM    pixels
        CROSS JOIN
                sphere
        CROSS JOIN
                origin
        CROSS JOIN LATERAL
                (
                SELECT  1 AS a, 2 * pixel_unit * (origin - center) AS b, (|(origin - center))^2 - radius^2 AS c
                ) q3
        CROSS JOIN LATERAL
                (
                SELECT b ^ 2 - 4 * a * c AS discriminant
                ) q4
        WHERE   discriminant > 0
        ),
        sphere_intersections AS
        (
        SELECT  x, y, pixel_unit, intersection_distance
        FROM    sphere_intersection_coefficients
        CROSS JOIN LATERAL
                (
                SELECT  t AS intersection_distance
                FROM    (VALUES (1), (-1)) q (sign)
                CROSS JOIN LATERAL
                        (
                        SELECT  (-b + SQRT(discriminant) * sign) / (2 * a) AS t
                        WHERE   discriminant > 0
                        ) q2
                WHERE   t > 0
                ORDER BY
                        t
                LIMIT   1
                ) q
        ),
        closest_intersections AS
        (
        SELECT  object, i, j, x, y, pixel_unit, intersection_distance
        FROM    (
                SELECT  q.*,
                        ROW_NUMBER() OVER (PARTITION BY x, y ORDER BY intersection_distance) rn
                FROM    (
                        SELECT  'sphere' AS object, 0 AS i, 0 AS j, x, y, pixel_unit, intersection_distance
                        FROM    sphere_intersections
                        UNION ALL
                        SELECT  'spike', i, j, x, y, pixel_unit, intersection_distance
                        FROM    spike_intersections
                        ) q
                ) q
        WHERE   rn = 1
        )
SELECT  *
FROM    closest_intersections
ORDER BY
        x, y

We see that some pixels show the intersection with the sphere, and some with the triangle 4 of the spike 0. When the intersection is with the spike, the distance is lower: the spike protrudes out of the sphere and the pixel ray hits it earlier than it would the sphere.

Reflections

To calculate the pixel intensity at each intersection we will be using the Phong reflection model.

This is a simple, yet quite realistic model which accounts for three kinds of reflection, each one adding to another.

• The ambient reflection is the same for all intersections. It's either there, if there is an intersection, or not. It simulates the scattered light from the walls, ceiling, other objects and stuff like this, which is always there and which is dim and uniform enough to bother with its exact calculation.
• The diffuse reflection, or Lambertian reflection, accounts for the light from the main source scattered, or diffusely reflected by the surface. By the Lambert's cosine law, the amount of scattered light which hits the pixel does not depend on the view point, but it does depend on the angle between the object's surface normal at this point and the source of light.
• The specular reflection is responsible for the "shiny" part of the reflection. Dr. Bui Tuong Phong, the author of the model, had noticed that some glossy surfaces, like marble or still water, behave somewhat in between the rough surfaces and perfect mirrors. They don't reflect all the light at the angle equal to the angle of incidence, but the closer the angle of reflection to the angle of incidence is, the more light they do reflect. A reflection of a lamp on a bowling ball looks larger and smoother than on a polished metal ball.

  The Phong model treats this component of reflected light as proportional to (R̅ ⋅ V̅)α, where R̅ is the reflection unit vector, V̅ is the pixel unit vector and α is the shininess constant. The shinier the object, the larger is this constant, the more sharp is the peak. For a perfect mirror, this component is infinite, all the light is reflected strictly at the angle equal to the angle of incidence, and it looks like a shiny dot on the object. For an object which is merely glossy, the reflection of the light source looks like a blurred spot

Unlike the previous term, this one does not depend explicitly on the surface normal and the direction to the light source, but it does depend implicitly on them, because R̅ depends on both these terms.

The formula for the intensity of light at each point is A + L * Il * (L̅ ⋅ N̅) + S * Is * (R̅ ⋅ V̅)α. A, L and S are reflection constants, I_l and I_s are light intensities for diffuse and specular reflection. They are different because the prevailing kind of reflection can depend on the light's color and polarization.

So, to calculate the intensity of reflected light at each pixel, we need to know the light intensities and the surface parameters (ambient reflection constant, Lambertian reflection constant, specular reflection constant and shininess). We also need to know the direction between the pixel intersection point and the light source. Finally, we need to know the direction of the vector normal to the surface at this point.

The surface parameters we will specify next to the rest of fields in the sphere and spike CTEs. For the light source, we will create another single-record CTE with the parameters (position and intensities).

As for the normal vectors in the intersection points, it makes sense to calculate them at the same time we calculate the intersections. They need to be calculated differently for the sphere and for the triangles.

Strictly speaking, there are two normal vectors to each surface point, differing in sign (because the surface has two sides). However, we can only see one side at each time. This means that we will have to select one normal vector from two, namely the one which gives the positive dot product with the point of view.

This is a problem for the triangles, but this is not a problem for the sphere because of its geometry. We will always have to select the normal vector pointing outside the sphere. A vector pointing inside the sphere would never give a positive dot product with the point of view, because this intersection would be the farthest one of the two, and would be blocked by another one. After all, we can never see the inside of a sphere if we are located outside!

For the sphere, the normal vector is just the unit vector from the center of the sphere to the intersection point.

For a triangle, the normal vector is a unit vector for a cross product of any two of its sides. The sign will have to be selected as described above.

The reflection unit vector can be calculated as 2 * (L̅ ⋅ N̅) * N̅ - L̅.

Let's calculate all these parameters and substitute them into Phong's formula. As before, we'll be only doing this for a small section of the grid.

WITH    origin AS
        (
        SELECT  (0, 0, 0)::VEC3 AS origin
        ),
        light AS
        (
        SELECT  (-5000, 5000, -5000)::VEC3 AS position, 0.7 AS specular, 0.6 AS lambertian
        ),
        pixels AS
        (
        SELECT  *
        FROM    GENERATE_SERIES(-5, 4) x
        CROSS JOIN
                GENERATE_SERIES(-5, -4) y
        CROSS JOIN LATERAL
                (
                SELECT  ||((x, y, 100)::VEC3) AS pixel_unit
                ) pixel_unit
        ),
        sphere AS
        (
        SELECT 120 AS radius, (0, 0, 300)::VEC3 AS center, 0.7 AS specular, 0.6 AS lambertian, 1 AS shininess
        ),
        spikes AS
        (
        SELECT  i, spherical.*
        FROM    (
                SELECT  40 AS spikes
                ) constants
        CROSS JOIN LATERAL
                GENERATE_SERIES(0, spikes - 1) AS i
        CROSS JOIN LATERAL
                (
                SELECT  ACOS(2 * i / spikes::DOUBLE PRECISION - 1) AS theta,
                        PI() * (3 - SQRT(5)) * i AS phi
                ) spherical
        ),
        spike_parameters AS
        (
        SELECT  0.3 AS height, 0.1 AS width, 2 AS specular, 4 AS lambertian, 4 AS shininess
        ),
        spike_origin_vertices AS
        (
        SELECT  1 AS j, (0, 0, 1)::VEC3 AS origin_vertex
        UNION ALL
        SELECT  i + 2 AS j, (0, 0, 1 + height)::VEC3 + (SIN(2 * PI() * i / 3), COS(2 * PI() * i / 3), 0)::VEC3 * width AS origin_vertex
        FROM    spike_parameters
        CROSS JOIN
                GENERATE_SERIES(0, 2) i
        ),
        spike_vertices AS
        (
        SELECT  i, j, vertex
        FROM    spikes
        CROSS JOIN
                spike_origin_vertices
        CROSS JOIN LATERAL
                (
                SELECT  ((origin_vertex).x * COS(theta) + (origin_vertex).z * SIN(theta), (origin_vertex).y, (origin_vertex).z * COS(theta) - (origin_vertex).x * SIN(theta))::VEC3 AS rotated_y
                ) y
        CROSS JOIN LATERAL
                (
                SELECT  ((rotated_y).x * COS(phi) - (rotated_y).y * SIN(phi), (rotated_y).y * COS(phi) + (rotated_y).x * SIN(phi), (rotated_y).z)::VEC3 AS rotated
                ) rotated
        CROSS JOIN LATERAL
                (
                SELECT  (sphere.radius * rotated) + sphere.center AS vertex
                FROM    sphere
                ) mapped
        ),
        spike_triangles AS
        (
        SELECT  i,
                sv1.j + sv2.j + sv3.j - 5 AS j,
                sv1.vertex AS one, sv2.vertex AS two, sv3.vertex AS three
        FROM    spike_vertices sv1
        JOIN    spike_vertices sv2
        USING   (i)
        JOIN    spike_vertices sv3
        USING   (i)
        WHERE   sv2.j > sv1.j
                AND sv3.j > sv2.j
        ),
        spike_intersections AS
        (
        SELECT  i, j, x, y, pixel_unit, intersection_distance, intersection, intersection_norm_unit,
                specular, lambertian, shininess
        FROM    pixels
        CROSS JOIN
                origin
        CROSS JOIN
                spike_triangles
        CROSS JOIN
                spike_parameters
        CROSS JOIN LATERAL
                (
                SELECT  two - one AS edge1, three - one AS edge2
                ) edges
        CROSS JOIN LATERAL
                (
                SELECT  pixel_unit ** edge2 AS p_vector
                ) normal
        CROSS JOIN LATERAL
                (
                SELECT  edge1 * p_vector AS determinant
                ) det
        CROSS JOIN LATERAL
                (
                SELECT  origin - one AS t_vector
                ) t
        CROSS JOIN LATERAL
                (
                SELECT  (t_vector * p_vector) / determinant AS u
                WHERE   determinant NOT BETWEEN -1e-8 AND 1e-8
                ) u
        CROSS JOIN LATERAL
                (
                SELECT  t_vector ** edge1 AS q_vector
                WHERE   u BETWEEN 0 AND 1
                ) q
        CROSS JOIN LATERAL
                (
                SELECT  (pixel_unit * q_vector) / determinant AS v
                ) v
        CROSS JOIN LATERAL
                (
                SELECT  (edge2 * q_vector) / determinant AS intersection_distance
                WHERE   v >= 0 AND u + v <= 1
                ) distance
        CROSS JOIN LATERAL
                (
                SELECT  origin + (pixel_unit * intersection_distance) AS intersection
                FROM    origin
                ) intersection
        CROSS JOIN LATERAL
                (
                SELECT  ||(edge1 ** edge2) AS intersection_norm_unit_oriented
                ) norm_oriented
        CROSS JOIN LATERAL
                (
                SELECT  intersection_norm_unit_oriented * -SIGN(pixel_unit * intersection_norm_unit_oriented) AS intersection_norm_unit
                ) norm
        ),
        sphere_intersection_coefficients AS
        (
        SELECT  *
        FROM    pixels
        CROSS JOIN
                sphere
        CROSS JOIN
                origin
        CROSS JOIN LATERAL
                (
                SELECT  1 AS a, 2 * pixel_unit * (origin - center) AS b, (|(origin - center))^2 - radius^2 AS c
                ) q3
        CROSS JOIN LATERAL
                (
                SELECT b ^ 2 - 4 * a * c AS discriminant
                ) q4
        WHERE   discriminant > 0
        ),
        sphere_intersections AS
        (
        SELECT  x, y, pixel_unit, intersection_distance, intersection, intersection_norm_unit,
                specular, lambertian, shininess
        FROM    sphere_intersection_coefficients
        CROSS JOIN LATERAL
                (
                SELECT  t AS intersection_distance
                FROM    (VALUES (1), (-1)) q (sign)
                CROSS JOIN LATERAL
                        (
                        SELECT  (-b + SQRT(discriminant) * sign) / (2 * a) AS t
                        WHERE   discriminant > 0
                        ) q2
                WHERE   t > 0
                ORDER BY
                        t
                LIMIT   1
                ) q
        CROSS JOIN LATERAL
                (
                SELECT  origin + (pixel_unit * intersection_distance) AS intersection
                FROM    origin
                ) intersection
        CROSS JOIN LATERAL
                (
                SELECT  ||(intersection - center) AS intersection_norm_unit
                FROM    sphere
                ) q3
        ),
        closest_intersections AS
        (
        SELECT  *
        FROM    (
                SELECT  q.*,
                        ROW_NUMBER() OVER (PARTITION BY x, y ORDER BY intersection_distance) rn
                FROM    (
                        SELECT  'sphere' AS object, 0 AS i, 0 AS j, *
                        FROM    sphere_intersections
                        UNION ALL
                        SELECT  'spike', *
                        FROM    spike_intersections
                        ) q
                ) q
        WHERE   rn = 1
        ),
        intensities AS
        (
        SELECT  object, i, j, x, y, intersection, intersection_norm_unit, intensity
        FROM    closest_intersections
        CROSS JOIN
                light
        CROSS JOIN LATERAL
                (
                SELECT  ||(light.position - intersection) AS light_unit
                ) q
        CROSS JOIN LATERAL
                (
                SELECT  2 * (light_unit * intersection_norm_unit) * intersection_norm_unit - light_unit AS reflection_unit
                ) q2
        CROSS JOIN LATERAL
                (
                SELECT  GREATEST(light_unit * intersection_norm_unit, 0) AS norm_term
                ) q3
        CROSS JOIN LATERAL
                (
                SELECT  0.125 +
                        closest_intersections.lambertian * light.lambertian * norm_term +
                        closest_intersections.specular * light.specular * CASE WHEN norm_term > 0 THEN GREATEST(reflection_unit * -pixel_unit, 0) ^ shininess ELSE 0 END +
                        0 AS intensity
                ) q4
        )
SELECT  *
FROM    intensities
ORDER BY
        x, y

We can see that the normal vector to the spike triangle is the same for all the points. This is expected for a flat surface. The intensity values differ ever so slightly for the points on the spike triangle. Also, the spike is much more shiny than the sphere, because of the way the constants are picked (for aesthetic reasons).

The picture

Now that we have all the pieces together, let's build the picture of the dreaded coronavirus.

Note that this query takes quite a time (about 3 minutes on my machine), but hey, that's ray tracing in SQL! As Dr. Kajiya had put it in one of his papers, "ray tracing isn't too slow; computers are too slow".

Here's the query:

WITH    origin AS
        (
        SELECT  (0, 0, 0)::VEC3 AS origin
        ),
        light AS
        (
        SELECT  (-5000, 5000, -5000)::VEC3 AS position, 0.7 AS specular, 0.6 AS lambertian
        ),
        pixels AS
        (
        SELECT  *
        FROM    GENERATE_SERIES(-70, 70) x
        CROSS JOIN
                GENERATE_SERIES(-70, 70) y
        CROSS JOIN LATERAL
                (
                SELECT  ||((x, y, 100)::VEC3) AS pixel_unit
                ) pixel_unit
        ),
        sphere AS
        (
        SELECT 120 AS radius, (0, 0, 300)::VEC3 AS center, 0.7 AS specular, 0.6 AS lambertian, 1 AS shininess
        ),
        spikes AS
        (
        SELECT  i, spherical.*
        FROM    (
                SELECT  40 AS spikes
                ) constants
        CROSS JOIN LATERAL
                GENERATE_SERIES(0, spikes - 1) AS i
        CROSS JOIN LATERAL
                (
                SELECT  ACOS(2 * i / spikes::DOUBLE PRECISION - 1) AS theta,
                        PI() * (3 - SQRT(5)) * i AS phi
                ) spherical
        ),
        spike_parameters AS
        (
        SELECT  0.3 AS height, 0.1 AS width, 2 AS specular, 4 AS lambertian, 4 AS shininess
        ),
        spike_origin_vertices AS
        (
        SELECT  1 AS j, (0, 0, 1)::VEC3 AS origin_vertex
        UNION ALL
        SELECT  i + 2 AS j, (0, 0, 1 + height)::VEC3 + (SIN(2 * PI() * i / 3), COS(2 * PI() * i / 3), 0)::VEC3 * width AS origin_vertex
        FROM    spike_parameters
        CROSS JOIN
                GENERATE_SERIES(0, 2) i
        ),
        spike_vertices AS
        (
        SELECT  i, j, vertex
        FROM    spikes
        CROSS JOIN
                spike_origin_vertices
        CROSS JOIN LATERAL
                (
                SELECT  ((origin_vertex).x * COS(theta) + (origin_vertex).z * SIN(theta), (origin_vertex).y, (origin_vertex).z * COS(theta) - (origin_vertex).x * SIN(theta))::VEC3 AS rotated_y
                ) y
        CROSS JOIN LATERAL
                (
                SELECT  ((rotated_y).x * COS(phi) - (rotated_y).y * SIN(phi), (rotated_y).y * COS(phi) + (rotated_y).x * SIN(phi), (rotated_y).z)::VEC3 AS rotated
                ) rotated
        CROSS JOIN LATERAL
                (
                SELECT  (sphere.radius * rotated) + sphere.center AS vertex
                FROM    sphere
                ) mapped
        ),
        spike_triangles AS
        (
        SELECT  sv1.vertex AS one, sv2.vertex AS two, sv3.vertex AS three
        FROM    spike_vertices sv1
        JOIN    spike_vertices sv2
        USING   (i)
        JOIN    spike_vertices sv3
        USING   (i)
        WHERE   sv2.j > sv1.j
                AND sv3.j > sv2.j
        ),
        spike_intersections AS
        (
        SELECT  x, y, pixel_unit, intersection_distance, intersection, intersection_norm_unit,
                specular, lambertian, shininess
        FROM    pixels
        CROSS JOIN
                origin
        CROSS JOIN
                spike_triangles
        CROSS JOIN
                spike_parameters
        CROSS JOIN LATERAL
                (
                SELECT  two - one AS edge1, three - one AS edge2
                ) edges
        CROSS JOIN LATERAL
                (
                SELECT  pixel_unit ** edge2 AS p_vector
                ) normal
        CROSS JOIN LATERAL
                (
                SELECT  edge1 * p_vector AS determinant
                ) det
        CROSS JOIN LATERAL
                (
                SELECT  origin - one AS t_vector
                ) t
        CROSS JOIN LATERAL
                (
                SELECT  (t_vector * p_vector) / determinant AS u
                WHERE   determinant NOT BETWEEN -1e-8 AND 1e-8
                ) u
        CROSS JOIN LATERAL
                (
                SELECT  t_vector ** edge1 AS q_vector
                WHERE   u BETWEEN 0 AND 1
                ) q
        CROSS JOIN LATERAL
                (
                SELECT  (pixel_unit * q_vector) / determinant AS v
                ) v
        CROSS JOIN LATERAL
                (
                SELECT  (edge2 * q_vector) / determinant AS intersection_distance
                WHERE   v >= 0 AND u + v <= 1
                ) distance
        CROSS JOIN LATERAL
                (
                SELECT  origin + (pixel_unit * intersection_distance) AS intersection
                FROM    origin
                ) intersection
        CROSS JOIN LATERAL
                (
                SELECT  ||(edge1 ** edge2) AS intersection_norm_unit_oriented
                ) norm_oriented
        CROSS JOIN LATERAL
                (
                SELECT  intersection_norm_unit_oriented * -SIGN(pixel_unit * intersection_norm_unit_oriented) AS intersection_norm_unit
                ) norm
        ),
        sphere_intersection_coefficients AS
        (
        SELECT  *
        FROM    pixels
        CROSS JOIN
                sphere
        CROSS JOIN
                origin
        CROSS JOIN LATERAL
                (
                SELECT  1 AS a, 2 * pixel_unit * (origin - center) AS b, (|(origin - center))^2 - radius^2 AS c
                ) q3
        CROSS JOIN LATERAL
                (
                SELECT b ^ 2 - 4 * a * c AS discriminant
                ) q4
        WHERE   discriminant > 0
        ),
        sphere_intersections AS
        (
        SELECT  x, y, pixel_unit, intersection_distance, intersection, intersection_norm_unit,
                specular, lambertian, shininess
        FROM    sphere_intersection_coefficients
        CROSS JOIN LATERAL
                (
                SELECT  t AS intersection_distance
                FROM    (VALUES (1), (-1)) q (sign)
                CROSS JOIN LATERAL
                        (
                        SELECT  (-b + SQRT(discriminant) * sign) / (2 * a) AS t
                        WHERE   discriminant > 0
                        ) q2
                WHERE   t > 0
                ORDER BY
                        t
                LIMIT   1
                ) q
        CROSS JOIN LATERAL
                (
                SELECT  origin + (pixel_unit * intersection_distance) AS intersection
                FROM    origin
                ) intersection
        CROSS JOIN LATERAL
                (
                SELECT  ||(intersection - center) AS intersection_norm_unit
                FROM    sphere
                ) q3
        ),
        closest_intersections AS
        (
        SELECT  *
        FROM    (
                SELECT  q.*,
                        ROW_NUMBER() OVER (PARTITION BY x, y ORDER BY intersection_distance) rn
                FROM    (
                        SELECT  *
                        FROM    sphere_intersections
                        UNION ALL
                        SELECT  *
                        FROM    spike_intersections
                        ) q
                ) q
        WHERE   rn = 1
        ),
        intensities AS
        (
        SELECT  x, y, intensity
        FROM    closest_intersections
        CROSS JOIN
                light
        CROSS JOIN LATERAL
                (
                SELECT  ||(light.position - intersection) AS light_unit
                ) q
        CROSS JOIN LATERAL
                (
                SELECT  2 * (light_unit * intersection_norm_unit) * intersection_norm_unit - light_unit AS reflection_unit
                ) q2
        CROSS JOIN LATERAL
                (
                SELECT  GREATEST(light_unit * intersection_norm_unit, 0) AS norm_term
                ) q3
        CROSS JOIN LATERAL
                (
                SELECT  0.125 +
                        closest_intersections.lambertian * light.lambertian * norm_term +
                        closest_intersections.specular * light.specular * CASE WHEN norm_term > 0 THEN GREATEST(reflection_unit * -pixel_unit, 0) ^ shininess ELSE 0 END +
                        0 AS intensity
                ) q4
        )
SELECT  STRING_AGG(dot, '' ORDER BY x) AS picture
FROM    pixels p
LEFT JOIN
        intensities i
USING   (x, y)
CROSS JOIN LATERAL
        (
        SELECT  ' .:—=+*#%@' AS palette
        ) palette
CROSS JOIN LATERAL
        (
        SELECT  SUBSTRING(palette, LEAST(FLOOR(COALESCE(intensity, 0) * LENGTH(palette))::INT + 1, LENGTH(palette)), 1) AS dot
        ) dots
GROUP BY
        y
ORDER BY
        y

You can view the queries here: https://github.com/quassnoi/explain-extended-2021

In the New Year, I wish all of us to stay safe, get vaccinated and not let the filthy beast affect our lives anymore!

Previous New Year posts:

• Happy New 2010 Year!
• Happy New 2011 Year!
• Happy New 2012 Year!
• Happy New 2013 Year!
• Happy New 2014 Year!
• Happy New 2015 Year!
• Happy New 2016 Year!
• Happy New 2017 Year!
• Happy New 2018 Year!
• Happy New 2019 Year!
• Happy New 2020 Year!