Checking the correct offset direction for Bezier curve offset

For a two dimensional spline it is possible to define a left and a right side, or let us call the directions counterclockwise and clockwise side (or maybe port and starboard would be appropriate). The side is easily identified as the cross product of the tangent in direction of point order crossed by the paper plane.

Image 1: Offset side in 2D.

The higher dimensional cases are more problematic, I will be examining the 3D space only. A free space curve has no inherent good second direction. Mathematicians have defined, what is essentially a hack, that the second direction is in the direction of curvature which has a mathematical definition. The only problem is that curvature is subject to change and can flip direction along the curve. This is not good.

So how to counter this? Well simplest is if a curve on a surface of a 2D plane or surface. In this case you can simply use the underlying primitives directions as reference and do a 2D case. This has the added benefit of not doing a somewhat random offset direction.

For the free curve case there are several approaches you can use, they all have the downside of needing to keep some extra data with you.

One is to calculate the location where the curvature is 0 and "cut" the curve at that point and continue with a flipped curvature direction. In case of a straight line just continue with last good direction. A variation of this is not to use curvature but a user defined vector and use that instead.

This can also be done numerically by sampling multiple points along and detect the flip by comparing nearby result and not letting it turn more than a certain degree untill you declare a flip (what your doing sound like a variation of this). One other is to propagate a vector along the curve. One way to do this is to specify a vector at each point and let the user flip them if he feels the need. Or interpolate vectors only at certain points on the curve.

Most applications simply just use the will not do it or bug out option. This has the benefit of being rather sane, as the curvature based bi-normal is a bit erratic to say the least. This way the user has to come up with a way to make the directional control with for example a ruled surface.

*  Bézier, NURBS or something other like a Catmul-Rom spline doesn't matter


It's generally a good idea to add noise like this when you're using a gradient, to avoid visible banding in the gradient, especially on smartphones. Often smartphone screens claim to be 24-bit colour but the panel itself is actually only 16 or 20 bit.

The chief difference between the example and your attempt to reproduce it is that your noise is in RGB space, while the example only adds luminance noise (i.e. brightness). You could achieve this by using HSV colours and only adding to the V channel, or with RGB colours and adding the same amount to each channel.


Rotate a cylinder from XY plane to given points

When you say $acos(\|v\|\|p\|)$ are you saying $acos(\|\hat{v} \cdot \hat{p}\|)$? If not the second, your angle will be some bizarre thing (In fact, you'll always get 90deg).
How are you applying the rotation? glRotate? Have you tried glRotate with the negative of your estimated angle?
Are you doing rotation/scaling/translation in the right order? Could you post the actual code that is executing if it is short? (If not, then a small runnable excerpt?)

Also: I'd be quite careful, as the cross-product-as-axis approach will produce goofy results when $v \parallel u$.


This is not a full answer, I would just like to share the knowledge I obtained by studying two of the papers mentioned in the question: Steerable Importance Sampling and Practical Product Importance Sampling for Direct Illumination.

Steerable Importance Sampling

In this paper they propose a method for sampling the product of the clamped cosine component and environment map lighting:

$$
L_{EM}\left(\omega_{i}\right)\left(\omega_{i}\cdot n\right)^+
$$

They make use of the fact that a piece-wise linear approximation of the product function can be relatively well expressed and partially pre-computed using the first nine spherical harmonic bases. They build this approximation on top of an adaptively triangulated EM and use it as an importance function for sampling.

They pre-compute and store approximation coefficients for each triangle vertex and also coefficients for computation of approximation integral over the triangle for each triangle. These coefficients are called vertex and triangle weights. Then they make use of the fact that is it possible to easily compute coefficients for an integral over a set of triangles just by summing the individual triangle weights without incorporating additional spherical harmonic bases. This allows them to build a balanced binary tree over the triangles where each node contains coefficients for computing approximation integral over the node's sub-tree triangles.

The sampling procedure consists of selecting a triangle and sampling its area:

A triangle is chosen by descending down the pre-built binary tree    with probability proportional to the sub-integral approximations. This costs $O\left(\log N_{\triangle}\right)$ on-the-fly computations of sub-integrals, each consisting of one inner product of clamped-cosine spherical harmonic coordinates with the pre-computed coefficients.
The chosen triangle surface is then sampled in $O\left(1\right)$ time in a    bi-linear fashion by a novel stratified sampling strategy proposed in    the paper.

To me, this looks like a promising technique, but the classical question with papers is how it will behave in the real life. On the one hand, there may be pathological cases when the EM is hard to approximate with triangulated piece-wise linear function, which can lead to an enormous amount of triangles and/or to poor sample quality. On the other hand, it can instantly provide a relatively good approximation of the whole EM contribution, which can be useful when sampling multiple light sources.

Practical Product Importance Sampling for Direct Illumination

In this paper they propose a method for sampling the product of environment map lighting and cosine-weighted surface reflectance:

$$
L_{EM}\left(\omega_{i}\right)f_{r}\left(\omega_{i},\omega_{o},n\right)\left(\omega_{i}\cdot n\right)^+
$$

The only pre-processing in this method is computation of a hierarchical representation of the EM (either mipmap or wavelet based). The rest is done on the fly during the sampling.

The sampling procedure:

Building an on-the-fly BRDF approximation: They first draw several BRDF importance samples and evaluate $f_{r}\left(\omega_{i},\omega_{o},n\right)\left(\omega_{i}\cdot n\right)^+$. From these values they build a quadtree-based piece-wise constant approximation of the BRDF, where each leaf of the tree contains exactly one sample.
Computing a product of the BRDF approximation and the EM: Multiplication is done at the BRDF quadtree leaves and averaged values are propagated to parents.
Product sampling: uniform samples are fed through the product tree using simple sample warping.

The procedure should generate relatively good samples at the cost of heavy pre-computation – they show that roughly 100–200 BRDF samples are needed for BRDF approximation to achieve the best sampling performance. This may make it suitable for purely direct illumination computations, where one generates many samples per shading point, but is most probably too expensive for global illumination algorithms (e.g. uni- or bi-directional path tracers), where you usually generate only a few samples per shading point.


How to produce simple 2D illustrations to accompany geometry answers?

I use a variety of tools depending on the particular needs of the moment. I have a bias toward free tools, as I don't have access to things like Illustrator or Mathematica.

For general vector art purposes I use the free Inkscape. It operates on .svg files and can also render them out as high-quality .pngs or other bitmap formats. For example, here's a diagram I made in Inkscape to explain UV mapping:

For 2D function graphs, I often use the online graphing calculator Desmos, which has the handy capability of letting you share graphs online so others can view and modify them. You can also create sliders that can be varied to play with curve parameters in real time. For example, here is a configurable sharpening function.

For more complex graphing tasks that involve a bit of programming, I use Maxima, a free computer algebra system that also has graphing features. I used Maxima to create the graphs in this answer about Gaussian blurs:

This was a bit too complex for Desmos, since it involved programmatically creating a piecewise approximation of a Gaussian.

For really odd graphing tasks that aren't handled well by anything else, I've used Python with the svgwrite package to programmatically generate diagrams into .svg files. These can then be touched up in Inkscape if desired and rendered to bitmaps using Inkscape or ImageMagick. I used that workflow to create the images in this blog post on depth precision, like this:

In this case I wanted to not only graph a function, but place tick marks programmatically and nonuniformly along both the axes and the graph, to illustrate discrete depth buffer values.

And finally, if nothing else suffices, I'll use Python with numpy and pillow to render images pixel-by-pixel. That's how I created the Perlin noise images in this answer (with the help of the noise package too, in this case).


Does angular diameter to FoV ratio equal object's image to total image ratio?

This is an interesting question, as intuitively we would think that the image size of the object will be 50% of the image width. However, upon some playing around with triangles, I found this to not be the case.
For example, take a look at the diagram below. Here we have two overlapping triangles, one with angular diameter of 50 degrees (representing the object's angular diameter), and one with angular diameter 100 degrees (representing the camera's FoV). As we can see, the object has a width of 1 unit as seen from the camera, however the image width is around 3 units - therefore in this case the horizontal image size of the object will be around 33% of the image width.

If you were to play around with different values for D and d, you would find that the angular diameter ratios would vary. There would be a way to calculate this ratio directly with some trigonometry, however I do not know of how much use that would be to you.
Hope this helps in some way :)


How can deal with batch rendering when the elements in the batch changed every frame?

Allocate a VBO with a reasonable amount of storage during startup, and update the data in it each frame. This will help the driver manage the memory efficiently. Don't destroy and recreate the buffer unless absolutely necessary; too much resource creation/destruction churn will put more pressure on the driver and can potentially lead to stalls.

To be more specific: for data that will be updated every frame, you probably want to initialize it using glBufferStorage with GL_DYNAMIC_STORAGE_BIT (or, if you're on an older OpenGL version, glBufferData with GL_DYNAMIC_DRAW). Then, to update it each frame, use glMapBufferRange with GL_MAP_INVALIDATE_BUFFER_BIT. This should be an efficient, fast path in the graphics driver.


For an introduction, you can give a look at: Ray Tracing from the Ground Up (Amazon link here).

It starts really from the basics, and provides simple implementations for the concepts that are gradually introduced.


Why is the scan line filling algorithm so seemingly over complicated?

The scan-line algorithm (as described on Wikipedia for instance) is concerned with generating the pixels in order, left-to-right and top-to-bottom, with each pixel needing to be touched only once. It was developed in the late 1960s, for devices with no framebuffer memory—so it has to generate each pixel just-in-time as it scans out to the display.

The constraint to generate the pixels exactly in order is a strong one. Now, it's possible to imagine a very simple and naive algorithm that doesn't need tables, buckets, or sorting to accomplish this. For example, for each pixel, you might just iterate over all the triangles, test whether the current pixel falls inside it, and keep track of the closest triangle that passes that test. (That's effectively ray tracing without any acceleration structure.)

It's pretty clear that this is inefficient as soon as you have more than a handful of triangles. The complicated data structures and stuff that show up in the classic scan-line algorithm are there to optimize this process. For instance, pre-bucketing the edges by Y-coordinate and maintaining an "active edge table" lets you quickly and incrementally identify the edges that affect each scanline as you move down the image. Keeping the edge intersections sorted by X-coordinate allows you to quickly generate the pixels left-to-right within each scanline, and keeping a Z-sorted list of active triangles as you scan left-to-right enables you to do hidden surface removal at the same time.

Incidentally, this algorithm is of mainly historical interest. It is very far from how modern GPU rasterizers or even software rasterizers work. Nowadays we pretty much always have a framebuffer and use double buffering, so we don't need to generate the pixels precisely in order anymore; consequently, we usually iterate over the triangles first, then find which pixels each triangle covers, rather than the other way around.

Moreover, the scanline algorithm is pretty inherently serial and doesn't make good use of today's highly parallel hardware. Modern approaches rely on multicore processing and SIMD to test many pixels against many triangles at once, rather than incrementally updating sorted lists of edges and suchlike.

For more on modernized software rasterization, see Fabian Giesen's articles Triangle rasterization in practice and Optimizing the basic rasterizer.


How does the Waves class in the book Introduction to 3D Game Programming with DirecX11 work?

This implementation represents the water as a regular height field grid. The function parameters m and n are the dimensions of the grid, dx is the regular distance between the water vertices, dt is the time step, speed is the speed of the movement and damping is a constant to ease disturbance as the water flows from the disturbance center.


Ray Tracing, why does the GGX make my rendering result darker and darker for each frame?

http://graphicrants.blogspot.com/2013/08/specular-brdf-reference.htm - here is reference with different geometry shading equations
http://jcgt.org/published/0003/02/03/paper.pdf - paper explaining how G term works and on what assumptions it is build.

Looking again at your code, I'm not sure (based on first link to Karis blog) if your G term is calculated correctly.
  float k = pow(roughness + 1, 2) / 8.f;
  k = pow(roughness, 4);

You're overwritting k with roughness^4, which doesn't fall into any Schlick equation that can be seen in Karis blog and if I recall correctly, I haven't seen that in any papers. Another note is, that you're passing ωo and ωi to your geometry function, but based on your comment that this is Schlick G and overall context, I think that you should use NoL and NoV (dot first pointToLight direction and then viewDir).
I was trying to do my own image, but I found [Blinn 1977] paper which sums up how Smith G model works (take a look at figure 1) - https://www.microsoft.com/en-us/research/wp-content/uploads/1977/01/p192-blinn.pdf#page=2&zoom=auto,-264,459
Some time ago, I was playing a little with different D, F, G equations - https://github.com/komilll/dxFramework/blob/master/dxFramework/Shaders/PS_BRDF.hlsl - take a look at PS_BRDF.hlsl and PS_BRDF_Helper.hlsl for further reference. I was basing my equations on Karis blog, some are modified base on other papers e.g. Frostbite and Disney diffuse are based on respective papers if I recall correctly.

I'm not sure if that's causing a problem, but it's worth fixing for sure. Also I don't think that you're using for summing frames
(previous_average * frame_count + current_color) / (frame_count + 1)

is correct, neither for uniform or importance GGX sampling. You should derive pdf in that case, it can be found in link that I've sent you before - http://cwyman.org/code/dxrTutors/tutors/Tutor14/tutorial14.md.html - it's generally nice tutorial, with good step by step introduction to DXR.


Correcting my "look-at" matrix so that it works on non-camera objects?

You don't seem to get the concept right:

A look-at camera matrix is created by defining two points, the "from" point (the dog eyes for example assuming you take the position of the eyes of the dog and average them to get the point in between) and a "to" point which is where the dog looks at. So if your dog is at (3,3,3) and that you want it to look at the some food which is at (-10, 100, 20), then you need to call the lookAt() method with these 2 points and it will create a matrix that will actually orient and position the camera properly.

Now you are getting confused (probably) because most cameras in most 3D packages are positioned at the origin (by default) and looking down the z-axis (by default). You can create such camera with a lookAt() function btw that would look like that:

lookAt(Vec3f(0), Vec3f(0,0,-1));

For example. Now check what the matrix this call will produce with the matrix of the default camera in Maya. You will notice they are different. The second is an identity matrix, the first one is not.

Now for example if you create a default camera in Maya it will be position at (0,0,0) and looking down the negative z-axis. If you look at the 4x4 matrix of this camera you will get identity. And maybe what you are trying to do is to create a 4x4matrix using the lookAt function and multiply this matrix by the Maya's default camera matrix in which case, your camera will point in the wrong direction (it will be flipped along the camera's local coordinate z-axis). That's because the matrix returns by lookAt properly orient the camera (as suggested in the first paragraph) but since the default camera in Maya is already oriented along the negative-z axis, your Maya camera will look the wrong way.

So in summary:

If you use the lookAt method in a renderer, it will work as you expect. There's no z problem.
if you try to use it on a camera that is already oriented along the negative z-axis and yet has identity matrix yes indeed you will be looking the wrong way. This can be fixed easily by flipping the camera back one more time along the z-axis. Simply create another matrix that scales the object/camera by (1,1,-1).

EDIT

The OP says (see comment):

  I did spot what looks like a typo: it says "(mind the direction of this vector: it is To − From not From − To)", but then you go on to write the line of code Vec3f forward = Normalize(from - to);, which contradicts what you were trying to clarify.

But a note in the reference says:

  The camera is aligned along the negative z-axis. Does it mean I need to rotate the camera by 180 degrees along the y-axis or scale it up by -1 along the z-axis? Not at all. Transforming a camera is no different from transforming any other object in a scene. Keep in mind that in ray-tracing, we build the primary rays as if the camera was located in its default position. This is explained in the lesson Ray Tracing: Generating Camera Rays. This is when we actually reverse the direction of the rays. In other words, the z-coordinates of the rays' directions at this point are always negative: the camera in its default position, looks down along the negative z-axis. Those primary rays are then transformed by the camera-to-world matrix. Therefore, there is no need to account for the default orientation of the camera when the 4x4 camera-to-world matrix is built.

So if you look at the reference they use by default a right-hand coordinate system, which means x-axis points to the right, y is up and z goes out of the screen (outward, going your direction when x-axis points to the right). So they build a matrix that "follows this convention" which makes sense. Then they say in the note "yes but wait, normally the camera points down the negative z-axis, so you should be doing to - from and not from - to, and they explain that in the case of ray-tracing the fact that the camera looks down the z-axis is performed when the rays are cast (by inverting their z-coordinate).

Vec3f rayDirectionInCamSpace{ x, y, -1 }; // invert here
normalize(rayDirectionInCamSpace);
Vec3f rayDirInWorldSpace = transform(rayDirectionInCamSpace, camToWorldMat);
pixelColor = trace(rayOrigWorldSpace, rayDirInWorldSpace);

So no it's not a typo. Now the reference you are pointing out maybe assume this matrix needs to be used in some OpenGL render in which there is no possibility to invert they ray direction when they cast. So they indeed use to - from instead which is the same as multiplying the cam-to-world matrix you would get from the reference by a matrix that flips the coordinates z-axis (multiply by (1,1,1-)).

So no typo, things are quite logical in both cases. The result you get is always constrained to a lot of small details, such as the handiness of the coordinate system and how is the matrix later used down your graphics pipeline (OpenGL pipeline, ray-tracer, renderer, animation package, etc. which may follow its own conventions, etc.).

EDIT 2

Note that in the reference your provide (screenshot of book) they do target - eye but then they use the negative of that in the matrix construction. Not sure why. Then also no idea why the multiply the eye by the vectors for the translation part of the matrix. It seems to be using row-order matrices for the axis of the coordinate system but then the translation part should be at mat[3][0], mat[3][1] and mat[3][2] and not mat[0][3], mat[1][3] and mat[2][3]. Not sure what book this is but seems to be terribly wrong and confusing.

EDIT 3

Code (that seems to be the only solution):

void lookAt(const Vec3f& from, const Vec3f& to, Matrix44f &V)
{
    static const Vec3f tmp(0, 1, 0);
    Vec3f forward = normalize(from - to);
    Vec3f right = crossProduct(normalize(tmp), forward);
    Vec3f up = crossProduct(forward, right);

    V[0][0] = right.x;
    V[0][1] = right.y;
    V[0][2] = right.z;
    V[1][0] = up.x;
    V[1][1] = up.y;
    V[1][2] = up.z;
    V[2][0] = forward.x;
    V[2][1] = forward.y;
    V[2][2] = forward.z;

    V[3][0] = from.x;
    V[3][1] = from.y;
    V[3][2] = from.z;
}

int main(int argc, char **argv)
{
    ...
    lookAt(from, to, camToWorldMatrix);
    worldToCamMatrix = camera.camToWorldMatrix.inverse();

    Matrix44f projectionMatrix;

    Matrix44f openGLWorldToCamMatrix = worldToCamMatrix;
    Matrix44f openGLProjectionMatrix = projectionMatrix;

    while (!glfwWindowShouldClose(window)) {
        glViewport(0, 0, width, height);

        glUniformMatrix4fv(projectionLocation, 1, GL_FALSE, (const GLfloat *)&openGLProjectionMatrix[0][0]);
        glUniformMatrix4fv(viewLocation, 1, GL_FALSE, (const GLfloat *)&openGLWorldToCamMatrix[0][0]);

        ...
        wwapBuffers(window);

    }

    return 0;
}


Performance impact during texture layer loading

From the way the problem only manifests when you start overwriting previously used texture tiles, I suspect you're running afoul of some driver heuristic that's attempting to manage the lifetime of the texture data on the GPU.

In general, when you update buffer or texture data from the CPU, the driver has to make sure it doesn't overwrite anything that might still be in use by the GPU (considering that the CPU and GPU run asynchronously to each other). There are various ways the driver can do this, such as internally double-buffering the texture, inserting sync points (for the CPU to wait on the GPU), or initially copying data to a temporary staging buffer and then later re-copying it to its final destination. All this stuff is under-the-hood and the driver makes its own decisions about what strategy to employ in any given situation, using some heuristics based on texture size, usage patterns, etc.

I can only guess at the details, but I suspect what's happening is that the driver is getting into a state where it thinks those previously-used texture tiles are still in use on the GPU somehow. Then when you start to overwrite them, it decides that it needs to make a second instance of the whole texture, and laboriously copies everything over before finally doing the updates you wanted—or something like that. That could be the many tiny draws you see in Nsight.

It's still a mystery why it's taking tens of milliseconds to do that, though. Maybe the textures are large enough that a second copy of them causes you to run out of VRAM? In that case, you'd get a lot of thrashing as things are spilled to main RAM. That's just a guess, though. As for why the bad behavior only happens for a few tiles and then goes back to normal, it could be that the driver heuristics are deciding to switch to a different strategy at that point.

Anyway, as for how to fix it, you might try using glInvalidateTexSubImage whenever you free a texture tile. This lets the driver know that you don't care what was in that texture tile before, which hopefully will convince it to just update the texture in-place rather than doing anything crazy with it.

Another thing to try is to make sure you don't try to re-use a texture tile in the very next frame after you free it, but give it a frame or two of "downtime" before re-using it. But it sounds like that's probably the case already, given that you prefer unused tiles before you start re-using them.

If that doesn't help, I'd recommend sending an email to devsupport@nvidia.com. The OpenGL driver team there can likely take your Nsight captures and help diagnose the issue, and fix it if it's due to a driver bug.


Deriving blur from real optical formulae

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

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

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

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


How to assign a texture to a single entry of a texture array?

glTexSubImage3D(GL_TEXTURE_2D_ARRAY, 0, 0, 0, 0, 894, 894, 0, color_format, GL_UNSIGNED_BYTE, data);

You're transmitting 3D data. But you told OpenGL that the 3D block of data has a depth of 0. So the number of pixels you're transmitting is... 0.

I have no idea why OpenGL doesn't provide an error when passing 0 for any of the sizing components (since that means you're not passing data). But that's your problem. You need to pass 1 for the depth if you're uploading a single 2D slice to the array.


Straight from the specification of lineTo:

  If the object's path has no subpaths, then ensure there is a subpath for (x, y).

And "ensure that there is a subpath" links to this bit:

  When the user agent is to ensure there is a subpath for a coordinate (x, y) on a path, the user agent must check to see if the path has its need new subpath flag set. If it does, then the user agent must create a new subpath with the point (x, y) as its first (and only) point, as if the moveTo() method had been called, and must then unset the path's need new subpath flag.

In short, if you call lineTo without an explicit beginPath, then a new path is created implicitly, with the given point as its starting point.


You can render with a raymarcher and using distance fields. For the hat you would define the center and radius of the cylinder, and would draw any march which falls within such circle, no matter what "height" you are in, it will be infinite. You can also add a restriction, like y > 0 to put a beginning to your hat.


Worley noise, also known as cellular noise, has the same property. It just as easily implemented as Perlin noise and easily extends to higher dimensions. Thus the slicing of 4D Worley noise will produce a 3D Worley noise. However, it is not necessarily a noise function but rather a texture function, producing cellular-like characteristics. With FBM applied to Worley noise some interesting patterns can be created.

You can find additional types of procedural noise functions in this survey paper.


DFG pops up in the family of microfacet based BRDFs. It is simply the product of three terms:

D : The microfacet distribution.
F : The fresnel coefficient.
G : The geometric attenuation between microfacets.

When someone says Cook-Torrance, they usually mean a microfacet BRDF where the distribution (D) is Beckmann, which I think is what the original Cook-Torrance paper was about.

GGX is really just a different distribution (D term). GTR is another one. The names are often used as a shortcut to mean a BRDF using that distribution.

Now all of these assume specular reflection (or refraction) on the microfacets, hence the fresnel term. If you take a similar idea but use diffuse (lambertian) microfacets, you get the Oren-Nayar BRDF. I think that one is based on a gaussian distribution of the microfacets but I'm not familiar with the details.


At which stage of the pipeline is the projection transform applied?

I know that in a typical GLSL shader you need to apply the modelview matrix to the vertices, so that you obtain a vertex in eye coordinates.

This is not quite correct, you need to apply the modelviewprojection matrix and thus transform your vertices into Clip Space Coordinates.
You can apply the ModelView Matrix in certain instances, e.g. if you want to calculate light in Eye Space. However, the gl_Position will then still not use the ModelView but the ModelViewProjection Matrix and your light/whatever- position will use the ModelView Matrix and be transferred to the fragment shader via the out variables.

For a more detailed discussion, you can look into the comments on this question, I had to be corrected about that myself (since not all literature is clear/correct about Normalized Device Space and Clip Space).


It's the normalization factor to make sure the BRDF always reflects the same amount of energy regardless of the value of the specular exponent n. Without that factor, changing the specular exponent changes the overall reflectance of the material.

Further reading: The Blinn-Phong Normalization Zoo | The Tenth Planet Blog


Rotating scene or complex object using center of AABB

Solution:

Instead of re-calculating center, pre-defined center is used for primitives. The pre-defined center is center of AABB if other value is not supplied. Averaging all center of mesh primitives gives mesh center and averaging all center of meshes gives scene center. After a mesh is transformed, center also gets transformed.

Now center is not moving after rotations. Using center of scene's AABB seems not good idea. Because each shape of child AABB is changing after transformations.


What is this vertical line contour effect in glitch art called?

The technique is called halftone.  Achieving continuity must be done through subsampling / super-sampling, but the essence of the algorithm goes like this.
Imagine walking across a horizontal line of the image accumulating energy from each input pixel as you go.  Once you have accumulated some threshold amount of energy in the accumulator, emit a single pixel, and then reduce the energy stored in the accumulator by a constant amount.
In code, that would be something like:
for (int y = 0; y < height; ++y) {
    double accumulator = 0.0;
    for (int x = 0; x < width; ++x) {
        accumulator += input[x, y];
        if (accumulator >= 0.0) {
            output[x, y] = 1.0;
            accumulator -= K; // try values like 5.0 or more for K.
        }
    }
}

The larger K is, the more energy must be accumulated before outputting a pixel, and hence, the more horizontal distance there will be between the lines.  If K is 1, that's actually a naïve approach to the dither algorithm.
Supersample to achieve continuity.  You can also emit discrete (x,y) points instead of rendering an output pixel directly and then render a new image by rendering a vertical curve passing through the nth point generated of each scanline.
Perform the operation discretely on separate colour channels to get what you described as the "chromatic aberration" effect.
The effect degenerates the more complexity you encounter in the image on a single line.  For example, after just two eyes on the same line, the lines are extremely wavy at the right of your example image.  Here is the effect exaggerated even more.


How to set flags and usage glnamedBufferStorage in new OpenGL 4.4+ (GL_STATIC_DRAW)?

The usage flag in BufferData was basically useless. There is no functional difference between each value and drivers used other ways to decide how to handle its data.

The flags parameter on the other hand actually changes what you can do with the buffer.


What physical properties are “lacking” to keep this 3D scene from looking like a real photograph?

Few things, but usually this is what it takes to make the difference:

1- the material reflection at his head, he is bald, yet the diffuse texture shows color difference where the hair is, this means he has a shaved head, not a natural bald, this should translate in reflection, take a look at his head, top right (top left for the image), the reflection is very smooth, no bumps at all, there should be some indication in reflection to the tiny hairs that you can see in the texture map.

2- the wrinkles in his face need to be emphasized, he looks like some kind of a worrier, yet he has very smooth skin like a 15 year old has, I'm not talking about deep scars, I am talking about skin texture and bumpiness due to rough way of life, skin bumps, moles, a minor infection that will affect both coloration and bumpiness.

3- his mouth is off, the size and positioning, I know that many people has small mouths, and a relatively large distance between the nose and the mouth, but that drives the attention of the viewers, it feels wrong, even if it was a photo of a real person, you can spare this kind of dragging attention to extreme shapes if reality is what you're after.

4- This is a minor thing, but it had a role to play, the background and outfit, they look fake, this drives the viewers automatically and subconsciously to feel that this is not real, try to put more effort in these, making the outfit less uniform, and the background to have more realistic pattern, with more hue variation.

5- And last, the eyes, there's no reflection, yet the lighting is low and the image is small, so I am not sure if this is due to lighting or to the material itself.

That all been said, all these remarks are minor, the artwork looks awesome, I don't think a CG artwork should look like a photo, a perfectly realistic CG artwork that looks exactly like a photo, will become just another photo to the viewer.


The output of an LCD or other display is continuous by virtue of its being a physical device, with pixels that have some size and shape, and emit some physical spectrum of light.

from Apple Retina Display by Bryan Jones

My guess as to what's meant by the comment that "the sampled representation is converted into a continuous one by the video hardware" is just that the input to the hardware is a sampled digital signal that drives the physical pixels, while the output is a continuous light field such as the above. So the "reconstruction filter", so to speak, is the physical pattern of light emitted when the display is turned on.

(Also note that displays often do additional digital processing of the input image before displaying it, e.g. to adjust brightness, contrast, color temperature, sharpness, and so forth. But those things wouldn't be called a reconstruction filter.)


Can I parallelise the scene for a raytraced image, rather than just the pixels?

Yes, such a thing is possible. There are a few different ways of doing it, but the basic idea is to split up the scene into chunks and assign these chunks to different machines. Each machine then traces rays locally within their assigned chunk.

The tricky bit is of course when rays leave a chunk, in which case machines need to synchronize and exchange rays that cross the local chunk boundaries. A recent paper from this year's EGSR, Distributed Out-of-Core
Stochastic Progressive Photon Mapping, describes how to do this efficiently for Progressive Photon Mapping, although other light transport methods might be feasible.

Although such distributed ray tracing is possible, it is certainly not the most efficient form of parallelization, and non-uniform lighting directly translates into non-uniform distribution of the workload. It is really only beneficial when the size of the scene far exceeds the memory of a single machine.


Tasks of a professional graphics programmer?

I don't know if this is exactly what you're looking for. I work tangentially in the film and TV industry. I don't work for a studio, but I work on the software that studios use for their productions. It's software that any user could buy if they wanted to, but it has been used in feature films, television shows, commercials, music videos, and more.

About half of my time is spent programming and debugging. This is the fun stuff. It involves implementing rendering and filters of various types, sometimes based on papers I've read from SIGGRAPH or cool web sites or book. I don't do very much user interface work, though there is some. I work with OpenGL and Metal mostly. I need to know image processing algorithms in addition to general computer science algorithms, like graph processing and traversing, for example. This includes dealing with source control, continuous integration, stepping through code to find and fix bugs, talking to QA about problems, etc.

The other half of my time is spent doing non-programming things. This is a variety of things like attending conferences to learn (things like the Game Developers Conference or SIGGRAPH), conventions to promote the company and/or product (things like NAB, IBC, etc.), answering tech support emails and writing sample code (we have an API so 3rd parties can extend our product), attending meetings (scrum, management updates, company-wide meetings, etc.), talking with marketing about how our new features work, meeting with customers to understand how they use the product, explaining to our documentation people how the product works so they can document it, etc.

This seems like a pretty similar breakdown to other programming jobs in other industries from talking to my friends who work in other areas.


Camera rotation effect on pitch, yaw, roll

Thanks for the suggestion to look up Euler angle rotations.  https://www.learnopencv.com/rotation-matrix-to-euler-angles/ is a good reference to answer this question.


Calculating Skeleton-Deformed Mesh Normals

You're right, standard skinning pipelines don't produce correct normals. You get similar problems when using translation on joints, which is very common in facial animation. The only way to get a correct result is to take the average of the surrounding face normals - this is not information you will be able to get in the vertex shader, so you'll have to use a geometry shader or run your skinning on CPU or a compute shader.


What is the fundamental reasoning for anti aliasing using multiple random samples within a pixel?

From a signal processing point of view, you're sampling a continuous-domain signal, and you need to filter it to get rid of frequencies beyond the Nyquist limit. It's that filtering that leads to integrating over the pixel area—or more generally, integrating over the support of your antialiasing kernel (which need not be a box).

Consider your rendering function that takes a sample point $x, y$ in screen space and returns back the color found at that point. (Let's ignore any issues of random sampling for the moment, and assume it returns a "perfect" rendered color for that specific point.) This function effectively defines a 2D continuous-domain signal. Or to put it another way, it defines an infinite-resolution image, since nothing prevents this function from having features at arbitrarily small scales. In terms of frequency domain: the function is not band-limited; it can include components of arbitrarily high spatial frequencies.

Now you want to convert it down to a finite number of pixels. Much like digitizing an audio signal, when you sample it, you'll get aliasing unless you first eliminate frequencies beyond the Nyquist limit imposed by the sampling rate. In other words, you have to get rid of features smaller than the pixel grid. To do this, you apply a low-pass filter. The ideal low-pass filter is the sinc function, but for various reasons of practicality we use other filters (that don't perfectly eliminate frequencies beyond the Nyquist limit, but they at least attenuate them).

Low-pass filtering is done by convolution. If $f(x, y)$ is the rendering function and $k(x, y)$ is a filter kernel, then we can define the low-pass filtered image as:

$$f_\text{filtered}(x, y) = \iint f(x', y') \, k(x' - x, y'- y) \, dx' \, dy'$$

Then the image is safe to sample, so the final pixel values are gotten by just evaluating $f_\text{filtered}$ at the pixel coordinates.

If $k$ is a box filter, which looks like $k = 1$ within the pixel box and $k = 0$ elsewhere, then this simplifies to just integrating $f$ over the pixel box. But as noted, the box filter isn't so great and there are better choices such as tent, bicubic, and Gaussian filters.

Anyway, now that we have an integral, we can use Monte Carlo for it and combine it with any other integrals we might like to do—for lighting, motion blur, and so on. We can even apply importance sampling to the $k$ factor in the integral, by generating samples for each pixel that are distributed around the pixel center according to $k$.


Help with wireframe artifacts on mesh procedural UV texturing

This looks like a side effect of bilinear filtering where it can sample texels neighbouring those that are at the edges of the used areas of the texture - the red/green boundary in your case.

One solution for this is to perform a "colour bleed" pass on the texture where for every 'empty'(Green) texel with none-empty (Red) neighbours you can output the average of all the none-empty neighbours, so what this will do is add a 1 pixel boundary to the red part of your image that matches its neighbours (In this case the average should be red)

This is a solution to a similar problem seen with alpha blended textures with sharp edges where colour bleeding can sometimes appear at the boundaries between semi-transparent and completely transparent pixels when magnifying the image. Most modern art packages will automatically colour bleed the RGB channels when saving RGBA images for this reason.


This is achieved via Barycentric Interpolation.

First, we find the barycentric coordinates of $P$. Barycentric coordinates represent how much weight each vertex contributes to the point, and can be used to interpolate any value which is known at the vertices across the face of a triangle.

Consider the 3 inner triangles $ABP$, $PBC$ and $PCA$.

We can say that the barycentric coordinate, or weight of the vertex $A$ on the point $P$ is proportional to the ratio of the area of inner triangle $PBC$ to the area of the whole triangle $ABC$.

This is intuitively evident if we consider that as $P$ approaches $A$ the triangle $PBC$ grows larger and the other two become smaller.

Also intuitively evident should be that the sum of the barycentric coordinates of a point inside a triangle always equals $1$. So, it is enough to find only two of the coordinates to derive the 3rd.

The method for computing the barycentric coordinates is:

$$\begin{aligned}
{Bary}_A &= \frac{(B_y-C_y)(P_x-C_x) + (C_x-B_x)(P_y-C_y)}{(B_y-C_y)(A_x-C_x) + (C_x-B_x)(A_y-C_y)}\\
{Bary}_B &= \frac{(C_y-A_y)(P_x-C_x) + (A_x-C_x)(P_y-C_y)}{(B_y-C_y)(A_x-C_x) + (C_x-B_x)(A_y-C_y)}\\
{Bary}_C &= 1 - {Bary}_A - {Bary}_B
\end{aligned}$$

The derivation and reasoning is explained in the wikipedia article.

Once you have coordinates, you can determine the texture coordinates of $P$ by interpolating the values at the vertices using the barycentric coordinates as weights:

$$P_{uv} = {Bary}_A \cdot A_{uv} + {Bary}_B \cdot B_{uv} + {Bary}_C \cdot C_{uv}$$

The reasoning is also explained very nicely in this presentation.

Also see this question for efficient methods of computation.


How do I use ColorSync Utility to convert between Display P3 and sRGB?

Assuming that "Display P3" refers to the DCI-P3 colour space, this covers a wider gamut of colours than sRGB, so that an exact match is not possible for P3 colours that fall outside the gamut of sRGB. The example colour in the question is the reddest red that DCI-P3 can produce. This is considerably redder than the reddest red that sRGB can produce, so the nearest available red is shown instead.

You are right that (1, 0, 0) in sRGB is not the same as (1, 0, 0) in P3, but it is the closest it is possible to get. If you enter (1, 0, 0) in sRGB, you should see that converting it to P3 gives a red that is not the extreme (1, 0, 0). Any colour in sRGB can be converted to P3, which covers all of sRGB and more. However, only some of the colours in P3 can be converted to sRGB, since sRGB doesn't cover all of those colours. For colours that sRGB doesn't cover, the conversion utility used appears to be giving the nearest approximation available within sRGB.

It is analogous to converting RGB to greyscale. Any shade of grey in RGB can be converted to greyscale and still look the same, but outside of the set of greys, any other colour will be changed when you convert to greyscale. So red will be converted to a grey that best approximates it.

In answer to your final sentence, if you want to take an sRGB colour and find the most saturated form of that colour that P3 can represent, then it is going to be less straightforward when you are dealing with a colour other than one of the primaries red, green, and blue. You need to consider whether you want the nearest hue amongst the maximum saturation colours in P3, or the most saturated of the exact hue match colours in P3. These are subtly different, because for a given hue its maximum saturation available within P3 is less than the maximum saturation available for the primaries, since the colour space is a triangle. You can compromise by deviating only slightly from the desired hue to achieve slightly more saturation, but the degree of compromise will depend on your purpose.

Even for the primaries, you need to ask whether you want to convert the greenest green sRGB can represent to the greenest green P3 can represent, or whether you want to take into account that those two colours are not quite the same hue. since sRGB and P3 don't use quite the same primaries (only the blue primary is identical). If you want the hue to remain constant, it might help to convert to an intermediary colour space that includes hue as a parameter.


How do we approximate the integral of the measurement contribution function in Metropolis Light Transport?

This approximation is typically done by running a bidirectional path tracer with a modest amount of samples per pixel. This means multiple paths per pixel to approximate the integral of $f_j$, the measurement contribution function. Veach's original MLT paper explains how it can be done in a way that eliminates start-up bias (see Section 5.1).


What is the difference between various space-filling curves?

The difference is how well a mapping preserves locality and how easy it is to encode/decode the keys. The paper "Linear Clustering of Objects with Multiple Attributes" by H V Jagadish says: "Through algebraic analysis, and through computer simulation, we showed that under most circumstances, the Hilbert mapping performed as well as or better than the best of alternative mappings suggested in the literature". On the other hand, z-order is a bit simpler to use, for example compare the various methods listed in Bit Twiddling Hacks for z-order and Wikipedia for Hilbert-order.

As for the applications, I think the main advantage in using space filling curves is that they map points from higher dimensional space to space of lower dimension. For example, they make it possible to window query for points using traditional B-tree database index. Again, on the other hand, the disadvantage is that one needs to know the bounds of the input in advance as it is difficult to "resize" the mapping later.

PS: "Z-curve" is the same as "Morton code".

PPS: Additional mappings include Peano curve and for applications see also Geohash.


Can't understand the Importance sampling GGX

In microfacet BRDFs, the half-vector is the same as the microfacet normal. The half-vector is exactly the required normal for a microfacet to reflect light from the incident ray to the outgoing ray, by the law of reflection. So, for given incident and outgoing directions, only those microfacets whose normals are aligned with the half-vector are "active".

Conversely, when importance-sampling a microfacet distribution, we wish to select a random microfacet normal, then reflect the incident ray through this normal to produce the outgoing ray. The incident and outgoing rays will then have the chosen microfacet normal as their half-vector.

For an isotropic BRDF, you will sample theta from the inverse CDF, and you must also sample phi, the azimuth angle about the normal, uniformly between 0 and 2π; then convert those into a unit vector using the usual spherical coordinates (as mentioned in the linked blog post), combined with a tangent-to-world transformation. For an anisotropic BRDF, both theta and phi would be supplied by sampling an inverse CDF.


How to compute normal of surface from implicit equation for ray-marching?

The usual way of doing this with raymarching is to define your surface as a signed-distance field and use finite differences to get the gradient of the distance function at the point you’re sampling. In other words, if you have a function map(p) that returns the signed distance value at a point p, the normal at p is given by:

float epsilon = 0.001; // arbitrary — should be smaller than any surface detail in your distance function, but not so small as to get lost in float precision
float centerDistance = map(p);
float xDistance = map(p + float3(epsilon, 0, 0));
float yDistance = map(p + float3(0, epsilon, 0));
float zDistance = map(p + float3(0, 0, epsilon));
float3 normal = (float3(xDistance, yDistance, zDistance) - centerDistance) / epsilon;

Note that you can improve the accuracy by taking an extra 3 samples on the other “side” of the center position and subtracting those from the first three, then dividing by 2× the epsilon value, but if your distance function is complicated then that can be a lot of additional work for not much visible benefit.

The sample function you provide isn’t a signed-distance function as-is, and I don’t know of a straightforward way to convert it into one, but there’s an article by Íñigo Quílez (kind of a legend in this field) here with various functions you can make use of, including one for an ellipsoid. There’s also a lot of reference material available on Shadertoy—search for “raymarching” and you’ll find plenty of examples.


Operator associativity and 4x4 matrices - performance question

As Olivier mentioned, it's not safe in general to apply algebraic optimizations to floating-point math, as floating point doesn't actually obey associativity rules due to roundoff error. $(a + b) - b$ will not be equal to $a$, if the intermediate result $(a + b)$ had to be rounded; nor will $a \cdot (b + c)$ be equal to $a \cdot b + a \cdot c$ in general, for the same reason.

That being said, shader compilers are typically very relaxed about floating-point math optimizations. They behave like using -ffast-math on a C++ compiler, where they make various simplifying assumptions—including that algebraic operations are associative—in order to optimize more aggressively. So it's not out of the question that a GLSL compiler could actually make the optimization you describe.

That said, GLSL compilers vary wildly in quality; there is no canonical GLSL compiler, every GPU vendor implements their own (for OpenGL at least—in the case of Vulkan there are a couple standardized options), and who knows what level of optimizations they might or might not apply under various circumstances. If you wish the code to be optimized in a particular way, better to write it explicitly that way.

One other point,

  the second version should (at least in theory) be much cheaper to compute

It's only cheaper if you're transforming just a single vertex. If you want to transform thousands of vertices by the same matrix, it will be cheaper to multiply all the matrices together up front; then each vertex only needs to be transformed by a single 4×4 matrix instead of several.


How to record N-body simulation as a video?

user1118321's suggestion is the correct one but i'll elaborate a bit.

You run your simulation giving it a fixed time step between frames.  This time step would be 1/30th of a second (33.3333 ms) for instance if you wanted a 30 fps movie.

You let your simulation calculations and your rendering run, taking as long as they need to, and then write the rendered frame out to disk, with the frame number in the name (like frame0.bmp, frame1.bmp, etc).

After you have all these frames, you can use ffmpeg to combine them into a video of the desired frame rate.

If you want audio in your video (like sounds when things collide), you can actually keep an audio timeline (what audio events happened on which frames, or at what points in time) while doing your simulation, and then after your rendering is finished, you can generate the audio stream that goes along with your video, which ffmpeg will gladly combine with your video.

I've used this on a few projects in times past and it's worked great.  I also wrote up some more details about this process a while back:

Recording lagless demo videos of a laggy game | The blog at the bottom of the sea


Refraction: given an incoming and transmitted direction, can I deduce the normal?

See equation (16) in Microfacet Models for Refraction through Rough Surfaces : $-(\eta_i w_i + \eta_o w_o)$ which you'll probably want to normalize. $\eta$ are the two indices of refraction and $w$ the two vectors (they use $i$ and $o$ in the paper). Also look at the left part of figure 7 on the same page to see where all the vectors are pointing.


Is it possible to draw the intersection of two squares of known dimension using only x,y coordinates and without an if statement?

Basic calculations
I'm assuming you want to calculate the bounds of the rectangle of intersection, so you can draw it later. As AlanWolfe mentioned, you will want to use the max and min functions. This is so you can retrieve the left-most, right-most, top and bottom squares. We'll call these placeholders right, bottom, left, and top (each of these will point to one of the two squares).
Now that you know these, you can start mathing. The top-left position of the intersection will be at (right.x_left, bottom.y_top). The bottom-right position of the intersection will be at (left.x_left + left.width, bottom.y_top + bottom.height).
Rendering the intersection
To draw this rectangle (unconditionally), first you will need a location and a dimension, because most draw functions require these.
Intersection.width = Intersection.x_right - Intersection.x_left

and the same for the height. Now suppose the squares don't overlap. This would mean at least one dimension is negative. We don't want anything to draw in this case, so we will set any negative dimensions to zero.
Intersection.width = max(Intersection.width, 0)

and the same for the height. And to wrap up, we render our rectangle:
fillRect(Intersection.x_left, Intersection.y_top, Intersection.width, Intersection.height)


Yes, your understanding it correct. The Laplace-Beltrami depends on the current state $x$ and you have to recompute $L$ if $x$ changes. Therefore you cannot write matix of $\Delta$ without knowing the state $x$.

To expand on solving $\frac{\partial x}{\partial t} = -\Delta x$:

Discretization in space turns $\Delta$ to $L$, discretization in time turns $\frac{\partial x}{\partial t}$ to $\frac{x_{n+1} - x_n}{h}$ and decide if we evaluate the right hand side at time $t_n$ or $t_{n+1}$

Option 1 = Forward Euler:
$$
x_{n+1} - x_n = h L_n x_n
$$
where $L_n = L(x_n)$, $L$ is considered as a function of the state $x$. As you have mentioned this method can be unstable if we take too big step $h$. Therefore, we often resort to the second option.

Option 2 = Backward Euler:
$$
x_{n+1} - x_n = h L_{n+1} x_{n+1}
$$
Because $L_{n+1}$ depends on $x_{n+1}$ the above equation is non-linear system, we cannot solve it with a simple linear solve.

The common approach is to Taylor expand $L$ around the point $x_n$, general Taylor expansion around point $x$ is
$$
L(x + \Delta x) = L(x) + \nabla L(x) \cdot \Delta x + \dots
$$
therefore for $x= x_n$ and $\Delta x = x_{n+1} - x_n$ we get
$$
L_{n+1} = L(x_{n+1}) = L(x_n + \Delta x) = L(x_n) + \nabla L(x_n) \cdot \Delta x + \dots
$$

The simplest thing to do is to take just the first term of the Taylor expansion, i.e. we replace $L_{n+1}$ with $L_n$. This way, we end up with a linear problem
$$
x_{n+1} - x_n = h L_{n} x_{n+1}
$$

You can take higher order terms or solve the Backward Euler equation with non-linear solver. However, in your case that is probably unnecessary and the headache of computing the gradient of $L$ is not worth it(I might be wrong here, I do not have hands on experience with solving this equation).


Render fragments with uvs outside (0,1) range as transparent

Assuming that alpha works properly in the rest of the code it can be solved on shader level:
out vec4 color;
in vec2 uv;
uniform sampler2D tex;

void main() {
    color = texture(tex, uv);
    color.a = 0.0;
    if (all(greaterThan(uv, vec2(0.0))) && all(lessThan(uv, vec2(1.0)))) {
        color.a = 1.0;
    }
}


How do you loop through the contents of a buffer in DX11?

There are a couple of issues.  First you need to make sure that veBuffer has the correct type, so that probably means you need to cast it to pointer to vertex:
const auto vertices = reinterpret_cast<const Vertex*>(reBuffer);

Second a range-based for-loop needs to know where the range ends and giving it a raw pointer does not provide that information, so you cannot use it here.  Instead you will need to know how many vertices there are and use a normal for-loop:
const size_t numVerts = 32;
for (size_t i = 0; i < numVerts; ++i) {
    const Vertex& vert = vertices[i];
    ...
}


Creating an "ink spread" effect with a glsl fragment shader

If you want to implement this with GLSL, the general concept that could be applied is multi-pass rendering. In a first render pass, draw some ink on the screen and write the result in a frame buffer object (FBO). Then in the second render pass, draw a screen-size quad and attach the texture of the FBO. In this second render pass you could use a GLSL shader that distributes the ink. In the example below I used a simple Erode and a Gauss filter but more complicated effects are possible here. The result is then used as the background for the first render pass that draws more ink, which is passed on to the second render pass, and so on...

I have generated the example above with an online tool called "GSN Composer", which uses GLSL shaders and WebGL. You find the corresponding node graph here:

https://www.gsn-lib.org/index.html#projectName=foruminkspread&graphName=erode

Currently, in this very simple implementation, the spreading of ink will never stop. To achieve this you would need to remember a per-pixel state (e.g. a "mass" or "reservoir", as you have mentioned in your question) and use this per-pixel state in the shader that performs the spreading. This state could be encoded in the texture output (e.g. in the alpha channel) or could be written into a second texture that is passed along the shader pipeline.


Fill the plane with pentagons as tightly as possible in a regular way

If I understand the general idea correctly, you are going to close the gaps which are appearing after further structure growth. I cannot do it algorithmically but geometrically it is quite easy to do. I use Illustrator here and manually cover the gaps with two additional shapes.

So in addition to the pentagon and flat rhombus, I will need a 5-fold star and another rhombus:

So after first steps we get such a structure and the filling will be:

(Initial center is marked red)

I don't fill all holes and I leave out the obvious ones just to make arrangements more visible.

Then next step will give this:

Here I do only halves of crack fillings, those are symmetrical. As seen those cracks indeed can be covered using same shapes without gaps and have self-similar patterns. How this can be achieved algorithmically, I don't know unfortunatelly. But this is interesting question, probably there are solutions somewhere, IDK.

Theoretically if one analyse several crack fillings for each step, one can guess a relatively
simple algorithm.


How does OBJ format determine vertex normals?

It is common for modeling software to compute the normals when they aren't included in the loaded file (regardless of the format). This is usually highly dependent on the settings chosen when the file is loaded. Vertex normals can be computed fairly easily. The software can also compute smooth or faceted normals so you can switch between the two.
Your example includes faces...the f 3 2 1 is a face. Each number refers to  one of the vertex values which are numbered starting from the top. v -2.22045e-16 -0.291306 -2.22045e-16 is vertex number 1 in this case. The software can easily use that information to compute normals.
A normal is computed from the cross product of the edge vectors. Stictly speaking the obj format doesn't compute the actual normals, it just provides the information that Meshlab needs so that Meshlab can compute the normals. which looks something along the lines of:
facted_normal = normalize( cross( 2-3, 1-3))
Where the numbers for 1, 2 and 3 are replaced with the actual values for those vertices. Other types of normals are usually computed based of that.
I loaded this file into Blender and it is well behaved, so odds are meshlab has some options turned on that are showing split normals or vertex normals or one of the many other possible variants out there.


`worldViewProjection` doesn't work for a single component

This was entirely my fault. I simply forgot to set the uniform locations before rendering!

I implemented a function setUniformLocations() which sets the uniform locations passed to it:

public static setUniformLocations(
    gl: WebGL2RenderingContext,
    uniforms: IUniform[],
    program: WebGLProgram): IUniform[] {

    const vao = gl.createVertexArray();
    gl.bindVertexArray(vao);

    for (const uniform of uniforms) {
        uniforms.find((u) => u.name === uniform.name).location = gl.getUniformLocation(program, uniform.name);
    }

    return uniforms;
}

I normally call this function in the constructor of the according class.

this._uniforms = GLHelper.setUniformLocations(this.gl, this._uniforms, this._program);

This time I forgot.


How to get the transform amount from a screen vector and a direction vector?

I would not recommend to correlate the extrusion direction with its screen projection because of the asymptotic behaviour in case the axis points toward the screen or the resulting extrusion crosses the focal point of the camera.

One consistent way of doing it would be computing the extrusion height independently of the axis direction:

Find a virtual construction plane of your operation which is parallel to the screen. For example, in the real world, this plane would be a plane of your computer monitor screen reachable by hand as you sit at your computer.

Usually, you have some kind of central point which serves as a center of screen camera rotation as you navigate your viewport. You can draw your construction plane through this point.

Another way is drawing your construction plane through the intersection of the mouse cursor and your object as you start your extrusion operation.

Either way, the construction plane must be parallel to the screen and in front of your camera view for this method to work.

If you are using an orthogonal projection (camera without a Field of View) then you can draw the construction plane through any point.

Compute the length of the projection of mouse motion on the construction plane your figured out in point (1). This is your desired extrusion height.
Provide a means of moving the construction plane back and forth along the view axis, if the current viewport has FOV. You don't need to display and the move the construction plane if your viewport uses orthogonal projection.

Using this method your extrusion operation height would be visually consistent with what user expects with mouse motion, and the user would not experience any asymptotic and singularity issues as with the axis projection method.


why is translating in 3D space the same as shearing in 4D space?

In a linear transformation system, your origin is always a fixed point, since 0*anything = 0. So imagine you have a cinema screen, and the origin is at the centre of the screen. Using linear transformations, you can rotate, scale or shear the image, what you can't do is move it, since you have a fixed point in the middle.

Now add a dimension, and move your origin out of the screen through this new dimension. Your origin is now at the projector, at <0, 0, 0>, and the centre of your screen is at <0, 0, 1>. So as well as the image transforms you had available before, you can now effectively move your image, by shearing through the space in between the projector and the screen.

Add another dimension and you have the 4D -> 3D case.


How to link Intel Open Image Denoise in Qt Creator and use it to denoise noisy Monte Carlo path tracing renderings?

All the previous answers/comments were correct. The weird thing about my case was the symlinks of the libraries in the build directory were not valid and were not pointing to the right files. The output of ls -lh /home/Warrior/Desktop/pathtracer/oidn/build was as follow:
total 160K
drwxrwxr-x 4 amir amir 4.0K Feb  4 00:09 apps
-rw-rw-r-- 1 amir amir  19K Feb  4 00:09 CMakeCache.txt
drwxrwxr-x 7 amir amir 4.0K Feb  4 00:09 CMakeFiles
-rw-rw-r-- 1 amir amir 6.5K Feb  4 00:09 cmake_install.cmake
drwxrwxr-x 3 amir amir 4.0K Feb  4 00:09 common
-rw-r--r-- 1 amir amir 3.6K Feb  4 00:09 CPackConfig.cmake
-rw-r--r-- 1 amir amir 4.0K Feb  4 00:09 CPackSourceConfig.cmake
-rw-rw-r-- 1 amir amir 103K Feb  4 00:09 Makefile
drwxrwxr-x 3 amir amir 4.0K Feb  4 00:08 mkl-dnn
-rw-r--r-- 1 amir amir 1.6K Feb  4 00:09 OpenImageDenoiseConfigVersion.cmake

But when I would go to the directory I would see the following:

What I did was instead of cloning the oidn repo somewhere else and building it, I cloned the oidn repo in /home/Warrior/Desktop/pathtracer/ and built it there. After that I could see the correct file names/symlinks via ls -lh as follow. Now my program compiles and things work fine. I just had to use LIBS += -L"oidn/build/" -lOpenImageDenoise when compiling my code:
total 51M
drwxrwxr-x 4 amir amir 4.0K Feb  9 10:14 apps
-rw-rw-r-- 1 amir amir  19K Feb  9 10:14 CMakeCache.txt
drwxrwxr-x 7 amir amir 4.0K Feb  9 10:19 CMakeFiles
-rw-rw-r-- 1 amir amir 5.6K Feb  9 10:14 cmake_install.cmake
-rw-rw-r-- 1 amir amir 2.5K Feb  9 10:17 color_ispc_avx2.h
-rw-rw-r-- 1 amir amir 2.5K Feb  9 10:17 color_ispc_avx512skx.h
-rw-rw-r-- 1 amir amir 3.1K Feb  9 10:17 color_ispc.h
-rw-rw-r-- 1 amir amir 2.5K Feb  9 10:17 color_ispc_sse4.h
drwxrwxr-x 3 amir amir 4.0K Feb  9 10:14 common
-rw-r--r-- 1 amir amir 3.3K Feb  9 10:14 CPackConfig.cmake
-rw-r--r-- 1 amir amir 3.7K Feb  9 10:14 CPackSourceConfig.cmake
-rw-rw-r-- 1 amir amir 2.6K Feb  9 10:17 input_reorder_ispc_avx2.h
-rw-rw-r-- 1 amir amir 2.6K Feb  9 10:17 input_reorder_ispc_avx512skx.h
-rw-rw-r-- 1 amir amir 3.1K Feb  9 10:17 input_reorder_ispc.h
-rw-rw-r-- 1 amir amir 2.6K Feb  9 10:17 input_reorder_ispc_sse4.h
-rw-rw-r-- 1 amir amir  41K Feb  9 10:14 libcommon.a
-rw-rw-r-- 1 amir amir  16M Feb  9 10:17 libdnnl.a
lrwxrwxrwx 1 amir amir   24 Feb  9 10:19 libOpenImageDenoise.so -> libOpenImageDenoise.so.0
lrwxrwxrwx 1 amir amir   28 Feb  9 10:19 libOpenImageDenoise.so.0 -> libOpenImageDenoise.so.1.3.0
-rwxrwxr-x 1 amir amir  34M Feb  9 10:19 libOpenImageDenoise.so.1.3.0
-rw-rw-r-- 1 amir amir  37K Feb  9 10:19 libutils.a
drwxrwxr-x 2 amir amir 4.0K Feb  9 10:17 local__core
-rw-rw-r-- 1 amir amir 104K Feb  9 10:14 Makefile
drwxrwxr-x 3 amir amir 4.0K Feb  9 10:14 mkl-dnn
-rwxrwxr-x 1 amir amir 181K Feb  9 10:19 oidnBenchmark
-rwxrwxr-x 1 amir amir  59K Feb  9 10:19 oidnDenoise
-rwxrwxr-x 1 amir amir 882K Feb  9 10:19 oidnTest
-rw-r--r-- 1 amir amir 1.6K Feb  9 10:14 OpenImageDenoiseConfigVersion.cmake
-rw-rw-r-- 1 amir amir 1.5K Feb  9 10:17 output_copy_ispc_avx2.h
-rw-rw-r-- 1 amir amir 1.5K Feb  9 10:17 output_copy_ispc_avx512skx.h
-rw-rw-r-- 1 amir amir 1.9K Feb  9 10:17 output_copy_ispc.h
-rw-rw-r-- 1 amir amir 1.5K Feb  9 10:17 output_copy_ispc_sse4.h
-rw-rw-r-- 1 amir amir 2.5K Feb  9 10:17 output_reorder_ispc_avx2.h
-rw-rw-r-- 1 amir amir 2.5K Feb  9 10:17 output_reorder_ispc_avx512skx.h
-rw-rw-r-- 1 amir amir 3.1K Feb  9 10:17 output_reorder_ispc.h
-rw-rw-r-- 1 amir amir 2.5K Feb  9 10:17 output_reorder_ispc_sse4.h
-rw-rw-r-- 1 amir amir 1.5K Feb  9 10:17 upsample_ispc_avx2.h
-rw-rw-r-- 1 amir amir 1.5K Feb  9 10:17 upsample_ispc_avx512skx.h
-rw-rw-r-- 1 amir amir 1.8K Feb  9 10:17 upsample_ispc.h
-rw-rw-r-- 1 amir amir 1.5K Feb  9 10:17 upsample_ispc_sse4.h
drwxrwxr-x 2 amir amir 4.0K Feb  9 10:17 weights


Estimating real world object distance from photo

Yes, theoretically, for the sphere no for other things in image. If it is possible to augment the situation then it is possible to measure anything you can track much more reliably. First, knowing the scale of a object that you can successfully track is a requirement to getting a scale for the image as the cameras lost all scaling info. But this requirement alone is not enough.

In practice you do not know the other settings of a camera (unless calibrated with stereo geometry first), all kinds of problems arise when you do real camera tracking. Such as real cameras lens not being ideal. By tracking more than 2 features on the scene you can automatically estimate the actual focal distances and camera distortions. More importantly you can get an error estimate on your data.

One way is to have two photos from different angles in a scene with enough depth variation (for a stationary object this is quite easy). You can then use stereo geometry to solve the depth many such tracker apps exist as part special effects packages and making your own is not entirely out of the question (fun link, fundamental matrix song). The second way is to have a known plane, the plane can be solved by knowing 4 points and their relative locations. In reality it is better use more points for error and lens correction estimates. A third way is to use structured light (a cloud of laser light points, or a line for example) to triangulate the point distances. Or you could use a light field camera.


How to handle specular highlights with reflection / total internal reflection in a ray tracer?

But, at this point, do i apply a specular highlight to the surface that i hit from the inside?

Unless you have a light source inside your object, there's no point in doing lighting on the inside surface as you'll hit the object again before you reach the light or environment.

If you want to support a light in the object, then by all means do your lighting but make sure to trace rays to check that you actually reach the lights without exiting the object.

Otherwise, that lighting should come from your (internal) reflected ray which will eventually hit the surface again, be refracted and get light from the outside.

  which ray should i use for the lighting - the refracted ray or the reflected ray, or both?

I think the confusion here is that when you write of "calculate a specular highlight" there is an assumed BRDF which may not make much sense for refraction. At least probably not without a few tweaks. What you really need is a BTDF (bidirectional transmittance distribution function).

If you want to keep it simple, I think you could adapt a Blinn-Phong model for refraction by using the half vector of equation (16) in Microfacet Models for Refraction through Rough Surfaces.


Yes, scene graphs can be non-tree graphs. Your question is a little vague, but assuming you're talking about something like the scene graph in rendering a frame of video, for say a feature film, then yes, there are often references between objects. For example, in a compositing application, you might have a piece of footage with an effect applied. The effect may have as one of its inputs another piece of footage which is also in the scene. (It might use the other image for blending, as a reflection map, as a diffuse or specular lighting map, a bump map etc.) So you could have a scene that looks like this:

scene contains:
    footage A contains:
        filter "Multiply Blend" which takes input:
            footage B
    footage B
    footage C

In something like a game, you might objects that refer to other objects that they can collide with:

scene contains:
    ground which can collide with:
        player
        ball
        goal
    player which can collide with:
        ground
        ball
        goal
    ball which can collide with:
        player
        ground
        goal

etc.


Importance sampling a cosine distribution

You can sample using absolutely any distribution you want, as long as you weight the results by dividing by the pdf of the sampled distribution. It will converge to the right answer (as long as the distribution is nonzero everywhere that you want to integrate). Different distributions will give different amounts of variance though. The trick with importance sampling is to find a distribution that minimizes variance while still being cheap to compute. A distribution similar to the shape of the function you're integrating, or equal to some factor of that function (so that it cancels when you divide), works best. In other words you try to guide the sampling toward areas that are more important to the result, hence, "importance" sampling.
So what the author is saying is the image on the left was sampled with $\omega_o$ (or maybe $\omega_m$, I'm not sure) drawn from a cosine hemisphere distribution, just as an example of a distribution that does not match the BRDF very well. Then they calculate the samples as $f(\omega_i, \omega_o) / p_{\cos}(\omega_o)$ where $p_{\cos}$ is the pdf for the cosine distribution. This is just to show you that it gives the "correct" result, but with a ton of variance. (Although it looks like there are some black areas on the left image that suggest it's not really converging to the correct answer...but let's assume that's a bug.)
Then the image on the right is sampled with (presumably) $\omega_m$ drawn from a distribution matching the $D(\omega_m)$ term from the BRDF. Then when you calculate the BRDF you can leave out the $D$ factor since it's canceled by the $D$ factor in the pdf. (The pdf probably still has some normalization and some conversion of probability from $\omega_m$ to $\omega_o$, so you still have to divide by those factors.)


What would be an example of decoding-intensive video?

You could try to play 1080p and 2160p videos, possibly 60fps in YouTube. If you right-click on the playback window and click on "Stats for nerds" you will get more information about used video and audio codecs and dropped frames, which can indicate that your hardware may be too slow. Currently the videos are mostly encoded as VP9 codec, but you can change settings to prefer AV1 codec, at https://www.youtube.com/account_playback. However please note that not all videos are available as AV1, but you can find some videos at AV1 Beta Launch Playlist.


Is there a objective reason for matrix naming conventions?

I think the naming order is intuitive because it is in reading order (left to right), e.g., worldViewProjection means that your point/direction is first multiplied by the world matrix, then the view matrix, and then the projection matrix. In this manner, you know the correct multiplication order by just reading the variable name and you do not have to think about mathematics at all. In your case, I think the problem is that your intuition seems to be focused on how you write the equations on paper or code.

The idea that a point is a column vector and matrices are multiplied on the left is also just a mathematical convention that you are used. You can use a row vector and multiply on the right by the transpose matrix, and the math is exactly the same, but would imply a reverse writing order in compound transformations. Some people just prefer one notation over the other.

I have read some old computer graphics books that use row vectors as points/directions, but it seems that most modern books have converged on using column vectors. IMHO, maybe this change during the years is one of the sources of confusion for so many new CG students.

For instance, if you choose to represent points as row vectors you would write something like this in your code: new_point = point * world * view * projection. On the other hand, if you use column vectors you would write this: new_point = projection * view * world * point. The important thing is that these matrices must not be the same, that is, the world matrix in the first code is the transpose of the world matrix in the second code, and the same is true for the others. As you can see, if you choose one math convention over the other, this will probably have influence on your naming conventions while writing your code.


What transformations should I apply to a model to produce this 3D effect?

Assuming your optimal viewing angle is parallel to the surface of the display and the pyramid is made from faces that are 45 degrees to its virtual (non-existant) base, it's actually just a simple non-transformed image (besides the reflection).

1:1 projection. No transformations. No scaling.


Barycentric coordinates for 5-sided polytope of two triangles and three trapezoids

I think I can do this using trilinear interpolation and treating two vertices on each side as merged and using my own volume calculations. I just need to interpolate some values given by the three points and a z. http://bmia.bmt.tue.nl/People/BRomeny/Courses/8C080/Interpolation.pdf


Rendering equation in terms of paths rather than directions

This is a very good question. There is a common misconception that Monte Carlo, or integration is applied "recursively" on the rendering equation. That is not what's happening. Numerical integration methods are tailored to problems of the form:

$$I = \int_{\Omega}{f(x)d\mu(x)} \approx \sum_{k=0}^{N-1}w(x_k)f(x_k)$$

Note that this is not the case for the rendering equation. Namely the fact that it is an equation, and the unknown radiance function $L$ shows up on both sides. In order to apply numerical integration techniques (such as Monte Carlo), one has to bring the rendering equation to the above form. Fortunately this is a well studied problem. The rendering equation is a specific case of a Fredholm integral equation of the second kind. If one can show that the scattering operator in that equation is a contraction $||T||<1$ (which holds if brdfs are energy conserving, there are cases where this is not required but that goes outside the scope of my answer), then we may perform a Neumann expansion (https://en.wikipedia.org/wiki/Liouville%E2%80%93Neumann_series).

The commonly used form of the rendering equation is the solid angle formulation (here $\sigma$ is the solid angle measure, $\sigma(\omega) = \sin\theta\,d\theta d\phi$):

$$L(x,\omega_o) = L_e(x,\omega_o) + \int_{\Omega}{f(\omega_o,x,\omega_i)L_i(x,\omega_i)\cos\theta_i\,d\sigma(\omega_i)}$$

In order to perform the expansion, I will rewrite this in the area formulation (here $\mu(x)$ is the Lebesgue area measure):

$$L(x_1 \rightarrow x_0) = L_e(x1 \rightarrow x_0) + \int_{M}{f(x_2 \rightarrow x_1 \rightarrow x_0)L(x_2 \rightarrow x_1)\cos\theta_{x_1}\frac{\cos\theta_{x_2}}{||x_2-x_1||^2}V(x_2,x_1)\,d\mu(x_2)}$$

Where $M$ is the set of all surface points in the scene, $V(x,y)$ is the visibility function which is $1$ if there is nothing between $x$ and $y$ and $0$ otherwise. If the normals of the surfaces at point $x_1$ and $x_2$ are respectively $N_{x_1}$ and $N_{x_2}$, then $\cos\theta_{x_1} = N_{x_1} \cdot \frac{x_2-x1}{||x_2-x_1||}$ and $\cos\theta_{x_2} = N_{x_2} \cdot \frac{x_1-x2}{||x_2-x_1||}$. As for the radiance, $L(x_2 \rightarrow x_1)$ gives the radiance arriving at $x_1$ from the direction of $x_2$. And the brdf notation relationship is: $f(x_2 \rightarrow x_1 \rightarrow x_0) = f(x_1 \rightarrow x_0, x_1, x_1 \rightarrow x_2)$ (that means $\omega_o = x_1 \rightarrow x_0$ and $\omega_i = x_1 \rightarrow x_2$).

I will rewrite the above for simplicity and conciseness as:

$$L(x_1 \rightarrow x_0) = L_e(x1 \rightarrow x_0) + L_r(x_1 \rightarrow x_0)$$

Let us now split the incoming radiance into direct illumination arriving at $x_1$ and indirect (at least one bounce) illumination. The direct illumination is obviously due to direct light source rays arriving at $x_1$:

$$L(x_1 \rightarrow x_0) = L_e(x1 \rightarrow x_0) + \int_{M}{f(x_2 \rightarrow x_1 \rightarrow x_0)(L_e(x_2 \rightarrow x_1) + L_r(x_2 \rightarrow x_1))\cos\theta_{x_1}\frac{\cos\theta_{x_2}}{||x_2-x_1||^2}V(x_2,x_1)\,d\mu(x_2)} =
L_e(x1 \rightarrow x_0) + \int_{M}{f(x_2 \rightarrow x_1 \rightarrow x_0)L_e(x_2 \rightarrow x_1)\cos\theta_{x_1}\frac{\cos\theta_{x_2}}{||x_2-x_1||^2}V(x_2,x_1)\,d\mu(x_2)} + $$

$$\int_{M}f(x_2 \rightarrow x_1 \rightarrow x_0)\Bigg[\int_{M}f(x_3 \rightarrow x_2 \rightarrow x_1)L(x_3 \rightarrow x_2)\cos\theta_{x_2}$$

$$\frac{\cos\theta_{x_3}}{||x_3-x_2||^2}V(x_3,x_2)\,d\mu(x_3)\Bigg]\cos\theta_{x_1}\frac{\cos\theta_{x_2}}{||x_2-x_1||^2}V(x_2,x_1)\,d\mu(x_2)$$

All I did here was separate the sum $L_e(x_2 \rightarrow x_1) + L_r(x_2 \rightarrow x_1)$ into two integrals, and then I additionally expanded $L_r$ into the corresponding integral using the recursive definition. You can expand this until infinity. For conciseness, I will rewrite the integration as an operator: $L_r = TL$. Now we can see the relationship to the Neumann expansion. Since the rendering equation can be written as:

$$L = L_e + TL$$

The solution is formally given as:

$$(I-T)L = L_e$$
$$L = (I-T)^{-1}L_e$$

Applying the Neumann expansion yields:

$$L = \sum_{i=0}^{\infty}T^iL_e$$

Note that each term of the sum is an increasingly dimensional integral. The first term is obviously just $L_e(x_1 \rightarrow x_0)$, the second term $TL_e$ is the integral for $L_e(x_2 \rightarrow x_2)$ that I wrote when I split $L_e$ and $L_r$, and so on. What each term gives you is the energy coming from a light source after $i$ bounces. The first term gives you the radiance emitted from $x_1$ towards $x_0$ (paths of length 0), the second term gives you the radiance that is due to direct illumination of $x_1$ scattered towards $x_0$ (paths of length 1), then you have the radiance due to lights two bounces away (paths of length 2) etc. With each bounce each new integral in the sum gains two dimensions (if we are just integrating over incoming directions). I would also like to emphasize that the integration variable $x_2$ in $TL_e$ is not the same $x_2$ as the one in $T^iL_e$ (this means that in the example with integration, $x_2$ in the first integral is not the same as in the second one - they are just integration variables). Another important fact is that within each sum each path starts at the film $x_0$ and ends at a light source ($L_e(x_2\rightarrow x_1)$ is non-zero only for points $x_2$ lying on light sources).

Now that we have a sum of integrals, we can apply our numerical techniques to each integral to estimate the sum. An obvious optimization, is reusing the samples used to compute $T^iL_e$ to compute $T^{i+1}L_e = T(T^iL_e)$. This yields a formulation where we can formally write out the sum of integrals as one integral with integration over all paths that start at the camera film and end at a light source. This also illustrates why Monte Carlo is a preferred technique: since we want to estimate an infinite dimensional integral, and Monte Carlo's convergence doesn't depend on the dimensionality unlike standard quadrature rules (there's also the point that it doesn't care about the smoothness of the integrand, but that's a double edged sword).

As you may have noticed, with the path integration formulation an exponential path count growth is not required. Splitting a path into multiple new ones is a technique known as splitting. While it is beneficial for high energy paths, on average with each bounce the energy gets lower due to attenuation, so this is in most cases counterproductive. On the other hand early path termination (for low energy paths) is more often than not beneficial in terms of efficiency (that's simply cutting off the sum at some point and not computing the remaining infinitely many terms, because you deem they have a contribution that is too low), and that's where Russian roulette comes in.

All of this is actually contained in Kajiya's paper on the rendering equation, even though he doesn't go into the details I went. He just refers to Rubinstein's book: https://dl.acm.org/citation.cfm?id=539488 (note - the first edition, later editions do not have the part he is referring to). In this book how to solve Fredholm integrals of the second kind is described (what I explained above, but more formally).

I hope my explanation was useful, and will continue being useful for future readers.


Looking at the video you provided, I would tend to think it's something along those lines:

A particle emitter, to which all threads are attached, is moving around. I haven't really figured what its motion is dictated by.
Each thread is made of particles flowing away from the emitter. The thread has a direction, used for the initial speed of the particles. Once created,  they are independent from the emitter.
The particles use a quad or a rectangle textured with a fuzzy ellipse shape (you can see them individually when a thread suddenly changes direction) oriented in the direction of the motion. Their color is based on the duration since their creation, and they use additive alpha blending (the saturation to white is typical). After a certain duration, they fade out and die, only to be recycled into new particles.
I suspect the whole thing might be in a vector field, or something similar blowing wind onto the particles, which would explain why the threads sometimes suddenly change direction.
There might be some feedback post-processing effect to smooth everything out.


Can somebody explain this Ray Tracing Function?

Good answer above.    (I like that the book adds in p then subtracts it.)

This is a hard-coded grey world as a "hello world" kind of ray tracer.   The 0.5 is the hard-coded reflectance.   Keep working through book and the object reflectance and max-recursion will be addressed.   Ray branching I avoid like the plague.   In cost benefit I am always a skeptic on ray trees versus ray paths.


Is modifying a texture (painting on it) considered a "state change"?

Updating an area of memory in the graphics device (a texture, buffer, and the like) is not quite the same as changing a rendering state.

What makes a render state change expensive is the amount of work the driver has to do to validate the new state(s) and reorder the pipeline. This will most likely also incur some synchronization between CPU and graphics device. However, the amount of data transferred between the devices should be small for a state change (probably just a few commands).

For a texture/buffer update on the other hand, the main cost lies in the data transfer itself. In theory, unless you're reading the texture data back to the CPU after the update, there should be no synchronization or pipeline stalls. However, another aspect should be considered: API overhead. Even if the amount of data you're sending to the graphics device is small, if you do that often enough, eventually the cost of communicating with the driver/device will become greater then the cost of the data transfer. That's another reason why batching is so important when optimizing a renderer.

So in your case, the best approach, it seems to me, would be to keep a system-memory copy of the texture that you update whenever new data arrives. Set a dirty flag and consolidate as much updates as possible into one glTexSubImage for the whole texture (or a large sequential portion of it). You can also play with Pixel Buffer Objects and try to do asynchronous data transfer to reduce pipeline stalls as much as possible. If you can implement some kind of double buffering, then you can write to one copy of the texture while the other is being rendered with. This tutorial explores that scenario. That's my intuitive approach, I'd try to reduce the number of API calls and "batch" the texture updates. That being said, this is very speculative, and you'd have to profile and compare that to other approaches, like doing several small updates, to know for sure which is the most performant in your usage case.

As a side note, this presentation by NVidia is also relevant and provides a lot of good insights: Approaching Zero Driver Overhead in OpenGL.


BRDF that interpolates between diffuse, glossy and reflective?

Any microfacet BRDF with a roughness parameter will do something like this—for example using the GGX or Beckmann NDFs. When the roughness goes to zero it becomes a mirror; as roughness increases the reflection will be wider, and in the high-roughness limit should be essentially Lambert I think. Roughness is not necessarily limited to [0, 1] as roughness values higher than 1 are possible, but typically roughness = 1 will be very close to Lambert already.
One thing to note about the roughness = 0 limit is that the NDF will become singular (a delta function), so will very likely break the shader code by producing infinite values or NaNs and such. So as a practical matter it might be necessary to special-case that, or clamp the roughness to some low minimum value like 0.0001.


How are normals assigned to vertices in different types of shading?

Flat shading uses the same normal for each vertex of the triangle. In this case the normal is constant, and one color is evaluated for the triangle.

For Gouraud and Phong shading normals are different per vertex and can represent a curved surface (even though the face is flat).

However, for Gouraud shading is evaluated per vertex, and the resulting color is interpolated. For Phong, the normals are interpolated first and the shading is evaluated per pixel.

Also note there's a bit of an overloaded term between the Phong Lighting/Reflection model and Phong Interpolation, Per-pixel shading which is what you're referring to here.

Given these 3 shading types, the next logical step is actually normal maps shading. This is basically where the normal is provided per pixel, and shading is evaluated per pixel.

And so you have the trend of doing more calculation at finer details / higher resolution.

Here's a table!

                 Normal         |   Lighting
Flat:            Per Triangle   |   Per Triangle
Gouraud:         Per Vertex     |   Per Vertex
Phong:           Per Vertex     |   Per Pixel
Normal Mapping:  Per Pixel      |   Per Pixel


What it means to "triangulate complex 3D objects" is not unambiguous.
Just one possible interpretation: You have a 3D polygon in space, and you
want to triangulate that. This is NP-hard:

  Barequet, Gill, Matthew Dickerson, and David Eppstein. "On triangulating three-dimensional polygons." Proceedings of the 12th Symposium on Computational Geometry. ACM, 1996. (ACM link.)


How is an ICC profile embedded or built into an image?

I realized that littlecms, suggested by joojaa in a comment above is very complete in dealing with the ICC profiles.

I did a few tests in their color translator program with different destination color spaces and the differences are great as you can see:

Original

sRGB-like

adobeRGB1998 (it seems it's the original color space)

appleRGB

And the most desirable results in my case are within the ProPhotoRGB color space:

linearProPhotoRGB

All of the above were following an automatic rendering intent but other intents bring slight variations (samples applied together with linearProPhotoRGB destination color space):

Absolute Colorimetric

Legacy Unadapted Absolute Colorimetric

I also realized that GIMP was actually restoring the original profile (not the opposite as I thought, i.e. instead of applying an embedded ICC profile, it was taking it off) and with this littlecms, now it does seem much simpler to "translate" color spaces than I thought.


OpenGL pipeline stages, when is the primitive assembly executed?

For a diagram to be comprehensive, it'd have to show way too much detail. It's just a simple tool.
Primitive assembly is the process of taking a sequence of vertices for a particular rendering primitive type, and generating a sequence of base primitives (points, lines, triangles) from them.
This process does indeed happen after vertex post-processing, as the input to the rasterizer is a sequence of base primitives. After all, you don't rasterize a triangle strip; you rasterize a triangle.
However, the input to a geometry shader invocation is a single base primitive, taken from the sequence of vertices provided up to that point. As such, a form of primitive assembly also must happen before the GS, to convert the vertex sequence into a sequence of base primitives.
This is also one reason why GS's are slow. If there is no GS, then this primitive assembly doesn't happen.
The output of the geometry shader is a bunch of vertex sequences, which still later undergo primitive assembly based on the output primitive type of the GS. A GS can output triangle strips, so there has to be primitive assembly to turn those strips into individual triangles.


How to store > 256 transparency values in textures?

A normal argb image has 8 bits per channel and therefor 32 bits per pixel, the size of an integer. That's why it's stored in an integer a lot of the times.

Floats are also 32 bits.

Something that you can do, is to have a separate image that has only 1 channel, for the transparency. If you make that a float type (or for 16 bit, another type that uses 16 bits instead of 32 bits), then it will take about the same room as an extra normal image.

After a quick google, I found out that most programs are limited to 8 bit per channel, but there should be programs that have 16 bit per channel.
This gives you at least more precision.

When loading in the normal 16 bit image (your texture, if you use 16 bits per channel), you could convert the normal rgb part to 8 bits, and load that in and then create a new image that has 1 channel that is 16 or 32 bits, and then store the alpha value in there.

This way it should be similar to having an extra image, and you do, but instead of 2 textures, you have the transparency in another image.

What would then also be possible to do, is to store the transparency in a normal rgba image (inside your program as a texture), as you were talking about, and then in the shader, you manipulate the rgba values to restore the transparency value. This way you wouldn't need a floating point image.

This is just a quick idea. As far as I know, this should work. (This is my first answer here! whoop! whoop!)

Hope it helps and good luck!

https://developer.mbed.org/handbook/C-Data-Types

(^for the data sizes)


Unity: Pass texture pointer from OpenGL to Unity

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

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

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

Best.


Shader code for putting an atmosphere over a planet

I believe that Space Engine as a base algorithm still uses Precomputed Atmospheric Scattering by E. Bruneton. It is released on github. So it is good starting point if you want to implement atmospheric scattering.

For beginner this algorithm is pretty advanced stuff written in C++11, and modern OpenGL, so such knowledge is required. In short, to avoid a cost of solving rendering equation for scattering in participating medium (single as well as multiple scattering) it uses precomputed lookup textures where it is possible and equations for Rayleigh and Mie scattering based on physical model.

There is also simpler demo computed in real time on shadertoy. This is probably better place to start. It involves raycasting a sphere and integrating Rayleigh and Mie scattering along a ray in real time. It is based on GPU Gems article, which provides some good information how it works.


Is there a term for this sort of graphics?

If the terrain is indeed represented by a 3D mesh as it looks like in the images, then I'm not aware of any specific terminology (and I dont believe there might be one).
I could say it was similar to a faceted shading model, but shading looks very rudimentary and lighting computations might not exist at all.
In case you intend to replicate these looks, I'm afraid you already hinted how to proceed very low-polygon non-textured meshes and no lighting.