Barycentric Coordinates of a Point on a Line

If you have a point $C$ on a line defined by a point $A$ and a point $B$, here is how you calculate the barycentric coordinates $(t_0,t_1)$:

$t_0 = \frac{sqrt((C_x-A_x)^2+(C_y-A_y)^2)}{sqrt((B_x-A_x)^2+(B_y-A_y)^2)}$

$t_1 = 1.0 - t_0$

In other words, you calculate the length of the vector from $A$ to $C$, and divide that by the length of the vector from $A$ to $B$.  This gives you $t_0$, and you can subtract that from 1.0 to get $t_1$ since barycentric coordinates always add up to 1.0.

A more mathish style version of the above is this:

$t_0 = \frac{\left\| C-A\right\|}{\left\| B-A\right\|}$

$t_1 = 1.0 - t_0$

But it's important to note that $C$ must be on the line defined by $A$ and $B$.  If that isn't the case, let us know and let us know what it is you are trying to do specifically.  For instance, maybe you want to project some other point $D$ onto the line $\overline{AB}$ and then find the barycentric coordinates of that projected point?


Ray tracing - partial derivatives for texture lookup

Yes, ray differentials are the way to go. The Paper by Igehy introduces them for the use case of filtered texture lookups. When generating the primary rays, you initialise the differentials to reflect the pixel footprint. As the ray progresses through the scene, you update the differentials at every bounce. When it comes to a texture lookup, you need to project the differentials to your s,t texture space (aka UV), typically by going to surface parameter space as an intermediate.

The case that you mention, the calculation of derivatives for a ray surface intersection, is implemented in pbrt:
Here the surface is described by a point p, normal n and two tangent vectors dPdu and dPdv. The resulting dudx, dudy, dvdx and dvdy are the differentials of the u, v parameters of the surface with respect to screen space x and y. Conversion from surface parameters to texture space happens here.


Drawing Directed Lines (with arrow endings)

I don't see anything wrong with your code except for possibly the normalization of the direction. If the length of the vector gets updated automatically, it will change after the division of the x coordinate by the length. What you probably want is

public void Normalize()
{
    float len = this.Length;
    m_x /= len;
    m_y /= len;
}

instead. This would also explain why the arrow endings look more off for horizontal directions than for vertical directions. If the x component of the direction is small, then a division of m_x by the length only leads to a small change in the length of the direction. (Try to change the order of the two divisions in your Normalize function; then it should be the other way around.)


How can I concentrate points in areas of higher curvature?

The idea I would try to apply would be the following: I make the example for the curve, but it should be straightforward for the application for the surface.

Let's say we have a curve $\gamma$ uniformly parametrized. Let's say the parameter of the curve is $s$. Your goal is to sample point corresponding to value of $s$ such that the curvature is high.

If you get the magnitude of the curvature $c$, this will be function of $s$ too. So, if you normalize the function $|c|$, you will get a probability distribution. If you get the integral of such distribution, you will have the cumulative distribution. Let's call this cumulative function $C(s)$.

The sampling problem from a distribution given by the cumulative function is well known, so basically once you have sampled a set of value $s_0,s_1,\dots,s_n$, such value will be related to the points of interest.

The application of this method to the case of surface should be straight, since basically you have a two dimensional cumulative distribution function, but the sampling problem is exactly the same.

Just to give some detail, it is basically sampling from a distribution given the cumulative function involves two steps:

take a random value in the interval $[0,1]$, let's say $k$
solve the equation $C(s) = k$.

This approach is exact, of course it is expensive, but if you like such approach you can work on optimization.


Can I determine the pixel geometry programmatically?

It appears that Microsoft has punted on this in Windows 7:

This is the method available in the control panel for selecting what layout ClearType uses.

Additionally, it seems that iOS and the Windows modern UI style de-emphasize subpixel antialiasing heavily, due to the prevalence of animations and screen rotations. As a result I expect the OS vendors to not spend a lot of effort trying to figure out the subpixel layout of every screen.


How to raytrace triangular faces when camera is not at origin?

Wow I'm stupid, the whole problem was I thought it was 1/t and it is actually 1-t in parametric form. Hope this helps someone anyway.


How can I utilize all cores of CPU In my ray tracer with std::thread?

This is not really a CG issue but you appear to have an integer division problem. Assuming mMaxThreads is an integer,  (threadIndex / mMaxThreads) is always 0 and ((threadIndex + 1) / mMaxThreads) is 0 for all threads but the last one, where it is 1. So one thread is doing all the work.


Special Texture Stretching without deformation

The Nine Patches algorithm is a really nice and powerful solution if you are using textures but I ended up by using a fully procedural solution.

I am just drawing a round box like this.

float udRoundBox( vec3 p, vec3 b, float r )
{
  return length(max(abs(p)-b,0.0))-r;
}

You can simply avoid the 3rd dimension if you don't need it.
This method comes from the Distance function article from Inigo Quilez.


Just to add to whychmaster's reply, in my experience, on a GPU, a floating-point operation, when used in quoting FLOPS (floating-point operations per second) benchmark figures will refer to Addition (equiv subtraction) and Multiplication.
The other standard IEEE operators, division + sqrt, are more computationally expensive and so unlikely to be included in the figure. (Note that, there is often HW support in GPUs for accelerating operations like 1/x, 1/sqrt(x) etc, perhaps to make them as (or nearly as) fast as muls & adds)
Having said that, it is very likely that a GPU will actually implement a Fused Multipy-Add operation, FMAD, i.e. do  A * B + C in "one step" as this is

more accurate than doing separate multiply and add instructions
is faster than separate ops
is about the same silicon cost
is a very common building block of CG operations (e.g dot products, convolutions)

I suspect that in quoted figures, though, an FMAD would be considered 2 operations.


Techniques for making tesselated shapes appear smoother?

GL_SMOOTH interpolates the vertex attributes of the primitive (triangle), it looks like your normals are setup per face, rather than per vertex.

If using normal per face (i.e. a single normal for the triangle), increasing the number of triangles will produce a less faceted look as there is more granularity when interpolating the normals over the surface of a given surface area given the increase in the number of triangles within that area. It is recommended to use normals per vertex for smooth surfaces when feature edges are not important (i.e. Not for cubes!).

If using normal per vertex (i.e. a normal for each vertex, resulting in each triangle having 3 vertices), each 3 vertices of a face will be slightly different if the neighbouring triangles are not coplanar to the triangle in question (i.e. not a flat surface). This will aid in the smoothness of the shading.

Normals per vertex can be easily calculated by averaging the three normals of the three vertices defining a triangle, however your surface needs to be defined with shared vertices for this. The exception to this is if you have a sphere or ellipsoid, in which case you can calculate a more precise normal easily by unitising vertex - sphereOrigin rather than averaging the three normals of the face.

On a side note: Are you using two lights in the scene of the screen shot you posted or one? If two it looks like the normals are currently not correct as the shading does not appear correct since there is a well-lit area of the surface on the bottom left region. With a single light, either the front (top right) of the sphere would be lit, or the bottom left, but not both.

Also note that simply increasing the power of the light simply over exposes the surface colour rather than improving the smoothing. You lose the surface definition.


Does prerendering actually help performance?

As far as I know, this sort of thing is mainly about shader compilation. One of the main reasons why a game may experience hitches the first time something renders is that the shaders necessary to render it haven't finished compiling yet, and the driver has to finish that work before the frame can proceed.

A little bit of background. When you write shaders in HLSL, you compile them with Microsoft's HLSL compiler (fxc.exe), which outputs a hardware-independent bytecode. When you call CreateVertexShader etc. in D3D, that bytecode gets passed to the underlying driver, which ultimately compiles it into the machine code for the specific GPU you're running on.

This internal compilation step isn't trivial—it's a full-blown optimizing compiler. So it can take a little while on a complicated shader. Because of this, when you call D3D's CreateVertexShader, the driver doesn't necessarily finish compiling the shader right then and there; it may add it to a queue, then return control to the application. Meanwhile, a background thread (perhaps more than one) consumes the queue and compiles the shaders. (The driver may also have a disk cache of previously-compiled shaders, so that it doesn't have to redo all this work every time the game is launched.)

This background-thread compilation means the main game thread can get on with other work, e.g. loading other data or building a scene graph or what-have-you, so it's generally a win for overall loading times. But the consequence of it is that the shaders may not all be finished compiling when the game starts rendering. Unfortunately, there is no way with DX11 and earlier APIs for the app to tell when the shaders are done. The best you can do is to render something with the shader—this effectively makes the app wait for the compiler to finish. That's why games sometimes do this offscreen pre-rendering business.

I don't have any hard numbers handy, but qualitatively, one can reduce the impact of shader compilation by front-loading it during loading. Create all your shaders first, then the compiler will run through those while you load everything else. However, with real-time streaming  (e.g. open-world games) it can still be quite difficult to ensure no hitches.

It's also worth noting that there are some cases where state changes can cause a shader recompile. Things like vertex buffer formats, render target formats, and shader linking details (e.g. re-using the same pixel shader with different vertex shaders, or with tessellation shaders) can potentially cause the driver to have to generate different versions of the shader under the hood. The details vary by vendor and hardware. This isn't necessarily a full from-scratch recompile—often it's just sticking on some prologue or epilogue code—but it's something to be aware of.

This is one place where the new APIs (DX12, Metal, and Vulkan) are in a good position to solve the problem. These APIs have "pipeline state" objects where all shader stages and all possibly relevant states are rolled up into a single object, and when you create one of these, the driver has all the information it needs to actually and completely finalize all the shaders. DX12 also makes shader caching explicit. Also, at least for NVIDIA DX12 drivers (I don't know about others), shader compilation is done fully inside the pipeline state creation method, and the app can use multiple threads itself in order to parallelize shader compilation and other work. So, the practice of offscreen pre-rendering to ensure shaders are finalized should no longer be necessary in DX12.


Is glDispatchCompute(a,b,c) any faster than glDispatchCompute(a*b*c,1,1)?

Broadly speaking, no, these won't really matter.
In specific cases, it could, but only because it might influence the order in which work groups are processed. And that would mean that it might influence cache behavior and the like, which could matter if neighboring work groups access related memories.
But even in such a case, which one is better would depend entirely on the nature of the shader's memory access patterns and the driver's decisions about ordering work groups. One shader might run slightly faster with a linear ordering, while a different shader might run faster with a multi-dimensional ordering. It all depends on what the shader is doing (and how the driver handles it), so you can't say that one is inherently better than the other.


When is a compute shader more efficient than a pixel shader for image filtering?

An architectural advantage of compute shaders for image processing is that they skip the ROP step. It's very likely that writes from pixel shaders go through all the regular blending hardware even if you don't use it. Generally speaking compute shaders go through a different (and often more direct) path to memory, so you may avoid a bottleneck that you would otherwise have. I've heard of fairly sizable performance wins attributed to this.

An architectural disadvantage of compute shaders is that the GPU no longer knows which work items retire to which pixels. If you are using the pixel shading pipeline, the GPU has the opportunity to pack work into a warp/wavefront that write to an area of the render target which is contiguous in memory (which may be Z-order tiled or something like that for performance reasons). If you are using a compute pipeline, the GPU may no longer kick work in optimal batches, leading to more bandwidth use.

You may be able to turn that altered warp/wavefront packing into an advantage again, though, if you know that your particular operation has a substructure that you can exploit by packing related work into the same thread group. Like you said, you could in theory give the sampling hardware a break by sampling one value per lane and putting the result in groupshared memory for other lanes to access without sampling. Whether this is a win depends on how expensive your groupshared memory is: if it's cheaper than the lowest-level texture cache, then this may be a win, but there's no guarantee of that. GPUs already deal pretty well with highly local texture fetches (by necessity).

If you have an intermediate stages in the operation where you want to share results, it may make more sense to use groupshared memory (since you can't fall back on the texture sampling hardware without having actually written out your intermediate result to memory). Unfortunately you also can't depend on having results from any other thread group, so the second stage would have to limit itself to only what is available in the same tile. I think the canonical example here is computing the average luminance of the screen for auto-exposure. I could also imagine combining texture upsampling with some other operation (since upsampling, unlike downsampling and blurs, doesn't depend on any values outside a given tile).


how to build the feature curves of this algorithm?

The way I read it, you have a 2D cubic Bézier that defines the projections of the ridge line (or river bed, or whatever) onto the X/Y plane. In Figure 4, these are the dark purple lines.

So that gives you the direction that the ridge line moves. Now you have to define the height of the ridge line. This is done by picking some points along the 2D Bézier and setting the following values at each of the chosen points:

h = height
r = radius (this appears more like a width at the ridge line than a curved radius)

You can then also set the angle of the terrain as it comes off the radius. These are:

`phi` and `theta`
`a` and `b` are the lengths of the slopes

Keep in mind that you don't have to set every value for every control point you pick. I assume any you don't set are defaulted to 0, or some reasonable default.

So for a hill, you would set h to the height of the hill at a few point. You would then set phi and theta to be horizontal, so yes, pi / 4. The length of a and b will decide how quickly the hill falls off as you move away from the ridge line.

It looks like the above 6 values are turned into another cubic Bézier curve that is perpendicular to the first one for the rasterization step in section 4 of the paper. h doesn't factor into the shape of the perpendicular Bézier because it's also perpendicular to it. So r, theta, phi, a and b are used to make a cubic Bézier. I believe that you can think about like this: the point on the ridge line where you define one of these control points is the link between 2 cubic Bézier segments - one to the left of the point and the other to the right. The images they show only show a single point between 2 Bézier curves, but the way I'm reading it there are actually 2 curves that look something like this:

So I think normally, you'd have 2 Béziers. One on each side of the control point they show in the paper. (Pardon my bad drawing.)  a would be the distance between control points 2 and 3 of the left curve. Theta would be the angle formed between the vertical and a of the left curve. r would be the distance between points 3 and 4 of the left curve. I'm unclear on how one would pick control point 1 in this scenario, though.

Ai and Ri have to do with the noise parameters and don't appear to me to be related to the Bézier curves at all.


UV partial derivatives of a cylinder shape primitive

I think it's easiest to get the tangent frame by writing the forward mapping from $(u,v)$ to $\mathbf{p}$:
$$
    \mathbf{p}(u,v) = \begin{bmatrix} R \, \cos (2\pi u) \\ R \, \sin (2\pi u) \\ vL \end{bmatrix}
$$
Then you can see that the derivatives are:
$$
    \frac{\partial\mathbf{p}}{\partial u} = \begin{bmatrix} -2\pi R \, \sin(2\pi u) \\ 2\pi R \, \cos(2\pi u) \\ 0\end{bmatrix} = \begin{bmatrix} -2\pi y \\ 2\pi x \\ 0\end{bmatrix}, \qquad
    \frac{\partial\mathbf{p}}{\partial v} = \begin{bmatrix} 0 \\ 0 \\ L \end{bmatrix}
$$
which matches the pbrt code.
I think your derivation went wrong by assuming that $\partial\mathbf{p}/\partial u$ was the component-wise reciprocal of $\partial u / \partial{\mathbf{p}}$. It's not so simple, because the latter object is really the gradient of $u$, where $u$ is now considered as a scalar field defined over all space (not just on the cylinder). To invert this properly we really need another coordinate running perpendicular to $u$ in the $xy$ plane, so that we have a locally invertible mapping; then we can take its Jacobian matrix and invert that. We can use the radius $R$ as such a coordinate. Then we have the inverse mapping,
$$
\begin{aligned}
    u(x,y) &= \frac{1}{2\pi} \arctan(y,x) \\
    R(x,y) &= \sqrt{x^2 + y^2}
\end{aligned}
$$
whose Jacobian is:
$$
    \frac{\partial(u,R)}{\partial(x,y)} =
    \begin{bmatrix}
        \tfrac{\partial u}{\partial x} & \tfrac{\partial u}{\partial y} \\
        \tfrac{\partial R}{\partial x} & \tfrac{\partial R}{\partial y}
    \end{bmatrix} =
    \begin{bmatrix}
        -y/(2\pi R^2) & x/(2\pi R^2) \\
        x/R & y/R
    \end{bmatrix}
$$
and the inverse of that is:
$$
    \frac{\partial(x,y)}{\partial(u,R)} =
    -2\pi R \begin{bmatrix}
        y/R & -x/(2\pi R^2) \\
        -x/R & -y/(2\pi R^2)
    \end{bmatrix} =
    \begin{bmatrix}
        -2\pi y & x/R \\
        2\pi x & y/R
    \end{bmatrix}
$$
from which you can then read off $\partial{(x,y)}/\partial u$ as the first column, which matches the result obtained for $\partial\mathbf{p}/\partial u$ above.


Ray Tracing Bounding Volume Hierarchies nodes

Bounding Volumes
A bounding volume is for example an Axis-Aligned Bounding-Box (AABB), a generaal Bounding-Box not aligned with the major axes or a Bounding-Sphere. They differ in

how much memory they take (e.g. a sphere takes only 4 floats, AABB 6 floats)
how good they enclose the object's geometry
how fast are the related computations (mainly intersection tests).

You will find an extensive overview at Wikipedia.
Hierarchy of volumes
So in the hierarchy of bounding volumes, instead of whole meshes you will store only their bounding volumes. You choose a single type, e.g. a bounding sphere, and will use it to build the whole BVH. At the bottom level each bounding volume will probably enclose only a single object. You will build the hierarchy up by enclosing pairs or smaller groups of volumes in a larger volume.
So using spheres you would enclose each object into a bounding sphere at first. Then at the higher levels of the hierarchy, you would always create a new sphere enclosing a few spheres from the previous level which are close to each other. You proceed until a single root sphere encloses all 3D objects in your scene.
There are different strategies to create BVH. But all have in common that complex structures like triangle meshes are hierarchically clustered using very simple entities.


You can do this by writing a renderer plugin for 3ds Max. This gives you the same kind of access as the built-in renderers, including a settings dialog, a preview window, and an interface where you can write the final frame of a full render (and 3ds Max will write it to the output file). The renderer plugin gets full access to the scene graph when the user hits "Render" (or when the preview window is in use). Any complex geometry (such as subdivision surfaces) will already have been converted to a triangle mesh for you.

PaulHK already linked the official Programmers' Guide online. You'll also find it in the help files that are installed with 3ds Max, along with sample source code for a trivial renderer plugin. Because your plugin needs to link against the 3ds Max SDK, you'll need to compile it with the same version of Microsoft Visual Studio that was used to compile the particular version of 3ds Max you're targeting. (This means that if you want to target multiple versions of 3ds Max, you might need multiple builds.) Again, how to do this is described in the Programmers' Guide.


GLSL shapes signed distance field implementation explanation?

Torus

A torus is defined by two parameters: the major radius, and the minor radius. The major radius (t.x) is the radius of the big ring (in red in the diagram), and the minor radius (t.y) is the radius of the circular cross-section. The x and y here are just indices into the vector: they're unrelated to the x and y axes.

The left side of the diagram shows the xz plane, and the right side shows a cross-section with y going horizontally.

p.xz is the projection of the point onto the plane of the torus, so length(p.xz)-t.x gives the distance from that projected point to the inside of the torus. (I've labelled it $r$ on the diagram.) p.y is the distance from the point to the plane of the torus, so the length of the vector ($r$, p.y) is the distance from the point to the inside of the torus. (Computing the length of a two-element vector is a handy way to use Pythagoras' theorem to find the hypotenuse of the triangle.)

Cone

Like the torus, the cone code uses the fact that the shape is symmetric about an axis: the z axis in this case. The cone is defined by a unit vector c that's normal to the surface of the cone, and the apex of the cone is at the origin. Similar to the above, length(p.xy) gets the distance from the point to the z axis, so that we can work in the plane where the cone is a triangle. vec2(q,p.z) is the projection of p into this plane. The dot product projects this vector onto the direction c: how far the point is along this direction tells us how far it is from the surface of the triangle.


Bicubic bezier surface from 4 bezier curves

Yes, its pretty standard stuff. Most 3D applications our there can do this. Since the nomenclature of 3d applications is not standardized, the tool has different names in different applications: Loft (Maya, 3DSMax), Blend (Creo, Catia...) and so on.

Image 1: Blend between 3 curves, not closed but same principle applies. Open curves are easier as they clearly define the boundary conditions. (software used is Creo because that's what I had open)
You usually want the first edge be from a surface edge because otherwise you will lose the tangency info. But nothing prevents you from using a natural tangent direction, you just lose control. Also you need to pay attention to how you intend the curves to connect (what point connects to what point).
Personally I prefer drawing the hull or using some of the network based tools. This way allow you to specify control curves in 2 directions, thus giving control of what to connect where and how it should work mid run.

Image 2: Curve network.
There are many ways to do this, infinite ways in fact. Also if you want to be compatible with other mainstream CAD apps use NURBS instead of beziers than its easier for you to write a exporter/importer if your data set is rich enough. Beziers are just special cases of NURBS after all so your code itself should not need to change all that much (fitting also gets easier).
What you do is simply fit a curve trough your points in the second direction this will then give you the corresponding hull in the other direction. In your case you could just make a Bezier patch for each segment and then try to deduce a natural tangent direction.
One way
Like I said theres a infinite ways of constructing these, you have an artistic freedom of doing whatever you see fit it's a over constrained problem. One way of constructing this would be to have a patch for each consecutive cross section. This is nice and tidy as it generalizes works well with an arbitrary number of cross sections.

Image 3: Example of deriving one row for 2 consecutive bezier spans.
First recognize there are 3 cases (ill cover 2):

Free edges,
Internal edges,
External constraint.

I'll cover 1 and 2. First the Internal edge as its simple you want it to be smooth so the tangent length and direction in both patches must match. If that's not true then its not sufficiently continuous. You could simply take that this is the furthest point from the surrounding patches control points and due to mean value theorem it should have a tangent direction perpendicular to these (Line A-C is parallel with D-E). The tangent length should be about 1/3 of the length of the distance to the adjoining vertices (A-B or A-C, in this case the shorter of the spans. It could be longer or average. though that could pose problems). Obviously if you have more spans just repeat for one span forward and after for more internal spans if needed.
For free edges find the point on the perpendicular line that is as far from the point as half distance between the points (length B-D is equal to length (A-B)/2 and length B-E is equal to length (B-D)/2). The free tangent points towards this point. The length of the tangent is as long as one third of corresponding point difference (length (A-B)/3 and (B-C)/3). Obviously the distances and points could be arbitrary but these are somewhat natural.
Repeat for next point sequence in line and next and next. You end up with 2 spans that means 2 times 16 points.  There are other ways but this is relatively easy to convert into a visual form without going into deeper maths. But this is pretty well known math you can find many ways to form these by searching in literature.


What is wrong with my metallic PBR shader?

What is wrong ? Almost everything. It is not a PBR material at all. PBR metallic materials do not use ambient term and do use generally microfacet BRDFs and ComplexFresnel.


Why are texture coordinates often called UVs?

This is not a definitive answer, but it is generally accepted that Ed Catmull introduced Texture Mapping in his 1974 thesis, "A SUBDIVISION ALGORITHM FOR COMPUTER DISPLAY OF CURVED SURFACES"
In that, he uses (U,V) to access the image data (see the page labeled 36 in the above)

MAPPING
Photographs, drawings, or any picture can be mapped onto bivariate patches. This
is one of the most interesting consequences of the patch splitting algorithm. It
gives a method for putting texture, drawings, or photographs onto surfaces....
...If a photograph is scanned in at a resolution of x times y then every
element can be referenced by u·x and v·y where 0<=u,v<=1. In general, one could think of the intensity as a function l(u,v) where I references a picture.

I believe this thesis also introduced the concept of the Z-Buffer (Page 32)


Does a gaussian blur based bloom/glow shader require gamma correction? Any other tips?

In short:

You should not gamma correct your glow map. In fact, you should do everything in a linear color space. At the very end, when doing any color grading (which is the very last step), you convert the final image to the right color space. The color space include that gamma correction.

In long:

Gamma correction is a step in the encoding of images. It is usually done so that we use more bits to store the shadows and less bits to store the highlights, since we see less detail in the highlights anyway. It is also part of some color spaces. These color spaces make sure the image looks the way it is supposed to. It is the final step, often part of color grading. Just before it is send out to the display, or saved to a file. This means that everything before should be in a linear color space.

Usually things would look like this:

Load textures and convert them to a linear space, so you do not double gamma correct..
Render the image.
Do the post effects.
Convert the image to the final color space with color grading. Probably sRGB Gamma. Converting to the color space includes the gamma correcting.

Color grading is not necessary, but optional.

When having everything at linear and doing the gamma correction at the very end, you will need a high precision for the color. Using floating point framebuffers is recommended. When using the default 8-bit framebuffer the whole way, you end up getting banding effects when gamma correcting.

This setup is not always done for games, to my knowledge, since most people will barely notice it or even care about it. In the film industry, VFX and animation, using the proper color space is more important.

When not doing it the 'right' way, the textures are not converted to linear and at the end, there is no gamma correction. Lights are supposed to have a quadratic falloff to be physically correct, but when gamma correcting it with a gamma of 2.0, it becomes a linear falloff. Lights can therefore get a linear falloff and look normal. As Noah Witherspoon pointed out, when blending and blurring you should use a linear color space, however this is not always done either.

The reasons for not doing it the 'right' way, is probably that it still looks fine. You do not need to do it this way, to end up with good graphics. It is always a matter of what the artist likes the most.


Do I need separate Graphics and Transfer commands?

Your question seems... very confused. The "challenge" you tried to do specifically stated to "use a different queue family specifically for transfer operations". But that's not what your code is trying to do; you used the same queue family.

The point of creating a "transfer queue" is to allow you to perform transfer operations without having to submit them as part of another batch. vkQueueSubmit cannot be called on the same queue from multiple threads. So if one thread builds some transfer work, but the other thread is doing some graphics work, if you only have one queue then the two threads must have some communication mechanism to submit that work.

If you have a transfer queue, that inter-thread communication is no longer necessary. Or at least, it is no longer necessary for the basic submission action. Obviously whichever batch of work that uses the results of the transfer commands will have to receive a semaphore to wait on from the thread that submitted the transfer work. But that is usually a less contentious inter-thread communication path, which can easily be done with a lockless queue of some kind.

The challenge in question is specifically talking about a queue family that only has the VK_QUEUE_TRANSFER_BIT flag set (or more specifically, does not have graphics or compute capabilities). Such transfer-only queues are usually tied to specialized hardware that (potentially) perform transfers faster than a graphics queue, but likely with extra limitations (as defined in VkQueueFamilyProperties::minImageTransferGranularity and other properties).

So basically, you're doing the challenge wrong.

Do keep in mind that Vulkan does not require implementations to provide a transfer-only queue. So if you want to take advantage of the presence of such a queue, you also have to take into account the possibility that no such queue exists.

Now, here's the thing: the term "queue family" is used because you can potentially get multiple queues from the same family. VkQueueFamilyProperties::queueCount tells you how many queues a particular family allows you to get.

So if a GPU does not provide a dedicated transfer-only queue, but you still want the asynchronous benefits of submitting transfer operations outside of non-transfer work, you can create a second graphics queue.

Of course, Vulkan also doesn't require implementations to provide more than one graphics or compute queue. So again, you have to have a synchronous fallback plan in case the hardware only provides you with one queue.


If you have the prerendered images, you can "simply" render them onto a grid of quads. This is done with texturing, i.e. you upload the image to your GPU as a texture and apply texturecoordinates to quads, such that when you render the quads, the textures are applied correctly to their surfaces.

You will however need to learn OpenGL (preferrably with GLSL), Direct3D with HLSL or even Vulkan, all of which will take some time. Depneding on where your goals lie, you may be better of going with a preexisting rendering engine and just learning like Unity or Unreal.

If you do want to learn one of the mentioned Graphics APIs, there are numerous tutorials on the web. For OpenGL I can recommend the youtube tutorials by ThinMatrix, for the other two APIs you'd need to ask someone else. Furthermore, there are splendid books online as well as offline. If you have some money to burn, the basics of rendering are covered in Real-Time Rendering and Mathematics for 3d Game Programming and Computer Graphics. In regards to OpenGL, there are the Online Tutorials of Learn OpenGL, as well as offline books in pdf and epub version.


Jim Blinn's book Jim Blinn's Corner: Notation, Notation, Notation has a couple of chapters which go through this in detail with all of the edge cases. This book is a collection of essays from his column in IEEE Computer Graphics and Applications. Highly recommended.


Why would a 4d texture be used for color conversion?

I think your quote refers to conversion between color spaces, rather than color grading. 4d textures could be used as lookup tables to convert color spaces like CMYK to some other space. Since CMYK has 4 components, a 3d texture wouldn't be enough.


Make Object created in a list blend in/dissolve out of a scene

It looks like there are a few steps to make this work:

Linearly interpolate the alpha of the mesh from transparent to fully opaque, or vice versa.

Note 1: It seems as if the alpha starts to fade a few seconds after the effect has started.
Note 2: The alpha of the "effect meshes" (covered in #2 and #3) should be changed independently of the original mesh. This will help produce the desired effect.

Generate (just the color buffer for) a mesh with identical vertices, but make the color blue or green. The brightness of each vertice's/triangle's color should be calculated using either random noise (i.e. using a random number generator), or Perlin Noise; or a combination of the two. For the flashing effect, it would be best to regenerate the mesh each frame (with a different noise seed, if using Perlin Noise). Next, calculate a weighted average between the colors of this "effect mesh" and the those of the original mesh, according to how much time has passed. For instance, assuming t is the percentage of the way through the effect that has passed; at t=0, make the color buffer 0*effect_mesh + 1*original_mesh; and at t=0.3, make the color buffer 0.3*effect_mesh + 0.7*original_mesh.
For the vertical line effect, think about ray-tracing all the vertices (using multiple vertical rays). If a vertex is hit by (or near) a ray, make that vertex's color lighter/darker. This effect would be best done with a weighted average also, according to the time (since the lines are stronger towards the beginning).

*** Although this "Star Trek" effect varies from scene to scene, you should order these effects; don't do them all at the same time.


Logarithmic spiral with equal vertex spacing, what equations?


Anti-aliasing / Filtering in Ray Tracing

There is a great paper from 2006 on this topic, Filter Importance Sampling. They propose your method 2, study the properties, and come out generally in favor of it. They claim that this method gives smoother rendering results because it weights all samples that contribute to a pixel equally, thereby reducing variance in the final pixel values. This makes some sense, as it's a general maxim in Monte Carlo rendering that importance-sampling will give lower variance than weighted samples.

Method 2 also has the advantage of being slightly easier to parallelize because each pixel's computations are independent of all other pixels, while in method 1, sample results are shared across neighboring pixels (and therefore have to be synchronized/communicated somehow when pixels are parallelized across multiple processors). For the same reason, it's easier to do adaptive sampling (more samples in high-variance areas of the image) with method 2 than method 1.

In the paper, they also experimented with a Mitchell filter, sampling from abs() of the filter and then weighting each sample with either +1 or −1, like @trichoplax suggested. But this ended up actually increasing the variance and being worse than method 1, so they conclude that method 2 is only usable for positive filters.

That being said, the results from this paper may not be universally applicable, and it may be somewhat scene-dependent which sampling method is better. I wrote a blog post investigating this question independently in 2014, using a synthetic "image function" rather than full rendering, and found method 1 to give more visually pleasing results due to smoothing high-contrast edges more nicely. Benedikt Bitterli also commented on that post reporting a similar issue with his renderer (excess high-frequency noise around light sources when using method 2). Beyond that, I found the main difference between the methods was the frequency of the resulting noise: method 2 gives higher-frequency, "pixel-sized" noise, while method 1 gives noise "grains" that are 2-3 pixels across, but the amplitude of noise was similar for both, so which kind of noise looks less bad is probably a matter of personal preference.


When should quaternions be used to represent rotation and scaling in 3D?

I want to start with misconceptions:

Modern GPUs (NVIDIA for quite a while, and AMD since Southern Islands) do not meaningfully support vector/matrix operations natively in hardware. They are vector architectures in a different direction: each component of a vector (x, y, z) are generally 32- or 64-valued, containing values for each element in a lane. So a 3D dot product is not usually an instruction, it is a multiply and two multiply-adds.

Additionally, counting primitive operations like multiply-add, transforming a vector by a quaternion is more expensive than transforming a vector by a matrix. Transforming a vector by a 3x3 matrix is 3 multiplies and 6 multiply-adds, and transforming a vector by a quaternion is two quaternion multiplies, each of which consist of 4 multiplies and 12 multiply-adds. (You can get less naïve than this—here's a writeup on a faster way—but it's still not as cheap as multiplying a vector by a matrix.)

However, performance is not always determined simply by counting the number of ALU operations it performs. Quaternions require less space than the equivalent matrix (assuming you are only doing pure rotation/scale), and that means less storage space and less memory traffic. This is often important in animation (which is conveniently also often where the nice interpolation properties of quaternions show up).

Other than that:

Matrices use more space because they support more operations. A 3x3 matrix can contain nonuniform scale, skew, reflection, and orthogonal projection.
Matrices can be naturally thought of as basis vectors, and easily constructed from those vectors.
Multiplying one quaternion by another (composing two rotations) is less operations than multiplying one matrix by another.


Why is the line from the camera to vanishing point parallel to the other parallel lines?

FYI every set of parallel lines will have their own vanishing point. The reason why there is usually only a single vanishing point (the most obvious one) talked about in an image is due to aesthetic choice.

A perspective projection can be seen as a point projection of a 3D scene onto a 2D plane. The way to determine where a point in the 3D scene ends up on the 2D plane is to draw a line between the point and the camera (the sight line) and see where it intersects the projection target plane.

To find the vanishing point of a set of parallel lines you can take a the limit of traveling on that line infinitely far. At the limit the sight line will be parallel to the line.


Does an object spanning percentage of the Camera's FoV correlate with the object's image percentage?

@ Physician --- We can use trigonometry and set-up a condition whereby the horizontal angle of view will be 100°
On a full frame 35mm camera format demotions 24mm by 36mm, mount a 15mm lens and horizontal angle of view will be 100°
On an APS-C compact, the frame size is 16mm by 24mm. Mount a 10mm lens and the horizontal angle of view is 100°
In other words, the angle of view of camera system intertwines focal length and format dimension.
For the full frame camera above:
The image plane falls 15mm downstream from the lens.  The Image triangle has a height of 15mm. The size of the angle at the apex of this triangle is 50°. We bisect this triangle creating two right triangles. The angle at the apex is 25°. The height of this triangle is 15mm. If you draw this triangle, the height is the adjacent side and the base is the opposite size. The tan of 25° = 0.4663. Thus Opposite / Adjacent = 0.4663. Thus Opposite side = 15mm X 0.4663 = 6.9946mm. The length of the image of the object will be 2X this value = 14mm.
For the APS-C format:  The image plane falls 10mm downstream. Such as lash-up delivers a 100° horizontal angle of view. The image triangle angles are the same as for the full frame. Bisecting the image triangle, we get two right triangles each with a 25° apex. The tan of this angle = 0.4663. The height of this triangle is the focal length of the lens = 10mm. 10mm X 0.4663 = 4.66mm X 2 = 9.33mm. The image of the object will be 9.33mm long.


Problem on understanding gamma correction

Almost all image formats store colors gamma-encoded, so if you write out those values to a file then that's what you'll see on opening it. The pixel pickers etc operate on the same values stored in the file. It's much less common to store linear-light color values in a file (one example is digital camera RAW files).

When you "see color C1" on the screen, you are seeing C1^gamma. Or rather, you are seeing however that particular monitor renders the RGB value C1, which is at least approximately described by C1^gamma for some gamma. The color you see does not differ between the paint program and the programmatically generated file because they're both sending C1 to the monitor, and the monitor is doing the same thing with it.

Re: your update, you wrote:

  if you open this file on a viewer you will see pixels with brightness half as white

Think about this carefully. How are you judging that the pixels are "half as white"? Are you measuring them with a photometer and finding that they emit half the luminous intensity of a full white pixel? Or are you looking at them and judging that they appear roughly halfway between black and white?

Human perception of brightness is very non-linear. A value that looks perceptually halfway between black and white is physically closer to 20% of the brightness of white, not 50%.

A good way to find a true 50% brightness is to create an image with a checkerboard of black and white pixels, then defocus your eyes slightly, or stand farther away from the monitor, so that the pattern blurs into a uniform grey. That grey is 50% brightness (because it's an average over 50% black pixels and 50% white pixels). You will find that this grey looks considerably brighter than an RGB (128, 128, 128) grey.

For example, open this image in a new tab to be sure it's appearing at 1:1 resolution (browsers resizing it do not necessarily respect gamma):

You will see that the checkerboard in the center looks a lot brighter than the 128 grey on the left, and should be pretty close to the 186 grey on the right. The 186 is encoding a 50% brightness at 2.2 gamma.

Re: update 2, you wrote:

  So with input of 186 the screen outputs 50% of white intensity? And the eye perceives it as 50% as bright as white?

Not quite—the eye perceives a 50% physical intensity as more than halfway between black and white. Maybe around 70% or so. Conversely, the grey that the eye perceives as halfway between black and white is more like a 20% physical intensity.

Actually, gamma-encoded values are a fairly close match to how we perceive brightness, which is convenient (although it's sort of a coincidence; the original reason for gamma was due to how CRT phosphors worked, not anything to do with perception). For a true perceptually uniform color space, you can look at CIELAB, which has a more complicated nonlinear transformation that's designed to mimic average human perception.

  is it the programmer responsibility to write gamma-corrected values into the file or the whatever API we use to save image file formats does it internally?

It's the programmer's responsibility. Almost all image file formats, libraries, and APIs are just passing around gamma-encoded values and do not do any transformations on them.


Yes, because refractive index can vary with wavelength. This is the origin of colored specular reflection in metals such as gold and copper; most other materials have essentially uncolored specular. At each wavelength, the specular reflectance at normal incidence is related to the refractive index at that wavelength according to the formula you mentioned.

In case you're interested, the site refractiveindex.info has data tables and graphs of refractive index versus wavelength for various materials. For example, here's gold:


Yes, the back buffer is for the GPU to write the frame that's in-progress, while the display controller sends the front buffer to the attached display. Swapping buffers doesn't really copy anything these days. The display controller has a register which tells it the address of the frame buffer to display. Swapping buffers simply puts the address of the old back buffer in this register, making it the front buffer. You should only swap buffers when the frame is complete (i.e. you've finished drawing all the objects) and you want it to be displayed. This frame will continue to be displayed until you swap buffers again.

In some configurations (e.g. in EGL it depends on what swap interval you set) swapping buffers will make the CPU wait (block) until the GPU has finished rendering that frame. In others, the CPU can continue to queue up draw calls for further frames.


How does the view-independence of radiosity increase the calculations involved?

While it may not entirely clear from the formulation of the wikipedia article, the author raises an important issue: In contrast to many other approaches, Radiosity needs to perform its calculations for all existing patches, not only the visible ones.
It is not the limitation to diffuse surfaces that increases the involved computations, but the fact that radiance is computed for all surfaces in the entire scene, not only the visible ones.

This is in strong contrast to other global illumination techniques like path-tracing, where radiance is computed only for visible samples. While the view-paths may still reach every point of the scene, there might be parts of a scene that are never reached by any view-rays/paths. Therefore, there are no computations at all. Comparing with local GI the "problem of viewpoint independence" its even more apparent.

On the other hand, as the Wikipedia article suggests, this may also be seen as a very useful property, since the computations do not need to be performed for each different viewpoint. This is not the case for most other techniques.


smallpt: when to stop diffuse ray scattering

p is gonna be zero for light sources. Light sources are defined by emission variable and color is zero for them. Hence the reflectance parameter calculated previous to this line p would be zero giving us obj.e


The key trick is the perspective divide, This divides the position by the w component.

If you swap z and w (using a identity matrix that has swapped 3rd and 4th rows) and then do a perspective divide you end up with a double pyramid with a angle between the side (FoV) of 90° if applied to a unit cube.

Once you have that transformation you can remap you can do further remapping pre and post the swap to make the corners of the input and output match up. The well known algorithm to create a perspective matrix has all those remaps inlined as an optimization.


Render with camera perspective off-center

Yes, you can use an off-axis projection matrix.
This is what I use in my code (note: I shift the centre upwards, not left as you do in your example.)

void camera_setAspectRatio(float aspect, float zNear, float zFar, bool offaxis)
{
        // create a projection matrix
        const float f = 1.0f / tanf(fovy_radians/2.0f);
        fovx_radians = 2.0f * atanf( aspect * tanf(fovy_radians/2.0f) );
        float* mout = proj.data;

        mout[0] = f / aspect;
        mout[1] = 0.0f;
        mout[2] = 0.0f;
        mout[3] = 0.0f;

        mout[4] = 0.0f;
        mout[5] = f;
        mout[6] = 0.0f;
        mout[7] = 0.0f;

        mout[8] = 0.0f;
        mout[9] = offaxis ? -0.25f : 0.0f; // off axis projection!
        mout[10] = (zFar+zNear) / (zNear-zFar);
        mout[11] = -1.0f;

        mout[12] = 0.0f;
        mout[13] = 0.0f;
        mout[14] = 2 * zFar * zNear /  (zNear-zFar);
        mout[15] = 0.0f;
}

Stereo Projection screens (not head-mounted) tend to use Off-Axis projection matrices as well. Paul Bourke has a great write up on these projections.


Is the microfacet GGX BSDF normally implemented as separate BRDF and BTDFs?

In current common material models, Diffuse and Specular material properties are handled separately.

Check out slide 43 on http://blog.selfshadow.com/publications/s2015-shading-course/hoffman/s2015_pbs_physics_math_slides.pdf

Basically use a diffuse BRDF (lambert, oren nayar) if the material's refraction is scattering is less than 1 pixel. Use refraction and subsurface models if the scattering is larger than a pixel. I.e. in a glass refraction case the light exit is much greater than 1 pixel away from the light entrance.

You may also be interested in a diffuse parameterization for GGX so that you can use share the same roughness term. https://www.gdcvault.com/play/1024478/PBR-Diffuse-Lighting-for-GGX


How can I account for specular interactions in shadow rays?

Shadow rays that hit a glass surface or similar are stopped because there's no direct path from the light source to the shading point—there's a piece of glass in the way.

for perfectly specular transmission/reflection, you can continue tracing the shadow ray because there's only one possible outgoing vector for a given incoming vector

That won't work, because the transmission/reflection will change the direction of the ray. The shadow ray was presumably created by pointing it at a light source, but after transmission/reflection, it's most likely not going to hit the light source anymore.
The idea of firing rays directly at light sources can be seen as a sampling strategy designed to exploit knowledge of which directions light is most likely to come from, and put your rays into those directions. If you're receiving light via a specular interaction rather than directly from the source, you don't know anymore which directions you should sample in.
I imagine for simple cases like a glass sphere, etc, you could algebraically solve for a light path that would refract toward a chosen point on a light source. Then you could trace shadow rays along that path. But when you get into geometry as complex as the dragon, there's no realistic way to calculate which directions the light is going to come in from.
When specular surfaces like glass or metal focus light onto a diffuse surface, that produces a caustic. Caustics are a notorious pain point for vanilla path tracing due to this problem—there is no way to know where to concentrate your samples to make your rays hit light sources, so the best you can do is fire rays in all directions (or sample from the BRDF) and hope that enough of them will find their way to a light source. Technically, path tracing can render caustics correctly, but it takes an extremely large number of samples to do it.
More practically, this issue can be solved by using photon mapping (which was invented precisely to deal with caustics), or some kind of path guiding methods like Metropolis and friends that can "discover" the caustics and dynamically concentrate samples around them.


Given two 3d vectors, find the rotation matrix between the two

You can construct a rotation matrix from an "axis", or 3 vectors. This is done by calculating 3 direction (normalized) vectors for the 3 axis of our new rotated coordinate system, they are forward, up and right vectors.
In your case let's say we have 2 vectors called v1 and v2.
we can produce a direction from them via (glsl psuedo code):
vec3 forwardVector = normalize(v2 - v1);

Then for our up vector we can pick the world-up axis (e.g. 0,1,0), although this may give you problems if the forward vector is very close to the up axis, in that case you can pick another arbitrary axis that's not nearby, like {0,0,1}. For the sake of simplicity we'll stick to the world up axis of {0,1,0}.
The last vector can be derived automatically from the other 2 vectors using a cross product.
so now we should have something like:
vec3 forwardVector = normalize(v2 - v1);
vec3 upVector = vec3(0,1,0);
vec3 rightVector = normalize(cross(forwardVector, upVector));

These 3 vectors can be directly plugged into a 3x3 matrix to form the rotation matrix.
r.x r.y r.z
u.x u.y u.z
f.x f.y f.z

where f=forwardVector, u=upVector and r=rightVector.
This is similar to an openGL "LookAt" matrix.
Be aware there are many rotation matrices that are valid solutions because 2 vectors cannot describe what the rotation is along the v1->v2 axis (That's what the 'up' axis fixes in this case).


Why do Vertex Buffer Objects improve performance?

In general, when you render an object in an immediate mode—issuing line drawing commands for instance—you build up a series of commands that you submit to the graphics card to draw. If you're drawing a lot of data, or drawing very frequently, you can waste a lot of time sending this data over and over again.

A vertex buffer allows you to produce a single object that you submit to the graphics card once. If you don't need to change your geometry, you can leave it on the graphics card and simply send the graphics card a request to draw that object. Since it avoids the copy every time you draw, there's much less overhead for each draw.

Note that using a vertex buffer object doesn't always provide a very significant speedup. If you're only drawing the object once per frame, and you're replacing the geometry in between each frame, then you don't get the benefits of avoiding the copy each frame.

Most of my experience comes from writing programs using graphics APIs like OpenGL, so someone who's messed with the backend of a graphics driver can probably provide a more detailed answer, but I hope this makes things a bit clearer.


Shear(?) matrix projection along Z (depth) axis

I belive what you are looking for is a scale Matrix, or actually it will end upp with as a shear matrix for you.
Usually they look like this
Sx  0   0   0
0   Sy  0   0
0   0  Sz   0
0   0  0    1

If you have no scaling, Sx, Sy, Sz represent the scaling in corresponding dimension.
So put the to 1 for no scaling.
But now you want to scale the height depending on the z value. So we want to modify it abit. If you would put Sx, Sy, Sz to the value one. And the put a higher value than one on Syz, then it would scale the height (Y), the futher away you are.
Sx  0    0   0
0   Sy  Syz  0
0   0   Sz   0
0   0   0    1


Are there tricks for getting proper sort ordering on particle systems?

First option should be to make all particles able to go through the same pipeline. Perhaps with an uber shader. That way you can batch them all.

Positive ieee floating point numbers can be sorted like unsigned integer. And there are O(n) algorithms to sort those for example radix sort.


When you are making an SDF, you commonly have infinitely detailed data (like, bezier curves making up a letter) that you are making into an SDF image of a specific size.

If you make one that is 64x64, the first mip will be 32x32. The right way to make that mip is to do the original process of creating the SDF texture from the source data, but do it at 32x32 instead of 64x64. Repeat for the other mips.

It will break down at lower resolutions though, but that is just the nature of how SDFs work. The more correct, but slower way is to do super sampling. Take multiple samples at different places in the pixel and average the result.


Can D3D12 resource barriers be used to synchronize between the compute and 3D engine?

My understanding would be that resource barriers are not at all about synchronization. They are about memory ordering, which is a related but separate concern. Basically, synchronization is about making sure that A happens before B. Memory ordering is about making sure that if A did X, someone how saw A happen will always also see X and no one will ever see X without first also having seen A happen.

In Direct3D terms, resource barriers are about making sure that if shader A writes to texture X, shader B running after shader A actually gets to see the result of shader A having written to texture X. Fences, on the other hand, are about making sure that shader A actually runs before shader B.

Direct3D resource barriers are just you stating in which way you are going to be using a resource. You cannot rely on just a resource barrier to automagically perform the necessary synchronization for you. The section you linked would seem to agree (emphasis mine):

  Although resource state is shared across all Compute and 3D queues, it is not permitted to write to the resource simultaneously on different queues. "Simultaneously" here means unsynchronized, noting unsynchronized execution is not possible on some hardware. The following rules apply.

  Only one queue can write to a resource at a time.
  Multiple queues can read from the resource as long as they don’t read the bytes being modified by the writer (reading bytes being simultaneously written produces undefined results).
  A fence must be used to synchronize after writing before another queue can read the written bytes or make any write access.


What is a tbuffer and what is their use case?

Keep in mind, that I'm still learning and I have no professional capacity in terms of graphics programming. So you might want to continue research with provided resources.

t*#*: A register for a texture buffer (tbuffer).[1]

To start with - texture buffer is what you can see in hlsl code as: Texture2D texName : register(t0) or as Texture2D texName : register(t0, space0). If I remember correctly, spaces were added in DirectX 12, maybe they're available in DX 11.3, but that's different discussion. Back to my point - you're using texture buffer every time you're binding texture to your shader and sending texture data to GPU. Therefore, using PSSetShaderResources [2] will be most trivial usage of texture buffer, espacially in DX11 code. You can use constant buffer to store pixel data and then send it to GPU, but why should you bother, when Texture2D or TextureCube and other data structures are already in hlsl?

cbuffer (DX11/DX12) forces you to pack data in 32-bit groups [3]. You might have experienced it, if not, I encourage you to experiment a little. In the past I didn't follow this rules, having that packed like this, which causes problem with incorrect data in your shader:
cbuffer buf : register(b0)
{
    float3 a;
    // Here is where problem occurs
    float2 b;
    float3 c;
};

Correct packing rule would be:
cbuffer buf : register(b0)
{
    float3 a;
    float padding_a;
    float2 b;
    float2 padding_b;
    float3 c;
    float padding_c;
};

Of course, I am using worst possible example here, because correct data is 8 bytes and you're sending additional 4 bytes (whole hlsl vector) as padding. But you get an idea. It might be easier to sort data with bigger buffers, especially if you're using uber buffer in your hobby project.

Additional problem in DX12 is that cbuffer has to be 256-byte aligned [4], unless you're using root constants, which allows you to send not aligned data, but they have other limitation. I was only talking about them with other developer, so I don't really have sources except [5]. There was something about size of buffer and that root constants are working magic in AMD drivers, but it isn't something that I really need right now, as I am focusing on learning, not extreme optimization of code.

Shader Model 5 allows you to use new data structures [6]. I wasn't working with compute shaders, but for DXR hlsl code, I am using RTTexture2D as UAV. However, I am passing indices, vertices and RaytracingAccelerationStructure as texture buffers. I didn't see anyone to use those three as cbuffer. I guess it's possible, at least for indices and vertices, but I don't see a reason to do so:
RaytracingAccelerationStructure SceneBVH : register(t0);
ByteAddressBuffer indices : register(t1);
StructuredBuffer<Vertex> vertices : register(t2);

Indices are simply 16 or 32 bit index for accessing vertex data. I don't see a reason to hassle with creating cbuffer for single array and having to deal with padding and providing correct size of buffer. However, I am really interested if there is a reason to do it sometimes.
Same with vertices - in rasterization, you're using structure which don't have to follow packing rules e.g.:
struct VertexInputType
{
    float3 position : POSITION;
    float3 normal : NORMAL;
    float3 tangent : TANGENT;
    float3 binormal : BINORMAL;
    float2 uv : TEXCOORD0;
};

So for raytracing, I am just using struct with same data and I'm accesing it through structured buffer, regsitered as texture buffer:
struct Vertex
{
    float3 position;
    float3 normal;
    float3 tangent;
    float3 binormal;
    float2 uv;
};
StructuredBuffer<Vertex> vertices : register(t2);

Regarding RaytracingAccelerationStructure - it's built-in hlsl type, no need to hassle with creating cbuffer.

I'm sorry for speaking from DX12, DXR perspective, but that's something that I'm currently working on, so it's a simplest way to provide examples from that area.
[1] https://docs.microsoft.com/en-us/windows/win32/direct3dgetstarted/work-with-shaders-and-shader-resources#read-from-the-constant-buffers
[2] https://docs.microsoft.com/en-us/windows/win32/api/d3d11/nf-d3d11-id3d11devicecontext-pssetshaderresources
[3] https://docs.microsoft.com/en-us/windows/win32/direct3dhlsl/dx-graphics-hlsl-packing-rules
[4] https://social.msdn.microsoft.com/Forums/sqlserver/en-US/78dec419-6827-40fc-a048-01c3cccb92ef/directx-10-11-and-12-constant-buffer-alignment?forum=vsga
[5] https://docs.microsoft.com/en-us/windows/win32/direct3d12/using-constants-directly-in-the-root-signature
[6] https://docs.microsoft.com/en-us/windows/win32/direct3d11/direct3d-11-advanced-stages-cs-resources


Is it possible to find out how many fragments made it through the stencil test?

One possible approach could be the use of Hardware Occlusion Query.

You can use the facts that, by specification, the Stencil Test is executed before the depth test, and only the fragments that pass the depth test are counted by the Occlusion Query.

A simple example (not tested) would be like:

    GLuint samples_query = 0;
    GLuint samples_passed = 0;
    glGenQueries(1, &samples_query);
    // Initialize your buffers and textures ...
    glEnable(GL_DEPTH_TEST);
    glEnable(GL_STENCIL_TEST);

    // Set up the values on the stencil buffer ...

    // Now we count the fragments that pass the stencil test
    glDepthFunc(GL_ALWAYS); // Set up the depth test to always pass
    glBeginQuery(GL_SAMPLES_PASSED, samples_query);
    // Render your meshes here
    glEndQuery(GL_SAMPLES_PASSED);
    glGetQueryObjectuiv(samples_query, GL_QUERY_RESULT, &samples_passed);
    // samples_passed holds the number of fragments that passed the stencil test (if any)

    // Release your resources ...
    glDeleteQueries(1, &samples_query);

Note that the call to obtain the number of samples will call forcibly the flush of the pipeline and wait for the query to finish.
If you need a more asynchronous approach you can query wether the occlusion query is done or not by using:

    GLuint query_done = 0;
    glGetQueryObjectuiv(samples_query, GL_QUERY_RESULT_AVAILABLE, &query_done);
    if (query_done != 0)
        // Your query result is ready
    else
        // Maybe check the next frame?


Check the last column in the table for the deciphered version of the existing Map Equations for easier comprehension.


Is there an analogue of OpenGL's ImageLoadStore in DirectX 12?

You should be able to use RWTexture3D.
The hlsl documentation page about RWTexture2D also has sample code that you might find useful.


How do I implement object picking, using OBB in OpenGL?

Ray vs OBB intersection works by shooting a ray towards the OBB's planes. This results in 2 intersections referred to as t1 and t2 (per axis), also referred to as tNear and tFar or tMin and tMax.

Checking whether the ray intersects with the box, can then be done by comparing the largest tNear against the smallest tFar. If the largest tNear is less than the smallest tFar then there is a intersection.


You can't see it in the image you uploaded because it's a JPEG, but PNG files which have a full alpha channel can do the trick you mentioned. It works with websites or apps where the "thumbnail" image is displayed on a white background and the "full" image is displayed on a black background. In the thumbnail display, white pixels and completely transparent pixels look the same, because of the white background; in the full image, black pixels and completely transparent pixels look the same.

Thus, you can have a white-on-black picture and a black-on-white picture in the same image, and the background hides one of the two images. You can easily see what's going on if you look at the image against a background of some other colour.


It is  common that the camera is pointing towards negative Z, in order to satisfy right hand rule. So you need to know the convention.

As for theinverse matrix of the camera matrix. Sort of yes. You can do it. However the camera matrix is a homogeneous space. While generally our other matrixes are affine spaces. If you take care of the distinction and do the proper camera divide, if tge camera needs one. Then yes you could do that.


Overlaying graphics on particle simulation

For rendering of gases, I think the usual approach is to simply render each particle as a tiny disc. Gases don't really coalesce into surfaces like liquids do, so this should produce acceptable results. You could perhaps apply a light blur over the gas layer afterwards to soften it and hide the fact that it is made of discrete elements.

Liquids, on the other hand, tend to coalesce together to form droplets and smooth surfaces, so you need to derive a surface from the particles somehow. One way to do this is to use Metaballs, which also display this behavior and can be tweaked to suit different liquids and particle densities. By interpreting each particle as a metaball, you will have an implicit equation representing your liquid surface. To render this implicit surface you can then use an algorithm like Marching Tetrahedra to convert it to triangles, or utilize Ray Marching to directly render it. (Ray marching can be easily done in realtime in a fragment shader these days.) You can also use this approach for gases if you want a somewhat softer look.


There is no such thing as a core profile in GLSL 1.30. The concept of API profiles was introduced with OpenGL 3.2 which happened after GLSL 1.30 was released. The syntax for the version specification in GLSL 1.30 is

#version <number>

Apart from that, I highly recommend that you check to confirm that GLSL 1.30 really is the maximum version you're supposed to be able to use on your system. It seems not unlikely that you just don't have a proper OpenGL driver installed. Even software implementations like MESA support more recent versions than that these days…


Hint:

The three coordinates can be handled independently.

For instance in $x$, you need to solve

$$x_0 + t (x_1-x_0) = 0$$
 and
$$x_0 + t (x_1-x_0) = x_{max}$$


Gaussian blurring a function defined over a surface mesh

A few different approaches

I'll consider a few variations on your specific request, since you mention efficiency and I suspect your specific request may be the least efficient. I'll also suggest ways of improving the efficiency without varying from your intended approach, so you can weigh up the alternatives.

Blurring the volume instead of the surface

If you want the distance metric to be 3D Euclidean distance instead of 2D Euclidean distance within the surface, then you could perform the blur on a regular 3D grid to which the scalar function you have in mind has been applied. Then you can use the final result of your Difference of Gaussians to calculate the scalar values at the vertices of your irregular mesh. This avoids having to take into account the mesh shape for the bulk of the calculation.

The 3D grid is likely to have a much larger number of vertices than the 2D mesh, but they will all be equally spaced and the large kernel blur can be achieved by repeated application of a small kernel blur taking into account only the 6 nearest neighbours, which will always be at a constant distance away. This approach involves potentially more calculation, but the ease with which the regular grid could be GPU accelerated may appeal.

This will give a different result from performing the blur on the mesh vertices using 3D Euclidean distance. For example, the 3D grid approach will be affected by distinctive regions of the 3D scalar function that are near but not on the mesh. This may be desirable or not depending on your specific purpose.

Using 2D distance instead of 3D distance

If you find that you need the distance metric to be 2D Euclidean within the surface, then you can get a good approximation to a larger kernel Gaussian blur by repeatedly applying a smaller kernel Gaussian blur. If there is not too much variation in the edge lengths within your mesh you may be able to choose a kernel size which allows for only including vertices one edge away at each iteration. This allows for only using single edge lengths to calculate the size of the contribution of a vertex, rather than calculating a 2D multi-edge distance.

3D distance using the surface without the volume

If you need the calculation to be precisely as described in your question, being calculated within the mesh rather than within the surrounding volume, but also using the 3D Euclidean distance, then using nearest neighbours and several iterations will not work. Unless the mesh is near to flat, the repeated application of a nearest neighbour blur will result in an approximation to the 2D Euclidean distance case, since the values will only be able to bleed from vertex to vertex, not directly along the shortest path as they would in a single pass. This will give less spread than would be achieved by a single pass that calculates the 3D distance to a vertex 10 edges away. (I have used 10 edges since you mention a 10-hop in your comment on the question.)

Implementing the blur in a single pass will mean calculating the 3D Euclidean distance between every vertex and every other vertex within a 10 edge radius. This will be expensive, but perfectly possible. Since you mention efficiency, consider that there are some redundancies you can eliminate provided you have sufficient available memory.

The two blurs that you produce prior to taking the Difference of Gaussians will use the same set of 3D distances up to the edge radius of the smaller kernel blur. If you can save these then you only need to calculate them once, rather than once per blur.

Also, each distance will be used twice per blur - once in each direction as the length from vertex A to vertex B is the same as the length from B to A. Caching/memoising these distances will avoid calculating them twice.

Effects from arbitrarily many edges away

If the surface curves such that some vertices which are many edges away are still near enough to affect each other in 3D distance, then rather than considering vertices within a certain number of edges away, you may need to consider all vertices within a certain 3D radius, regardless of how long the path via edges. In this case you can consider fewer vertices by using space partitioning, choosing a specific method which suits the mesh.

If you don't want parts of the surface which approach each other to influence each other, then you probably want a 2D distance metric rather than the 3D one.

If you have a wide range of different edge lengths then you may find the same problem of not being able to define a set number of edges to traverse, even if the mesh is fairly flat. Again you may need to define a 3D distance instead of a number of edges, and consider all vertices that lie within that radius.


Strange sphere reflection effects in ray tracing project. I cannot figure out what causes them

The issue was caused by an incorrect calculation of the reflection direction vector.
With D ray direction and N the normal vector:

R = D - 2 * dot(D, N) * N

The issue was caused by calculating the components of R as follows:

R[i] = D[i] - 2 * (D[i] * N[i]) * N[i]

It took me a while to find the mistake because this produced a correct reflection with the horizontal plane.


According to Polygon-Mesh-Processing, section 7.2.2. paragraph "Error quadrics"

  Error quadrics have major advantages regarding memory consumption and computational efficiency: each vertex $x$ stores a symmetric 4x4 matrix $Q$ only, and the error can efficiently be computed in constant time as $x^T Q x$ - no matter how many planes are associated to vertex $x$ during decimation. Because of this, error quadrics are one of the most frequently employed techniques in mesh decimation.

The author of the reference I mentioned specifically refers to the paper you mention in your question. My interpretation is that more recent methods might be more sophisticated but resource-wise quadric errors are still better.

Also I'd add is relatively easy to implement anyway.


Transforming a ray from camera space to world space

Now, we can write this direction vector in homogeneous form and apply the camera to world transformation. This should give us the direction vector for the ray in world space.

That is the problem. When transforming a direction, you must ensure the w coordinate is 0, so that it doesn't get translated by your camera-to-world matrix. What you're doing in the second method is that you're transforming the direction vector as a point, and then subtracting the camera's position, which effectively undoes the translation.


In my answer I relate to WebGL 2 / Opengl ES 3.0.

As far as I know ARB_shader_image_load_store is not supported. But there is discard available for fragment shaders in WebGL 2! So what you can do is to

make a depth pass with a straight-forward fragment shader and
then in a second pass test gl_FragCoord.z against the value of the depth texture.

It's up to you, if it's worth the additional texture lookup and rendering the scene twice.


How is Anti Aliasing Implemented in Ray Tracing?

Raxvan is completely right that "traditional" anti aliasing techniques will work in raytracing, including those that use information such as depth to do antialiasing.  You could even do temporal anti aliasing in ray tracing for instance.

Julien expanded on Raxvan's 2nd item which was an explanation of super sampling, and showed how you'd actually do that, also mentioning that you can randomize the location of the samples within the pixel but then you are entering signal processing country which is a lot deeper, and it definitely is!

As Julien said, if you want to do $N$ samples per pixel, you can break the pixel up into $N$ evenly distributed sample points (on a grid basically) and average those samples.

If you do that, you can still get aliasing though.  It is better than NOT doing it, because you are increasing your sampling rate, so will be able to handle higher frequency data (aka smaller details), but it can still cause aliasing.

If you instead take $N$ random samples within a pixel, you are effectively trading aliasing for noise.  Noise is easier on the eyes and looks more natural than aliasing, so is the preferred result usually.  I believe it's even provably the ideal situation with higher sample counts but don't have more info on that ):

When you use just "regular" random numbers like you'd get from rand() or std::uniform_int_distribution, that is called "white noise" because it contains all frequencies, like how white light is made up of all other colors (frequencies) of light.

Using white noise to randomize the samples within a pixel has the problem that sometimes your samples will clump together.  For instance, if you average 100 samples in a pixel, but they ALL end up being in the upper left corner of the pixel, you aren't going to get ANY information about the other parts of the pixel, so your final resulting pixel color will be missing information about what color it should be.

A better approach is to use something called blue noise which only contains high frequency components (like how blue light is high frequency light).

The benefit of blue noise is that you get even coverage over the pixel, like you get with a uniform sampling grid, but, you still get some randomness, which turns aliasing into noise and gives you a better looking image.

Unfortunately, blue noise can be very costly to compute, and the best methods all seem to be patented (what the heck?!), but one way to do this, invented by pixar (and patented too i think but not 100% sure) is to make an even grid of sample points, then randomly offset each sample point a small amount - like a random amount between plus or minus half the width and height of the sampling grid.  This way you get a sort of blue noise sampling for pretty cheap.

Note that this is a form of stratified sampling, and poisson disk sampling is a form of that too, which is also a way of generating blue noise:
https://www.jasondavies.com/poisson-disc/

If you are interested in going deeper you probably will also want to check out this question and answer!

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

Lastly, this stuff is starting to stray into the realm of monte carlo path tracing which is the common method for doing photorealistic raytracing.  if you are interested in learning more about that, give this a read!

http://blog.demofox.org/2016/09/21/path-tracing-getting-started-with-diffuse-and-emissive/


Algorithms to Remove High Frequency Noise from Path Tracing

There are, and I am looking forward to seeing the specifics of other answers, but one way to deal with this is to not have the noise (or as much noise) in the source data to begin with.

The noise is coming from the fact that there is high variance in the rendering - the number of samples you've taken haven't converged enough to the actual right answer of the integral, and so some pixels are too high/bright and some are too low/dim (in each color channel).

The problem is this: If you use white noise random numbers to do your sampling, you may get samples clumping together like the image below.  Given enough samples, it will converge, but it will take a while before it gives good coverage over the sampling space.  Find a region of empty space in the image below (like in the lower right) and imagine that there was a small, bright light there and that the scene was dark everywhere else.  You can see how not having any samples there is going to make a problem for rendering.

Alternately, you could sample at even intervals like the below, but that will give you aliasing artifacts instead of noise, which is worse.

One idea is to use low discrepancy sequences and do quasi monte carlo integration (https://en.wikipedia.org/wiki/Quasi-Monte_Carlo_method).  Low discrepancy sequences are related to blue noise, which has only high frequency components.  Going these routes, you get faster convergence of $O(1/N)$ instead of $O(\sqrt{N})$.  These give better coverage of the sample space, but since there is some randomness (or random like qualities) to them, they don't have the aliasing issues that regularly spaced sampling does.

Here is a "jittered grid" where you sample on a grid, but use small random offsets within a cell size.  This was invented by pixar and was under patent for a while but is no longer:

Here is a common low discrepancy sequence called the Halton sequence (basically a 2d version of Van Der Corpus)

And here is a poisson disc sampling, using Mitchel's best candidate algorithm:

More information, including the source code that generated these images can be found here: https://blog.demofox.org/2017/05/29/when-random-numbers-are-too-random-low-discrepancy-sequences/


SSBO get block of data instead of iterating

yes, and no. Loading the floats into the SSBO is loading them for the shader so yes you still have to iterate over them. And GPU's generally are good at latency hiding so whether or not you read floats or vec4's the GPU should be able to read them in just as fast. But processing them as vec4's instead of 1 float at a time will be faster since the code is taking advantage of the GPU's vector processing.


I'm not going to dive into much details about affine transformations and such, you are better off reading up on these concepts from good graphics books.

I think The matrix you need is this.

$\begin{bmatrix}
0.5 & 0 & 0 & 4\\
0 & 0.5 & 0 & 6\\
0 & 0 & 0.5 & 2\\
0 & 0 & 0 & 1\\
\end{bmatrix}$

The reason why we need 4x4 matrices is because of affine transformations. For more details on why we have 4x4 matrices and why are we using 4D representation you can check my answer here. Briefly speaking, Vectors representing direction have an extra 0 at the end since translations don't affect these sort of vectors and Vectors representing points have an extra 1 at the end. This makes our vectors 4D so we can apply 4x4 transformation matrices. Since here we are concerned with points we will append an extra 1 to each vectors.

Let's begin by how did I derive it. First of all we can see that they provided the origin of the $UVW$ frame in terms of $IJK$ frame. This is $[4i, 6j, 2k]$. In other words the vector $[0,0,0]$ in $UVW$ frame must correspond to $[4i, 6j, 2k]$ in $IJK$ frame. Let's forget the 0.5's in the diagonal for now. We have,

$\begin{bmatrix}
0 & 0 & 0 & 4\\
0 & 0 & 0 & 6\\
0 & 0 & 0 & 2\\
0 & 0 & 0 & 1\\
\end{bmatrix}
*
\begin{bmatrix}
0\\
0\\
0\\
1\\
\end{bmatrix}
=
\begin{bmatrix}
4\\
6\\
2\\
1\\
\end{bmatrix} $

Next the condition is if we move 1 unit in $IJK$ we move 2 units in $UVW$ or in other words if we move 1 unit in $UVW$ we move 0.5 units in $IJK$. This means if we move from origin to $[1,1,1]$ in $UVW$ we moved half a unit in the $IJK$ frame in each of the axis. Testing the original matrix now we have,

$\begin{bmatrix}
0.5 & 0 & 0 & 4\\
0 & 0.5 & 0 & 6\\
0 & 0 & 0.5 & 2\\
0 & 0 & 0 & 1\\
\end{bmatrix}
*
\begin{bmatrix}
1\\
1\\
1\\
1\\
\end{bmatrix}
=
\begin{bmatrix}
4.5\\
6.5\\
2.5\\
1\\
\end{bmatrix} $

When constructing matrices this way, you want the translation factor in the last column. This factor is given to you by telling you the origin of the $UVW$ frame in terms of $IJK$ frame. This means the $UVW$ frame is that much displaced/translated when viewed from $IJK$ frame.

Secondly any scaling factor goes along the diagonal like above. This factor was given to you by telling you "If you move $X$ units in 1 frame, you move $Y$ units in the other".
Which is what we want I guess. Is that what you need?


How to sample 3D points to visualize a B-spline surface?

Based on your linked question, the problem appears to be that you don't get enough samples in highly curved areas of the surface. You could increase the sample density everywhere, but then you oversample the flatter areas and it becomes too slow.

So, you need to increase the sample density locally in areas where it's needed. A conceptually simple algorithm for this is adaptive refinement. You would start out with a uniform sample grid in UV space. Then, detect grid squares in which the curvature is too high for the grid spacing (above some threshold), and subdivide them, adding additional samples in between the existing ones. This process can be repeated recursively on the new, smaller grid squares, subdividing as finely as the surface requires. You can tune the curvature threshold to balance fidelity versus performance.

Since B-splines (NURBS) are based on rational functions, I would guess that the maximum curvature in a region can be found (or at least bounded) analytically. This would be ideal, if possible, since it would avoid the problem of missing a small high-curvature area if none of the initial grid points happen to land in it. However, in a pinch the curvature can be estimated by looking at the sample points you've evaluated. For instance, if $p_1, p_2, p_3$ are three sample points collinear in UV space, they should be close to collinear as well, if the surface is adequately sampled. So, you could for instance evaluate
$$\frac{\|p_2 - p_1\| + \|p_3 - p_2\|}{\|p_3 - p_1\|}$$
where the double bars mean vector length. If this is greater than some threshold value such as 1.1 or 1.5 (note that 1.0 would be perfectly flat), then you would subdivide the area between these points.


Ray tracing - tangent space for a point on a sphere

The tangent space is spanned by the tangent to the point and the bitangent (which is orthogonal to both tangent and normal).

So you need to calculate the tangent which is achieved by calculating the cross-product of the ray-direction and the normal. $T = N \times DIR$
The resulting vector will be orthogonal to the normal and thereby be the tangent.

Now calculate the cross-product of the tangent and the normal $BT = T \times N$ to create a vector orthogonal to both. This vector is the bitangent.

Tangent $T$ and bitangent $BT$ span a plane which is the tangent-space of your intersection-point.


Opengl black screen for a simple translation

We fixed it ! it was glm::mat4(); no initializing to a identity matrix, glm::mat4(1.0f); works fine, thank you for helping me !


Properties of the image reconstruction filter in rendering

I am replying to the last edit of the question.

You can show that $\|h_j * f \|^2_2 \leq \|h_j\|^2_1\|f\|^2_2$ from Young's theorem. Provided that  $\|h_j\|_1 \leq 1$ you have: $\|h_j * f \|^2_2 \leq \|f\|^2_2$. Additionally for $g_j$ I believe you want to clamp the indicator function to $\{h_j > 0 \wedge f > 0\}$. Then $\int_{\{h_j>0 \wedge f>0\}} a\leq \int_{\{f>0\}}a$ (since the integrands are always non-negative: $a\geq 0$). This basically proves that the inequality holds for the left terms of the quadratic expansions, the middle terms are both negative, so those may be ignored. The only thing left to prove is: $(\frac{1}{\mu(\{h_j>0 \wedge f>0\})}\int h_jf)^2 \leq (\frac{1}{\mu(\{f>0\})}\int f)^2$, which is the case if $h_j \leq 1$ everywhere. Basically, as long as your filter $h_j$ preserves the minimum and maximum of the image, you're good (this is not an unreasonable requirement for a filter, note however, that there are filers that do not preserve min-max).

Edit: My bad, the middle terms may not be ignored (at least the one on the right of the inequality). Which means you are left to prove that:

$$\frac{(\int h_jf)^2}{\mu(\{h_j>0 \wedge f>0\})} \geq \frac{(\int f)^2}{\mu(\{f>0\})}$$, which will be a problem.

I guess one has to use tighter bounds to relate the mixed term with the squared terms.

Edit 2:

Let's consider the expression as a whole then:

$$\int \left(h_jf - \frac{\int h_jf}{\mu(\{h_j > 0 \wedge f > 0\})} \right)^2 \leq \alpha \int \left(f - \frac{\int f}{\mu(\{f > 0\})} \right)^2$$
$$ 0 \leq \alpha\int (f)^2 - \int (h_jf)^2 +\alpha\frac{(\int f)^2}{\mu^2(\{f > 0\})} - \frac{(\int h_jf)^2}{\mu^2(\{h_j > 0 \wedge f > 0\})} + \frac{2(\int h_jf)^2}{\mu(\{h_j > 0 \wedge f > 0\})} - \alpha\frac{2(\int f)^2}{\mu(\{f > 0\})}$$
$$ 0 \leq \alpha\left(\int (f)^2 - \frac{(1-2\mu(\{f > 0\}))(\int f)^2}{\mu^2(\{f > 0\})}\right) - \left( \int (h_jf)^2 - \frac{(1-2\mu(\{h_j > 0 \wedge f > 0\}))(\int h_jf)^2}{\mu^2(\{h_j > 0 \wedge f > 0\})}\right) $$

$$0 \leq \left(\int(\sqrt{\alpha}f)^2-\int (h_jf)^2\right) - \left(\frac{\alpha(1-2\mu(\{f > 0\}))(\int f)^2}{\mu^2(\{f > 0\})}  - \frac{(1-2\mu(\{h_j > 0 \wedge f > 0\}))(\int h_jf)^2}{\mu^2(\{h_j > 0 \wedge f > 0\})}\right) $$

$$0 \leq \left(\int(\sqrt{\alpha}f)^2-\int (h_jf)^2\right) - \left(a \left(\int \sqrt{\alpha}f\right)^2 - b\left(\int h_jf\right)^2\right)$$

Note also that $\left(\int u\right)^2 \leq \int u^2$. Then with an appropriate $\alpha$ and $h_j$ bounded it should hold. You may want to study the $a,b$ tems also.


sdl2: Compiler throws "Undefined Reference"

The problem is definitely a lack of linker flags.
It can see that you want to include <SDL/SDL2.h>, which is a sort of promise that said functions contained in that header have been built and will be available; but when your linker looks to find the actual machine code, it needs a little more than the header. It needs the library itself.
For SDL2, you would usually use something like:
gcc -g main.c -o ../ludo -lSDL2

On Linux, this is more-or-less the same command, but for complete definitions it also instructs the compiler to look inside /usr/lib/x86_64-linux-gnu/libSDL2.a (or wherever it is on your system). Note that the lib prefix on the file name is assumed, and the extension is irrelevant.
Do that, and your program should compile. Note that some SDL extensions, like Image and Mixer, require their own libraries. An SDL-specific shortcut to ensure that you've got everything is to use the handy pkg-config, like so:
gcc -g main.c -o ../ludo $(pkg-config --cflags --libs sdl2)

Which will be interpreted by the terminal as
gcc -g main.c -o ../ludo -D_REENTRANT -I/usr/include/SDL2 -lSDL2

This is less generic however, a bit more paranoid, and a bit of a cheat; but perhaps it will help you.
The thing to remember with C is that a function needs to be declared to be used in code, and that's compiler-end; but to be built into an actual stand-alone program it also has to be defined, which is where the post-compiler linker comes in.


Tracing shadow rays in a scene with mirrors

It depends how you define easy and what kind of constraints you have. The general case of this is rendering caustics but that's probably not what you're looking for if real time is your target.

If your mirrors are always flat as in your demo and you only want to support a single bounce, the easiest I can see would be to reflect your light sources on the other side of the mirror's plane, before you start a frame. Let's call each reflection a virtual light source. Then for each of these, test if the ray between whatever you're shading and the virtual light source hits inside the mirror (simple plane/ray test). If it does hit then you trace the actual shadow ray to the mirror and the rest of the shadow ray from the mirror to the actual light source.

This obviously won't scale well to many mirrors. I think it could be adapted to 2-3 bounces if the mirror and light count is low. The number of virtual lights could quickly grow out of hand though.


In short, performance reasons are why they aren't programmable.

History and Market

In the past, there used to be separate cores for vertex and fragment processors to avoid bloated FPU designs. There were some mathematical operations you could only do in fragment shader code for instance (because they were mostly only relevant for fragment shaders). This would produce severe hardware bottlenecks for applications that didn't max out the potential of each type of core.

As programmable shaders became more popular, universal units were introduced. More and more stages of the graphics pipeline were implemented in hardware to help with scaling. During this time, GPGPU also became more popular, so vendors had to incorporate some of this functionality. It's still important to note though that the majority of the income from GPUs were still video games, so this couldn't interfere with performance.

Eventually a big player, Intel, decided to invest in programmable rasterizers with their Larrabee architecture. This project was supposed to be groundbreaking, but the performance was apparently less than desired. It was shut down, and parts of it were salvaged for Xeon Phi processors. It's worth noting though that the other vendors haven't implemented this.

Attempts at Software Rasterizers

There have been some attempts at rasterization through software, but they all seem to have issues with performance.

One notable effort was an attempt by Nvidia in 2011 in this paper. This was released close to when Larrabee was terminated, so it's very possible that this was a response to that. Regardless, there are some performance figures in this, and most of them show performance multiple times slower than hardware rasterizers.

Technical Issues with Software Rasterization

There are many issues that were faced in the Nvidia paper. Here are some of the most important issues with software rasterizers though:

Major Issues

Interpolation:
The hardware implementation generates interpolation equations in specialized hardware. This is slow for the software renderer since it had to be done in the fragment shader.
Anti-aliasing:
There were also performance issues with anti-aliasing (specifically with memory). Information regarding the sub-pixel samples must be stored on-chip memory, which isn't enough to hold this. Julien Guertault pointed out that the texture cache/cache may be slower with software. MSAA certainly has issues here because it overflows the cache (the non-texture caches) and goes into memory off of the chip. Rasterizers compress data that is stored in that memory, which also helps with performance here.
Power Consumption:
Simon F pointed out that power consumption would be lower. The paper did mention that custom ALUs are in rasterizers (which would reduce power consumption), and this would make sense since fragment and vertex processing units in the past used to have custom instruction sets (so likely custom ALUs as well). It certainly would be a bottleneck in many systems (e.g., mobile), although this has implications beyond performance.

Summary

TL;DR: there are too many inefficiencies that software rendering can't past, and these things add up. There are also many larger limitations, especially when you are dealing with VRAM bandwidth, synchronization problems, and extra computations.


Trying to understand environmental cubemaps

Your confusion comes from fact that in some basic tutorials like one you point cubemap which is used for applying lighting is global and static - it does not come from actual geometry but from infite sky box which moves with camera position so it is valid everywhere.

And there is nothing mentioned about second type - local cubemaps which are actually computed using local geometry and are valid only in point of their placement.

This can be approached using multiple local cubemaps and connecting them with local objects situated around them but what is problematic is that it generates seams between them especially for high frequency reflections, which you probably would like to test.

For now you can generate few local cubemaps - each one centered at one of your small objects and used on it. They would move with objects if you plan to change their position. They would take into account your local geometry and, if you render it, your global skybox or skydome. You can precompute these cubemaps each time you change the position of your objects or do it every frame.

If you want to make your whole sponza reflective with high frequency, that's harder problem but one idea is to voxelize whole sponza with different mipmap voxel levels and then ray trace it as in voxel cone tracing technique. Voxelization could be done as precomputation but it is generally memory intensive.


In the paper you referenced $f$ refers to the fraction of water covered by foam that modifies optical properties on the water’s surface. It is explained on page 256 in: "Oceanic Whitecaps: Their Role in Air Sea Exchange Processes" (1986), by E. Monahan and G. MacNiocaill.

  How could I generate the waves for a water body with whitecaps & foam and variable intensity of the waves?

There are more techniques for rendering oceans than there are oceans to render. Some are based on speed (with reasonable results) and some are based on accuracy. On the fourth page of "Deep Water Animation and Rendering" is a discussion of rendering foam, with the prior pages covering other calculations.

Maya has a number of preset demos that can be altered or you can start from scratch. See this YouTube tutorial "Crafting the Wave: How to Render an Ocean in Maya".

Which of those references is more applicable depends on: if you simply want to understand more about it, if you are writing your own software, or if you want to render realistic oceans without necessarily knowing all the calculations behind it.