Best practice for managing many shaders and merging multiple shaders into one?

I target mobile/tablet devices with my product
  this is less of a concern with desktop GPUs.
  unity and unreal allowing full flexibility on adding shaders & materials at-will

Well, Unity and Unreal can also be used to make games targetting mobile/tablet devices, right?
Do you think you can throw at them all the shaders & materials you want, and expect them to magically output a performant application every time?
Try to do it and you will find that they also struggle at some point.

This can also happen on desktop ; the newest smartphones are totally faster than my oldest laptop.

  What is the best practice for managing multiple (dozens) of shaders in a scene with many materials?
  how many parameters are too much?

The best practice is to test on the real hardware that you are targetting. If you target a wide range of mobile devices, then you have to decide the lowest spec that your product will officially support.
Then, it is a matter of finding a good balance between performance (few materials/shaders) and visual quality (more variety/complexity).
If you are comfortable with the process, you can prepare a number of performance profiles based on the type of device that will run your product (lower-end, higher-end?), with better quality for the higher end.

  Since shaders do have a lot of parameters populated from the model
  material/mesh and sent down to the fragment shader, I fear that
  unifying the shaders and passing too many parameters in the VBO each
  frame/vertex/fragment, would be really bad for performance.

You meant Uniform Buffer Object, not Vertex Buffer Object, right?
I expect that the size of the uniform data that you pass to the GPU is not that big?
However if your shaders are complex, it could take time for the GPU to execute them.

Here again it is a balance, between the amount of parameters that you pass, and the complexity of the shader, that you have to find.
You just can try different solutions on real hardware, note down the performance and aspect of each one, and somehow compare them ; to find the optimal solution for your product.


In Monte Carlo integration, the samples $x_1, x_2, \ldots x_N$ are independent, identically-distributed random variables. This implies they all have the same expectation value. The derived quantities $f(x_i)/p(x_i)$ also all have the same expectation value. So the expectation of a sum of $N$ of these is the same as $N$ times the expectation of any one of them:
$$E\left[ \sum_{i=1}^N \frac{f(x_i)}{p(x_i)} \right] = N \cdot E\left[\frac{f(x_1)}{p(x_1)}\right]$$
(I picked $x_1$ arbitrarily.)

Then, the expectation value of a continuous random variable is found by integrating over its probability measure:
$$E\left[\frac{f(x_1)}{p(x_1)}\right] = \int_{\text{domain of }p(x)} \frac{f(x)}{p(x)} \, p(x) \, dx$$


Is it good practice to use all the available texture units?

Modern hardware doesn't really have the concept of texture binding points as exposed by OpenGL. Rather, the shader unit uses a descriptor (which is just some kind of fat pointer) which can potentially address any texture as long as it's resident in video memory. This is what makes things like bindless textures possible. So the large amount of "texture units" available in current implementations is simply trying to handwave this now-irrelevant part of the API.

That said, the recommended way of implementing what you're trying to do (avoiding re-binding of textures) is to use texture arrays, which allow you to dynamically index into a set of textures in the shader, as long as they all have the same format and size. This presentation contains more details about these and other techniques for reducing driver overhead when rendering with modern OpenGL: Approaching Zero Driver Overhead


Since the cube is a convex object, you could first find the silhouette edges of the cube, by using dot product between cube face normals & view vectors. Then you can sort the edges by the smallest screen y-coordinate for each edge and rasterize the cube a scanline at the time by filling the span between two active edges for each scanline. For each scanline you progress the two edges by adding the edge slopes (x1-x0)/(y1-y0) to the current x-coordinate of the edges.


The signature of glfwGetProcAddress is:

GLFWglproc glfwGetProcAddress(const char *);

It's a function that takes a string and returns a GLFWglproc. which is defined as

typedef void(* GLFWglproc) (void);

Essentially, that means glfwGetProcAddress is a function that takes a string and returns a pointer to function with no parameters and no return value.

GLADloadproc is a type alias defined as follows:

typedef void* (* GLADloadproc)(const char *name);

Essentially, this means that this is the type of a function that takes a string and returns a void pointer (this is a type of pointer that can be safely converted to from every other type of pointer).

in the expression

gladLoadGLLoader((GLADloadproc)glfwGetProcAddress)

Inside, there is the expression

(GLADloadproc)glfwGetProcAddress

Where we have a type in parenthesis followed by an identifier. This is known as a type cast, it means we are taking some value of some type and converting it to some other type.

In this case, we take our GLADloadproc function that takes string and returns pointer to function and we 'disguise' it as a function that takes string and returns void*. This works because we can convert any pointer to void*, even in function pointer return types.

This is all to end up passing it to gladLoadGLLoader, which expects a value of type GLADloadproc.

TL;DR

gladLoadGLLoader takes a GLADloadproc but we want to give it a GLFWglproc (*)(const char *). To be able to do this, we force the type system to take it anyways by using a type cast.


Why do identical polygons not z-fight in OpenGL?

There are explicit rules in the OpenGL specification about invariance guarantees. For the purpose of this discussion, these rules boil down to the following.

Given a particular series of vertex processing shaders (VS, tessellation, and GS), there are some set of input values, uniforms, and the like which lead to the computation of the final gl_Position value for the primitives in that process. If all of these values are binary identical (including the way the vertex format leads to any VS inputs), then you're guaranteed that gl_Position will have binary identical values. And if two primitive of the same type have binary identical positions, the rasterizer is required to generate binary identical fragment positions.

And thus, no z-fighting.

Now, that's kind of restrictive, as it only works with the same shaders. You can use different shaders, so long as:

all shaders involved get binary identical values on expressions leading to the computation of gl_Position,
all shaders involved use the exact same expressions (and I mean exact; a + b isn't the "same expression" as b + a, no matter what mathematics says) to compute gl_Position, and
gl_Position and any inter-shader in/out variables leading to its computation are qualified with the invariant qualifier in all shaders. This tells the compiler not to optimize certain things away in expressions leading to the computation of such values.


The effect is called chromatic aberration; the way it’s usually done in games, and what’s being done there, is by reading the red / green / blue channels of an image at different offset positions. As with your other question (you’re looking to get into the NFT market, I’m guessing?), there are most likely plugins for After Effects that will let you create something along these lines, but the original author probably used a tool like Processing or Shadertoy.


Typically you will need VRAM to store your frame buffer (memory used to store pixels for use with the display and for rendering). This would usually be 4 bytes per pixel so a 1080p frame buffer is 1920x1080x4 bytes (1 byte each for R,G,B components + 1 padding byte). Your card is going to use 2 of these for double buffering so it can use one for display output and the other for constructing the next frame. Any increase in resolution will require an increase in VRAM usage. If HDR is used then you're looking at 16 bytes per pixel so you can use floating point accuracy. Certain types of post-processing special effects may also need additional rendering buffers which are some division of the screen resolution so will consequently use more VRAM if using a higher resolution.

After that the majority of VRAM is going to be allocated for textures which for a typical game is the biggest user of VRAM. Then mesh data needs storing in VRAM which usually consists of vertices storing information like geometry position/normals/colours/transformation matrices/etc. Finally compiled shaders will need to occupy VRAM but they should be low down on memory usage compared to textures and mesh data. There will be other uses for VRAM too depending on how the game was engineered so this is not an exhaustive list of what VRAM is used for.

In short: Display resolution will change the amount of VRAM required but not by a large amount. For a video game most VRAM is from textures which are usually defined by the quality settings of the game and are not affected by display resolution. The 12gb card you mentioned will allow higher texture resolution but with your choice of 1080p or 4k display resolution making little difference in VRAM usage. However for 4k it makes sense to use the highest texture quality.


I can find two possible reasons for the image not converging.

#1. Every sample is the same

For every sample, you generate random rays. You do that when you shoot the ray through a pixel (for anti-aliasing and DoF) and when you sample the hemisphere (for a new indirect bounce). The problem would be if for every sample, it would generate the same direction, then every sample is the exact same colour and you would not be doing any converging.

Most of the times when generating random numbers on the GPU, we generate them in advance and save them to a texture. When we need a number, we just sample the texture. Of course, for every pixel you want a different random value, and you would then use something like the thread id or the pixel coordinate to sample the texture. Then, every frame must also be different, so you would also use the frame counter to offset the thread id. But, notice how there is nothing to make sure that every sample also has a new random number. If for every sample you generate the same random direction, then every sample is the exact same colour and you would not be doing any converging.

This can often easily be fixed by keeping track of how many random numbers you
already have generated and use that to offset every lookup. Basically, have a struct that holds the u and v values for the texture lookup. You initialize it with the thread id or pixel coordinate together with either the frame count or just a random number you generate on the CPU and pass it to the kernel as a global offset to the random values. Make sure that the struct is initialized in the very beginning of the kernel (where you initialize total_color). Then when you need a random number, you pass a reference to that struct to the method that generates a random number, it then looks up the random number in the texture and then it offsets the u and v values in the struct so that the next time you will get a different number.

Of course, I have not seen how you generate the rays (both on shooting a ray through a pixel and sampling a hemisphere), so I cannot say for sure if this is the problem, but it does seem very much like it.

#2. Not properly averaging multiple frames

  I thought taking 64 samples 2 times, and averaging them would be analogous to taking 128 samples once.

This is actually true. You average samples to get an image.
$$IMG1 = \frac{sample1 + sample2 + sample3 + samples4}{4}$$
Which you then do again for the second image.
$$IMG2 = \frac{sample5 + sample6 + sample7 + samples8}{4}$$
Which you then combine.
$$RESULT = \frac{IMG1 + IMG2}{2}=$$
$$\frac{\frac{sample1 + sample2 + sample3 + samples4}{4}+ \frac{sample5 + sample6 + sample7 + samples8}{4}}{2}=$$
$$\frac{\frac{sample1 + sample2 + sample3 + samples4 + sample5 + sample6 + sample7 + samples8}{4}}{2}=$$
$$\frac{sample1 + sample2 + sample3 + samples4 + sample5 + sample6 + sample7 + samples8}{8}$$
As you can see, it would indeed be like sampling twice as much.

However, your piece of code does not properly combine two images.

float4 prev_color = read_imagef(outputImage, sampler, pixel);
color += prev_color;
color /= 2.0f;
write_imagef(outputImage, pixel, color);

When we look at what you are doing.
$$outColor = \frac{newColor + oldColor}{2}$$
Where $outColor$ is the colour saved in the image.
$newColor$ is the colour you just calculated.
$oldColor$ is the color of the previous frame.

We see that it does what it is supposed to do, it averages both frames.

Of course, when the third frame comes around, $outColor$ becomes $oldColor$ as you swap the buffers of outputImage and inputImage. Thus, we can substitute it.
$$outColor = \frac{newColor + \frac{prevColor + oldColor}{2}}{2}$$
Where $prevColor$ was the color calculated the previous frame.

Although, this does not seem that good.
$newColor$ contributes 50%, while $prevColor$ and $oldColor$ only contribute 25%. This is a problem. You want every frame to contribute equally because then you are just adding up the samples. This means that you need to compensate for this.

There is an extremely simple solution to this. Multiply $oldColor$ by the amount of frames excluding the current frame and to then divide by the amount of frames including the current frame. you would then get this.
$$outColor = \frac{newColor + 2*\frac{prevColor + oldColor}{2}}{2 + 1}$$ Which would cancel out the divide and give you.
$$outColor = \frac{newColor + prevColor + oldColor}{3}$$ (The $+1$ is because we multiplied by 2 at the top, basically making it a total of 3 parts, and thus we need to divide by 3.)

If you would be calculating the fifth frame you would do it like this.

float4 prev_color = read_imagef(outputImage, sampler, pixel);
color += prev_color * 4.0f;
color /= 5.0f;
write_imagef(outputImage, pixel, color);

If you would be calculating the tenth frame you would do it like this.

float4 prev_color = read_imagef(outputImage, sampler, pixel);
color += prev_color * 9.0f;
color /= 10.0f;
write_imagef(outputImage, pixel, color);

This does mean that you need to pass through a counter which tells you how many frames you already have calculated. Most progressive or interactive path tracers actually only do one sample per frame (make sure that the ray through the pixel is randomly jittered), since it would take less time to calculate one frame, so that you see changes much earlier and often it is fine if you edit your scene using only one sample since you still get enough information.


Temporal by definition will ghost. I’d limit past frame sampling, apply weight < 0.5, limit movement etc.


Why use image load / stores instead of FBOs

Why would you use Images instead of Framebuffers?

Image Load/Store is not intended to replace framebuffers. It's a feature that allows (relatively) arbitrary reading and writing to memory. So you would use images when you need to arbitrarily read and write to memory. You would use framebuffers when you're rasterizing primitives for rendering purposes.

The fact that you can build scenarios where you would use the latter in place of the former does not mean that this is the purpose of the feature.

  What advantages do they offer?

Again, it's not meant to replace framebuffers. But the comparison is pretty simple.

A particular FS invocation cannot:

Choose what framebuffer texel(s) to write to. An FS invocation will only ever write outputs to the framebuffer texel(s) that map to its current fragment location. gl_FragCoord is an input value; it is non-modifiable. An FS's outputs will either be discard or written to the framebuffer texel(s) corresponding to gl_FragCoord.xy.

Yes, through point rasterization, you can select which fragment is being generated. But it isn't the fragment shader itself that's deciding this; it's whatever process wrote the vertex to the buffer that it's being read from. That is not under the direct control of the FS invocation itself.
Choose to write to multiple distinct framebuffer texel(s) locations. A fragment cannot cover multiple texels.
Choose which output attachments in the framebuffer to write to. The decision of which fragment shader outputs will be written is made by the write mask state, which the FS has zero control over. An FS can choose to be discarded, but that discards all of its outputs.
Choose to read from the current framebuffer texel(s) it is writing to. Blending is done by a fixed-function stage; the FS has zero control over this. The closest you can get is to employ some texture barrier gymnastics, but that's not something the FS can directly control.
Choose to read from framebuffer texel(s) other than the one it's writing to. Even with texture barrier gymnastics, you cannot read from other people's texels.

Now granted, image load/store doesn't exactly make read/modify/write passes easy, due to incoherent memory access rules.

Of course, image load/store cannot handle (easily):

Depth testing. Or at least, not in any way that is nearly as efficient as the rasterizer. The read/modify/write nature of the depth test makes this exceedingly tricky with image load/store. With a regular depth buffer however, this is trivial.
Simple blending. Again, you can do it, but it's really hard, and far less efficient than the normal rendering process.


It is not necessarily affine. An affine matrix in homogeneous coordinates has a form like:

$$\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ 0 & 0 & 1 \end{bmatrix}$$

(assuming you use a column vector convention). Here, the upper-left 2×2 submatrix is the linear part, and $(a_{13}, a_{23})$ is the translation vector of the affine transform.

If the lower row of the matrix has some values other than $[0, 0, 1]$, then it is in general a projective transform, not an affine one.


Why is rotating this prism in iOS Metal not rendering correctly?

The problem was in the shader, and I had the order of multiplication wrong with the MVP matrix.

So this:

out.pos = v * mvp_matrix;

should be:

out.pos = mvp_matrix * v;

For whatever reason, I'm used to dealing with row vectors and had myself convinced the problem was due to the matrices themselves and the clip region..


Extract visible vertices from a 3D geometry model

Idea A:
Draw an invisible mesh that will occlude the points we don't want.

Create a mesh from the point cloud.
Render that mesh to a depth buffer but not to the color buffer.
Render the point cloud using a depth test "closer or equal".

This approach should give the expected result, but the problem with it is the first part, which is not trivial at all.
Idea B:
Approximate the mesh with fat particles, that are big enough to hide points in the background.

Pick a size that is roughly the typical distance between two points.
Use that size as a particle diameter.
Render the point clouds using such large particles, to a depth buffer but not to the color buffer. Maybe offset that point cloud a little, so it's a bit farther from the camera.
Render the point cloud using points this time, and using a depth test "closer or equal".

This method should be a lot simpler to implement than A, because there is no need to triangulate anymore. The result will not be perfect though: some points that should be visible might disappear, and some that shouldn't be visible will remain visible. Tweaking the particle diameter and offset will yield different results.
Idea C:
As a very cheap trick, just render the point cloud with a fog function. The farther points are from the camera, the less visible they are. This method won't hide any point, but still help reading the volume.


How do I replace a mesh's texture with vertex colours?

You can sample the texture in your vertex shader and set the color appropriately. You don't say what system you're using, but for something like an OpenGL vertex shader, it would look something like this (this is a little old school, so it may need updating for modern GL) :

uniform sampler2DRect texture;
varying vec4 color;
void main ()
{
    color  = texture2DRect(texture, gl_MultiTexCoord0.xy);

    vec4    newVertex   = vec4(gl_Vertex.xyz, 1.0);

    gl_Position = gl_ModelViewProjectionMatrix * newVertex;
}

Your fragment shader would then just set the color of the fragment to the interpolated color coming from the vertex shader. Something like this:

varying vec4 color;
void main()
{
    gl_FragColor = color;
}


The constructor of hyperboloid in pbrt-v3

They must define it as a surface of revolution of a line segment.
From the comments, they are probably the a and c coefficients of the canonical hyperboloid equation. With a surface of revolution, a == b so no need to store b separately.
It's clearly related to figuring out ah and ch from the line segment. I can't say more precisely without working out the math involved. I'll leave that to someone with more free time.


What's the difference between irradiance and fluence/radiant exposure?

I can only cite what I have learned in my lecture on Global Illumination techniques which was unfortunately some time ago:

Radiant Power : The amount of energy emitted by a light source in unit time. denoted by $\Phi$ and is measured in Watt which equals to Joules per second. (This does not specify any area!)
Irradiance: The irradiance denotes the incident radiant energy on a surface. It is defined as $E=\frac{d\Phi}{dA}$ - This quantity denotes the incoming energy per area and is measure ind $\frac{W}{m^2}$ Note here: radiant energy is measured in Watt which is Joules per second!
Radiant exitance: This is the reverse direction of Irradiance and is the energy an area emits in unit time and has the same definiton and units as Irradiance only with reverse direction $M = B = \frac{d\Phi}{dA}$. This is also often called radiosity B.
Radiance: Is the power per projected area per solid angle. It is defined as $L=\frac{d²\Phi}{d \omega * dA * cos\theta}=[\frac{W}{sr*m^2}]$ where $\omega$ denotes the solid angle, $\theta$ denotes the angle between emitting surface normal and the direction the energy is traveling to. So this is not attenuated by the distance but by the angles.

In general as I learned the distance between emitting and receiving surface is only important for Radiance because the projected area will be smaller if they are farther apart.

The time you are missing in your definitions seem to come from the conversion of Watt to Joules Per Second.


Using extensions in Android OpenGL ES 2.0

I have encountered devices that advertise certain extensions via glGetString(GL_EXTENSIONS) and do not expose the entry points. I have also encountered the opposite, the extension was not in the extension strings but the entry points were there and it worked. It depends on whether the vendor has disabled any features in the driver from the GPU provider.

You will need to check both the extension strings and check if the entry point is not NULL via eglGetProcAddress.

Apart from that, each extension should conform to the Khronos registry of extensions so you can include gl2ext.h directly from there if you wanted.


There are two general ways the graphics driver can implement Map().

One way is to set up a virtual memory mapping that makes the resource's actual memory (could be system memory or GPU memory) visible as part of the application's address space. If it does happen to be GPU memory, any writes you do would typically be buffered up and sent across the PCIe bus to the GPU, using write combining. (Even on a unified-memory system, mapped resource memory probably has some different caching behavior from "regular" memory, since CPU and GPU caches aren't coherent with each other.)

The other way the driver can implement Map() is to allocate memory and create a local copy of the resource, let your app edit it, then copy the modified memory back to the resource when Unmap() is called. This might need to be done if the actual memory layout of the resource on the GPU is "weird" somehow, like having a nonstandard tiling scheme or some such.

The copying approach is clearly more expensive due to the extra copies required, so I'd expect the driver to use the memory mapping approach whenever possible.


How can I generate procedural noise on a sphere?

I'd consider just going with 3D noise and evaluating it on the surface of the sphere.

For gradient noise which is naturally in the domain of the surface of the sphere, you need a regular pattern of sample points on the surface that have natural connectivity information, with roughly equal area in each cell, so you can interpolate or sum adjacent values. I wonder if something like a Fibonacci grid might work:

I haven't chewed through the math to determine how much work it would be to figure out the indices of and distance to your four neighbors (I don't even know if you end up having four well-defined neighbors in all cases), and I suspect it may be less efficient than simply using 3D noise.

Edit: Someone else has chewed through the math! See this new paper on Spherical Fibonacci Mapping. It seems that it would be straightforward to adapt it to sphere noise.

If you are rendering a sphere, not just evaluating noise on the surface of a sphere, and are fine with tessellating your sphere to the resolution of your noise lattice, you can create a geodesic grid on the surface of the sphere (a subdivided icosahedron, usually):

Each vertex of the sphere can have a randomly generated gradient for gradient noise. To get this information to the pixel shader (unless you want straightforward interpolation like value noise), you may need a technique like this article's wireframe rendering with barycentric coordinates: do unindexed rendering, with each vertex containing the barycentric coordinates of that vertex in the triangle. You can then read from SV_PrimitiveID (or the OpenGL equivalent) in the pixel shader, read the three noise gradients from the vertices based on what triangle you are on, and use whatever noise calculation you like using the interpolated barycentric coordinates.

I think the most difficult part of this method is coming up with a scheme to map your triangle ID to three samples in order to look up the noise values at each vertex.

If you need multiple octaves of noise or noise at a finer resolution than your sphere model, you may be able to do a coarse geodesic grid with vertices and do a few levels of subdivision in the pixel shader. i.e. from the barycentric coordinates, figure out which subdivided triangle you would be in if the mesh was further tessellated, and then figure out what the primitive ID and barycentric coordinates would be for that triangle.


Can I use QueryPerformanceCounter to measure single draw or compute call performance?

That diagram includes a CPU block titled "wait for GPU". I do not see the part of your code that includes an equivalent command.

Furthermore, even if it did include glFinish/glClientWaitSync/etc, there is no guarantee that the CPU time in any way relates to how much time the GPU spent processing it. glFinish might block on a mutex, thus ending the thread's timeslice. Which means that the thread will start up sometime after that mutex is released. How long after? Who knows, but your QueryPerformanceCounter call will include that time.

If you want to measure GPU time, you should use tools that actually measure GPU time. For example, timer queries. You bracket some commands with glBeginQuery(GL_TIME_ELAPSED, queryObject) and glEndQuery(GL_TIME_ELAPSED). Then you can read the result after some time has elapsed. Like maybe next frame.

And unlike glFinish, this will not absolutely murder your performance ;)


understand a pseudo angle implementation

The "pseudo-angle" this function defines runs from 0 to 1 counterclockwise around the circle, with 0 at the negative x-axis, and each 25% of the range mapping to one quadrant (so 0.25 at the negative y-axis, 0.5 at the positive x-axis, etc).
First, the expression dx / (math.Abs(dx) + math.Abs(dy)), can be seen as normalizing the vector in the 1-norm or taxicab norm as opposed to the usual Euclidean norm. In other words, it's projecting the vector onto the unit diamond:

After that, it looks at the sign of y and remaps x accordingly. For the lower half of the diamond (negative y), it maps x linearly from [−1, 1] to [0, 0.5]. For the top half it maps x from [1, –1] to [0.5, 1] (reversing the direction of the mapping this time). So, as you start from the point (−1, 0) and go counterclockwise around the diamond, the pseudo-angle result will increase from 0 to 1.
It is probably used because the system wants to sort vectors by their angle (or in this case, incorporate the angle in a hash apparently), but it doesn't need the actual Euclidean angle which would require an expensive atan2 function to calculate. It just needs something monotonically related to the angle but cheap to calculate, and this pseudo-angle fits the bill.


What is the difference between radiance and irradiance in BRDF

First of all, irradiance at a certain point of a surface is the density of radiant flux (power) per unit of surface area, while radiance at a certain point of a surface in a certain direction is the density of radiant flux per unit of surface area and unit of solid angle. Vaguely speaking, irradiance is the amount of light incoming to a certain point from possibly all directions, while radiance is the amount of light incoming to a point from a single direction. Why on Earth do we talk about irradiance from a direction then?
The thing is that the BRDF definition assumes that there is radiance $L_i\left(\omega_i\right)$ incoming to point $x$ through an infinitesimal cone $\mathrm{d}\omega_{i}$ around the direction $\omega_i$ with solid angle size denoted $\mathrm{d}\sigma\left(\omega_{i}\right)$, which then generates some irradiance on the surface. This "partial" irradiance is then denoted $\mathrm{d}E\left(\omega_{i}\right)$. This incident light can be partially absorbed by the surface and partially scattered in all directions. The out-scattered contribution of $\mathrm{d}E\left(\omega_{i}\right)$ to radiance outgoing in direction $\omega_o$ is denoted as $\mathrm{d}L_{o}\left(\omega_{o}\right)$.
It can be shown experimentally that the amount of outgoing radiance $\mathrm{d}L_{o}\left(\omega_{o}\right)$ is proportional to the incoming irradiance $\mathrm{d}E\left(\omega_{i}\right)$ and the BRDF is then defined simply as the ratio between the two:
$$
f_{r}\left(\omega_{i}\rightarrow\omega_{o}\right) = \frac{\mathrm{d}L_{o}\left(\omega_{o}\right)}{\mathrm{d}E\left(\omega_{i}\right)}
$$


Perlin Noise, surface deformation NOT along the normal vector

You can do this, and the results can be interesting, but they’re pretty far from looking like realistic terrain. Here’s a plane deformed with 3D simplex noise (Perlin doesn’t look significantly better):

The issue is that there’s no volume to the surface, per se, and nothing preventing it from passing through itself. Deforming it only along the normal vector mostly avoids that. For a surface that looks like terrain, with actual caves, you’ll probably need some kind of voxel system like the one described in this GPU Gems article.


Compute Speed of Spheres in 3D animation

I'd use the euclidean distance between two 3D positions. Velocity is distance over time. Time is a single frame in this case.

$$\sqrt{(x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2}$$

If you only move the camera, it is simply that frame of reference that moves, so all spheres have a relative velocity equal to the camera's velocity.


How to get the length of paths in an given font/character

The way I would approach this is to get the outline curves of the font and then break each one down into some fixed number of steps (like, say, 100), and make a series of very small straight lines.

So if you've got some bezier curve, you loop t from 0 to 1 in 100 steps and plug it into the equation:

double x = calculateBezier(0.0, c0.x, c1.x, c2.x, c3.x);
double y = calculateBezier(0.0, c0.y, c1.y, c2.y, c3.y);
double dist = 0.0;
for (double t = 0.01; t <= 1.0; t+= 0.01)
{
    double nextX = calculateBezier(t, c0.x, c1.x, c2.x, c3.x);
    double nextY = calculateBezier(t, c0.y, c1.y, c2.y, c3.y);
    double nextDist = sqrt((nextX - x) * (nextX - x) + (nextY - y) * (nextY - y));
    dist += nextDist;
}

Where calculateBezier() is a function that calculates a point on a bezier curve given the parameter t and control points c0 through c3:

double calculateBezier(double t, double c0, double c1, double c2, double c3)
{
    double t1 = 1.0 - t;

    return (t1 * t1 * t1) * c0 + 3.0 * (t1 * t1) * c1 +
            3.0 * t1 * (t * t) * c2 + (t * t * t) * c3;
}

And just add the distances of each part of the curve until you're out of curves.

Note that this assumes cubic bezier curves. For quadratic, the calculateBezier function would calculate a quadratic one (like described here).


How does a computer upscale 1024x768 resolution to 1920x1080?

In essence a image is  a group of point samples (read A pixel is not a little square3). When you transform or scale the image you need to resample it. So what you do, theoretically, is take the point samples and convert them into a continuous function. Then you sample that continuous function and reconstruct a signal. So, there are two or three different phases here.

Converting the samples to a continuous function, (function reconstruction).
(Transforming)
Re-sampling the signal

Note that none of these steps have a fixed form. In practice, when optimized, it's impossible to tell that there are steps. The transformation does not really have to be simple it could be mapping the shape into a spiral etc.

Image 1: A 1-D signal reconstructed by different filters.

In practice, quite a bit of lore is known about signal reconstruction in the signal processing field. Designing these filters and choosing the right one is a art form on its own. But, in essence, the choice of filter is a tradeoff between blurring and ringing. Of course the algorithm also has other qualities such as how many instructions it takes to implement and how fast and how much memory it needs etc. Which can be very important in realtime or embedded applications.

Image 2: Overview of entire process.


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


Should ray throughput always decrease when bouncing?

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


How can I get a signed distance (SDF) from a mesh?

Generally it boils down to finding the nearest element. Computing distance to it. Checking the normal to tell if inside or out. There are a number of edge cases. Also you usually compute the "angle weighted psuedonormal" as its called.
Here's the usual reference "Signed Distance Computation using the Angle
Weighted Pseudo-normal" by J. Bærentzen.
There are many such libraries with the implementation if you don't want to write it yourself. I personally use ImplicitPolyDataDistance function in VTK. Also a fan of LibIGLs. But there are many lighter weight ones to do it as well, https://github.com/Magnus2/MeshSdf.
It's also worth mentioning that you can compute SDFs in other less accurate ways for more performance. Namely, another popular approach is to voxelize/binarize a shape and compute the SDF of the binary image instead. The distances at the boundary of the shape won't be exact but it is a good estimate and often faster.
EDIT: Also worth mentioning narrow bands and hierarchical SDFs. Those are very popular as well. Instead of computing SDF values for every element. Say in an image. We can compute SDF values only near the boundary of the shape. With an octree you can compute it at each vertex and interpolate to get inbetween estimates. There are special methods for computing these as well and libraries such as OpenVDB.


Polarising reflection and refraction for a wrap around water surface

The most commonly suggested method seems to be Mueller calculus, which boils down to tracking the Stokes parameters of a light ray to represent the polarization of light transmitted along that ray. A ray might be unpolarized—Stokes parameters of (1, 0, 0, 0)—or it may be circularly or linearly polarized in various directions, which is a property of light in the aggregate. At the surface light is scattered according to the polarization and the Stokes vector is propagated by multiplying it by the Mueller matrix of the surface.

Here's a writeup by Toshiya Hachisuka about ray tracing while tracking light polarization. It seems like a good introduction, and there are several references that seem promising. The article argues for direct tracking of the polarization state of the ray: instead of an aggregate representation, individually tracking the direction and frequency of the two harmonic oscillations of a given light ray. This may have the disadvantage that you need more samples to reproduce the polarization effects accurately, but it may be able to reproduce more effects (in the article, thin-film interference).


Can I have lower color bit depth in sample buffer than target framebuffer?

What does that even mean?

It means that multisampling is resolved.

A multisample framebuffer contains multiple samples for each pixel. The combination of samples represents the color of that pixel. The process of converting from multiple samples to a single color is multisample resolving.

By blitting from a multisampled buffer to a non-multisampled buffer, you are telling OpenGL to resolve each pixel's worth of samples into a single sample. Typically, this is done by adding all the sample values up and dividing by the number of samples, but implementations can take weighted averages or something else, depending on the precise nature of each sample.

Generally speaking, it's best that you let the implementation do its job.

  Is there a way to do it in a cross implementation defined manner?

You can program your own resolution by rendering with a multisample texture to a non-multisampled texture. In the FS, you read each fragment's worth of pixels, doing whatever operation you would like to combine them.

But that's probably not going to be as fast as letting the hardware do the multisample resolve itself.


Can V-Sync avoid tearing? I think it can't. Am I wrong?

V-sync means that the GPU won't switch buffer until the current frame is done showing.

When you tell to swap the buffer with V-Sync enabled the actual swap doesn't happen until V-Blank, a small window between frames where nothing is being sent to the screen.

The reason I think you cannot see the update tear is that it moves across the screen during the update. While a non v-sync tear stays in the same place on the screen for 16 milliseconds.


Closed form solution to a problem of tangent circles

First let us solve the radius of your circle. You can see that the centers of the circles form a triangle. There are 8 possible triangles, of which 4 are mirror solutions across the line that connects the centers of your original circles.  The other solutions are permutations of inside tangent and outside tangent connections.

Image 1: All possible solutions

So knowing the quadrant is hardly enough for the choice of solution.

Solving the triangles

So for eaxmple let the circle be on the inner outside of both circles (the sign changes for the other solutions). Now you need to solve the equation:

$$
 D^2=(x+r_1)^2+(x+r_2)^2-2(x+r_1)(x+r_2)\cos(\beta)
$$

Now this can be solved in a computer algebra system and gives 2 solutions with numerous corner constraints (do it yourself it's important for robust code and obviously you need 4 solutions to check.)

$$
x=\frac{1}{2} \left(\pm\sqrt{\frac{\cos (\beta) (r_1-r_2)^2-2 D^2+(r_1-r_2)^2}{\cos (\beta)-1}}-r_1-r_2\right)
$$

Once you have x solved you can solve the center. This can be done in many ways but I would proceed by solving angle $\alpha$ and some vector rotation, normalization and scalar multiplication. Same calculation can be used to find the intersection point and the tangent is just this vector rotated by 90 degrees.

You could also use an alternate formulation with 2 angles and a vector sum of $a+b-c=0$, with 2 unknowns and 2 equations.


What formula or algorithm can I use to draw a 3D Sphere without using OpenGL?

Scratchapixel has a nice tutorial on writing a basic rasterizer here. Also, you could use the projection algorithm here to get the position of the vertices in screen space, then use Bresenham's algorithm or DDA to draw lines in between. If you want to fill them too you can use scanline (you can find it on Wikipedia).

For ellipsoids, you can either just turn them into a triangle mesh, or the approach here might work, although I haven't tried it.

The Z-Buffer algorithm is very straightforward, just calculate the distance from the point to the camera and store that somewhere for each pixel, then make sure that that's less than whatever was already there.


How can we convert a probability density according to solid angle to a density according to area?

I don't understand what's the problem. We want to make a change of variables so that instead of integrating over the range of directions $d\omega_i$ in the upper hemisphere we integrate over patches of area $dA$ of a region $R$.

This whole thing is given in section 5.5 as Hubble said but let me post it again.

We know by definition that,

$\omega = A *  cos(\theta)/r^2$

Now if you differentiate this, you get

$d\omega = dA * cos(\theta) / r^2$

where,

$Jacobian = cos(\theta) / r^2$

as Hubble said. I don't know what are you trying to do by calculating Jacobian through determinant but as far as I know, that Jacobian determinant is for multi-variables and involves partial derivatives. For change of variables involving a single variable the Jacobian can be gotten by simply derivating as I did.

Now this might seem the inverse but wait a moment. What you posted in the question is a relation between the PDFs.

Consider a single sample of the rendering equation integrated over the domain of $\omega_i$

$L_s = BRDF * L_f * (\omega_i . n) / p(\omega)  $

After changing the domain by above equations this would become

$L_s = BRDF * L_f * (\omega_i . n) (\omega_i . n_r)\, / \,  (p(A) * r^2)  $

$L_s = BRDF * L_f * (\omega_i . n) * JN\, / \,  p(A) $
where $Jacobian = JN$

This means,
$ 1 / p(\omega) =  JN/p(A)$

$ p(A) = p(\omega) * JN$

$ p(A) = p(\omega) * cos(\theta)/r^2$

This is what you see above. Simple as ****.


Why integrate over a hemisphere (and not a sphere) to solve the rendering equation?

The form of the rendering equation that uses only the BRDF ($f$ in your example, often called $f_r$) and integrates over one hemisphere does not account for transmission.

When adding in transmission, it's fairly common to add a second integral over the opposite hemisphere, using a different BTDF function (bidirectional transmission distribution function). This is equivalent to an integral over the full sphere of directions with a BSDF function, but since that function would usually have to be defined as a piecewise function, writing it as two integrals can be more straightforward.


How do you apply a normal map to a 3d mesh?

There are numerous approaches to setting up tangent bases on a mesh, and unfortunately, no totally universal standard for how they are calculated.
Tangents are based on the mesh's UV mapping, so that the tangent vector points (at least roughly) along the U axis in texture space, and the bitangent along the V axis. This means each triangle in a mesh has its own tangent basis (as long as the UV mapping is nondegenerate). This often gets converted into a smoothed tangent basis per vertex, much like you can create smoothed normals by setting each vertex normal to an average over the normals of the surrounding triangles. Different tools and engines might do the smoothing in different ways, might handle UV seams and mirroring differently, might orthonormalize the basis per vertex or per pixel or not at all, etc. Here is a nice article to read if you want to know more. And as noted there, there is a method for tangent spaces called mikktspace that is gradually gaining widespread support and might become a de facto standard.
Object/character models like the one you linked usually have normal maps baked out from a content creation tool (e.g. based on an ultra-high-detail source mesh); ideally, you need the tangent basis you use in your renderer to match the tangent basis used in the tool to do the bake. Sometimes there will be vertex tangents stored in the mesh, but the .obj format doesn't support tangents so you would need to see if the model has tangents in one of the other formats, like .fbx.


A website devoted to prescanned photogrammetry is :Quixel.com They did (at least some of) the assets for the version 5 unreal engine demo.
I have gotten excellent results with software from Agisoft but it costs real money.
For good free software try 3DF Zephyr
This software will just generate the mesh, they will provide only minimal navigation of the result which is outside the scope of generation.
To navigate the mesh you can use something like windows 3D viewer.
Or better yet, and way more fun, write your own!


Why does the 1/r² term appear with point sources?

The concept of a point source is an approximation. Physically, light sources are extended objects and emit light from every point on their surface; but when you're far enough away (i.e. the distance to the source is large compared to its size) it's useful to approximate it as a point source.

You can get the $1/r^2$ law out of it as follows. Imagine a spherical area light with some radius $r_\text{light}$, and you're looking at it from a distance $r$ away. Then, we can approximate the solid angle that it subtends from your point of view as the area of a circle of radius $r_\text{light}/r$ (just using similar triangles). This area will be $\pi (r_\text{light}/r)^2$, so it's proportional to $1/r^2$.

Note that this approximation becomes exact in the limit $r_\text{light}/r \to 0$, i.e. when the light source is very far away or very small. It breaks down if the source is too large or too close.

If the source emits a constant radiance from every point on its surface, then when you integrate over solid angle in the rendering equation, you'll get a total irradiance proportional to $1/r^2$. In order to approximate it as a point source in a renderer, we skip the integration and just add an irradiance proportional to $1/r^2$ directly.


Should the alpha channel be gamma corrected

We know that in PNG,BMP,etc... the pixel value stored is not in the linear RGB space.

This is not necessarily true. You can store whatever color space you want into an image, it doesn't even need to be colors (such as normal maps).

The alpha channel is generally linear. The alpha channel doesn't get displayed, but it is generally a non-color term used for transparency (or whatever else). Because they don't need to display on a monitor, there's no reason to store in in gamma space. If you did, you would unnecessarily lose precision at the lower end of the alpha values. Normals maps follow a similar line of reasoning, as explained very well by Julien Guertault.


calculating size of rectangle which fully obscures a sphere

How long should $i$ and $j$ be? By similar triangles, they should have length $R\times (1-\frac{R}{\lVert V\rVert})$.

I believe this is where you went wrong. If you draw the 2D version of this on a sheet of paper, it will be obvious that the edge of your triangle which touches the side of your square goes through the sphere. More so as the camera comes closer.

That line needs to be tangent to the sphere. So your large similar triangle will have its right angle at the contact point with the sphere, not at the sphere's center as you did. If my trigonometry is not too rusty, this gives

$\frac{R}{\sqrt{\lVert V\rVert^2-R^2}} \times \left(\lVert V\rVert-R\right)$


Factor of 2 converting radiance to irradiance

I don't really understand what are you doing. I think your first equation shouldn't have a cos factor in it.

We have the relation,

$BRDF = dL_r / dE_i$

That is the brdf is the ratio of reflected radiance to incoming "irradiance". Re-arranging this gives us,

$ dL_r =  BRDF * dE_i$

For diffuse surfaces we know,

$BRDF = \alpha/\pi$

Substituting in above equation we have,

$dL_r = \alpha/\pi * dE_i$

Irrespective of what angle we view, the radiance is the same. Now for a diffuse surface with albedo = 1 we have what you said,

$dL_r = dE_i/\pi$

Now in order to measure Irradiance, we can use the relation,

$E = \int_{2\pi} L \cos(\theta) d\omega$

$E = \displaystyle\int_{\phi = 0}^{2\pi} \displaystyle\int_{\theta = 0}^{\pi/2} L \cos(\theta) \sin(\theta) d\theta d\phi$

$E = E/\pi * \pi$

$E = E$

So you got back the same irradiance.


Since this is linear interpolation it boils down to solving a linear equation
$y = a + b t$
You can substitute the following $y = y_s$, $a = y_1$(intercept) and $b = y_2 - y_1$ (slope) to get the following
$$y_s = y_1 + t(y_2 - y_1)$$
$$y_s - y_1 = t(y_2 - y_1)$$
$$\frac{y_s - y_1}{y_2 - y_1} = t$$
The same thing happens for points where you can write
$$D = A + (B-A)t$$
$$ (D-A) = (B-A)\cdot t$$
$$ (D-A)^2 = ((B-A)\cdot t)^2$$
$$ \sqrt{(D-A)^2} = \sqrt{((B-A)\cdot t)^2}$$
$$ \sqrt{(D-A)^2} = \sqrt{(B-A)^2 \cdot t^2}$$
$$ \sqrt{(D-A)^2} = \sqrt{(B-A)^2} \cdot\sqrt{ t^2}$$
$$ |D-A| = |B-A| \cdot t$$
$$ t = \frac{|D-A|}{|B-A|}$$


You might want to check matplotlib-cpp. You'll find a "funny-looking xkcd-styled example" in the README.
I also saw xkcd related entries on matplotlib's Python documentation. So, it should be available in Python too. However, I did not use these libraries myself. I just stumbled upon matplotlib-cpp while I was looking for a plotting library for C++.


Double rasterization work on triangles edges.

That's why the seam between two triangles can be rasterized twice.

The seam won't be rasterized twice because of the OpenGL or D3D fill/rasterisation rules.

What I think Humus is alluding to is that a fragment shader may be executed multiple times on pixels spanning the shared edge, which just implies wasted effort.

The moral of the story is "avoid small or long/thin triangles if you want efficient rendering".


It seems like you're asking two things.  I can't really speak technically about JBU, but I can give an overview of the necessary concepts and bilateral filtering generally.  You'll probably need to find more details yourself, but this should give a coherent structure to start from.

Fixing "Image"s
Many image-processing people view filtering as either something to be done as a post-process, or as a way to scale images.  The mathematician's perspective is more accurate, descriptive, and complex.
When you have an image file, you don't have a real image.  You have pixel "values", which are strictly valid only at the center of pixels.  Roughly speaking, each pixel has an area that takes on that color.  Mathematically, what you're doing is using a nearest-neighbor reconstruction filter to reconstruct the true image from pixel values.  So an "image" is really a set of samples and a reconstruction filter, typically nearest-neighbor.
Resampling
When you upscale an image, what you're really doing is resampling.  How do you do this?  You take your signal, and generate a new set of samples from it.  But for an image, I just said we don't have the real image; we have to reconstruct it first.  So a typical resampling procedure looks like this:

Take your pixel values.
Reconstruct an image between the values, which live at pixel centers, by using reconstruction filtering.
Measure a new set of samples from the reconstructed image, possibly spaced differently.
Throw away the old samples and use the new ones.

The key point is that the reconstruction filter used is the "filter" used for resampling.  For example, a "bicubic upsample" would use a bicubic reconstruction filter.
It turns out that if you use a sinc function for the reconstruction filter, and you don't shrink the image, you get (theoretically) the same image back.  In fact, if you upscale the image and then downscale it again, you'll get the same answer.  This isn't generally true for any other filter.
Most filters approximate a sinc filter to more or lesser degrees.  Deviations result in image loss.
Filtering
If you use the same before/after sampling rate, you do something called "filtering".  Formally, filtering is resampling without a change in sample rate.  Filtering is a special case of resampling.
The only possible purpose to filtering it to throw away information (e.g. a Gaussian) or possibly also to distort or add spurious information (e.g. a Sobel).
Bilateral Filtering
You can use many different filters for filtering.  Bilateral filtering is one of them.  I don't believe it's considered state-of-the-art for any purpose now (though I don't know what is; probably some newfangled machine-learning thing I expect).
A bilateral filter is a non-linear filter, meaning it doesn't have a nice Fourier representation.  Conceptually, it is like a Gaussian, but the source colors themselves also factor into a pixel's "distance" from others.  This tends to smooth flat areas while preserving edges.
Bilateral Upscaling
By now, I hope you see where I'm going with this: any filter you can use for filtering, you can use for upscaling or downscaling, which are really the same thing (because all of it is resampling).
If you've implemented a bilateral filter, you should be able to rather simply implement a bilateral filter resampler - just as, if you've implemented a Gaussian filter, you should be able to implement a Gaussian filter resampler.
The only difficulty may be the definitions.  Unfortunately, the distinction between filtering and resampling is systemically muddied to occasional non-existence in the image-processing literature.  For implementing a bilateral-type filter, I would rely heavily on its similarities to Gaussian resampling.  IIRC (I haven't read the paper in years) bilateral filtering is posed in terms of pixels, not samples.


how many watts are in the SPD of the color white?

There is no unique answer or approach and it's not simple either. I'm afraid that once you know the absolute spectral distribution of light that falls onto your camera, you've only computed half the story. What happens now depends on what observer you are trying to simulate, and to what accuracy.

I can highly recommend Moving Frostbite to Physically Based Rendering 3.0 from SIGGRAPH 2014. It contains many details especially section 4.3 and 5.1.

If you are trying to simulate a camera you need to answer some other questions for yourself:

What ISO am I simulating?
What is my shutter speed?
What is my aperture?

With the answers to that you can correct how much light falls onto your sensor and how sensitive it is to that light. Perhaps you don't have to correct for aperture if you already simulate this for depth-of-field.

Then knowing that you need to know how that sensor responds to light of different wavelengths. Perhaps you can find some spectral responsitivity diagrams for the camera sensor you are interested in simulating? Then you can divide by integrate for each sensor type (R, G, B) to get RAW exposure, which still needs processing to turn into sRGB (which involves choosing a whitepoint). It's not simple.

  I'm not an expert in color, and it's possible that my undestanding of the situation is wrong. Does it make sense to ask, what absolute(i.e. not relative) SPD of radiance should the pixel get(watts per steradian per meter cubed), so that its RGB value is exactly white (1, 1, 1) ?

Absolutely not. If you are indoors and look at a piece of paper, that is white, right? Now if you move outdoors, a much greater amount of light hits your eye, yet you still view it as white.

Assuming you view your image on a monitor, there's a lot of color spaces involved: the absolute spectral distribution, the camera, the stored image, the monitor used to display the image and your eyes. They all have different ideas of what 'white' means, and they all have different settings to tweak it: brightness, whitepoint, exposure, viewing conditions, etc.


What is the difference between a Sampler and an Image variables in GLSL?

Image and Samplers differ in many ways. Unlike Samplers, images can only be addressed with an integer coordinate that ranges from  (0,0) to (x,y), where x and y are the width and height of the image. They support random access reads like texelFetch, and support no filtering as you're accessing raw texels. The most interesting thing with images is that they support read and write access, in Open GLES there are limitations, wherein an image may only be specified as readonly or writeonly, not both, in desktop GL though you may specify an image as read/write so that you can both read and write to that image in a single pass.

Using a sampler, if you wanted to write a color into it you'd have to use an FBO, then a second pass to consume that. With images, you can do this all in one pass by writing to the image and reading from it, you'd have to handle the synchronisation yourself of course.
Images are very useful for image processing etc.


Can microfacet BRDF (including ggx) represent perfect mirror?

If this was correct, $\frac{D(H)G(L,V,H)}{4|N⋅V|}$ should be a dirac delta function. But I don't think so.

Actually, you're close to your answer - you just are trying to find out, what the assumption was in the very beginning - if you're talking about a specular BRDF. The term  $F\frac{DG}{4|N⋅V||N⋅L|}$ only works, if you start with the precondition, that you're dealing with a perfect mirror, otherwise this would be an integral itself. Have a look at PBR Diffuse Lighting for GGX+SmithMicrosurfaces by Earl Hammon Jr., where explains how to get to the specular BRDF. Especially have a look at slide 29:

Microfacet BRDF is a perfect mirror

I.e., light reflects if and only if $m = H$

Mathematically, BRDF is a scaled dirac delta $\delta_m(H, m)$

So the perfect mirror is your model in the beginning$^1$ - but only for the specular BRDF!

If you are trying to see, if including the diffuse BRDF into your BRDF is a perfect mirror, then this assumption does not hold any more and thus you don't have the dirac delta function.

$^1$ Note that he speaks about a microfacet being a perfect mirror. However, if all your microfacets are perfectly aligned, then your microgeometry "equals" your macrogeometry


Can I write histogram equalization, blurring, sharpening grayscale image etc.. as shader programs with openGL?

Yes you can. There are basically two approaches:

You draw everything to a framebuffer object. Then you apply image processing fragment shaders to resulting texture of the previous draw call, where the shaders take the pervious texture with a sampler.

or if you have a particular image and you would like to have more control on the numeric precision, the number of neighbours available to you while working on pixel etc:

You apply compute shaders directly to image.

For the first approach, see Learn OpenGL, for compute shaders, see Anton's tutorial


What kind of image processing was used in the Facebook "Celebrate Pride" filter?

The top Zuckerberg image looks like you first convert your background Zuckerberg image to grayscale. Then you can add a second layer with the rainbow texture but use the 'Screen' blending mode, this will cause it to colourise the layers beneath it, but in such a way the whites are preserved. You could also try 'Multiply' mode too if the results are too bright.

The second one with the facebook logo looks like alpha blending. This can be done using a default layer in photoshop but setting the layer opacity to some midway value, i'd guess around 75% opacity. In this version the rainbow texture is the background, and the layer contains the 'F' logo.


How different can the number of pixels in a straight line be to its real length?

The worst case is when the line is at a 45-degree angle; then its length (in pixel-sized units) is a factor of $\sqrt{2}$ times the number of pixels, or about 40% longer.

You can see this if you imagine starting with a vertical line of say 50px, then gradually sliding the endpoint away horizontally. Note that the number of pixels drawn doesn't change! You still have exactly one pixel turned on in each row that the line crosses. The pixels will adjust their horizontal positions but no new pixels will turn on, until the line surpasses 45 degrees from the vertical. The same is true if you start with a horizontal line and slide the endpoint vertically.

The number of pixels in the line effectively only looks at the vertical or horizontal distance, whichever is larger, and ignores the smaller dimension.
Mathematically, the number of pixels is $\max(|\Delta x|, |\Delta y|)$ while the Euclidean distance is $\sqrt{\Delta x^2 + \Delta y^2}$.


I fully understand the confusion (I was also a bit confused when I saw it). So, I decided to look at the ACES implementation of the ST 2084 EOTF. https://github.com/ampas/aces-dev/blob/v1.0.3/transforms/ctl/lib/ACESlib.Utilities_Color.ctl#L402
From there it is very clear: You perform the transfer function on the R, G and B values separately and not by first calculating the luminance. The confusion most likely comes from the fact that we are talking about two different kinds of luminance.
We have relative luminance, which is the luminance we use when talking about colours and can be calculated with $Y=0.2126R+0.7152G+0.0722B$
We also have luminance for which the SI unit is cd/m2 and basically means how bright something is.
You, and I when I first read your post, were thinking about relative luminance, but the Wikipedia article was talking about just luminance. In this case, it talks about the luminance of red light, green light and blue light separately (since a camera captures it separately and a monitor produces it separately). I'm pretty sure the ST 2048 spec is more clear about this, but that's behind a paywall.


How to work around the pointer limitations of GLSL?

I'll get to the specifics of your question presently, but first it is important to discuss the broad thesis of your idea: porting C++ code as is to GLSL.
This is a bad idea. The execution model of shaders is fundamentally different from the way a normal program (even a multi-threaded program) works. Attempting to port such code from one execution model to another is not a good idea. Your best bet is to try to make the algorithms involved fit the nature of the rendering pipeline. That is, you need to understand how rendering works, then try to apply the mechanics of rendering to the problem you're trying to solve.
For example, in regular CPU code, you might implement a ray-tracer by iterating over a bunch of locations on the screen and doing ray-scene intersection tests for each of them. Doing that iteration on the GPU in GLSL code is wrong. Instead, you should map rendering concepts to the needs of a ray-tracing algorithm. You might have each fragment shader do a single ray/scene processing step, with the fragment shader inputs being a direction from the camera. There would be no explicit loop over the screen's pixels; the loop happens implicitly in the rasterizer. The rasterizer generates a fragment shader invocation for each pixel-sized area of a rendered object. You would simply render a single quad that covers the whole screen.
Things like that are how you need to think of these kinds of problems in a GPU-centric computation model.

GLSL does not have pointers. But it does have shader storage buffer objects, which allow you to define arbitrarily large regions of memory that the shader can read from. A "pointer" in this model is just an index into an array of data structures. You have to do the indexing explicitly in GLSL code.
Similarly, any form of "late binding" can just be a switch statement based on some value in a data structure. The different casees define which function to call.


Mirror Reflections: Ray Tracing or Rasterisation?

There are a couple of special cases where mirror-like reflections can be rendered efficiently using rasterization techniques, and these are commonly used in games, although they don't work for the general case.

Planar reflections

If the reflecting surface is flat or reasonably close to flat, the reflected image can be rasterized in an separate rendering pass, by applying a reflection transform and clip plane to the scene. The reflection is then composited into the final image using a shader (or, in older games, stencil buffer masking). Also, a small amount of distortion or blurring might be applied to the image if the reflecting surface isn't intended to appear perfectly flat and smooth.

This is very commonly used in games for water reflection in lakes or the ocean, as long as any waves aren't too large. This page from the UE4 docs shows a number of examples of planar reflections used for water. The same technique is also often used for mirror reflections in glossy floors, or for actual mirrors.

The main strength of planar reflection is that it automatically gives you correct perspective and occlusion in the reflected image. The main limitation, aside from being limited to planar surfaces, is the performance cost of rendering an extra copy of the scene (with full lighting, etc) for each reflection plane. For instance, having multiple water surfaces at different elevations (e.g. in a boat lock) would require another rendering pass for each elevation. For this reason, usually when planar reflections are used, they're set up in such a way that only one reflection plane is visible at a time.

Cubemap reflections

Another special case is when the reflecting object is reasonably small and, ideally, not too close to other objects around it. In that case, you can approximate reflections by rasterizing a cubemap around the object. Then it can be used as an environment map during shading, by sampling the cubemap along the reflection vector calculated in the shader.

Many games use pre-rendered cubemaps to render static environment reflections on characters and dynamic objects; somewhat less commonly, you can also dynamically re-render a cubemap every frame to get dynamic reflections of other objects. This might be used on the cars in a racing game, for instance.

The advantage of cubemap reflections is that they can be applied to any shape of object, not just planar ones. However, since the cubemap is rendered from a single point, the reflections have increasingly incorrect perspective and occlusion the farther you get from that point. This can be partly addressed by parallax correction, and by blending between multiple cubemaps based on location, but ultimately cubemaps just don't give you the accuracy of either planar reflections or raytracing. Moreover, re-rendering a dynamic cubemap every frame is quite expensive, as it requires 6 extra rendering passes (one for each face of the cube).

Hybrid techniques

There are also cases where rasterization is used to build an initial data structure for ray tracing / casting to traverse into.

Screen-space reflection can be understood this way, as it involves first rendering the screen using ordinary rasterization techniques, then treating the Z-buffer as a heightfield and marching rays through it.
Voxel cone tracing is a similar case, where the scene is voxelized using rasterization, and then rays can be marched through the voxel grid. This perhaps doesn't qualify as "mirror-like" reflections since in practice the voxel grid isn't high-resolution enough for that, but it could be done in principle.

So yes, depending on what degree of approximation you can accept and what constraints you can work within, there are quite a few cases where mirror-like reflections can be achieved using rasterization, or hybrid rasterization/raytracing techniques.


Vulkan: using gl_PrimitiveID without geometryShader feature

It seems likely that the mechanism by which the gl_PrimitiveID is generated is, in some way, based on the hardware that makes geometry shaders possible.

OpenGL ES 3.1 contains compute shaders but not geometry shaders (or tessellation). But it doesn't have gl_PrimitiveID. By contrast, ES 3.2 contains GS's and gl_PrimitiveID. The same goes for desktop OpenGL 3.1 vs. 3.2 and its GS functionality. Having a GS and having the ability to identify primitives in a sequence seem to be part of the same general functionality.

  In my pipeline I used gl_PrimitiveID to fetch per-face data from global ssbo and to render face ids to attachment for later reuse.

Without the hardware functionality, there is no general way to compute a count of primitives. There are only specific solutions for specific circumstances.

For example, if you are using non-index rendering of triangle lists, then you can use gl_VertexID / 3 to compute a primitive ID. Obviously that's a pretty tight restriction, but it would work in that circumstance.

You could also provide each vertex with a per-vertex face ID. The VS would output this value using flat interpolation. By carefully matching the face ID with the provoking vertex for the primitive topology, you can effectively assign a face ID to a face. This works best with indexed triangle lists, as you can play around with the order of indices in a triangle more than you can with strips. You basically have to decide which vertex maps to which face, and make sure that that face uses that vertex as the first one.


Why does opengl perform clipping in clip coordinates?

The space after clip space is normalized device coordinate space, which is obtained by dividing clip.xyz by clip.w. If W is zero... oops.

Oh sure, you could put a conditional test to see if W is zero or very close to zero. Or you can just do the clipping math in clip space and never have to worry about it. No post-clipping vertex can have a W of zero, because the clip box for each vertex is based on being in the closed range (-W, W). And if W is 0, then the closed range is an empty set, and thus the vertex is not on that range.

Also, negative values of W are outside of the clipping space. The closed range (-W, W) is inverted where W is negative, but the meaning of the range itself is not. Consider a W of -1; the range becomes (1, -1). There are no numbers that are simultaneously greater than 1 and less than -1. Therefore that space is empty.

And yet, division by a negative W can still land NDC points in the [-1, 1] NDC-space range, where they will not be clipped if you didn't clip them beforehand.


It's not too difficult to set up your app to stream frames directly to ffmpeg. This blog post by Miles Macklin shows how, and I've used this method in my own apps.

You spawn ffmpeg as a child process from your app, and at the end of each frame you just read back the GL framebuffer and write its contents to ffmpeg's stdin. The relevant code bits are excerpted here.

Initializing ffmpeg:

// start ffmpeg telling it to expect raw rgba 720p-60hz frames
// -i - tells it to read frames from stdin
const char* cmd = "ffmpeg -r 60 -f rawvideo -pix_fmt rgba -s 1280x720 -i - "
                  "-threads 0 -preset fast -y -pix_fmt yuv420p -crf 21 -vf vflip output.mp4";

// open pipe to ffmpeg's stdin in binary write mode
FILE* ffmpeg = _popen(cmd, "wb");

Writing a frame:

glReadPixels(0, 0, width, height, GL_RGBA, GL_UNSIGNED_BYTE, buffer);
fwrite(buffer, sizeof(int)*width*height, 1, ffmpeg);

You can adjust the quality by altering the -crf option in the command line (lower numbers are better quality). Don't forget to change the resolution (-s) and framerate (-r) to match what you want, as well.

Depending on the quality, how fast your PC is, etc, it might not quite work in real-time. It's a good idea to set up your app to use a fixed per-frame timestep when you're in video capture mode, so the resulting animation doesn't depend on how much slowdown there was in the capture process.


I made a visualization a little like yours years ago, and used a similar technique, because it's the natural approach to the problem. After I spent a few years writing drivers for modern GPUs, I wouldn't use this approach at all, because it creates a dependency from one frame to the next, which introduces pipeline stalls and/or false sharing.

Today, I'd keep the path of the particle as part of its state, and render the whole trail each frame, instead of reusing the old FBO. It's conceptually a little more difficult, and you have to be careful that the way you draw the past path is exactly the same each frame. For example, if you represent the path with some kind of cubic spline, you can't just move the final control point forward, because it'll change the whole curve slightly. A circular buffer of (x, y) positions would be much easier to get right.

What you get in exchange for this extra complexity is (a) reduced frame time, memory use, and power consumption (which might be important if you're on mobile), and (b) you won't have your current problem at all. Maybe the frame time advantage won't be as big in WebGL, because you'll end up doing extra work in Javascript, which is slow, but memory bandwidth is more of a bottleneck on most systems. Just make sure that the co-ordinates you store are in the space of your simulation, and transform them into screen co-ordinates when you draw them, just like you would with a 3D scene.


What are the fConstX in Spherical Harmonics?

I'm not an expert on Spherical Harmonic (SH), but I'll try to answer.

SH functions are basis functions. That means that they encode a function $f$ this way by a given basis function $B_i$:

$$\int_{}^{} f(x) B_i(x) dx=fConst_i$$

So what we can do is to use the $fConst_i$ to reconstruct an approximation of the function $f$ this way:

$$\sum_{n=1}^{i} fConst_iB(i)$$

In your case they could look like this:

$$
fConst1 = 0.429043;
$$
$$
fConst2 = 0.511664;
$$
$$
fConst3 = 0.743125;
$$
$$
fConst4 = 0.886227;
$$
$$
fConst5 = 0.247708;
$$

These coefficients approximate diffuse lighting.


TL;DR: Your $G1$ formula is wrong.

Just to avoid confusion, I am assuming the isotropic version of the BRDF, the Smith microfacet model (as opposed to the V-cavity model), and the GGX microfacet distribution.

According to Heitz 2014, the masking/shadowing term $G1$ is

$$
\chi^{+}\left(\omega_{v}\cdot\omega_{m}\right)
\frac{2}{1+\sqrt{1+\alpha_{o}^{2}\tan^{2}\theta_{v}}}
$$

and according to Walter 2007, the formula is

$$
\chi^{+}\left(\frac{\omega_{v}\cdot\omega_{g}}{\omega_{v}\cdot\omega_{m}}\right)\frac{2}{1+\sqrt{1+\alpha^{2}\tan^{2}\theta_{v}}}
$$

Where $\omega_{m}$ is the microfacet normal direction (halfway vector), $\omega_{g}$ is the main (geometric) normal direction (normal), $\omega_{v}$ is the incoming or outgoing direction, $\alpha$ is the isotropic roughness parameter, and $\chi^{+}\left(a\right)$ is the positive characteristic function, or the Heaviside step function (equals one if $a>0$ and zero otherwise).

As you can notice, the half vector $\omega_{m}$ is used only to make sure the $G1$ is zero if the geometrical configuration is forbidden. More precisely, it makes sure that the back surface of the microsurface is never visible from the $\omega_{v}$ direction on the front side of the macrosurface and vice versa (the latter case is meaningful only when also refractions are supported). If the calling code guarantees this, then you can obviously omit this parameter. That is probably the reason why they did so in Unity.

Your implementation, on the other hand, uses the half vector to compute the cosine of the $\omega_{v}$ direction with respect to the microfacet, which leads to computation of something else than the presented formulae.

If it is of any help, then this is my implementation of the $G1$ factor:

float SmithMaskingFunctionGgx(
    const Vec3f &aDir,  // the direction to compute masking for (either incoming or outgoing)
    const Vec3f &aMicrofacetNormal,
    const float  aRoughnessAlpha)
{
    PG3_ASSERT_VEC3F_NORMALIZED(aDir);
    PG3_ASSERT_VEC3F_NORMALIZED(aMicrofacetNormal);
    PG3_ASSERT_FLOAT_NONNEGATIVE(aRoughnessAlpha);

    if (aMicrofacetNormal.z <= 0)
        return 0.0f;

    const float cosThetaVM = Dot(aDir, aMicrofacetNormal);
    if ((aDir.z * cosThetaVM) < 0.0f)
        return 0.0f; // up direction is below microfacet or vice versa

    const float roughnessAlphaSqr = aRoughnessAlpha * aRoughnessAlpha;
    const float tanThetaSqr = Geom::TanThetaSqr(aDir);
    const float root = std::sqrt(1.0f + roughnessAlphaSqr * tanThetaSqr);

    const float result = 2.0f / (1.0f + root);

    PG3_ASSERT_FLOAT_IN_RANGE(result, 0.0f, 1.0f);

    return result;
}


Finding the nearest neighbor of polygons

You're looking for a nearest neighbor search: https://en.wikipedia.org/wiki/Nearest_neighbor_search
We've solved similar problems in 3D space using spatial subdivision (kd tree). Depending on polygon count and uniformity, a linear search or hash grid would be options too. A linear search, while O(N^2) can still be faster than setting up extra data structures for small a small number of polygons/vertices. A hash grid should be faster to update than a kd tree if your scene is dynamic, and I would expect it to be more efficient when the scene is uniform in vertex density.


Perspective correct interpolation z-buffer

Yes, that's correct. Perspective-correct interpolation works by (for some quantity $u$ to be interpolated) calculating $u/z$ and $1/z$ at each vertex, linearly interpolating those values in screen space, then calculating $u = (u/z) / (1/z)$ at each sample point. This is done by the GPU hardware behind the scenes.


How to handle a parametric equation in raytracing?

Many years ago I worked on a ray tracer that handled parametric surfaces, so this is unlikely to be state of the art, but, IIRC, I used a combination of interval arithmetic with (binary?) subdivision and Newton-Rhapson.

The interval arithmetic + subdivision constructed (conservative) bounding boxes which could be used for intersection rejection. I think I may also have used interval arithmetic on the 1st derivative to help determine when it was safe to launch into Newton-Rhapson. Typically the Newton steps converged very rapidly.

You can also pre-dice your parametric surfaces and put them into an acceleration structure (Voxel grid/BVH) to speed up the process.

This work was based on/inspired by Daniel Toth's 1985 paper "On Ray Tracing Parametric Surfaces"


High variance in path tracing specular objects

If your bunny is purely specular, then sampling the light directly at the shade point would give no contribution since the specular BSDF is a delta BSDF. It generally evaluates to zero for any direction other than the mirror direction. If it was a glossy BSDF, then it might be possible that the pdf value could be very small so that the monte-carlo estimator $\frac{f(X)}{P(X)}$ just explodes (which is generally the source of fireflies i.e having a very small pdf and the function is large).

If it was a mirror direction, then you should just see a black bunny for direct illumination and just the reflection of the surrounding scene for multiple bounces. Multiple bounces from the specular bunny should not be cause fireflies, because if you're sampling based on cosine term at the matte surface, the estimator becomes $\frac{L_i(X)f_r(X)\cos(\theta_i)}{\text{pdf}(\theta_i)}$, where $\theta_i$ is the sampled direction using BSDF sampling at the matte surface. The pdf and the $\cos(\theta_i)$ should cancel out each other perfectly. However since you're using a point light, you can never sample from the BSDF since the light is a delta light and needs an explicit connection. Hence the pdf for the estimator would be 1 all the time. The only way for large values to creep up are in the numerator. So I would suggest you to check the values by increasing bounce one by one to see when the problems start creeping in.


In a physically based BRDF, what vector should be used to compute the Fresnel coefficient?

In Schlick's 1994 paper, "An Inexpensive Model for Physically-Based Rendering", where they derive the approximation, the formula is:
$$F_{\lambda}(u) = f_{\lambda} + (1 - f_{\lambda})(1 - u)^{5}$$
Where

So, to answer your first question, $\theta$ refers to the angle between the view vector and the half vector. Consider for a minute that the surface is a perfect mirror. So:
$$V \equiv reflect(V')$$
In this case:
$$N \equiv H$$
For microfacet-base BRDFs, the $D(h_{r})$ term refers to the statistical percentage of microfacet normals that are oriented towards $H$. Aka, what percentage of the incoming light will bounce in the outgoing direction.
As for why we use Fresnel in a BRDF, it has to do with the fact that a BRDF by itself is only a portion of the full BSDF. A BRDF attenuates the reflected portion of light and a BTDF attenuates the refracted. We use the Fresnel to calculate the amount of reflected vs. refracted light, so we can properly attenuate it with the BRDF and BTDF.
$$BSDF = BRDF + BTDF\\
$$
$$\begin{align*}
L_{\text{o}}(p, \omega_{\text{o}}) &=  L_{e}(p, \omega_{\text{o}}) \ + \  \int_{\Omega} BSDF * L_{\text{i}}(p, \omega_{\text{i}}) \left | \cos \theta_{\text{i}} \right | d\omega_{\text{i}} \\
 &= L_{e}(p, \omega_{\text{o}}) \ + \  \int_{\Omega} BRDF * L_{\text{i, reflected}}(p, \omega_{\text{i}}) \left | \cos \theta_{\text{i}} \right | d\omega_{\text{i}}  \ + \  \int_{\Omega} BTDF * L_{\text{i, refracted}}(p, \omega_{\text{i}}) * \left | \cos \theta_{\text{i}} \right | d\omega_{\text{i}}
\end{align*}$$
So, in summary, we use $D$ to get the percentage of light that will bounce in the outgoing direction, and $F$, to find out what percentage of the remaining light will reflect/refract. Both these use $H$, because that is the surface orientation that allows a mirror reflection between $V$ and $V'$


Description of the process involved in displaying digital images

Starting from reading the data from a file to displaying the image on the screen.

As was pointed out in the comments, this is an extremely broad question, and depends on what operating system you're using, and which tool sets for that OS you're using. That said, I think it is possible to talk generally about the process with much of it (if not all) being similar from machine to machine, and OS to OS.

I am puzzled regarding which parts of this process are handled by the OS and which are handled by the GPU.

I'll try to point out which parts are handled strictly by the OS, your program, and the graphics hardware, and I'll even provide C-like pseudo code. My pretend code will not be efficient, simply illustrative, and I'll be pretending to read this example image file.

For example, how does a 1920x1080 image eventually get displayed on screens with higher or lower display resolutions?

I think some of my pretend pseudo code will illustrate that, but the short answer is that you will manually resize your image to better fit the screen (or any other drawing area) inside your program.
The Long Answer
The first step is to open the image file. It doesn't matter what the format of the image is, from the view of your program it'll simply be a list of bytes. Some languages have tools and libraries that can be used to open files in a generic (OS independent) way, but often this step is OS specific. However, even though the task of opening and/or reading the file is dependent on the operating system, once your program has copied the content into memory it's just bytes for your program to deal with however you desire.
File imgFile = OS_OPEN_FILE("sample.ppm"); // Ask the OS for a file.
byte buffer[imgFile.size()];               // Set aside memory to hold the file's content.
for(int i = 0; i < imgFile.size(); ++i) {  // Start a reading loop.
    buffer[i] = OS_READ_BYTE(imgFile);     // Ask the OS for all bytes in the file, one by one.
}

Once you have access to all of the bytes in the entire image file it is up to you, the programmer, to know what each byte means, then to make your program read and interpret them correctly. For example, the portable bitmap format starts with a couple of bytes letting you know if the picture file is color, gray scale, or black and white. There are bytes that let you know what size the picture is, bytes that should be ignored because they're comments intended for humans to read, and finally bytes that are the dots of color to be sent to the screen. It is up to the program to sort through all of that.
int bufSize = imgFile.size();
for(int i = 0; i < bufSize; ++i) {  // loop through all the bytes read earlier.
    if(buffer[i] == 'P' && buffer[i+1] == '3') {  // Handle the image header.
        /* Do something */
    }
    else if(buffer[i] == '#') {  // Handle comments
        /* Do something */
    }
    else if(buffer[i] == /*...*/) {
        /* Do something */
    }
   /* ... etc. */
}

The next task is OS specific. You'll need to somehow ask the operating system to give your program access to a "drawing surface". Again, there are libraries such as SDL or SFML that provide surfaces in a generic way, but they are hiding the surface the system gave them behind an easy to use wrapper.
Often the drawing surface given by the OS is just a big block of bytes and you can simply copy data into it. Again, it is up to you to know which byte does what (like, which bytes represent red, green, or blue), but the task itself is usually just to copy from one chunk of memory into another. This is known as a Block Image Transfer, or a BLIT, and it's at this point that you may be handling memory inside the graphics card. Ultimately, it's up to the OS to decide where the drawing surface keeps the bytes that represent colors, but you can usually request that those bytes be on the GPU.
Surface sampleImage = OS_GET_GPU_DRAWING_SURFACE();
if(sampleImage.good()) {
    /* Copy in bytes from buffer[], sorting or converting data as needed. */
}

Depending on what tools you're using, and what task you're performing, a drawing surface you request from the OS can be the screen itself, and other times the screen is something that you can't touch directly but can only send surfaces to. Even when you can draw directly to the screen you will want to draw to a hidden surface, then blit the finished work over to the screen to prevent flickering and other visual problems. This is called double buffering, and mostly applies to the screen because it's always visible to the user.
Generally you will have many small drawing surfaces that represent things like, buttons, pieces of terrain in a game, letters, or whatever, and you'll use those to assemble a larger image on a bigger (usually screen sized) drawing surface that you will finally send to the screen as one big block.
Just to reiterate, these steps are all being done by your program, but the way your program has to do them is OS specific, and much of the memory being handled at this late stage is probably on the GPU.
// Ready the screen and a drawing area.
Surface * screen = OS_GET_SCREEN_SURFACE();
Surface * screenBuffer = OS_GET_GPU_DRAWING_SURFACE();

// Make the image 800% bigger!
sampleImage.resize(8.0f);
// Copy the image to the screen sized buffer, centering it.
// blitTo(surface, x, y)
// I'm assuming that in this pretend code images draw from their top left corner.
sampleImage.blitTo(screenBuffer,
                   screenBuffer.width()/2 - sampleImage.width()/2,
                   screenBuffer.height()/2 - sampleImage.height()/2);
/* ... etc. */
/* Do more image drawing work. */

// Finally, display all our complex behind the scenes work to the user.
// This is called a buffer swap.
Surface * temp = NULL;
temp = screen;
screen = screenBuffer;
screenBuffer = temp;
temp = NULL;

The Short Answer
In summary...
Task your program needs to do (what's primarily in charge):

Open a file (OS)
Read file (OS & Program)
Interpret data (Program)
Get a drawing surface (OS & GPU)
Fill the drawing surface with image data. (Program & GPU)
Assemble your final picture. (Program & GPU)
Gain access to the screen (OS)
Send your final picture to the screen. (OS and/or Program)


Path tracing: How to ensure we are sampling a direction vector within the visible hemisphere of a surface?

If it’s on the wrong side of the normal, you don’t have to throw it away—negating it will give you a vector that’s in the visible hemisphere.
Answering your comment, to get a vector that’s within some angle θ of the normal, this should work (GLSL):
vec3 direction(vec3 normal, vec2 randomValues, float maxTheta) {
    // pick an orthogonal tangent vector using the method described here: http://lolengine.net/blog/2013/09/21/picking-orthogonal-vector-combing-coconuts
    vec3 tangent = (normal.x > normal.z) ? vec3(-normal.y, normal.x, 0.0) : vec3(0.0, -normal.z, normal.y);
    // and a bitangent vector orthogonal to both
    vec3 bitangent = cross(normal, tangent);

    float phi = 2.0 * PI * randomValues.x;
    float theta = randomValues.y * maxTheta;
    float sinTheta = sin(theta);
    return tangent * cos(phi) * sinTheta + bitangent * sin(phi) * sinTheta + normal * cos(theta);
}


I dont understand how a GPU computes and branches, could you please help me fix this shader preformance

Rather than trying to speed up your shader by fiddling with the code, can I suggest you make use of the GPU's inherent ability to rapidly fill polygons? This will allow it to eliminate the vast bulk of the cases where there is no need to even make a decision.

What I suggest you try is something like the following: Create, say, a heptagon (i.e. the blue dashed outline (note: it shares some of the boundary with the green but I just noticed this is not particularly obvious in the diagram)) that starts from the original location and extends to the new position but lies inside the extremes of the new location's circle.
The shader for this is then utterly trivial.

Then there can be one or two cases depending on how much you want to optimise. If the tail extends outside of the end position circle, you draw the green-dashed nonagon with your decision-based shader. This contains vastly fewer pixels.
If the tail doesn't extend outside, then you'll need to mirror this polygon to get the other side of the circle. Alternatively, rather than make the decision, you could always mirror it and let the Z-test eliminate unneeded fills.

Your examples also have the have the capsules shown as opaque. If this is the case, you can also make gains with, say, TBDR-systems (a blatant PowerVR plug :-) ) or those with early-Z, provided you supply all your opaque "tail" portions first and then follow with all parts requiring decisions.

The shapes of course were just examples. You may find it's better to use more segments in your outlines in order to reduce the regions requiring the complex shader.


Will the gaussian kernels in this fragment shader be computed for every fragment?

It depends. On a desktop, if the Gaussian function is always called with the same arguments, the shader compiler will probably optimize it so it's evaluated only once. On a mobile platform though, the compiler has less time to do optimizations and might miss it. Either way, you can confirm by looking at the generated assembly.

Note that Aras Pranckevičius of Unity wrote a GLSL optimizer to address the problem of sub-optimal compilers on mobile: https://github.com/aras-p/glsl-optimizer

That being said, regarding your sample code, I would definitely at least expose (1.0 / sqrt(2.0 * 3.141592 * deviation)) as a uniform.


I'm not entirely sure, but I guess it might have to do with the notions of intrinsic vs extrinsic curvature. The Laplace–Beltrami operator on a manifold is, I think, intrinsic in the sense that it can be constructed from the "internal" geometry within the manifold itself, and doesn't depend on any embedding in an ambient space. The Steklov operator, on the other hand, looks like it's written explicitly in terms of an embedding, so would be extrinsic in that sense.


How to convert decomposed transformations between different coordinate systems?

Maybe not the best solution but one that should work would be to multiply them by matrix like
\begin{bmatrix}-1&0&0\\0&1&0\\0&0&-1\end{bmatrix}
And multiply by this your mvp matrix to convert meshes vertices (this could be done offline).
Quaternion can be change to matrix form, multiply and save as quaternion.
You also need to change culling in OpenGL because triangle order will be inverted.
I did this in one project and it works fine. If I will find this project I will write something more about it.
//EDIT
This matrix will invert your Z-axis and X-axis
//EDIT2
To apply this solution to each to component you need to:

Change translation x and z component to -x and -z.

Rotation in matrix form multiply by matrix above and change it to quaternion.

Scale leave how it is.
4. Iterate through all meshes and multiply each of their vertices by this matrix.

//EDIT3
Instead of matrix you can just multiply rotation by quaternion (0, 0, 1, 0) (in order (w, x, y, z)) and new quaternion will look like this: (-y,z,w,-x). This operation (and matrix too) just rotate object by 180 degree in y-axis.


How to subdivide faces with a half edge DS?

No one answered so I will. The solution is to do the subdivision in 2 passes, first generating new triangles as I showed in the question.

Then on the second pass you flip them.

As described there:
https://github.com/cmu462/Scotty3D/wiki/Loop-Subdivision


Is there a tool capable of drawing a triangular linear gradient fill?

Most 2D graphics programs are able to do linear gradients with arbitrary orientations. If you don't mind a little work, it is possible to set this up to imitate the 2D linear interpolation across a triangle.
Set up two layers:

Pick two of the colors, say red and green, and set up a linear gradient between those two colors along the red-green edge.
In a second layer, set up a linear gradient between blue and transparent-blue (i.e. the same color with alpha = 0). Orient it perpendicular to the red-green edge.