Just turning the comments into a fully-fledged answer as it may prove useful to others.

IMHO, the alpha channel should always be linear and treated as such.

You described blending a black quad over a white background using alpha=0.95, and mapping the result back to non-linear space. We'd expect a resulting, physical intensity value of 0.05.

If the non-linear space is sRGB, then a linear value of "0.05" would correspond to an sRGB value of "0.248" which, when mapped to 8-bit, is around 63 (or hex, 3D)

If your display is reasonably calibrated, you can verify this by creating, say, a 5x4 black image, setting one pixel to white, and then tiling that to a suitable size, e.g 200x200 pixels. Then draw a, say, 50x50 square in the middle and fill it with 0x3D3D3D. Here's an example:

If you take a few steps back from your monitor and squint, the centre grey region and the dithered outer should be almost indistinguishable, demonstrating a correct physical model.

The problem you are describing appears to be related to the fact that our visual system is non-linear, much as our hearing is too (though probably not with the same curve). For example, the analogue volume control on a hifi amplifier will use a logarithmic potentiometer rather than a linear one .

You thus 'simply' need to change your per-frame alpha non-linearly to get a perceptually even fade out.


Why did tessellation come to be a prominent feature?

The main purpose of tesselation is to increase the resolution of the mesh, while only transferring a small amount of triangle data around. In addition, tessellation allows you to dynamically change the LOD of the mesh, so you can optimize your shader calls.

So, we can pass the GPU, say, 3000 triangles, and have it tesselate it to 300000 triangles. We are essentially trading storage space/bandwidth for compute power. Since GPUs have lots and lots of compute, and memory is limited/slow, this is a pretty good tradeoff.

As for "not being heavily used in real-time graphics". I somewhat disagree. Many AAA games have been using tessellation for a long time. That said, tessellation is hard to get right. Done wrong, tesselation can lead to lots of problems, such as cracking, or over tessellating to sub-pixel triangles, which destroys your fill rate.

Combined, these problems make it difficult to implement good and fast tessellation for your everyday person. Thus, most of the uses of tessellation you see nowadays are in AA or AAA games, game engines, and offline rendering tools.

That said, there is active research going on that is trying to better utilize the tessellation hardware, and make it easier to use. For example: Efficient GPU Rendering of Subdivision Surfaces using Adaptive Quadtrees

In the end, tessellation is a great feature, but is difficult to get right


PBR missing rim darkening on non metallic rough surfaces

The problem is in the diffuse term, which can be seen by making the specular portion of IBL not be added into the result.  The diffuse only render will not have the darkening, and of course, the specular can't make the diffuse darker. (Note, fresnel vs no fresnel is not the issue here).

The core problem is that lambertian diffuse is being used, which is view independent. There are other view independent diffuse terms which give results closer to the ground truth.

Two example diffuse terms are Burley (Disney) and Oren-Nayar.

Burley is talked about here: https://google.github.io/filament/Filament.md.html#toc4.5

An approach used in Titanfall is described here:
https://www.gdcvault.com/play/1024478/PBR-Diffuse-Lighting-for-GGX


Particle & Fluid Simulation: Where to learn?

PBRT handles only the rendering part, not the simulation part. There are several books out on game physics that handle fluid simulation, the one on my bookshelf is David Eberly's "Game Physics".


What is the accepted method of converting shininess to roughness, and vice versa?

As you already note, there is no clear cut interpretation/conversion for these values. I think it is even much worse: Depending on your BRDF and internal limitations (like having defined exponents ranging from 2-2048) the interpretation is completely different. Like suggested in the comments, it might be the best to render a series with different values and fit a conversion curve until the value looks intuitive.

A few examples I was able to find some blog posts that mention something about that topic:

At Dontnod entertainment they use a "perceptual linear distribution".
Sébastien Lagarde acknowledges the problem in this blog post and writes a few notes on that.
Brian Karis uses a squared roughness values in this Microfacet BRDF overview. This illustrates also nicely how differently  roughness is used in different Normal Distribution Functions. Blinn-Phong power is here defined as 2/roughness^4 - 2.
Frostbite uses a squared remapping. ie. Roughness =  (1 − Smoothness)^2
Details about it and their entire material system is explained in section 3.2 of Sebastien Lagarde's writeup.
This blog post suggests to define the roughness for a Beckmann distribution from the roughness alpha as:


Composition of decomposed transformations

Yes you can, provided that all the nodes' scales are only uniform scaling—the same scaling factor along all three axes. If there is nonuniform scaling (different factors along different axes) then things get more complicated, and the composed transformation generally can't be decomposed cleanly into scale, rotation, and translation.

One way to accomplish this is to convert all the individual nodes' transforms into 4×4 matrix form, multiply the matrices down the tree from root to leaf, then decompose the resulting matrix (see this answer for how to do the decomposition).

However, if you have all the transforms in TRS form (translate, rotate, scale—assuming column vector convention) and you want the final transform in TRS form, you can also combine them directly.

The final scale is just the product of all the individual nodes' scales
The final rotation is just the composition of all the individual nodes' rotations (in quaternion form for example)
The final translation is made by summing the translation vectors, but transforming each vector by the composite scale and rotation of all the nodes above it in the tree. For example if you have a transformation $T_1 R_1 S_1 T_2 R_2 S_2 T_3 R_3 S_3$ (numbering from the root down to the leaf), then the final translation is $t_1 + (R_1 S_1) t_2 + (R_1 S_1 R_2 S_2) t_3$.


Why are RGB tertiary colors so similar to RGB primary colors?

Likely a gamma correction issue. Tertiary colours are some permutation of (1, 0.5, 0) in RGB, right? So the issue arises because without applying gamma, 0.5 appears less than half as bright as 1, so the blended colour skews towards the primary (1,0,0) rather than the secondary (1,1,0). The image below has linear colour on top, gamma corrected on the bottom, you can see the gamma corrected colours appear much closer to a halfway tone.


The usual way of doing this is to work in a color space that has luminance as one of the color components. For example the Y'CbCr color space has a luminance channel (you can think of this as similar to brightness which is a poorly defined term) and 2 chrominance or color channels. There are other color spaces that separate colors in this way, too, such as HSV, HLS, Lab, etc.

There are conversions to go between RGB and these various color spaces. The article linked above shows some of the conversions for going between Y'CbCr and sRGB, for example.

You don't say what type of graphics library you'll be using, so it's hard to give any more detail than what I've given you here. But I've done this in OpenGL, Core Graphics, and others. It shouldn't be too difficult.

Edited to add: I should point out that HSV defines value or "brightness" not as luminance but as the average of red, green, and blue, and as such does not track very will with our intuitive sense of brightness. HLS also suffers from some issues (discontinuities, bad but popular reference implementations) and is generally not recommended by the color scientists that I know.


$E_n$ is illuminance (or irradiance), $L_i$ is luminance (or radiance) in direction $\omega_i$ and $d\omega$ is differential solid angle.

TL;DR: $d\omega_i$ is the area of a pixel in the cubemap projected onto the unit sphere.

And then a longer explanation (:

Solid angle is essentially the projected area of an object on a unit sphere and has unit of steradians (sr). Entire unit sphere has surface area of $4\pi$ and hemisphere has the area of $2\pi$, which comes up a lot in lighting calculations due to integrals over spheres and hemispheres respectively.

If you would want to calculate solid angle for an object and point $p$, you would project the object on the unit sphere at $p$ and calculate the area of the projection. Further the object is from $p$ smaller the solid angle gets. As an example, we can calculate solid angle for Sun as seen from the Earth. We know that Sun has about 0.52 degrees subtended angle from Earth, so we can calculate that the solid angle is:
$$2\pi*(1-Cos(0.52/2))\simeq0.00006469sr$$

So the differential solid angle is just an infinitesimally small solid angle. If you would numerically integrate over a sphere ($4\pi$ steradians) by taking $N$ uniform samples on the sphere, your $d\omega$ would be $4\pi/N$. With integrals $N$ goes to infinity, while with numerical approximation of the integral $N$ is some "sufficiently large" number.

However, in this paper you refer to they perform integration by rendering to a cubemap at $p$ and then iterating through all the rendered pixels. But because these pixels are rendered onto planes (cubemap faces), their projection on the unit sphere isn't constant, so $d\omega$ needs to be weighted by the position in the rendered image, which they explain in your second quote:

  Computing the solid angle dωi subtended by a particular texel is a
  little more involved. Texels near the center of the texture subtend a
  somewhat larger angle, while texels near the corner take up a smaller
  solid angle when we map the cube texture onto a sphere.


Fresnel equation with incident angle, n and k as input

@PaulHK's answer is correct I'm sure, here's a bit of a check to show that the IOR() function is calculating the reflection coefficients for $s$ and $p$ polarizations then averaging the two assuming unpolarized incident light.

The reflection is for a single interface, and at least at normal incidence the results reduce to the simpler $(n_2-n_1)^2/(n_2+n_1)^2$. With this script you can check against other on-line calculators using non-zero incident angle and nonzero values for $k$.

I cleaned up the python a bit (indents) and added the plot. You can see when exiting from glass to air, the reflectivity goes to 100% at the critical angle for total internal reflection:

def IOR(n,k):
    theta_deg = 0

    fresnel = []

    while theta_deg <= 90:
        theta = math.radians(theta_deg)

        a = math.sqrt((math.sqrt((n**2-k**2-(math.sin(theta))**2)**2 +
                                 ((4 * n**2) * k**2)) + (n**2 - k**2 -
                                            (math.sin(theta))**2))/2)

        b = math.sqrt((math.sqrt((n**2-k**2-(math.sin(theta))**2)**2 +
                                 ((4 * n**2) * k**2)) - (n**2 - k**2 -
                                            (math.sin(theta))**2))/2)

        Fs = (a**2+b**2-(2 * a * math.cos(theta))+
              (math.cos(theta))**2)/(a**2+b**2 +
                            (2 * a * math.cos(theta))+(math.cos(theta))**2)

        Fp = Fs * ((a**2+b**2 -
                    (2 * a * math.sin(theta) * math.tan(theta)) +
                    (math.sin(theta))**2*(math.tan(theta))**2)/(a**2+b**2 +
                    (2 * a * math.sin(theta) * math.tan(theta)) +
                    (math.sin(theta))**2*(math.tan(theta))**2))

        R = (Fs + Fp)/2

        fresnel.append((R, Fs, Fp))

        theta_deg += 1
    return fresnel

import math
import matplotlib.pyplot as plt

n1, n2  = 1.00, 1.45  # into glass, one side only
n2ovrn1 = n2/n1

Rin, Fsin, Fpin = zip(*IOR(n2ovrn1, 0))

print Rin[0]
print (n2ovrn1-1.0)**2 / (n2ovrn1+1.0)**2  # check

n1, n2  = 1.45, 1.00  # outof glass, one side only
n2ovrn1 = n2/n1

Rout, Fsout, Fpout = zip(*IOR(n2ovrn1, 0))

print Rout[0]
print (n2ovrn1-1.0)**2 / (n2ovrn1+1.0)**2  # check

plt.figure()
plt.plot(range(91), Rin, '-k', linewidth=2)
plt.plot(range(91), Rout, '-k', linewidth=2)
plt.plot(range(91), Fsin, '-k')
plt.plot(range(91), Fsout, '-k')
plt.plot(range(91), Fpin, '--k')
plt.plot(range(91), Fpout, '--k')
plt.show()


Best way to clear a framebuffer image without drawing?

You are looking for vkCmdClearColorImage,

It does exactly what it says on the tin, it clears an image to a specified color. This command uses the Transfer stage as far as barriers are concerned. You will need to add one between the clear and the vkCmdDispatch.

There is also a vkCmdClearDepthStencilImage if you need it.


Ray tracing implicit surfaces using Sturm's theorem

This is a rather late reply but it might interest others who come up to this question.
Your question assumes that your ray tracing platform does not support double precision. This is unlikely.
If you are using CUDA or OpenCL they both support double precision vectors.
If you are using Vulkan things are a little more complicated due to specification saying: "The precision of double-precision instructions is at least that of single precision."
Though not a ray tracing platform OpenGL also supports double precision vectors from OpenGL 4.0 and on.
If you are tracing on CPU, chances are you can define a double precision vector data type, if it is not already defined.

If, in the unlikely case that your ray tracing platform does not support double precision, or you are obliged to use single precision in any way, the only thing that I can think of is using pseudo-remainder sequences for rem- i.e. euclidean division.
Particularly you might want to check out Subresultant pseudo-remainder sequence from the wikipedia section in the link. I am adding the pseudocode of the algorithm in here:
r_0 = a
r_1 = b
for (i:= 1;r_i not equal to 0; i := i+1)
    d_i := degree_of(r_{i-1}) - degree_of(r_i)
    g_i := leading_coefficient_of(r_i)
    if i = 1 then
        b_1 := (-1)^{d_1 + 1}
        p_1 := -1
    else
        b_i := -g_{i-1} (p_i)^d_i
        p_i := ((-g_{i-1})^d_{i-1}) / (p_{i-1}^{d_{i-1}-1}
    end if
    r_i := rem((g_i)^{d_i+1} r_{i-1}, r_i) / b_i
end for


Is the division of the coordinates x/w,y/w,z/w the last operation that happens on the GPU?

The division you are referring to is generally called "The Perspective Divide". You can google it and get information overload on helpful documentation. The division itself happens after the geometry shader before the fragment shader as part of rasterization. And is part of the process that creates the fragments that the pixel shader is run on. It is also an important part of vertex attribute interpolation. All of OpenGL, Vulkan, DirectX perform the perspective divide in the same place which is done by the GPU's hardware.
The most adopted floating point standard is IEEE 754 and yes virtually all vendors implement and follow these standards. I highly recommend (slowly) wrapping your mind around this standard to get a deeper understanding of how floats are processed. Here is a good first factoid to memorize (if you want): The standard put the sign bit for floating point numbers in the same place as the sign bit for integers to make conversion between and understanding both easier.
Quantization in this context refers to the accuracy of floating point numbers. There are many numbers that can not be fully represented in floating point math on the GPU, such as 1/3. Quantization is the process of representing that number as closely as possible within the limitations of the hardware. But minor difference in hardware implementations that are still standards compliant can and often will produce slightly different results.
Rather then focus on the limits of floating point numbers or attempting to write code that operates at the edge of these limits, coders generally use techniques that avoid the issue altogether like utilizing hardware offset capabilities to side step the issue completely is one method. (but is situation specific, like when applying decals)


How to describe a position of a point w.r.t the position and orientation of 3 other points

FYI, the 3 other points lie on a plane

Of course they do. Any set of three points defines a plane (except in the degenerate case joojaa mentions, where they lie on many planes). If the points are $p$, $q$, and $r$, then $(q - p) \wedge (r - p)$ is the normal to the plane; call it $n$.

This plane defines a co-ordinate system. If you make it the XY plane, then $z = a \dot\ n$. You can also project $a$ onto the plane to get the other two co-ordinates: $x = (a - p) \dot\ (q - p)$ and $y = (a - p) \dot\ (r - p)$.

Note that the co-ordinate system I've described here is pretty rubbish: it isn't even orthogonal. It's good enough to uniquely identify a point, as long as you're careful about all the signs.

But if you want an orthonormal basis, you can construct it from this one. First, instead of $q - p$ and $r - p$, normalize those vectors. Then $n$ will also have length 1, and the $x$ axis will also be normal. The last problem is the $y$ axis, which is not orthogonal to the $x$ axis. You can get a better $y$ axis by crossing the other two and normalizing, and while you're doing that you might as well also choose the order to ensure the basis is right-handed. Your three basis vectors are as follows:

$$
\begin{eqnarray}
\hat{z} & = & \frac{n}{|n|} \\
\hat{x} & = & \frac{q - p}{|q - p|} \\
\hat{y} & = & \hat{z} \wedge \hat{x}
\end{eqnarray}
$$


I personally like glad. The have an online generator as well http://glad.dav1d.de

It will generate exactly what you need for your platform.


Calculating Material/Texture Placement from Bitmap to 3D Model

For each visible voxel side you need to iterate over all triangles and find if ray from voxel side normal intersects that particular triangle. Having intersection point you can then calculate barycentric coordinates of this point within intersected triangle. And with barycentric coordinates you can calculate UV of intersection point and therefore sample color from original texture. If you need more than one pixel per voxel side then you can cast more rays.


Screenspace Normals - Creation, Normal Maps, and Unpacking

How are screenspace normals created, and is this step before or after using normal maps or bump maps?

They are created after using normal maps. In deferred rendering, you write to the various buffers (diffuse, normal, depth, etc.) in fragment shaders. By this time, the normal maps would have already been read.

  If it's done before using normal maps, how are normal maps going to be affected by the screenspace-ness of the normals,...

If you are referring to normal maps on models, normal maps can stay the same. You can compress normal maps as well if you want, but this becomes irrelevant to what you do for writing to the normal buffer.

  ...and if after, how can you justify using a Model-View matrix on all fragments rather on vertices? Isn't that a lot more calculations?

Your model's normals are computed per vertex and interpolated through rasterization to be per fragment. The normal maps are processed per fragment, but these are just texture lookups, and these can be applied to the camera-space normals produced from the previous pipeline steps. Therefore, a view matrix is no longer necessary.

  Finally, how are they unpacked? I realize you can get the blue component using pythagoras, but how are they returned to world space?

The idea behind 2-component normals versus 3-component normals is that the vector magnitude isn't often taken into account for normals. Also, normals facing away from the camera generally aren't useful. With two components, you can use model a hemisphere facing the camera. You unpack them the same way you pack them.

Think about normals as positions. When converting from world space (x, y, z) to screenspace (x, y, depth), but the depth does not effect where on the screen a mesh gets displayed. Convert the 3D vector to a vector in terms of screenspace (through view matrix multiplications).

Then you can discard the z-component (aligned with the camera) as long as you normals are unit vectors (because you can derive z from the other two vectors).


How to modify Perlin (not simplex) noise to create continental like terrain generation

Perlin noise is just a base block, not very interesting by itself. You don't need to modify it, but to combine and filter it in interesting ways. Look at how to make fractal Brownian motion (fBm) with it for example, which combines octaves based on few parameters to get a richer texture.

The question of terrain rendering is a difficult one and a topic of active research. Many generated terrains have this artificial and uninteresting look. This is exactly what the paper "Modeling Real-World Terrain with Exponentially Distributed Noise" tries to address, by proposing a formula based on statistics measured from actual terrain, which yields more convincing shapes.

Another source of valuable information is Iñigo Quilez's website, where he mentions procedural generation a lot, including tricks to fake erosion.


Is there a general rule with regard to using the inverse of the "camera" to world when using a homogeneous transform matrix?

Is there a general rule with regard to using the inverse of the “camera” to world when using a homogeneous transform matrix?

No, there is no general rule. You should always think about which approach is the best one for your specific use-case.

Nothing else is ever done with the local-to-world matrix - other than taking its inverse. So I'm wondering if it would be easier just to define the world-to-local matrix directly?

Don't forget, that the source you linked is a tutorial. It has a limited scope and avoids confusing the reader with additional information he might not understand at the time of reading. That said, if I have an isolated look at the rendering process, I can't think of a reason why you absolutely need the cam-to-world matrix there because you want to transform everything into camera space before continuing with other transformations (projection, clipping, etc.). Even though it seems to make no sense to transform something out of camera space just to transform it back later, there still might be some reasons to do exactly this. For example, if some renderable object is attached to the camera coordinate system. Using an extra shader or a branch statement just for objects defined in the camera system usually costs you much more performance than two obsolete transformations.
If we now consider other operations, not just rendering, you might find some useful applications of the camera to world matrix, for example, the movement of a camera that always moves into view direction or orthogonal to it (free camera mode in an FPS). The target point of the movement can be easily determined in camera space. It is just a vector of the movement speeds in each principle direction. Use the camera to world matrix to transform the target point into world space and you have your new location (if no collisions occur).
Lastly, another reason to take the "inverse-route" might be that it is just easier, less work, or more convenient. All objects, including the camera, usually have a world space coordinate and orientation. The transformation to world space is always the same series of operations. Scale-rotate-translate. Since it so common and frequently used, you probably have an optimized function that does not multiply 3 individual matrices but calculates all the entries of the resulting matrix directly. You can still do the same for the world-to-local transformation, but it is less frequently used and if you only need it for your world-to-camera transformation it might be just easier to calculate the inverse since every linear algebra library you might use should provide this function. If you do this not too frequently each frame, it shouldn't affect performance in any notable way.
As I said in the beginning, you should always focus on what your program needs. If you have no use for the camera-to-world matrix other than to calculate the world-to-camera matrix, you can definitely skip this and calculate the world-to-camera matrix directly. Just ask yourself, if it is worth it to add some extra lines of code for that.


How is orthographic projection used in computer graphics technically classified as a projection?

I think it is more of a terminology shared specifically in the computer graphics community. The transformation $P$ results in a vector one-step-ahead from getting the actual projection. From a 3D position transformed with $P$, we can obtain both

the 2D-projected position in normalized device coordinate (NDC) (by obtaining xy component)
the depth (by obtaining the z component) in NDC

This so-called projection matrix is designed specifically to satisfy the needs in the computer graphics pipeline. The key difference between an ordinary N-D-to-(N-1)-D projection is that the projection matrix results in 2D positions with "depth", which is required for obtaining a depth map used for many purposes in the computer graphics pipeline.


In order to understand Russian Roulette, let's look at a very basic backward path tracer:

void RenderPixel(uint x, uint y, UniformSampler *sampler) {
    Ray ray = m_scene->Camera.CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);

    // Bounce the ray around the scene
    for (uint bounces = 0; bounces < 10; ++bounces) {
        m_scene->Intersect(ray);

        // The ray missed. Return the background color
        if (ray.geomID == RTC_INVALID_GEOMETRY_ID) {
            color += throughput * float3(0.846f, 0.933f, 0.949f);
            break;
        }

        // We hit an object

        // Fetch the material
        Material *material = m_scene->GetMaterial(ray.geomID);
        // The object might be emissive. If so, it will have a corresponding light
        // Otherwise, GetLight will return nullptr
        Light *light = m_scene->GetLight(ray.geomID);

        // If we hit a light, add the emmisive light
        if (light != nullptr) {
            color += throughput * light->Le();
        }

        float3 normal = normalize(ray.Ng);
        float3 wo = normalize(-ray.dir);
        float3 surfacePos = ray.org + ray.dir * ray.tfar;

        // Get the new ray direction
        // Choose the direction based on the material
        float3 wi = material->Sample(wo, normal, sampler);
        float pdf = material->Pdf(wi, normal);

        // Accumulate the brdf attenuation
        throughput = throughput * material->Eval(wi, wo, normal) / pdf;

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.org = surfacePos;

        // Reset the other ray properties
        ray.dir = wi;
        ray.tnear = 0.001f;
        ray.tfar = embree::inf;
        ray.geomID = RTC_INVALID_GEOMETRY_ID;
        ray.primID = RTC_INVALID_GEOMETRY_ID;
        ray.instID = RTC_INVALID_GEOMETRY_ID;
        ray.mask = 0xFFFFFFFF;
        ray.time = 0.0f;
    }

    m_scene->Camera.FrameBuffer.SplatPixel(x, y, color);
}

IE. we bounce around the scene, accumulating color and light attenuation as we go. In order to be completely mathematically unbiased, bounces should go to infinity. But this is unrealistic, and as you noted, not visually necessary; for most scenes, after a certain number of bounces, say 10, the amount of contribution to the final color is very very minimal.

So in order to save computing resources, many path tracers have a hard limit to the number of bounces. This adds bias.

That said, it's hard to choose what that hard limit should be. Some scenes look great after 2 bounces; others (say with transmission or SSS) may take up to 10 or 20.

If we choose too low, the image will be visibly biased. But if we choose too high, we're wasting computation energy and time.

One way to solve this, as you noted, is to terminate the path after we reach some threshold of attenuation. This also adds bias.

Clamping after a threshold, will work, but again, how do we choose the threshold? If we choose too large, the image will be visibly biased, too small, and we're wasting resources.

Russian Roulette attempts to solve these problems in an unbiased way. First, here is the code:

void RenderPixel(uint x, uint y, UniformSampler *sampler) {
    Ray ray = m_scene->Camera.CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);

    // Bounce the ray around the scene
    for (uint bounces = 0; bounces < 10; ++bounces) {
        m_scene->Intersect(ray);

        // The ray missed. Return the background color
        if (ray.geomID == RTC_INVALID_GEOMETRY_ID) {
            color += throughput * float3(0.846f, 0.933f, 0.949f);
            break;
        }

        // We hit an object

        // Fetch the material
        Material *material = m_scene->GetMaterial(ray.geomID);
        // The object might be emissive. If so, it will have a corresponding light
        // Otherwise, GetLight will return nullptr
        Light *light = m_scene->GetLight(ray.geomID);

        // If we hit a light, add the emmisive light
        if (light != nullptr) {
            color += throughput * light->Le();
        }

        float3 normal = normalize(ray.Ng);
        float3 wo = normalize(-ray.dir);
        float3 surfacePos = ray.org + ray.dir * ray.tfar;

        // Get the new ray direction
        // Choose the direction based on the material
        float3 wi = material->Sample(wo, normal, sampler);
        float pdf = material->Pdf(wi, normal);

        // Accumulate the brdf attenuation
        throughput = throughput * material->Eval(wi, wo, normal) / pdf;

        // Russian Roulette
        // Randomly terminate a path with a probability inversely equal to the throughput
        float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
        if (sampler->NextFloat() > p) {
            break;
        }

        // Add the energy we 'lose' by randomly terminating paths
        throughput *= 1 / p;

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.org = surfacePos;

        // Reset the other ray properties
        ray.dir = wi;
        ray.tnear = 0.001f;
        ray.tfar = embree::inf;
        ray.geomID = RTC_INVALID_GEOMETRY_ID;
        ray.primID = RTC_INVALID_GEOMETRY_ID;
        ray.instID = RTC_INVALID_GEOMETRY_ID;
        ray.mask = 0xFFFFFFFF;
        ray.time = 0.0f;
    }

    m_scene->Camera.FrameBuffer.SplatPixel(x, y, color);
}

Russian Roulette randomly terminates a path with a probability inversely equal to the throughput. So paths with low throughput that won't contribute much to the scene are more likely to be terminated.

If we stop there, we're still biased. We 'lose' the energy of the path we randomly terminate. To make it unbiased, we boost the energy of the non-terminated paths by their probability to be terminated. This, along with being random, makes Russian Roulette unbiased.

To answer your last questions:

Does Russian Roulette give an unbiased result?

Yes

Is Russian Roulette necessary for an unbiased result?

Depends on what you mean by unbiased. If you mean mathematically, then yes. However, if you mean visually, then no. You just have to choose you max path depth and cutoff threshold very very carefully. This can be very tedious since it can change from scene to scene.

Can you use a fixed probability (cut-off), and then redistribute the 'lost' energy. Is this unbiased?

If you use a fixed probability, you are adding bias. By redistributing the 'lost' energy, you reduce the bias, but it is still mathematically biased. To be completely unbiased, it must be random.

If the energy that would be lost by terminating a ray without redistributing its energy is eventually lost anyway (as the rays to which it is redistributed are also eventually terminated), how does this improve the situation?

Russian Roulette only stops the bouncing. It doesn't remove the sample completely. Also, the 'lost' energy is accounted for in the bounces up to the termination. So the only way for the the energy to be 'eventually lost anyway' would be to have a completely black room.

In the end, Russian Roulette is a very simple algorithm that uses a very small amount of extra computational resources. In exchange, it can save a large amount of computational resources. Therefore, I can't really see a reason not to use it.


Is it better to blur the input or output of an edge detection shader for noise reduction?

Standard blur removes high frequency content from the signal, whereas edge detection usually look into high frequency to detect edges. Be careful on how much blurring to apply to ensure that you don't lose desirable edges. The goal of blurring is to perform noise reduction, so the best would be to come up with a model of the noise present in your images and what your desired edges looks like (their frequency, orientation, what caused them: shading on a smooth object, shadows, etc.) and tailor your blur filter to remove the noise but not the edges.

Seems like you want to implement an edge detector into a shader. Be sure to check what are the sources of noise in your image. For example, if you are displaying pure lambertian surfaces, chances are that there is virtually no noise in your data, so you don't need to apply a smoothing filter before detecting edges.

Also, instead of using a standard gaussian blur, you could look into edge-preserving filters like the bilateral filter. It is much slower, but is good to remove noise without altering too much the edges.

If you are performing filtering and edge detection through convolutions (the standard way of doing it), it doesn't matter if you do one before another, as convolutions are commutative, meaning $a*b = b*a$, where $*$ is the convolution operator, $a$ being the edge detector and $b$ the smoothing kernel.

On non-associative operations (like the bilateral filter), you must ask yourself the question "what is the goal of this operation?" Basically, blurring before the edge detection will diminish the image noise while blurring after will give you blurred edges (potentially caused by noise) instead of sharp edges. The same edges will be found, though, but the blurred one will be less precise. It depends on what you want to achieve: it could be interesting for art, but it just removes interesting information for signal analysis.

Not sure what "applying an edge to a color buffer" means; is the goal to perform edge enhancement?


You don't need to render to texture to achieve an effect like this. Use stencil buffer.

The visibility polygon can be drawn using GL_TRIANGLE_FAN, if you start from the middle:

At the beginning of each frame:

Clear the stencil buffer (and enable stencil test)
Set the stencil function to always pass, and to replace the current value in the buffer to a greater than zero number (for example one)
Draw the visible part as a GL_TRIANGLE_FAN. If you have multiple characters/lights, draw the visibility polygon for each.

Now the stencil buffer contains 1 at visible locations and 0 where the character can't see.

Once you have a stencil buffer like this, you can control which objects you want to draw using glStencilFunc. Do not write to the stencil buffer, just use the stencil test to discard fragments which are visible/invisible.

If you want to draw a background which is visible everywhere, use GL_ALWAYS. For objects which are only displayed in the visible area use GL_EQUAL and 1. Similarly you can draw something (for example a dark quad) just in the invisible parts of the scene with GL_EQUAL and 0.

If it is not clear, how to use stencil buffers, check out this tutorial


Fuse detailed road mesh with regular grid terrain mesh

I think the best way to tackle this problem would be to make an algorithm that works on a single triangle at a time.

You'd look at each terrain triangle individually that the road intersected and process it.

Looking at one triangle individually, you would have to chop that triangle into one or more meshes to make it contour to the road.

Here are some of the cases to consider in the image below.  I'm betting there is some triangle slicing algorithm that would be suited to this, and I'm not sure of a name off hand, but the main contribution I'm proposing here is that you deal with triangles individually.

When cutting the triangle on the x/y plane to make a hole for for the road, you would also want the terrain to take on the z value (height) of the road, for those vertices that touch the road.  This will make the terrain rise up to meet the road.

At this point, you'd have a terrain mesh which did conform to the road but it would be ugly.

Instead of having a nice even looking hill that went up to the road at reasonable slopes, you'd have some triangles that went very sharply up to the road, and others that had a very flat slope.

What you would need to do now is do some kind of "relaxation" algorithm to address this.

Basically, if something is too steep, because there wasn't enough triangle to make the height adjustment, you need to make your whole triangle closer in height to the road and bleed off some of the height adjustment to neighboring triangles.

If something was not steep enough for your liking, because there was a lot of triangle to make a small height adjustment, you may want to break the triangle up into more polygons where it gets closer to the road and have that smaller area handle the change in elevation.

During this second step, when you added a vertex to an edge between two triangles you would also need to adjust the triangle that shares the edge to prevent cracks from showing up in the terrain.

The details are a bit vague, but hopefully the two step approach will help you tackle the problem.

Make the terrain mesh adapt to the road, without caring about how it looks.
Refine the terrain mesh to make it appear to meet the road in a more organic man made way.


Best practices for VRAM management in dynamic environments (e.g. games)

For consoles (and D3D12/Metal/Vulkan), you know everything you need to know to develop an allocation strategy. You know:

How much memory, of each kind, is available.
How big each kind of resource will be. Textures of a particular format and size take up X bytes, require Y alignment, etc.

Given such information, you can devise an allocation strategy for data uploading. Open-world games use streaming of fixed-sized blocks of data. Each sector of the game contains resources that take up X number of bytes, and therefore, artists have X number of bytes of storage to work with. It's their job to fit all of the stuff needed for that region of the game into that much space. There will also be stream blocks for larger collections of regions, or other data that gets shared between all game regions.

For other games, some of them load everything at level-load time. Some do hard-loads in the middle of "levels" (see many Source-engine games). Others stream things more dynamically, but even there, artists are restricted to using resources that total up to X number of bytes.

On consoles, such sizes can be known a priori, but with Vulkan/Metal/D3D12, these have to be queried for a particular piece of hardware. So in the latter case, you have to be more reactive to what's available, but all of the information is available to the user.

For other APIs/platforms (OpenGL/D3D pre-12), things are rather more complicated. Consoles and the lower-level APIs give you direct access to memory allocations, which gives you more flexibility. You can use a piece of memory for vertex data at one point, then put an image into it later on. Memory is divorced from its current use, so it's much easier for an application to choose how to recycle memory.

For GL/D3D-pre-12, things are difficult. A buffer object is a buffer object, and the memory currently associated with it cannot later be used as a texture. Well, it can, but you cannot explicitly ask that to be done. The most you can do is deallocate the buffer, allocate the image, and hope that the buffer's memory was recycled.

For cases where you're doing dynamic loading in fixed-sized blocks, you basically have to give your artists more limitations. Instead of limiting the artists by memory, you have to spell out the limitations by object So you'd say that a stream block can have "X numbers of textures of Y formats and Z size, and W bytes of vertex data". Artists can't sacrifice 4 textures to allow one texture to be 4 times the size, or sacrifice a texture to get more vertex data, or sacrifice some vertex data to get more texture space.

That is what you have to do if you don't trust the implementation's memory allocator to effectively recycle resources. For hard-load cases, it's generally OK to just delete all of the objects and expect the implementation to recycle effectively. But even here, it's a good idea to make sure that you don't try to allocate all of the available storage; leave the implementation a healthy percentage for fragmentation and so forth.

Of course, these kinds of things are a big part of the reason why IHVs have application-specific optimizations: they know up-front which resources were allocated for which purpose, and the driver makers can design allocation schemes that work with these assumptions.

Which is why Vulkan et. al exist: so that application developers don't have to rely on IHVs to work around problems in the API.


What are the differences between GPU instancing and standard instancing?

Unity provides a pretty clear description of GPU instancing. It is exactly the same as instancing. Employees of Unity even spell this out by initially arguing against it years ago before finally implementing it.

Before this, Unity solely relied on batching, which was a technique that combined meshes together into one to reduce draw calls.


Visualizing the Lane-Riesenfeld Algorithm

The Lane-Riesenfeld algorithm subdivides the control polygon of a B-spline to create a new control polygon with the same limit spline.  It's made up of two steps: first, duplicating all of the control points $P_i$ into $P^\prime_{2i}$ and $P^\prime_{2i+1}$; then, moving each point to the midpoint between it and the next point, so $P^\prime_i \rightarrow \frac{1}{2} P^\prime_i + \frac{1}{2}P^\prime_{i+1}$.  These steps are illustrated in the figure:

On the left is an initial control polygon.  In the middle, I've duplicated the vertices.  On the right, the vertices have moved halfway to the next vertex.  Notice that only half of the vertices move in the first movement step: this is because $P^\prime_{2i} = P^\prime_{2i+1},$ but $P^\prime_{2i+1} \neq P^\prime_{2i+2}.$  These eight vertices are a refined control polygon for the linear B-spline defined by the initial four vertices.

Now, we can do a second movement step (without another duplication):

On the left, we've moved each vertex to the midpoint between it and its neighbour; note that all of the vertices move this time (since none are in the same position).  On the right, we have drawn the polygon with these eight vertices.  This is a refined control polygon for the quadratic B-spline defined by the initial four vertices.  You may also recognize this as the same polygon you get by Chaikin corner-cutting, which also gets you the quadratic B-spline.

Now we can perform more duplicate-move-move steps to further refine the polygon, and thus more closely approximate the quadratic B-spline curve:

Here's an animation of this process:

If, instead, we do a third movement step without a duplication step (that is, one duplication followed by three successive movements) we get a refined control polygon for the cubic B-spline defined by the initial four vertices:

An animation of this process:

In general, doing $k$ movement steps after each duplicate gives us the refined polygon for the $C^k$ B-spline.

EDIT: Added animations.


Do the alignement and declaration order of uniforms matter?

I went looking for an answer, so I downloaded AMD's shader analyzer to view the assembly produced when compiled for GCN. In the assembly below vector registers are v# and scalar registers are s#.

It would appear that the uniforms even vector uniforms are passed into the shader as separate scalars, so a vec3 would use 3 scalar registers. The bit I found confusing was the v0 to v4, I'm not sure if v0 is a full 4 float register or a single float in a register, with a full vector register spanning v0 to v3. One way or another it didn't appear to change between the two versions so I can assume the definition order didn't affect the assembly.

http://amd-dev.wpengine.netdna-cdn.com/wordpress/media/2013/07/AMD_GCN3_Instruction_Set_Architecture.pdf

#version 450

uniform vec2 xy;
uniform vec2 zw;

out vec4 v;

void main(){
    v.xy = xy;
    v.zw = zw;
}

shader
  asic(VI)
  type(VS)

  v_cndmask_b32  v0, s0, v0, vcc
  v_mov_b32     v0, 0
  v_mov_b32     v1, 1.0
  exp           pos0, v0, v0, v0, v1 done
  s_andn2_b32   s0, s5, 0x3fff0000
  s_mov_b32     s1, s0
  s_mov_b32     s2, s6
  s_mov_b32     s3, s7
  s_mov_b32     s0, s4
  s_buffer_load_dwordx2  s[4:5], s[0:3], 0x00
  s_buffer_load_dwordx2  s[0:1], s[0:3], 0x10
  s_waitcnt     expcnt(0) & lgkmcnt(0)
  v_mov_b32     v0, s4
  v_mov_b32     v1, s5
  v_mov_b32     v2, s0
  v_mov_b32     v3, s1
  exp           param0, v0, v1, v2, v3
end

#version 450

uniform vec2 xy;
uniform float z;
uniform vec2 zw;

out vec4 v;

void main(){
    v.xy = xy;
    v.zw = zw;
    v.w += z;
}

shader
  asic(VI)
  type(VS)

  v_cndmask_b32  v0, s0, v0, vcc
  v_mov_b32     v0, 0
  v_mov_b32     v1, 1.0
  s_andn2_b32   s0, s5, 0x3fff0000
  exp           pos0, v0, v0, v0, v1 done
  s_mov_b32     s1, s0
  s_mov_b32     s2, s6
  s_mov_b32     s3, s7
  s_mov_b32     s0, s4
  s_buffer_load_dword  s4, s[0:3], 0x10
  s_buffer_load_dwordx2  s[6:7], s[0:3], 0x00
  s_buffer_load_dwordx2  s[0:1], s[0:3], 0x20
  s_waitcnt     expcnt(0) & lgkmcnt(0)
  v_mov_b32     v0, s4
  v_add_f32     v0, s1, v0
  v_mov_b32     v1, s6
  v_mov_b32     v2, s7
  v_mov_b32     v3, s0
  exp           param0, v1, v2, v3, v0
end


Are they just a way to tell Vulkan to use a certain feature implemented by a driver?

Well... yes. Vulkan extensions are simply additional features that Vulkan implementations may provide if they so choose to. They add new functions, structs, and valid enumerators to the API, and they can change some of the behavior of existing functions.

You have to detect whether the implementation supports these extensions before creating an instance/device. And you have to explicitly ask for such extensions when you create the instance/device.


Why does this gl_FragDepth calculation work?

What you are missing is, that in OpenGL's NDC space (i.e. clip space after division by w) all 3 coordinates are in the range $[-1,1]$. So ndc_depth is in the range $[-1,1]$, while your computation assumes a range of $[0,1]$.
So let's take your computation but add an initial step to map $[-1,1]$ to $[0,1]$:
float ndc_depth = Pclip.z / Pclip.w;
float win_depth = (ndc_depth + 1.0) / 2.0;
gl_FragDepth = gl_DepthRange.diff * win_depth + gl_DepthRange.near;

If we do some transformations (and realizing that gl_DepthRange.diff is just convenience for gl_DepthRange.far - gl_DepthRange.near), we get
gl_FragDepth = gl_DepthRange.diff * (ndc_depth + 1.0) / 2.0 + gl_DepthRange.near;

gl_FragDepth = (gl_DepthRange.diff * ndc_depth + gl_DepthRange.diff) / 2.0
               + gl_DepthRange.near;

gl_FragDepth = (gl_DepthRange.diff * ndc_depth + gl_DepthRange.diff +
                2.0 * gl_DepthRange.near) / 2.0;

gl_FragDepth = (gl_DepthRange.diff * ndc_depth +
                gl_DepthRange.far - gl_DepthRange.near + 2.0 * gl_DepthRange.near) / 2.0;

gl_FragDepth = (gl_DepthRange.diff * ndc_depth +
                gl_DepthRange.far + gl_DepthRange.near) / 2.0;

So, to answer your specfic questions:

That's just a result from simplifying the calculation based on the fact that gl_DepthRange.diff is nothing else than gl_DepthRange.far - gl_DepthRange.near.

Because the depth range in NDC space is $[-1,1]$, which has a width of $2$.

It doesn't. These values have nothing to do with the near and far values you used to construct your projection matrix, in fact you could use any arbitrary projection matrix (or none at all). They are the depth range values set with glDepthRange, which have a totally different meaning (okay, let's say slightly different but related).
Basically, and special situations notwithstanding, the projection matrix defines where you near and far clipping planes are in eye-space, i.e. which eye-space range is mapped to the $[-1,1]$ NDC depth range. Then the depth range defines what these NDC near/far values of -1 and 1 are mapped to in the actual depth buffer. So the names "near"/"far" do relate to the near and far clipping planes. But both transformations are basically independent of each other, using the general NDC (i.e. the normalized device space) as a mediator between both stages. This way glDepthRange more or less just extends what glViewport does onto the z-dimension.


How can I draw a tube on basis of position and angle between two 3d points?

The simplest solution is to use a lookat matrix. This is a matrix calculated from a "eye" point, a "target" point and an up vector.

The resulting matrix will then make the Z axis point from the eye to the target and the y axis in the general direction of up.


Always rotate in the direction of mouse movement

Transformations, and especially rotations, may have different visual effects, depending on the order they are applied to the target object. For example, using a rotation $M_R$ and a translation $M_T$ to the object's vertices $v_i$, $M_R \cdot M_T \cdot v_i$ will usually1 result in a different pose than $M_T \cdot M_R \cdot v_i$ does.

90 degree rotation and translation along the x axis applied to the box. Left: translation first ($M_R \cdot M_T \cdot v_i$). Right: rotation first ($M_T \cdot M_R \cdot v_i$).

The key here is to keep in mind to which coordinate system the transformation is applied. Assuming that your to be rotated object initially is in the origin without transformations applied, the first rotation will always yield the desired effect. This is because you're applying the first transformation to the object.

But once your object is already transformed (rotated, translated, whatever) by $M_1$, technically you can already apply the next transformation $M_2$ in two places. Namely,

$M_1 \cdot M_2 \cdot v_i$ or
$M_2 \cdot M_1 \cdot v_i$

You describe that

  so if I rotate the object around the x axis by 180° and then move my mouse left, the object will rotate to the right instead of left (since the y axis is now inverted).

This is, because the second rotation is applied to the object itself, which is case 1 of the above. But you want to have the second case, so you need to apply the rotations the other way around. That way, the new rotation will be applied to the already rotated object, i.e. the outer coordinate system.

Ideally, as long as the user holds the mouse button down you adjust that transformation $M_2$. Calculate the distance the mouse moved from the point the user pressed the button, recalculate the whole rotation ($M_2$) and drop whatever you calculated in the time between the mouse-down event and the previous frame. So for the ongoing rotation, pretend that the user moved the mouse instantly from the start to the current position.

Once the mouse button is released, you can finalize the transformation and multiply it into $M_1$, which is the identity matrix in the beginning. (Assuming you always want to rotate around the world's center.)

1 There are some special cases where the result is the same, e.g. if the translation direction matches the rotation axis.


What's the difference between Clipping and Culling?

When it comes to graphics, the term "culling" by itself doesn't really have a single definition or location within the pipeline. The term is sometimes applied to primitives and it is sometimes applied to fragments. It is sometimes applied to faces based on orientation relative to the camera, and it is sometimes applied to primitives based on various parameters. So it's not useful to talk about "culling" as though it is a single thing; there are a lot of things called "culling".

By contrast, "clipping" with regard to the graphics pipeline is almost always used to refer to what happens to the vertices of primitives after any programmable shading, but before rasterization. Clipping happens differently for different kinds of primitives, but the basic idea is to take a primitive and split it into one or more primitives, such that all of the new primitives are entirely within the clip-space volume.

If a primitive is partially outside of the clip-space boundary, the process of clipping creates one or more new triangles that are entirely within that boundary. The clip-space boundary, through the viewport transform, defines everything which is visible with the current rendering state, so if something is outside of the clip-space boundary, it is not visible. Well, in theory; in practice, there are ways to draw stuff outside of the viewport (point size, line width, multisampling, etc), but nevermind those now.

Now this describes the abstract model of how OpenGL and similar specifications talk about rendering. What actually happens in the hardware is often a different matter. Because clipping has the potential to generate many primitives from one (and thus make the pipeline choke a bit, since it has to run more rasterization cycles), hardware tends to avoid performing actual clipping unless absolutely necessary. So the notion that it is a performance optimization is not really true.

Note that a primitive which is entirely outside of the clip-space boundary gets clipped into nothingness. This particular case is often called "culling" the primitive, which is why the terms sometimes get used in similar places. Point primitives are either rasterized or culled; they don't ever undergo real clipping (since there is only one vertex). But overall, this process of reducing primitives to those within the clip-space boundary is referred to as "clipping", not "culling".

Clipping also allows users to define some user-defined clipping parameters, which clip primitives based on values provided by the previous vertex processing stages. You can use these to implement clipping planes, such that primitives are clipped at the intersection of an invisible plane. All parts of a primitive on one side of the plane are visible, but not on the other side.

Things which are referred to as "culling" formally within rendering APIs:

Face culling. Discarding triangle primitives based on their winding order relative to the viewport.
Vertex-based culling (a relatively new feature), whereby entire primitives can be culled based on vertex processing parameters.

The article you're talking about seems to mostly equate "culling" with visibility culling techniques. These are computations implemented by user code which chooses to render or not render an object based on whether it is visible. These can range from the very simple frustum culling test (ie: asking if the object is even close to the visible region) to complex BSP systems or portal visibility culling (which attempt to detect if objects are entirely behind static geometry). But these are not part of the GPU's rendering pipeline, even if GPU processes are sometimes used to do such testing.

Overall, that article gets more wrong than it does right, so it's best to just pretend you never read it.


Does this 3D rendering method exist? How is it called? Is there any C++ algorithm for it?

Yes, such rendering exists. It is called a parametric surface and there are several methods you can render such entities. In this case you are talking of bezier patches.

Turn the object into enough triangles. The triangilation is called tesselation or meshing depending on what literature you read.
There are lots of applications that do this allready. Most cad applications work this way and many graphics apps like blender can work
this way.

raytrace or raymarch the surface in a shader.

Microdice the surface

etc.

The reason this is not often done is that the on line tesselation is for most intent and purposes just a hindrance to asset developpers. You get the same with subdivision schemes and its just easier. The primary reason to drag the bezier (or nurbs, or t-splines...) is that unless you need to find better intersection of underlying primitives then it is often just a waste of resources. Although subdividing your mesh  via a subdivision surface is allready more useful


The answer in the screenshot is wrong on two counts. First, you're correct that in the middle matrix the −1 should be on the third row, not the second. The other error is that the $u, v$ basis vectors are normalized, but not orthogonalized, so $M$ is not in general a rotation matrix; therefore its inverse is not simply its transpose. (The reflection doesn't actually need an orthonormal basis, so you could fix the answer either way: by using the full inverse for $M$, or by fully orthonormalizing $u$ and $v$.)

The rule of thumb for setting up a local coordinate system (I assume orthonormal) for a given object is to set some of the axes to be lined up with the object, then choose arbitrary perpendicular vectors for any remaining axes.

For a point, there's nothing to line up with, so you can use any arbitrary axes.
For a line, you line up the $u$ axis with the line, then there's nothing to constrain the $v$ axis, so you can choose any arbitrary vector perpendicular to $u$. That's what the $v=u\times(0,0,1)^T$ thing was getting at: it's a convenient formula to get some vector perpendicular to $u$ if you don't otherwise care what direction it is. Except that this won't work if $u$ happens to be parallel to the z-axis (the cross product returns zero then), so you actually need to detect that and switch to a different formula, such as $v=u\times(1,0,0)^T$, in that case.
For a plane, you line up the $u$ and $v$ axes with the plane. So, you can pick any arbitrary $u$ that lies in the plane, then find a $v$ that lies in the plane and is perpendicular to $u$.

Addendum on orthonormalization

Suppose you have a couple of 3D vectors $u, v$ that are non-parallel (linearly independent) but otherwise arbitrary, and you want to generate a 3D orthonormal basis $u', v', w'$ from them. There are two equivalent ways to do it.

Use the Gram–Schmidt process. This works as follows: first, normalize $u$. Then, for each additional vector, subtract off its projection onto the previous vectors, and normalize what's left over. For the 3D case this looks like:
$$ \begin{aligned}
u' &= \text{normalize}(u) = \frac{u}{|u|} \\
v' &= \text{normalize} \bigl( v - (u' \cdot v)u' \bigr) \\
w' &= u' \times v'
\end{aligned} $$
Note that the notation $(u' \cdot v)u'$ means take the dot product of $u'$ with $v$ (giving a scalar), then multiply that scalar by $u'$. This gives the vector component of $v$ that is parallel to $u'$. When this is subtracted from $v$, the remaining component is perpendicular to $u'$.

Finally, $w'$ is set to the cross product of $u'$ and $v'$ to ensure the basis is right-handed. Since $u', v'$ are already orthonormal, we don't need to normalize $w'$—it comes out already normalized. (This isn't strictly part of the Gram–Schmidt process.)

By the way, the Gram–Schmidt process works for making an orthonormal basis in any number of dimensions, not just two or three.
The other way, which only works in 3D, is to use cross products.
$$ \begin{aligned}
u' &= \text{normalize}(u) \\
w' &= \text{normalize}(u' \times v) \\
v' &= w' \times u'
\end{aligned} $$
Here, $w'$ has to be normalized because $u'$ and $v$ are not orthonormal. But then $u'$ and $w'$ are orthonormal, so we can calculate $v'$ from them without a normalize. This gives the same results as the previous process, so it's up to you which you prefer to use.


What is shader warming, and why does Direct3D seem to avoid it?

I hadn't heard the phrase 'shader warming' before, but from some googling around, it seems that some games/engines use this term for a method of forcing the shaders to be fully compiled during level load, as explained in this answer, to avoid hitches during gameplay.

In Metal/Vulkan/D3D12, it's true that shaders are precompiled to a bytecode form, but this is hardware-independent bytecode; the shaders still need to be further compiled to the GPU-specific machine code, which (in those APIs) is done by the graphics driver when a pipeline state object is created. If there are a lot of shaders, that can potentially take awhile. The fully compiled shaders can also be cached, so this only needs to be done the first time you run the app (at least until you update the driver, which could invalidate the cached results).

So my guess about what's going on in your app is that the Metal/Vulkan renderers are calling it 'shader warming' when they go through and create all the pipeline state objects they're going to need—which entails compiling all the shaders. Meanwhile, in the Windows renderer, assuming it's D3D11 or earlier (not D3D12), perhaps the app is relying on the driver to compile shaders in the background, which lets it appear to load faster—at the cost of potentially having hitches a bit later if it tries to use a shader that's not finished compiling yet.


In order to translate a vtkCellLocator, you have only to translate a bounding box of root octant. However this bounding box cannot be publicly accessible. So it is necessary to create a new class inheriting from a vtkCellLocator and provide access functions to get and set it directly. The following is a simple example.

#include <vtkTransform.h>
#include <vtkTransformPolyDataFilter.h>
#include <vtkSmartPointer.h>
#include <vtkSphereSource.h>

#include <vtkCellLocator.h>
#include <vtkObjectFactory.h>

class vtkMyCellLocator : public vtkCellLocator
{
public:
  static vtkMyCellLocator *New();

  vtkSetVectorMacro(Bounds,double,6);
  vtkGetVectorMacro(Bounds,double,6);
};
vtkStandardNewMacro(vtkMyCellLocator);

int main(int, char *[])
{
  auto sphereSource = vtkSmartPointer<vtkSphereSource>::New();
  sphereSource->Update();

  // Create the tree
  auto cellLocator = vtkSmartPointer<vtkMyCellLocator>::New();
  cellLocator->SetDataSet(sphereSource->GetOutput());
  cellLocator->BuildLocator();

  double translate[3] = {10.0, 0, 0};
  auto transform = vtkSmartPointer<vtkTransform>::New();
  transform->Translate(translate);

  auto transformFilter = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
  transformFilter->SetTransform(transform);
  transformFilter->SetInputConnection(sphereSource->GetOutputPort());
  transformFilter->Update();
  cellLocator->SetDataSet(transformFilter->GetOutput());

  double bounds[6];
  cellLocator->GetBounds(bounds);
  for (int i = 0; i < 3; ++i)
  {
    bounds[2*i] += translate[i];
    bounds[2*i+1] += translate[i];
  }
  cellLocator->SetBounds(bounds);

  double testPoint[3] = {2.0, 0.0, 0.0};

  //Find the closest points to TestPoint
  double closestPoint[3];//the coordinates of the closest point will be returned here
  double closestPointDist2; //the squared distance to the closest point will be returned here
  vtkIdType cellId; //the cell id of the cell containing the closest point will be returned here
  int subId; //this is rarely used (in triangle strips only, I believe)
  cellLocator->FindClosestPoint(testPoint, closestPoint, cellId, subId, closestPointDist2);

  std::cout << "Coordinates of closest point: " << closestPoint[0] << " " << closestPoint[1] << " " << closestPoint[2] << std::endl;
  std::cout << "Squared distance to closest point: " << closestPointDist2 << std::endl;
  std::cout << "CellId: " << cellId << std::endl;

  return EXIT_SUCCESS;
}


How do CAD programs render non-polygonal 3D objects?

Depends on the CAD program. But yes basically they do just what you describe they triangulate the model into a mesh/lines on demand and then display that. How they chose to do this depends on the CAD, most modernized codebases probably use buffers. But older cads may in fact use the old immediate mode draw calls.

This said it is possible that the triangulation happens on the GPU. Though i doubt it is implemented that way in many CAD software. There are also other reasons for the CAD to triangulate the meshes, such as intersection testing, CAM path planning and export, that has nothing to do with display as such and would result in the need of implementing both a software and a gpu triangulator doubling workload or and making maintenance difficult.

Most CAD systems allow you to access the triangulated representation, through their API as well as export. But if all else fails you could just attach a debug dll to see what data the CAD dumps to the display card. So for a definitive answer ask your CAD software to tell it to you.


How does the rendering equation incorporate shadowing

Light that is blocked will mean that $L(p*, -w_i)$ under the integral is 0 plus how much light the blocking object itself reflects.

In other words the shadowing is embedded in the incoming light function.


Energy Conservation for Blinn-Phong BRDF

First off, this is probably not the best way of testing energy conservation, because not all incident inclinations are visible in the image. Most of the environment is also black, so a less "peaky" BRDF such as Lambert would come out very dark with no way of telling whether it's energy conserving or not.

What is more commonly used is a so-called "furnace test", where you place a single sphere into a white environment map. When you render this sphere using the BRDF of interest, it should blend into the environment perfectly if the BRDF is energy conserving (make sure no color channels are clipping). Because it is a sphere, you will be able to see the BRDF at all incident inclinations. Additionally, the sphere is convex, so lacking any other objects in the scene, no reflected ray will be occluded, and every pixel of the sphere in the image will be an integral over the BRDF at a different inclination.

If you did this test, you would also notice that your Phong BRDF is most likely not energy conserving. The commonly used normalization term of $\frac{n+1}{2 \pi}$ (or $\frac{n+2}{2 \pi}$, depending on who you ask) only holds at normal incidence. At other inclinations, and particularly at grazing angle, up to half of the reflected rays will land below the surface and are usually terminated (because transmitting through a reflective surface does not make sense).
The correct normalization term that accounts for the part of the Phong lobe that is clipped below the hemisphere is very difficult to compute. Jim Arvo derived a closed-form solution for integer exponents using double-axis moments, summarised in the Global Illumination Compendium, paragraph 31a) and 31b) (page 17/18). These formulas are not very practical for use in a renderer though and I don't suggest implementing them.

It should be clear that if a simple model such as Phong is already this difficult to normalize, then Blinn-Phong is even more difficult. First of all, the $\frac{n+1}{2 \pi}$ normalization term you mention is only correct if you use Blinn-Phong as a microfacet distribution inside a full microfacet BRDF (e.g. Torrance-Sparrow or Walter). If you directly use the Blinn-Phong BRDF, then a more correct normalization term is $\frac{(n+2)(n+4)}{(8 \pi 2^{-n/2} + n)}$, derived in this short article by Fabian Giesen (second page) and also listed on the page you mention. However, this normalization term is again only correct at normal incidence and will not be exact at other inclinations. From what I can tell, no exact normalization term exists.


Resolved by adding precision highp sampler2D to my shaders.

About default precision settings on ES - https://www.khronos.org/files/opengles_shading_language.pdf page 36.


Efficient methods to rasterize a scene of connected triangles

It's an interesting question, because the advice changes over time. Having said that:

On a GPU, by far the most efficient way, as Nathan Reed said in the comment, is to use the rasterisation hardware because that is that it is designed for.
On a CPU, one of the most efficient ways is to more-or-less mimic what rasterisation hardware does.

GPU vendors are understandably quiet about the specifics of their hardware design, but the basic idea hasn't changed for a few decades.

Conceptually split the screen into tiles, e.g. of 8 by 8 pixels.
Traverse the tiles which overlap with a triangle using some variant on Pineda's method.
For each tile, compute a bit mask to find the exact pixels (or samples, if you're subsampling) in the bucket which are inside the triangle, using half-plane calculations. In GPU speak, the samples in a tile which are inside the triangle are called a fragment.
Do whatever processing is needed on the fragment.

On a modern CPU, fragment processing (including the half-space calculations) can be made extremely efficient by using vector instructions.


Vertices of a regular polygon given the incircle radius

what is the relation between radius of the in-circle and circum-circle of a polygon?

That is $cos ( \frac{2\pi}{n}*\frac{1}{2} ) = cos ( \frac{\pi}{n})$

The triangle with edges from the center to the middle of a side, from the center to an adjacent corner and half the side connecting them is a right angle triangle. The angle of the point at the center $ =\frac\pi 2$.

Then it's simply the cosine to get the ratio between adjacent and hypotenuse.


How to calculate the view vector in Phong Lighting Model?

EDIT

While the basic principle of vector calculation below is correct, as Sebastián Mestre pointed out, in lighting you would use the vector from the lighted fragment to the eye position. Therefore, you would in fact use $v = -(point - eye) = eye - point = -point$, so you were actually correct in your "negative" assumption. Still (what might have been more a typo than a misunderstanding), it is the negative position of the lighted fragment, rather than the negative eye (which would still be the origin).

As Sebastián also said in a comment to you, you would usually normalize the view vector after you have calculated it. Sebastián is the MVP.

Calculating a vector $v$ from point $A$ to point $B$ is subtracting the values of $B$ from $A$:

$v = \left(\begin{array}{c}x\\y\\z\end{array}\right) = B - A = \left(\begin{array}{c}b_x - a_x\\b_y-a_y\\b_z-a_z\end{array}\right)$

So in your case, you need to subtract the $eye$ point from the point you're calculating light ($pos$) for:

$
v = point -eye
$

Your $eye$ is the origin, so remaining is only the $point$. Therefore: no, it is not the negative position of the lighted point, but the positive position.


It's not that hard. If you have just planar or angular light sources, you can think of them as one light source split into multiple chunks and the only thing to deal with is how to sample this multi-light and how to compute the PDF of the resulting samples.
Picking probability
First, you need to setup the picking probability $P(l)$ for each light source $l$. The picking probabilities can be any non-negative numbers, you just have to make sure that are always non-zero if the contributions of the respective lights are non-zero. However, the closer the probabilities are to the actual (relative) light contribution, the better the overall performance of your Monte Carlo estimator will be.
Contribution estimates may be:

$1$: naive, uniform probability; good to start with; used by PBRT
$power$: better
$power / distance^2$: even better; closer light sources have greater contributions
…and other improvements: make zero if not visible at all, multiply by the cosine of the inclination angle, etc.

Note that to get the light-picking probabilities, you need to normalise the contribution estimates to the sum of all of them:
$$
P(l)=
  \frac
    {ContribEst_l}
    {\sum_{a\in\mathrm{Lights}}
      {ContribEst_a}}
$$
You could even use more different picking strategies if you find they work just for certain type of scenes and combine them using multiple importance sampling technique…
Sampling
When sampling the lights, you just select one light $l$ with the appropriate probability $P(l)$ and then sample the light surface/angle with its own sampling strategy (with PDF $p_l$). The resulting PDF $p$ will be:
$$
p(x)=P(l) * p_l(x)
$$
Evaluating PDF
When computing the light sampling PDF for given direction (as needed for multiple importance sampling), you just find the light source at the direction, evaluate the PDF of the lights own strategy and, again, multiply it by the probability of picking this particular light.
Conclusion
Now you have a sampling strategy which gives you PDF which can be used in the Monte Carlo estimators (e.g. $f(x)/p(x)$) and the resulting estimator will be unbiased.


Why is the CMYK space four sided in chromaticity charts

There are two aspects to this:

The traditional chromaticity graph is actually a two-dimensional slice of a three-dimensional volume.  The slice usually used is one that matches up (more or less) to the fully-saturated primaries of an RGB-based colorspace; the full RGB colorspace is a tetrahedron.
CMYK is a four-primary colorspace, with primaries that are not linearly independent.  This makes the full volume a complex shape.  Further, at least one (black) and possibly all four of the primaries don't lie on the plane of the 2D chromaticity graph.

The CMYK space is four-sided because it's a 2D slice of a complex 3D shape that (probably) doesn't lie on any of the surfaces of that shape.


Multi-view stereo (MVS) is the general term
given to a group of techniques that use stereo correspondence as their
main cue and use more than two images.

A quote from 'Multi-View Stereo: A Tutorial' by Yasutaka Furukawa and Carlos Hernández.
So to paraphrase: We have a set of images that is larger than two, and use them in a pairwise manner by applying techniques that finde stereo correspondences to reconstruct the object shown in them.
If you want to know more about MVS algorithms that paper I linked seems to be pretty good and at least Furukawa is a name that I know and have read papers from, so he seems to know the topic.


Explanation
With a non-linear scale, you apply different weights for each pixel (or whatever unit you are using). You can use Euclidean distance towards the closest pixels in both directions to determine the weights.
Example
Normalization
For example, in the image below, the red image is the downsampled image (3x3) and the black image is the original image (4x4). The dots represent the sample points (center of the pixels).
You can align the two images by treating each pixel coordinate as a set of two values (x and y) between 0 and 1 (inclusive). Pixels should not be treated as square, so you should treat each pixel as a point rather than a square (when aligning the two images).

In the above image, the four corners would all carry the weight of the original image's four corners, for example. The points in the middle would carry different weights depending on the difference between the four neighboring points.
With numbers:
This is a hypothetical area to sample in larger images (not related to the image above). The distances are between original pixels are normalized from 0 to 1 to simplify the math.
The black dots represent the sample points (pixels) of the original image. The larger red dots represent the sample points of the downsampled image. The red dot in the top left is closer to the black dot in the top left than the other three surrounding dots, so that pixel will have a heavier weight for determining the color of the downsampled pixel.

Horizontal distance from:

Right 2 original (green) pixels = 0.1 pixels
Left 2 original (green) pixels = 0.9 pixels

Vertical distance from:

Top 2 original (green) pixels = 0.1 pixels
Bottom 2 original (green) pixels = 0.9 pixels

Weight for final pixel:
$$weight(0, 0) = (1.0 - 0.1)(1.0 - 0.1) = 0.81$$
$$weight(0, 1) = (1.0 - 0.1)(1.0 - 0.9) = 0.09$$
$$weight(1, 0) = (1.0 - 0.9)(1.0 - 0.1) = 0.09$$
$$weight(1, 1) = (1.0 - 0.9)(1.0 - 0.9) = 0.01$$
Subtract from 1 both the horizontal and vertical distance (because closer pixels should carry more weight). The scale factor in this example is 0.75 in each direction.
Scaling Up vs. Scaling Down
Scaling up in the same  process as scaling down. By normalizing both the original image and the new image to between (0, 0) and (1, 1), you can do either with the same process. As an exercise, in the first figure, if the red points were from the original image and the black points were from the enlarged image, you should be able to see that the process is the same.


Yes, the statement in the W3C doc is correct. You can tell because the background term includes premultiplication by the background alpha αb. If it were non-premultiplied, it wouldn't have an αb in it.


Your actual result is correct in terms of what you can achieve in a single-pass with a fragment shader. It calculates the distance of the fragment's world position to the line segment. That distance is different for each fragment so the stretched black ellipsoid is the correct output.

A fragment shader has only local information about the current fragment being processed. You can't do more with that. Your expected result requires some global knowledge about the scene, since you are searching for the minimum distance across the whole surface.

Without any thoughts on optimization, the naive solution path should involve contextual computation of local minima of the fragment-to-linesegment distance. In order to get such context, you need to compute the distance field first and then access it to find the local minima in a separate pass.

Compute the raw distance field in the first pass (same as your example but no radius subtraction at the end.
Mark the pixels with local minima in the second pass (e.g. by setting the alpha channel to 1). This will be screen-space pass. The best is to have a rgba float texture at output with rgb being fragment's world position and a being 1 for local minima and 0 else.
Paint the black circles/spheres in a third pass (also screen-space).

Having the information about the local minimum of your distance function enables you to output the right shape for around the minima in the last pass. In pseudocode:

vec4 worldPos = tex2D(distance_texture, frag.x, frag.y);
if (worldPos.a > 0.5) //is minimum
    return vec3(0,0,0);
else
{
    //Sample the screen-space neighborhood in larger and larger circles
    // until a fragment with distance.a == 1 is found
    // or until all fragments are more than radius away (in world coords).

    //Minimum world space distance of any sampled fragment.
    float minDist = 0.0;
    int iter = 1;
    //Search for a local minimum until all fragments are more than radius away,
    // since then even if a local minimum gets found it will no longer produce
    // the black shading.
    while (minDist < radius && minDist >= 0.0 && iter < Math.Max(screenWidth, screenHeight))
        minDist = sampleNeighborhood(iter++, worldPos);

    if (minDist < 0.0) //found a local minimum
        return vec3(0,0,0);
    else
        return vec3(1,1,1); //or any other computed shading

    //I am aware that the function call above checks central fragments several times.
    //That should be of course optimized in production code.
}

float sampleNeighborhood(int neigh, vec4 worldPos)
{
    float minDist = float.MaxValue;
    for(int x = Math.max(0, frag.x - neigh); x <= Math.Min(screenWidth, frag.x + neigh; ++x)
        for(int y = Math.max(0, frag.y - neigh); x <= Math.Min(screenHeight, frag.y + neigh; ++y)
        {
            float p = tex2D(distance_texture, x, y);
            float d = distance(p.xyz, worldPos.xyz);
            if (p.a > 0.5 && d <= radius)
                //Return negative distance to indicate a local minimum was found.
                return -1.0;
            else
                minDist = Math.min(minDist, d);
        }
    return minDist;
}

Finding the local minima in the second pass is much easier, since you only need to sample the local 3x3 neighborhood.

I hope that helps you in reaching the desired rendering!


Which perspective projection matrix to use

Even though you can find multiple slightly different formulations for the perspective matrix on the internet, they all do more or less the same thing. They project everything inside a space with the shape of a truncated pyramid (your view frustum) into a space with the shape of a quad and a defined value range (clip space).
Most differences are related to how you define the shape of your frustum. For example, your first matrix uses a height, a width, and an angle together with z-far and z-near value. The second one uses just an x- and y-FOV and the z-far and z-near values to describe the shape of the frustum. This explains the differences in the upper left section.
So what about the lower right section? One big issue here is, that you can choose to use column or row vectors for your vertices/points. Column vectors are multiplied on the right-hand side and row vectors on the left-hand side. The choice you make determines if you have to transpose (switch columns and rows) your matrix or not. The matrices you presented are just targetting different vector representations. I can't tell you without testing/rereading which one is meant for which representation, but if you transpose one of them, they start looking quite similar.
There are still some other minor differences like varying signs, different factors, or a missing addition.
The derivation of your second matrix can be found here.
I can't tell you exactly where those small differences in the first matrix come from, but you can find some possible explanations in this link (Have a look into the "Are There Different Ways of Building this Matrix?" section).
Another good read is this link here.
Which matrix is "correct" depends on how your Graphic API defines its coordinate systems and it's clip space (For OpenGL - you can take the one from the last link I provided). However, there is no "best" choice because they all do more or less the same. Unless you have no special needs for your program (depth precision) or other fancy ideas I can't think of now, just pick the one suited for your API with the inputs (FOV vs width+height+angle) that you like most and start having fun. ;)


Volume rendering is not like ray-tracing, it is like "scene rendering". i.e. there exists several algorithms to render volumes.

One close to ray-tracing is ray-marching, and has may variants. The simplest:
for every pixel in the screen plane, trace a ray starting from the eye point to the screen pixel location, and advance along the ray by constant steps. At each voxel calculate the pixel color and transparency, blend it to the pixel value, and so on.

Blend =

$C_{pix} += T_{pix}*C_{vox}*(1-T_{vox})$
$T_{pix} *= T_{vox}$

with $C_{vox},T_{vox}$ the color and transparency of the current voxel, and $C_{pix},T_{pix}$ the cumulated color and transparency of the pixel.


Main techniques to draw a point at non integer coordinates?

Instead of having a binary on/off per pixel you can instead make a heat map. The more points fall within the pixel the brighter the pixel is. Essentially create a 2D histogram which you then renormalize to get the grayscale for the output.

If you want to blur the points then you can move a portion of the brightness to the next pixel based on the remaining fraction after multiplying to the resolution. If it the coordinate is (0.55313, 0.236831) and your resolution is 1000x1000 then 0.87*0.17 goes to 553, 236, 0.87*0.83 goes to 553, 237, 0.13*0.17 goes to 554, 236 and 0.13*0.83 goes to 554, 237. This is distributed linearly but there are other options.


The effect in the photo is very close to a simple scale-bias per pixel. After a bit of tweaking, I found that applying the transformation: $x' = 0.77x + 38$ to the raw pixel values (as bytes) gives something quite close to your output:

In this case the effect of the scale-bias is to reduce the contrast (scale factor < 1) while keeping the overall brightness about the same.

Example Python code (with numpy and pillow) to do this transformation:

import numpy as np
from PIL import Image
a = np.array(Image.open("input.png"), dtype=np.float32)
b = a * 0.77 + 38
Image.fromarray(np.uint8(np.rint(b))).save("output.png")

For learning how to implement many basic image processing operations, I've found these Image Processing Operator Worksheets from the University of Edinburgh helpful.


How would I create a virtual cabling system?

The shape you’re trying to draw is called a catenary: it’s the shape that a cable/cord of constant density takes when supported at each end. You’ll have to do some research to find a parametric equation for its shape—this page has a start, though it doesn’t let you substitute in the endpoints so you’ll need some additional work there.

Once you have an equation that gives you a point along the curve as a function f(t) = (x,y), where t is in the range [0,1], you can draw it the same way you would any other parametric curve: step along the function from 0 to 1 in small increments, drawing individual line segments from the value of f at the previous increment.

For nicer-looking cables, you can follow that process with a wide line in a dark color followed by a narrow line in a lighter one to give the appearance of shading, draw a shadow with the same curve offset vertically in partially-transparent black, etc.

As to the patch ports themselves… you can make grids pretty easily by using an image for a single port and drawing it at regular intervals across the space you want to cover. To figure out which port a given screen position corresponds to (e.g. if the user clicks on a particular one), get that position, divide it by the size of the interval you drew the ports at, and take the ceiling of that value. In other words, if you’re drawing a port every 20 pixels and the user clicks on (42, 27), dividing each coordinate value by 20 gives you (2.1, 1.35), which ceils to (3, 2)—i.e. the third column and second row.


Ray tracing geometries with same position, translation and scale

This is more of a comment then an answer:
I have resolved this problem in a ray tracer using Constructive Solid Geometry(CSG). Which was done by adding operators that allowed us to create objects that are the union, difference, and intersection of a combination of objects. Here is a wiki page on the subject.
CSG sounds difficult to implement, but it is actually surprisingly straight forward.
But, this system was added for other reasons, it just turned out that we could "subtract" the light from the box (using your post as an example), then the intersection for the light itself falls out naturally since it can no longer hit the box where the two intersect.
This really is just reiterating what everyone has been saying in the comments, every solution is either done by creating real geometry that avoids the problem (like CSG) or use some sort of hack. I posted the CSG answer because it resolves the problem so naturally if it is available.


Techniques for rendering high polycount models

The straightforward solution is to simple not render the high poly model at all. Have a lower detail model that you can switch to once the model is too far away for it to make a difference. There is no point in rendering 20 triangles that all end up in the same pixel.

Next option is to partition the model and then cull the individual parts. You can also make each (set of) patch have different Levels of detail.

This partition can take the form of patches that are mostly flat. That way you can use backface culling on the whole group by taking the dot product of the average normal and comparing it to the look vector, if it is too large only the backsides will be showing and there is no point in rendering that patch.


Multiple Lights Shading - One Pass vs Multiple Passes

It's completely possible to do all pixel lights in the fragment shader (or, say, do 4 at a time) by looping over an array. However, this comes at a great performance cost: You're going to be calculating lighting for every light on every fragment of the scene, even if that geometry isn't actually affected by the light's influence.

I suspect that is reason why Unity doesn't do so: it probably culls the objects and submits only geometry that might be influenced by the respective light on each pass, meaning that shading for a light is only evaluated in the fragments where it may have an effect. This is a trade off between processing less fragments but duplicating the geometry cost.


What is an index buffer and how is it related to vertex buffers?

Could someone ELI5 to me, what is an index buffer and how is it
  related to vertex buffers

Your vertex buffer contains the X and Y coordinates of 5 vertices. They are:

index |  X  |  Y
  0   | 0.0 | 0.0
  1   | 1.0 | 0.0
  2   | 0.0 | 0.6
  3   | 1.0 | 0.6
  4   | 0.5 | 1.0

Your index buffer contains information about which lines to draw between these vertices. It uses the index of each vertex in the vertex buffer as a value.

Since you are drawing lines, each pair of consecutive values in your index buffer indicate a line segment. For example, if the index buffer starts with 0, 2, it means draw a line between the 0th and 2nd vertices in the vertex array, which in this case would draw a line going from [0.0, 0.0] to [0.0, 0.6].

In the following graphic each pair of indices is color coordinated with the line it indicates:

Similarly, if you were drawing triangles instead of lines, you would need to supply an index buffer where each 3 consecutive values indicate the indices of three vertices in the vertex buffer, e.g.

0, 1, 2,
2, 1, 3,
2, 3, 4,


Self intersection artifacts in 3D image raytracing

I want to know if there is a better way to check for self intersection, other than checking the original position?

Since cubes are convex solids, a self-intersection check shouldn't introduce artefacts in your case, while it would in a more general ray-tracer (see my answer about the shadow terminator problem for more about that). In general, it's a good idea for your ray object to have a reference of some kind to the object it came from, as well as the kind of ray it is (primary, shadow, reflection, refraction). Normally the reference should be to the exact polygon, not just the whole object; in your case, it should be the particular cube. You could identify it by the co-ordinates of the cube, if you're using spatial hashing to store your voxel grid.

Then, in your intersection test, if your ray is a shadow ray, you should check whether the origin cube is the same as the cube being intersected, and discard that intersection if so.


Instanced Stereo Rendering vs. Multiple Command Buffers

Yes, there are advantages. You can render different command buffers to different queues.

Modern NVIDIA hardware, in particular, offers a full sixteen separate queues that are capable of rendering independent graphics operations. By putting each eye in a separate queue, you may be able to use the available hardware resources more efficiently.

Then again, it may not. It all really depends on how things work out down in the lowest levels of hardware.

  Regardless, the issue with synchronization that I'm concerned about would be if one eye finishes considerably faster than the other eye. It seems that with stereo instanced rendering, this would be less likely to happen since the workload is more divisible.

It would basically never happen with instanced rendering. But you say that like it's a good thing.

If one eye happens to have more scene complexity than the other eye, then that means instanced stereo rendering will be wasting a lot of vertex processing on vertices for the "weak" eye that will never be seen. Whereas if you built the command buffers separately, you would only submit the stuff that matters for both scenes.

If you submit both CBs to the same queue, the total time spent processing those CBs can be less than the time spent processing the larger CB, since it isn't computing a bunch of triangles that won't be culled. If you submit the two CBs to different queues, how that balances out depends entirely on how well the hardware can adapt to the workload it is given. If it can adapt well, then you'll be fine.

So for this particular case, the multiple CB setup may represent a more efficient use of hardware resources. Granted, this particular case is rather rare.

Synchronization for this case is really not a problem. If one eye is so complicated relative to the other that you drop a frame, then it doesn't matter which form you're using to render. You're going to drop a frame.

With Vulkan, rendering entirely separate stuff to multiple queues isn't difficult. Indeed, the only real synchronization you might need would be to synchronize the presentation of the swapchain images, which is easily done with semaphores.

  The other issue I see where synchronization could cause a performance hit is that you must wait for both threads to complete sending their respective command buffer before sending it again.

No, you don't. If you're reusing CBs, you can always use the VK_COMMAND_BUFFER_USAGE_SIMULTANEOUS_USE_BIT, which allows you to submit a CB that's already in use. And if you're not reusing CBs (that is, your CBs aren't static), you don't care; just double-buffer your command pools so that you aren't trying to reset a CB that's being consumed.

The latter is the general standard method for rendering dynamic scenes with Vulkan.


Spectral path tracing - image color/brightness incorrect

The problem lies mainly in CIE1931XYZ::tristimulusValues() function, where you normalize the resulting color to the luminance of your illuminant which causes that directly observed light source has luminance 1, but everything else is much darker. That is a nice thing to do if you just want to visualize colours of various reflectance spectra under a given illumination, but is probably not the best thing to do in a global illumination renderer.

To be more specific, you cannot expect a reflectance profile to be rendered always completely white just because it reflects all incident energy (constant spectral reflectance equal to 1 everywhere). Neglecting the material directional properties (BSDF, etc.), the properties of the light reflected from a point depend both on reflectance spectral profile of the material, and on the spectral and angular properties of the incident light. The incident light therefore depends not only on the SPD of the illuminant used for your light source, but also on its size and distance from the shaded point. The points closer to the light source will be brighter, while distant points will be darker.

I would suggest removing the luminance normalization in your Spectrum-to-XYZ conversion completely and then adjusting the overall brightness of your scene by multiplying the normalized SPD of the used illuminant (D65 in your case) with some brightness parameter. Whether you multiply it before rendering, during rendering or during the XYZ conversion step is completely up to you. The choice shouldn't affect the result.


The generations of Intel's integrated graphics

You want to look up "Intel Graphics Developer Guides", they are quite thorough.  Here's a link to the latest:

https://software.intel.com/sites/default/files/managed/22/6d/6th_gen_graphics_api_dev_guide_dec_2015.pdf

And here's the architectural Cliff Notes for the past five generations of Intel GPUs:

"In several important respects, the 6th gen Intel Core processors' graphics hardware architecture is similar to that of 5th generation Intel® CoreTM processors"

"The graphics hardware architecture of the 5th gen Intel Core processors is similar to that of the 4th gen processors."

"The graphics hardware architecture in the 4th gen Intel Core processor is similar to that in the 3rd gen processor."

"The new generation of microarchitecture, codenamed Ivy Bridge, provides another jump in functionality and performance over Sandy Bridge microarchitecture."

"With the introduction of the Intel microarchitecture codenamed Sandy Bridge, the graphics processor has moved onto the same die as the CPU, and is now referred to as “processor graphics”. In addition, processor graphics has enjoyed numerous architectural improvements that yield significant performance improvements over previous generations of Intel® integrated graphics parts."