Can the framebuffer in OpenGL jargon or its color buffer end up being the one and the same memory location as the one scanned by a display controller?

Can it be (eg. if I run a video game full-screen) that the compositor could simply declare the application's off-screen buffer the same as its 'screen buffer'

Yes. This is one of the things that "fullscreen exclusive" mode would do for you: it enables the application's swap buffers (front/back default framebuffers) to be directly scanned out by the GPU's display engine, without needing to go through a compositor first, saving a little bit of performance.
Nowadays there is also support for "borderless windowed fullscreen" (see Fullscreen Optimizations for Windows 10 for instance, though I expect other OSes have similar functionality). This is a mode where you still go through the compositor, but the compositor has a shortcut method when it detects that your app is covering the full screen with no other windows rendering on top of it. In that case it can skip the compositing steps and scan out directly from your swap buffers, while retaining the ability to seamlessly switch back to compositing mode if another window should pop up (such as an overlay, or if you Alt+Tab away from the game, etc.)

Does the compositor designate the 'off screen buffer' for the application

Also yes. When you set up an OpenGL context (or any other graphics API context) your app will request the number and format of swap buffers it wants (subject to the choices offered by the OS/driver); then the OS/driver will give you the buffers, and you render into them. You can create your own framebuffers as well, such as depth/stencil buffers or other offscreen buffers for whatever purpose, but ultimately your app will need to put the image you want to show up on the screen into the provided swap buffers. The buffers will be set up so the compositor can also access them. When you do a glSwapBuffers call, you're informing the OS/compositor that you finished a new frame and asking it to please show that on the screen (either by scanning it out directly or by including it in the next compositor frame).


How to convert from Object space into World space (exercise from 3D Math Primer book)

For others who might have this question, here a perhaps more straight forward numeric answer (I'm now working myself on this book as well). Consider that:

We are given object space coordinates that need to be transformed into upright and subsequently world space coordinates.
We know that a vector $\vec{v}$ can be represented as a linear transformation of its axises in relation to the set of basis vectors $\vec{p}$, $\vec{q}$, $\vec{r}$, when these are orthonormal (perpendicular to each other and with length $l = 1$)
The equation to represent this is thus: $\vec{v} = v_{x}p + v_{y}q + v_{z}r$
Given the basis vectors we can computed this for each object reference using the above equation. So, the first example would be:

$\vec{u} = -1 \begin{bmatrix}0.866\\ 0\\ -0.500 \end{bmatrix} + 2 \begin{bmatrix}0\\ 1\\ 0\end{bmatrix} + 0 \begin{bmatrix}0.500\\ 0\\ 0.866\end{bmatrix} = \begin{bmatrix}-0.866\\ 2\\ 0.500\end{bmatrix}$

Then, to compute this in world space, you simply need to add your world space reference to the upright vector. Meaning:

$\vec{w} = o + u$

For the first example in the book that would be:

$\vec{w} = \begin{bmatrix}1\\ 10\\ 3\end{bmatrix} + \begin{bmatrix}-0.866\\ 2\\ 0.500\end{bmatrix} = \begin{bmatrix}0.134\\ 12\\ 3.500\end{bmatrix}$


You can do it by only rendering a flat plane and then making a grid in fragment shader. This is actually a nice exercise.

To get these two color values, you can use color selector, which is available in any graphic software for example MS Paint

Some code hints:

vec2 coord = worldPos.xz / stepSize; //square size in world space
vec2 frac = fract(coord); //fractional parts of squares
//interpolation, grad is smoothness of line gradients
vec2 mult = smoothstep(0.0, grad, frac) - smoothstep(1.0-grad, 1.0, frac);
vec3 col = mix(col1, col2, mult.x * mult.y);


Progressive Path Tracing with Explicit Light Sampling

There are multiple areas in path tracing that can be importance sampled. In addition, each of those areas can also use Multiple Importance Sampling, first proposed in Veach and Guibas's 1995 paper. To better explain, let's look at a backwards 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);
    SurfaceInteraction interaction;

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

        // The ray missed. Return the background color
        if (ray.GeomID == INVALID_GEOMETRY_ID) {
            color += throughput * m_scene->BackgroundColor;
            break;
        }

        // 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 emission
        if (light != nullptr) {
            color += throughput * light->Le();
        }

        interaction.Position = ray.Origin + ray.Direction * ray.TFar;
        interaction.Normal = normalize(m_scene->InterpolateNormal(ray.GeomID, ray.PrimID, ray.U, ray.V));
        interaction.OutputDirection = normalize(-ray.Direction);

        // Get the new ray direction
        // Choose the direction based on the bsdf
        material->bsdf->Sample(interaction, sampler);
        float pdf = material->bsdf->Pdf(interaction);

        // Accumulate the weight
        throughput = throughput * material->bsdf->Eval(interaction) / pdf;

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.Origin = interaction.Position;

        // Reset the other ray properties
        ray.Direction = interaction.InputDirection;
        ray.TNear = 0.001f;
        ray.TFar = infinity;

        // Russian Roulette
        if (bounces > 3) {
            float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
            if (sampler->NextFloat() > p) {
                break;
            }

            throughput *= 1 / p;
        }
    }

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

In English:

Shoot a ray through the scene
Check if we hit anything. If not we return the skybox color and break.
Check if we hit a light. If so, we add the light emission to our color accumulation
Choose a new direction for the next ray. We can do this uniformly, or importance sample based on the BRDF
Evaluate the BRDF and accumulate it. Here we have to divide by the pdf of our chosen direction, in order to follow the Monte Carlo Algorithm.
Create a new ray based on our chosen direction and where we just came from
[Optional] Use Russian Roulette to choose if we should terminate the ray
Goto 1

With this code, we only get color if the ray eventually hits a light. In addition, it doesn't support punctual light sources, since they have no area.

To fix this, we sample the lights directly at every bounce. We have to do a few small changes:

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);
    SurfaceInteraction interaction;

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

        // The ray missed. Return the background color
        if (ray.GeomID == INVALID_GEOMETRY_ID) {
            color += throughput * m_scene->BackgroundColor;
            break;
        }

        // 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 this is the first bounce or if we just had a specular bounce,
        // we need to add the emmisive light
        if ((bounces == 0 || (interaction.SampledLobe & BSDFLobe::Specular) != 0) && light != nullptr) {
            color += throughput * light->Le();
        }

        interaction.Position = ray.Origin + ray.Direction * ray.TFar;
        interaction.Normal = normalize(m_scene->InterpolateNormal(ray.GeomID, ray.PrimID, ray.U, ray.V));
        interaction.OutputDirection = normalize(-ray.Direction);

        // Calculate the direct lighting
        color += throughput * SampleLights(sampler, interaction, material->bsdf, light);

        // Get the new ray direction
        // Choose the direction based on the bsdf
        material->bsdf->Sample(interaction, sampler);
        float pdf = material->bsdf->Pdf(interaction);

        // Accumulate the weight
        throughput = throughput * material->bsdf->Eval(interaction) / pdf;

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.Origin = interaction.Position;

        // Reset the other ray properties
        ray.Direction = interaction.InputDirection;
        ray.TNear = 0.001f;
        ray.TFar = infinity;

        // Russian Roulette
        if (bounces > 3) {
            float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
            if (sampler->NextFloat() > p) {
                break;
            }

            throughput *= 1 / p;
        }
    }

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

First, we add "color += throughput * SampleLights(...)". I'll go into detail about SampleLights() in a bit. But, essentially, it loops through all the lights, and returns their contribution to the color, attenuated by the BSDF.

This is great, but we need to make one more change in order to make it correct; specifically, what happens when we hit a light. In the old code, we added the light's emission to the color accumulation. But now we directly sample the light every bounce, so if we added the light's emission, we would "double dip". Therefore, the correct thing to do is... nothing; we skip accumulating the light's emission.

However, there are two corner cases:

The first ray
Perfectly specular bounces (aka mirrors)

If the first ray hits the light, you should see the light's emission directly. So if we skip it, all the lights will show up as black, even though the surfaces around them are lit.

When you hit a perfectly specular surfaces you can't directly sample a light, because an input ray has only one output. Well, technically, we could check if the input ray is going to hit a light, but there's no point; the main Path Tracing loop is going to do that anyway. Therefore, if we hit a light just after we hit a specular surface, we need to accumulate the color. If we don't, lights will be black in mirrors.

Now, let's delve into SampleLights():

float3 SampleLights(UniformSampler *sampler, SurfaceInteraction interaction, BSDF *bsdf, Light *hitLight) const {
    std::size_t numLights = m_scene->NumLights();

    float3 L(0.0f);
    for (uint i = 0; i < numLights; ++i) {
        Light *light = &m_scene->Lights[i];

        // Don't let a light contribute light to itself
        if (light == hitLight) {
            continue;
        }

        L = L + EstimateDirect(light, sampler, interaction, bsdf);
    }

    return L;
}

In English:

Loop through all the lights
Skip the light if we hit it

Don't double dip

Accumulate the direct lighting from all the lights
Return the direct lighting

Finally, EstimateDirect() is just evaluating $BSDF(p, \omega_{\text{i}}, \omega_{\text{o}}) L_{\text{i}}(p, \omega_{\text{i}})$

For punctual light sources, this is simple as:

float3 EstimateDirect(Light *light, UniformSampler *sampler, SurfaceInteraction &interaction, BSDF *bsdf) const {
    // Only sample if the BRDF is non-specular
    if ((bsdf->SupportedLobes & ~BSDFLobe::Specular) != 0) {
        return float3(0.0f);
    }

    interaction.InputDirection = normalize(light->Origin - interaction.Position);
    return bsdf->Eval(interaction) * light->Li;
}

However, if we want lights to have area, we first need to sample a point on the light. Therefore, the full definition is:

float3 EstimateDirect(Light *light, UniformSampler *sampler, SurfaceInteraction &interaction, BSDF *bsdf) const {
    float3 directLighting = float3(0.0f);

    // Only sample if the BRDF is non-specular
    if ((bsdf->SupportedLobes & ~BSDFLobe::Specular) != 0) {
        float pdf;
        float3 Li = light->SampleLi(sampler, m_scene, interaction, &pdf);

        // Make sure the pdf isn't zero and the radiance isn't black
        if (pdf != 0.0f && !all(Li)) {
            directLighting += bsdf->Eval(interaction) * Li / pdf;
        }
    }

    return directLighting;
}

We can implement light->SampleLi however we want; we can choose the point uniformly, or importance sample. In either case, we divide the radiosity by the pdf of choosing the point. Again, to satisfy the requirements of Monte Carlo.

If the BRDF is highly view dependent, it may be better to choose a point based on the BRDF, instead of a random point on the light. But how do we choose? Sample based on the light, or based on the BRDF?

Why not both? Enter Multiple Importance Sampling. In short, we evaluate $BSDF(p, \omega_{\text{i}}, \omega_{\text{o}}) L_{\text{i}}(p, \omega_{\text{i}})$ multiple times, using different sampling techniques, then average them together using weights based on their pdfs. In code this is:

float3 EstimateDirect(Light *light, UniformSampler *sampler, SurfaceInteraction &interaction, BSDF *bsdf) const {
    float3 directLighting = float3(0.0f);
    float3 f;
    float lightPdf, scatteringPdf;

    // Sample lighting with multiple importance sampling
    // Only sample if the BRDF is non-specular
    if ((bsdf->SupportedLobes & ~BSDFLobe::Specular) != 0) {
        float3 Li = light->SampleLi(sampler, m_scene, interaction, &lightPdf);

        // Make sure the pdf isn't zero and the radiance isn't black
        if (lightPdf != 0.0f && !all(Li)) {
            // Calculate the brdf value
            f = bsdf->Eval(interaction);
            scatteringPdf = bsdf->Pdf(interaction);

            if (scatteringPdf != 0.0f && !all(f)) {
                float weight = PowerHeuristic(1, lightPdf, 1, scatteringPdf);
                directLighting += f * Li * weight / lightPdf;
            }
        }
    }

    // Sample brdf with multiple importance sampling
    bsdf->Sample(interaction, sampler);
    f = bsdf->Eval(interaction);
    scatteringPdf = bsdf->Pdf(interaction);
    if (scatteringPdf != 0.0f && !all(f)) {
        lightPdf = light->PdfLi(m_scene, interaction);
        if (lightPdf == 0.0f) {
            // We didn't hit anything, so ignore the brdf sample
            return directLighting;
        }

        float weight = PowerHeuristic(1, scatteringPdf, 1, lightPdf);
        float3 Li = light->Le();
        directLighting += f * Li * weight / scatteringPdf;
    }

    return directLighting;
}

In English:

First, we sample the light

This updates interaction.InputDirection
Gives us the Li for the light
And the pdf of choosing that point on the light

Check that the pdf is valid and the radiance is non-zero
Evaluate the BSDF using the sampled InputDirection
Calculate the pdf for the BSDF given the sampled InputDirection

Essentially, how likely is this sample, if we were to sample using the BSDF, instead of the light

Calculate the weight, using the light pdf and the BSDF pdf

Veach and Guibas define a couple different ways to calculate the weight. Experimentally, they found the power heuristic with a power of 2 to work the best for most cases. I refer you to the paper for more details. The implementation is below

Multiply the weight with the direct lighting calculation and divide by the light pdf. (For Monte Carlo) And add to the direct light accumulation.
Then, we sample the BRDF

This updates interaction.InputDirection

Evaluate the BRDF
Get the pdf for choosing this direction based on the BRDF
Calculate the light pdf, given the sampled InputDirection

This is the mirror of before. How likely is this direction, if we were to sample the light

If lightPdf == 0.0f, then the ray missed the light, so just return the direct lighting from the light sample.
Otherwise, calculate the weight, and add the BSDF direct lighting to the accumulation
Finally, return the accumulated direct lighting

.

inline float PowerHeuristic(uint numf, float fPdf, uint numg, float gPdf) {
    float f = numf * fPdf;
    float g = numg * gPdf;

    return (f * f) / (f * f + g * g);
}

There are a number of optimizations / improvements you can do in these functions, but I've pared them down to try to make them easier to comprehend. If you would like, I can share some of these improvements.

Only Sampling One Light

In SampleLights() we loop through all the lights, and get their contribution. For a small number of lights, this is fine, but for hundreds or thousands of lights, this gets expensive. Fortunately, we can exploit the fact that Monte Carlo Integration is a giant average. Example:

Let's define

$$h(x) = f(x) + g(x)$$

Currently, we're estimating $h(x)$ by:

$$h(x) = \frac{1}{N} \sum_{i=1}^N f(x_i) + g(x_i)$$

But, calculating both $f(x)$ and $g(x)$ is expensive, so instead we do:

$$h(x) = \frac{1}{N} \sum_{i=1}^N \frac{r(\zeta, x)}{pdf}$$

Where $\zeta$ is a uniform random variable, and $r(\zeta, x)$ is defined as:

$$r(\zeta, x) =
\begin{cases}
  f(x), & 0.0 \leq \zeta \lt 0.5 \\
  g(x), & 0.5 \leq \zeta \lt 1.0
\end{cases}$$

In this case $pdf = \frac{1}{2}$ because the pdf must integrate to 1, and there are 2 functions to choose from.

In English:

Randomly choose either $f(x)$ or $g(x)$ to evaluate.
Divide the result by $\frac{1}{2}$ (since there are two items)
Average

As N gets large, the estimate will converge to the correct solution.

We can apply this same principle to light sampling. Instead of sampling every light, we randomly pick one, and multiply the result by the number of lights (This is the same as dividing by the fractional pdf):

float3 SampleOneLight(UniformSampler *sampler, SurfaceInteraction interaction, BSDF *bsdf, Light *hitLight) const {
    std::size_t numLights = m_scene->NumLights();

    // Return black if there are no lights
    // And don't let a light contribute light to itself
    // Aka, if we hit a light
    // This is the special case where there is only 1 light
    if (numLights == 0 || numLights == 1 && hitLight != nullptr) {
        return float3(0.0f);
    }

    // Don't let a light contribute light to itself
    // Choose another one
    Light *light;
    do {
        light = m_scene->RandomOneLight(sampler);
    } while (light == hitLight);

    return numLights * EstimateDirect(light, sampler, interaction, bsdf);
}

In this code, all the lights have an equal chance of being picked. However, we can importance sample, if we like. For example, we can give larger lights a higher chance of being picked, or lights closer to the hit surface. You just have to divide the result by the pdf, which would no longer be $\frac{1}{\text{numLights}}$.

Multiple Importance Sampling the "New Ray" Direction

The current code only importance samples the "New Ray" direction based on the BSDF. What if we want to also importance sample based on the location of lights?

Taking from what we learned above, one method would be to shoot two "new" rays and weight each based on their pdfs. However, this is both computationally expensive, and hard to implement without recursion.

To overcome this, we can apply the same principles we learned by sampling only one light. That is, randomly choose one to sample, and divide by the pdf of choosing it.

// Get the new ray direction

// Randomly (uniform) choose whether to sample based on the BSDF or the Lights
float p = sampler->NextFloat();

Light *light = m_scene->RandomLight();

if (p < 0.5f) {
    // Choose the direction based on the bsdf
    material->bsdf->Sample(interaction, sampler);
    float bsdfPdf = material->bsdf->Pdf(interaction);

    float lightPdf = light->PdfLi(m_scene, interaction);
    float weight = PowerHeuristic(1, bsdfPdf, 1, lightPdf);

    // Accumulate the throughput
    throughput = throughput * weight * material->bsdf->Eval(interaction) / bsdfPdf;

} else {
    // Choose the direction based on a light
    float lightPdf;
    light->SampleLi(sampler, m_scene, interaction, &lightPdf);

    float bsdfPdf = material->bsdf->Pdf(interaction);
    float weight = PowerHeuristic(1, lightPdf, 1, bsdfPdf);

    // Accumulate the throughput
    throughput = throughput * weight * material->bsdf->Eval(interaction) / lightPdf;
}

That all said, do we really want to importance sample the "New Ray" direction based on the light? For direct lighting, the radiosity is affected by both the BSDF of the surface, and the direction of the light. But for indirect lighting, the radiosity is almost exclusively defined by the BSDF of the surface hit before. So, adding light importance sampling doesn't give us anything.

Therefore, it is common to only importance sample the "New Direction" with the BSDF, but apply Multiple Importance Sampling to the direct lighting.


Volume Rendering - Dataset not rendered correctly

Sorry to post such a silly solution but apparently the problem was reading 16 bit uint data into 8bit uint array. This was the reason for inexplainable noise, interleaving and other artifacts.

After I explored more about what the medical datasets actually represent and what flavors they come in (8,16,32bit) I managed to solve the problem.


glMultiDrawArraysIndirect does not work when specifying an offset

The problem was my class handling the population fo the draw indirect commands.

Each time I added a new command to the indirect buffer I also incremented the baseInstance, so the cube had baseInstance = 0 and the sphere had baseInstance = 1.

The reason why it doesn't render the sphere is because glMultiDrawArraysIndirect pretty much evaluates as glDrawArraysInstancedBaseInstance(mode,  cmd->first,  cmd->count,  cmd->instanceCount, cmd->baseInstance); (as per the OGL wiki).

The drawID is calculated as (gl_InstanceID / divisor) +baseInstance (link), which means that the drawID is out of range when trying to render with an offset, as the sphere now appears as the first element in the indirect buffer.


Efficiently sampling specific surfaces using MIS in path tracing

I will not answer your questions directly, but will try to change your mindset a bit. I hope got your intentions reasonably well...

The usual way the path tracing is implemented is that the reflection integral is split into direct and indirect illumination and they are estimated separately.

The direct illumination can be estimated using more sampling strategies which are combined together using the multi-sample multiple importance sampling method by Veach. Typically, there is just BSDF sampling and light sampling used, but you can also add things like light portals sampling and other more sophisticated stuff. I guess this is where you are headed.

The indirect illumination is usually estimated by just one sampling strategy (typically BSDF sampling) because you want to avoid exponential growth of recursive paths count. That's the reason why path tracing has poor performance for indirect illumination reflected/refracted through glossy/specular materials.

A nice thing about the split integral is that you can share samples between the direct and indirect integral. If you have two integrals over the same domain and have a sampling strategy that can be used for both, then you can generate a sample to estimate the first integral and then use the same sample for the second integral as well. There is no problem with correlated samples here because we have two separate computations. In practice you generate a BSDF sample with some PDF, shoot a ray in that direction and use the intersection point to estimate both direct and indirect integral separately.

PS: I guess that for indirect illumination estimator you could possibly use the one-sample multiple importance sampling method by Veach to incorporate the light portals and other sampling strategies, but I have some doubts about the efficiency of the resulting estimator, especially when reusing the samples for the direct illumination integral.


You already have a 2D parametrisation, don't you? One of the dimensions is the longitudal axis (in mm?) and the other is the circumferential axis (in degrees). The only problem I see is when you have rectangles wrapping around the 0/360-degree boundary. One workaround for that would be to duplicate each rectangle (or just the ones on the boundary) so that you get one copy on each side of the boundary. Does that help? And if you need to have a distance, and not an angle, for the circular dimension, that is readily available as $l = r a \pi/180$, where $r$ is the radius and $a$ is the angle in degrees.


adobe illustrator art brush, how does it work?

These rely on a transformation from global coordinates to path local coordinates, where the local axes are along and out (or more often u and v). In essence this all works on a principle similar to how offset works.

Image 1: Depiction of local axes.

Now the first one ($A$) is a sweep of an oval. Which is technically the most challenging of the four. In essence You offset edges directionally and cap them with a fixed directional oval. Although this has all kinds of interesting edge cases. A lot of mathematical lore is known about sweeping shapes as this is somewhat hard to do in 3D and doing so with an arbitrary shape is useful for the cad crowd. At it simplest its just sampling many times, which allows for change in shape along line, but quite sophisticated approaches are available.

Image 2: Oval sweep.

The second is simply (bombing) a shape in random places in the new coordinate system. Mathematically this is pretty simple. Generate a number in along direction then find that point. Then generate a second number and multiply that by the normal direction vector of your path then add this to the point.

The third and fourth just do this for a lot of points read in form a Cartesian template object. Again you should familiarize yourself with how offset and graduated offset is done. You can find a lot of references but one of my favorite online resource for Béziers is in the Pomax Bézier primer. Its just lot of small offset lines to transform from.


Can gl_WorkGroupSize be set outside the compute shader?

Core OpenGL has no mechanism for this. The ARB_compute_variable_group_size extension provides such a feature. It doesn't allow for indirect compute operations though.

It doesn't seem to have much AMD support though (outside of open source Linux drivers). And Intel support is nil.


As long as no pixel on the texture is used twice, you can render the geometry and display the uv-coordinates (texture coordinates). Usually they are a combination of red and green (2d).
When using a brush, you only need to read out the uv-coordinates you've hit and color the texture at position (uv) in the desired color.
Hint:
paint onto a second texture, so you'll not destroy the original texture color.
Drawbacks:
Sometimes the texture resolution of different triangles within a geometry vary. So at some faces it will look nice, and on others it will be pixelated.
As said before, if the mesh (geometry) uses the same part of a texture twice (two faces aiming at the same texture coordinates), and you draw on one of them, the other will be painted as well.


File format for swept profile with changing normal

Well realistically your choices are either IGES or STEP. IGES is slightly simpler but you will not successfully write either format without buying the standard (which in case of step is actually so many pages that by the time you have read them a year has passed, and lots of money is wasted on paper). You could also use one of the nonstandard kernel internal formats like Parasolid as its also commonly supported. Most other formats are polygon data so not suitable.

I would explore an alternative option to writing your own exporter, buy a commercial or take a free CAD kernel and let the kernel (open CASCADE) do the exporting for you.

On the other hand this is not hard to do inside your CAD application. Which you will need anyway to verify the result. Most CAD applications have a COM API's (don't know about NX most likely, but both creo and catia do have this) so you can edit, even do the generation from a text editor via python if you so wish. This is by far the least amount of effort (and it's usually less work to set up than using just the kernel).

Image 1: Swept with spine angle $60 \cos(4 t)+110$ and polar trajectory with $r = \cos(4 t)+3$ in Creo. All input from program injection.

So the approach depends on why you want to do what you do. And what your budget/time constraint is.

PS: NURBS (the underlying format of CAD data) solutions are not necessarily analytic, but parametric fits to analytical solutions. Sometimes these are exact sometimes not so exact as not all functions can be fitted to rational spline (though many can). Even when they can fitting is a more general solution.


Outside of GL_R11F_G11F_B10F, 3-channel formats are not allowed for image load/store. And you'll note that rgb32f is not in the list of valid format qualifiers for GLSL variables, nor are any other 3-channel formats (except r11f_g11f_b10f).


CUDA/OpenCL vs OpenGL compute shader for Path tracer

To first clear your confusion around the terms:

GPGPU stands for General Purpose computing on GPUs

CUDA is the specific NVIDIA API to perform GPGPU only on their hardware

OpenGL is a graphics specific API and is vendor independent
OpenCL is a parallel programming compute API and is vendor independent

Compute Shaders are a way to perform general purpose computations within
a rendering system but without being part of a fixed rendering architecture.
They exist for OpenGL (https://www.khronos.org/opengl/wiki/Compute_Shader)
but also for other graphics APIs like DirectX (DirectCompute) etc.

OpenCL is not only built for GPUs and can also be used to write code
for example to run massive parallel computations on supercomputers.

NVIDIA develops their own ray tracer based on CUDA and GPU computing.
The project is called Optix (https://developer.nvidia.com/optix)

To your performance question. It depends, as always on your setup and use case.
Some people did already the effort to measure differences.
One paper for example is this:
https://arxiv.org/ftp/arxiv/papers/1005/1005.2581.pdf

Another one is this:
https://dl.acm.org/citation.cfm?id=2066955

Or one thesis about it:
https://liu.diva-portal.org/smash/get/diva2:909410/FULLTEXT01.pdf

In their comparisons did CUDA often show good results. But when they for example used image processing algorithms that can be parallized very efficiently from the GPU then the comparison is nearly worthless compared what you might do in your path tracing code.

If you do not use any "classic" rendering, which you wanna combine your path traced results, I see not why you should use compute shaders for a pure path tracer. If platform independence is important OpenCL might be interesting for you since you then can take the code and run it on a complete other hardware for comparison. If you will stick to NVIDIA GPUs anyways you will likely get better performance using directly CUDA.


What explains the high specularity of metals?

Warning: I am not a physicist.

As Dan Hulme already explained, light can't travel through metals, so dealing with IOR is a lot more... complex. I will answer why that happens and how to calculate the reflection coefficient.

Explanation: Metals are filled with free electrons. Those electrons react to external fields and reposition until electrostatic equilibrium is met (the electric field is zero inside a conductor in electrostatic equilibrium). When electromagnetic waves hit a metallic surface, the free electrons move until the field that they create cancels the field of the incoming wave. Those electrons grouped together radiate a wave going out nearly the same as the one that hit the surface (i.e. with very low attenuation). How much is attenuated depends on the material properties.

From this explanation it is clear that conductivity is a key part of the high reflection coefficient on metals.

Math-wise, what you are missing is the complex index of refraction. On good conductors, such as metals, the complex term of the IOR is relevant and key for explaining this phenomena.

Practically, in rendering, achieving good metal parameters is more visual based. Artists adjust to their preference until it looks believable. Often you see a metalness parameter with specific handling for materials marked as metal.

Involved answer:

The complex index of refraction can be seen if we use Ohm's Law $J = \sigma \vec{E}$, which holds for conductors, on the Ampère-Maxwell equation using sinusoidal waves $\vec{E} = e^{i\omega t}$:

$$
\vec{\nabla} \times \vec{H} = \sigma\vec{E} + \frac{\partial \vec{D}}{\partial t} = \sigma \vec{E} + i\omega \epsilon \vec{E}
$$
$$
 = i\omega \left( \epsilon - i \frac{\sigma}{\omega} \right)\vec{E} = i \omega \epsilon_m\vec{E}
$$

Note how we can interpret that whole term as a complex permittitivity $\epsilon_m$ and that $\sigma$ is the conductivity of the material.

This affects the IOR, as its definition is given by:

$$
  n' = \sqrt{\frac{\epsilon_m}{\epsilon_0}} = \sqrt{\frac{\left(\epsilon - i \sigma / \omega\right)}{\epsilon_0}} = n_{\text{real}} + in_{\text{img}}
$$

This shows how $n'$ can be complex. Also, note how very good conductors have a relevant complex term, as $\sigma \gg \epsilon_0 \omega$. Since it would take a lot, I will skip some steps with a reference, page 27: it can be shown that, since $\sigma \gg \epsilon_0\omega$, (we are dealing with $\omega$ of the visible spectrum):
$$
  n_{\text{real}} \approx n_{\text{img}}
$$

and reflection from metals with normal incidence, from a medium with IOR $n$, given that $n' \gg n$:

$$
  R = \frac{(n_{\text{real}} - n)^2 + n_{\text{img}}^2}{(n_{\text{real}} + n)^2 + n_{\text{img}}^2} \approx 1
$$

Agreeing that a good conductor is, in general, a good reflector.

The famous Introduction to Electrodynamics from Griffiths, pages 392-398, explains this and a lot more in a similar fashion.


Confusion between usages of linear RGB and sRGB

The monitor interprets its input data as encoded in sRGB color space. This means that to actually display them, i.e. to get brightness of the subpixels it does the transformation $C_{\mathrm{sRGB}}\to C_{\mathrm{linear}}$. This transformation is basically application of gamma of $\approx2.2$ to the value, so that output brightness will be superlinearly dependent on input data. Since non-extremal input values are $<1$, they'll all be darker than with linear response.

Now if you have rendered an image in linear RGB color space, you need to compensate the transformation the monitor does. Thus you convert the image to sRGB, i.e. apply the $C_{\mathrm{linear}}\to C_{\mathrm{sRGB}}$ transformation. This transformation is inverse of what the monitor does: it applies a gamma of $\approx1/2.2$ to the values, so all non-extremal values will be larger than original.

In the OP I was coming to the wrong conclusion because of a mistake when trying to match the output of monitor with my input. To match it correctly one needs to either apply $C_{\mathrm{sRGB}}\to C_{\mathrm{linear}}$ transformation to input data to simulate what the monitor does with them, or apply backward transformation to the measured output — to undo monitor's processing.


Gradient descent (Not ML) on arbitrary meshes

Maybe you can find a heuristic to detect when it gets stuck (e.g. gradient magnitude is too small, or it stays on the same triangle or returns to a recently visited triangle too many times in a row, or it hasn't moved at least X distance over the last N steps, etc) and just take a random step in some direction, e.g. to some randomly chosen neighboring triangle? That might help shake it out of the difficult region and get back to smooth descent.
For an even lower-tech solution, if you don't mind some randomness you could just inject occasional random moves all the time (without attempting to detect whether it's stuck). Kind of like stochastic gradient descent in ML.
If the gradient you're following is a geodesic distance, I think it can have saddle points, but not local minima (except for the global minimum), is that right? If that's a constraint you can rely on, you could try to use a second-order method, finding an estimate of the local Hessian and looking for eigenvectors with negative eigenvalues, which should be directions of principal curvature that will accelerate the gradient descent, and in particular, will pull you out of saddle points. I'm sure there must be a literature on such methods although I don't know much about it.


What is the PDF for path tracing in the paper "Learning the light transport the reinforced way

I have to answer my question.
Consider that all patches have the same value. Then choosing a patch uniformly and selecting a random point uniformly on that patch is as the same as sampling the hemisphere uniformly. If $p_s=Pr(Q_k \in Q)=1/n$, which n is the number of patches, I can conclude that $p_\omega=p_s*(n/2\pi)=1/2\pi$ .

And if patches don't have the same value, then $p_s=Pr(Q_k \in Q)=q_k/(q_1+...+q_n)$ and
$p_\omega=p_s*(n/2\pi)=(q_k/(q_1+...+q_n))*(n/2\pi)$


It seems that nothing was wrong with the code I posted.
I deleted my complete build and rebuild it, and now it seems to work.
I don't see anything on screen but i think that error is unrelated.

I can see the matrix in the nsight debugger.


Do you have any "input data" from which you constructed the matrix (like fov, near plane distance, etc.)? Or is the matrix itself the only data you have?

EDIT:
Ok, so if you want to know the size of the near plane while you have viewProj matrix available do this:

a) compute viewProj^-1

b) vertices of the near plane in clip-space have coordinates (in OpenGL):
(-1, -1, -1, 1)
(-1,  1, -1, 1)
( 1, -1, -1, 1)
( 1,  1, -1, 1)

and on the far plane:
(-1, -1, 1, 1)
(-1,  1, 1, 1)
( 1, -1, 1, 1)
( 1,  1, 1, 1)

It's a unit cube spanning (-1, -1, -1) and (1, 1, 1) coordinates with w=1

c) take each of those vertices, mutiply by viewProj^-1 and finally divide all resulting components by the w-component. That gives you the frustum's vertices in world space. From this you can easily find world space planes constraining the view frustum or the bounding box or whatever you like.

I think it's not possible to extract left/right/bottom/top/near/far from the matrix directly. But since you can compute the world space coordinates of frustum's corners (as shown above) you can easily find left/right/bottom/top/near/far.


Can you do GPGPU on integrated graphics like Intel HD 630 on Core i3?

CUDA is only available on nvidia cards.

The open technology is OpenCL: https://es.wikipedia.org/wiki/OpenCL

Another alternative is to use Compute Shaders. DirectX, OpenGL, and Vulkan all support compute shaders


It's hardware dependent, but as it stands, you shouldn't expect a performance drop on modern dedicated GPUs. Typically they have a pipeline state cache which fits such a small number of states easily.

The only case when that cache even becomes relevant is if the workload performed in that state is insufficient to saturate the shaders, so execution of subsequent drawcalls overlaps. At that point you can eventually exceed the size of the cache if you have too many different pipeline states and effective utilization drops.

Except that when you do, you will also have trouble combining all the states into one, due to a limit on texture bindings per state and stage, so you would still end up with multiple states. Depending on the API you will have different limits, lower bound being set e.g. by OpenGL 3.x spec requiring 16 texture binding slots per shader, but realistically it doesn't even go much higher than 40-64 binding slots per shader even in high end hardware.

At that point you are better off using a texture array instead - the level of indirection associated with that doesn't cause any significant performance penalty any more, on modern architectures. You "only" have to group your textures by size and tiling at that point.

The alternative is to bake the individual textures into a a single, large texture, which does still allow a certain level of mipmapping, but only down to a level at which pixels from different source textures don't blend. Also you loose the ability to use native tiling, but have to do it manually in shader code.


OpenGL: Why glTextureSubImage2D crashes?

Pixel transfer operations are governed by more than just your parameters to a pixel transfer function. They are also governed by parameters passed to the glPixelStore functions.

Of particular interest to you is GL_PACK_ALIGNMENT. This tells OpenGL what the byte alignment for each row of data is in your array. By default, this value is four, not one.

This means that OpenGL will believe that the second row of your image starts at this address: &myData[1028], since 1025 is not aligned to 4 bytes. And that means it will run off the end of your allocated buffer eventually.

If you are using data with 1-byte alignment for rows, you must inform OpenGL of this by calling glPixelStorei(GL_PACK_ALIGNMENT, 1);. Note that you should not do this by default; do it only when dealing with byte aligned data.


Leaving the actual derivation of mathematics aside, I'll try to give a general description of how can we combine the NDF with importance sampling in spherical coordinates.
What is NDF ?
The problem is the following. Intiutively, we need the NDF for "zooming in" from the realm of macrosurfaces to that of microfacets.
Mathematically it is needed for the following equivalence:
$$dA_h = NDF(h) dw_{h} A$$.

$A$: flat macrosurface patch

$dw_{h}$: is the area of a portion of a unit spherical cap defined by the orientation $h$

$dA_{h}$: is the total area that is being approximated by the orientation $h$

see here for more.
These terms are understood better when we considere the following equivalence:
$$d\Omega = 4\pi (\frac{dS}{A})(\hat{r} \cdot \hat{n})$$
You can read this equivalence as the following:
One can approximate a solid angle subtended by a flat surface whose area equals $dS$, orientation $\hat{n}$, and distance from the viewer $r$ where $\hat{r} = \frac{\vec{r}}{r}$ and $A = 4\pi r^2$. The $\vec{r}$ is a position vector of the surface area $dS$ with respect to the center of the sphere. You can read more about it in wikipedia
The point here is that whenever we consider approximating solid angles with projected surfaces, we need that direction and distance term.
Whether we consider solid angles over microfacets or macrosurfaces don't make a difference. From this perspective, you can think of $NDF$ as "regulation" term for two different granularities involved in approximating the same structure.
Since we are modeling the same structure from two different granularities, we can't "touch" properties that coincides with physical/geometric properties of the structure, but we can "touch" the regulation term, hence different ways of computing ndf are provided by different researchers.
Whatever it maybe though, it has to satisfy the following equivalence for unit sphere:
$$1=\frac{1}{A}\int_{S} (n \cdot h) dA_{h} = \int_{\Omega} NDF(h) (n \cdot h) dw_{h}$$
Where S is the domain of the area of the macrosurface patch. This reads simply as the following: The area we are calculating as a spherical cap using solid angle bits, and the area we are calculating as a macrosurface patch using macro surface patch bits has to be the same.

What is Importance Sampling ?
Here is the second bit of information necessary to understand the procedure.
Let X be a random variable, ie a function whose signature is $X: M \to R $, where $M$ is the sample space and $R$ is the set of real numbers. Expected value $E[X]$ is given under the probability distribution $P$ as the following $E[X;P] = \frac{1}{n} \sum_{i=1}^n x_i$ where n is the total number of samples drawn from the distribution and $x_i$ is sample generated from $X$ using the distribution P.
The precision of the estimate calculated this way depends on the variance of the random variable $X$ under the probability distribution $P$. If we have infinite number of samples, then its okay to use low precision estimates, but we never do. In fact, it would be great if we have a small number of samples, with a more precise way of computing the estimate.
Now think of the following, we have some random variable X, that outputs values according to some distribution, and we want to understand how X works under that distribution. Wouldn't it be better to select the input values for X so that the graph made by the output values look more like the actual graph of the X under distribution P. This is what we are trying to do with importance sampling. The key here is to understand that though we may not know how X maps values in its domain to its codomain, we may know how the probability distribution for choosing samples works. Hence if we choose a function that is similar to random variable's actual probability distribution, we can significantly reduce variance.

How they interact ?
Now coming to your question. We have seen that in order to use importance sampling we need to use function that resembles to the actual probability distribution used for sampling the spherical cap area. So if our function is $f=NDF(h) (n \cdot h)$, we need use a function that is similar to the actual probability distribution of this function to sample values from it in order to get an estimate with higher precision.
Since we are dealing with points on unit sphere, we are dealing with spherical coordinates. The vectors associated with these points are also in spherical coordinates. Now spherical coordinates are defined using two angles and a distance to the center of the sphere. In unit sphere the distance is 1. This leaves us with
the two angles when considering vectors.
The probability distribution associated with $f$ then has two parameters, $\theta$ and $\phi$ due to using spherical coordinates.
Now the probability density function (pdf) of $f$ is given in here. From that we calculate its cumulative distribution function (cdf) which is also given in here.
We need the cdf because of the following property: If $T$ is a continuous random variable with cumulative distribution function $cdf(T)$, then the random variable $Y = cdf(T)$ has a uniform distribution on [0, 1]. Meaning that the inverse of $cdf(T)$ whose input would come from a uniform distribution on [0, 1] would output a value according to the probability distribution of $T$. More on that on wikipedia as suggested by peter.
Coming back to our point. We are looking for function closer to the actual probability distribution of $f$. Why ? Because it is the actual probability distribution that "causes" $f$ to output $x_i$ (sample generated in accordance with the probability distribution of random variable as stated in importance sampling section), since the probability distribution of $f$ takes two parameters, the angles of spherical coordinate, we need to adapt pdf and cdf to deal with two parameters.
It turns out, we can sample one of the angles uniformly, hence
phi = X_i.x * 2 * PI given that X_i.x is a random number between [0,1].
Then we only need to deal with pdf and cdf of a single parameter.
If we can find out manually the inverse of cdf of that parameter, we can map a random value between [0, 1] to a value according to correct distribution of that parameter. The piece of code that deal with X_i.y represents that manually found inverse of cdf.
You can read about the involved mathematics from:

Here

Here. This gives the full derivation

and Here


World coordinates, Normalised device coordinates and device coordinates

World coordinates are at the bedrock of your game world. The 3D positions of all objects are ultimately specified in world space, either directly or through a node hierarchy. The ground, buildings, trees, and other stationary things will be fixed in world space. Gameplay calculations and physics will be done in world space (perhaps with some local recentering for precision purposes, if the world is large). The axes of world coordinates are often taken to represent compass directions, such as X = east, Y = north, Z = up.

Local coordinates are attached to an object (there is one local coordinate space for each object). The axes may represent something meaningful to the object, for instance X = forward, Y = left, Z = up for an object such as a character, vehicle, gun, etc. that has an inherent orientation. As the object moves around, the relationship between local and world space (as expressed by a transformation matrix) will change. For instance if you flip your car upside down, its local Z axis ("up" in local space) will now be pointing "down" in world space.

Camera or view coordinates are local coordinates attached to the camera. This is still a 3D coordinate system, without any projection or anything, but with the axes aligned to match the screen orientation: usually X = right, Y = up, Z = backward. The transformation from world to view space is often known as the "view matrix".

Clip space coordinates are the coordinates output by a vertex shader: coordinates to which the projection matrix has been applied, but not the perspective divide. This is a 4D (homogeneous) space. (World, local, and view space are 3D with an implicit w = 1.) It's so named because this is the space in which view frustum clipping and culling takes place (conceptually, at least).

Normalized device coordinates, also commonly known as "screen space" although that term is a little loose, are what you get after applying the perspective divide to clip space coordinates. The 3D coordinates now represent the 2D positions of points on screen, with X and Y in [−1, 1], together with the depth within the depth buffer range, Z in [0, 1] for D3D or [−1, 1] for OpenGL. The axis orientation is X = right, Y = up, and Z can be either forward or backward depending on the depth buffer configuration.

Device coordinates are 2D pixel coordinates within the render target, with (0, 0) in the upper-left corner, X = right, Y = down. One arrives at device coordinates by applying the viewport transform to the normalized device coordinates; the viewport controls at what pixel offset and resolution the image appears in the render target. One thing that's important to note is that device coordinates are not integer values; they aren't indices of pixels. They're a continuous 2D coordinate system, using pixel-sized units, but in which fractional (subpixel) values are perfectly valid. Pixel centers lie at 0.5 offsets (like 0.5, 1.5, 2.5, ...) in this coordinate system.


For the lightning, I recommend using a midpoint displacement algorithm. You start with a line segment between any 2 points A and B (this works in either 2D or 3D). Calculate the midpoint of the segment AB. Now move that point a random amount in the direction perpendicular to the line segment AB and call it point C. Replace the original segment AB with 2 line segments AC and CB. Repeat with the new line segments until you are satisfied with how it looks. In my experience you'll need somewhere between 5 and 10 iterations depending on the resolution you're working at. You can add a glow to the resulting geometry to make it look more realistic.

For sparks you need to make a particle system. These are pretty fun. You generally have a point (or other shape) that generates particles every few milliseconds. That is, particles are spontaneously formed on the point or boundary of the shape periodically. Each particle has various properties such as location, velocity, color, size, lifetime, etc. When animating, you can update the velocity and position based on various forces such as gravity, wind, or some sort of convection current or flow field. After a particle has reached its "lifetime" it can disappear, or fade out, and a new particle will be generated in its place at the point or on the boundary of the shape. The particles themselves can be drawn as single pixels, small shapes, or images. There's an interesting wikipedia article on particle systems.


The box.obj file has no vertex normals, and by default Mitsuba will generate smooth normals for OBJ files that don't specify their own normals. This creates the magnification effect: the box with smooth normals forms a convex lens!

By adding this line to the box object in the scene file:

<boolean name="faceNormals" value="true"/>

I got results that look equivalent to yours:


UV mapping assigns fixed positions in the texture to each vertex. You'll have to distort the image so that it maps to a constant area in "model size" so to speak. I found the following other question that will probably give you a solution: https://blender.stackexchange.com/questions/22872/unwrap-uv-map-proportionally-to-object-size

It won't give disjointed parts like in your image, but I think it will still provide what you're looking for.


projecting 3d onto 2d plane is not showing up correctly

Your translate_coords_to_center(0,0) function should be done after you have done the projection from 3d to 2d.
e.g.
for i in range(N):
        x2 = cx + radius * math.cos(PI* angle * i / 180)
        y2 = cy + radius * math.sin(PI* angle * i / 180)

        # Do a 2d projection
        x2 = K* x2 / cz
        y2 = K* y2 / cz

        sx, sy = translate_coords_to_center(x2, y2)

        pixel(sy, sy)

while running:
    circle(0, 0, cz, 50)
    cz += 1
    pygame.display.flip()
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False

The effect you are seeing is caused by 2 things.

Because a constant x/y offset is incorrectly added in 3d world space the array of circles will appear to the right and below of the camera's position.
Because the 0,0 origin of view/camera space is not the centre of the display, instead it's at the top-left (Adding the window centre coords AFTER projection fixes that)


What is the reason that a boolean operation between two almost identical 3D triangular meshes fails?

Without debugging the routines with the offending data, it's impossible to say conclusively. However given the examples you have given it sounds like the libraries are falling foul of the coplanar triangle cases when evaluating the meshes. This is certainly an Achilles heel of many CSG implementations out there (far to many in my experience). In many cases, both theory and implementations seem to fail to address these coplanar cases in many CSG libs.

Certainly the fact that it crashes or fails to produce the desired result implies that the libraries you use are buggy or incomplete in this regard.

  The libraries would operate normally in any other case of Boolean operations between two different meshes.

When there are no coplanar triangles between manifolds, the result is well defined and since CSG operations use triangle/triangle intersections fundamentally, if the case of coplanar triangle intersections is not accounted for, everything should be fine and dandy when there are no coplanar triangle tests between meshes. Of course this is not what you require, but makes sense that this operates normally in these scenarios.

  The 2 almost-identical meshes I wanted to use were clean and water-tight. I even used an example shape which I subtracted from itself. It crashed too.

This is certainly not correct and I would say a bug. A crash implies something is very wrong and either the coplanar triangle support is buggy, limited, or simply not implemented at all.

I suspect the latter unless you can get a result out with something simple such as diff'ing one cube from another... if you can get a result in some sandboxed cases, then it could be any of the former possibilities.

To see if your CSG routines can handle the coplanar cases, subtracting one cube from another should produce nothing (which as you said it fails). You could also attempt to subtract one cube from another which has a single face coplanar. e.g. Two cubes, both of length length. One positioned at x - length should produce a cube which is a closed manifold. i.e. Cube A (1x1x1) positioned at 0,0,0... and cube B (1x1x1) positioned at 1,0,0.

  These 2 similar meshes are originating from the same mesh, but with a small difference, like a dent or some noise. I tried to introduce stronger noise or bigger changes, but it still would not work if the 2 meshes had their biggest part identical.

Probably just one triangle in a mesh is coplanar to another triangle in the other mesh which you are operating on, and thus failing. This may be a way to test by creating two meshes with only one triangle from each mesh being coplanar to another in the other mesh. If this fails then It might indicate that this corner case of coplanar intersection is not catered for at all. As said above, if it does work, then this may point to buggy behaviour at much higher level when dealing with coplanar triangles.

  By enlarging 1 of the 2 meshes by more than 0.1%, or by transforming/moving(see attached picture) it by a similar amount, the algorithm starts working. But by that point, the shapes are too different to be useful.

I presume these libs work on floating point numbers, in which case this sounds like the FP equality tolerances are kicking in and evaluating 'almost' coplanar triangles as coplanar, and then subsequently failing to result in the CSG operation. Of course this is problematic for you if the tolerances are high since the input geometry is expected to be 'different enough' to not fall within these tolerances.

Unfortunately, there is nothing you can really do which I think you would find satisfactory since you are effectively at the mercy of the accuracy of the implementation.

If you suspect the libs can handle coplanar triangles between meshes, and its just failing in many cases, you might need to find out how you can kludge the input geometry to not fall foul of the shortcomings that may or may not be present for these corner cases in the implementations:

Remeshing the input geometry may help, although this may result is varying quality of input data in comparison to your measured data, which kind of defeats the purpose and will result in loss of accuracy wrt the original input meshes, but may at least give a result.
Re-tessellating triangles while retaining the mesh features may help. This involves inserting Steiner points into a given triangle and re-tessellating into more triangles. This will result in larger meshes (mem footprint, not geometrically speaking) but may help the routines to resolve coplanar intersections. Unlike the previous suggestion, It will retain the accuracy of the topology of the original meshes though.
Change the FP tolerance (if applicable) that the routines use which may allow you to keep meshes a similar size. However the problem is still there, but this may be satisfactory for your needs allowing you to not have to scale as much. Certainly from a visual standpoint this would be satisfactory I would have thought, however from a numerical point of view this may not be. However this does not fix the issue, rather attempt to circumvent it.

Given that the majority of the meshes will perhaps be similar, if the vertices are what are similar (i.e coplanar triangles have the SAME vertex values, although that screen shot you posted suggests otherwise) you could perhaps simply diff the meshes yourself, finding all the vertices that are not equal when comparing and then tessellating these. This is far more complex though as you would need to keep track of the edges around the different vertices (these edges will use vertices which are equal to others in the other mesh, so not a quick and dirty fix, rather a slightly different approach avoiding the libs you use) so that you can tessellate the different vertices into a meaningful manifold. I wouldn't recommend this as it is quite a specialised case, and may not be generalised enough for your needs.

Ultimately, if you can't circumvent the issue present with the current libs. You may have to find another. I certainly wasn't able to find many that matched my needs with a correct lic that was applicable to work, however I think there are certainly libs out there that would cater for this. Look for CAD libs as this is an important operation in many CAD programs.

A few years ago I spent some time testing and evaluating various different C++ CSG libs for our needs at work as we required a fairly robust implementation for production. While there were various success, they all failed at some point, and coplanar triangles were usually the culprit. It's quite frustrating as the answer is actually quite straight forward and simply requires the Boolean operation to be performed in 2D between the coplanar triangles and then testing if the resulting split triangles are within the 'other' triangle. Why this rather crucial point is missed in a lot of implementations and theory is beyond me.
/rant


Can I raytrace gravitational lensing using only a single point source of gravity?

As you have probably observed, in general, you cannot solve in closed-form the equations for ray trajectories around multiple distorting objects.

The standard approach for such things is just to discretize it.  This looks like Eulerian integration.  Just step your ray a little bit toward the object, calculate gravity from all sources and bend it, then step it more, und so weiter.  This is called raymarching.

Since you're interested in problems on an astrophysical scale, you can assume the ray is linear far enough away from the object, with minimal error.

Here's a scene I rendered a while ago using a homebrew photonmapper (no QMC) using this raymarching technique.

I present: relativistic photon mapping: black hole in a Cornell box!


Correctness or otherwise of the procedure for generating an orbiting camera in WebGl

Your math looks correct to me. Your terminology is a little off - technically what you are creating here is just the 'view' matrix rather than the 'modelView'. If you're just drawing a single sphere at the origin then it doesn't make a difference, but normally the modelView is unique for each object in a scene - it's the object's model-to-world transform matrix multiplied by the view matrix you describe above.


This looks like the VkFormat being selected for the VkSurfaceFormat, most likekly VK_FORMAT_R8G8B8A8_SRGB. Try VK_FORMAT_R8G8B8A8_UNORM for the VkFormat with a corresponding colorspace (like VK_COLOR_SPACE_SRGB_NONLINEAR_KHR). That should give the dark color you are expecting.


We tackle this problem by computing, like shaders do, the color of a certain pixel with position $y$ at timestamp $t$, namely $I_{x, y}(t)$. All we want to know if that pixel must draw a bar on timestamp $t$. So simple as that, we start with

float horizon = 0.5;

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
   // Normalized pixel coordinates (from 0 to 1)
   vec2 uv = fragCoord / iResolution.xy;

   float position = uv.y;
   float time = iTime;

   // everything above our horizon we make white
   if (position > horizon - 0.01) {
      fragColor = vec4(1.0);
      return;
   }

   // check if pixel is on the bar, and if so color it white.
   float color = isOnBar(time, position);

   fragColor = vec4(color, color, color, 1.0);
}

Next step is to determine whether that pixel actually is on the bar.
This problem can be solved with the rules of gravity. We visualize the bars falling down from the horizon. For any pixel $I_{x, y}$, we can say the pixel is on the first bar if

$h_{bar}(t) = I_{y}(t)$.

We know however that bars fall with gravity from the horizon, so:

$h_{bar}(t) = h_{horizon} - gt^{2}$.

Hence, the falling time $ t_{\text{fall}} $, when the bar hits the pixel, can be calculated with

$ t_{\text{fall}} = \sqrt{ \frac{h_{horizon} - I_{y}(t_{\text{fall}})}g } $.

However, this was just our first bar. Any other time a bar hitting that pixel, is when $
(t - t_{\text{fall}}) \bmod a = 0 $, where $ a $ is our falling frequency.

e.g. first bar hits a pixel after 3 seconds, and every 2 seconds another bar falls, then for instance on the 7th second a bar again hits our pixel: $ (7 - 3) \bmod 2 $ is indeed $0$.

How much room we give when computing that modulo close to zero, is the width of our bar. Right now, we can incorporate the widening of the bars when they drop.

float gravity = 0.1;
float timeBetweenBars = 0.5;

float isOnBar(float time, float position) {
   float timeOfHit = sqrt((horizon - position) / gravity);
   float moddedTime = mod(time - timeOfHit, timeBetweenBars);

   // the lower the pixel, the more room we give to catching the bar with this pixel.
   if (moddedTime < 0.1 * (1.0 - position)) {
      return 1.0;
   };
   return 0.0;
}

Here is the code on ShaderToy


Why do games only allow certain resolutions?

Note that in your Dark Souls example, all the resolutions listed are the exact same aspect ratio.  So the question is, why doesn't it support arbitrary aspect ratios?

The answer is that it costs time and effort to support arbitrary aspect ratios, and if that money spent doesn't translate to increased profit, then it doesn't make sense to do.

This answer covers the 2d layout issues, which includes extra time and effort to create a versatile layout that can be resized so that elements do not overlap, stretch, or pixelate in an undesired manner, and that the overal UI maintains the best user experience.
The developer would have to ensure that the window doesn't get too large or too small.  Does it make the game any better if the user can drag the window size to be smaller than 800 x 450?  probably not.
you have to define the behavior of the camera when the aspect ratio changes.

For example, take the following scene:

If we allow the user to change the horizontal size, we can decide to "zoom out" horizontally so that the same area still shows up horzontally but now we can see more vertically.

Or, instead, we can keep the same zoom level but just clip the viewing area more:

Note, that if you don't specifically handle the changing aspect ratio, the default graphics pipeline behavior will be to stretch the image:

These decisions take time and effort to design, code and test, but management has to decide if it actually makes the user experience better enough to be worth the time.  Apparently, most games have decided its not worth the effort.  I'd even go further and suggest that allowing users to screw up the best layout of the game is a bad feature.


Vectorizing to polylines instead of paths

If you look at chapter 8 of the SVG Specification, it describes how to parse a path element. The short version is that you'll want to find the d attribute of the element. That element should be a string describing the curve. It will contain the following commands:

m - Moveto command
l - Lineto command
c, s, q, & a - Curve commands
z - close path

Simple lines will just be combinations of m and l sequences. The curves will be either circular arcs, quadratic bezier or cubic beziers.

To convert a Bezier curve into line segments, you can do the following:

const float kNumPolyLines = 10; // <- change this to make more lines per bezier segment
const float kDeltaT = 1.0 / kNumPolyLines;
Point c0;
Point c1;
Point c2;
Point c3;
GetBezierControlPoints(&c0, &c1, &c2, &c3); // <- get the 3 or 4 control points for the next bezier segment
Point prevPoint = c0;
Point currentPoint;
for (float t = kDeltaT; t <= 1.0; t += kDeltaT)
{
    currentPoint = CalculateBezier(c0, c1, c2, c3, t);

    // Now you have the next line segment in your poly line.
    // It goes from prevPoint to currentPoint. Store it or draw it,
    //or whatever you need

    // Get ready for the next iteration
    prevPoint = currentPoint;
}

You can find the formula for both quadratic and cubic bezier curves on Wikipedia.


Environment map importance sampling (beginner)

From all other implementations I have seen, you generally sample twice, once in the Y dimension to select a row and then once in the X dimension to selct the pixel from this row.
I am not sure if this is necissary, maybe you're one dimensional sampling is fine.
I think you also need to multiply your pdf so that it is relative to solid angle (so it's total sums to 4pi).
You should be able to do this by dividing by (2 * Pi * Pi * sinTheta).
PBRT is probably the best textbook to help you. You can find it online here:
http://www.pbr-book.org/
The source code is also available. Here is the code for the environment light:
https://github.com/mmp/pbrt-v3/blob/aaa552a4b9cbf9dccb71450f47b268e0ed6370e2/src/lights/infinite.cpp
and the code for the Distribution2D sampler:
https://github.com/mmp/pbrt-v3/blob/aaa552a4b9cbf9dccb71450f47b268e0ed6370e2/src/core/sampling.h


Is there a way to compile part of a shader into a spirv?

What you want is generally not possible for OpenGL's SPIR-V support. SPIR-V compiles to shader objects, and shader objects theoretically allow you to link multiple objects of the same stage in a program. However, OpenGL's usage of SPIR-V requires that each SPIR-V shader be a complete shader stage with its entrypoint function.
So the requirement of shaderc mirrors the requirements of OpenGL.
The best you could do is combine multiple SPIR-V binaries into a single binary before handing it off to OpenGL. But that's a non-trivial process, since you basically have to renumber every reference in at least one of them, and combine any references that map to the same type.


Sharp Corners with Signed Distance Fields Fonts

EDIT: Please see my other answer with a concrete solution.

I have actually solved this exact problem over a year ago for my master's thesis. In the Valve paper, they show that you can AND two distance fields to achieve this, which works as long as you only have one convex corner. For concave corners, you also need the OR operation. This guy actually developed some obscure system to switch between the two operations using four texture channels.

However, there is a much simpler operation that can facilitate both AND and OR depending on the situation, and this is the principal idea of my thesis: the median of three. So basically, you use exactly three channels (ideal for RGB), which are completely interchangeable, and combine them using the median operation (choose the middle value out of the three).

To accomodate anti-aliasing, we don't work with just booleans, but floating point values, and the AND operation becomes the minimum, and the OR becomes the maximum of two values. The median of three can indeed do both: if a < b, for (a, a, b), the median is the minimum, and for (a, b, b), it is the maximum.

The rendering process is still extremely simple. The entire fragment shader including anti-aliasing can look something like this:

int main() {
    // Bilinear sampling of the distance field
    vec3 s = texture2D(sdf, p).rgb;
    // Acquire the signed distance
    float d = median(s.r, s.g, s.b) - 0.5;
    // Weight between inside and outside (anti-aliasing)
    float w = clamp(d/fwidth(d) + 0.5, 0.0, 1.0);
    // Combining the background and foreground color
    gl_FragColor = mix(outsideColor, insideColor, w);
}

So the only difference from the original method is computing the median right after sampling the texture. You will have to implement the median function though, which can be done with just 4 min/max operations.

Now of course, the question is, how do I build such a three-channel distance field? And this is the tricky part. The most obvious approach that I took in the beginning was to perform a decomposition of the input shape/glyph into three components, and then generate a conventional distance field out of each. The rules for this decomposition aren't that complicated. Firstly , the area with at least 2 out of 3 channels on is the inside. Then, if you imagine this as the RGB color channels, convex corners must be made of a secondary color, and its two primary components continue outward. Concave corners are the inverse: Two secondary colors enclose their common primary color, and the wedge between where both edges continue inward is white. I also found that some padding is necessary where two primary or two secondary colors would otherwise touch to avoid artifacts (for example, in the middle stroke of the "N" in the picture).

The following image is an example decomposition generated by the program from my thesis:

This approach however has some drawbacks. One of them is that the special effects, such as outlines and shadows will no longer work correctly. Fortunatelly, I also came up with a second, much more elegant method, which generates the distance fields directly, and even supports all of the graphical effects. It is also included in my thesis and so is also over a year old. I am not going to give any more details right now, because I am currently writing a paper that describes this second technique in detail, but I will post it here as soon as it's finished.

Anyway, here is an example of the difference in quality. The texture resolution is the same in each image, but the left one uses a regular texture, the middle one uses an ordinary distance field, and the right one uses my three-channel distance field. The performance overhead is only the difference between sampling an RGB texture versus a monochrome one.


Difference between isnan and comparison to self

IEEE-754 defines the isNaN property, which determines if a value is NaN. However, it is also defined to be a non-signaling operation; calling isNaN on an sNaN will never result in a signaled failure. By contrast, equality testing may signal when given sNaN values.

Now, GLSL explicitly states that signaling doesn't matter:

  Support for signaling NaNs is not required and exceptions are never raised.

But of course, there's also the fact that isnan is a vector function, while !(x == x) only works on scalars. So if you want to test if any of the elements of a vec3 are NaN, you can do any(isnan(some_vec)), while the equality-testing equivalent would be even more obtuse and verbose: any(not(equal(some_vec, some_vec))).

And all of that ignores the simple usability reasons for having isnan. !(x == x) is an idiom which must be learned, while isnan is an easily looked up function whose name tells you want it does. After all, most users of floats tend to forget about NaN, so such an expression seems like a tautology, while "isnan" explicitly reminds such people that NaNs are a real thing and they are reading code where that fact is important.

Also, as previously stated, isNaN is a property defined by IEEE-754. It is also a function provided by the standard libraries of C99, C++11, Java, C#, Swift, Go, JavaScript, and I can keep going, but I think my point is clear.


Does premultiplied alpha give order independent transparency?

Premultiplied alpha itself does not give you order independent transparency, no.

This page talks about how it can be used as part of an order independent transparency solution however: http://casual-effects.blogspot.com/2015/03/implemented-weighted-blended-order.html

Other benefits of premultiplied alpha include:

better mipmaps for textures that contain alpha: https://developer.nvidia.com/content/alpha-blending-pre-or-not-pre
Being able to do both transparency and additive blending within a single texture without needing to change render states or issue multiple draw calls: http://amindforeverprogramming.blogspot.com/2013/07/why-alpha-premultiplied-colour-blending.html
Being able to represent more realistic lighting situations, such as the fact that being non opaque doesn't always mean transparency: http://distantsoulsdev.blogspot.com.ar/2013/02/premultiplied-alpha-online-interactive.html
Being able to more naturally / correctly / effortlessly composite images that have alpha channels.
Being able to get better results from image processing (blurs / convolution). No black halo: https://stackoverflow.com/questions/4854839/how-to-use-pre-multiplied-during-image-convolution-to-solve-alpha-bleed-problem
Premultiplied alpha also allows transparent objects to be drawn front to back instead of back to front, possibly saving pixel writes.


Fast approximation of spherical texture coordinates

Incidentally I implemented something similar for similar hardware some 20 years ago (:

IIRC, I calculated 1D look-up-table using acos for parallel projection of a cylinder and just scaled it for each scanline to reduce computation cost. The LUT just gives you x-coordinate offset to the texture and you add constant offset for each scanline for the rotation.

For y-axis the input texture can be preprocessed to have proper spherical distortion applied and you can just linearly scale the texture vertically. Though it's not that much extra computation to do it run-time either.

This doesn't give you proper perspective projection but might be enough for your purposes.

Edit: Didn't notice you said you can't afford LUT. You could do the table look-up every N pixels and linearly interpolate inbetween. This was a common method in the past to do affordable perspective correct texture mapping for software rasterizers to avoid div/pixel


This type of artifact is a tell-tale sign that you've messed up your texture image format at the byte level. It's hard to tell what exactly is wrong, but it's something along the lines of bad pitch or a component layout mismatch.

Have a read of this wiki page and try again. glTexImage2D should work, keep trying until it does.


Are the benefits of a Metal early depth test limited to one draw call

I don't know what Apple have done with Metal per se, but as the iOS devices have PowerVR-derived GPUs, normally I would think that doing as much opaque-only drawing "up front" would usually be beneficial, but if, as in your case, only 5% of the screen is covered by opaque geometry, I'm not sure how much the software overhead will affect things.

  The reason I am unsure is that I know the second you call discard on a fragment you remove the possibility of Metal preforming early depth test which is allegedly quite an optimization. What I dont know is if discards disable early depth test on draw call 1 do they also disable the optimization for draw call 2's geometry?

With PowerVR devices, a depth test can "always" be done as the first step in the pipeline, so any geometry that was submitted earlier in the scene can cull later things prior to fragment shading.

With alpha-tested (AKA "punch-though"), the Z-test can thus be done against previous geometry, but there may be a cost of then doing the texturing/shading and feeding that back to update the depths for subsequent geometry.

If you aren't doing alpha-test/discards, and are doing your own sorting, then this feedback is not required (but obviously there may be more shading/texturing work to perform).

(Such a pity, though, that we removed the OIT feature from the early PowerVR hardware)

  Does the fragment stage happen after each draw calls vertex phase or after all of them?

I'm not sure what you are asking, but for a given scene, all the vertex shading (across multiple calls) is performed**, the geometry buffered, and only then rendering commences on the entire scene.

**There are some caveats here but that probably doesn't affect you.

  Do you think I should explore sorting the triangles and using alpha blending on draw call 1? Again this would be a negligible CPU cost would it actually help the GPU?

Are you asking "should I profile/time the results and see what is best for my use case?".  Probably "yes".


When I implemented real-time area lighting, there were two documents I kept referring to: "Moving Frostbite to PBR" by Sebastien Lagarde and "Real Shading in Unreal Engine 4" by Brian Karis. Also "Lighting of Killzone Shadow Fall" by Michal Drobot is an interesting read on the topic. None of these documents deal with area shadows though but focus on energy conserving approximation of the lighting integral. While there are some code snippets, these documents lack quite a bit of details of how to implement area lighting though, so you need to do fair amount of research on your own. For example I posted a question and an answer here related to implementing rectangular area lights.
A very valuable tool I used was to implement a Monte Carlo reference model to easily compare different approximations to the ground truth. Some people also use Mitsuba renderer to validate their implementations. I think this ground truth reference is extremely important so that you know if there are mistakes in your implementation and also to see where/how these approximations fail and if you can somehow try to improve them. Implementing Monte Carlo integration is also quite a good exercise to understand the underlying problem.


Importance sampling GGX NDF - fireflies and bright final result

I've made a few changes in my code, so let's start with basic image that I had problem we that we'll state as a problem at the beginning:

After changing part of the code from:
float weight = abs(dot(wo, wm) / (dot(normalTS, wo) * dot(normalTS, wm)));
if (dot(normalTS, wi) > 0.0f && dot(wi, wm) > 0.0f)

To this (i.e. moving part of calculations to world-space):
float weight = abs(dot(wo, wm) / (dot(triangleNormal, incomingRayDirWS) * dot(normalTS, wm)));
if (dot(triangleNormal, rayDirWS) > 0.0f && dot(wi, wm) > 0.0f)

I am getting better results. Here are series of images with different settings:
1024 spp, 4 bounces

5000 spp, 4 bounces

5000 spp, 1 bounce

Notice that 4 bounces version is having more fireflies and noise. Long specular paths might be a reason of introducing noise. Also, 5000 spp are getting darker than 1000 spp. Reason is that path is probably introducing no changes, therefore average color is getting darker. You can see it in the last line of main post. I'm going to put it below for reference:
if (dot(triangleNormal, rayDirWS) > 0.0f && dot(wi, wm) > 0.0f)
    throughput *= F * G * weight;
else
    throughput *= 0;

Last thing. Here is image generated with better method ( https://schuttejoe.github.io/post/ggximportancesamplingpart2/ ). 1024 spp, 4 bounces:

Overall, I got rid of weird noise and artifacts. Algorithm got closer to better method and I think it can be treated as fixed.


Programmatically generating bin files for GLTF models

The .bin file that goes with a .gltf only includes mesh data (position, normal, uv, other vertex attributes as needed).  It does not normally include any images (unless you're using the base64 export "glTF Embedded", but don't use that for this, use "glTF Separate" instead).

If you look in the JSON text of the .gltf file, you'll see an images array, where image files are called out by uri (filename, if local, or relative URLs for the web).

Let's look at the Damaged Helmet sample model for an example.  The JPG images are called out in the .gltf file, separately from the .bin file.  The data from the .bin file is used by buffers, which in turn are used by bufferViews, accessors, and then mesh primitives.  It is exclusively vertex mesh data, no images.  This is a separate chain from how the JPGs are being used.  The JPGs are called out by images, which are referenced by textures which in turn are used by materials, which are referenced by mesh primitives.  For a more graphical version of this explanation, take a look at the Overview Chart for glTF.

So when exporting from Blender, choose glTF Separate from the combo box in the export options panel.  This will produce a .gltf file, a .bin file, and may copy several images (be sure to export to a different folder from your source files, so nothing gets overwritten).  To change the images, just update the uri fields in the images array inside of the .gltf file, and keep the same .bin file around (the one that came out of Blender).

That way, you can swap out images with other images, but keep the same mesh data.


If I have an nVidia graphics card on Linux, do I need Mesa or FreeGLUT or anything like that?

The Mesa3D project is a software implementation of OpenGL; that is, an implementation that does all the rendering on the CPU with no GPU involvement at all. It's designed in two layers: a backend which renders each frame from some internal state, and a frontend which implements the OpenGL API and sets up the state for the backend to use. Most OpenGL implementations have a similar two-layer architecture, so that the same backend can be used by an OpenGL frontend, an OpenCL frontend, a Direct3D frontend, etc.

The open-source NVidia drivers (and most open-source drivers) use Mesa's frontend to provide the OpenGL API and a backend specific to that GPU. The OpenGL API is quite complicated to implement, so it makes sense for vendors to share an existing front-end in this way. However, closed-source drivers do not use Mesa's frontend and have their own closed-source frontend.

If you install and use a closed-source GPU driver, it will include its own libgl, which doesn't use Mesa at all. You don't need to install anything else to run OpenGL applications. (You might need an extra "-dev" package to get the header files needed to compile OpenGL applications.)

If you install and use an open-source GPU driver, it will use Mesa, but the packages in your distribution should have the dependencies set up correctly, so that if an external Mesa package is needed, it will be installed automatically.

Either way, you shouldn't have to worry about whether the OpenGL implementation you're using is based on Mesa or not.


Very efficient/fast marching cube implementation

This seems to be eluding to a Marching Cubes LOD algorithm such as: Place the entire volume in one giant cube. Break that volume into NxNxN cubes. And continue doing so until the cubes are at the finest granularity needed for the density function about 8 level works. Each volume then responds with either: No voxels - stop processing that sub volume OR Yes there are voxels -- process the sub volumes and keeping doing this until:  the voxel level is reached. Finally set a LOD - such as equating a screen space distance range to each subvolume level and treat the subvolume for that LOD as the final voxel. Cracks are a major issue with this algorithm and why you won't see many implementation of it.

There is also Eric Lengyel's Transvoxel algorithm: https://transvoxel.org

and here is a link an implementation of it
https://github.com/stoyannk/voxels


Get screen size and then resize image accordingly

The OS does not, in many or most cases, know the physical screen size! Often the only thing you can glean is the make and model, sometimes not even that. So if you do not want to have a database then the best you can do is ask your user.

There are also some interesting corner cases. What do you do if the monitor is in fact a projector? What do you do in multi monitor setups? What if the image is halfway between 2 monitors? What do you do if the operator is using pixel zoom for accessibility reasons? What do you do if the image is shown on 2 monitors at the same time? What do you do if the monitor is huge, like one used in a conference room and your image turns out to be very small in comparison (like 4 by 4 pixels)? ...


The notation $\delta_{i,j}$ is the Kronecker delta, a notation commonly used in physics. It's defined as:
$$\delta_{i,j} \equiv \begin{cases}1,&i=j\\0,&i\neq j\end{cases}$$
So, as you suspected, it's essentially a shorthand for the identity matrix.

The notation $\mathbf{S}_i \in \mathrm{R}^s$ means that each $\mathbf{S}_i$ is an $s$-dimensional vector, so the index $j$ labels the components of this vector. The term "indicator" probably refers to the concept of an indicator function, which is a function that labels some set or point with the value 1 while mapping everything else to 0. So, the "$i$-th spring indicator", in this context, is a vector that has a 1 in the $i$-th component and 0 everywhere else. That's exactly the Kronecker delta, i.e. the identity matrix.

The way this "indicator vector" gets used in equation 12 in the paper is in the form of an outer product with another vector, $\mathbf{A}_i \in \mathrm{R}^m$. The two vectors are multiplied together in the form $\mathbf{A}_i \mathbf{S}_i^{\mathrm{T}}$, which is an outer product. The result is an $m \times s$ matrix with the contents of $\mathbf{A}_i$ in the $i$-th column, and zeroes everywhere else.


Procedural generation of biological models

I am playing with L-Systems and with some thorough understanding of how they work I think you could manage to get something useful out of them. You can see Paul Bourke's site for a sneak peek of what they do and you can read the  book The Algorithmic Beauty of Plants by Przemyslaw Prusinkiewicz and Aristid Lindenmayer. I am going to try and explain what L-Systems do in case the links die, L-Systems are a way of creating simple rules that define the topology of plants/organisms/fractals by simple axioms, this axioms can be rewritten in terms of other rules that can create incredibly organic behaviours, originally this systems were abstract in the sense that they lacked geometric meaning, but over the years people have figured out how to transform this topologies into geometry through interpreting the rules as "directions" much like turtle graphics. I have seen this systems replace basic geometric rules such as "draw a line between two points" into create a complex building with "A or C characteristics" so I think you could define your interpreter in terms of biologically accurate creatures.


Path tracer - multi layered materials and importance sampling

Wenzel Jakob et al presented a framework for layered Materials at SIGGRAPH 2014. Section 6.2 explains importance sampling. If you prefer code over equations, the method is implemented in the Mitsuba renderer.


Is there a way to do sphere caustics analytically?

You might find this paper useful—its authors describe a way to make reasonably plausible caustic effects just by modulating their shadows with the surface normal and refractive index. It won’t get you the more complex/interesting hot-spot shapes that you’d get with photon mapping, but it could work well depending on your scene.

  To extend stochastic shadow maps to mimic caustics, we decrease shadowing where light strikes a transmissive surface at normal incidence. We increase shadowing where the refractor is at a glancing angle and thus has a low Fresnel transmission factor.
  Specifically, during shadow map rendering, we increase ρ(λ) where the magnitude of the dot product of the surface normal nˆ and the vector to the light ωˆi is close to 1, and decrease ρ(λ) where the dot product is close to zero. We also increase this effect with the relative refractive index η.


Heat Method (Crane et Al) How do we pick u?

From my reading, u seems to be a suitable function that approximates heat flow on the surface of a manifold?

$u$ is the function that describes how your quantity behaves/evolves in a certain field. In the paper, the quantity is the temperature or heat flux, I guess. However, most of the time there is no analytical solution/formula for $u$. This is where methods like Finite Elements (FEM) come into play. By discretizing your field, you can piecewise approximate your function $u$.
In your case, you would use your mesh, which is already a discretization of your surface. Your elements are the triangles and you need to define how the nodal quantities are interpolated inside each triangle. --- Here, linear interpolation is probably the way to go. Otherwise, you need to remesh your geometry or introduce additional nodes for higher-order approximations.
Then you have to assign to each Node/Vertex an initial value $u_0$ as written in the answer of gilgamec. Afterwards, you build and solve your finite element system and get the nodal distribution of $u$ that actually solves your equation or system of equations. The finer your mesh gets, the better the solution. Higher-order interpolations will also help with accuracy.
So $u$ or its nodal values are what you are actually looking for as lightxbulb said in his comment. It's your unknown quantity.
If this doesn't help, you might want to read some literature about the finite element method. Can't tell how helpful the following links are, but a short glimpse looked promising. You will see, that they use $u$ all over the place. So I hope one of them will help you:

A gentle introduction to the Finite Element Method
PE281 Finite Element Method Course Notes
Introduction to the finite element method
finite element analysis by hand

I also had a link to a good online tutorial similar to the last link I provided that helped me a lot in understanding the fundamentals. If I find the link, I ll add it to my answer.

Found the link I was referring to. Unfortunately, it is in german:

FEM Handrechnung


Is it possible to render red / blue 3d from one image and a depth buffer?

Yes, this is essentially the same problem that occours in paralax mapping. What you basically have is a colored height field. That needs to be rendered from a second view.

There are several ways in which you can approach this. You could cheat and just shift pixels by their depth and ignore occlusion. Or you could just render a textured height map as if every pixel would be opaque. Or you could use any of the available paralax mapping methods to eliminate occlusion.

How well this works depends on the. A scene multiple transparent surfaces wouldnt work very well without some deep map.


Enter/Exit of unit cube representable as 4 numbers?

3D object surface is 2D domain and can thus be parametrized with a 2D function. In case of a unit cube you could parametrize the surface for example by mapping each 6 faces to a 2D image strip. So if you know two points of the ray that intersects the cube, you can define it with 4D function using the parametrization