How to calculate camera distance so that after projective transformation a quad is filling the whole screen

Fortunately your scenario is rather simple, your camera is on a line directly perpendicular above and centered on the quad you're concerned about. So what you want is for the quad to fill the whole screen, I assume (more or less) exactly the whole screen, i.e. the size of your projected quad in window space is your entire viewport.

The mathematics of this are rather easy if you understand how the perspective matrix is actually constructed and what it does and if you know some basic trigonometry. The matrix constructed with glm::perspective takes as viewing angle the field of view along the y-axis and is based on goold old gluPerspective. And if we take a look at the actual matrix, we can see that it transforms a $1$ high object sitting on the view axis at distance $1$ from the camera to a height of $\cot \frac{\alpha}{2}$ ($\alpha$ being your view angle). If you move it further away it obviously gets smaller and if you increase its height, it gets larger, so it's $\frac{h}{d}\cot\frac{\alpha}{2}$ for camera distance $d$ and object height $h$. We want that height to be equal to $1$ (the height, from the center, of the view plane in NDC space).

So all you need to know in addition to that is the size of your square in the y-direction. We take half that height (because the square is centered) and multiply that by $\cot\frac{\alpha}{2}$ to get the distance it has to be from the camera. So if your square vertices have the y-coordinate +/-H and given your other values:

$$Y-X = H\cot22.5°$$
$$Y = X+H\cot22.5°$$
$$Y \approx X+2.41421H$$


Where is the shown display image actually stored and accessed?

For modern GPUs it's commonly stored in VRAM which for discrete GPUs are separate from the CPU memory. This usually means there is a penalty associated with letting the CPU access this as the data typically has to be copied across some bus. Mobile devices tend to use shared RAM for both GPU and CPU.


Optimal memory access when using lookup tables on GPU?

The best place to put a look up table for a GPU compute shader depends on the size of the lookup table, and the frequency/coherency of access. In your case (you mentioned 4kb), shared local memory would likely be best (assuming you do not need this memory for other purposes in the same kernel). This memory has different names in different APIs, but is the same architectural thing and follows the same performance guidelines:

CUDA: threadgroup shared memory
DirectCompute: groupshared memory
OpenCL: local memory
Metal: threadgroup memory
OpenGL: shared memory

Storing the lookup table in global memory as a read-only buffer may perform just as well, depending on the cache size(s) of the particular GPU you're running on.

Note that I'm presuming this is a read-only lookup table. A read-write lookup table is a completely different beast, and you don't have any good options there.


Why do I need memory barrier when doing atomic operations?

glMemoryBarrier is for when you want to make writes from one rendering command visible to reads from a subsequent rendering command.

That is, you can do all the atomics you like during a rendering command. But without the barrier, there's no guarantee that any of those atomic actions will be visible to later commands. Yes, they will still be atomic, but this only ensures that they will each execute in some order and each operation will fully execute before the next. What they won't ensure is that operations which OpenGL defines has happening "later" will execute before earlier operations.

Sometimes, you actually don't care about the order of the operations; being atomic may be enough. And sometimes, you do care.

Consider the case of a "linked-list" style blending buffer. Each pixel in the image is an integer that represents one of the items in the list. Each FS invocation will use an atomic increment to get a unique identifier, then swap that identifier with the current value in the image. It writes its pixel data to the entry in an SSBO, along with the identifier it swapped out of the image.

When doing that atomic increment on the unique identifier, you don't care if FS invocations from two different rendering commands execute out of order. You simply need the operation to be atomic. Similarly, the swap doesn't have to be in order; it just needs to be atomic. As long as they all execute, and the atomic operations do their jobs, then they will be fine.

Now, when it comes time to do blending based on this data, that is when you need to be able to read all of the data that was previously written. And that rendering operation requires glMemoryBarrier.

Of course, you're not necessarily using atomic operations to read these values.


Does smooth lighting work with Gouraud shading on single triangles?

After thinking about it for some days I came up with a proof sketch (which is hopefully correct) that it is possible to use Gouraud shading on disconnected triangles, if they share an edge with the same normals as in the picture. The simple idea of this proof is to check that there is no discontinuity between the intensity of the triangles by approaching from each side to the edge. So given two triangles $T_1: (P_1,P_2,P_3)$ and $T_2: (P_1', P_2', P_3')$ which share an edge and let's assume w.l.o.g. that this edges are constructed by $(P_2,P_3)$ and $(P_2', P_3')$. Because Gouraud shading only computes the light intensity at the vertex positions and interpolate in between them, it's save to assume that the intensity at $P_2$ is the same as on $P_2'$, same goes for $P_3$ and $P_3'$. Because it's also save to assume that the values inside the triangles are getting computed correctly so it's only necessary to show that there is no discontinuity between the triangles. To interpolate inside the triangle it's rational to use Barycentric coordinates so that the intensity at any point inside the triangle is given by
\begin{align*}
 I  &= aI_2   + bI_3   + (1 - (a+b)) I_1\\
 I' &= a'I_2' + b'I_3' + (1 - (a'+b')) I_1'
\end{align*}
By now going from $I_1$ respectively $I_1'$ to the edge, a sequence can be constructed to archive this, by setting $a_n = \frac{c}{n+1}, b_n=\frac{d}{n+1}, c + d = 1$ and analogously for $a_n', b_n'$. Then the intensity is
\begin{align*}
 I_n  &= \frac{c}{n+1} I_2  + \frac{d}{n+1} I_3  + (1 - \frac{c+d}{n+1}) I_1\\
 I_n' &= \frac{c'}{n+1}I_2' + \frac{d'}{n+1}I_3' + (1 - \frac{c'+d'}{n+1}) I_1'
\end{align*}
If one now approaches the edge with:
\begin{align*}
 &\lim\limits_{n\rightarrow 0^+} I_n = cI_2 + dI_3\\
 &\lim\limits_{n\rightarrow 0^+} I_n' = c'I_2 + d'I_3'
\end{align*}
and using the requirement that $I_n$ and $I_n'$ share the same position $P_n$ if they are on an edge, which can be constructed by
\begin{align*}
 P_n &= \frac{c}{n+1} P_2  + \frac{d}{n+1} P_3  + (1 - \frac{c+d}{n+1}) P_1\\
 P_n &= \frac{c'}{n+1}P_2' + \frac{d'}{n+1}P_3' + (1 - \frac{c'+d'}{n+1}) P_1'
\end{align*}
with the fact that $P_2 =P_2'$ and $P_3 = P_3'$ it follows:
\begin{align*}
 &\lim\limits_{n\rightarrow 0^+} P_n = cP_2 + dP_3\\
 &\lim\limits_{n\rightarrow 0^+} P_n' = c'P_2 + d'P_3
\end{align*}
and $d = 1 - c$
\begin{align*}
 &\lim\limits_{n\rightarrow 0^+} P_n = cP_2 + (1-c)P_3\\
 &\lim\limits_{n\rightarrow 0^+} P_n' = c'P_2 + (1-c')P_3
\end{align*}
this results in $c = c'$ and therefore this gives the final result:
\begin{align*}
 \lim\limits_{n\rightarrow 0^+} I_n = cI_2 + (1-c)I_3 = c'I_2' + (1-c')I_3' = \lim\limits_{n\rightarrow 0^+}I_n'
\end{align*}
Thereby is possible to use Gouraud shading over disconnected triangles as in my case, because the intensity on the edge is the same on both triangles.

What actually happened in my case was that the normals differs slightly (thank you @Rahul for the tip) which therefore results in this strange look. I fixed them by averaging the normals over a small radius to ensure
the proper set up of the normals.


Triangle - Why is it the basic building block for complex shapes

Triangles are always flat, unlike a quad which can be made in a way that is bent / not all vertices on the same plane.

If you were to interpolate say vertex normals across each scanline, on a 4-edge bent quad it will render with horizontal (scanline) artefacts. If you made that same bent quad out of triangles, it will render correctly as each triangle has a different normal slope.

That's one reason why triangles are preferable for scanline rasterising


Radiosity does not account for specular reflections (i.e. it only handles diffuse reflections). Whitted's ray-tracing only considers glossy or diffuse reflection, possibly mirror-reflected. And finally, Kajiya's path-tracing is the most general one [2], handling any number of diffuse, glossy and specular reflections.

So I think it depends on what you means by "ray-tracing": the technique developed by Whitted or any kind of "tracing rays"...

Side-note: Heckbert [1] (or Shirley?) devised a classification of light scattering events which took place as the light traveled from the luminaire to the eye. In general it has the following form:

L(S|D)*E

"L" stands for luminaire, "D" for diffuse reflection, "S" for specular reflection or refraction, "E" for eye, and the symbols "*", "|", "()", "[]" come from regular expressions notation and denote "zero or more", "or", "grouping", "one of", respectively. Veach [3] extended the notation in his famous dissertation by "D" for Lambertian, "S" for specular and "G" for glossy reflection, and "T" for transmission.

In particular, the following techniques are classified as:

OpenGL shading: EDL
Appel's ray-casting: E(D|G)L
Whitted's ray-tracing: E[S*](D|G)L
Kajiya's path-tracing: E[(D|G|S)+(D|G)]L
Golar's radiosity: ED*L

[1] Paul S. Heckbert. Adaptive radiosity textures for bidirectional ray tracing. SIGGRAPH Computer Graphics, Volume 24, Number 4, August 1990

[2] The Siggraph 2001 course "State of the Art in Monte Carlo Ray Tracing for Realistic Image Synthesis" says the following: "Distributed ray tracing and path tracing
includes multiple bounces involving non-specular scattering such as E(D|G)*L. However, even these methods ignore paths of the form E(D|G)S*L; that is, multiple specular bounces from the light source as in a caustic."

[3] Eric Veach. Robust Monte Carlo Methods for Light Transport Simulation. Ph.D. dissertation, Stanford University, December 1997


Join two bezier curves so that the result is two-times continuously differentiable

For 2 bezier curves to connect and be continuous at all, the last control point of the first curve (C in your case) must be the same as the first control point of the next curve (D). So c3 has to equal d0. That gives you C0 continuity or positional continuity.

The next step is to make sure that the tangents where the 2 curves meet are the same. That will ensure c1 continuity. (That is continuity of the first derivatives of the 2 curves.) To do this, the slope must be equal at both points. In order to do that, you need to ensure that the line from c2 to c3/d0 is continuous with the line from c3/d0 to d1.

I don't know of a good way to describe C2 continuity, but this Stack Overflow answer describes it better than I could. That answer should also explain the calculation you need for your proof.


Maybe something quick (to set up, no guarantees of actual speed) and dirty.

Whatever your method was, it probably boiled down to a vector function
\begin{align}
\begin{bmatrix} x,\, y,\, z\end{bmatrix} = F\Big([s,\, t,\, u] \,\, \Big| \,\, V_{ijk} \, : i,j,k = 1..N\Big)
\end{align}  where $V_{ijk} \in \mathbb{R}^3$are the vertices of the three dimensional lattice.

Basically, you currently have the new coordinates $[x,\, y,\, z]$ of a point, and you want to find the original ones $[s,\, t,\, u]$. I hope you still have the coordinates of the lattice vertices $V_{ijk}$. To retrieve the initial coordinates $[s,\, t,\, u]$ you have to solve for $[s,\, t,\, u]$ the system of equations
\begin{align}
\begin{bmatrix} x,\, y,\, z\end{bmatrix} = F\Big([s,\, t,\, u] \,\, \Big| \,\, V_{ijk} \, : i,j,k = 1..N\Big)
\end{align} where $[x,\, y,\, z]$ is given and the points $V_{ijk}$ are given
For simplicity, I am going to write
$$\begin{bmatrix} x,\, y,\, z\end{bmatrix} = F\big([s,\, t,\, u] \,\, \big| \,\, V \,\big)$$
To solve it, use a type of Newton's method.
Set a small enough number $h$.

Set your initial $s_0=x, t_0=y, u_0=z$ and set up your error equal to $err = 1$ (or whatever).

While $err \geq h^2$ execute the following iteration:

Assume you have a current approximation $[s_n,\, t_n,\, u_n]$;
Calculate the following three three-vector (as rows)
\begin{align}
&\Delta_s F_n = \frac{1}{h} \Big( \, F\big([s_n+h,\, t_n,\, u_n] \,\, \big| \,\, V \,\big) - F\big([s_n,\, t_n,\, u_n] \,\, \Big| \,\, V \,\big) \, \Big)\\
&\Delta_t F_n = \frac{1}{h} \Big( \, F\big([s_n,\, t_n+h,\, u_n] \,\, \big| \,\, V \,\big) - F\big([s_n,\, t_n,\, u_n] \,\, \Big| \,\, V \,\big) \, \Big)\\
&\Delta_u F_n = \frac{1}{h} \Big( \, F\big([s_n,\, t_n,\, u_n+h] \,\, \big| \,\, V \,\big) - F\big([s_n,\, t_n,\, u_n] \,\, \Big| \,\, V \,\big) \, \Big)
\end{align}
Set the following three by three matrix
$$M_n = \begin{bmatrix} \Delta_s F_n \\  \Delta_t F_n \\
   \Delta_u F_n\end{bmatrix}$$ (the vectors are the ones computed above but transposed and are therefore vector-columns);
Calculate the inverse matrix $M_n^{-1}$;
Calculate the update
$$[s_{n+1},\, t_{n+1},\, u_{n+1}] = [s_n,\, t_n,\, u_n] + \Big(\,[x,\, y,\, z] - F\big([s_n,\, t_n,\, u_n] \,\, \Big| \,\, V \,\big)\,\Big)\, M_n^{-1}$$
(all of this is written as vector matrix multiplication, where you are working with $1 \times 3$ row-vectors and multiplying such$1 \times 3$ row-vectors with the $3 \times 3$ matrix $M_n^{-1}$ on the right)
Calculate the error $$err = \Big\|\, [x,\, y,\, z] - F\big([s_{n+1},\, t_{n+1},\, u_{n+1}] \,\, \Big| \,\, V \,\big)\, \Big\|^2$$
Go back to the beginning of the iteration to check whether $err \geq h^2$

The result of this iteration is an approximation of the coordinates $[s,\, t,\, u]$. I hope this works.


Computing heat diffusion creates weird results

I messed up I am doing a double negative and that is causing the problem:
const float laplace_coeff = -0.5f * ((1.f / tan(alpha)) + (1.f / tan(beta)));
This should be positive.


This is perfectly possible

Although the difference may not especially noticeable, I would expect sampling taking into account the exact pixel geometry to give a slightly more accurate image. You just need to offset your pixel centres per colour component according to the (average) location of the subpixel(s) of that colour. Note that not all pixel layouts have a one to one correspondence between pixels and sub pixels.

For example, penTile RGBG has twice as many green subpixels as red and blue, as this image from Wikipedia shows:

I am not aware of any technical reason that would prevent this being used to produce arbitrary full colour images. In fact a colourful scene will have less noticeable colour artefacts than black on white text, which makes the colour differences hardest to camouflage.

Fonts are rendered on demand

The relevant difference between rendering a raytraced scene and rendering fonts is that fonts tend to be rendered on demand, and can take into account the screen being used. Contrasting with this, a raytraced scene is often prerendered and then displayed on many different types of screen (with different pixel geometry). For example, displaying your raytraced image on a webpage will prevent tailoring it to a specific monitor type.

If you were designing a realtime raytracing program and you had access to check the pixel geometry of the monitor, then you could raytrace to the specific sub pixel layout. However, offline raytracing that produces a still image can only be tailored to a single type of pixel geometry, which will then make the image look worse on any other pixel geometry. You could work around this by rendering a set of different images and choosing the appropriate one when it is later displayed on a particular type of monitor.

There is unlikely to be a long term benefit

So there is no reason you couldn't develop sub pixel rendering for a raytracer, but it means taking into account an aspect of the monitor that is not always known. Another thing to bear in mind is that you will be developing this software for a shrinking market. Sub pixel rendering is useful for screens that have relatively low resolution. As more and more screens (even mobile screens) are approaching such high resolution that the human eye cannot detect the difference made by sub pixel rendering, your work is likely to be more of theoretical interest than practical use.


What is the use of epsilon in the Möller–Trumbore intersection algorithm?

The epsilon is there to counter floating-point accuracy problems.  Note that the code checks that the algorithm is between $-\epsilon$ and $\epsilon$; that is, close to zero.  If the determinant is zero, then the ray lies in the plane of the triangle and can be safely said to not intersect it (as the alternative, an edge-on intersection, should be caught by neighbouring triangles in the mesh).

The Cramer's rule computation of the intersection point then divides by the determinant, and here's where the floating-point accuracy comes in.  If the determinant is merely very close to zero, not exactly equal to it, then its reciprocal will be very large.  We multiply the other determinants by this reciprocal, so small floating-point errors in their values will be blown up into significance and make the difference between an intersection and a non-intersection.

I haven't tested it, but I suspect that leaving out the epsilon will cause intersections on the silhouette of a mesh (i.e. with triangles close to parallel with the ray direction) to become very nondeterministic.  This would result in nearby rays either missing or hitting the triangles unpredictably and giving you speckling.

The value of the epsilon depends on the floating-point precision you're using. It should be small enough that only very grazing intersections are culled, but large enough that it avoids the speckling problem.  It can be determined by experiment, or by just eyeballing it.  The implementation you linked to uses $\epsilon = 10^{-7}$, which is a reasonable choice for single-precision arithmetic; you may want to use $10^{-9}$ or $10^{-10}$ for double-precision computations.


Using original vertices coordinates after clipping

A vertex consists of not just the position, but all of the values passed from the vertex processing stage to the rasterizer. So if a vertex gets clipped, the new vertices generated in that process must generate new values for those vertex processing outputs. It does this using the same math that it uses to generate the new vertex position.


Why is the color in the cube being weirdly swapped?

It is likely that your specific problem arises from a singularity in atan when m_position.z goes to 0. Rather than try to debug that, I'll point out that you don't need to transform your camera position into an angle to then modify it: rather, just create a rotation from your newPhi and newTheta angles, in the form of a rotation matrix or quaternion, and then apply it to your vector, directly rotating it without doing the perilous round-trip to an angle.


Deciphering Affine/Projective Transformation Code

Your understanding of the matrix structure in Q3 is correct. This code just does not construct a matrix explicitly and the matrix multiplication is applied implicitly. I think this part might cause your confusion.
Instead of deciphering the code, I would rather derive the transform and compare it with the code.

The affine (6 degrees of freedom) and projective matrix (8 dof) map points to points. So you can build a system of equations for each pair of points in different frames. You will need to solve this system for the components of the matrix (Q1).

Let $M$ be the transform matrix, then

$$
\begin{align}
M (0,0,1) &= (X_1,Y_1,1) \\
M (1,0,1) &= (X_2,Y_2,1) \\
M (0,1,1) &= (X_3,Y_3,1) \\
\end{align}
$$
Each of above lines gives you two equations. Now you have 6 equations and the affine matrix can be derived.
For the projective matrix, we need extra two equations using the point $(X_4,Y_4)$:

$$
\begin{align}
H_4 &= M (1,1,1) \\
\frac{H_4}{H_4.z} &= (X_4,Y_4,1)
\end{align}
$$
Note that I have to divide the homogeneous part for the transformed point.

Solve this linear system you can get an expression for your matrix $M$. Perform $M (f_x,f_y,1)$ you will see how it transforms the point. I believe this solves your Q2, Q3, Q5 and Q7.

For Q4, I think it is just a convention of pixel position and your understanding should be correct.

Appendix: maple code to verify the answer

with(ListTools):
# affine matrix
m:=<m11,m21,0|m12,m22,0|m13,m23,1>;
# Build equations
E:=[convert(m.<0,0,1> - <X1,Y1,1>,list),convert(m.<1,0,1> - <X2,Y2,1>,list), convert(m.<0,1,1> - <X3,Y3,1>,list)];
# E=[[m13 - X1, m23 - Y1, 0], [m11 + m13 - X2, m21 + m23 - Y2, 0],[m12 + m13 - X3, m22 + m23 - Y3, 0]]
# Solve E=0
sol:=eliminate(Flatten(E),indets(m));
# Plug the solution back to m:
subs(sol[1],m);
                           [-X1 + X2    -X1 + X3    X1]
                           [                          ]
                           [-Y1 + Y2    -Y1 + Y3    Y1]
                           [                          ]
                           [   0           0        1 ]

# projective matrix
p:=<m11,m21,a|m12,m22,b|m13,m23,1>;
# Compute transformed points
U0:=p.<0,0,1>;
U1:=p.<1,0,1>;
U2:=p.<1,1,1>;
U2:=p.<0,1,1>;
U3:=p.<1,1,1>;
# Build equations
E:=[convert(U0/U0[3]-<X1,Y1,1>,list),convert(U1/U1[3]-<X2,Y2,1>,list),convert(U2/U2[3]-<X3,Y3,1>,list),convert(U3/U3[3]-<X4,Y4,1>,list)];
# Sove E=0
solp:=eliminate(Flatten(E),indets(p));
# copy their code for xx and yy:
xx:=((X4-X1)*(Y3-Y1)-(Y4-Y1)*(X3-X1))/((X2-X1)*(Y3-Y1)-(Y2-Y1)*(X3-X1));
yy:=((X4-X1)*(Y2-Y1)-(Y4-Y1)*(X2-X1))/((X3-X1)*(Y2-Y1)-(Y3-Y1)*(X2-X1));
# Verify our a and b, compare ours with their expressing using xx and yy
subs(solp[1],a);
subs(solp[1],b);
simplify(subs(solp[1],a)-(yy-1)/(1-xx-yy));
# It is 0 as expected
simplify(subs(solp[1],b)-(xx-1)/(1-xx-yy));
# It is also 0


Energy redistribution sample re-integration

It sounds like you yet have to adjust the importance function of the Markov Chain according to the way you assign paths to each of the two estimators, path tracing and ERPT.
If you think about it in the framework of MIS, you have two estimators with binary weights (according to your threshold). Now, the simplest (but inefficient) option to combine them would be to always run ERPT and path tracing for all paths, applying the weight to the splats only (this weighting is what I understand you do by separating them into different framebuffers).
In your case, however, you only enter the ERPT chain for the paths that you want to further explore, therefore you do not actually enter the chain in a stationary distribution proportional to the measurement contribution. This breaks ERPT, as the startup bias removal no longer works correctly. To fix it, you need to additionally account for your path selection strategy in the importance function such that it matches the initial distribution of paths entering ERPT: The MIS weight of the ERPT estimator needs to become a part of its targeted distribution, preventing energy from running off into paths that would never actually trigger an ERPT run. Thus, you should be able to make the startup bias removal work with your path selection strategy.


My solution to this ended up being not sorting at all. Manipualting an already initialized HE is a nightmare, so instead I built the HE with this property from the get go.
The initialization algorithm is pretty much:

Create a new half edge.
Check if the pair of the half edge exists.
If it exists, you are allocated at index of pair + 1
If it does not exist, increment an offset by 2 and store yourself at offset - 2 (i.e. allocate 2 slots in the edge array and store yourself in the even slot).

Doing it this way before setting next and prev and company guarantees that at any one time the HE is properly initialized.


How can I draw a filled circle on a grid with the less rectangle as possible?

To get a proper circle, you can use a transformation matrix to transform your cell coordinates into equally spaced cartesian coordinates. In the case of your 50 x 100 pixel cells, your transformation matrix would probably either double the width or halve the height. So something like this:

| 0.5 0.0 |
| 0.0 1.0 |

If you want to draw a circle of radius r at a center point (cx, cy), you could iterate over the cells from (cx - r, cy - r) to (cx + r, cy + r) and test if each one is in the circle, and if so, fill it in, otherwise, don't. Something like this pseudocode:

double rSquared = r * r;
for (int row = cy - r; row < cy + r; row++)
{
    for (int col = cx - r; col < cx + r; col++)
    {
        double x = col * transform[0][0] + row * transform[0][1];
        double y = col * transform[1][0] + row * transform[1][1];
        double distSquared = x * x + y * y;
        if (distSquared <= rSquared)
        {
            // inside the circle, fill in this cell (col, row)
        }
        else
        {
            // outside the circle, don't fill in this cell (col, row)
        }
    }
}

As for drawing rectangles, it's not clear what you're trying to accomplish there. There's no need to use rectangles when drawing a circle.


How does a GPU process a task by using multiple ALUs?

Architecture Background

GPU don't schedule to individual ALUs, but rather groups of them. Terminology will vary based on GPU company (NVidia, ATI, and Intel), but I will stick with NVidia terminology since they all work somewhat similarly at the high level.

NVidia has a concept of CUDA cores, each of which contains an ALU and/or FPU, which make up collections called Streaming Multiprocessors (SM). The Fermi architecture has 32 CUDA cores per SM, for example.

GPUs operate in a SIMD-like manner. This means that unlike a CPU, which traditionally has had one execution unit per instruction, GPUs have many execution units per instruction (CUDA core). These groups of execution are called warps.

Answer

GPUs contain scheduling units that attempt to reach maximal utilization. Warps are scheduled in a way to utilize as much of the GPU as possible. One way to do this is to use a round-robin approach. Latency is something that a GPU struggles with compared to a CPU, so scheduling a process as soon as possible is often desirable. Like an OS (which controls CPU scheduling), however, simplicity and hiding latency are still desirable traits to an extent, although other tradeoffs may be deemed more important.

Your workload will determine how scheduling is done. If you run a CUDA kernel or shader on 320 pieces of data (threads), you will need at least 10 warps to be scheduled (this depends on the type of data as well). If you have 1 thread that needs to be calculated, oh well, you still need to schedule one warp to do the calculation.

Of course, what ratchet freak said is true about GPU architectures and algorithms being proprietary, especially when compared to CPUs.


I know im a bit late, but hopefully this can help other people. you can add detail by using Fractal Brownian Motion. There is a great article here that i used as a guide to make my own version in c#. Also I would check to make sure your vars are ints.


What is the basic difference between Radar Screen, CRT and Oscilloscope displays?

At many places in CG hardwares and algorithms you have the choice between vector description and raster description. Oscilloscops, radar, lasers shows, plotting tables, laser cuts directly follow lines and curves; these are the vectorials devices. LCD, inkjet and laser printers, stereolithographic 3D printing, follow pixels and voxels; these are the raster devices. The CRT, but also some laser shows, use the oscilloscope technology to implement raster: the "curve" is a repeated Z narrowly sweeping the screen line after line, with intensity modulation making the "pixels". Radars are indeed between this and pure vectors since rotation angle is swept. Details about various forms of radar display and oscilloscope is available on the wikipedia page.

( About representations, the vector/raster dichotomy corresponds to inkscape/illustrator vs gimp/photoshop, 3D meshes vs 3D voxels, list of lobes vs angular or SH decomposition, Lagrangian animation vs Eulerian animation, etc. )


Jitter is supposed to make edges more irregular; that's how it masks aliasing in the first place. I didn't realise this when I made the question, but the noise in path tracers comes from undersampled global illumination; because each sample travels in a random direction, it's possible for a visually-significant number of samples to reach unrepresentative areas (e.g. for rays at the bottom of a curve to mostly hit the walls of the curve rather than travel upward towards), and the dissonance between these and the actual image function creates the fuzzy variance you get in path-tracing renderers.

That sort of noise is completely different to the pseudo-one-dimensional edge variance you expect to get from ray jitter at one sample/pixel. Single-sample ray jitter essentially means pulling edges randomly closer to or further away from their actual tendency in order to hide the high-frequency patterns introduced by aliasing; it looks nothing like "fuzzy" GI noise because (at least at the single-sample scale I was thinking about when I made this question) it's one-dimensional (variance about an edge) rather than two-dimensional (variance over a surface).

So, this question is invalid. It assumes that jittered edges are supposed to be fuzzy (as in GI noise), but forgets that edges are one-dimensional and that the jitter can't magically expand them to two dimensions and spread them over some domain (at least, not without spreading the aliasing as well (multi-pixel jitter)). Jitter doesn't produce fuzzy noise because it isn't supposed to; it's only supposed to vary the relative distance of points "on" an edge from the actual edge positions given by the image signal.


It's not 100% clear what the author means here, but I'll choose to interpret "screen coordinates" as "pixel coordinates". These would be related to the projected points by a 2D coordinate transformation.

You're correct that projection is done by dividing each component of a point by its $z$ component. (That's not quite true—actually, we usually use homogeneous coordinates, which have four components $xyzw$, and divide by the $w$ component, which has previously been set to eye-space $z$ by the projection matrix. This frees up the $z$ component of the vectors to be used for the z-buffer.) The projected points $q_i$ are then in "normalized device coordinates" or NDC space, which has a range of [−1, 1] on each axis.

To calculate the resulting pixel coordinates on the screen, we then need to perform a 2D scale-and-translate operation, to map the [−1, 1] range to the desired range of pixels. For example, if we want a 640×480 viewport, we could calculate
$$x_i = 640 \cdot \frac{q^x_i + 1}{2} \\ y_i = 480 \cdot \frac{q^y_i + 1}{2}$$
You might also need to flip the $y$ axis in this transformation, for instance if the NDC space is Y-up but the pixel coordinates are Y-down.


An algorithms for covering a 2d polygon with a predetermined number of rectangles?

(Not an answer but my comments are too large for a "comment")

This is more computational geometry than graphics. First of all, it would help if you could give an example (a drawing) which describes your input and your expected output. Are you working with low count polygons, or precise coastal maps? Do you want a real-time solution or is it off-line? It could be a big help if you would describe what you want to obtain with these rectangles covering your polygon. Also, what are your needs concerning the solution(s), for instance do you want the area of the rectangles to be as least as possible?

Intrinsically this is NP-hard, this means there is a vast amount of possible solutions when the size of input data rises. Think of a circle covered by a square divided in two not necessarily equally sized rectangles. There is an infinite set of solutions already there concerning where they are split. Then just think of the amount of possible solutions for optimally packing a star shape into five or so rectangles, it's simply "not possible" if you want it done on a computer, in finite time.

So, you will be looking for a heuristic, that gives you a solution but not necessarily an optimal one. This means that you need to figure out what you want to accept being a solution.

Edit: One variant of the problem at hand, that I didn't think of, could be that you already have N rectangles (size is given up front, they might be of the same size) that you wish to cover a polygon with. It's not simpler, also it could be impossible, for instance if the total area of the rectangles is less than the area of the polygon. So: Please describe your problem at bit more :)


Modelling Young's double slit experiment

It is in fact possible to augment a ray tracer to make it capable of simulating wave effects. The paper Rendering Wave Effects with Augmented Light Field
 describes one way to do it. In summary, they introduce a framework called Augmented Light Field that allows them to model wave effects with a ray-based representation. In this framework, rays can carry negative radiance in addition to positive radiance; intuitively, rays carrying negative radiance can "subtract" light from surfaces that it cannot reach because of interference effects, for example.

I haven't implemented it personally, and I can't speak for its implementation complexity, but I believe it is a good starting point.


How does a GLSL Sampler2D read a Vec2 UV coordinate to determine UV range

It's hard to say without seeing the shader in question. My guess is that every animation frame in the atlas has the same size. This fixed size allows mapping any UV coordinates to the proper animation frame's UV extents.

For example, with 16 frames in a texture:

0 0.25 0.5 0.75 1.0
\  \   |   /  _/
 +--+--+--+--+
 | 0| 1| 2| 3|
 +--+--+--+--+ 0.25
 | 4| 5| 6| 7|
 +--+--+--+--+ 0.5
 | 8| 9|10|11|
 +--+--+--+--+ 0.75
 |12|13|14|15|
 +--+--+--+--+ 1.0

This code should do the trick to convert UV coordinates (range [0, 1]) to a given frame's UV extents:

vec2 uv = ...
uint frameId = ... // [0, 15]
const uint rowFrameCount = 4; // implies frameUVWidth = 1 / 4 = 0.25
const vec2 frameUVSize = vec2(1.0 / rowFrameCount, 0.25);

// Find the column and row offset
uvec2 frameOffset = uvec2(frameId % rowFrameCount, frameId / rowFrameCount);

// Offset to the proper column
uv += frameOffset;

// Scale the coordinates in the range of a single frame (here [f_0, f_0 + 0.25])
uv *= frameUVSize;

Some results:

uv = vec2(0, 1); frameId = 0; uv_out = vec2(0, 0.25);
uv = vec2(0, 0.5); frameId = 3; uv_out = vec2(0.75, 0.125);
uv = vec2(0.5, 0.5); frameId = 14; uv_out = vec2(0.625, 0.875);

In truth however, this kind of texture trickery is less useful these days since array textures are now widely available with better results and simpler indexing.


Camera Model View transformation in Ray Tracing

The space at which you transform your vertices is completely up to you, because it depends on what algorithms and kind of effects that you are trying to achieve.

As of my personal experience, I usually shoot rays in world space because eventually we all need some sort of "world-space" acceleration data structure, such as a space-partition tree, that gathers all objects into a single space. As this data structure is probably going to be pre-computed for static objects, world space seems to be the most preferred one, as the vertices representation never change when you move the camera around.

When it comes the time to actually test the ray-object intersection, I transform the ray direction onto object space by using the inverse of the model matrix. This makes it possible for all intersection tests to be implemented as if the object were axis-aligned. Thus, you can easily use AABBs (axis-aligned bounding boxes) for coarse intersection testing, even when your models are rotated in world space. This is also a good optimisation for dynamic objects, as the AABB does not change in object space.


How to count aspect ratio of an image when you have dimensions?

You need to take your width and height and make them into a fraction, then reduce it. To do that, you'll need to find the greatest common divisor of the numerator and denominator, and divide them both by that.

Now take the new fraction you got by dividing both the numerator and denominator by the greatest common divisor and do it again. Repeat this process until the greatest common divisor is 1. At that point you have your irreducible fraction and you can print it.

Another way to do it is to use Euclid's Algorithm.

Example:

Let's look at your top example, 1,000 x 1,200

1,000 / 1,200 = ?

We can see by examination that they both end in "00" so they're both divisible by 100. So start there:

1,000 ÷ 100 = 10
1,200 ÷ 100 = 12

Next, we see that 10 and 12 are both even, so they're both divisible by 2:

10 ÷ 2 = 5
12 ÷ 2 = 6

Finally, we see that 5 and 6 share no common factors besides 1. (This is because 5 is prime, so its only factors are 5 and 1. 6 has the factor 2 and 3, as well as 1 and 6.) So we have reached the end. The reduced fraction is 5 / 6.

One More Thing

You might also want to be able to detect cases of irrational ratios, such as the golden ratio, which is quite popular in visual arts. Such a ratio will not be possible to divide into small integers, so checking whether width / height is within some small distance of the golden ratio may be sufficient.


Why does skinny triangle is avoided in triangulation algorithm?

We don't want to avoid skinny triangles as much as we want to avoid really long triangles.

Long triangles make interpolation of properties less accurate because the points in which the data is stored are farther apart. This is equivalent to – when using a 2160p (4K UHD) screen – the preference of watching a 1080p video in fullscreen mode to watching the same video but in 720p, also in fullscreen mode. The 720p video is streched 3 times when the image is rescaled to fullscreen while the 1080p video is only streched 2 times; hence the pixels in the 720p video are streched farther apart compared with the pixels in the 1080p video and the interpolation of the colors in the pixels become less acurate (simply because the 720p video can't contain as fine details as the 1080p video can), which is another way of saying that the 1080p video looks better than the 720p video when viewed in fullscreen on a 2160p screen.

However, the skinnier our triangles are, the longer they have to be in order to cover the same area, unless we use more triangles in total. The former is undesirable because of above mentioned reason. The latter is usually undesirable too since it increases the required memory consumption as well as the computational load on the algorithms that work with them since there are more triangles to loop through. We usually therefore want to avoid skinny triangles.


Euler Angles - Gimbal lock, why non-orthogonal axes

OK I think I know what your problem is, gimbal lock does not really lock anything from a mathematical point if view, only for certain operations. See euler suggested that euler angles can define all space angles. And that is true. Its just that euler does not say that the space is uniform. Highly mathematical people also say that the euler space has impossible solutions since some positions are overlapping, but in practice this is not really a problem.

Now since you can reach any angle, and you can always unwind your angle you can reach any angle in the body. Provided that your ready to unravel the situation and start over for a computer program doing animation and a mathematician this seems reasonable.

However, this is not what most people are doing. Instead if you try to do some physical things like calculating the distance between 2 positions, the energy needed to rotate, use a external force to move the system or make a delta move in the gimbals space then this non uniformity kicks in. Simply the space is not uniform the same way that a real space is at the poles the distance to a position next to you in real space can be surprisingly far away in Euler rotation space.

So when somebody is euler locked it does not mean they can not move anywhere they just can not move in the direction they would want to move in that particular orientation of the gimbal, and close to it. Its a bit like living in the suburbs and wanting to reach your neighbors house (that is 100 m away) but being forced to go trough the center of the city (15 km away). Energetically or time wise highly wasteful.

This is a problem in 3D graphics only if you choose to interpolate in Euler angles. But for a control engineer or a mechanical engineer this can be  a nightmare, since you can suddenly have a discrepancy between what your mechanism and your control logic is saying.

Now it is possible to construct the gimbal so that the pole does not really happen, just not with the orthogonal axes the space is still not uniform but it does not experience the extreme case. However it might just be easier to use some other formulation like quaternionions, or axis angle instead. (quaternionions aren't really practical if your building a real mechanism though, its possible but very convoluted)

PS: I would draw an image, its just that drawing something the wraps perfectly upon itself and is a 3d field is a bit difficult to do.


If you really need this level of debugging (spoiler: you don't, more on that later), AFAIK only NVIDIA's NSight can provide it. However, in order to remain interactive, you probably need two machines - one to actually run your shader, and the other to debug it - because stopping a GPU for debugging will - quite understandably - freeze the display of the computer, making it unusable.

One step down from that is RenderDoc and its ability to swap shaders in a capture on the fly, as well as step through Direct3D bytecode and inspect its registers (i.e. a virtual machine for the DXBC, it's most often sufficiently detailed to actually find shader bugs).

And what most graphics programmers do most days is simply what I like to call "your mom's debugger", which on the CPU means sprinkling your code with print statements with the intention of later reviewing the log; but in shaders, you just stuff your shader outputs with debug values instead and inspect those values. There are also ways to "print text" from shaders (for instance, this ShaderToy).


Regular voxel grids and RTX accelerated ray tracing

As you said, the RTX Turing architecture comes with wired primitive-ray intersection, (to be more specific, triangle-ray intersection). The BVH is built by specifying the Bounding Box program to OptiX, the signature of which is:

RT_PROGRAM void my_boundingbox_build_program( int, float result[6] )

As you can guess, the result must contain the minimun (3 first components) and maximum (3 last components) bounds of the bounding box. You can use this to build virtual bounding boxes even if you have no geometry (for instance, a chunk of 16x16x16 voxels).

When you perform rtTrace , the OptiX API will detect the collision through the BVH which it has build with the help of your bounding box program. In that momment, it will call the Intersection program, which you can specify and which signature is:

RT_PROGRAM void my_intersection_program( int prim_index )

Where you will report any possible intersection within your 16x16x16 chunk of voxels (obviously you can ignore the prim_index, which is only useful when using geometry primitives).

Now, if you would use their special structure RTgeometrytriangles , the Turing RTX architecture would be able to make use of the wired primitive intersection, which is where performance boost would be gained, and you would not need to implement your own intersection program.

To sumarize:
Since you are not using the specific RTgeometrytriangles , you will not access the intersection performance boost from the Turing architecture.


One method is to do projection mapping. I understand that most 3D rendering applications can do that. I don't have much experience with them, but I have done it using AfterEffects using the Camera Mapper in Buena Depth Cue.

Algorithmically, it's achieved by projecting the geometry back onto the image to generate texture coordinates. So one way to do that is as follows:

Line up the image of the scene in front of the corresponding geometry
For each triangle in your geometry, project the 3 corners onto the image
The coordinates where the projection hits your image are the texture coordinates you need for that triangle


You are on the right track but what you need to do is to calculate u/w and v/w, and also 1/w for each vertex, which you interpolate linearly in screen space in your rasterizer. Then for every pixel you divide the interpolated u/w and v/w coordinates with the interpolated 1/w to get perspective correct uv-coordinates for the pixel.

The same applies to all the vertex attributes, e.g. if you need to get perspective correct color or position of vertices per pixel, then you interpolate c/w and p/w and divide those values by 1/w per pixel. Even if you don't implement software rasterizer this can be quite useful for example in ray marching of screen-space reflection implementation. I have seen people doing the ray marching for SSR in world space because they don't know how to properly interpolate position in screen space.

What old software rasterizers used to do is to perform this 1/(1/w) division only every X pixels and linearly interpolate values in-between because the division was relatively expensive operation to perform every pixel.


Are textures automatically Morton-coded when loaded to VRAM?

In vulkan, like a lot of other implicit operations, this transition is made explicit. Images can be in one of 2 tiling layouts: optimal and linear. Linear images have a fixed layout while optimal is with implementation-defined space-filling curve. (There is no way currently to figure out which exact curve is used).

To fill a optimally tiled image you need to copy from another image or from a vkBuffer (which has the data in linear format). The swizzle is done as part of the copy. But with the requirement that the data source must not be changed during the operation. Because other operations are allowed be read from the copy operation's source.


3D scene construction and change of coordinates

I assume the term object description refers to the set of 3D vertices a mesh is made of. And I also assume reference frame actually means coordinate space. For a book on computer graphics, it sure uses some very uncommon words.

In computer graphics, we often have to transform a set of vertices from one coordinate space to another. For example, an artist may have modelled a rubber ducky, aligned in such a way that it looks down the positive x-axis and it uses centimeters for its units.

When we place the ducky in our world (let's say it uses meters as the primary unit), we scale, rotate and translate it so that it neatly floats on the water in a bath tub.

Then we have our camera, which can move freely around in our scene with its own translation and rotation parameters. And finally there is the notion of projection, which takes into account the field of view and the aspect ratio of our frame to project the 3D point onto a 2D image.

Our job is to calculate exactly where in our 2D image the 3D vertex ends up, if at all. And for this, we will have to go through the sequence of all our coordinate spaces: from object to world to camera (a.k.a. view) to projection.

So that's how most graphics programmers visualise this process and hence why the book also follows this sequence: it takes the object and follows it through its transformations until it ends up in our 2D image.

It seems you prefer it the other way around: you take the camera as a fixed object and rotate the world around it. But you have to admit this isn't very logical, as the camera sure isn't fixed and the world does not revolve around it.


How to translate mouse movement to a rotation gizmo?

Map the mouse position to the circle by drawing an imaginary line between the mouse and the center of the gizmo on the plane of the current circle. Do this for the start and current positions.

Then you can easily get the quaternion rotation around the center from the start to the end by doing normalize(quat(1+dot(v_1, v_2), cross(v_1, v_2))).

This works because the dot product is the cos of the angle between the vectors and the cross is the vector perpendicular to the 2 vectors (aka the rotation angle) and has length of the sin of the angle between the vectors. This means that you can create the double quaternion out of them. Then you can nlerp that double quaternion with the unit quaternion to get the rotation you need.

However with that information you can do an arcball by projecting the mouse onto the ball that is the gizmo and using the start and end vectors that way.


I'm not sure what happens there. Looking at the rendering of "considerations" in your example, the "cons" looks like what you get when you do crude sub pixel anti-aliasing, without applying the necessary "fix up" to make it look good. But then the "der" in the same word looks like there's no sub-pixel anti-aliasing.

You can find a good article about how good looking sub-pixel anti-aliasing works here https://www.grc.com/ct/ctwhat.htm and for the "fix up" part here https://www.grc.com/ct/freeandclear.htm and here https://www.grc.com/ct/cttech.htm

What it comes down to is: if you just compute the R/G/B intensities by calculating the covered area for each sub-pixel, the resulting color fringes are very ugly/irritating. To fix it you have to apply a horizontal blur/low-pass filter to the sub-pixel image. The result is an image that's less sharp but also has far less saturated colors. I.e. something like your "Check you" example.


"Importance Sampling of Many Lights with Adaptive Tree Splitting" - paper: participating media

Someone else feel free to correct me as I am not 100% certain, but from what I can gather:
They are describing all the normalized vectors v, which point at the ray between points v1 and v2.
They are doing this using one parameter angle φ.
They are essentially interpolating between two vectors o1 and o2, with the equation v = o1 cosφ + o2 sinφ
o1 is simply the normalized v1. o2 is a vector at right angles with the o1 along the plane of the ray. (I think this would be calculated by cross product of o1 and N from the paper).
Because this angle can sweep greater than the length of the ray defined by v1 and v2 they also define a limit Bmax so that φ always stays within the bounds of the ray.
I've tried to annotate the diagram below (but I think I made a mistake and o2 would have an angle create than v1).

They also discuss it briefly here: http://www.aconty.com/pdf/importance-sampling-lights-slides.pdf
Might be worth mesaging the orinal authors to see if this is correct.


How can I detect edges between different colours of the same brightness?

I think you can convert from the RGB space to the HSV one, or whatever color space has the HUE in a single channel.

Take the HUE channel, and make the edge detection on that one.

Here a simple Matlab script to achieve the result.

I = imread('image.png');
hsv = rgb2hsv(I);
hue = hsv(:,:,1);
edges = edge(uint8(hue),'sobel');
imshow(edges);


OK, I have figured out this problem.

For the texture coordinates, there are two kind of texture coordinates. One is normalized texture coordinates, which is in the range [0, 1]. The other is texel space, which is in the range [0, size), where size is the size of the texture. For texelFetch(), the texel space is used.

For lod, it means the level of detail in the mipmap. We can simply use 0 for the base level (the original size).

Thus, an example of using texelFetch() should be:

vec4 texelValue= texelFetch(texture2D, ivec2(gl_FragCoord.xy), 0);


The issue is that on your system, writing out "\n" is writing out 0x0D0A, whereas on my system, it's only writing out 0x0A. If you explicitly write out 0x0A instead of "\n" I think it will work for you.

Here's a hex dump on macOS Mojave:

50 36 0A 35
31 32 20 35
31 32 0A 32
35 35 0A
FF 00 00 // <- pixels start here
FF 00 00
FF 00 00
FF 00 00
FF 00 00
FF 00 00


How do I use barycentric coordinates to interpolate vertex normal?

Thats the way your defining the vectors. See your defining vectors $V_1$ and $V_2$ as pointing outwards from the point $P_2$.

Now then $denom$ is the square of the area of the parallelogram of formed by completing the triangle with a duplicate of itself. Another way to look at this is that its the 2 tines the square of the are of the triangle.

Now then the calculation:

(d11 * d20 - d01 * d21)

Is the square area of the parallelogram of $P_2$-$T$-$P_3$

Now if you divide this with denom, then you have calculated the total area relative area of triangle $P_2$-$T$-$P_3$. So the formula (d11 * d20 - d01 * d21) / denom is the triangle size in relation to whole triangle size of:

This is the weight that the opposing corner exerts on the point. which in other words it is the barycenter weight related to point $P_1$.

The next line is simply the relative area of triangle $P_2$-$T$-$P_1$ or:

Which is opposite to $P_3$ being the barycentric coordinate of $P_3$.

Lastly you calculate the rest which is the remaining triangle which is opposite of $P_2$ using 1 - bary[0] - bary[1].

Note: You could have chosen the first vectors differently for a different order.


What is the algorithm used in pbrt-v3 for triangle intersection?

The permutation is described in section 3.6.2, page 158-159.

A reason behind this permutation is that a shear matrix (which involves dividing z component of ray direction) is computed to transform from world space to ray triangle intersection coordinate space.
Permuting the coordinates is to avoid divide-by-0 error by mapping nonzero component to +z.


How do I handle uniforms in a deferred renderer?

I'd say that the fundamental problem here is your approach. You basically want the outside renderer to throw a bunch of random stuff at the wall, then let the internal rendering system sort it out without having any idea what any of that stuff actually means.

Here's what I mean. Consider the case of skinned meshes vs. unskinned meshes. The internal rendering system has no idea which ones are skinned or not; all it knows is that some draw commands use a particular program and a large UBO/SSBO, while others use a different program and much smaller UBOs.

And that kinda works... to a point. And that point is when you're dealing with state values that are no longer encapsulated in objects.

Uniform values are program state. And you might reasonably have two different objects that use the same program have the same value for that uniform. They might have different values for other uniforms, but there would be different variance in which objects use which values.

But here's the thing: in pretty much every case where this happens, someone knows that it happened. And the code that knows that is the higher level code. After all, if two objects have common programs and common uniform values, they probably have a bunch of other stuff in common too. So odds are good that the higher level code knows when such things would happen.

By using this "throw stuff against the wall" approach, by separating the business end of your renderer from the source of that rendering data, you destroy the semantic information you could have used to detect such things. That is, it is the higher level system that should handle the "sort by state change" part, since it can do a much better job.

However, given that you've effectively designed yourself into a corner, I'd say that the simplest approach is to develop a way to encode uniform values, rather than relying on callbacks.

Uniforms can be scalars, vectors, or matrices (arrays and structs use multiple uniform locations). So you need the equivalent of a "variant" type that can store all of such data, of the types supported by GLSL. And then, for each object, you walk through each uniform location in that object's data block and set that value. Sure, "visitation" is not the most instruction-cache friendly approach, but it's a lot better performing than calling a std::function. At least with visitation, you're executing the same block of code each time.

So what you'd have is a (possibly dynamic) array of uniform data, where each uniform datum is a location, a type, and the variant data. I'll leave the storage optimization of the "variant" as an exercise for the reader.

Generally speaking however, non-block uniform usage should be fairly minimal, dealing only with values that change frequently while rendering with a program. So even in that case, things shouldn't be too bad.


Why is the symbol for solid angle a small omega in the definition of the BRDF?

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


Actual vs Perceived Brightness of RGB Colour

Two different effects are causing the observation mentioned in the video.

On one side, the vast majority of screens have a non linear response: if the RGB value is half as much, the emitted light intensity is not half as much. This behavior originally comes from cathode ray tube (CRT) displays, which produced different light intensities by varying the voltage used to generate the electron beam. From the article:

  CRTs have a pronounced triode characteristic, which results in significant gamma (a nonlinear relationship in an electron gun between applied video voltage and beam intensity).

As CRT became ubiquitous, systems assumed this behavior to be present and took it into account. Nowadays CRT have mostly been replaced by LCD technology, which emulate this behavior for compatibility.

Gamma is usually considered to be around $2.2$, which means the light intensity follows a function $I=x^{2.2}$. This is where the 22% mentioned in the video comes from: ($0.5^{2.2} \approx 0.218$).

On the other side, the human vision is not linear either. According to Wikipedia, a typical sunset on a clear day is about 400 lux, and clear sky daylight is about 10000 to 25000 lux. Yet the latter doesn't feel over 25 times as bright as the former.

In the end, the non linear response of screens is a good thing, because it gives more precision where it is needed. Note that it is not logarithmic though, but closer to a square function. If you want to read more on the topic, you can look for articles on "gamma correction".


Is there a way to perform calculations of Mandelbrot set using only integer numbers?

Probably the easiest way to get an arbitrary precision Mandelbrot set using full-precision integers is to combine two integers into a rational number $p/q$. These numbers are added by finding the common denominator:
$$ \frac{m}{n} + \frac{p}{q} = \frac{mq + np}{nq} $$
and multiplied by just multiplying the respective numerators and denominators:
$$ \frac{m}{n} * \frac{p}{q} = \frac{mp}{nq} $$
Then a single point $\vec{p}$ would be represented in rational coordinates by
$$ \vec{p} = \left( \frac{n_x}{d_x} , \frac{n_y}{d_y} \right) $$
and the iteration for the point would proceed
\begin{eqnarray*}
\vec{z} &\rightarrow& \vec{z}^2 + \vec{p} \\
\left( \frac{a_x}{b_x} , \frac{a_y}{b_y} \right) &\rightarrow&
\left( \frac{a_x^2}{b_x^2} - \frac{a_y^2}{b_y^2}
, \frac{2 a_x a_y}{b_x b_y} \right) +
\left( \frac{n_x}{d_x} , \frac{n_y}{d_y} \right)\\
&=&
\left( \frac{a_x^2 b_y^2 - a_y^2 b_x^2}{b_x^2 b_y^2} +
\frac{n_x}{d_x}
, \frac{2 a_x a_y}{b_x b_y} + \frac{n_y}{d_y} \right) \\
&=& \left(
\frac{(a_x^2 b_y^2 - a_y^2 b_x^2)d_x + n_x b_x^2 b_y^2}
     {b_x^2 b_y^2 d_x},
\frac{2 a_x a_y d_y + b_x b_y n_y}{b_x b_y d_y}
\right).
\end{eqnarray*}
Easy this may be, but I can't claim it'll be fast; near the boundaries of the Mandelbrot set, the numerators and denominators in the iteration will get pretty big fairly soon, even if you reduce the rationals by dividing by the greatest common denominator.


Why BMP file format last line of image is first line on file?

Here's a quote from Petzold:

    So, in DIBs, the bottom row of the image is the first row of the file,
    and the top row of the image is the last row in the file. This is
    called a bottom-up organization. Because this organization is
    counterintuitive, you may ask why it's done this way.

    Well, it all goes back to the OS/2 Presentation Manager. Someone at
    IBM decided that all coordinate systems in PM—including those for
    windows, graphics, and bitmaps—should be consistent. This provoked a
    debate: Most people, including programmers who have worked with
    full-screen text programming or windowing environments, think in terms
    of vertical coordinates that increase going down the screen. However,
    hardcore computer graphics programmers approach the video display from
    a perspective that originates in the mathematics of analytic geometry.
    This involves a rectangular (or Cartesian) coordinate system where
    increasing vertical coordinates go up in space.

    In short, the mathematicians won. Everything in PM was saddled with a
    bottom-left origin, including window coordinates. And that's how DIBs
    came to be this way.

  Source: Charles Petzold, Programming for Windows 5th Edition, Chapter
  15.

from
https://stackoverflow.com/questions/8346115/why-are-bmps-stored-upside-down


Can you batch vertex data and still sort by depth?

Drawing back-to-front is not optional with blending (unless you use an order-independent blending technique), so getting the visual effect you actually want has to trump performance.

But drawing front-to-back is purely a performance thing. Specifically, it makes it possible for early-Z to save you fragment processing time, as earlier objects are in front of later ones. So it's a GPU performance improvement.

Batching is a performance improvement for the CPU, in that you're spending less CPU time doing state change stuff.

So you need to decide which is more important: CPU performance or GPU performance. With deferred rendering, saving FS computation time also means saving memory bandwidth.

But the general construction of your scene, what you're actually drawing, can help inform which is more important.

Oftentimes what happens is that you render static terrain first, perhaps using your visibility system to order what you render coarsely front-to-back. Then you render your dynamic objects. Because dynamic objects usually require a lot of state changes, you won't necessarily draw all of them in a single batch (though there are techniques you can use to hide even texture changes).

In any case, the point is that there is a tradeoff. And the one you pick depends on exactly what you're rendering. A top-down game like StarCraft or DotA would never bother with front-to-back sorting; you just render the terrain after the characters. A tight-quarters FPS by contrast would probably really need to employ such techniques, particularly for terrain.


Could some give an explanation or hint about this kind of equation? $\left(- \sqrt{X^{2} + Y^{2}} + 1\right) \cos{\left (2 \pi X + \phi \right )}$

The R part generates a cone where the point is at (0, 0, 1) and it spreads out below that. It meets the x-y plane at the unit circle:

The cos component defines an extruded cosine wave along the X-axis:

So the final equation takes the cosine wave and multiplies its amplitude by the cone. It will have full amplitude at the origin, and decrease to 0 at the unit circle on the x-y plane. Past that it becomes negative:

The phi component shifts the sine wave forward or backward along the x axis.


Need help reasoning about video oscillation >60hz

Yeah, that makes sense. Most flat panel monitors on the market have a 60Hz refresh rate! So you are not going to be able to flicker faster than 60Hz due to technical limitation.

Were you to do this with a electron sweeping oscilloscope, or a LED lamp connected to a PWM source fast enough then no problem (you can find this in arduino samples). It would work.

You can also emulate the effect by rendering multiple frames and then blending them. But you should take careto do correct kind of blending as the monitor output is not nesseserily linear!


How will the z-buffers have the same values even if polygons are sent in different order?

What your teacher means is that at the end of the frame, once all polygons in the scene have been drawn, the z-buffer will have the same values. This is because the z-buffer keeps the minimum Z within each pixel, over all the polygons drawn to that pixel, and the minimum doesn't depend on the order.

If you look at the z-buffer in the middle of the frame, when some polygons, but not all, have been drawn, then naturally the contents of the z-buffer will depend on what's been drawn so far, which would depend on the order in general.

With regard to the color buffer, the same things are true if all the Z values within each pixel are distinct. However there's one case where the results do depend on the order: when multiple polygons have the same Z value at a given pixel. This could happen if the same polygons are rendered twice with different colors, or some such. Then, the behavior depends on the setting of the depth test (i.e. glDepthFunc in OpenGL). If it's strictly "less", then the first polygon drawn that has the minimum Z will "win" and its color will end up in the color buffer at the end. But if the depth test is set to "less-equal", then each polygon that has the minimum Z will overwrite the one before it, and the last one drawn will end up in the color buffer at the end. So in that case, the order of drawing polygons, and the state of the depth test, do make a difference.


How do we know if a ray intersected a triangle in one in a million

If you want to know all triangles that are hit by the ray, you have to check all of them. However, that doesn't mean that you have to check each of them individually. Usually, you apply some sort of space partitioning to your scene and group related objects in a bounding volume hierarchy. Then it gets a lot more efficient because you can drop multiple triangles with a single check.
For example, consider a scene with a rectangular area. You might subdivide this area into a 10x10 grid. Now check which parts of the grid are hit by the ray. All objects located in grid cells that aren't hit by the ray can be discarded and you only need to check the remaining objects. Now let's say that one of the hit objects is a human. Before you test the individual triangles, you narrow down the affected area in multiple substeps by checking against its bounding volume hierarchy. The traversal of the hierarchy might look like this: upper body hit -> left arm hit -> hand hit -> thumb hit. So with 4 checks, you can discard all triangles of the humans' mesh that do not belong to the thumb.
In case you only want to find the first triangle that is hit by the ray, you can also sort all hit cells and objects by distance before performing detailed tests. Then you start your checks at the closest object and continue until you get the first hit.


How to capture frames in a MacOS MTKView for recording to a video

I've seen no answer to my question. But after lots of experimentation and further study of Apple documentation for clues I have figured out how to capture GPU generated frames.

First, I gave up on the attempt to have the shader write to the screen and an offscreen buffer in the same render pass. The odd (and undocumented as far as I can tell) behavior of MTKView in which the onscreen view occupies only a portion of the original area of the screen. I think this a bug in MTKView. It seems to act as if each attachment should appear somewhere on the screen. MTKView is not checking to see if the attachments are all on the screen.

Since MTKView calls the draw() function 60 times per second I decided to alternate. One draw() goes to the screen and the following draw() goes to my off screen buffer. This gives me a 30 Hz frame rate for the captured video frames, which is what I wanted anyway. I can watch the screen and see what will appear in the video.

For the video capture draw() calls I now define the MTLTextureDescriptor in this way:

let offScreenPixelFormat = MTLPixelFormat.rgba32Float   // 32Float is convenient for AVFoundation video recording
let offScreenBufferDescriptor = MTLTextureDescriptor.texture2DDescriptor(pixelFormat: offScreenPixelFormat,
          width: offW, height: offH, mipmapped: false)
offScreenBufferDescriptor.usage = [.shaderWrite, .renderTarget]
offScreenBufferDescriptor.storageMode = .managed
offScreenBufferDescriptor.resourceOptions = [.storageModeManaged]
offScreenBufferDescriptor.textureType = MTLTextureType.type2D
offScreenTextureBuffer = device.makeTexture(descriptor:    offScreenBufferDescriptor)

It turns out this still did not work. I got nothing back from the GPU. That is, I got nothing back until I added this code to the end of the commandBuffer:

  let blitCommandEncoder = commandBuffer.makeBlitCommandEncoder()
  blitCommandEncoder?.synchronize(resource: offScreenTextureBuffer!)
  blitCommandEncoder?.endEncoding()

This causes the GPU memory version of the offScreenTextureBuffer to be copied back to the CPU memory version. I figured this out by looking at the Apple documentation. Like all Apple documentation, it doesn't tell you why you might want to synchronize. But you definitely do need to do this. I have seen several online code examples where this is NOT done and they claim to get data back from the GPU.
I tried one example code without the synchronize command and nothing was returned. So I assume the requirement to synchronize is a new feature and the online code is just out of date. It would be nice if the Apple documentation made it clear this is now a requirement.

With these changes I can retrieve GPU frames. I have created a short video using GPU-generated frames.

If anyone is interested in seeing a complete code sample let me know.


I'm trying to get a 2D screen position so I unproject a point from 3D into 2D screen coordinates by that following function:

Not sure if Ogre does it implicitely, but using a 4x4 matrix to transform a 3D point to get the correct 2D coordinates (as it would be rendered with the same matrix),  requires a transition from homogeneous coordinates by dividing the result by w.

In other words, you must multiply the matrix by a 4-component vector (x,y,z,1.0f) instead of a 3-component vector. And then the resulting vector will look like (X,Y,Z,W) and you need to divide by W to get the position like it would be rendered.

Like this:

// Transform world coordinates.
Ogre::Vector4 transformedWorldPosition = worldViewProjection * Ogre::Vector4(worldPosition,1.0);

//Divide by w
transformedWorldPosition = transformedWorldPosition / transformedWorldPosition.w;

// Calculate screen coordinates from world coordinates.
Ogre::Vector2 screenPosition;
screenPosition.x = (int) ( ( ( transformedWorldPosition.x + 1 ) / 2.0 ) * screenRes.width + 0.5 );
screenPosition.y = (int) ( ( ( 1 - transformedWorldPosition.y ) / 2.0 ) * screenRes.height + 0.5 );


If you're using Python, you can use a library such as matplotlib for this. It can take a numpy array and plot it as a standard chart using the function plot or as an image using imshow. The resulting image/plot can be displayed in an interactive window, or saved to a file, or it shows inline in your Jupyter notebook if you're using that.

There are other (better?) Python chart plotting libraries but the functionality is basic and IMO matplotlib is easier to get started with (perhaps just because it's older and more established).

You can also do this using OpenCV, if your input is an OpenCV image. Its plotting functions are named similarly to the matplotlib ones. OpenCV is available in Python, C++, and a handful of other languages.


There are  a couple of ways to answer this question: an algebraic way and a geometric way.

Algebraically, we can identify the units that the BRDF must have by looking at its place in the rendering equation. The classic rendering equation is:

$$L_\text{outgoing}(\omega) = L_\text{emitted}(\omega) + \int_\Omega L_\text{incoming}(\omega') \, f_\text{BRDF}(\omega, \omega') \, (n \cdot \omega') \, d\omega'$$

The output value on the left is a radiance, so the result of the integral must also be a radiance. The integrand contains a radiance multiplied by a solid angle $d\omega'$, so something else in the integrand has to cancel out that factor of solid angle. The $n \cdot \omega'$ factor is dimensionless, and the only other thing there is the BRDF—so to make the whole thing work out, the BRDF must have units of inverse solid angle. Equivalently, the BRDF can be seen as a ratio of radiance to irradiance, since they differ by a factor of solid angle in the denominator of radiance.

Another way to see it is that the BRDF plays a role similar to a probability density. If you look at how probability densities work, they have units inverse to the volume of their domain. For instance, a 1D probability density has units of inverse length (probability per unit length, but probability itself is dimensionless), a 2D one has units of inverse area, and so on. The BRDF acts much like a probability density defined on the hemisphere, giving a likelihood for a photon coming in from a given direction to be reflected into some other direction. So, like any other probability density on a spherical domain, it has units of inverse solid angle.

Geometrically, we can get right down to brass tacks and take apart what's going on in the integral in the rendering equation. Recall that an integral means to subdivide the domain into tiny pieces and sum the integrand over all the pieces (in the limit as the pieces get infinitesimally small). Let's look at one such piece. The integrand should result in an infinitesimal amount of radiance $dL$, since we're going to sum over many pieces to arrive at a finite outgoing radiance. So a single infinitesimal piece of the integral looks like:

$$dL = L_\text{incoming} \, f_\text{BRDF} \, (n \cdot \omega') \, d\omega'$$

If we regroup the factors a bit, the combination $L_\text{incoming} \, (n \cdot \omega') \, d\omega'$ calculates the irradiance on the surface due to the light coming from the infinitesimal solid angle $d\omega'$. Since it's arriving from an infinitesimal amount of solid angle, it produces an infinitesimal irradiance $dE$.

$$dL = f_\text{BRDF} \, dE$$

or

$$f_\text{BRDF} = \frac{dL}{dE}$$

So the BRDF acts as a proportionality constant between the infinitesimal irradiance arriving at the surface from an infinitesimal solid angle, and the infinitesimal outgoing radiance generated thereby. It couldn't be a ratio of radiances, because we have a finite incoming radiance, and we need an infinitesimal outgoing radiance if we want to sum up many pieces of the integral and get a finite result. To make that happen, the BRDF would have to be infinitesimal-valued, which...isn't a thing, in standard mathematics. :)

I hope some of this helps. There are a variety of equivalent ways to look at this problem, as with so many things in math and physics.


Compute the interpolation factor between the SLERP of 2 quaternions, given those quaternions and the SLERP result

We can do this mostly by measuring angles between the quaternions, with a bit of sign flip to get the orientations right. Angles between normalized quaternions are given by the arccos of their dot product (as a 4D vector).
$$
t = \frac{\arccos(q_1 \cdot q)}{\arccos(q_1 \cdot q_2)} \ \text{sign}\bigl(q \cdot [q_2 - (q_1 \cdot q_2) q_1]\bigr)
$$
The first factor is measuring the angle from $q_1$ to $q$ as a fraction of the angle from $q_1$ to $q_2$. This ensures that $t$ will be 0 at $q_1$ and 1 at $q_2$. The other factor ensures that $t$ increases as you move from $q_1$ toward $q_2$ and goes negative on the other side. The quantity in brackets is the projection of $q_2$ on the (hyper)plane perpendicular to $q_1$, then we're dotting $q$ with that to decide which side of the circle we're on.
Note that the point where $t$ swaps from negative to positive will be the point opposite $q_1$ on the circle—not a point centered across from $q_1$ and $q_2$ as your diagram has it.
(Also, it's worth noting that this boils down to essentially the same thing as the logarithm solution in the other answer. Their factor $\check p/\|\check p\|$ corresponds to our sign factor: all the quaternions in the problem are in the same plane and so will have $\check p$ along the same axis; the only thing that can differ between the normalized $\check p$s is the sign.)


What are the outcodes that would be assigned to these points?

The diagram on the Wikipedia entry for Cohen-Sutherland Clipping explains it pretty well:

        left    central right
top     1001    1000    1010
central 0001    0000    0010
bottom  0101    0100    0110

So in your example, p1 would be left, central, so its code  would be 0001.

p2 would be central, bottom, so 0100

p3 would  be  central, central, so 0000

p4 would be central, top, so 1000


What is Tessellation in computer graphics

Tesselation is a technique that allows you to generate primitives (triangles, lines, points and such) on the graphics-card. Specifically, it lets you repeatedly subdivide the current geometry into a finer mesh.

This allows you to load a relatively coarse mesh on your graphics card, generate more vertices and triangles dynamically and then have a mesh on the screen that looks much smoother.

Most of the time this tesselation is done anew in each single frame and this could be the reason that your frame-rate drops once you enable this.

Tesselation is done in multiple stages, and it is done AFTER the vertex shader.

The terms for each stage varies based on the API. In DirectX, it is the Hull Shader, Hardware Tessellation, and the Domain Shader. In OpenGL, they are called the Tessellation Control Shader, Tesselation Primitive Generation, and the Tessellation Evaluation Shader.

The first and the last stage is programmable, the actual tesselation is done by the hardware in a fixed function stage.

In the Tesselation Control shader you set the type and number of subdivisions.

Then the hardware tessellator divides the geometry according to the Control Shader.

Lastly, the Tessellation Evaluation Shader is called for each newly generated vertex. In this shader, you set the type of primitive you want to generate and how the vertices are spaced and many other things. This shader can also be used to do all sorts of per-vertex calculations just like a vertex shader. It is guaranteed to be called at least once for each generated vertex.

If you want to do any further work on the primitives, you can add a Geometry shader.


What is the most efficient line algorithm using a shader program

Here is a (mostly) minimal version of an SDF for a line segment using a signed distance field as suggested in the comments. I take no credit for this code, it is your function pieced together with IQ's SDF function done in the style you listed above.
Be sure to go to IQ's website, it is one of the best learning resource out there for this material. Though it can be hard to digest till you get all the basics down. IQ's Index Page.
float udSegment( in vec2 p, in vec2 a, in vec2 b )
{
   vec2 ba = b-a;
   vec2 pa = p-a;
   float h =clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 );
   //return length(pa-h*ba); // use actual length

   // or use length squared
   vec2 sqrd = pa-h*ba;
   sqrd = sqrd * sqrd;
   return sqrd.x+sqrd.y;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
   vec2 p = ((2.0*fragCoord-iResolution.xy)/iResolution.y)*1.4;

    // give it some spin
    vec2 v1 = iResolution.xy * vec2(cos(iTime),sin(iTime));

    vec2 v2 = vec2(0,0); // set the start point
    float th = 0.004; // Controls line thickness

    float d = udSegment( p, v1, v2 ) - th; // compute sdf

    vec3 col = vec3(1.0) - sign(d);
    fragColor = vec4(col,1.0);
}

If all this is doing is drawing a line, then there is no need to compute the actual length, so don't...just return lenght squared. I made an edit to the code above that skips the sqrt, but line thickness is now squared.
Bresenham's algorithm was intended to be used for drawing individual pixels and is often implemented directly in the hardware of GPU's. So it drives which pixels to draw in its own loop and tells the hardware which pixel to draw.
ShaderToy is drawing the entire screen, so every pixel is covered probably by drawing a fullscreen triangle, then calling the mainImage function in the fragement shader.
To implement Bresenham's algorithm in ShaderToy the code would need to determine if a particular fragment was actually on the line and in so doing we lose any advantage gained by this clever line drawing technique.
The most efficient way to draw a line on a GPU is not to use ShaderToy at all but write our own bit of code that tells the hardware to just draw a line, this is where graphics API's like OpenGL come into play. Which ultimately would call a fragment shader but unlike ShaderToy, the fragment shader would only be called for those fragments which are actually on the line and our fragment shader would then be reduced to something ultra simple like:
void main(){ return vec4(1,1,1,1); }

Which would result in a pure white line being drawn on screen.


Confusion about how inverse bind pose is actually calculated and used?

To get a vertex in root space into C space, I would have to do (CBA) * v.

Well, yes. But that's not actually what you want to do in skeletal animation. You have it reverse.

It's the other way around. You assume you have a point in C space and want to compute its position in root/global space, in order to pump it further through the transformation pipeline, which gets increasingly more global with each transformation until your point ends up in eye-space (and further). That's why A doesn't transform from root space into A space, it transforms from A space into root space. So $ABC$ transforms from C space into root space.

Now your mesh vertices already are in root space, that's where the inverse bind pose comes into play, which transforms a point from root space into C space. So you basically transform your point from root space to C (or whatever joint) space with the inverse bind matrix that's only dependent on the joint and constant over the whole animation and then from that joint space back into root space with the animation-dependent joint transformation.

So C's transformation matrix for frame $t$, that transforms a vertex from unanimated root space, which is how your mesh is defined, into animated (with respect to C's joint angle) root space would be $(A_RB_RC_R)^{-1}A(t)B(t)C(t)$, where $M(t)$ is the local joint matrix at animation time $t$ and $M_R$ is the local joint matrix corresponding to the rest pose, the pose in which the mesh is defined, and thus $(A_RB_RC_R)^{-1}$ is C's inverse bind pose.


Reading neighbouring pixels during object draw call

In general you can't safely read from neighboring pixels in the render target you're currently writing to. You can sample from them, but there's no guarantees about how other pixels are scheduled on the GPU, so it's a race condition—you might get the neighboring pixel value either before or after modification, and it might change frame to frame due to unpredictable scheduling differences.
You can read from the same pixel you're about to write to, but only if your geometry is set up to ensure that nothing is overlapping in screen space (e.g. a full-screen pass). If other triangles in the same render pass are overlapping the same area on screen, then you have a race condition again, similar to the neighboring-pixels case. (You can use fragment_shader_interlock or equivalent functionality to address this case, but it is not universally supported across all GPU vendors, and it can also be slow.)
The depth and stencil buffers, as well as fixed-function blending (additive, alpha, etc) can get around this because they have special-purpose hardware that sorts the fragments in submission order in order to apply those simple fixed operations.
I believe the way stochastic transparency usually works is to randomize the set of pixels output and use spatial/temporal filtering to clean up the noise afterward. So you wouldn't have such fine-grained control like "within each 2x2 block, there will be N pixels of this, and M pixels of that". It would be more probablistic, like "on average, X% of the pixels will be this, and Y% of the pixels will be that", with the exact outputs being randomly chosen each frame.


Why negate z when constructing projection matrix OpenGL

I'm assuming your scene is constructed based on right-handed coordinates.
If you are using OpenGL, yes. If you are using Direct3D, no.
The projection matrix maps [-n, -f] into [0, 1]. This weird property comes from the fact that the eye coordinates are right-handed, but the clip space (NDC) uses left-handed coordinates. Hence the implementation of ortho().
Direct3D uses a z-axis that goes outwards the screen. Implementations for Direct3D projection matrices will not produce negated z when coupled with your application.


Why are depth buffer values clamped to (0, 1)?

(Just converting comment into an answer)
It seems you might be asking one of two possible questions.

If you meant "Why doesn't the Depth buffer contain camera space Z"
then it's perhaps worth reading Why do GPUs divide clip space Z by W, for position?

If, OTOH, you just meant "why is it restricted to [0,1)?", then that's probably because old systems represented the Z-buffer depth with fixed point numbers, hence a natural [0,1) range.
 Of course, later hardware arrived that could use floating-point depths; IIRC Dreamcast's depth representation was an arbitrary positive float representing $\frac{1}{w}$.


Modify camera matrix so that it renders a "chunk" of the frustum

It depends a little how you construct the matrix. I assume that you use the function perspective_matrix from your library, which actually works pretty much like good old gluPerspective). But what this actually does is nothing else than use the more general function frustum_matrix, which in turn creates a common arbitrary view frustum matrix, like good old glFrustum did (from the (simplified) linalg source code):

mat4 perspective_matrix(T fovy, T aspect, T n, T f)
{
    T y = n*std::tan(fovy / 2), x = y*aspect;
    return frustum_matrix(-x, x, -y, y, n, f);
}

And this general frustum matrix function is actually all you need for creating arbitrary assymmetrical frustum matrices, like the ones you want. The paramters it takes are exactly the left, right, bottom, and top boundaries of the near plane, i.e. the x and y boundaries of the frustum on the z = near (or -near) plane. So all you actually need to do is take these left/right/bottom/top parameters and cut those ranges into whatever sub-intervals you need.

It depends a little on what information you actually have. If you have the parameters you construct your perspective matrix from (i.e. the fovY and the aspect) you can just replace the above procedure with something more elaborate that computes the frustum parameters from these two and then scale-biases them into the corresponding intervals.

If all you have is an already constructed perspective matrix, you might have to retrieve the parameters from that matrix first, which shouldn't be too difficult when looking at how the matrix is constructed. If you already have a general assymmetric frustum matrix, it might get a little more elaborate. But it always really all just comes down to tweaking the $M_{11}, M_{13}, M_{22}, M_{23}$ elements of the matrix a little with some scaling and offseting that I'm sure you can figure out when looking at how these matrices are constructed. The rest can stay as they are (provided you don't need chunking in the z-direction either, which would be used for different techniques than what you are trying to do).


There is algorithm called as delaunay triangulation which does triangulation or points it's comparably simple to understand . But it is very slow. If you don't want to implement your own algorithm from scratch please have look at pcl
http://pointclouds.org/documentation/tutorials/greedy_projection.php
Also have look at CGAL
http://www.cgal.org.


Why do T-junctions in meshes result in cracks?

lhf's answer is good from the perspective of tessellation, but these can occur with simpler triangle mesh use cases.

Take this trivial example of three, screen-space triangles, ABC, ADE and DBE...

Although point E was, mathematically, intended to be exactly on the line segment AB, the pipeline won't be using fully precise values, such as rational numbers (e.g. https://gmplib.org/). Instead, it will likely be using floats, and so some approximation/error will be introduced. The result is probably going to be something like:

Note that all of the vertices may have inaccuracies. Although the example above shows a crack, the T-junction may instead result in overlap along the edge causing pixels to be drawn twice. This might not seem as bad, but it can cause problems with transparency or stencil operations.

You might then think that with floating-point the error introduced will be insignificant, but in a renderer, the screen-space vertex (X,Y) values are nearly always represented by fixed-point numbers and so the displacement from the ideal location will usually be much greater. Further, as the rendering hardware "interpolates" the line segment pixel-by-pixel with its own internal precision, there is even more chance it will diverge from the rounded location of E.

If the T-junction is "removed" by, say, also dividing triangle ABC into two, i.e. AEC and EBC, the problem will go away as the shifts introduced by the errors will all be consistent.

Now, you might ask why do renderers (especially HW) use fixed-point maths for the vertex XY coordinates? Why don't they use floating-point in order to reduce the problem? Although some did (e.g. Sega's Dreamcast) it can lead to another problem where the triangle set-up maths becomes catastrophically inaccurate, particularly for long-thin triangles, and they change size in unpleasant ways.


How do I interpret the format of the .flame files used by apophysis?

There is some information about the XML file format used for .flame files. There are a number of forks of this program written in different languages, the original site is by definition the place to go for answers; unless you are using a fork - the original  seems to have been abandoned years ago.

Coefs is 6 floats which defines the coefficients for the affine part of the function, in row order (a d b e c f). Xform has an interpolation_type which is a string describing the method for interpolation of xform positions: options are "log", "linear", "old", etc. (a list of a few dozen options is here). There can be multiple xforms, each applies a transformation to the previous state.

The spherical transform is complex, there's a 19 page .docx tutorial with 9 pages of examples at: "Understanding Spherical Transforms: Apophysis Tut" by UniversalKinase (direct link to .ZIP file). Another tutorial is "Apo-Tuto Sphericals" by Fiery-Fire.


You need to learn and understand STEP data exchange format (ISO 10303-21) which is the fundamental exchange format which the IFC schema is built on top of. This is no easy feat though, especially when comparing to a file format such as OBJ. OBJ files are designed to define geometry (and other information via associated MTL) and so are optimised for this purpose. IFC on the other hand, is designed to exchange many aspects of AEC including building, services, structural, furnishings, energy model etc and so is far more complex with the geometry definitions being only a subset of the information supplied in an IFC. It would be easier to use an existing IFC library to do this, and even then, an understanding of the structure of the IFC entities and their hierarchy is still required to traverse the relationships (inverse attributes) between the object definitions in an IFC file in order to reconstruct the desired geometry.

OBJ files explicitly define the vertices, normals and faces that represent a mesh, however while there are IFC entities which can represent meshes in a similar fashion (namely IfcConnectedFaceSet), most geometry defined in an IFC file should be considered more like a series of instructions/commands to recreate the geometry. For this reason you would also require a geometry library (one which can do things like extrusions, Boolean/CSG operations, curves, NURBS etc). For example, one of the main things that is represented in many IFC files would be the walls, doors, windows etc of a building, and while these could be represented as explicit tessellated data or breps, in my experience these tend to be represented as an IfcExtrudedAreaSolid (https://standards.buildingsmart.org/IFC/RELEASE/IFC2x3/TC1/HTML/ifcgeometricmodelresource/lexical/ifcextrudedareasolid.htm) which supplies the direction and length to extrude a path/outline curve (the curve being defined by an IfcProfileDef which is an attribute of IfcExtrudedAreaSolid, https://standards.buildingsmart.org/IFC/RELEASE/IFC2x3/TC1/HTML/ifcprofileresource/lexical/ifcprofiledef.htm). This gets even more complicated when there is an opening (i.e window recess or door opening in the wall) which again is related (through an inverse attribute) to the geometry representing the opening to be removed through a geometric Boolean operation (CSG). As you can see it starts to get complicated fairly quickly and isn't just a case of reading in raw vertices, normals and faces etc.

All geometry in an IFC file is defined by objects which are sub types of IfcRepresentationItem (http://dcom.arch.gatech.edu/aci/cic/schema/ifcgeometryresource/lexical/ifcrepresentationitem.htm) and so that would be the place to start to understand what types of geometry definitions are used by IFC. Other items within the model (namely furnishings) are far more likely to be represented by raw tessellated data (IfcFacetedBrep, https://standards.buildingsmart.org/IFC/RELEASE/IFC4_1/FINAL/HTML/schema/ifcgeometricmodelresource/lexical/ifcfacetedbrep.htm) which is geometry that would be more familiar to you, and there is a chance you can isolate these definitions within in an IFC file and concentrate on processing that (you would still need to understand STEP format, namely the inverse attributes involved, in order to faithfully generate the mesh geometry). Also this geometry could be either explicitly located (in world space), but in many cases is also positioned through IfcPlacements (https://standards.buildingsmart.org/IFC/RELEASE/IFC2x3/TC1/HTML/ifcgeometryresource/lexical/ifcaxis2placement3d.htm) locally, so you need to build up the hierarchy of transformations based on these definitions.

In short, if you wanted to RYO IFC parser to grab geometry, it would be a steep learning curve and require you to understand many other concepts of STEP/IFC outside just knowing the geometry types and transforms. Easier to use an IFC library coupled with a geometry library (although many IFC libs do come with geometry libs to do this for you), but even then, depending on the lib used, you may need to have a basic understanding of the building hierarchy involved in order to correctly place and reconstruct much of the geometry.

Also note, while many of the geometry concepts remain the same through the various versions of IFC (IFc2x3, IFC4, IFC4 ADD1 etc), in some cases there are subtle differences in the schema in relation to geometry object definitions. However, much of the version difference in IFC is do to with the building hierarchy, and associated information for the various disciplines in AEC.