Radiance at a nearly perpendicular angle is almost infinite?

It's only more intense if $d\phi$ remains constant. However, for light source with constant luminance the flux is a function of $cos(\theta)$, i.e. depends on the light's angle of incidence with the surface and cancels out the term.

Edit: I feel the need to clarify this a bit because you say that:

  thus light is somehow more "intense"

It's not that the "light is somehow more intense", but rather that the light has to be more intense to have same flux in grazing angles as when it's parallel to the surface. If you think of flux as particles hitting to this differential surface, the number of particles hitting the surface is smaller more perpendicular the light gets to the surface. So you have to throw more particles at the surface (i.e. increase the radiance) for the same flux until when the surface is completely perpendicular to the light, none of the particles hit the surface, thus the radiance goes to infinity.


What effects does path tracing capture that recursive ray tracing does not?

Generally speaking, path tracing removes a number of assumptions that ray tracing makes. Ray tracing usually assumes that there is no indirect lighting (or that indirect lighting can be approximated by a constant function), because handling indirect lighting would require casting many additional rays whenever you shade an intersection point. Ray tracing usually assumes that there is no such thing as non-glossy reflection, because while it is fairly easy to handle glossy/mirror reflection (you reflect the ray and continue with the raycast), if you have non-glossy specular reflection you have a problem very much like indirect lighting: you need to trace a bunch more rays and then do lighting wherever those rays hit. Ray tracing also usually assumes that there are no area lights, because it is very cheap to cast a single ray towards a light to determine if an intersection point is in shadow, but if you have an area light you generally need to cast multiple rays towards different points on that light to determine the amount that the intersection point is shadowed.

And as soon as you need to break any of those assumptions with ray tracing, you now have a sampling problem: you need to decide how to distribute the many additional rays you need to trace in order to get an accurate result that doesn't have noticeable aliasing.


The variance for the Monte Carlo estimator

I don't have that book to check the context of this, but from the equations you posted, yes, it looks like you're right. The $1/N$ factor should be applied to both terms. That agrees with the formula for variance from statistics, which is $E[X^2] - E[X]^2$.


Confused about prefiltering environment maps (Manson 16)

Your first quote is referring to "Split-sum approximation" presented in "Real Shading in Unreal Engine 4" by Brian Karis, and also referred in the paper [Kar13]:

$$\frac{1}{N}\sum_{k=1}^N \frac{L_i(l_k)f(l_k,v)cos\theta_{l_k}}{p(l_k,v)}\approx \bigg(\frac{1}{N}\sum_{k=1}^N L_i(l_k)\bigg)\bigg(\frac{1}{N}\sum_{k=1}^N \frac{f(l_k,v)cos\theta_{l_k}}{p(l_k,v)}\bigg)$$

In split-sum approximation the incident radiance (cubemap) is pre-convolved with GGX kernel using for example importance sampling, and normal-incident specular reflectance ($F_0$ or "specular color" you refer to) is multiplied with another pre-convolved 2D look-up-texture (referred as "an auxiliary specular color map") shown below.

Next about your second question of importance sampling being inefficient due to band-limited signal - It simply means that the incident radiance is sampled with finite resolution (cubemap resolution) and thus higher frequencies have been filtered out. But because importance sampling biases samples to important regions of the function (peak of the GGX BRDF) this region may be sampled in unnecessarily high frequency thus wasting samples.


Compute queues in general don't necessarily mean you can now do 2x dispatches in parallel. A single queue that fully saturates the compute units will have better throughput. Multiple queues are useful if one queue consumes less resources (shared memory or registers), then secondary queues can then overlap on the same compute unit.

For real-time rendering this is especially the case with things like shadow rendering that are very light on compute/shaders but heavy on fixed function hardware, thus freeing the GPU scheduler to run the secondary queue async.

Also found this in the release notes. Don't know if it's the same issue, but could just be that CodeXL is not great. I would expect that it might have not have the best instrumentation for which dispatches are in flight.

https://developer.amd.com/wordpress/media/2013/02/AMD_CodeXL_Release_Notes.pdf

  For an application that performs simultaneous asynchronous data transfer and kernel execution,
  the timeline shown in the Application Trace session view will not show these operations being
  overlapped. This is because the driver and hardware force these operations to be synchronous
  while profiling. (333981)


Is there a way I can make this easier for openGL to draw?

You only need to create the VBOs and VAO once and then you can reuse them every frame.

void init(){

    glGenVertexArrays(1, &VAO1);
    glBindVertexArray(VAO1);

    glGenBuffers(1, &VBO1);
    glBindBuffer(GL_ARRAY_BUFFER, VBO1);
    glBufferData(GL_ARRAY_BUFFER, sizeof(data1), data1, GL_STATIC_DRAW);
    glBindBuffer(GL_ARRAY_BUFFER, VBO1);
    glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, 0);
    glEnableVertexAttribArray(0);

    glBindBuffer(GL_ARRAY_BUFFER, 0);

    glGenBuffers(1, &VBO2);
    glBindBuffer(GL_ARRAY_BUFFER, VBO2);
    glBufferData(GL_ARRAY_BUFFER, sizeof(data2), data2, GL_STATIC_DRAW);
    glVertexAttribPointer(1, 1, GL_FLOAT, GL_FALSE, 0, 0);
    glEnableVertexAttribArray(1);

    uniformLoc_centerX = glGetUniformLocation(program, "centerX");
    uniformLoc_centerY = glGetUniformLocation(program, "centerY");
    uniformLoc_sizeX = glGetUniformLocation(program, "sizeX");
    uniformLoc_sizeY = glGetUniformLocation(program, "sizeY");
    uniformLoc_size= glGetUniformLocation(program, "size");
}

void draw(){
    glUseProgram(program);
    glBindVertexArray(VAO1);

    glUniform1f(uniformLoc_centerX, centerX);
    glUniform1f(uniformLoc_centerY, centerY);
    glUniform1f(uniformLoc_sizeX, scrwidth);
    glUniform1f(uniformLoc_sizeY, scrheight);
    glUniform1f(uniformLoc_size, dotSize);

    glDrawArrays(GL_POINTS, 0, dot count);
}

If the data ever changes then all you need to do is:

glGenBuffers(1, &VBO1);
glBufferData(GL_ARRAY_BUFFER, sizeof(data1), data1, GL_STATIC_DRAW);
glBindBuffer(GL_ARRAY_BUFFER, VBO2);
glBufferData(GL_ARRAY_BUFFER, sizeof(data2), data2, GL_STATIC_DRAW);


For an nVidia only solution you can use floating point atomic add intrinsics (like NvInterlockedAddFp32) Unlocking GPU Intrinsics in HLSL | NVIDIA Developer

I tried this on 80.000 vertex mesh and it's quite fast (something like 1 or 2 ms on a GTX980M, if I remember correctly)

Just beware of compiling your shaders in release for the intrinsics to work (due to nvidia bug/limitation)

Also beware of vertex splits (due to UV discontinuities for example), you'll have to handle them or else you'll have unwanted hard edges at UV seams.


How to visualize density data as volume using OpenGL?

What you're looking for is "volumetric rendering", here is one resource with code that takes data in the form of images of the layers:
Getting started with Volume Rendering using OpenGL
A quick hack would be to use a grid of particles/"billboards", and make them transparent based on the density.


Getting contrast color to a given color in RGB or ARGB or HEX

I read this in a magazine article probably 25 years ago. I don't recall the full derivation, but they said to take each component (red, green, and blue) and invert the high bit of the color. In a C-like programming language, you could do:

 component = (~component & 0x80) | (component & 0x7F);

This can be simplified to:

component ^= 0x80;

for  a single component or to:

component ^= 0x80808000;

if you want to do 4 components at once (and don't affect the alpha).

This assumes your components are 8-bit values. This avoids the problems you get with the middle values if you simply invert them. In that case 127 becomes 128, and vice-versa. But with the above solution, 127 becomes 255 and 128 becomes 0.


Non-uniform rational B-spline (NURBS) Basics

B-Splines and Beziers are parallel inventions of more or less the same thing. Where Beziers try to start from the idea of fitting tangents. B-Splines start with the idea of basis functions. NURB Splines (or the rational part in fact) are just generalizations of B-Splines so you can describe accurate conic sections*, as they are of special interest in engineering.

First let us start with a simple NURB Spline terminology. The rationale of these curves is a bit different than for Beziers. First there is the concept of a span. A span would roughly be equivalent to a whole Bezier spline except in nurbs you can have any number of spans.

Image 1: One cubic NURBS span. This is a bit atypical in formulation

Each span is formed by the curve degree + 1  control points**. Each curve can consist of any number of points. Each consecutive span reuses the points form previous span by dropping one point and taking one point more in the list. So making more complex curves is as easy as just appending more points to the curve.

NOTE: The images curves are a bit atypically parametrized, ill explain what this means in the next section. When i take the concept of knots up. This is just a easier way to explain how the curves glue together.

Image 2: 2 cubic spans after each other, each span uses 4 points. together they form one curve. They share most points with each other.

By now we have probably answered the 2 question about adding complexity. But I would like to add that this scheme ensures better continuity than a bezier curve. Additionally you can make the point array that forms the hull cyclic. Forming a closed curve.

Image 3: A closed cubic NURBS surface has as many spans as it has points. Each color is one span.

Parametrization

Up until this point one could just say that the stringing together the spans is a trick just as "sewing" the Bezier curves. But there is a difference. The curve is parametrisized along its length. So the curves are not separate they do not interpolate form 0 to 1 on each span like Beziers do. Instead the underlying curve has a cusomizable parameter range. The parameter is stored in something called a knot, and each knot can have a arbitrary increasing value in the sequence. So you can parametrize the entire curves u range to 0 - 1 or 0 to 12. The parametrization also does not have to be uniform.

This parametrization changes how the curve is shaped. Why would this be useful? Well you can adjust the tension along the curve for one. Or you could encode the length of the curve into to the U parameter. One peculiar use is to make the NURBS curve act like a Bezier curve either fully or just partially (bezier like in the ends but not in the middle for example).

Image 4: Same points different knot sequences. The green NURBS curve corresponds to a Bezier curve that has a parameter range of 0-2 instead of 0-1

Ok so what are the knots? They are simply the ranges of the basis functions. Since the cubic b-spline with 4 points has 4 interpolating functions it needs 8 knots. Only areas where 3 functions overlap and sum up to 1.0 can a line be drawn.

Image 5: 2 different basis functions, a bezier like and a uniform segment parametrisation, spread to 0-1 range.

And now we have mostly described the answer to question 1. The range is not defined you can stretch the basis functions as you see fit. And finally the knot vector simply produces the parameter ranges for the basis functions. Theres still one more thing that governs the shape of the curve and that is the weight vector. But that another story to be told elsewhere.

* This rational in this case means that a NURBS curve does not have to be a polynomial, as you cant describe a circle with polynomials.

** One can define other types of points.


Perlin Noise with Smooth function vs Lerp

It looks like the way you're using smoothstep isn't quite right.

With lerp, the first two parameters are the endpoints of the output range, and the third is a 0–1 value specifying how far to interpolate. It looks like you're trying to use smoothstep the same way, but it doesn't work like that: the first two parameters to smoothstep are endpoints of the input range, and the third is a value in the input range; the output of smoothstep is always in 0–1. It's almost the inverse of lerp in that regard. (See also the GLM docs for smoothstep).

To use smoothstep here, you'd want to do something like this:

glm::lerp(bottomLeftDot, bottomRightDot, glm::smoothstep(0, 1, sx));

This applies smoothstep to the sx value (which is between 0 and 1). The result, which is another 0–1 value, is then used to lerp between the two "dot" values.


Generating and Combining Spectral Colors

Your way of calculating XYZ functions is probably the most efficient way to go about calculating accurate colors from a spectrum. It is standard practice afaik, for examples the books Physically Based Rendering (3rd) and Real-Time Rendering (3rd) both use this method.

You can add the colors in RGB space, but only if you convert from sRGB to linear RGB first. Otherwise you need to take into account, that sRGB sums will not lead to correct colors. The blogpost Adventures with Gamma-Correct Rendering by Naty Hoffman is a good read regarding this, topc:

  Computing shading in sRGB space is like doing math in a world where 1+1=3.

This is the problem you have when adding colors in sRGB.

As to whether $XYZ$ can be summed, I think so. If we look at the definition of the coordinates (for $S$ is the spectral function and $x, y, z$ are the CIE functions):

$X = \frac{1}{\int_\lambda y(\lambda) d\lambda} \int_\lambda S(\lambda) x(\lambda) d\lambda$

Now two $X$ coordinates $X_1, X_2$ would have the same definition, only the spectrals are different ($S_1, S_2$).

Thus you have

$X_1 + X_2 = \frac{1}{\int_\lambda y(\lambda) d\lambda} \int_\lambda S_1(\lambda) x(\lambda) d\lambda + \frac{1}{\int_\lambda y(\lambda) d\lambda} \int_\lambda S_2(\lambda) x(\lambda) d\lambda\\
= \frac{1}{\int_\lambda y(\lambda) d\lambda} \left(\int_\lambda S_1(\lambda) x(\lambda) d\lambda + \int_\lambda S_2(\lambda) x(\lambda) d\lambda\right)\\
= \frac{1}{\int_\lambda y(\lambda) d\lambda}\int_\lambda S_1(\lambda) x(\lambda) + S_2(\lambda) x(\lambda) d\lambda \\
= \frac{1}{\int_\lambda y(\lambda) d\lambda}\int_\lambda (S_1(\lambda)  + S_2(\lambda)) x(\lambda) d\lambda\\
=X$

The same holds for other coordinates.

Therefore, your best option is probably to sum your XYZ coordinates and convert them to sRGB space in the end.


In convolution, two mathematical functions are combined to produce a third function. In image processing functions are usually called kernels. A kernel is nothing more than a (square) array of pixels (a small image so to speak). Usually, the values in the kernel add up to one. This is to make sure no energy is added or removed from the image after the operation.

Specifically, a Gaussian kernel (used for Gaussian blur) is a square array of pixels where the pixel values correspond to the values of a Gaussian curve (in 2D).

Each pixel in the image gets multiplied by the Gaussian kernel. This is done by placing the center pixel of the kernel on the image pixel and multiplying the values in the original image with the pixels in the kernel that overlap. The values resulting from these multiplications are added up and that result is used for the value at the destination pixel. Looking at the image, you would multiply the value at (0,0) in the input array by the value at (i) in the kernel array, the value at (1,0) in the input array by the value at (h) in the kernel array, and so on. and then add all these values to get the value for (1,1) at the output image.

To answer your second question first, the larger the kernel, the more expensive the operation. So, the larger the radius of the blur, the longer the operation will take.

To answer your first question, as explained above, convolution can be done by multiplying each input pixel with the entire kernel. However, if the kernel is symmetrical (which a Gaussian kernel is) you can also multiply each axis (x and y) independently, which will decrease the total number of multiplications. In proper mathematical terms, if a matrix is separable it can be decomposed into (M×1) and (1×N) matrices. For the Gaussian kernel above this means you can also use the following kernels:

$$\frac1{256}\cdot\begin{bmatrix}
1&4&6&4&1\\
4&16&24&16&4\\
6&24&36&24&6\\
4&16&24&16&4\\
1&4&6&4&1
\end{bmatrix}
=
\frac1{256}\cdot\begin{bmatrix}
1\\4\\6\\4\\1
\end{bmatrix}\cdot\begin{bmatrix}
1&4&6&4&1
\end{bmatrix}
$$

You would now multiply each pixel in the input image with both kernels and add the resulting values to get the value for the output pixel.

For more information on how to see if a kernel is separable, follow this link.

Edit: the two kernels shown above use slightly different values. This is because the (sigma) parameter used for the Gaussian curve to create these kernels were slightly different in both cases. For an explanation on which parameters influence the shape of the Gaussian curve and thus the values in the kernel follow this link

Edit: in the second image above it says the kernel that is used is flipped. This of course only makes any difference if the kernel you use is not symmetric. The reason why you need to flip the kernel has to do with the mathematical properties of the convolution operation (see link for a more in depth explanation on convolution). Simply put: if you would not flip the kernel, the result of the convolution operation will be flipped. By flipping the kernel, you get the correct result.


Converting Vertex and Index buffers to a File Format

No idea about DirectX, but in OpenGL the buffers are just a block of binary data so you can easily load and save them.

In any 3D framework, changing the object position to achieve a different relative arrangement of a scene is not done through the view matrix (that represents the camera position and viewing direction) but through the model matrix of each model. Note that the vertex and index buffers stay unaffected by matrix transformations, that is the main benefit of such pipelines.

Perhaps you could do your study with less effort using a 3D game engine, to list just a few:

https://godotengine.org/
https://www.unrealengine.com/
https://unity.com/

I am not aware that any of them would support loading the raw buffers straight away, but they can load OBJ files. OBJ is one of the oldest text-based 3D model formats. It has a simple structure, so you can easily convert your buffer data to an OBJ file.


Is there some book about how to design a 3D math lib?

Books: http://www.essentialmath.com/book.htm is the one.

Best algos on the real 3D topics (and some tricky 2D like image resize and premult-alpha ), explained and with source most of the time.

The website is very nice too, the "tutorial" full of gems

Adding http://www.geometrictools.com/index.html which is the most complete source code repository on the topic, with loads of very precise algos implemented.


Is there different between 'if' and 'ternary operator' in GLSL?

The difference is obvious. The ternary operator will return one of the two alternative values it is given; that's what it is for. An if condition that doesn't have an else block won't do anything if the condition is false.
The ternary operator version will set color to axesColor or background. But by that point in your code, you've gone past 3 other conditions that may have changed the value of color from its initial value of background. But the ternary operator doesn't care. It will overwrite whatever is in color with either axesColor or background.
What you probably wanted was color = abs(coord.y) < 0.005 ? axesColor : color;, so that the value will be preserved if the condition is false.


By "bounces by the BRDF", I assume you mean picking random directions in the hemisphere and weighting by the BRDF, then averaging over those samples. I guess you're accumulating one sample per frame, so then you divide by the current number of frames.

To be theoretically "correct", you should also be weighting those samples by $2\pi$, as that's the solid angle of the hemisphere you're integrating over. That's not too important when you just have a single sampling method, but when you're going to combine multiple ones it's easier to be sure everything's consistent if you have these factors right.

Then, to incorporate explicit direct light sampling, you have to convert from the light-sampling probability density to the solid angle probability density at the receiver point. If the light samples are uniform over the light source's area, then in place of the simple $2\pi$ weight, you'd use this weight factor instead:
$$\frac{A}{r^2} (N_\text{light} \cdot -L)$$
where $A$ is the area of the light source, $r$ is the distance between the receiver point (current path point) and the chosen light source point, $N_\text{light}$ is the normal of the light point, and $L$ is the unit vector from the path point, toward the light point.

This factor accounts for how dense the sampled light source points will be in solid angle of the receiver point, varying both with distance and with the angle on the light source's surface (the points will bunch up in areas where the source surface slopes away from the receiver, such as around the edges of a spherical area light).

After weighting by this factor, the direct light samples are in the same solid angle domain as the BRDF samples, so you can just add them together and divide by the total number of samples—which would now be twice the number of frames, if you're doing one BRDF sample and one light sample per frame.

Also, an important thing to get right is not double-counting the light's illumination. Since you're explicitly sampling the lights now, you have to make sure not to include the emissive color if you do happen to randomly hit a light while doing the BRDF samples! Otherwise you'd end up with the lights appearing twice as bright as they should.

A refinement of this is that you may wish to apply explicit light sampling only to the diffuse part of the BRDF, and not the specular part. A narrow specular highlight is better handled through importance-sampling the BRDF instead. This implies that when doing explicit light samples, you would only evaluate the diffuse component of the BRDF, and then if you hit a light source through a specular reflection ray, you would include its emissive color there.


can half of a picture be processed and viewed?

Yes, most formats work this way! For example JPEG can be decoded up to the point where the chunks are available:

Image 1 Scaled version of original JPEG picture at (1 048 198 bytes) and corrupted one (610 922 bytes) side by side. Original image by Maruf Mostafa image available here

This is easy for you to test take an image and just simply delete some bytes out of it (like with a text editor). Some files are easy just open them, like the JPEG in question that just opens up in many editors without a hitch.

PNG required some slight special handling to uncompress. It did't open in Photoshop or windows image preview. But was still relatively straightforward to open as the web browser happily opened even partial PNG files to the point it was defined to.


Q1: how to represent path tracing as the first one...
Your materials have brought enough information, the first formula, quoted from your material, presents as

The path integral formulation of light transport formalizes the idea
by writing the camera response as an integral over all light transport
paths of all lengths in the scene, which given path encompasses

the “amount” of light emitted along the path,
the light carrying capacity of the path, and
the sensitivity of the sensor to light brought along the path.

Can't interpret the formula better
Q2: Should I generate paths for different lengths and do this procedure over and over...
Q3: ...randomly generate paths with different lengths and sum their estimate up since each length...
Quoted from your material,

The Monte Carlo integration procedure consists in generating a
‘sample’, that is, a random x-value from the integration domain, ...

I can't say the second one is a formula, it is just the Monte Carlo estimator of the path integral formulation of light transport or it is just a method for programming the formulation; or course you could have other methods.


Before we begin, let us differentiate between two things:

The shape of a sphere,
and the topology of a sphere.

A NURBS surface can make the shape of a sphere. In a typical configuration, it will be 'open' at the poles. That is the mathematical function of the surface does not wrap over the pole in that it is not a true sphere (it has the shape of a sphere). Topologically the sphere is in this typical configuration a cylinder.

It is possible to leave the other direction open too in which case it is topologically a plane. Several CAD applications choose this approach. A torus is also possible if you allow a thin no volume sliver at the center of your sphere.

Image 1: Turning cylinder to sphere. Note this is not a topological sphere it is still a cylinder as the top is open (even though infinitesimally small).

  Is the sphere homeomorphic to one of the above surfaces?

No, but you can still have a spherical shape even if it does not satisfy the topology condition of mathematics.

Why only 3 topological families?

Simply, a NURBS surface has only 4 possible configurations of wrapping around the parameter space:

It does not wrap at all. Topology: Plane.
It wraps around the U direction. Topology: Cylinder.
It wraps around V direction. Same as above. Topology: Cylinder.
It wraps around both U and V. Topology: Torus.

A wrap is always periodic so it goes from - direction to + direction. It can not arbitrarily connect (on a mathematical level).


Weird glColor3f results in framebuffer objects

You need to set the colour back to white before drawing the framebuffer to the screen. Otherwise, the whole framebuffer gets multiplied with red.

glBindFramebuffer(GL_FRAMEBUFFER, 0);

glColor3f(1, 1, 1); // <------- add this (doesn't have to be exactly here)

glBindTexture(GL_TEXTURE_2D, fbo_tex[0]);
glMatrixMode(GL_PROJECTION); glLoadIdentity(); glOrtho(0, SCRWIDTH, SCRHEIGHT, 0, -1, 1);
glMatrixMode(GL_MODELVIEW); glLoadIdentity();
{
    glBegin(GL_QUADS);
    glTexCoord2f(0, 0); glVertex3i(0, 0, 0);
    glTexCoord2f(0, 1); glVertex3i(0, SCRHEIGHT, 0);
    glTexCoord2f(1, 1); glVertex3i(SCRWIDTH, SCRHEIGHT, 0);
    glTexCoord2f(1, 0); glVertex3i(SCRWIDTH, 0, 0);
    glEnd();
}
glBindTexture(GL_TEXTURE_2D, 0);


In what order to apply lighting and translucency?

When blending multiple layers, physically the "right" thing to do is calculating lighting on each layer separately, then composite the lit layers together. This way, for instance, you can have a translucent sprite standing in a spotlight, in front of a dark background, and the sprite will be well-lit while the background visible through it stays dark. Or conversely, you can have a dark sprite silhouetted against a brightly lit background.

If you blended the sprite with the background before lighting, it wouldn't be clear which lighting environment you should use if the sprite and background have different ones.


how to put a pixel somewhere in a window(directx)

there is no direct way to place a pixel on a window I am afraid.
An indirect way would be to draw a small quad equal to the size of a pixel in relation to your window
and essential y d3d 11 and 12 are most useful to draw geometry
So I recommend you start from there and take your creativity ahead


Your example images do not make the desired effect very apparent, but the opposite of a barrel distortion is a pincushion distortion:

I found a fisheye shader on ShaderToy which produces the following image:

By slightly modifying the shader code:

// This shader is based on cafe's fisheye shader:
// https://www.shadertoy.com/view/ll2GWV

vec3 checker(vec2 uv){
    return vec3(abs(floor(mod(uv.x*10.,2.))-floor(mod(uv.y*10.,2.))));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = fragCoord.xy*2. / iResolution.xy-vec2(1.);

    //------------------------------------------------
    // To use in Godot, port this section:
    //------------------------------------------------

    // I picked these somewhat arbitrarily
    const float BARREL = -1.0;
    const float PINCUSHION = 0.1;

    float effect = PINCUSHION; // Set effect to either BARREL or PINCUSHION
    float effect_scale = 1.0;  // Play with this to slightly vary the results

    /// Fisheye Distortion ///
    float d=length(uv);
    float z = sqrt(1.0 + d * d * effect);
    float r = atan(d, z) / 3.14159;
    r *= effect_scale;
    float phi = atan(uv.y, uv.x);

    uv = vec2(r*cos(phi)+.5,r*sin(phi)+.5);

    //------------------------------------------------
    // end relevant logic
    //------------------------------------------------

    fragColor = vec4(checker(uv),1.);
}

We can easily obtain the opposite effect:

If this is not the exact effect you're looking for, I recommend playing around with the code a bit yourself. Change some signs from + to -, add scale factors and explore the possibilities. You may even stumble across a more interesting effect than what were you looking for.


If you are writing your own Vertex/Fragment Shader there is another additional possibility to do this. It is much more complicated but might be useful for you and/or others. Additionally it speeds up the whole drawing process since it uses only one call to a draw command. The maximum number of Viewports is defined by GL_MAX_VIEWPORTS and is usually 16.

Since OpenGL 4.1 the method glViewportArrayv exists. It can define an array of Viewports. The Viewports which are created by this method have an index assigned.

This index can be used in the Vertex Shader to set the Viewport to which the scene is rendered. (You have to include the "ARB_shader_viewport_layer_array" extension in the shader code)

In your case I would suggest to do the following:

Store the 4 camera matrices in a Shader Storage Buffer (Array of mat4) to have them in the Vertex Shader
Use indexed drawing, for example :glDrawElementsInstanced
use the build in gl_InstanceID of the Vertex Shader to access the camera matrix array
set the build in variable output variable gl_ViewportIndex in the Fragment Shader to the gl_InstanceID. (for details see ARB_fragment_layer_viewport)


How much precision do I need in my G-Buffer?

First of all, you don't need position in the G-buffer at all. The position of a pixel can be reconstructed from the depth buffer, knowing the camera setup and the pixel's screen-space xy position. So you can get rid of that whole buffer.

Also, you don't ordinarily need tangent vectors in the G-buffer either. They're only needed for converting normal maps from tangent space, and for parallax mapping; these would be done during the G-buffer fill pass (when you have tangents from the mesh you're rendering), and the G-buffer would only store normals in world or view space.

Material properties like colors, roughness, and metallic are usually just 8-bit values in the G-buffer, since they're sourced from 8-bit textures. Same for AO.

Height is also not needed in the G-buffer unless you're going to be doing some kind of multi-pass blending that depends on it, but if you do need it, 8 bits is probably enough for that too.

Normals can be benefit from being stored as 16-bit values rather than 8-bit. Half-float is okay, but 16-bit fixed-point is even better, as it gives you more uniform precision across all orientations (half-float is more precise near the axes and loses some precision away from them). Moreover, you can cut them from 3 components down to 2 using octahedral mapping.

So, at the end of the day, a minimal G-buffer might look like:

Material color + metallic: RGBA8
Octahedral world-space normal + roughness + AO: RGBA16

and that's all! Only 12 bytes per pixel.

Alternatively, you could use an RG16 buffer for the normals, and move roughness + AO into a separate 8-bit buffer. That would give you some room to grow should you eventually need more G-buffer components of either 8-bit or 16-bit sizes.


3D projection that increases objects' size as they become more distant

After a quick google of 'Inverted perspective' I found out that you have five different names for it; Reverse perspective, inverse perspective, inverted perspective, divergent perspective and Byzantine perspective. https://en.wikipedia.org/wiki/Reverse_perspective

I think that if you want to use it in a render engine using a projection matrix, in the perspective divide, you use multiply by w instead of divide. This is because in a perspective projection, the divide by w is what actually gives the perspective. In a orthographic projection w is always 1, meaning that there is no actual perspective.

For a ray tracer you probably just need to scale the image plane up of course, but probably instead of having the origin of the ray to be the center of the camera, you have some kind of second image plane that is larger than the actual image plane. If both are the same size, the rays are straight and you have an orthographic projection, if the actual image plane is smaller then you have your Reverse perspective.

Quick extra note: If you have proper lenses, you might be able to get this effect in real life. Yes, you would need to have a lens that is a few meters in diameter, but I do not think it is impossible! :P

I hope this will help you and that I have not made any mistakes. Good luck with whatever you are doing!


Is BVH faster than the octree/kd-tree for raytracing the objects on a GPU?

First, it is important to observe that there are basically two different "classes" of methods for building accelerations structures for ray tracing: 1) space subdivision methods (e.g. BSP, kD-Tree, Octree, etc.); and 2) object subdivision methods (e.g. BVH). While space subdivison methods work by recursively subdividing the space with planes, without sticking to the geometry embedded into that space, object subdivison methods work by subdividing the geometry recursively into smaller and smaller parts, wrapping each part with a, usually tight, volume. More on the basics of space and object subdivision that can be found in the PBRT v3 online book.
Each method (space or object sudivision) has its own pros and cons. In the case of space subdivision methods, since the subspaces do not overlap, we can usually traverse the structure in front-to-back, or back-to-front, order more easily. When a ray is traversing such structure, as soon as it hits a surface, we can stop the traversal. This usually leads to faster traversal schemes. Several software renderers take advantage of the traversal efficiency given by the space subdivision schemes. On the other hand, space subdivision schemes may be more intricate to implement (usually you have to tweak some epsilons) and may lead to deeper trees. Also, they do not like very much dynamic geometry. If the geometry encoded into a space subdivision acceleration structure changes, we usually have to rebuild the acceleration structure from scratch.
Object subdivision methods have quite different characteristics. Since the object is subdivided with volumes, and these volumes may overlap, the traversal is traditionally slower. We cannot, for instance, stop traversing a BVH as soon as a ray finds an intersection with a surface. Since the volumes overlap, we may need to check for potential intersections with nearby primitives before quitting the traversal. On the other hand, it may be easier to implement a BVH because we do not have to split the object parts with planes. Also, BVHs usually generate shallower structures (which may eventually compensate the slower traversal). However, one of the most interesting aspects of the BVHs is that they are dynamic geometry-friendly. If geometry changes (but not much, actually), we can simply locally adjust the size and position of the corresponding bounding volume (i.e. refitting). These adjustments may cause the need to adjust the parent volumes, a procedure that may culminate in a chain reaction that may reach the root node of the BVH. All in all, assuming that we have a reasonably balanced BVH, these operations will be ~O(log n), which is really fast and cool.
The above paragraph describes the traditional characteristics of the BVHs. There is currently lots of research going on on the topic, and several of its ineficiencies have been significantly reduced over time. It should be noted, for instance, that still in 2007, Wald, Boulos and Shirley have designed a very efficient BVH-based acceleration structure, for interactive ray tracing on GPUs, that rivals with space subdivision schemes such as kD-Trees. This work was very influential.
Actually there are several reasons that make the BVHs the preferred acceleration data structure for interactive ray tracing o GPUs. According to Wodniok and Goesele, the main factors that have turned the BVH into a popular acceleration structure in the context of GPU ray tracing are its "memory footprint (...) efficient empty space cut-off, fast construction, and simple update procedure (...) while (...) offering similar (...) performance as kD-Trees". Regarding fast BVH reconstruction, it may be interesting to take a look at the LBVH, by Lauterbach et. al., the HLBVH, by Pantaleoni
and Luebke, and the one by Karras and Aila, that are capable to rebuild BVHs iteractively. Overlapping impacts the traversal procedure efficiency because, when a ray hit the overlap, it has to be tested against the two child nodes. SBVHs, by Stich et. al. is a hybrid structure that employs both spatial and object subdivisions such that overlap is reduced. The cons, according to the paper, are the fact that SBVH does not handle dynamic scenes anymore, and that it may increase memory consumption up to 30%. The SBVH building procedure, also, is not as fast as for the regular BVHs.
The current NVidia RTX GPUs support the traditional BVH, the SBVH and Karras/Aila acceleration structure. As pointed out and detailed by @DirectX_Programmer, in his answer below, these structures are actually combined, by the RTX drivers, into a more flexible arrangement composed by two levels.


Given a background image B and a composie image of A OVER B, can we somehow recover image A?

Just to confirm, assuming non-premultiplied alpha, at one point you performed the following, per pixel operation:
$$ Result_{RGB} = Src_{A} * Src_{RGB} + (1.0 - Src_{A})* Background_{RGB} $$

You say you've lost the $Src_{ARBG}$ data,  but have retained the $Background_{RGB}$ and the $Result_{RGB}$, but want to recover the $Src$ data?

Unless you can make some other assumptions, you might be out of luck. It's a rather underconstrained problem.
Jim Blinn notes in "Fun with Premultiplied Alpha" that there are some commercial systems that are used to 'de-bluescreen' but they may require the operator to identify, say, what is 'Source' and 'Background'.

As an extreme example, you may have a pixel in the result which exactly matches the background but the source pixel may also have been the same colour.

If, on the other hand, $Result$ and $Background$ don't match, then you know that the source is involved but if the $Src$, $Result$, $Background$ RGB values are all collinear, then there could be multiple possible solutions.

Having said this, there may be some heuristics that perhaps could be used, but it'd require some thinking.
Perhaps, pixels which differ and are not collinear could be used to estimate $Src_A$, with some level of 'certainty', and this could be 'filtered' to neighbouring 'uncertain' pixels and the process repeated.
Non-constant pixels in the destination with a local neighbourhood which all match the background might be assumed to have an alpha of 0.0.


It's generally referred to as a "state machine". It has, as you noticed, a bunch of switches and parameters. All of them are globally accessible to any code that has access to the current context. Your program sets them and then issues draw calls. Anything you don't change, will remain in that state until you do change it.

It is also a Hardware Abstraction. The same OpenGL routines should have roughly the same effect regardless of the hardware it's running on. If you upload 3 vertices for a triangle and then draw them, that should produce similar results on an ATI, Nvidia, or Intel GPU, for example. You don't need machine-specific instructions for each one.


I think you may be focusing on the wrong thing here. Moving the data from an SSBO to a uniform buffer might give some speedup, sure, but that is a micro-optimization.
I would search for algorithmic optimizations first. As noted in the comments, processing all 200 records for every vertex is both a lot of data to read, and a lot of computation to do if there's a lot of vertices. Also note that in the course of rendering, vertices will sometimes be processed more than once as they are used by multiple triangles (the vertex cache helps with this, but it can't eliminate all duplicate processing). So you will be re-doing this expensive computation multiple times for the same vertex, sometimes.
Does every vertex really need to process every data element? I see that your data elements have a position and radius. If they represent some kind of mesh deformation that happens within a localized volume of space, then it might be beneficial to store those in an acceleration structure. This would allow each vertex to efficiently narrow down its processing to only those elements that affect it. For example, you could create a regular grid spanning the bounding box of your mesh, and in each grid cell store a list of indices into the main data array for the elements that intersect that grid cell. The creation of this grid could be done in a compute shader earlier in the frame, and should be quite fast. Then, when drawing the mesh, the vertex shader would look up which grid cell the vertex is in, and process those data elements. This is very similar to how we do tiled/clustered shading for light sources in many game engines.
If spatial subdivision like that doesn't get the job done (for example if many of the elements are of large radius, so that most/all of the vertices really do need to process most/all of the elements), then are the elements' influences on the vertices parallelizable in any way? For instance, are you summing a perturbation generated by each element, to arrive at a total displacement for the vertex, or something along those lines? If so, you could run compute shader passes to pre-calculate the displacements and store them into a buffer to be read later by the vertex shader. The compute passes could be set up with work groups that distribute jointly over the vertices and the work elements (the exact distribution and mapping of compute threads would have to be experimented with). Each work group might, for instance, load a small chunk of vertices and a chunk of displacement elements into shared memory, calculate the displacements, then store them out to the per-vertex displacement buffer using atomic adds. Or, they might store out to various intermediate buffers (to avoid atomic contention) and have another set of compute shaders that sum up the intermediate values to obtain the final values. It would take some experimentation to see how to do it most efficiently, but you could likely get some big speedups here.


How fast should I expect to render the silhouette of a 20k triangle model?

OpenGL (scanline rendering) is always going to be way faster than ray tracing.

Ray tracing computes pixels one by one using physics (at least a good 90% of the time it's going to be real physics). Although some of it can now be done with a GPU, it is not really helping as much as one would expect. Calculating speculative lightning is not something you can easily do in parallel (ray 1 will bounce from object A to object B, ray 2 will bounce from object A to object C... you can't parallelize the computations on B and C.) The GPU helps with the matrix math, though.

In comparison, OpenGL is a total fake in comparison. It renders triangles really fast with shades coming from textures and fast gradient computations from your various light sources.


The most data I've seen is on the relative expense of various state changes is from Cass Everitt and John McDonald's talk on reducing OpenGL API overhead from January 2014.  Their talk included this slide (at 31:55):

The talk doesn't give any more info on how they measured this (or even whether they're measuring CPU or GPU cost, or both!).  But at least it dovetails with the conventional wisdom: render target and shader program changes are the most expensive, uniform updates the least, with vertex buffers and texture changes somewhere in the middle.  The rest of their talk also has a lot of interesting wisdom about reducing state-change overhead.


How to convert a thread ID into Screen Space Coord in an OpenGL Compute Shader?

The problem is actually in your fragment shader:

color_out = texture(source_image, gl_FragCoord.xy);

The texture() function accepts normalized coordinates which range from 0.0 to 1.0. The gl_FragCoord built in contains window coordinates which range from (0.0, 0.0) to (window width, window height).

To fix this, change the fragment shader to this:

color_out = texture(source_image, gl_FragCoord.xy / vec2(1280.0, 720.0));

That is assuming your window is 1280x720, change as necessary.


reducing image size increases file size?

JPEG is a lossy format and depends on both the 2d frequency components of the image and the user specified quality level.

It is possible that down-scaling an image can increase the higher frequency components and potentially result in an increase in file size, typically higher frequency components a encoded using a higher number of bits compared to lower frequency components. An example would be to imagine shrinking a checkerboard texture, the smaller it gets the more it shifts into higher frequency.

It may also be possible that your original image was encoded with a low quality factor and your new image has a high quality factor so it can also gain extra bytes from that. On a more technical level this may be because of a mismatch in quantization bands leading to higher number of DCT coefficients being needed compared to the original JPEG. The quality parameter is not standardized and different applications have different ideas of what '90% quality' means.

Another problem is that because you are using JPEG as a source image it will introduce ringing artifacts into the image, which then need more bits to encode when trying to encode it.

Unfortunately there is no easy fix for this apart from reducing quality.


Remove clearly not visible vertices from a polygon

I solved it myself. Here is a result:

The example is out of this post: here


How many 3D textures does OpenGL support

As gllampert pointed out in the comments the value is hardware dependent. You can retrieve it with glGet, using GL_MAX_COMBINED_TEXTURE_IMAGE_UNIT. You can find how different hardware performs here.

However, in OpenGL 3 there is a lower bound of at least 48 simultaneously used textures, no matter which type. [source]


Outcodes in Cohen-Sutherland Clipping Algorithm

The outcodes just represent four boolean flags as a bitfield. Codewise it is easier to move a bitfield around than four booleans. Functionality-wise it becomes very easy to check if any bits are set (just check if zero), and combining outcodes is as simple as a bitwise OR.

On a side-note: Since bit-operations are fast, it is not a performance issue that the individual bits need to be extracted (e.g. outcode & BOTTOM).


When spawning the refraction/transmission ray, don't adjust the ray origin by adding the normal: adjust it by subtracting the normal. This moves the ray origin inside the object, so it won't immediately self-intersect. Your current problem is caused by all your refraction rays immediately hitting the surface.

Don't forget that for accurate refractions, the ray needs to be able to intersect with the inside or back face of the object as well. If you're intersecting back-facing surfaces, you need to flip the normal before shading or adding/subtracting it to the origin.


Will sRGB color space get replaced by wider gamut

Displays are already moving towards wider gamuts. I believe that newer iPhones, for example, are all P3 displays. So there is no doubt that sRGB usage will be replaced by wider gamuts in the future.

There are really 2 important questions here:

Which gamuts will people end up using?
When can/should we move our apps to using wider gamuts?

The answer to the first question will depend on the market. People have been pretty happy with sRGB for a few decades. P3 is better, but most consumers probably can't tell the difference. There are much wider gamuts like ACEScg, Rec. 2020, etc. Will it be affordable to move to those? If it is, will consumers care?

The answer to the second question depends entirely on what your target audience is using. If you can guarantee they're using devices with P3 displays (say, if you require the iPhone X's IR sensor), then you should move to P3. If your users are high-end users creating output for professional use, it's probably wise to move to P3. If your users need to create content that looks reasonably good on end-user devices like web pages rendered on aging machines, then it is probably a waste of effort to move away from sRGB.

Note that there is a cost to moving to a wider gamut. You will need to do significantly more work to ensure that all of your processing is done in the right color space because input images to your app will likely come in a variety of gamuts and you now have to convert them all to your working space, process them, and  then convert to the output space. That's not cheap to do and not easy to get right.


Offset.xy is definatelly off because it just return XY coordinates from screenspace:

Clamping solved 'small copies of model' on the model itself, thanks!

sampleDepth actually looks like this and indeed is not in viewspace and it may really explain those little models now:

Using viewspace (offset = mul(projection, offset) ---> offset = mul(view, offset)) solves depth problem but final buffer is clear white. No effects are actually visible. What next step do you suggest? If depth comparing is correct now, what might be off?

UPDATE:
Comparision is working but in a wrong way. Below line results in just fully white buffer screen:

occlusion += (sampleDepth >= sample.z + bias ? 1.0 : 0.0);

Changing it to:

occlusion += (sampleDepth >= sample.z - bias ? 1.0 : 0.0);

Actually gives us comparision but it is really bold and doesn't include only holes and creases but whole model. Changing bias to small values like 0.001f doesn't help:

Update #2:
Directly comparing gives result presented below. However it doesn't seem satisfying enough to use as it seems to present rather edges than creases and holes:

occlusion += (sampleDepth == sample.z ? 1.0 : 0.0);


Taken from Wikipedia's article on Cohen-Sutherland:

  Both endpoints are in different regions: in case of this nontrivial situation the algorithm finds one of the two points that is outside the viewport region (there will be at least one point outside). The intersection of the outpoint and extended viewport border is then calculated (i.e. with the parametric equation for the line), and this new point replaces the outpoint. The algorithm repeats until a trivial accept or reject occurs.

So yes, since your new first point is inside the viewport, but your second point is outside, the situation is non-trivial (i.e. no "go render it all" or "go discard it all"), it will repeat with the second point. Afterwards, both points will be inside the viewport, therefore it is a trivial draw and no further action is needed.


Is there ongoing research on Markov Chain Monte Carlo algorithms (other than Metropolis-Hastings) for use in ray tracing?

Overview

Here is a short overview of the most used space representations, MLT variants and mutation strategies for these MLT variants. As you can see, there are quite some papers dating back to 2017 (e.g., three papers explore combining the Path Space and the Primary Sample Space by jumping back and forth between the two).

Path Space (PS) representation, Metropolis Light Transport (MLT)

Veach E., Guibas L. J.: Metropolis Light Transport. In Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Techniques (New York, NY, USA, 1997), SIGGRAPH ’97, ACM Press/Addison-Wesley Publishing Co.
Veach E.: Robust monte carlo methods for light transport simulation. Stanford University PhD thesis, 1997.
Van de Woestijne J., Frederickx R., Billen N., Dutré P.: Temporal Coherence for Metropolis Light Transport. Eurographics Symposium on Rendering - Experimental Ideas & Implementations 2017.

[Example implementation available in Mitsuba]

Primary Sample Space (PSS) representation, PSSMLT

Kelemen C., Szirmay-Kalos L., Antal G., Csonka F.: A Simple and Robust Mutation Strategy for the Metropolis Light Transport Algorithm.

Extended Primary Sample Space (EPSS) representation, Multiplexed MLT (MMLT)

Hachisuka T., Kaplanyan A. S., Dachsbacher C.: Multiplexed Metropolis Light Transport. ACM Trans. Graph. 33, 4 (July 2014).

Note that MMLT improves upon PSSMLT.

[Example implementation available in pbrt-v3]

Manifold Exploration MLT (MEMLT)

Jakob W., Marschner S.: Manifold Exploration: A Markov Chain Monte Carlo Technique for Rendering Scenes with Difficult Specular Transport. ACM Trans. Graph. 31, 4 (July 2012).
Jakob W.: Light Transport On Path-Space Manifolds. Cornell University PhD thesis, 2013.

[Example implementation available in Mitsuba]

Half Vector Space representation, Half Vector Space Light Transport (HSLT)

Kaplanyan A. S., Hanika J., Dachsbacher C.: The Natural-constraint Representation of the Path Space for Efficient Light Transport Simulation. ACM Trans. Graph. 33, 4 (July 2014).
Hanika J., Kaplanyan A., Dachsbacher C.: Improved half vector space light transport. In Computer Graphics Forum (2015), vol. 34, Wiley Online Library.

Note that HSLT is more general than MEMLT: MEMLT is an optional add-on mutation for MLT, whereas HSLT replaces the complete Path Space representation.

Gradient Domain representation, GDMLT

Lehtinen J., Karras T., Laine S., Aittala M., Durand F., Aila T.: Gradient-domain Metropolis Light Transport. ACM Trans. Graph. 32, 4 (July 2013).
Manzi M., Rousselle F., Kettunen M., Lehtinen J., Zwicker M.: Improved Sampling for Gradient-domain Metropolis Light Transport. ACM Trans. Graph. 33, 6 (Nov. 2014).

Note that there are variants for path tracing and bidirectional path tracing exploiting the Gradient Domain representation as well.

Path Space <-> Primary Sample Space representations

Reversible Jump MLT (RJMLT)

Bitterli B.: Informed Choices In Primary Sample Space. Master’s thesis, Swiss Federal Institute of Technology Zurich, 2015.
Bitterli B., Jakob W., Novák J., Jarosz W.: Reversible Jump Metropolis Light Transport Using Inverse Mappings. ACM Transactions on Graphics 37, 1 (Oct. 2017).

Charted MLT (CMLT)

Pantaleoni J.: Charted Metropolis Light Transport. ACM Trans. Graph. 36, 4 (July 2017).

Fusing State Spaces

Otsu H., Kaplanyan A. S., Hanika J., Dachsbacher C., Hachisuka T.: Fusing State Spaces for Markov Chain Monte Carlo Rendering. ACM Transactions on Graphics (Proc. of SIGGRAPH) 36, 4 (2017).

Genetic operators in Path Space and Extended Primary Sample Space, Genetic MLT (GMLT)

Delabie M.: Genetic Operators in Metropolis Light Transport. Master’s thesis, KU Leuven, Belgium, 2018.

Note that the Path Space is represented via 3D path vertex positions instead of directions to facilitate crossovers (compared to MLT).

Post Scriptum

  Its variants differ only in the used "mutation strategy".

Different mutation strategies and space representations are used. Furthermore, GMLT replaces and extends Markov states representing single paths to complete path populations (Evolutionary Monte Carlo, Evolutionary Markov Chain Monte Carlo <> Markov Chain Monte Carlo).

Some references to Evolutionary Monte Carlo:

Liang F., Wong W. H.: Evolutionary Monte Carlo: Applications to Cp Model Sampling and Change Point Problem. STATISTICA SINICA 10 (2000).
Liang F., Wong W. H.: Real-Parameter Evolutionary Monte Carlo With Applications to Bayesian Mixture Models. Journal of the American Statistical Association 96, 454 (2001).


Geometric interpretation of this bilinear interpolation equation?

If we look at the two terms with only x or only y they fully describe what happens when you go right or up from the starting point.

The xy-terms does not come into play at all.

When you want to go "into the square" or along the upper or the right edge then the xy-term comes into play. It cancels out the x/y-terms and introduce the P11-value that it should obtain in upper-right.

No full explanation here, just nudging along towards insight!


Which technique is fast to generate real-time reflections?

Usually modern game engines use Screen Space Reflections which is some kind of screen space raymarching.

Basically, a simple way to achieve this is to shoot rays from the G-buffer. This means you use your depth buffer as a height field by raymarching any ray you cast. Then sample the normal at this position and take the final composed lighting from you previous frame and project it.

Although it has the usual drawbacks associated with any screen space method. This means you only have what is visible in front of you (screen space) to compute your reflections, so anything outside the screen would require tricks to have them get reflected on your current frame (e.g. you can combine this with environment mapping reflections and such).

You can find a few articles about that and presentations on the Internet.


I think I got it!

Because $cos(\theta)$ integrates to $\pi$ over the hemisphere (and not $2\pi$). And the incoming light is multiplied by $cos(\theta)$ (and the BRDF).


Using instanced arrays for dynamic objects

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


Jittering and dithering are both techniques of adding noise to reduce visible artefacts (such as banding) in an image. They solve different kinds of artefacts so they are used in different situations. Jittering moves sample positions in space to reduce artefacts caused by regular sampling. Dithering changes the way colours are rounded (when reducing precision), replacing banding with noise.

When is jittering used?

A use case for jittering is when generating eye rays in a ray-tracer. Without jittering, you might generate one ray in the centre of each pixel. But then if you trace those at a regular pattern, such as a grille of bars each one pixel wide, you might find that every ray hits a bar, making the grille solid; or that every ray hits a gap, and the grille isn't seen at all.

Instead, you can add a small stochastic displacement to each ray direction, so that each ray goes through a different part of its pixel. That way, some rays will hit the bars of the grille and some will hit the gaps, making a noisy image instead of a completely wrong image.

When is dithering used?

The normal use of dithering is when rounding colours (intensity values) from a high-precision format (such as floating-point numbers) to a low-precision format (such as 8-bit integers). Say you have a smooth gradient you want to represent in an 8-bit image. If you round each number separately, then starting from the black end of the gradient, you get a lot of black pixels until they get light enough to round to 1/256; then these all have the same value until you get to 2/256; and so on. You end up with regular bands where all the pixels have the same value.

Instead, when you round the first number down, you can let it "donate" the lost intensity to another nearby pixel. That way, then at the black end of the gradient, maybe every tenth pixel is 1/256 and the other nine are 0. The ratio of black to nearly-black pixels increases gradually, until eventually all the pixels are 1/256, and then every tenth is 2/256, and so on. Instead of visible bands, the overall intensity increases smoothly along the gradient - but the price of this is that pixels which should be the same are slightly different, adding noise.

Dithering is used a lot in printing, where the colour resolution is quite poor (i.e. the steps between colours are big) but the spatial resolution (how small you can make each dot of colour) is very fine. It used to be used in the days of 256-colour displays for making colours which couldn't otherwise be represented, and it's still used to avoid banding when displaying (for example) a 16-bit-per-channel image on an 8-bit display.

Is dithering the only way?

Because traditional dithering algorithms involve giving information from each pixel to its neighbours, they don't parallelise well, so they're not suited to modern hardware. Another popular technique is to deliberately add noise to the source image. The range of the noise is less than one step of the destination format (e.g. the noise ranges from -1/512 to 1/512), so it can't be seen directly, but it means that each pixel will be rounded differently. Like true dithering, it replaces any banding with subtle noise in the output image.

Can we combine them?

There's no reason you can't use both in the same image, but as they solve different problems, that's not really "combining" them. It's like asking if we can combine hoses and axes. A fire truck has both, and the firefighters might use both at the same incident, but one is for extinguishing fires and the other is for freeing trapped people.

An example of when you might use both is if you ray-trace an image into a floating-point framebuffer, using jittered ray positions to avoid sampling artefacts, and then reduce that to an 8-bit image using dithering to better preserve the colours.


Controling order of execution of a (compute) shader on an array of vertices (data) in Vulkan

Vulkan is not really suited to at least some of the task you're trying to accomplish. It has no mechanism to declare from within a shader that some additional work should be invoked. Or at least, not directly.

The most effective way to execute this on the GPU would be to compute each level of all of the trees as a distinct dispatch operation. So you process all root nodes first, then the first level for each tree, then the second level, etc. Each compute shader invocation will need to do the processing for its node and store that computed data.

This also requires that you store your tree data in a fairly scrambled way. All of your trees need to live (on the GPU side) in the same general array, preferably with all nodes for a particular level stored contiguously with one another. This makes processing easier, but it makes adding new nodes harder.

The biggest performance problem will likely be the dependency issue. Each layer depends on the execution of the previous layer, which means you won't get much parallelism between layers. You're effectively inducing a stall between each set of layers. If you're executing thousands of nodes, you'll saturate your GPU's execution units just fine; it's the first few layers when you have few nodes that will have less performance than they could.

Now, you may be able to avoid this performance problem by using a different shader for the first couple of layers. Such a shader would have each invocation process, say, the first 3 levels of the tree. On the one hand, this mitigates synchronization costs. But on the other hand, each shader invocation has to read and process a lot more data. And reading/writing a bunch of data from multiple invocations can be expensive, compared to having a more pipelined operation.

So you should profile it and see what works best.

Oh, and FYI: OpenGL compute shaders are just as capable of this as Vulkans.


Phong model: why no multiplication by N dot L in specular term?

The $N \cdot L$ term (I've seen "geometry term" thrown around; not quite accurate, but good enough for here) is part of the rendering equation.  What you're looking at here is an estimate of the inside of the rendering equation—that is, the incoming radiance i times the geometry term dot(n,l) times the BRDF.
The simple answer for why the Phong BRDF doesn't include the geometry term is "Well, it probably should, but the Phong model predates the invention of the rendering equation by more than a decade.".

Phong, B.T.: Illumination for computer generated pictures.  Communications of the ACM 18(6), 311–317 (1975).
Kajiya, J.T.: The rendering equation.  In: Proceedings of the 13th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH ’86, pp. 143–150.  ACM, New York, NY, USA (1986).

This deficiency was fixed in what is now called the Modified Phong BRDF:

Lewis, R.R.: Making Shaders More Physically Plausible.  Computer Graphics Forum 13(2), 109–120 (1994)

At the end of the day, whether the term is present or not isn't an issue of correctness per-se.  Not having the geometry term when you might expect it just means that the BRDF corresponds to a different surface reflection that doesn't foreshorten much at glancing angles.  Indeed, many other BRDFs have terms that look like reciprocal geometry terms, so that they cancel with the actual geometry term in the rendering equation.
In code, the BRDFs look like:
vec3 f_phong    = rho_s/(TWO_PI/(a+1.0)) * pow(max(0.0,dot(r,v)),a) / dot(n,l);
vec3 f_modphong = rho_s/(TWO_PI/(a+2.0)) * pow(max(0.0,dot(r,v)),a);

And so when multiplied by the light and geometry term in the rendering equation (i.e., i * f * dot(n,l)), you get:
vec3 spec_phong    = i * rho_s/(TWO_PI/(a+1.0)) * pow(max(0.0,dot(r,v)),a);
vec3 spec_modphong = i * rho_s/(TWO_PI/(a+2.0)) * pow(max(0.0,dot(r,v)),a) * dot(n,l);

The obvious should be noted, however: while the Phong BRDF is simple to understand and implement, it is of limited accuracy on real-world materials.  It's not terrible, the way many people think it is, but it does lack important features that can be important on many materials.  It is not the state-of-the-art, and depending on your needs you may like to look elsewhere (GGX microfacet, etc.) for more modern BRDFs.


Incorrect values written in stencil buffer

OK, so this not really an answer, rather just sharing what I have learned.

To clarify, shape of the object has nothing to do with this problem at all. Be it a triangle or rectangle, and wherever you click on that object, this method should work (the 'But' for this part is explained below).

I think stencil buffer values are associated with each pixel on screen, and it does not directly maps to the OpenGL window space. Initially I was using a window size that is smaller than my screen size and on top of that it was not full screen. Each time I started (debugging from VS 2015), the window appeared not exactly at (0, 0) on the screen. Some of the screen space is also consumed by the window bar at the top. I think this is the reason for getting erratic values when reading back the data from stencil buffer.

So as a simple solution to this problem what I did is to make the window full screen, and everything worked like charm. As I use GLFW, here is how I have done it:

glfwCreateWindow(Width, Height, "Title", glfwGetPrimaryMonitor(), nullptr);

Just used glfwGetPrimaryMonitor(), while creating the window. And that's it.

And yes, it works for triangles too.


Issue with declaring inputs to geometry shader?

The issue was that I declared out vec3 gModelPosition[] as an array in my geometry shader.

It should have only been declared as a vec3 (i.e. out vec3 gModelPosition)


Are non manifold meshes problematic for Virtual Reality?

No, VR is no different to flat graphics applications in this respect. If your VR experience has physics, CSG (Boolean geometry), or any of the things that would have problems with non-manifolds in a flat application, then it'll have the same problems. If your application wouldn't have problems with non-manifolds as a flat application, it won't have problems as a VR application.

The only real difference is that it's harder to stop the camera being inside objects in VR, since it's attached to the user's head. In a flat application, you might override the camera position or adjust the camera controls to avoid the camera moving through objects in the first place. In VR, those techniques don't really apply, so you might want to have an inside-outside test, so that you can react if the user puts their head inside an object (e.g. by blanking the view completely and/or asking them to move). The inside-outside test would obviously have problems with non-manifolds. It's the same problem as in a flat application, but in a flat application you might not be doing the test at all.


Why is the back of a perspective frustrum larger than the front?

It might help to think of it this way:

Both the near plane and the far plane are the size that will fit onto the screen you are viewing on. The further away something is, the bigger it can be and still fit on the screen. Close up, a playing card can fill the whole screen, but in the distance, an entire building can fit on the screen.

Consider the example image:

If a playing card exactly fills the near plane, then you will be able to fit 4 playing cards on the far plane. This reflects the fact that when a playing card is twice as far from the camera, it appears to have half the height and width.

Think of the near plane and the far plane as simply the area that will fit on your screen. The further away you point your camera, the more will fit on the screen. This matches up with the fact that things look smaller in the distance - this means you can fit bigger things on the screen if they are further away.


What are forms of affine transformations?

Note: I have answered before the edit from trichoplax and I thought you were searching for other transformations other than the one you mentioned. The informations below are still useful so I will keep the answer here, but it does not directly answer your question.

Affine transformations (surprise!) map affine spaces to affine spaces. An affine space is substantially a vector space where you can establish an origin and define points as tuples of their coordinates. This is far from a formal definition, if you are interested in one I can edit the answer.

Now an affine transformation $T$ transform points into points and must preserve affine combinations:

$T(\lambda_1 P_1 + \lambda_2 P_2 + ...  +  \lambda_n P_n) = \lambda_1 T(P_1) + \lambda_2 T(P_2) + ... + \lambda_n T(P_n)  $

Given $\sum_{i}^{n} \lambda_i = 1 $

However these transformations can't be arbitrary as the following must be preserved:

Parallelism is preserved. This means that if you transform parallel lines they remain parallel after the affine transformation.
Relative ratios are preserved. This means that if you have $R = (1-\beta) P + \beta Q$  then $T(R) = (1 - \beta) T(P) + \beta T(Q)$
Addition between vector and points are preserved. Meaning that $T( P + \vec{v} ) = T(P) + T(\vec{v})$

With these properties in mind you can come up with a very big number of affine transforms yourself. A couple of obvious ones other the ones you mentioned are:

Translation

$$T_vp =
\begin{bmatrix}
1&0&0&v_x\\0&1&0&v_y\\0&0&1&v_z\\0&0&0&1\end{bmatrix}
\begin{bmatrix}
p_x+v_x\\p_x+v_y\\p_z+v_z\\1\end{bmatrix}=p+v$$

That moves a point into a specific direction by a specific amount specified by a displacement vector.

Shearing

image from wikipedia

That is a transform that displace all points in a given direction by an amount that is proportional to their perpendicular distance to a line parallel to that direction.

For example the transform matrix for an horizontal shear in a 2D space is given by:

$$
\begin{bmatrix}x'\\y'\end{bmatrix} = \begin{bmatrix}x+my\\y\end{bmatrix} =
\begin{bmatrix}1&m\\0&1\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}$$

Again, these two are just example, the important information you should really keep is the definition (and properties) of an affine transform; with that in mind it shouldn't be too hard to recognize an affine transform.  Also note that combining affine transforms will give you an affine transform!


There is no relation between direction's direction and the world axis. And that's fortunate otherwise it would mean your camera is not a free view, it's some kind of a axis bound camera, which has its usages but most likely nothing in your mind right now.

The default camera matrix is looking at +z when everything is Identity, that's most surely where you got confused. Know that this is purely a convention, but sticking to a widely adopted convention allows for easy compatibility with libraries like glm.

the direction vector goes outward, so its normalize(lookat - campos).

The view matrix can be constructed directly with TBN vectors in rows 0 1 2 :)
(that's your RUF)

Point 2. There is no translating of the camera to the origin, that's just how you chose to view (pun) it. There is translation of world objects toward the origin. Because the rasterizer will work in device coordinates (NDC) which is not configurable, so having a roaming camera (that travels) makes it necessary to indicate its translation in world as part of the last row (row 3), and in reverse since it's not a world matrix representing the camera position, it's the matrix that will bring back points into the view space. (by the way, incidentally meaning the actual camera world matrix is the inverse of the view matrix).

oh, that was the last phrase of pt 2.


Short answer:

Importance sampling is a method to reduce variance in Monte Carlo Integration by choosing an estimator close to the shape of the actual function.

PDF is an abbreviation for Probability Density Function. A $pdf(x)$ gives the probability of a random sample generated being $x$.

Long Answer:

To start, let's review what Monte Carlo Integration is, and what it looks like mathematically.

Monte Carlo Integration is an technique to estimate the value of an integral. It's typically used when there isn't a closed form solution to the integral. It looks like this:

$$\int f(x) \, dx \approx \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{pdf(x_{i})}$$

In english, this says that you can approximate an integral by averaging successive random samples from the function. As $N$ gets large, the approximation gets closer and closer to the solution. $pdf(x_{i})$ represents the probability density function of each random sample.

Let's do an example: Calculate the value of the integral $I$.

$$I = \int_{0}^{2\pi} e^{-x} \sin(x) dx$$

Let's use Monte Carlo Integration:

$$I \approx \frac{1}{N} \sum_{i=1}^N \frac{e^{-x} \sin(x_i)}{pdf(x_{i})}$$

A simple python program to calculate this is:

import random
import math

N = 200000
TwoPi = 2.0 * math.pi

sum = 0.0

for i in range(N):
    x = random.uniform(0, TwoPi)

    fx = math.exp(-x) * math.sin(x)
    pdf = 1 / (TwoPi - 0.0)

    sum += fx / pdf

I = (1 / N) * sum
print(I)

If we run the program we get $I = 0.4986941$

Using separation by parts, we can get the exact solution:

$$I = \frac{1}{2} (1 − e−2π) = 0.4990663$$

You'll notice that the Monte Carlo Solution is not quite correct. This is because it is an estimate. That said, as $N$ goes to infinity, the estimate should get closer and closer to the correct answer. Already at $N = 2000$ some runs are almost identical to the correct answer.

A note about the PDF: In this simple example, we always take a uniform random sample. A uniform random sample means every sample has the exact same probability of being chosen. We sample in the range $[0, 2\pi]$ so, $pdf(x) = 1 / (2\pi - 0)$

Importance sampling works by not uniformly sampling. Instead we try to choose more samples that contribute a lot to the result (important), and less samples that only contribute a little to the result (less important). Hence the name, importance sampling.

If you choose a sampling function whose pdf very closely matches the shape of $f$, you can greatly reduce the variance, which means you can take less samples. However, if you choose a sampling function whose value is very different from $f$, you can increase the variance. See the picture below:

Image from Wojciech Jarosz's Dissertation Appendix A

One example of importance sampling in Path Tracing is how to choose the direction of a ray after it hits a surface. If the surface is not perfectly specular (ie. a mirror or glass), the outgoing ray can be anywhere in the hemisphere.

We could uniformly sample the hemisphere to generate the new ray. However, we can exploit the fact that the rendering equation has a cosine factor in it:

$$L_{\text{o}}(p, \omega_{\text{o}}) =  L_{e}(p, \omega_{\text{o}}) + \int_{\Omega} f(p, \omega_{\text{i}}, \omega_{\text{o}}) L_{\text{i}}(p, \omega_{\text{i}}) \left | \cos \theta_{\text{i}} \right | d\omega_{\text{i}}$$

Specifically, we know that any rays at the horizon will be heavily attenuated (specifically, $\cos(x)$ ). So, rays generated near the horizon will not contribute very much to the final value.

To combat this, we use importance sampling. If we generate rays according to a cosine weighted hemisphere, we ensure that more rays are generated well above the horizon, and less near the horizon. This will lower variance and reduce noise.

In your case, you specified that you will be using a Cook-Torrance, microfacet-based BRDF. The common form being:

$$f(p, \omega_{\text{i}}, \omega_{\text{o}}) = \frac{F(\omega_{\text{i}}, h) G(\omega_{\text{i}}, \omega_{\text{o}}, h) D(h)}{4 \cos(\theta_{i}) \cos(\theta_{o})}$$

where

$$F(\omega_{\text{i}}, h) = \text{Fresnel function} \\
G(\omega_{\text{i}}, \omega_{\text{o}}, h) = \text{Geometry Masking and Shadowing function} \\
D(h) = \text{Normal Distribution Function}$$

The blog "A Graphic's Guy's Note" has an excellent write up on how to sample Cook-Torrance BRDFs. I will refer you to his blog post. That said, I will try to create a brief overview below:

The NDF is generally the dominant portion of the Cook-Torrance BRDF, so if we are going to importance sample, the we should sample based on the NDF.

Cook-Torrance doesn't specify a specific NDF to use; we are free to choose whichever one suits our fancy. That said, there are a few popular NDFs:

GGX
Beckmann
Blinn

Each NDF has it's own formula, thus each must be sampled differently. I am only going to show the final sampling function for each. If you would like to see how the formula is derived, see the blog post.

GGX is defined as:

$$D_{GGX}(m) = \frac{\alpha^2}{\pi((\alpha^2-1) \cos^2(\theta) + 1)^2}$$

To sample the spherical coordinates angle $\theta$, we can use the formula:

$$\theta = \arccos \left( \sqrt{\frac{\alpha^2}{\xi_1 (\alpha^2 - 1) + 1}} \right)$$

where $\xi$ is a uniform random variable.

We assume that the NDF is isotropic, so we can sample $\phi$ uniformly:

$$\phi = \xi_{2}$$

Beckmann is defined as:

$$D_{Beckmann}(m) = \frac{1}{\pi \alpha^2\cos^4(\theta)} e^{-\frac{\tan^2(\theta)}{\alpha^2}}$$

Which can be sampled with:

$$\theta = \arccos \left(\sqrt{\frac{1}{1 = \alpha^2 \ln(1 - \xi_1)}} \right) \\
\phi = \xi_2$$

Lastly, Blinn is defined as:

$$D_{Blinn}(m) = \frac{\alpha + 2}{2 \pi} (\cos(\theta))^{\alpha}$$

Which can be sampled with:

$$\theta = \arccos \left(\frac{1}{\xi_{1}^{\alpha + 1}} \right) \\
\phi = \xi_2$$

Putting it in Practice

Let's look at a basic backwards path tracer:

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

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

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

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

        // We hit an object

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

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

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

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

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

        // Shoot a new ray

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

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

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

IE. we bounce around the scene, accumulating color and light attenuation as we go. At each bounce, we have to choose a new direction for the ray. As mentioned above, we could uniformly sample the hemisphere to generate the new ray. However, the code is smarter; it importance samples the new direction based on the BRDF. (Note: This is the input direction, because we are a backwards path tracer)

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

Which could be implemented as:

void LambertBRDF::Sample(float3 outputDirection, float3 normal, UniformSampler *sampler) {
    float rand = sampler->NextFloat();
    float r = std::sqrtf(rand);
    float theta = sampler->NextFloat() * 2.0f * M_PI;

    float x = r * std::cosf(theta);
    float y = r * std::sinf(theta);

    // Project z up to the unit hemisphere
    float z = std::sqrtf(1.0f - x * x - y * y);

    return normalize(TransformToWorld(x, y, z, normal));
}

float3a TransformToWorld(float x, float y, float z, float3a &normal) {
    // Find an axis that is not parallel to normal
    float3a majorAxis;
    if (abs(normal.x) < 0.57735026919f /* 1 / sqrt(3) */) {
        majorAxis = float3a(1, 0, 0);
    } else if (abs(normal.y) < 0.57735026919f /* 1 / sqrt(3) */) {
        majorAxis = float3a(0, 1, 0);
    } else {
        majorAxis = float3a(0, 0, 1);
    }

    // Use majorAxis to create a coordinate system relative to world space
    float3a u = normalize(cross(normal, majorAxis));
    float3a v = cross(normal, u);
    float3a w = normal;

    // Transform from local coordinates to world coordinates
    return u * x +
           v * y +
           w * z;
}

float LambertBRDF::Pdf(float3 inputDirection, float3 normal) {
    return dot(inputDirection, normal) * M_1_PI;
}

After we sample the inputDirection ('wi' in the code), we use that to calculate the value of the BRDF. And then we divide by the pdf as per the Monte Carlo formula:

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

Where Eval() is just the BRDF function itself (Lambert, Blinn-Phong, Cook-Torrance, etc.):

float3 LambertBRDF::Eval(float3 inputDirection, float3 outputDirection, float3 normal) const override {
    return m_albedo * M_1_PI * dot(inputDirection, normal);
}


Map a texture onto a hyperbolic triangle

My guess is that to get the texture to look right, you'll have to subdivide the interior of the triangle as well, and approximate the non-linear UV mapping within it.

Currently, it looks like you're subdividing around the edges of the triangle, and forming a fan of smaller triangles between the edge and the incenter. This would be fine if you were just rendering the triangle in a solid color. But when applying a texture, you're getting discontinuities at the boundaries of those subdivisions, because the texture is just stretched linearly across each subdivided triangle and not correctly approximating the hyperbolic mapping along the radial direction.

You'll need to subdivide along both axes, something like this:

and position all the vertices appropriately in screen space and UV space to approximate the hyperbolic coordinates in the interior of the triangle. If you subdivide finely enough, this should produce the illusion of a continuously curving texture mapping.


My question is, how should I (...) rotate normals of my mesh.

Normals are rotated the same way as the rest of your mesh/point cloud. But make sure that you use a "pure rotation matrix" for this (no scaling, no translation). Since rotations preserve lengths, you do not need to renormalize your normals (usually - floating-point errors might accumulate).

This sentence of you

  Note that this rotation also needs to keep object in unit cube, where it was before running this code.

lets me guess, that you are actually asking how to rotate an object and its normals and not just the normals. In python, you can use the packages SciPy and NumPy for that.

The Rotation package (documentation) of SciPy helps you to generate an arbitrary rotation matrix. Then you can use NumPys matmul function to rotate your points and normals.

Here is a short example script that rotates a single point 20 degrees around the z axis:

from scipy.spatial.transform import Rotation as Rot
import numpy as np

p = np.array([1, 1, 0]).transpose()
rot_mat = Rot.from_euler(seq="z", angles=20, degrees=True).as_matrix()

p_r = np.matmul(rot_mat, p)

print(p)
print(p_r)
print(np.linalg.norm(p))
print(np.linalg.norm(p_r))

output:

[1 1 0]
[0.59767248 1.28171276 0.        ]
1.4142135623730951
1.414213562373095

  Info: Depending on your SciPy version you might need to replace as_matrix with as_dcm

As you can see, apart from some minor floating-point error, the length (last 2 outputs) did not change. Furthermore, I chose the example point, so that it gets outside the unit cube. To address this problem, you need to find the maximum absolute value inside the transformed data. You can use amax and amin to get the maximum and minimum value (see also here) if you put all your coordinates into one large 2d array (a matrix). Then compare the absolute values of the amax and amin results and pick the larger one. Divide all your coordinates (not the normals!) by this value, to scale your model back into the unit cube. Alternatively, you can create a uniform scaling matrix

$$
\begin{bmatrix}
s&0&0\\
0&s&0\\
0&0&s
\end{bmatrix}
$$

with $s$ being the reciprocal of the maximum absolute value. Multiply this matrix with all points.

To sum it up:

calculate the rotation matrix with SciPy
Multiply the rotation matrix with all points and normals
find the maximum absolute value in the rotated data
scale your coordinates with the found maximum value

Remark: If you put your points and normals into 2 matrices, you can apply the rotations and scaling matrices without a loop, which makes it much faster.


Why does opengl use 4d matrices for everything?

The bottom row allows you to create perspective foreshortening. That is, it makes lines that are getting further away appear to converge. When arranged this way, we call this a perspective projection matrix.

There are other ways to arrange a projection matrix where that foreshortening doesn't happen. For example in an orthographic projection. This graphic shows some different projections that are possible by  changing the values in the last row of a 4x4 matrix.

Mathematically, the difference between a perspective projection and an orthographic projection is that the last row of a perspective matrix has a division in the z component, whereas the orthographic projection has a 0 in the z component.


Is placing z value of vertex in w enough to achieve perspective projection in OpenGL?

The projection matrix distorts the view frustum (the volume the camera can see) into a unit cube. So everything with all coordinates in the range -1 to 1 after projection is potentially visible, and everything else can be clipped.

To make things further away look smaller, we need to divide by their Z distance. We can't put this value directly into the matrix since it is different for each point we want to project, so instead we place each point's Z into its projected W, giving (X, Y ,Z ,1) -> (X', Y', Z', Z), known as clip coordinates. This is then divided by the new W to give (X'/Z, Y'/Z, Z'/Z, 1),  known as normalized device coordinates. On the GPU this division happens automatically after the vertex shader and before clipping and rasterization.

The general form of a projection matrix is

A 0 0 0
0 B 0 0
0 0 C D
0 0 1 0

The A and B terms relate to the field of view, higher values give a narrower view. Specifically, B is the cotangent of half the vertical view angle, and A is B/aspect ratio.

The C and D terms are bias and scale terms which map the near plane to -1 (0 in DirectX), and the far plane to 1. These are used for writes to the depth buffer to determine which objects occlude others. The depth mapping is actually a function of 1/Z, rather than Z, since this gives better precision close to the near plane, and allows perspective-correct interpolation.

Also note that in OpenGL, W' is -Z rather than Z, since we are using a right-handed coordinate system, with negative Z pointing into the screen.


Resolution issue of generating a PNG icon for macOS menu bar

The issue is the difference between pixel and point sizes. On macOS and other Apple platforms, the user interface is measured in points, which don’t necessarily map one-to-one to pixels; on a Retina screen like the one you’re using, 1 point is 2 pixels. Since there’s no way for the system to know whether your image’s scale is the same as the screen’s, it’s assuming that each pixel in your larger image is meant to be one point.

One way to resolve this would be to explicitly set the size property on your NSImage to the point size your image is supposed to be drawn at, i.e. 19×19; you’d have to handle loading different images for 1x and 2x screens on your own, though, which is inconvenient. A better way would be to use an asset catalog in your project, which will let you specify assets for 1x and 2x screens and allow you to load the appropriate image by name.

As Simon F pointed out, future questions in this vein probably belong on Stack Overflow—it’s more about API usage than graphics concepts.


Rendering equation - why unsolvable directly?

I'm sadly not able to add a comment to the answer above (not enough reputation), so I will do it like this.

I'd like to point out that what Dragonseel describes is simply an integral equation (specifically a Fredholm equation of the second kind). There are many such equations which do have an analytic solution; even some forms of the rendering equation have one (e.g. the solution of a white furnace can be given using a simple convergent geometric series, even though the rendering equation is infinitely recursive).

It is also not necessary to bias the estimated solution by bounding the number of recursions. Russian Roulette provides an useful tool for giving us an unbiased solution for an infinitely recursive rendering equation.

The main difficulty lies in the fact that the functions for reflectance (BRDF), emitted radiance and visibility are highly complex and often contain many discontinuities. In these cases there often is no analytic solution, or it is simply unfeasible to find such a solution. This is also true in the one dimensional case; most integrals lack analytic solutions.

Finally I'd like to note that even though most cases of the rendering equation do not have analytic solutions, there is a lot of research in forms of the rendering equation which do have an analytic solution. Using such solutions (as approximations) when possible can significantly reduce noise and can speed up rendering times.


How does state of the art real time hair rendering work?

Sorry for the bad quality of my answer. I do not have access to a computer currently and editing from my phone is not a straightforward task. In particular I would love to be able to paste images.

I would say that the main challenges of simulating hair are:

replicating their very specific reaction to lighting (as a material)
replicating their volumetric property (as a geometry)
animating them in a realistic manner

Here is a list of references I gathered on the subject, ordered chronologically (mainly about the rendering part):

JAMES T. KAJIYA, TIMOTHY L. KAY, Rendering Fur With Three Dimensional Textures, 1989
JEROME LENGYEL, EMIL PRAUN, ADAM FINKELSTEIN, HUGUES HOPPE, Real-Time Fur over Arbitrary Surfaces, 2000
STEPHEN R. MARSCHNER, HENRIK WANN JENSEN, MIKE CAMMARANO, Light Scattering from Human Hair Fibers, 2003
ARMIN BRUDERLIN, SUNIL HADAP, TAE-YONG KIM, NADIA MAGNENAT-THALMANN , ULRICH NEUMANN, YIZHOU YU, STEVE WORLEY, Photorealistic Hair Modeling, Animation, and Rendering, Siggraph 2003 Course
THORSTEN SCHEUERMANN, Hair Rendering and Shading, ShaderX3, 2004
THORSTEN SCHEUERMANN, Hair Rendering and Shading, ATI / GDC 2004
THORSTEN SCHEUERMANN, Practical Real-Time Hair Rendering and Shading, ATI / Siggraph 2004
THORSTEN SCHEUERMANN, Practical Real-Time Hair Rendering and Shading, ATI / Siggraph 2004 (slides)
MARTIN KOSTER, JORG HABER, HANS-PETER SEIDEL, Real-Time Rendering of Human Hair using Programmable Graphics Hardware, 2004
HUBERT NGUYEN, WILLIAM DONNELLY, Hair Animation and Rendering in the Nalu Demo, GPU Gems 2 (Chapter 23), 2005
LENA PETROVIC, MARK HENNE, JOHN ANDERSON, Volumetric Methods for Simulation and Rendering of Hair, Pixar 2006
CHRISTOPHE HERY, RAVI RAMAMOORTHI, Importance Sampling of Reflections from Hair Fibers, Pixar 2007
SARAH TARIQ, LOUIS BAVOIL, Real-Time Hair Rendering on the GPU, Siggraph 2008.
ARNO ZINKE, CEM YUKSEL, ANDREAS WEBER, JOHN KEYSER, Dual Scattering Approximation for Fast Multiple Scattering in Hair, Siggraph 2008
IMAN SADEGHI, HEATHER PRITCHETT, HENRIK WANN JENSEN, RASMUS TAMSTORF, An Artist Friendly Hair Shading System, 2010
EUGENE D’EON, GUILLAUME FRANCOIS, MARTIN HILL, JOE LETTERI, JEAN-MARIE AUBRY, An Energy-Conserving Hair Reflectance Model, 2011
XUAN YU, JASON C. YANG, JUSTIN HENSLEY, TAKAHIRO HARADA, JINGYI YU, A Framework for Rendering Complex Scattering Effects on Hair, 2012
JIAWEI OU, FENG XIE, PARASHAR KRISHNAMACHARI, FABIO PELLACINI, ISHair: Importance Sampling for Hair Scattering, 2012
SARAH INVERNIZZI, On Physically Based Hair Rendering, 2013
WOLFGANG ENGEL, Hair Rendering in Tomb Raider, 2013
EUGENE D’EON, STEVE MARSCHNER, JOHANNES HANIKA, Importance Sampling for Physically-Based Hair Fiber Models, 2013
TIMOTHY MARTIN, WOLFGANG ENGEL, NICOLAS THIBIEROZ, JASON YANG, AND JASON LACROIX, TressFX: Advanced Real-Time Hair Rendering, GPU Pro 5, 2015
LEONID PEKELIS, CHRISTOPHE HERY, RYUSUKE VILLEMIN, JUNYI LING, A Data-Driven Light Scattering Model for Hair, Pixar 2015

I did not read all of these material yet, but here is a little of what I remember (I am more interested by the real-time solutions) :

Kajiya-kay laid the basis for realistic hair / fur rendering
Marschner proposed an improved model, and noted that hair has two distinct specular highlights (one white due to direct reflection, the other colored and sparkling due to reflection inside of the hair fiber)
Scheuermann proposed approximations to simulate the behaviour described by Marschner, approximations suited for real-time (using tricks such as two separate specular highlights, artist-designed specular offset and exponent, noise textures). He also proposed using three rendering passes (opaque, transparent back and front) of polygons to simulate the "volumetric" qualities of hair
Most of the games I heard of seem to use Scheuermann approaches for hair shading. However, this is not an energy-conserving solution. Nevertheless, I think d'Eon attempts to fix this problem in the 2011 paper.
An other topic of interest is the use of tessellation, that has also been used recently to replace the polygonal approach and simulate hair strands individually.

The offline rendering world attempts to solve more advanced challenges. I think most of them are described in the references above.