Dealing with blown-out highlights in PBR I've been writing some custom shaders based around Unity's Standard BRDF, which is uses GGX distribution term and Smith visibility term. It looks great at either end of the smoothness scale (1 and 2 in the pic below). However, with the smoothness set to about 50%, things start to get a little ugly (3). Increasing the metalness makes it worse (4), and strongly coloured specular gets worse still, as different colour channels blow out at different rates(5).

This seems to be caused by specular term returning >1 around the highlights, probably caused by Unity's hack of multiplying specular by pi (instead of dividing diffuse) as shown in the code snippet below.

// Specular term
// HACK: theoretically we should divide diffuseTerm by Pi and not multiply specularTerm!
// BUT 1) that will make shader look significantly darker than Legacy ones
// and 2) on engine side "Non-important" lights have to be divided by Pi too in cases when they are injected into ambient SH
half roughness = PerceptualRoughnessToRoughness(perceptualRoughness);
half V = SmithJointGGXVisibilityTerm (nl, nv, roughness);
half D = GGXTerm (nh, roughness);
half specularTerm = V*D * UNITY_PI; // Torrance-Sparrow model, Fresnel is applied later

HDR rendering helps a lot, but unfortunately our project is on mobile and using LDR. I'm making a car paint shader with a clear coat, so the base coat will likely be semi-rough and semi-metallic in a lot of cases which seems to be a worst-case-scenario for this BRDF. Can anyone suggest a way to get saner specular highlights? Either some kind of rolloff for out-of-range specular values, or a different BRDF altogether?

Using physically based BRDFs only makes sense if your entire pipeline is built for physical units - the extreme range of values can't be displayed properly without some form of tone mapping. You didn't include that part of the code but from the looks of it I'd say you're doing a simple clamp() followed by linear->sRGB conversion, which causes the bad highlights.

Since you're not ready to run your entire frame in HDR, you could apply something simple like a sigmoid operator to just the specular component. This may look a bit off next to other BRDFs but should make those highlights much more pleasant.


Equation for camera projection I am trying to find an equation for a camera projection.
The goal is to map the point $P(P_x, P_y, P_z)$ from the world coordinates onto the window coordinates $Q(Q_x, Q_y)$.

The eye of viewing is located a distance D from the $\hat{e_1}$-$\hat{e_2}$ plane, and in the direction of -$\hat{e_3}$.
$\hat{e_1}$, $\hat{e_2}$ and $\hat{e_3}$ are unit vectors.

What is the equation to map P from world coordinates to window coordinates?
$\vec{Q}(Q_x, Q_y) = f(\vec{P}, \hat{e_1}, \hat{e_2}, \hat{e_3}, \vec{v}, \vec{c})$
with all the vectors on the right side described in terms of world coordinates, eg

$\vec{P} = (P_x, P_y, P_z)$

$\vec{C} = (C_x, C_y, C_z)$

$\vec{e_1} = ({e_1}_x, {e_1}_y, {e_1}_z)$

$\vec{e_2} = ({e_2}_x, {e_2}_y, {e_2}_z)$

$\vec{e_3} = ({e_3}_x, {e_3}_y, {e_3}_z)$

This is pretty straightforward matrix math and can be found on Wikipedia. It's quite elaborate, though, if you want to learn the math behind it and the other types of projection that exist. I will encourage you to checkout either Real-Time Rendering or Computer Graphics, Principles and Practice.


Why is the symbol for solid angle a small omega in the definition of the BRDF? In many resources the BRDF is defined like $$f_{r}(\omega_{i},\omega_{r}) = \frac{\mathrm{d}L_{r}(\omega_{r})}{L_{i}(\omega_{i})\,\cos{\theta_{i}}\,\mathrm{d}\omega_{i}}.$$ The index $r$ might be replaced by $o$ but overall the equation looks something like this. However, in some other resources I see the definition $$f_{r}(\theta_{i},\phi_{i};\theta_{r},\phi_{r}) \equiv \frac{\mathrm{d}L_{r}(\theta_{i},\phi_{i};\theta_{r},\phi_{r};E_{i})}{L_{i}(\theta_{i},\phi_{i})\,\cos{\theta_{i}}\,\mathrm{d}\omega_{i}}.$$ The second definition seems more precise to me as it – besides other things – clearly states that the BRDF in its purest form is a 4-dimensional function of two directions and clearly differentiates between the direction $(\theta_{i},\phi_{i})$ and the solid angle $\omega_{i}$.
It seems to me as if the first definition mixes up the direction of incidence and the solid angle in the direction of incidence. The symbol $\omega_{i}$ stands for the direction – at least I hope so as the BRDF is not a function of solid angles – to economize on writing, I understand that, but what is $\mathrm{d}\omega_{i}$ then? It must be the solid angle but $\omega_{i}$ is a direction, not a solid angle. Is $\mathrm{d}\omega_{i}$ simply a common way to express the solid angle in a specific direction? May I have misunderstood something?
Another thing I've noticed is that the solid angle in the definitions of radiance, luminance and basically all other radiometric and photometric quantities that deal with a solid angle is symbolized with $\Omega$, i.e. a big omega. Actually, as far as I can tell, solid angles are usually symbolized with $\Omega$. Why is it different in the definition of the BRDF? Am I missing something?

The use of spherical coordinates like $(\theta, \phi)$ is kind of an implementation detail that isn't an essential part of the definition of a BRDF, or other functions defined over a spherical domain generally.
It's a lot like how in vector math, we will often express a point or vector by a single symbol like $v$, understanding that "under the hood" it's going to be expressed as coordinates $(v_x, v_y, v_z)$ in some coordinate system. Writing out all the coordinates in every equation would make notation very verbose and would often obscure the vectorial relationships we really care about.
Similarly, with functions defined over a sphere or hemisphere, we will often express a direction by $\omega$, understanding that when we get to actually writing code this will stand in for some coordinates on the sphere, but the exact choice of representation is not relevant to the math.
The symbol $\mathrm{d}\omega$ represents a differential quantity of solid angle—this is an abuse of notation for sure, but it serves the same purpose of abstracting away the implementation details; in spherical coordinates that would be $\sin(\theta) \, \mathrm{d}\theta \, \mathrm{d}\phi$ for example, but it's kind of assumed as background knowledge that you know how to integrate things over a sphere in whatever coordinate representation you're using.
Differentials as seen in calculus sort of aren't actual values themselves, but are things that can be either integrated over (in which case they specify the measure of the integration), or you can take a ratio of two differentials, which defines a density function (kind of a generalized derivative).
The symbol $\Omega$ is most often used to represent the hemisphere in my experience, not solid angles generically. You'll see something like $\int_\Omega \text{...stuff...}\, \mathrm{d}\omega$ to mean integrate that function over the normal-facing hemisphere with respect to the solid angle measure (i.e. area on the unit sphere).


How do I convert a hit on an infinite plane to UV coordinates for texturing in a Raytracer? I am trying to convert a ray hit on an infinite plane, defined by an origin and a normal vector, into UV coordinates, so I can find the appropriate texel at that point.

Code I have is close to functional, however, when the plane has a normal of (0, 0, -1) the result is a plain white plane.

The code that calculates the UV coordinates is pasted below:

// Calculate UV coordinates for the texture
vec3 u = vec3( normal.y, -normal.x, 0 ).normalized();
vec3 v = normal.cross( u );

h.u = u.dot( h.coordinates );
h.v = v.dot( h.coordinates );

h is a Hit struct that stores relevant data from a ray-primitive intersection, coordinates is the vec3 where the intersection occurred.

The current situation is shown in the following image; The floor plane properly shows the texture, however, the back wall does not.

However, when I change the normal of the back wall ever so slightly, the texture appears

In image 1 the normal of the back plane was (0, 0, -1), while it was (0.0001, 0, -1) in the second image.

I have been trying to figure this one out, but I am hitting a dead end.

If your plane has a normal of $\begin{pmatrix}0 & 0 & z\end{pmatrix}^T$, then your computation

vec3 u = vec3( normal.y, -normal.x, 0 ).normalized();
vec3 v = normal.cross( u );

will result in u and v both being $\begin{pmatrix}0 & 0 & 0\end{pmatrix}^T$.

A more general approach would be, for example, to compute the cross product of your normal with each of the three base vectors and use whichever result happens to be the longest as your direction for u:

vec3 computePrimaryTexDir(vec3 normal)
{
    vec3 a = cross(normal, vec3(1, 0, 0));
    vec3 b = cross(normal, vec3(0, 1, 0));

    vec3 max_ab = dot(a, a) < dot(b, b) ? b : a;

    vec3 c = cross(normal, vec3(0, 0, 1));

    return normalize(dot(max_ab, max_ab) < dot(c, c) ? c : max_ab);
}

  …

    vec3 u = computePrimaryTexDir(normal);
    vec3 v = cross(n, u);


Triangulation of vertices of an ellipsoid I have the parametric model of an ellipsoid like this:

\begin{cases}
x=a_1\,cos \,\eta \:cos\,\omega  & -\pi\leq\omega\leq\pi \\ = _2\,cos\,\eta\:sin\,\omega \\ z=a_3\,sin \,\eta & -\frac{\pi}{2}\leq\eta\leq\frac{\pi}{2}
\end{cases}

and

\begin{cases}
n_x=\frac{1}{a_1}\,cos \,\eta \:cos\,\omega  & -\pi\leq\omega\leq\pi \\n_ = \frac{1}{a_2}\,cos\,\eta\:sin\,\omega \\ n_z=\frac{1}{a_3}\,sin \,\eta & -\frac{\pi}{2}\leq\eta\leq\frac{\pi}{2}
\end{cases}

First of all I do not understand the last system of equation what it represents.
 Shouldn't  parameterisation of an ellipsoid be the first system alone?

I have to create a polygon mesh (in webgl) from parametric model of ellipsoid. I can choose dimensions of semiaxes (namely values of $a_1$ $a_2$ and  $a_3$) and the sampling step of 2 angles (the parameters).  In this way I'm able to create in javascript with a little of code several vertices of ellipsoid. But I can't imagine how to realise the triangulation. Namely given the vertex buffer, I have to create the index buffer but in what way do I choose the indexes of the vertices to tessellate in a correct way ellipsoid surface? Is there a clever/smart way to proceed?

EDIT:

I implemented  that myself, in javascript (the code works):

 function degToRad(degrees)
     {
        return degrees * Math.PI / 180.0;
     }

     function calculateVNI()
     {
     var indices = [];
     var vertices = [];
     var normals=[]
     var a1 =10;
     var a2 =5;
     var a3 =7;
     //var pi = Math.PI;

     var step_eta   = 4; //step in degree (choose a divisor of 360)
         var step_omega = 2; //step in degree (choose a divisor of 180)

     var eta,omega;
     var eta_rad,omega_rad;
     var index=0;
     var x=0;
     var y=1;
     var z=2;

     for(eta = -90 ; eta <= 90 ; eta+=step_eta) //I could do the cycle in rad but I can have problems in condition of termination of the cycle (float value)
         {
            eta_rad = degToRad(eta);

            for(omega = -180 ; omega <= 180 ; omega +=step_omega , index+=3)
            {
               omega_rad = degToRad(omega)
               vertices[index+x]= a1 * Math.cos(eta_rad) * Math.cos(omega_rad);
               vertices[index+y]= a2 * Math.cos(eta_rad) * Math.sin(omega_rad);
               vertices[index+z]= a3 * Math.sin(eta_rad);

                   normals[index+x]= vertices[index+x] /(Math.pow(a1,2.0)) ;
               normals[index+y]= vertices[index+y] /(Math.pow(a2,2.0));
               normals[index+z]= vertices[index+z] /(Math.pow(a3,2.0));
            }
         }

             //Now triangolation. I believe that it needs a counterclockwise specification of vertices (else I have back-face)
     //each triangle is specified with 3 indices (gl.TRIANGLES).

             var num_samples_omega = 360/step_omega ;   //Actually there is another sample that coincides with the first sample.
     var num_samples_eta   = 180/step_eta ;

             var i,j;
     for(i=0,index=0 ; i<num_samples_eta ; i++) //il primo e l'ultimo campionamento coincidono con un solo punto ma va bene uguale
     {
         for(j=0 ; j<num_samples_omega ; j++,index+=6)
         {
            //Si definiscono ad ogni ciclo gli indici dei 2 triangoli del poligono formato da 4 campioni (due per ciascuna sezione di ellisse adiacente)
             indices[index]=i*num_samples_omega +i+j; //I could use index++ but like this is more clear.
                             indices[index+1]=(i+1)*num_samples_omega +2+i+j;
             indices[index+2]=(i+1)*num_samples_omega +1+i+j;

             indices[index+3]=i*num_samples_omega +i+j;
             indices[index+4]=i*num_samples_omega +1 +i+j;
             indices[index+5]=(i+1)*num_samples_omega +2 +i+ j;
         }
     }

  //alert("num vertici"+vertices.length+"\n"+"indice massimo:"+ indices[index-1])

 }

If you sample the two parameters $\eta$ and $\omega$ with steps $d\eta$ and $d\omega$, then you'll get a grid of points $v_{ij} = f(i\;d\eta,j\;d\omega)$.  Any four adjacent points will define a quadrilateral.  To get triangles, you just have to split each quad in two by a diagonal.

So in the example, you'd split the quadrilateral $\{v_{00},v_{01},v_{11},v_{10}\}$ into two triangles $\{v_{00},v_{01},v_{10}\}$ and $\{v_{01},v_{11},v_{10}\}$.  You'll want to make sure the points come in the right order so the triangle orientations match (or you'll end up backface-culling half of your triangles.)

Your second equation, for $n_{xyz}$, is finding the normal of each vertex, which you can use for shading.


Refraction in a Ray Tracer: What do with an intersection within the medium? I am currently following a tutorial at Scratchapixel.com on refraction

Here is the refract function:

Vec3f refract(const Vec3f &I, const Vec3f &N, const float &ior)
{
    float cosi = clamp(-1, 1, dotProduct(I, N));
    float etai = 1, etat = ior;
    Vec3f n = N;
    if (cosi < 0) { cosi = -cosi; } else { std::swap(etai, etat); n= -N; }
    float eta = etai / etat;
    float k = 1 - eta * eta * (1 - cosi * cosi);
    return k < 0 ? 0 : eta * I + (eta * cosi - sqrtf(k)) * n;
}

1) I understand that if there is no total internal reflection, first we must find the refraction direction (ior is the object's index of refraction):

Vec3f refractionDirection = refract(dir, hitNormal, isect.hitObject->ior).normalize();

2) We then compute the refraction ray origin (which is the point where the ray will leave the object?)

Vec3f refractionRayOrig = outside ? hitPoint - bias : hitPoint + bias;

3) Finally to get the refracted ray color we cast a ray using the refractionDirection & refractionRayOrig

refractionColor = castRay(refractionRayOrig, refractionDirection, objects, lights, options, depth + 1);

So far I understand how a glass block (or any other transparent object which undergoes refraction) bends the light towards the normal when entering the object and then away from the normal when exiting but I am confused as to how refraction works or is computed when there is another object in the medium as shown in the image below.

So when the ray hits the object from the inside do we use a BRDF to calculate the surface properties and then that ray leaves the glass block with the refraction origin and direction?

this is an interesting question (and I am actually an author on Scratchapixel so I can maybe help on that one)).

Things go as follows:

you cast the primary ray into the scene
the ray hits the glass which is a refractive-reflective/transparent material
you compute and cast two rays from the point of intersection: a reflective ray and a refractive ray
if the refractive ray hits an object contained in the transparent object (the glass) then you need to shade that object too. So if it's a plastic material, compute the result of that shader (including the lighting, so loop over all lights and add their contribution, etc.)
the color at the intersection point can finally be computed as a mix between the reflection and the refraction color (where as the mix is driven by the fresnel effect).

Now that works but this is essentially wrong.

Why? Because in fact light rays too are bent by say glass and water. So in fact when you cast a shadow ray from a point on the pen in the direction of the light, this ray direction is in fact wrong, since you should account for refraction and reflection (the ray should be bent when it goes from one medium to another in this case from water to air) but this is impossible to do with this algorithm/approach.

Though in general this gives acceptable results so no one cares much but if your goal is really a physically plausible image then this approach is wrong.

A solution to this maybe is to use photon mapping, where photons of light are cast from the light into the scene in a pre-pass. Thus when a light-photon hits the glass/water surface it will bent due to the law of refraction and eventually it will be "deposited" onto the pen. This is a bit hard to explain here but hopefully there will be a lesson on Scratchapixel about this in the future. Hope it helps.

EDIT

> Thanks for clarifying, I had one other issue though. In the source code provided on the website for that lesson there doesn't seem to be code present to spawn a ray when the light ray leaves the glass cylinder. Is this unnescessary or is it actually already present? If it is could you point out where in the code?

Hum it seems like my explanation wasn't clear if you ask this question. When the refractive ray hits the surface of the pen you need to compute the color of the pen at that intersection point (where the refractive ray intersected the pen's geometry). Now the pen is made out of wood for example which for the sake of simplicity is a diffuse material only (no specular). So in order to compute the color of the pen at the intersection point, you simply execute the traditional illumination loop:

Vec3f computePenColor(const Vec3f& objectColor, const Vec3f& intersectionPoint, normalAtIntersection) {
    Vec3f shadedPointColor = 0;
    // loop over all lights in the scene to cadd their contribution
    for (i = 0; i < scene.lights.size(); ++i) {
        Light currentLight = scene.lights[i];
        // let's say this light source is a point light source
        Vec3f lightDirection = currentLight.position - intersectionPoint;
        // this cast a shadow ray from the shaded point to the light position
        // if this ray intersects an object along the way then this point is in the shadow of this light and the function returns true
        bool shadow = castShadowRay(lightDirection.length(), lightDirection.normalize());
        // this point is not in the shadow of the light so add its contribution
        if (!shadow) {
            // the object is diffuse so apply cosine law only
            shadedPointColor += max(0, normalAtIntersection.dot(lightDirection)) * currentLight->intensity * currentLight->color;
        }
    }
    return shadedPointColor * objectColor;
}

So maybe as you can see more clearly from this example, is that when I say "cast a ray from the point of intersection to the light" I mean cast a shadow ray. This shadow ray helps to determine if the shaded point is in the shadow of that light. If not, then you can add the contribution of that light to the shaded point.

Now you can also see that the direction of this shadow ray is computed as the position of the point on the pen (the intersection point or shaded point) to the light position (we assume a point light source here for simplicity). So this is a straight line between the point on the pen and the light.

But this where things are wrong. This shoudn't be a straight line, since when the "shadow ray" (which is also a light ray but we go from the point on the pen to the light rather than from the light to the pen, but in essence they are the same thing) exits the water volume it should be bent due to refraction (and part of it should be reflected to). You can compute this new direction if you wish, but then if the shadow ray is bent when it leaves the volume of water then it won't travel in the direction of the light anymore (as showed in the pic below), so you can't use it as a shadow or light ray any longer.

So this algorithm that consists of looping over the lights and casting shadow ray can't be used if you really want to account for light ray bending due to reflection and refraction. Though as I said few people care about that because in reality, ignoring this fact, doesn't make much of a visual difference in 99% of cases. And as I pointed out, the solution to get it right if really needed is something like photon mapping where photon of lights are cast from the light into the scene into a pre-pass, which allows to follow the paths of light rays into the scene as they are being bent by refractive/reflective surfaces. But photon mapping is a rather complex algorithm to implement.


A recent approach for subsurface scattering I read "A Practical Model for Subsurface Light Transport" and "A Rapid Hierarchical Rendering Technique for Translucent Materials". If I understand correctly, the former does not consider the case when an object is lit behind and we can see the light is going through it.

The latter can achieve this but there is a preprocessing stage which is not suitable for my CUDA path tracer and I want it to be as natural as possible. So, do you have any paper in mind which is not obsolete or would you suggest me to solve the full radiative transfer equation?

As mentioned in the comments, I would highly suggest starting with Full Volumetric Scattering. This is two fold:

Since you are doing path tracing, adding volumetrics isn't super difficult.
Fully understanding how full volumetric scattering works will be a great basis for understanding the estimations. In addition, it can provide great "references" to see if your estimations are performing well / correctly.

With that in mind, below is a basic introduction to how to implement Full Volumetric Scattering in a backwards path tracer.

To start, let's look at the code for a backwards path tracer with only reflection, no transmission / refraction:

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);
}

In English:

Shoot a ray through the scene.
Check if we hit anything. If not we return the skybox color and break.
If this is the first ray, or we just bounced off a specular surface, check if we hit a light. If so, we add the light emission to our color accumulation.
Sample the direct lighting.
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.
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.

For more details on Direct Light Sampling, see this answer.

BSDFs

To add transmission, we first need to upgrade some of our BRDF materials to full BSDFs, with a BTDF function.

For example, a BSDF for an ideal specular dielectric (ie. glass) is:

class IdealSpecularDielectric : public BSDF {
public:
    IdealSpecularDielectric(float3 albedo, float ior)
        : BSDF(BSDFLobe::Specular, albedo),
          m_ior(ior) {
    }

private:
    float m_ior;

public:
    float3 Eval(SurfaceInteraction &interaction) const override {
        return m_albedo;
    }

    void Sample(SurfaceInteraction &interaction, UniformSampler *sampler) const override {
        float VdotN = dot(interaction.OutputDirection, interaction.Normal);
        float IORo = m_ior;
        if (VdotN < 0.0f) {
            IORo = 1.0f;
            interaction.Normal = -interaction.Normal;
            VdotN = -VdotN;
        }

        float eta = interaction.IORi / IORo;
        float sinSquaredThetaT = SinSquaredThetaT(VdotN, eta);
        float fresnel = Fresnel(interaction.IORi, IORo, VdotN, sinSquaredThetaT);

        float rand = sampler->NextFloat();
        if (rand <= fresnel) {
            // Reflect
            interaction.InputDirection = reflect(interaction.OutputDirection, interaction.Normal);
            interaction.SampledLobe = BSDFLobe::SpecularReflection;
            interaction.IORo = interaction.IORi;
        } else {
            // Refract
            interaction.InputDirection = refract(interaction.OutputDirection, interaction.Normal, VdotN, eta, sinSquaredThetaT);
            interaction.SampledLobe = BSDFLobe::SpecularTransmission;
            interaction.IORo = IORo;
        }

        if (AnyNan(interaction.InputDirection)) {
            printf("nan");
        }
    }

    float Pdf(SurfaceInteraction &interaction) const override {
        return 1.0f;
    }
};

The interesting part of the code is Sample(). This is where a ray will choose whether to reflect or refract. How we do this is really up to us, as shown in this answer. However, an obvious choice would be to sample based on the Fresnel equations.

In the case of an ideal Specular Dielectric, there are only two potential output directions: perfect refraction or perfect reflection. So we evaluate the Fresnel and randomly choose to refract / reflect with proportion equal to the Fresnel.

Media

Next, we need to talk about media, and how they affect light bouncing around. As explained by Nathan Reed in this answer:

  The way I like to think about volume scattering is that a photon traveling through a medium has a certain probability per unit length of interacting (getting scattered or absorbed). As long as it doesn't interact, it just goes in a straight line unimpeded and without losing energy. The greater the distance, the greater the probability that it interacts somewhere in that distance. The interaction probability per unit length is the coefficient $\sigma$ that you see in the equations. We usually have separate coefficients for scattering and absorption probabilities, so $\sigma = \sigma_s + \sigma_a$

  This probability per unit length is exactly the origin of the Beer-Lambert law. Slice a ray segment into infinitesimal intervals, treat each interval as an independent possible place to interact, then integrate along the ray; you get an exponential distribution (with rate parameter σσ) for the probability of interaction as a function of distance.

  You can technically choose the distance between events however you want, as long as you correctly weight the path for the probability that a photon can make it between two adjacent events without interacting with the medium. In other words, each path segment within the medium contributes a weight factor of $e^{-\sigma x}$, where $x$ is the length of the segment.

  Given this, a usually good choice for the distance is to importance-sample it from the exponential distribution. In other words, you set $x = -(\ln \xi)/\sigma$

In summary:

A ray travelling through a medium has some probability to:

Be absorbed
Be scattered

Since this is Monte Carlo Integration, we are free to choose the interaction distance however we want. However, a good choice is to importance sample from an exponential distribution similar to Beer-Lambert.

I chose to implement this with a scattering coefficient, and letting Beer-Lambert take care of the absorption coefficient.

class Medium {
public:
    Medium(float3 absorptionColor, float absorptionAtDistance)
        : m_absorptionCoefficient(-log(absorptionColor) / absorptionAtDistance) {
        // This method for calculating the absorption coefficient is borrowed from Burley's 2015 Siggraph Course Notes "Extending the Disney BRDF to a BSDF with Integrated Subsurface Scattering"
        // It's much more intutive to specify a color and a distance, then back-calculate the coefficient
    }
    virtual ~Medium() {
    }

protected:
    const float3a m_absorptionCoefficient;

public:
    virtual float SampleDistance(UniformSampler *sampler, float tFar, float *weight, float *pdf) const = 0;

    virtual float3a SampleScatterDirection(UniformSampler *sampler, float3a &wo, float *pdf) const = 0;
    virtual float ScatterDirectionPdf(float3a &wi, float3a &wo) const = 0;

    virtual float3 Transmission(float distance) const = 0;
};

A Non-scattering medium is one where a photon travels straight through, but it is attenuated along the way according to Beer-Lambert:

class NonScatteringMedium : public Medium {
public:
    NonScatteringMedium(float3 color, float atDistance)
        : Medium(color, atDistance) {
    }

public:
    float SampleDistance(UniformSampler *sampler, float tFar, float *weight, float *pdf) const override {
        *pdf = 1.0f;
        return tFar;
    }

    float3a SampleScatterDirection(UniformSampler *sampler, float3a &wo, float *pdf) const override {
        return wo;
    }
    float ScatterDirectionPdf(float3a &wi, float3a &wo) const override {
        return 1.0f;
    }

    float3 Transmission(float distance) const override {
        return exp(-m_absorptionCoefficient * distance);
    }
};

An Isotropic Scattering Medium allows the ray to have scattering (reflection) events while the ray is traveling through it. The ray can be reflected in any direction in the sphere around the interaction point (hence Isotropic).

The distance between reflections is random, based on an exponential "scatter" probability.

class IsotropicScatteringMedium : public Medium {
public:
    IsotropicScatteringMedium(float3 absorptionColor, float absorptionAtDistance, float scatteringDistance)
        : Medium(absorptionColor, absorptionAtDistance),
          m_scatteringCoefficient(1 / scatteringDistance) {
    }

private:
    float m_scatteringCoefficient;

public:
    float SampleDistance(UniformSampler *sampler, float tFar, float *weight, float *pdf) const override {
        float distance = -logf(sampler->NextFloat()) / m_scatteringCoefficient;
        // If we sample a distance farther than the next intersecting surface, clamp to the surface distance
        if (distance >= tFar) {
            *pdf = 1.0f;
            return tFar;
        }

        *pdf = std::exp(-m_scatteringCoefficient * distance);
        return distance;
    }
    float3a SampleScatterDirection(UniformSampler *sampler, float3a &wo, float *pdf) const override {
        *pdf = 0.25f * M_1_PI; // 1 / (4 * PI)
        return UniformSampleSphere(sampler);
    }

    float ScatterDirectionPdf(float3a &wi, float3a &wo) const override {
        return 0.25f * M_1_PI; // 1 / (4 * PI)
    }
    float3 Transmission(float distance) const override {
        return exp(-m_absorptionCoefficient * distance);
    }
};

At high scattering distances, the material behaves almost like a NonScatteringMedium, because it has a very small probability of scattering before it travels through the medium.

At low scattering distances, the material behaves like wax or jade.

As the scattering distance gets lower and lower, we have a higher and higher chance to scatter while traveling through the medium. This presents a few issues:

The number of bounces per ray is going to explode.
If Russian Roulette kills the ray, or if you go over "maxBounces", you will "lose" any color contribution from inside the medium.

Thus, you need to set maxBounces to some high number, and rely on Russian Roulette to terminate rays in the general case. In addition, IsotropicScatteringMedium is extremely inefficient for representing mostly opaque materials. To get better performance, we should use a non Isotropic phase function, like Henyey-Greenstein or Schlick. These bias the scattering in a certain direction. So we could set them to have high back-scattering, thereby lessening the issue of "lost" rays.

Putting it all together

So, with this new information, how do we modify the integrator to understand BSDFs and media?

First, we need to start tracking what medium we are currently in, and its Index of Refraction (IOR).

SurfaceInteraction interaction;
interaction.IORi = 1.0f; // Vacuum
Medium *medium = nullptr; // nullptr == vacuum

Then, we break the integrator up into two parts: transmission and material interaction.

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

    float3 color(0.0f);
    float3 throughput(1.0f);
    SurfaceInteraction interaction;
    interaction.IORi = 1.0f; // Air
    Medium *medium = nullptr;
    bool hitSurface = false;

    // Bounce the ray around the scene
    uint bounces = 0;
    const uint maxBounces = 1500;
    for (; 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;
        }

        // We hit an object
        hitSurface = true;

        // Calculate any transmission
        if (medium != nullptr) {
            float weight = 1.0f;
            float pdf = 1.0f;
            float distance = medium->SampleDistance(sampler, ray.TFar, &weight, &pdf);
            float3 transmission = medium->Transmission(distance);
            throughput = throughput * weight * transmission;

            if (distance < ray.TFar) {
                // Create a scatter event
                hitSurface = false;

                ray.Origin = ray.Origin + ray.Direction * distance;

                // Reset the other ray properties
                float directionPdf;
                float3a wo = normalize(ray.Direction);
                ray.Direction = medium->SampleScatterDirection(sampler, wo, &directionPdf);
                ray.TNear = 0.001f;
                ray.TFar = infinity;
                ray.GeomID = INVALID_GEOMETRY_ID;
                ray.PrimID = INVALID_PRIMATIVE_ID;
                ray.InstID = INVALID_INSTANCE_ID;
                ray.Mask = 0xFFFFFFFF;
                ray.Time = 0.0f;
            }
        }

        if (hitSurface) {
            // 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);
            interaction.IORo = 0.0f;

            // Calculate the direct lighting
            color += throughput * SampleOneLight(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;

            // Update the current IOR and medium if we refracted
            if (interaction.SampledLobe == BSDFLobe::SpecularTransmission) {
                interaction.IORi = interaction.IORo;
                medium = material->medium;
            }

            // Shoot a new ray

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

            // Reset the other ray properties
            ray.Direction = interaction.InputDirection;
            if (AnyNan(ray.Direction))
                printf("bad");
            ray.TNear = 0.001f;
            ray.TFar = infinity;
            ray.GeomID = INVALID_GEOMETRY_ID;
            ray.PrimID = INVALID_PRIMATIVE_ID;
            ray.InstID = INVALID_INSTANCE_ID;
            ray.Mask = 0xFFFFFFFF;
            ray.Time = 0.0f;
        }

        // 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;
        }
    }

    if (bounces == maxBounces) {
        printf("Over max bounces");
    }

    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 are currently in a medium (not in a vacuum).

Sample the medium for a distance.
Evaluate the medium's transmission and accumulate the throughput.
If the scatter distance is less than the distance from the ray origin to the next surface, make a scatter event.
Sample the medium for a scatter direction.

Check if we hit a surface (ie, we didn't have a scatter event).

If this is the first ray, or we just bounced off a specular surface, check if we hit a light. If so, we add the light emission to our color accumulation.
Sample the direct lighting.
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.
If we refracted into the material, update the IOR and medium for the next bounce.

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.


OpenGL: Purpose of transform matrix in Vertex shader where transform*vertex I'm trying to understand OpenGL by using Processing OpenGL ES.

If I have an example Cube to render.

I've read that a vertex shader will be run once for each Vertex. The fragment shader will be run (at least) once for each pixel.

In a vertex shader I have:

#define PROCESSING_COLOR_SHADER

uniform mat4 transform;

attribute vec4 vertex;
attribute vec4 color;

varying vec4 col;

void main(){

    gl_Position = vertex;

    col = color;
}

In a fragment shader I have:

gl_FragColor = col;

When ran, it creates a simple painted cube of the colour defined in the Processing sketch. However, if I change the vertex shader to...

gl_Position = vertex

(instead of transform*vertex) it draws a wireframe cube in exactly the same place. What exactly does the transform matrix do?

The transform matrix will be a uniform which is easier to update and upload to the gpu (just a single glUniform call) than applying the transform to all vertices on the cpu and uploading the new vertices. Especially with large models.


Computing a rotation: complex numbers vs rotation matrix A 2D vector can be rotated by an angle $\theta$ using the rotation matrix:
\begin{bmatrix}
\cos(\theta) & -\sin(\theta) \\
\sin(\theta) & \cos(\theta)
\end{bmatrix}

Or, it can be rotated by multiplying the vector by the complex number $c$:

$$c = \cos(\theta) + i\sin(\theta)$$

Is there any meaningful difference between these two methods? I tested both in MATLAB, and they seem to run at the same speed.

On a related note, is there some spatial transformation that complex numbers can do but matrices cannot?

Both methods end up doing the same calculations when you break it down.

Rotating a vector $u$ with a matrix:
$$\begin{bmatrix}\cos\theta & -\sin\theta \\ \sin\theta & \cos\theta\end{bmatrix} \begin{bmatrix}u_x\\u_y\end{bmatrix} =
\begin{bmatrix}u_x \cos\theta - u_y\sin\theta \\ u_x \sin\theta + u_y \cos\theta \end{bmatrix}$$

Rotating a vector $u$ using complex numbers:
$$\begin{aligned}
(\cos\theta + i\sin\theta)(u_x + iu_y) &= u_x\cos\theta + iu_y\cos\theta + iu_x\sin\theta - u_y\sin\theta \\
&=(u_x\cos\theta - u_y\sin\theta) + i(u_x\sin\theta + u_y\cos\theta)
\end{aligned}$$

I wouldn't expect either one to be appreciably faster or slower than the other, since they all end up doing the same set of basic operations, i.e. 4 multiplies and 2 adds.

Conceivably, complex numbers could be faster when you have a large number of rotations stored in an array, because they have only two components instead of four and therefore more of them fit into each cache line.

  On a related note, is there some spatial transformation that complex numbers can do but matrices cannot?

No. Matrices are more general than complex numbers. Any complex number $z$ can be represented by a matrix as:
$$\begin{bmatrix}\text{Re}(z) & -\text{Im}(z) \\ \text{Im}(z) & \text{Re}(z) \end{bmatrix}$$
This corresponds to rotation by the phase of $z$ combined with scaling by the magnitude of $z$. Complex numbers can only represent rotation and uniform scaling. Matrices can represent those, but also nonuniform scaling and shearing. Another way to see it is that complex numbers have only two degrees of freedom, while 2×2 matrices have four degrees of freedom.


Troubleshoot half vector sampling from a distribution (cook-torrence, trowbridge-reitz, etc) I am learning ray tracing and the mathematics behind it.  I have a working monte-carlo raytracer I am experimenting on.  I have gotten past pure diffuse BRDFS, area lights, acceleration structures, etc, and now I'm working on properly sampling the next ray direction from an anisotropic distribution.  I believe I have interpreted the D,G,F factors correctly, and have working tests that verify those function's outputs.

I am having trouble with drawing a "wi" sample from the half vector distribution function.   The purpose of the function is to give you a randomized half vector in shading space based on the "incoming" wo ray, which is the eye ray, and a roughness parameter.  As I graph out the function outputs, it seems like the
half-vector sampling function begins to rotate the spherical coordinates about the wrong axis (I know this makes no sense, but I'm just observing the graphed results).

Knowing that situation, I was hoping some people could post some pointers to good academic articles, open source code examples, or other educational material that could help me investigate and educate myself on the behavior of functions that generate half-vectors from distributions.  I have found plenty of the original source papers from Cook-Torrence, etc that explain their algorithms, but I cannot find any that explain in plain language HOW to generate the half vector for the eye ray to reflect about.

Note: I have also extensively explored PBRT and its source code, all of which has been incredibly enlightening, and I have seen their sample_wi function, i even went so far as to directly copy and translate it to webgl for troubleshooting my functions.  I know i must be doing or interpreting something wrong - I just can't find enough material out there to help me figure out what.

Here is a visualization of only the half vector normals.  They have been put into a range of [0,1] for viewing in the RGB spectrum.  I have also boosted the contrast of the image to make the error more easily visible.  You can clearly see that instead of randomly orienting themselves about the surface normal, they are instead correlating themselves radially about the camera eye-vector.

UPDATE:
here is a webgl "proof" i put together to isolate and illustrate my problem.  What should happen is the color should be fairly blended together, interpolated rather nicely.  Instead I get this obvious hard shift in the distribution throughout the 4 quadrants: http://rdtests.ml3ds-test.com/microfacettesting.html

The functions referenced by main function here can all be viewed directly in the source code of the page I linked to:

        vec3 finalColor = vec3(0.);

        // create a changing "wo" vector across screen space (simulates perspective ray change)
        vec2 screenUv = gl_FragCoord.xy / u_resolution;
        // fire rays into positive z (forward) space from "camera"s
        float z = 1.;

        vec3 rd = normalize(vec3(screenUv, z));

        // remap 0,1 to -1,1
        vec3 wo = (rd * 2.) - 1.;
        // wo is outward facing
        wo *= -wo;

        // convert wo to spherical coordinates
        vec3 Nt, Nb;
        // simulate hitting an infinite plane, facing us like a wall
        vec3 normal = vec3(0., 0.,-1.);
        CoordinateSystem(normal, Nt, Nb);

        // put wo into the shading space coordinate system of the "wall" normal
        vec3 shadingWo = WorldToShading(wo, Nt, Nb, normal);

        // sample half vector distribution from trowbridge-reitz
        // make random input a constant to allow error to be more clearly seen
        vec3 wh = Dist_TR_Sample_wh(shadingWo, vec2(.2,.1), .5, .5);

        // Reflect wo about the half vector normal
        vec3 wi = Reflect(shadingWo, wh);

        // put wi back into cartesian coordinate system
        vec3 wiWorld = ShadingToWorld(wi, Nt, Nb, normal);

        // now remap wi to a viewable color: 0,1
        finalColor = normalize((wiWorld + 1.) / 2.);

        gl_FragColor = vec4(finalColor, 1.);

The general idea for sampling half vector based distributions is that you generate $H$ and then compute $w_i$ by reflecting $w_o$ about $H$. This is so $H$ will be the half vector of your $w_i$ and $w_o$ pair. It is standard reflection:
$$w_i = -w_o + 2(w_o\cdot H)H $$

How you generate $H$ depends on the specific distribution. Generally, it is done in polar coordinates, with the angle from your distribution's center being picked using some specific function. Then the azimuth will be a uniform distribution, unless your distribution is anisotropic, in which case it gets more complicated.

For example, with Cook-Torrance, it goes something like:
$$\tan \theta = m \sqrt{ - \log \xi_1 }$$
$$\phi = \xi_2 2\pi$$

Then, given an orthonormal basis made of vectors $N$, $a$, $b$:

$$H = a\sin\theta\cos\phi + b\sin\theta\sin\phi + N\cos\theta$$

As for problems in your code, replacing this:

vec3 rd = normalize(vec3(screenUv, z));

// remap 0,1 to -1,1
vec3 wo = (rd * 2.) - 1.;
// wo is outward facing
wo *= -wo;

with:

vec3 wo = normalize(vec3(screenUv*2.0-1.0, -z));

Gives me a much nicer plot. Your remapping was incorrect as it was done on a vector which includes z and is no longer in 0,1 for x and y because of the normalization. This explains why the center was not in the center of the frame and the whole thing was leaning. Removing wo *= -wo; (a transformation which makes no sense) and using -z instead fixes what I think you were calling the hard transition in the 4 quadrants. That was caused by squaring the components of the vector.

The remaining glitch in the center appears to be numerical issues. Fixing that will require reworking the code to better handle some cases (likely where you have the 0.0008 constant for a sin and the 0.9999 for a cos).


Deriving blur from real optical formulae Looking at the Wikipedia page for Circle of Confusion, and at such posts as this one, I completely understand how to calculate the values for near, far, focal planes and hyperfocal distance. I am also aware of the Zeiss formula for my desired format.

What I'm trying to find out is how to take real camera parameters (focal distance, aperture, focal length) and get a blur radius or similar usable value to configure a blur with them.

Currently, I'm simply using this function:

float getBlurSize(float depth, float focusPoint, float focusScale) {
    float coc = clamp((1.0 / focusPoint - 1.0 / depth)*focusScale, -1.0, 1.0);
    return abs(coc) * MAX_BLUR_SIZE;
}

All real-world inputs are available to the shader - Aperture f-Number, focal distance, format dimensions (35mm in this case), maximum acceptable CoC from the Zeiss formula, focal length etc., as well as depth and colour. Currently, I am applying the calculations and blur in post.

First we can calculate the physical diameter of CoC in the image plane, given the lens parameters. This equation is from Wikipedia – Circle of confusion:
$$
c = {|S_2 - S_1| \over S_2} {f^2 \over N(S_1 - f)}
$$
where the variables are:

$c$: the physical CoC diameter in the image plane
$S_1$: focal distance (the distance at which a subject would be in perfect focus)
$S_2$: actual distance to the subject (the point whose CoC we are calculating)
$f$: the focal length of the lens
$N$: the f-number at which the lens aperture is set.

Then, divide this CoC diameter by the film or sensor frame size to get the blur diameter in texture coordinates. You can additionally multiply it by your rendering resolution to get the blur diameter in pixels if desired.

As an example calculation, suppose we render a 1920×1080 image, as if through a 85mm lens at f/1.4, on a 35mm camera (frame size 36×24 mm), and focused at 3m while our subject is 5m away. Then we can calculate:
$$
{|5\text{m} - 3\text{m}| \over 5\text{m}} {(85\text{mm})^2 \over 1.4(3\text{m} - 85\text{mm})} {1920\text{px} \over 36\text{mm}}
$$
(Don't forget to convert units, since we are mixing meters and millimeters here!) I get a blur diameter of ~38 px, so the radius would therefore be ~19 px.


Do the alignement and declaration order of uniforms matter? In the section 6.4 Constant Buffers of the book Practical Rendering & Computation with Direct3D 11 (pages 325, 326) it is mentioned:

  By default, the HLSL compiler will attempt to align constants such that they don't span multiple float4 registers. [...] The packing for an HLSL constant buffer can also be manually specified through the packoffset keyword.

I assume a similar rule will apply to the OpenGL equivalent, Uniform Buffer Objects, since they map to the same hardware feature.

What about vanilla uniforms though? What are the rules that apply when declaring uniforms?

uniform vec2 xy; // Can we expect the compiler to pack xy
uniform vec2 zw; // into a same four component register?

uniform vec2 rg;
uniform float foo; // Will this prevent from packing rg and ba?
uniform vec2 ba;   // If so, will foo eat up a full four components register?

If the compiler can do such optimizations, how good are they? Can we explicitly tell the compiler to pack or not, and when should we?

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

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

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

#version 450

uniform vec2 xy;
uniform vec2 zw;

out vec4 v;

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

shader
  asic(VI)
  type(VS)

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

#version 450

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

out vec4 v;

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

shader
  asic(VI)
  type(VS)

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


What's more optimal - sampling one channel at a time or all channels at once? I'm going through and optimizing some of my glsl code for a project, and have come across a scenario several times that I am not sure how best to optimize.

In summary - Is it faster to sample a texture once with .rgba vs 4 times with simply .r at the end?

Here's what I have now:

mat4 getLightProjectionMatrix(int lightIndex){
    mat4 m;

    m[0][0] = texelFetch( projMatricies , ivec2( 0 , lightIndex ) , 0 ).r;
    m[0][1] = texelFetch( projMatricies , ivec2( 1 , lightIndex ) , 0 ).r;
    m[0][2] = texelFetch( projMatricies , ivec2( 2 , lightIndex ) , 0 ).r;
    m[0][3] = texelFetch( projMatricies , ivec2( 3 , lightIndex ) , 0 ).r;

    m[1][0] = texelFetch( projMatricies , ivec2( 4 , lightIndex ) , 0 ).r;
    m[1][1] = texelFetch( projMatricies , ivec2( 5 , lightIndex ) , 0 ).r;
    m[1][2] = texelFetch( projMatricies , ivec2( 6 , lightIndex ) , 0 ).r;
    m[1][3] = texelFetch( projMatricies , ivec2( 7 , lightIndex ) , 0 ).r;

    m[2][0] = texelFetch( projMatricies , ivec2( 8 , lightIndex ) , 0 ).r;
    m[2][1] = texelFetch( projMatricies , ivec2( 9 , lightIndex ) , 0 ).r;
    m[2][2] = texelFetch( projMatricies , ivec2( 10 , lightIndex ) , 0 ).r;
    m[2][3] = texelFetch( projMatricies , ivec2( 11 , lightIndex ) , 0 ).r;

    m[3][0] = texelFetch( projMatricies , ivec2( 12 , lightIndex ) , 0 ).r;
    m[3][1] = texelFetch( projMatricies , ivec2( 13 , lightIndex ) , 0 ).r;
    m[3][2] = texelFetch( projMatricies , ivec2( 14 , lightIndex ) , 0 ).r;
    m[3][3] = texelFetch( projMatricies , ivec2( 15 , lightIndex ) , 0 ).r;

    return m;
}

Here's what I'm considering doing:

mat4 getLightProjectionMatrix(int lightIndex){
    mat4 m;

    vec4 A = texelFetch( projMatricies , ivec2( 0 , lightIndex ) , 0 ).rgba;
    m[0][0] = A[0];
    m[0][1] = A[1];
    m[0][2] = A[2];
    m[0][3] = A[3];

    vec4 B = texelFetch( projMatricies , ivec2( 0 , lightIndex ) , 0 ).rgba;
    m[1][0] = B[0];
    m[1][1] = B[1];
    m[1][2] = B[2];
    m[1][3] = B[3];

    vec4 C = texelFetch( projMatricies , ivec2( 0 , lightIndex ) , 0 ).rgba;
    m[2][0] = C[0];
    m[2][1] = C[1];
    m[2][2] = C[2];
    m[2][3] = C[3];

    vec4 D = texelFetch( projMatricies , ivec2( 0 , lightIndex ) , 0 ).rgba;
    m[3][0] = D[0];
    m[3][1] = D[1];
    m[3][2] = D[2];
    m[3][3] = D[3];

    return m;
}

Does the compiler know how to optimize both of these into the same thing in the end? or would it be more efficient to pack more data into a texture fetch no matter what?

Thanks!

It should be significantly faster to sample one RGBA value from one location in a texture than four R values from different locations in the texture.


BVH uses aabb for a sphere in Ray Tracing: The Next Week While studying Peter Shirley's Ray Tracing:The Next Week, I came across an issue that I cannot quite figure out:
Stepping through the random_scene() case of the code, found out that once a leaf of the BVH tree is reached, the aabb of the object is checked for a hit; even when the object is a sphere.
The check for a hit against the bounding_box of a sphere (aabb) is necessary but not sufficient; there are parts of the bounding_box which is not a part of a sphere.
Shouldn't it invoke the specific hit() function of the object when the leaf is reached?
Yet the code works somehow.

It does indeed do exactly what you described. Have a close look at how the bvh is constructed, it is a passed a list of concrete shapes that implement hittable, e.g. xy_rect, sphere, etc. The base cases for the recursive bvh construction create leaf nodes that directly reference the concrete shapes that make up the scene.
bvh_node::bvh_node(
    const std::vector<shared_ptr<hittable>>& src_objects,
    size_t start, size_t end, double time0, double time1
) {
    auto objects = src_objects; // Create a modifiable array of the source scene objects
    ...
    size_t object_span = end - start;

    // **create leaf nodes, concrete shapes referenced here**
    if (object_span == 1) {
        left = right = objects[start];
    } else if (object_span == 2) {
        if (comparator(objects[start], objects[start+1])) {
            left = objects[start];
            right = objects[start+1];
        } else {
            left = objects[start+1];
            right = objects[start];
        }
    } else {
       // create internal bvh_node
       ...
}

When bvh_node::hit is called and the ray/aabb test passes, the left and right child nodes have their hit functions called recursively. Eventually when this recursion bottoms out the concrete shapes referenced by the leaf nodes will have their hit functions called as you would expect.
bool bvh_node::hit(const ray& r, double t_min, double t_max, hit_record& rec) const {
    if (!box.hit(r, t_min, t_max))
        return false;

    bool hit_left = left->hit(r, t_min, t_max, rec);
    bool hit_right = right->hit(r, t_min, hit_left ? rec.t : t_max, rec);

    return hit_left || hit_right;
}


Pathtracing line and point primitives? So I'm in the midst of implementing a path-tracer. I have chosen to use Assimp library to load models and scene information.

According to it, it will provide mesh data and meshes may contain point, line, triangles and polygons primitives. It gives the option to triangulate polygons but I don't know what to do with lines and points. Should I ignore them?

This leads to a bigger question that came in my mind. Can we ray-trace points and lines? I wanted to use this kind of thing once to show my Bounding boxes in wire frame mode but didn't know how I would trace only the boundaries.

I'm guessing this has to do with providing a specific line/point width? Else they are too thin for intersection.

If we can render lines and points, then should I ignore them or not in my path tracer? Will I be missing out on details etc? Don't know what are they used to represent in 3D models as I have only seen triangles/polygons.

I'm guessing this has to do with providing a specific line/point width? Else they are too thin for intersection.

Exactly. True points and lines have no spatial extent: a point is infinitely small, a line is infinitely thin. You can't intersect those with a ray. But a point with a size is just a sphere, and a line with a size is a cylinder, and those are classic ray-casting primitives, so of course you can intersect them.

The only difficulty is if your path-tracer currently only has one geometry primitive (the polygon or polygon mesh): adding new primitives might be hard for you, if you've structured it for just one. But even in that case, instead of intersecting the points and lines directly, you can just turn them both into polygon meshes inside the renderer, before casting any rays. If you know how to polygonise a sphere and a cylinder you already know what to do.


Upscaling 720p to 1080p vs upscaling 768p to 1080p Suppose I can run some graphics in 1366x768 equally well as 1280x720, and I have a 1920x1080 monitor.

There are more pixels in 1366x768 than in 1280x720.  However, in 1280x720, a 2x2 grid of pixels fits exactly within a 3x3 grid of physical pixels in my monitor.

So which should look better, and is there any theoretical basis for which one looks better?  Does the 3:2 ratio mean anything at all?

With a good-quality upscaling filter (bicubic, for example), there's no particular importance to the 3:2 ratio. Starting from a higher-resolution source will produce a better-looking final image.

It's true that if you're programming the upscaler from scratch, making one for a fixed ratio like 3:2 is much simpler than making one that handles arbitrary ratios. But image-processing software, monitors, etc probably have arbitrary-ratio upscaling built in anyway.


Is there a way to do sphere caustics analytically? I have a Shadertoy which renders the image below using ray tracing.

I'd like to add caustics, but am aware that in general they are a hard solution, especially without bidirectional ray tracing or photon mapping.

Does anyone know if there is a way to analytically calculate caustics, or at least some of the caustics (or plausible caustic effects) in a situation like this?

There is refraction going on here, not just transparency, in case that makes a difference.

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 η.


What can be reduce in peripheral circle and increase in foveated cirle? Foveated rendering, is definitely one of the newest fields of real-time graphics rendering. The main goal of foveated rendering is reducing rendering cost to increase frame-rate without perception difference to the user from full-rendering. Eye-tracker is a great help here, as technology has exponentially evolved over years.
When we use an eye tracker, can gradually decrease the rendering qualities in peripheral zones. The decreased graphics parameters could be, for example, the color quality. Spatial frequency decrease (e.g. text) in the peripheral does not help much in special applications as decrease such quality will increase the reading time.  Could you give me some ideas about the graphics parameters which can be decreased from the angular distance of the gaze point?

Actually it's not that new of an idea -- it goes back to at least 1996.  I strongly suggest reading Microsoft Research's Foveated 3D Graphics.  (click 'View publication' and read the full pdf)
In a nutshell, human visual acuity is different near the center of vision compared to the peripheral.  The middle of the eye is packed more tightly with cone receptors making it more sensitive to detail and color.  In contrast, humans are quite sensitive to motion in their peripheral vision.  Evolutionarily speaking, it's a good idea to maybe take a look at that tiger sneaking up on you from the side.  Your motion perception will trigger a shift in your gaze so you can make sense of exactly what's trying to eat you.
The most obvious thing that's different is your ability to perceive detail (angular size).  Numerous things can be adjusted with the rendering to economize on detail.

Resolution can be reduced.
Texture sizes can be reduced.
Geometric (polygonal) detail can be reduced

Color perception is reduced in the periphery as well, therefore:

Texture maps with lower bit-depth can be used.
Precision of color calculations can be reduced.

CAUTION: There are nuances to making effective use of all of these techniques.  You have to be especially careful not to change the silhouette of the objects because the outer visual region is actually more sensitive to the motion differences that may be perceived by an edge jumping around.  It may also be tempting to turn off antialiasing, but artifacts introduced by crawling pixels are tantamount to motion and draw the eye to these artifacts.  Anti-aliasing sometimes needs to be increased negating the effects of foveated rendering.  A blended blurring technique is used that allows seamless transition between regions rendered at high-resolution and low-resolution. It's a complex subject -- read the paper.


Using instanced arrays for dynamic objects With the help of instanced arrays we can specify an array of transforms and render a number of objects with a single render call. I'm wondering if using instanced arrays make sense in situations where the transforms of the objects need to be updated on each frame? In that case we would need to reupload the buffer with transforms on every frame. Is this practical or should we use individual rendering or gl_Instance_ID (if the number of instances is not too high) instead?

If you need to update data, then you need to update data regardless of the mechanism you use to render the object. If you have 3000 objects to be rendered, and each frame they will be rendered with a different matrix transform, then you will have to send 192,000 bytes of memory (64-bytes-per-matrix) every frame. Whether explicitly through a write to buffer object storage or implicitly through setting uniform values, 192,000 bytes of data are going to have to find their way from CPU-space to GPU-space.
So long as you use reasonable buffer object streaming techniques, performance due to sending data shouldn't be a problem relative to other ways to send the data. Especially for fairly small data like per-model transforms. Of course, good buffer streaming does require having two buffers, but that's a memory storage tradeoff, not a runtime performance cost.


How to get assembly code from GLSL Shader I am new to optimizing shader code. Now that I have written an "optimzation", I want to make sure that it really speeds things up and is not just an unreadable version of what I intend to do. For this, I wanted to look into the assembly version of my shader, but I couldn't really find out how to do that.

I use CodeXL and APITrace, but the former doesn't have any Assembly view as far as I could see and the latter didn't work for my program for some reason. Most other debugging tools don't work (properly) because I program using OpenSceneGraph, meaning there is some Legacy OpenGL Code.

Therefore my question: How would I get to seeing assembly code of my shaders?

In unextended OpenGL that's indeed a little difficult, if possible at all. There is really no actual notion of "GLSL assembly code", or not one that's remotely vendor-independent, let alone a way to retrieve that.

There are a few extensions that might help you get that information, though, albeit in a not entirely intuitive or even platform-dependent way (but well, you're talking about micro-optimizations, so we're firmly in platform-dependent territory anyway).

First of all, there's GL_ARB_get_program_binary (core since 4.1), which let's you retrieve a platform-dependent binary version of the program (after linking). At least on NVidia hardware (maybe others, too, but I wouldn't know) this binary blob also contains ASCII assembly listings of all the contained shaders in NVidia-extended ARB syntax (continuously refined in the GL_NV_gpu_programX extensions). Those listings are quite easy to find in the binary when opening it with a normal text editor and the syntax is quite straight forward if you know your way around shaders and general assembly programming principles. As said, other vendors might offer similar listings inside their binaries.

Then there's also external programs that might compile your shader into some intermediate form, the most platform-independent of those being SPIR-V, for which Khronos-approved external compilers exist. Being an intermediate language this is somewhat of a "high-level assembly language" with an "assembly-like" structure but with many high-level abstractions that builds a bit of a compromise between the easy interpretability of assembly and the easy programability and platform-independence of a high-level language (think of Java bytecode).

With the GL_ARB_spirv extension (core since 4.6) you can directly use (externally precompiled) SPIR-V programs in OpenGL in place of GLSL shaders. I don't know if there is currently a way to let your driver compile GLSL into SPIR-V and then retrieve that. Otherwise some external SPIR-V module created from a reference compiler would not necessarily be equal to the machine code that your driver's actual GLSL compiler produces. It might also be that such an intermediate representation might already be too high-level for analysing micro-optimizations anyway.

All in all, and especially if you just want to use it for debugging/analyzing purposes on a dedicated machine, trying to look for a whatever-natured assembly listing in your driver's binary blob retrieved with glGetProgramBinary might be the easiest and most "close-to-the metal" way to get some insight into the compiled product. I use it occasionally just out of interest.

But of course no answer to this question would be complete without a warning that you shouldn't make too much out of it either. As on the CPU, generally trusting your compiler is a good default approach and you shouldn't let micro-optimizations like this guide your practice too much, let alone rely on any of that guidance being remotely portable (let alone assume that the assembly listed in the binary is actually what's uses as machine code). But if you want to know what your driver might make of it, it can be quite interesting and helpful, especially when deciding between two equally meaningful ways to write some function that might have unexpected differences in the "assembly".


How should I generate the view position for a light at infinity when creating depth maps? I'm working on implementing Bruneton's cloud rendering paper, which requires one to render a pair of depth maps, similar to how shadow mapping works. In this case, my light source is the sun (or possibly the moon), which I represent in my rendering engine as a normalized vector pointing in the direction of the sun or moon, which is essentially located at infinity from the current camera position.

Since the light source doesn't have a traditional "position," I'm having a little trouble determining how best to create the viewpoint transformations for the source.

Since I would also essentially use the same transformation for sun/moon shadowing, it'd kill two birds with one stone.

So my question is thus:
How can I generate a view transform for generating depth maps of area light sources located at infinity?

Instead of using a perspective projection you would use a orthographic projection.

Then the trick is to position the bounding box to in front of the normal camera.

An additional option is to add a skew/shear operation so it maps the light direction to the vertical and keeps the horizontal ground plane horizontal.

This avoids leaving a large dead zone underground near the light source.


Z-Buffer algorithm vs Painter's algorithm? In hidden surface removal, how does the Z-buffer algorithm work vs. painter's algorithm?

In the painter's algorithm, you first sort all graphics elements
on depth (deepest first) and
then one-by-one fully paint them into the image on top of each other.
That way, deeper elements are obscured by less deep element.
(Intersecting graphics element require special attention.)
In the depth-buffering algorithm,
you store the current depth of each pixel in the image
(in addition to other pixel info),
and fully paint each graphics element pixel-by-pixel,
where an image pixel is updated
only if its depth is more than that of the new pixel.

So, there is a time-memory tradeoff:

The painter's algorithm requires an extra sorting step (costs time),
but avoids storing additional (depth) information per pixel (saves memory).
Depth-buffering avoids sorting (saves time),
but stores extra (depth) information per pixel (costs memory).


Why inverse of an allignment matrix is the same as its transpose? $A_{v}^{-1}$ = $A_{v}^{T}$
$A_{v}$ = Allignment matrix that alligns vector v with z axis
$=$\begin{bmatrix}
    \frac{\lambda}{|v|}       & \frac{-ab}{\lambda|v|} & \frac{-ac}{\lambda|v|} &   0 \\
    0 & \frac{c}{|v|} & \frac{-b}{|v|} &   0 \\
    \frac{a}{|v|}       & \frac{b}{|v|} & \frac{c}{|v|} &   0 \\
0 & 0 & 0 & 1
\end{bmatrix}
here, $\lambda = \sqrt{b^2+c^2} $
$|v| = \sqrt{a^2+b^2+c^2} $
And vector $v = (a, b, c)$

To align vectors, you are applying a rotation.

If you consider a simple rotation matrix, e.g. rotation of angle $\theta$ around Z axis, then you will trivially see that its inverse matrix, a rotation of $-\theta$, is the transpose,  since $sin(-\theta)=-sin(\theta)$ and $cos(-\theta)=cos(\theta)$.

An arbitrary rotation matrix, R, can be constructed by multiplication of other rotations, e.g $R = A \cdot B$.  If the inverse of A and B exist, then $R^{-1} = B^{-1} \cdot A^{-1}$.
Similarly we know  $(A \cdot B)^T = (B^T \cdot A^T)$.

Put these together and you have your answer.


When and who introduced the idea of a scene graph? Simple question: is there a source that precedes the idea of presenting graphical application or 3D scene in a scene graph?
I came across Hierarchical geometric models for visible surface algorithms which introduces the idea of having the geometry in easily searchable and cullable format. So I wonder if this is the father of the idea of having scene graph, or has there been something published earlier with this idea.

Not 3D as such, but it seems in 1963 Ivan Sutherland's seminal "Sketchpad" (and wikipedia) had instancing. eg. See "Chapter VI RECURSIVE FUNCTIONS" of the 2nd reference.
IIRC somewhere there's a video of it in action - it made me annoyed as to how some of today's vector drawing tools don't seem to be able to do some of what sketchpad appeared to support.
Update: @Joojaa has kindly located a video   or this


Rendering same vertex-array in different modes here is my vertex array I am using:

vertices = [
    -1.0,  1.0,  0.0,
     0.0,  1.0,  0.0,
    -0.5,  1.0,  0.0,
     0.5,  1.0,  0.0,
     1.0,  1.0,  0.0,
     1.0,  0.5,  0.0,
     0.0,  0.5,  0.0,
     0.5,  0.5,  0.0,
    -0.5,  0.5,  0.0,
    -1.0,  0.5,  0.0,

    -1.0,  0.0,  0.0,
    -0.5,  0.0,  0.0,
     0.0,  0.0,  0.0,
     0.5,  0.0,  0.0,
     1.0,  0.0,  0.0,

    -1.0, -0.5,  0.0,
     0.0, -0.5,  0.0,
    -0.5, -0.5,  0.0,
     0.5, -0.5,  0.0,
     1.0, -0.5,  0.0,
     1.0, -1.0,  0.0,
     0.0, -1.0,  0.0,
     0.5, -1.0,  0.0,
    -0.5, -1.0,  0.0,
    -1.0, -1.0,  0.0,
];

In a html-select I want to select a rendering mode and use it inmy drawScene-method:

var renderModeSelector = document.getElementById("renderModeDropDown");
var selectedRenderModeValue = renderModeSelector.options[renderModeSelector.selectedIndex].value;

switch (selectedRenderModeValue) {
    case "0":
        gl.drawArrays(gl.POINTS, 0, squareVertexPositionBuffer.numItems);
        break;
    case "1":
        gl.drawArrays(gl.TRIANGLE_FAN, 0, squareVertexPositionBuffer.numItems);
        break;
    case "2":
        gl.drawArrays(gl.LINE_LOOP, 0, squareVertexPositionBuffer.numItems);
        break;
    default:
        console.log("No rendering mode selected.");
}

Case '0' is supposed to draw the points which works fine

Case '1' should draw a filled rectangle, with the shader below. Unfortunately it is messed up, see also in the picture below.

Fragment shader:

precision mediump float;

varying vec4 vColor;

void main(void) {
gl_FragColor = vColor;
}

Vertex-shader:

attribute vec3 aVertexPosition;
attribute vec4 aVertexColor;

uniform mat4 uMVMatrix;
uniform mat4 uPMatrix;

varying vec4 vColor;

void main(void) {
gl_Position = uPMatrix * uMVMatrix * vec4(aVertexPosition, 1.0);
vColor = aVertexColor;
gl_PointSize = 5.0;
}

And finally in case '2' it is supposed to draw only the outlines, which also does not work as expected, see below.

So I wonder if it is possible to meet my requirements while using one vertexPositionBuffer. Should I maybe create another vertexPositionBuffer or can I just change mine "on the fly", before rendering?

Your shaders are fine, the problem is topological, that is the way the vertices are connected to create triangles.

For cases 1 & 2: You cannot really draw a grid of points in this mode (GL_TRIANGLE_FAN) as vertex 0 (the top left one) will be used in every triangle, as well as other issues. You should be using Triangle/Line strips, although these aren't so straightforward as some vertices are shared between each row.

For example, imagine you have vertices arranged like so:

  0 1 2

  3 4 5

  6 7 8

You would draw strips by sending the above vertices in the following order:

  0,3,1 (outputs triangle 0,1,3)

  4,    (outputs triangle 1,3,4)

  2,    (outputs triangle 1,2,4)

etc...

This you can do in code by scanning across 1 row of vertices and including the vertex on the next row for every other vertex.

e.g.

int numColumns = 5;   // As per your example..
int numRows = 5;
int vertexOffset = 0;
int vertexSize = 3;          // 3 floats per vertex
int nextRow = numColumns * vertexSize;

// This loop is for each row-strip, which requires 2 vertex rows..
for(int j = 0; j < numRows - 1; j++) {

    // Inner loop sends vertices for each row-strip
    for(int n = 0; n < numColumns; n++, vertexOffset += vertexSize) {
       glVertex3fv(vertices + vertexOffset);            // Current row
       glVertex3fv(vertices + vertexOffset + nextRow);  // Row below
    }
    // This is special, we send a duplicate vertex to restart a new strip for the next row (This can be omitted if j == numRows - 1)
    glVertex3fv(vertices + vertexOffset - vertexSize + nextRow);
}

If you are concerned about sending duplicate vertices you could use an additional index buffer, as that will allow you to use the exact same vertex buffer you used for the GL_POINTS version.


Dynamic Ray-Triangle Intersection I am working on a small simulation software which has a 3D Renderer that is using OpenGL. I render objects to the screen as triangle meshes. To be able to select the triangles on the screen, I have implemented the Möller–Trumbore ray-triangle intersection algorithm. The algorithm works fine. However, I have another problem. My triangles are not static in the world. When I change the location of the triangle on the screen, I have to change the input vertices of the triangle in the algorithm but I don't know how to do it mathematically. Is there a way to use this algorithm for dynamic objects?

There are 2 ways to go about intersecting the triangle. Let the vertices of the triangle have positions $v_1, v_2, v_3$. Let the ray have origin $o$ and direction $d$. Let the model (4x4) matrix be $M$.
To find the new vertex coordinates one extends the positions with a 1 (to allow for translations) and multiplies by the model matrix. Let $u_i = (v_{i,x}, v_{i,y}, v_{i,z}, 1)$ then $w_i = Mu_i$. The resulting vertex positions are: $v_i' = (w_{i,x}, w_{i,y}, w_{i,z})$.
The other option is to transform the ray with the inverse matrix $M^{-1}$ and intersect with the non-transformed triangle. To achieve this extend $o$ with a 4th coord of 1 (to account for translation) and extend $d$ with a 4th coord of 0 (to ignore translation) then multiply both with $M^{-1}$:
$$o' = M^{-1}(o_x, o_y, o_z, 1)$$
$$d' = M^{-1}(d_x, d_y, d_z, 0)$$
Drop the 4th coordinate of $o'$ and $d'$ then intersect with the triangle formed by $v_1, v_2, v_3$.


Bitpacking into buffers with webgl (shadertoy) I'm working on a webgl pixel shader which is writing to a 16 bit floating point buffer (each color channel r,g,b,a is a 16 bit floating point number).

I'm trying to figure out if there are any decent ways to pack the following info into those 16 bits per color channel:

normal (x,y,z where z is always negative)
uv coordinates (x,y between 0 and 1)
material index (an integer. 0 to 3 would be plenty)
distance (a positive float)

I could live without distance if needed but the others are required.

I know that I can store just x,y of the normal and calculate z as needed since its always negative (the normal is in screenspace).

That still leaves 5 values to store in four 16 bit floats: normal.xy, uv.xy, material index.

Does anyone have any good methods for making this happen? It would be nice to have the distance too, otherwise I am limited to directional lighting instead of getting point lights as well.

This is in shadertoy if that context changes any answers.

I'll be sampling this buffer with nearest neighbor sampling, so you don't have to worry about it being stored in a way that interpolates well.

Thanks!

Assuming you are not going to filter the buffer using HW texture filtering: since you only need 4 values for the material index, pack it in the sign bits of uv x and y.


Problem with Basic Diffuse Lighting So recently I have been learning how to implement Diffuse Lighting in DirectX and I have done so with the following code:

float4 PS(VS_OUTPUT input) : SV_TARGET
{
    input.normal = normalize(input.normal);

    float4 diffuse = ObjTexture.Sample(ObjSamplerState, input.TexCoord);

    float3 finalColor;

    finalColor = diffuse * light.ambient;
    finalColor += saturate(dot(light.dir, input.normal) * light.diffuse * diffuse);

    return float4(finalColor, diffuse.a);
}

which has produced the following results:

Before I progress any further with specular lighting I wanted to clean up my code a little so I re-implemented this with a function and few changes.

float3 DirectDiffuseBRDF(float3 diffuseAlbedo, float nDotL)
{
    return (diffuseAlbedo * nDotL);
}

float4 PS(VS_OUTPUT input) : SV_TARGET
{
    input.normal = normalize(input.normal);

    float4 diffuseAlbedo = ObjTexture.Sample(ObjSamplerState, input.TexCoord);

    float nDotL = saturate(dot(input.normal, light.dir));

    float3 diffuseLighting = DirectDiffuseBRDF(diffuseAlbedo, nDotL);

    return float4(diffuseLighting, diffuseAlbedo.a);
}

Which produces:

As you can see the modified code has it's diffuse lighting a little different compared to the original. What am I missing from the original code to get it back to normal?

Your ambient lighting contribution is missing from the second one (:


Screen Space anti aliasing? I have the following challenge. I am given a black and white image with rasterized lines.

Given this input, I want to extend the edges to cover more than one pixel and antialias the edges.
My first attempt was to pass a morphological dilation followed by a gaussian filter. But that straight up did not work, all i got is for the image to look blurry.
I instead tried combining both effects into a single shader and go the following:

This is decent, but you can still see a lot of the original aliasing of the edges.
Assuming that I only have the input image, are there any image based / screen space based AA techniques I can use to make the edges look smooth?

The easiest way to deal with this would be to provide thickness for the edges in the continuous setting. That is, make your edges out of solid capsules/cylinders, then you would not have this issue. Technically, this is is neither a supersampling nor a filtering technique, but rather a reformulation of the problem in the continuous setting. Another equivalent reformulation is considering only the vertices and then your reconstruction function being thick edges between those (using a thick line algorithm).
This doesn't seem to be your setting however? Since you mentioned that you only have the raster (aliased) image already, in which case you have already introduced the error so it is hard to get rid off. If you can "access" the continuous model it is a lot easier to fix. If you cannot, then I could recommend using a different reconstruction filter than the dilation that you are using. For instance you can do a thickening based on distance to pixels (sdf) with a step function. This should remove the issues you are having, at least for the example you showed.


compression of gif animation I have a parchment animation at singingdays.co.uk that is a rather large file of about 80mb which obviously is taking a long time to download. I have managed to compress it down to about 60mb using convert and decreasing the colourmap for all the images down to 256 colours, but would like to compress it more without losing the quality. If I could somehow only dispose of part of the background, so that I'm only disposing the part of the background where the parchment is unfurling and keep the top part after the rolled up animation and shadow has passed I think that would help, but must admit that it is boggling my head trying to work it out, if you can even do that with convert, let alone figure it out. If you can do it with convert could you work out the command line and pass it on? If not could you recommend something else? I am working on linux.

GIF animation compressors normally do this automatically, using transparency in a given frame to avoid storing what has not changed from the previous frame.

The reason it's not working in your case is that your input is bad. The flat part of the parchment is not really static. It is moving slowly upwards, at a rate of about 1 pixel every 20 frames. Fix your animation so the non moving parts are perfectly identical between frames and you should get a much smaller file.

You would also get a smaller file if you avoided the noise in the unfurling part (most visible in the shadow). Such noise barely compresses at all. This looks like it came out of some 3D software so crank up the quality settings until the input is noise free.

If dithering, which is really added noise, is used when producing the images, make sure it does not change from frame to frame. Or just turn it off entirely. You might get a little banding but the file will be smaller. You'll need to decide which you prefer.


Unity: Pass texture pointer from OpenGL to Unity I'm working on a scene in Unity where I have textures loaded at runtime. As a result of that, the scene is being slowed down (even with co-routines), so I thought maybe I could move the texture loading task to a plugin and have it run as a parallel thread..

Is there a way to load a texture with OpenGL in a Unity native plugin, and return a pointer back to the scene for further use?

I am not sure if doing this is right or wrong or actually needed in practice. However, I was curious once & wanted to learn about Unity Native plugins, so i decided to do exactly what you need.

Basically, I wrote a native plugin, which loads image data using C++ & stb_image library, it then passes byte[] to Unity C# side where it gets applied on a quad. There are many things you need to consider if it's a cross platform app, I was able to run it on both PC & Android during the experiment.

Let me know if this is what you want to achieve & I can pass on the code to you.

Best.


Should ray throughput always decrease when bouncing? In the accepted answer for the question What is Importance Sampling?, in the example code there is this line:

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

If the PDF is between 0.0 and 1.0, the division will potentially increase the throughput. I do not think it makes sense. What am I missing?

Yes, throughput doesn't have to be in range [0, 1]. That's a reason why you have to use exposure settings, because for differente accumulation time, brightness of an image might differ. Nathan Reed described it in details here - https://computergraphics.stackexchange.com/a/10244/10129
Your final results and number of fireflies and bright surfaces, depends on PDF that you'll be using. By using Normal Distribution Function PDF, you'll get good results, because NDF is the most determining factor of BRDF shape. Relate to this blog post for details - https://schuttejoe.github.io/post/ggximportancesamplingpart1/

Back to throughput range. Assuming that you're sampling based on NDF, you'll get fireflies in case of having values in specific range, when light source is shining from very low angle and you're looking directly at the object. Therefore, your denominator might be much greater than nominator and result in small PDF causing fireflies:
$$
p_i(w_m, w_o) = \frac{D(w_m)(w_g \cdot w_i)}{4|w_g \cdot w_m|}
$$

Since 2014, Heitz was working on BSDFs, presenting good results in same year by sampling visible normals [2]. In 2017, he published much simpler, faster and better paper, improving his method [3]. There are also 2015 and 2016 papers in related field, but I haven't read them, but it might be beneficial if you're interested in BSDF. I would recommend checking his first 2014's work (1) where he explains physical assumption on which BRDF is based. It's worth knowing, what NDF and Geometry Shadowing exactly means in BRDF.
Also, there is a lot of good reading in pbrt book, especially:

http://www.pbr-book.org/3ed-2018/Light_Transport_I_Surface_Reflection/Sampling_Reflection_Functions.html#MicrofacetBxDFs
http://www.pbr-book.org/3ed-2018/Light_Transport_I_Surface_Reflection/Path_Tracing.html
http://www.pbr-book.org/3ed-2018/Monte_Carlo_Integration/Importance_Sampling.html
For more advanced method, I think it'd be worth reading about Russian Roulette and more efficient Path Tracing algorithms like Bidirectional Path Tracing, Metropolis Light Tracing and Photon Mapping methods

(1) Eric Heitz - Understanding the Masking-Shadowing Functionin Microfacet-Based BRDFs http://jcgt.org/published/0003/02/03/paper.pdf
(2) Eric Heitz, Eugene d’Eon - Importance Sampling Microfacet-Based BSDFs using the Distribution of Visible Normals https://hal.inria.fr/hal-00996995v1/document
(3) Eric Heitz - A Simpler and Exact Sampling Routine for the GGXDistribution of Visible Normals https://hal.archives-ouvertes.fr/hal-01509746/document


Where should the ray origin be put in perspective camera ray tracing? Where should the ray origin be put in perspective camera ray tracing?

The following gives a picture of the ray tracing in the perspective camera setting:
https://www.cs.cornell.edu/courses/cs4620/2011fa/lectures/06raytracingWeb.pdf

and says that the $e$ is the ray origin, which is the same for all rays.

However, I don't understand how $e$ is supposed to be selected.
Can it be the zero vector? Or perhaps it is supposed to be the camera center?

I don't get what you don't understand. All rays originate from the camera so It's basically the start of the camera's pyramid or more technically the viewing frustum.

Depends on the system you are using. Maybe your camera is fixed at the origin and only the world rotates/moves. In this case, usually the camera is fixed at the origin looking in either the +Z or -Z direction depending on if it's left or right handed coordinates system. So the ray origins will be fixed at (0,0,0).

If your camera moves, however, the ray origins  also move to where you place the camera.


Easing curves on the GPU? I was thinking that I could pass a timestamp in (from requestAnimationFrame) to the vert shader, and just have my easing curves in the GPU instead of using (for example) Tween.js from the outside.

Is this a good idea? Why or why not?

EDIT: After learning more GL and thinking about it, it doesn't seem like I want every shader instance to run a tween that I may only use for a handful of vertices. Can we have only certain shaders run tween calculations if the vertex being animated is the one of the current shader?

This is a pure optimisation problem. It depends on many factors like:

What you need to do with the tween
How many different tweens you need to calculate
How much signaling between systems you need
Whether you need to use that data in your CPU bound process or not.
How much time you have to write the code
What your tool set supports easiest
...

See if you only need one value to interpolate then your CPU is more than likely faster (unless you're running some weird embedded system). Even if you move a bunch of points you can still calculate the tween of that one position and do the needed cluster transformation in the GPU, like bone animation.

The GPU only starts to win against your processor if you need to calculate many truly separate animation curves in parallel. Or you're in an otherwise performance bound system like running JavaScript in a WebGL application. If you do use the GPU then it's most likely best to tween all vertices, at least unless your object is somehow massive. Since initiating a new draw call is costly and maintaining many states is more complex. Its just that those vertices that don't move would just move from starting point to starting point.


Range of integers returned by glGenBuffers? I'm developing a general Mesh class with the help of LearnOpenGL. In the example he creates a class which assumes that you want to generate a VAO, VBO, and an EBO but I don't usually use EBOs so I only want to generate an EBO under certain circumstances. So if I decide not to generate an EBO then I want to set my EBO to a number which cannot be generated by glGenBuffers.

// Draw Mesh
glBindVertexArray(this->VAO);

if(EBO not in range generated by glGenBuffers)
    glDrawArrays(GL_TRIANGLES, 0, numVertices);
} else {
    glDrawElements(GL_TRIANGLES, numIndices, GL_UNSIGNED_INT, 0);
}

glBindVertexArray(0);

My question is what is the range that a buffer object's name can be in?

According to the OpenGL spec, any integer other than zero can be a potential name. Zero is the only value that's guaranteed not to be an object name.
Relevant quote from the spec:

2.6.1.1    Name Spaces, Name Generation, and Object Creation
Each object type has a corresponding name space.  Names of objects are represented by unsigned integers of type uint.  The name zero is reserved by the GL; for some object types, zero names a default object of that type, and in others zero will never correspond to an actual instance of that object type.
Names of most types of objects are created by generating unused names using commands starting with Gen followed by the object type.  For example, the command GenBuffers returns one or more previously unused buffer object names.