When unsetting a VAO, should you also unbind the associated VBOs?

No, you do not need to unbind the VBO's. VBO's are only added to the state of the VAO when that VAO is bound. So any previously bound VBO's will not alter the state of the VAO that you are going to bind. The state of a bound VAO can be changed with the following calls:

glEnableVertexAttribArray
glDisableVertexAttribArray
glVertexAttribPointer
glBindBuffer(GL_ELEMENT_ARRAY_BUFFER)

You can bind and unbind VBO's without a currently bound VAO and even share them amongst VAO's. The VBO is simply a pointer to a data container that can be used for multiple purposes.


One strategy mentioned in Brian Karis's talk about TAA is neighborhood clamping. The general idea is that, for the previous frame's pixel to be valid, its color should be in the color range found in the neighborhood (say 3x3 pixels) of the current pixel this frame.

This rejects history from changing light conditions, which is probably what you want anyway if you don't want moving shadows to produce ghosting.

(Animated textures, depending on the speed of the animation, could also be handled with a motion vector, if you have a predictable UV mapping or can guess reasonably well.)


How to draw an orthogonal graph from its edges and vertices data?

libCola

cola.js (a.k.a. WebCola) is a JavaScript based rewrite of libcola which works well with D3.js

http://marvl.infotech.monash.edu/webcola/

Gridified Layout example:

http://marvl.infotech.monash.edu/webcola/examples/dotpowergraph.html

OGDF

Open Graph Drawing Framework

Have an orthogonal layout

  Orthogonal: Ortholayout with modules for the steps of shape and
  compaction, ClusterOrthoLayout

http://www.ogdf.net/doku.php/ogdf:features

Check projects that uses OGDF

http://www.ogdf.net/doku.php/project:external

HOLA

Human-like Orthogonal Layout Algorithm

https://vimeo.com/150509722

Kieffer, Steve, Tim Dwyer, Kim Marriott, and Michael Wybrow. "Hola: Human-like orthogonal network layout." Visualization and Computer Graphics, IEEE Transactions on 22, no. 1 (2016): 349-358.


How does graphics api like OpenGL determine which triangle is back face to cull?

It's actually fairly simple (and in essence 'equivalent' to what you said above). If you represent triangle ABC's projected XY coordinates as a matrix...

$$M_{ABC} =   \begin{bmatrix}
    X_A & X_B & X_C \\
    Y_A & Y_B & Y_C \\
    1 & 1 & 1 \\
    \end{bmatrix}
$$

...then the sign of the determinant e.g. $X_A (Y_B - Y_C)+X_B(Y_C-Ya)+X_C(Y_A-Y_B)$ will allow you to decide how to cull.
Further if the determinant is zero, then the on-screen triangle has zero area and can clearly be removed.

Another reason to use this notation is that the adjoint of $M_{ABC}$  gives you equations for the edges.


What is wrong with it seems to be explained on the next two slides: it leads to non-energy-conserving results, where the apparent brightness of the volume changes depending on the scattering coefficient.

My read of it is that the "wrong" integration code implicitly assumes a constant transmittance over the extent of each voxel. But the trouble is, when the scattering coefficient is high, the transmittance will fall off significantly within a single voxel. For instance, if you have voxels of size 10 cm, but the scattering coefficient is 1 / (5 cm), then the transmittance falls off by exp(10 cm / 5 cm), or about a factor of 7.4 over the length of the voxel!

Therefore, according to slide 28, a better approach is to explicitly account for the transmittance falloff over the voxel length, while still assuming constant inscattered light. Fortunately the integral has a simple closed-form solution. So, they're saying to evaluate that formula per voxel in the integration loop, in place of the former "wrong" one.

For completeness: in the code from slide 27, replace this line:

accumScatteringTransmittance.rgb += scatteringExtinction.rgb *
    accumScatteringTransmittance.a;

with this:

accumScatteringTransmittance.rgb += scatteringExtinction.rgb *
    accumScatteringTransmittance.a *
    ((1.0 - transmittance) / scatteringExtinction.a);

Here, accumScatteringTransmittance.a is the transmittance at the point the ray enters the voxel (i.e. the resultant transmittance of all the voxels in front of it), and then the ((1.0 - transmittance) / scatteringExtinction.a) factor accounts for the varying transmittance within the voxel.


Same specular lighting on both tris in a quad (for "low-poly" "flat shading" look)?

If I understand you correctly, the problem is that the provoking vertex isn't the same vertex between the two triangles in a quad, due to using triangle strips. You could switch to using an indexed triangle list, and set up indices in such a way that the provoking vertex is the same in each quad. For instance, if one quad is vertices {0, 1, 2, 3} (CCW) then you might write {0, 1, 2, 3, 0, 2} in your index buffer so that 2 is the provoking vertex for both triangles. This should still be pretty efficient.

Another possibility (more involved, but more flexible) is to use the built-in gl_VertexID variable in your vertex shader to figure out which quad you're in (eg int quadID = gl_VertexID / 4), and use an appropriate shading location for that quad—the centroid of its vertices, for instance. However, for this to work, I think you would have to "un-share" the vertices so that each quad has its own 4 unique vertices (the two triangles in a quad can still share). So you would have to use primitive restarting, or again switch to indexed non-strip draws. You could also have per-quad data in an SSBO that you manually fetch using your computed quad index.


Why the Bresenham's algorithm does not work for those two points?

Like Rahul said that algorithm is only in the case where $0 \lt \Delta y \lt \Delta x$.

So you need to adjust the coordinates to fit within that.

There are 8 cases:

$0 \lt \Delta y \lt \Delta x$ the normal base case
$0 \lt \Delta x \lt \Delta y$ swap x and y except for the drawPixel call.
$\Delta x \lt 0 \lt \Delta y$ and $|\Delta y < \Delta x|$ use case 1 but negate $\Delta x$ and instead of counting from x0 to x1 you count from x1 to x0
$\Delta x \lt 0 \lt \Delta y$ and $|\Delta x \lt \Delta y|$ use case 2 but negate $\Delta x$ and instead of counting from x0 to x1 you count from x1 to x0

The remaining cases have $\Delta y \lt 0$ so swap the 2 points and it will fall within the previous 4 cases.


How to render an object that recieves shadows but does not cast them in a variance shadow mapping system?

Not rendering a shadow caster could work but your mileage can vary. When variance is "large" a VSM is not able to exactly localize the occluders, which causes light leaks. Not rendering a large receiver (let say the floor of a room, a road, etc..) can dramatically increase the variance of the depth samples within a region of the shadow map that include the missing receiver(s) (e.g. the receiver depth is replaced by something much more distant, like the far plane).

On the other hand, if you skip rendering in the VSM some little object that is unlikely to ever be the most distant receiver in the scene, then you have some hope that things will work out (compatibly with the limitations of variance shadow maps). If you can afford the cost I suggest you to use something more robust like exponential variance shadow maps (see this work: https://software.intel.com/en-us/articles/sample-distribution-shadow-maps).


An Affine Transform is a Linear Transform + a Translation Vector.

$$
\begin{bmatrix}x'& y'\end{bmatrix}
=
\begin{bmatrix}x& y\end{bmatrix}
\cdot
\begin{bmatrix}a& b \\ c&d\end{bmatrix}
+
\begin{bmatrix}e& f\end{bmatrix}
$$

It can be applied to individual points or to lines or even Bezier curves. For lines, it preserves the property that parallel lines remain parallel. For Bezier curves, it preserves the convex-hull property of the control points.

Multiplied-out, it produces 2 equations for yielding a "transformed" coordinate pair $(x', y')$ from the original pair $(x, y)$ and a list of constants $(a, b, c, d, e, f)$.
$$
    x' =  a\cdot x + c\cdot y + e \\
    y' =  b\cdot x + d\cdot y + f
$$

Conveniently, the Linear transform and the Translation vector can be put together into a 3D matrix which can operate over 2D homogeneous coordinates.

$$
\begin{bmatrix}x'& y'&1\end{bmatrix}
=
\begin{bmatrix}x& y&1\end{bmatrix}
\cdot
\begin{bmatrix}a& b &0\\ c&d&0 \\ e&f&1\end{bmatrix}
$$

Which yields the same 2 equations above.

Very conveniently, the matrices themselves can be multiplied together to produce a third matrix (of constants) which performs the same transformation as the original 2 would perform in sequence. Put simply, the matrix multiplications are associative.

$$
\begin{matrix}
\begin{bmatrix}x''& y''&1\end{bmatrix}
& = &
\left(
\begin{bmatrix}x& y&1\end{bmatrix}
\cdot
\begin{bmatrix}a& b &0\\ c&d&0 \\ e&f&1\end{bmatrix}\right)
\cdot
\begin{bmatrix}g& h &0\\ i&j&0 \\ k&m&1\end{bmatrix}
\\
& = &
\begin{bmatrix}a \cdot x + c \cdot y+e & b \cdot x + d \cdot y+f &1 \end{bmatrix}
\cdot
\begin{bmatrix}g& h &0\\ i&j&0 \\ k&m&1\end{bmatrix}
\\
&=&
\begin{bmatrix}g(a \cdot x + c \cdot y+e) + i( b \cdot x + d \cdot y+f) + k
\\
h(a \cdot x + c \cdot y+e) + j( b \cdot x + d \cdot y+f) + m
\\1 \end{bmatrix}^T
\\
&=&
\begin{bmatrix}x& y&1\end{bmatrix} \cdot
\left(
\begin{bmatrix}a& b &0\\ c&d&0 \\ e&f&1\end{bmatrix}
\cdot
\begin{bmatrix}g& h &0\\ i&j&0 \\ k&m&1\end{bmatrix}
\right)
\\
&=&
\begin{bmatrix}x& y&1\end{bmatrix} \cdot
\begin{bmatrix}ag+bi& ah+bj &0\\ cg+di&ch+dj&0 \\ eg+fi+k&eh+fj+m&1\end{bmatrix}
\end{matrix}
$$

Alternatively you can consider a few basic transform types and compose any more complex transform by combining these (multiplying them together).

Identity transform

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

Scaling

$$\begin{bmatrix}
    S_x & 0 & 0 \\
    0  & S_y & 0 \\
    0 & 0 & 1 \\
\end{bmatrix}$$

*Note: a reflection can be performed with scaling parameters $(S_x, S_y) = (-1,1)$ or $(1,-1)$.

Translation

$$\begin{bmatrix}
    1 & 0 & 0 \\
    0 & 1 & 0 \\
    T_x & T_y & 1
\end{bmatrix}$$

Skew x by y

$$\begin{bmatrix}
    1  & Q_x & 0 \\
    0 & 1 & 0 \\
    0 & 0 & 1
\end{bmatrix}$$

Skew y by x

$$\begin{bmatrix}
    1  & 0 & 0 \\
    Q_y & 1 & 0 \\
    0 & 0 & 1
\end{bmatrix}$$

Rotation

$$\begin{bmatrix}
    \cos\theta & -sin\theta & 0\\
    \sin\theta & \cos\theta & 0\\
    0  &  0 &  1
\end{bmatrix}
$$

[Note I've shown the form of Matrix here which accepts a row vector on the left. The transpose of these matrices will work with a column vector on the right.]

A matrix composed purely from scaling, rotation, and translation can be decomposed back into these three components.


benefit of perlin noise over value noise

The benefit of perlin noise is the overall distribution of frequencies. Since value noise uses simple values that are interpolated, there is a higher chance, that a row of several values only differs a little. The consequence is, that some regions of your picture may contain little changes and some regions a lot of changes.

By using gradients you are reducing this effect because the interpolation is not done by value but instead calculated between tangents. Now it's more difficult to have a flat curve (both tangents must be collinear).

Source: As noted by Martin Ender the question was already posted on a different StackExchange community: see this Math.SE post.


Correct Specular Term of the Cook-Torrance / Torrance-Sparrow Model

I would trust Pharr and Humphreys on this. Equation 2 also agrees with the SIGGRAPH Physically Based Rendering course notes, as well as with equation 20 in the Walter et al paper that introduced the GGX distribution.

I've read somewhere that there is an error in the original Cook-Torrance paper that led them to miss the factor of 4 in the denominator, which was corrected in subsequent papers. I couldn't find a reference to this with a quick search though (if anyone knows one, please feel free to note it in the comments).

As for the factor of π, it might appear or not, depending on conventions. Sometimes it's factored into the normal distribution function D. For instance, if you look in the Walter et all GGX paper section 5.2 where they give the equations for several D functions, you can see they all have a π in the denominator. Note that this implies that the Lambertian BRDF should have a π in the denominator as well.

In real-time graphics, the π is often left out, in which case we can interpret it as having been factored into the light colors. Either way is fine, as long as you're consistent about either putting the π in or leaving it out of all the BRDFs you use.


No, it's not. The way many of the 3D sensors work is by projecting an infrared pattern onto the surface and measuring how it distorts. But with a 2D image, the pattern will simply be projected onto the flat 2D image, not onto the objects in the scene. So the 3D sensor will only sense a flat card.

Other methods work by combining 2 images taken with different cameras at different angles. This will have the same result as above. The 2D photo will be interpreted as a 2D surface by the combination of the 2 cameras.


How to decide which child of BVH to visit based on sign of ray's direction vector?

In your example of a BVH cell partitioned along the x axis, if the ray is going left to right (increasing x), then the cell with the lower x comes first; but if it's going right to left (decreasing x), then the cell with the higher x comes first. So, you can look at the sign of the x component of the direction vector to determine which cell to test first.


Proper architecture for rendering history of streamed data (points)

My advice would actually be to stick with solution 1 unless and until something more is needed. It's the simplest and easiest to understand, and it should be quite feasible to make it perform well.

On the one hand, it may seem wasteful to re-filter the points and re-upload a bunch of data to the VBO every frame. I guess that this is why you say it is a "bad practice". But consider the total amount of data involved. An order-of-magnitude estimate shows that processing this amount of data every frame is not much work for a modern computer.

The Point struct you quote looks like it's about 20 bytes (assuming that p is 3 floats and t is one float or int). With 50,000 points, the total size of the data is about 1 MB. Even if your actual data struct is a bit larger, still the total size of the data is only a few MB—a small amount by today's standards. The time needed to process this amount of data on the CPU, and transfer it to the GPU, is on the order of ~0.1 ms. (Assuming ~20 GB/sec CPU bandwidth, and ~8 GB/sec PCIe transfer bandwidth.) So, this is an entirely reasonable amount of data to re-process each frame, and in fact you could handle many times more without stressing your system at all.

So, just keep things simple and efficient for memory access:

Store the points in a big flat array (maybe using a ring buffer to allow new incoming points to overwrite old ones).
Each frame, scan the points using a tight loop on the CPU, and copy the ones that pass the filter into the VBO.

For best efficiency, probably you should use GL_DYNAMIC_DRAW usage when you create the VBO, then use glMapBufferRange each frame and copy the filtered points directly into the mapped buffer.
Use GL_MAP_WRITE_BIT | GL_MAP_INVALIDATE_BUFFER_BIT for the access parameter to glMapBufferRange; this will allow the GL driver to perform under-the-hood double-buffering to stream the data to the GPU efficiently. For more information on what's going on here, see Buffer Object Streaming.

Then unmap the buffer and issue the drawing commands, using a vertex shader to transform the points to screen-space as usual.


I'm not entirely sure if I understood the question, but it seems like you want to implement environment map illumination in your ray tracer (you seem to be thinking in real-time, though). I will approach the general solution for offline renderers without assuming anything about the scene except the environment sphere.

This is generally done by considering your environment as a light source. Let's say that you use a sphere with a radius big enough to contain your scene. Then, you sample it as you would sample an area light source. You can think of this sphere as a textured light source.

The contribution of the direct lighting on the radiance leaving a point $x$ on the direction $\Theta$ would be:
$$
L_{\text{direct}}(x \rightarrow \Theta) = L_{\text{lights}}(x \rightarrow \Theta) + L_{\text{env}}(x \rightarrow \Theta)
$$

Where the environment map contribution is:

$$
L_{\text{env}}(x \rightarrow \Theta) = {\int_{\Omega_x} L_{\text{tex}}}(x \leftarrow \Psi)f_r(x, \Theta \leftrightarrow \Psi)V(x, \Psi) \cos(\Psi, N_x)\, d\omega_{\Psi} \qquad (1)
$$

Note that, here, $\Psi$ are directions to the sphere, $V(x, \Psi) = 1$ if and only if the sphere is visible from that direction and $V(x, \Psi) = 0$ otherwise, $f_r$ is the BRDF and $L_{\text{tex}}$ is the color on the texture.

Directions $\Psi$ can be generated with a parameterization of the sphere. If you want something simple, then this is enough, just treat it as a textured area light and emit whatever color is on the texture on the sampled point. If you want something more elaborate, consider the rest of the text.

With this said, environment mapping can be inefficient if not implemented with care.

Also, I am assuming you only implemented area lights that emit uniformly, so it was okay to sample its area uniformly. However, it is clear that the sphere is a light source where some parts contribute more than others since the colors vary on the texture. This, combined with the fact that the sphere occupies the entire solid angle around a point (instead of the hemisphere that you usually sample), skyrockets the variance of the estimator of $(1)$. The visibility term $V(x, \Psi)$ further complicates things regarding variance, depending on your scene (if it very exposed or not).

In order to sample and estimate $(1)$, you have to consider:

Parameterization of the sphere: You need to decide how to parameterize your sphere. You need to have a space where you can generate random points and map to a point on the sphere to get a direction $\Psi$ and the color from the texture. You can just use typical the latitude-longitude parameterization.

For more elaborate, better distributed parameterizations: paraboloid parameterization [page 4]; concentric-map parameterization.

PDF $p(\Psi)$ used for sampling the sphere: It is possible to reduce the variance that I mentioned earlier. You can importance sample only using the term $\cos(\Psi, N_x)$ (i.e. cosine lobe sampling), but even so it is possible to do better than that. If the BRDF $f_r$ is simple, you can sample the product of the BRDF with $\cos(\Psi, N_x)$.

There are other methods. In [Kollig, Keller; 03] a PDF is built from the 2D array of the texture itself. This way, it is possible to get more samples from directions from where the map contributes more. A PDF from the product of this with the BRDF would yield the best results considering a compromise between the sampling cost and actual variance reduction.

For a more elaborate and adaptive solution, check [Agarwal et. al, 03].


Is there way for polygon decomposition by thickness?

It is a bad idea to try to abuse clustering algorithms for this. In particular k-means will not cut vertically, nor consider rectangles.

What you need to do is a trivial axis-aligned corner detection.it does not get much simpler than that!

You begin on the left, and whenever the location of the topmost pixel changes, a rectangle is complete.


When performing the unfolding the way you want, after you have unfolded one triangle in the plane, the next one should be one of the three triangles, adjacent to the unfolded one, i.e. the next triangle to be unfolded shares a common edge with the unfolded one.
Assume you have a triangle $ABC$ in three space with vertices  $$A = [A_1, A_2, A_3], \,\, B = [B_1, B_2, B_3],\,\, C = [C_1, C_2, C_3]$$ which is one of the triangles from the 3D mesh.
If $ABC$ is the first triangle you are starting the unfolding with, then:

calculate $l_{AB} = ||A - B|| = \sqrt{(A_1 - B_1)^2 + (A_2 - B_2)^2 + (A_3 - B_3)^2}$;
set the 2D vertices $a = [0,0]$ and $b = [l_{AB}, 0]$ as well as $$l_{ab} = l_{AB};

Else if $ABC$ is a next triangle to be unfolded in the plane, i.e. you have already unfolded one of $ABC$'s adjacent (neighboring) triangles, therefore one of $ABC$'s edges is already unfolded in the plane. Say that this edge is $AB$ and it is unfolded in the plane as the edge with vertices $a = [a_1, a_2]$ and $b = [b_1, b_2]$, where $l_{ab} = ||a-b|| = \sqrt{(a_1 - b_1)^2 + (a_2 - b_2)^2} = ||A-B|| = l_{AB}$.
Either way, for the triangle $ABC$ you already have the edge $AB$ unfolded as $ab$ in the plane, so that $l_{ab} = l_{AB}$ (i.e. you already know $a, \, b$ and $l_{ab}$). The question is, how do you find the coordinates of the third vertex $c = [c_1, c_2]$ so that $l_{bc} = l_{BC}$ and $l_{ca} = l_{CA}$? We impose the latter two conditions because two triangles that have pairwise equal edge-lengths are congruent (isometric, i.e. identical copies of each other).
For your information, I will write vectors like $B - A = [B_1 - A_1, \, B_2 - A_2, \, B_3 - A_3]$ and $C - A = [C_1 - A_1, \, C_2 - A_2, \, C_3 - A_3]$ and $b-a = [b_1 - a_1, b_2-a_2]$ and will use dot product $\cdot$ and cross product $\times$ between vectors. You already know what $||B-A|| = l_{AB}$ and $||C-A|| = l_{CA}$ are and that $l_{AB} = l_{ab} = \|b-a\|$.

Calculate $$s = \frac{\|(B-A)\times(C-A)\|}{l_{ab}^2}$$

Calculate $$c = \frac{(B-A)\cdot(C-A)}{l_{ab}^2}$$

Calculate both pairs of coordinates
\begin{align}
&c'_1 = a_1 + c(b_1-a_1) - s(b_2 - a_2)\\
&c'_2 = a_2 + c(b_2-a_2) + s(b_1 - a_1)
\end{align}
and
\begin{align}
&c''_1 = a_1 + c(b_1-a_1) + s(b_2 - a_2)\\
&c''_2 = a_2 + c(b_2-a_2) - s(b_1 - a_1)
\end{align}

You have two candidate points $c' = [c'_1, c'_2]$ and $c'' = [c''_1, c''_2]$

Perform a check whether $c'$ is on the same side of the edge $ab$ as the  triangle unfolded right before $ABC$.

If it is not, then $c = c'$
else if it is, then $c = c''$
As a result of this procedure, you have constructed the planar unfolding $abc$ of triangle $ABC$ and the next triangle from the 3D mesh to be unfolded (if any), should be the one that shares either edge $BC$ or $CA$ with $ABC$. Either way, for that next triangle, you already know the unfolding $bc$ of the edge $BC$ and the unfolding $ca$ of $CA$ (whichever is relevant), so you can apply the same procedure for that next triangle. And so on and so forth, until you have unfolded all triangles from the mesh.
Have in mind that there are 3D meshes that do not posses unfolding that do not self overlap, no matter what you do. For such meshes there no greedy algorithms or any other algorithm that can do it.


Possible methods to deal with artifacts when drawing billboards to a cubemap

Change the projection to always be the same no matter which face is being rendered.

The simplest way is you use the direction from billboard to the camera instead of the direction the camera is facing to rotate the billboard.


How to voxelize geometry using textures and shaders?

A fragment shader executes for a 2D rasterized position. When you output the normalized position it's still being projected onto the 2D xy image plane onto a render target. The render target is 2D, so you can't use the fragment shader output to achieve this.

You need to bind a 3D texture as a read/write texture in order to store values.

The other thing is how you're expecting to voxelize the geometry, are you trying to voxelize the surfaces? If so you can use a pixel shader but you need to render all faces (disable back face culling).


What actually happens in frame buffer when an image or document is scrolled?

In the olden days of text-mode terminals, the "frame buffer" was a buffer of text, not an image. Scrolling worked in the way you described. The buffer was larger than would fit on the screen, and to "scroll" it would change an offset in the display controller, so it would display different lines. This only makes sense when the whole framebuffer is scrolled at once, which never happens in a windowing system.

Nowadays, it depends on the application. A word processor will rasterise the text separately each frame: there's no off-screen buffer at all. It wouldn't make sense to store an image of the whole document, because the memory use would be huge, and the whole thing would have to be recreated every time the text changed.

On the other hand, an image editor will probably store the whole image, possibly as a texture in GPU memory. A program on the GPU copies the visible section of the image to the application's frame buffer, each frame. This can be done by simply having the view parameters (offset and zoom level) as uniforms and updating them as the view is changed. This is actually very fast even though it's essentially a copy, because GPUs are good at drawing textures efficiently.

A web browser might use either technique, or something more sophisticated. There was a trend in mobile browsers a couple of years ago to render in tiles. The structure of the page would be rasterised to a set of off-screen tiles: enough to more than cover the screen. A program on the GPU draws all the visible tiles into the correct place in the application's framebuffer, clipping the edges. That way, small scrolls or zooms can be performed by changing the offsets of the tiles, without having to rasterise the page again. When you scroll too close to the edge of the tiles, a new row of tiles can be rasterised in a separate thread.

This tile approach is fast to display, uses a fixed number of fixed-size textures, and allows a trade-off between responsiveness and scroll speed. It also allows lower-resolution tiles to be used if the page is too slow to rasterise during a fast scroll. It does require a lot more sophistication in the way you rasterise the page to the tiles, to know which parts of the page can affect which tiles.

Whichever technique the application uses for managing off-screen content, it's responsible itself for producing a new frame each time. The windowing system or compositor doesn't give it a way to specify a window of a larger texture to use without copying.


Attribute Location in Multiple Shader Programs

Each VAO, and shader for that matter, has a separate space of attribute locations. It's perfectly fine to have multiple shaders use the same attribute locations. Indeed, it's a good idea to do so if those attributes conceptually mean the same thing to both shaders.

A VAO can be used with a shader if that VAO's attribute locations match those the program expects. Well, technically they don't have to match. It's OK if the VAO provides an attribute that a shader doesn't use. It's more complicated if a VAO doesn't provide an attribute for a location that the shader expects.

So whatever it is that generates your shader prefix strings and VAOs should not be based on what has been generated for other shaders/VAOs. So if you're going to have attributes increase numerically starting from 0, then this should be the case for each independent shader. It should not be based on some constantly increasing global.


Based on the article, it seems like, when it adds "directional lighting information", the program tries to infer a location of a light source and base the normals on the colors with the added information of the directional lighting from the pixel to the inferred light source. I feel like I'm making no sense right now, so here's an example using their pictures:

As you can see, on the far right of the normal map, it's very pink, just like how the sample image is completely white. I would assume this is the case because the inferred light is directly on top of it, with all the normals (in sort of a curved fashion, resulting in the gradient) pointing right at it. In contrast...

This normal map isn't trying to figure out a light source, and instead just kind of uses the sample texture as a heightmap of sorts.
Late reply, but hope this helps! :)


Unsized arrays aren't implemented in HLSL. You can pass fixed-size arrays into functions like
void foo(float2 x[2][4])
{
    ...
}

(Note the brackets go after the parameter name, not after the type—just like declaring a variable.)
However, having different array sizes that can be passed into the same function won't work in HLSL the way it would in C/C++. For variable-size data, better to use a buffer resource, such as a constant buffer, structured buffer, etc (depending on what you're doing with this data). Or if it's in groupshared memory, you'll have to specify a maximum size and allocate the array for that.


What method is used for baking grayscale curvature maps

This is how I ended up doing it:
Instead of using vertex colors like most of the solutions I found online, I calculated the angle of each edge, and then used a shader to draw the lines in UV space. This Also gives me control over how thick the lines should be.
In some cases the angle calculated isn't right, but that should be easy to fix.

Here is the project: https://github.com/Jummit/curvature-baker


I'm not an expert on path tracing, maybe someone with more experience can come up with a better answer, but I'll have a go at answering your questions.

Firstly, for the purposes of a college-project scale path-tracer, VRAM capacity shouldn't really be an issue. A vertex with position, normal, tangent and UV fits in 48 bytes, so with a GB of buffers you can hold ~20 million vertices, not counting the space for your index buffers and kd-trees of course. Even my crappy 5-year-old midrange laptop card has 2GB, and modern cards have significantly more, so i doubt you're gonna exhaust VRAM with model data.

As for the local memory question, short answer is no. Group shared memory for shader cores is tiny, you shouldn't be allocating more than a few KB per shader group, so there's no way you'd get your acceleration structures in there.

One of the main issues with path tracing is coherence, which is what makes it difficult to do on GPU. Since ray directions are chosen at random, each group of shaders processing a nearby group of pixels can end up sampling radically different parts of the scene, and hence memory. This can introduce a lot of branching and cache misses if not handled carefully, and negate a lot of the advantages of running on GPU.

The current state of the art for production path tracers is probably Disney's Hyperion. This is a CPU path tracer that runs on absolutely enormous scenes (think San Fransokyo in Big Hero 6), that are potentially too big to fit into even main RAM, and need to be streamed in from disk on demand ('out-of-core' rendering). They mitigate the coherence problems by first generating rays en-masse, then sorting them into bundles with similar position and direction, so that consecutive ray solves are likely to touch the same part of the scene and utilize the cache. Even so, they're bound by memory and I/O speeds, not compute power. Good write up on it here including a technical PDF.

For your second question, that is more or less the way it's done. You just accumulate the light from consecutive passes onto your colour buffer, and divide by the number of passes to normalize - you probably want to use 32-bit float colour channels for this. Naive path tracers (which is probably what you're building to start with) just choose a truly random sample distribution and weight each sample by its probability. Smarter approaches bias the sampling towards the more significant directions then remove the bias when summing the color. For example, for a diffuse surface, you could (a) pick truly random directions in the hemisphere and weight each sample by $N.L$ or (b) sample in a cosine distribution (more samples around the normal and less around the horizon) and apply full weight to each sample.

If you haven't seen it already, Ray Tracey's blog has some good tutorials on optimizing path-tracing for GPU. Good luck with your project.


How to convert from glBegin() and glEnd() to VBOs?

Back when I did the transition from immediate mode to vertex buffer object (VBO) in my own code base, I have found it easier to proceed with these steps:

From naive immediate mode code...

The original code was something like this:

glBegin(GL_TRIANGLES);
glNormal3f(nx, ny, nz);
glVertex3f(px, py, pz);
// [...] arbitrary long sequence of such function calls.
glEnd();

1. To array based immediate mode code...

First, keep using the immediate mode, but make sure that all values passed to the vertex functions come from an array. This seems to be already the state of the code described in the question.

float vertexData[] = { ... }; // Can be all together, or split by attribute type, as you prefer.
glBegin(GL_TRIANGLES);
for (int i = 0; i < numberOfVertices; ++i)
{
    // The indices are going to depend on the data layout.
    glNormal3f(vertexData[0], vertexData[1], vertexData[2]);
    glVertex3f(vertexData[3], vertexData[4], vertexData[5]);
}
glEnd();

2. To vertex array version...

Second, you can replace the loop by a handful of calls with vertex arrays.

glEnableClientstate(GL_Vertex_Array);
glVertexPointer(3, GL_FLOAT, 6 * sizeof(float), vertexData);
glNormalPointer(GL_FLOAT, 6 * sizeof(float), vertexData + 3);
glDrawArrays(GL_TRIANGLES, 0, numberOfVertices);
glDisableClientState(GL_VERTEX_ARRAY);

3. To vertex buffer version

Finally from the vertex array based version, the transition to VBO is minimal:

// An initialization part will create and populate the vertex buffer:
GLuint vboId = 0;
glGenBuffers(1, &VBOid);
glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, VBOid);
glBufferData(GL_ELEMENT_ARRAY_BUFFER, 6 * sizeof(float), vertexData, GL_STATIC_DRAW);

// A draw part will bind the vertex buffer and issue the drawcall.
// It's almost the same code as before, except:
// - We bind a VBO
// - We pass (0+offset) as the last parameter for the gl*pointer().
glBindBufferARB(GL_ARRAY_BUFFER_ARB, VBOId);

glEnableClientstate(GL_Vertex_Array);
glVertexPointer(3, GL_FLOAT, 6 * sizeof(float), (void*)(0));
glNormalPointer(GL_FLOAT, 6 * sizeof(float), (void*)(0 + 3));
glDrawArrays(GL_TRIANGLES, 0, numberOfVertices);
glDisableClientState(GL_VERTEX_ARRAY);

And that's it. Both vertex array and VBO allow to work directly with vertices or with vertex index, so be careful about possible confusion there.

P.S.: I have written this answer without double checking the code, so it's probably ridden with typos and not meant for copy-pasting. I invite you to read the documentation, understand the API and refer to other examples when you implement it. At first, VBO can be really annoying to get to work properly.

I have used the same distinction when answering this other question about performance of VBO.


How does Primary Sample Space Metropolis Light Transport work?

PSSMT operates directly on the space of random numbers that generate valid light paths. As such, mutations in the unit hypercube lose their physical interpretation since they do not have direct knowledge of the actual light path constructed. Recent research in rendering has shown that it is possible to bridge the gap between path space (that acts directly on the paths) and primary sample space (that acts on uniform variates) to some extent. See Bitterli 2018, Pantaleoni 2017 and Otsu 2017 if that interests you.
I am not entirely sure what you mean here, but typically an initial path of length $k$ will be formed using path tracing, and thus random numbers will be used along the way ($k+2$ if you are using BDPT). PSSMLT will mutate this path by acting only on the random numbers. Path throughput is indeed independent of the pixel but you only care about paths actually reaching the eye. At the end of the day, you are importance sampling paths according to their contribution to the final image, and luminance is what is commonly used.


Shader cannot move texture x, y location

What you do is shifting the texture coordinates. This will move the texture on your polygon, and depending on your settings it will handle the borders which are now probably wrong. however, this will not move your polygon - if you want to move an entire object, you should not do that within the shader. you do that via modifying your model matrix.

I assume you have a rectangle, and you want to move it by 0.2 along the x axis. Then you would edit the correct element of your model matrix for the rectangle.

//depending on your setup, you may need [0][3] rather than [3][0]
modelMatrix[3][0] += 0.2;
//now upload the matrix to your shader
//the entire object will be moved by 0.2 units

If that is what you wanted, then you need to remove the fragment shader position calculation.

On a sidenote: I had problems with using += in a shader the other day, so it might just be that this doesn't compile. Try using position = position + 0.2; instead.


How can I send custom bitmap to frame buffer?

This answer uses GDI and a little bit of MFC. The MFC bits aren't necessary and just there for my convenience. You can hit the Windows API directly if you need to.

In your window class have some data members for the bitmap. I currently use MFC's CBitmap for convenience:

CBitmap             bitmap_;
unsigned char*      bitmapBits_;

bitmapBits_ is a pointer to the bitmap data that you can write directly into.

Before rendering, you'll need to create the bitmap using CreateDIBSection() and get the pointer to the bits of the bitmap:

BITMAPINFO bmi;
memset(&bmi.bmiHeader,0,sizeof(bmi.bmiHeader));
bmi.bmiHeader.biSize = sizeof(BITMAPINFOHEADER);
bmi.bmiHeader.biWidth = widthOfTheImage;
bmi.bmiHeader.biHeight = heightOfTheImage;
bmi.bmiHeader.biPlanes = 1;
bmi.bmiHeader.biBitCount = 24;
bmi.bmiHeader.biCompression = BI_RGB;

void* bitmapBits;
bitmap_.Attach(CreateDIBSection(0,&bmi,DIB_RGB_COLORS,&bitmapBits,0,0));
if(bitmap_.m_hObject == NULL)
    AfxMessageBox("CreateDIBSection() FAILED!\n");
bitmapBits_ = reinterpret_cast<unsigned char*>(bitmapBits);

When writing to the bitmap bits, you first need to determine the row size (see here).

// Compute the row size of a DIB.
// (See https://en.wikipedia.org/wiki/BMP_file_format#Pixel_storage and ComputePitch() from atlimage.h)
const unsigned int rowSize = ((24 * imageWidth + 31) / 32) * 4;

Then you can write the colour information into the bitmap's data. Note this is a bottom-up DIB; memory earlier in the buffer represents scanlines lower in the image. Another way of putting it is that the the first element of the bitdata relates to the bottom-left. Colour components are in BGR order.

unsigned char* dibPtr = &bitmapBits_[pixelY * rowSize + pixelX * 3];
dibPtr[0] = ???; // blue between 0 and 255
dibPtr[1] = ???; // green between 0 and 255
dibPtr[2] = ???; // red between 0 and 255

Sorry if the previous part is already known to you, as you asked simply how to get it painted to the screen. Here's my simplified paint routine from my MFC application:

void RendererWnd::OnPaint()
{
    CPaintDC dc(this);

    CDC dcMemory;
    dcMemory.CreateCompatibleDC(&dc);
    CBitmap* oldBmp = dcMemory.SelectObject(&bitmap_);
    dc.BitBlt(0,0,imageWidth,imageHeight,&dcMemory,0,0,SRCCOPY);
    dcMemory.SelectObject(oldBmp);
}

There are going to be many other ways of doing this. If you're doing a long render and want to update the window periodically then you'll probably want to have a timer in your window. When this timer fires you want to request a repaint of your window using InvalidateRect(). You'll need to synchronise access to the bitmap data between your rendering thread(s) and the thread that's doing the painting. This will ensure the painter can access the bitmap when it's not being written to. I use a CRITICAL_SECTION for this.


Just to comment on the 0.16*Specular^2 term that was mentioned in the comments by Karim:

Frostbite only remaps the specularity for their internal purpose so that they can pack gemstones specularity into the 8bit channel. They just lose some precision in favor of more variation. The specularity value is the original one when rendering.

Check the Disney BRDF viewer which is a awesome start for PBR, since it contains the whole BRDF from Disney.
BRDF Explorer | Walt Disney Animation Studios

Also, if you are looking at a better Fresnel function than

vec3 Light_F(in vec3 Specular, in float VH) {
    // Fresnel Schlick
    return Specular + (1-Specular) * pow(1-VH, 5.0f);
}

Use this instead since it takes roughness/gloss into account:

vec3 Light_F(in vec3 Specular, in float VH, in float gloss)
{
    return Specular + (max(vec3(gloss), Specular) - Specular) * pow(1-VH, 5.0f);
}


Mimicking Blender's "roughness" in my path tracer

For explicit light sampling: yep, you just evaluate the BRDF for that incoming direction and the output direction back toward the camera. There's probability involved in the case of an area light: you have to randomly choose a point on the light source, and include a factor to convert from probability density over light source area to a density over solid angle at the receiver point. It looks like you've got that already (that's the $\cos \theta' / (p(\mathbf{x}')\|\mathbf{x}-\mathbf{x}'\|^2)$ in the formula you mentioned).

For sampling the reflection ray, you could just pick a random ray from the hemisphere and weight it by the BRDF, the same as you'd do for diffuse. But that won't be very efficient; it'll take a huge number of samples to converge to a low-noise result.

What you really want to do is importance-sample the BRDF. This means choosing a random ray not uniformly over the hemisphere, but using a probability distribution that matches the shape of the BRDF (as much as possible). Then you'll automatically get rays that are clustered around the sharp reflection peak in the BRDF when the roughness is low, and spread out more over the hemisphere when the roughness is high.

So how do you do it? The paper Microfacet Models for Refraction through Rough Surfaces by Walter et al (which popularized the GGX distribution in graphics) handily gives importance-sampling formulas for several common BRDFs in section 5.2. Their approach is to importance-sample the normal distribution of the BRDF, then leave the Fresnel and geometry factors to be factored into the path weight.

For example, their GGX importance-sampling formula is (equations 35–36):
$$\begin{aligned}\theta_m &= \arctan \left( \frac{\alpha_g\sqrt{\xi_1}}{\sqrt{1-\xi_1}} \right) \\ \phi_m &= 2\pi\xi_2\end{aligned}$$
where $\alpha_g$ is the GGX roughness parameter, and $\xi_1, \xi_2$ are a pair of uniform random numbers in [0, 1]. The outputs $\theta_m, \phi_m$ are the polar coordinates (relative to the surface normal) of the sampled microfacet normal.

So to use this, you'd choose random inputs $\xi_1, \xi_2$, then evaluate the above equations, and use the resulting normal to compute the reflection ray, as well as the Fresnel and geometry factors from the BRDF.


I don't think having a equation for flat shading is meaningful as such. Mainly because the formula for flat shading is same as the formula for the model the flat shading is emulating. The difference is that flat shading is operating per primitive rather than per point or per vertex.

So there can in fact exist multiple formulas that work differently but are still all flat shading.


Raytracing in OpenGL using Compute Shaders

So i searched a lot after this and I think it was my confusion on FrameBuffer Objects.

I thought you could use FBO's just like a default FrameBuffer and display the texture image attached to it but you can't. It's only used for offscreen rendering. So while you can use rendering commands to draw something to a "texture image" attached to it, you can't "display the image" by making it default framebuffer or something like that.


What subdivision algorithm advances have occurred since Catmull-Clark?

More of an extended comment than an answer:

What do you mean by "optimization and accuracy"? Do you mean computational efficiency for some particular application, like ray tracing, physical simulation, CAD modeling, ....?

Note that the idea of "accuracy" for a subdivision scheme is not well-posed. Different schemes will have different limit surfaces, but there is no canonical way of declaring one limit surface to be "more accurate" than any another. One can pose some constraints on the type of limit surface one desires, but these constraints are again very application-dependent: one person might ask for $G^n$ everywhere, the next will complain because this precludes preserving sharp creases.

Catmull-Clark (and Loop, for triangle meshes) remains popular because of its simplicity, which in many cases outweigh its weaknesses (no handling of sharp features; loss of regularity at extraordinary vertices). Countless alternative schemes (which may or may not be improvements over Catmull-Clark, depending on the specific application) have been proposed -- if you have a specific application in mind, with specific requirements, we may be better able to help you navigate your options.


Conversion from cubic catmull-rom spline to cubic b-spline

Quoting the comments above for context:

Just to confirm, are you asking, given a set of $N$ CatRom control points,
$$\{CR_0, CR_1, CR_2, CR_3 ... CR_{n-1}\}$$
forming a piecewise curve, what is the equivalent $N$ points, $$\{B_0,B_1,B_2...B_{n-1}\}$$ for a matching piecewise uniform cubic bspline?

Exactly. That is what I'm looking for. I think conversion to cubic bezier's are easy. But I would like to know if this is even possible/doable with bsplines

I think it's a bit complicated: Using the approach you described above, i.e. $$ V_{targetformat} = M_{targetformat}^{-1} . M_{sourceformat}.V_{sourceformat}$$
you can easily go from a sequence of $n$ control points of Catmull-Rom to  $(n-3)$ sets of 4 Uniform BSpline control points, i.e. a total of $4n-12$ control points , but I can't see that this will allow you to stitch them together to obtain just $n$ points...
Assuming I've got the correct matrices for the 'typical' Catmull-Rom,
$$M_{CR}=\frac{1}{2}\begin{bmatrix}
-1 & 3 & -3 & 1\\
2 & -5 & 4 & -1 \\
-1 & 0 & 1 & 0 \\
0 & 2 & 0 & 0\\
 \end{bmatrix}$$
and uniform B-Spline:
$$M_{UB}=\frac{1}{6}\begin{bmatrix}
-1 & 3 & -3 & 1\\
3 & -6 & 3 & 0 \\
-3 & 0 & 3 & 0 \\
1 & 4 & 1 & 0\\
 \end{bmatrix}$$
then $M^{-1}_{UB}.M_{CR}$ should be
[UPDATE: Fixed errors in coefficents]
$$\frac{1}{6}\begin{bmatrix}
7 & -4 & 5 & -2\\
-2 & 11 & -4 & 1 \\
1 & -4 & 11 & -2 \\
-4 & 13 & -26 & 11\\
 \end{bmatrix}$$
It seems to me that, in any batch of four points, every output vertex thus depends on every input vertex, so it seems the neighbouring points can't be shared.
Having said this, I think perhaps there might be a way of getting an approximate solution (but I've not tried rendering it myself to check it so take with a grain of salt)
We know that a Catmull-Rom spline defined by $(CR_A, CR_B, CR_C, CR_D)$, passes through $CR_B$ and $CR_C$ at t=0 and t=1 respectively.
Similarly, using $\begin{bmatrix}t^3 & t^2 & t & 1\end{bmatrix}.M_{UB}$, we should see that at t=0 and t=1, a B-spline passes through $\frac{B_A + 4 B_B + B_C}{6}$ and $\frac{B_B + 4 B_C + B_D}{6}$.
First consider trying to go from BSpline to CR. If we had a 'loop' formed by N control points one would map from BSpline to CR with...
$$\begin{bmatrix}CR_0\\CR_1\\ ...\\CR_{n-1} \end{bmatrix} = M_{NxBtoCR}.\begin{bmatrix}B_0\\B_1\\ ...\\B_{n-1} \end{bmatrix}$$
where
$$
M_{NxBtoCR}= \frac{1}{6}.\begin{bmatrix}
4 & 1 & 0 & 0 & ... & 0 & 1\\
1 & 4 & 1 & 0 & ... & 0 & 0\\
0 & 1 & 4 & 1 & ... & 0 & 0\\
: & : &   &   & ... & : & :\\
1 & 0 & 0 & 0 & ... & 1 & 4\\
\end{bmatrix}
$$
If we invert that matrix to get $M_{NxBtoCR}^{-1}$, one should then be able to map the other way.
Using a maths package such as Maple, one finds that each row of $M_{NxBtoCR}^{-1}$ consists of the same set of values, all non-zero, which are progressively offset in each row. If we ignore the smaller terms, it suggests that $$B_i \approx 1.732 CR_i -0.464 (CR_{i-1}+CR_{i+1}) + 0.1243 (CR_{i-2}+CR_{i+2}) - 0.033(CR_{i-3}+CR_{i+3}) $$
There's probably some nice algebraic expression for this but I hope this might "do the trick".
Update Those weights don't quite sum to 1.0 - you'll need to normalise them slightly.
Obviously, this has assumed the curve is a loop. It's probably possible to adjust the weights for the end points but I'll leave that as an exercise for the reader.


OpenGL GLSL - Sobel Edge Detection Filter

The best results strongly depend on your use case. They also depend on what effect you want to achieve. Sobel is just an edge detection filter: the edges will depend on the input signal, choosing that input signal is up to you.

Here you are using the color image as an input, and the filter rightfully detects faint edges in the blue gradient, while the edges of the cube get interrupted where its color is too close from the background color.

Since I assume your program is also responsible for drawing the cube, you have access to other information that you can feed to your Sobel filter. For example depth and normals are good candidates for edge detection. The albedo before lighting could be used to. Test with different inputs and decide which ones to use depending on the results you get.
Regarding your question about how to combine the edge information, I suggest you filter g a little before using it. Then you can use it to interpolate between the original color and the wanted edge color.

For example you could try something like this:

    float g = sqrt(pow(gx, 2.0)+pow(gy, 2.0));

    // Try different values and see what happens
    g = smoothstep(0.4, 0.6, g);

    vec3 edgeColor = vec3(1., 0., 0.2);
    color = vec4(mix(diffuse, edgeColor, g), 1.);

Addendum

To use depth or normals, you would need to save them in a texture if that's not done already. When you create the frame buffer for your regular rendering pass, you can attach various textures to it (see glFramebufferTexture2D) and write other information to them than just the color of the scene.

If you attach a depth texture (with GL_DEPTH_ATTACHMENT), it will be used automatically for depth. If you attach one or several color textures (with GL_COLOR_ATTACHMENTi), you can write to them by declaring several outputs to your fragment shaders (this used to be done with gl_FragData; either way, see glDrawBuffers).

For more information on the topic, look up "Multi Render Target" (MRT).


How can I bind my texture in the correct orientation?

Binding the image is unrelated to its orientation. A shot in the dark is that your texture coordinates are either reversed or computed with the orientation in the wrong direction. You could just rotate the image as a work around. The main issue here is whether the upper left corner is your origin, or the lower left corner is the origin. But neither of these explain the image being sideways. At any rate, look at your texture coordinates, how the image is being copied into memory, how is is being passed to opengl, and the texture coordinates.


I think that it is important to classify between different types of animations here. One rough classification could be

Rigid body animations
Vertex skinning animations
Others...

Rigid body animations

This refers to animations where the actual geometry data does not change. For example, you could have a geometric object, like a unit cube, and apply (usually affine) transformations to it: You can rotate, translate or scale this object. This is usually accomplished by some sort of scene graph: Each node of this scene graph can have information about the rotation, translation and scaling that should be applied to the objects. In general, this is something that can be expressed with a 4x4 transformation matrix.

The key point is that you don't have to store "animated geometry data". You only have to store the geometry data once, in its original state. In addition to that, you store the key frames of the transformation.

In your case, if you wanted to do a rigid body animation, you could store the geometry in an OBJ file, and additionally store the key frame information in an arbitrary form, for example

TIME  TRANSLATION
0.0   ( 0, 0, 0)
0.1   (10, 0, 0)
0.2   (10,10, 0)
0.3   ( 0,10, 0)
0.4   ( 0, 0, 0)

This could describe a movement of the object along a rectangular path, taking 0.4 seconds. Fortunately, you don't have to create your own file format for that. There already are solutions for storing this sort of animation data. For example, the COLLADA format already supports animations like this. A more compact representation of such an animation can be stored in a glTF asset. For example, it can represent animations like this one, from the glTF sample models repository:

Vertex skinning animations

This technique refers, roughly speaking, to cases where you have an "initial" geometry, and you can define "rules" of how this geometry should be deformed based on the pose of a skeleton. As the name suggests, this is mainly intended for character animations, where you define a skeleton consisting of the "bones" of the character, and the geometry defines the "skin" of the character. Vertex skinning is also supported by different file formats like COLLADA or glTF. The effect is shown in this image, again, from the glTF sample models:

As already mentioned by joojaa in his answer, implementing vertex skinning according to the specifications of these file formats may be a bit tricky. Additionally, a considerable amount of the skinning computations is usually done by the vertex shader. So for a ray-tracer, you would have to move these computations into your application code.

Others...

This refers to what you already alluded to in the question: There may be cases where the animation is not a simple transformation of rigid objects, and where you cannot define any "rules" of how the vertices should be transformed. For example, when the vertex positions are computed by a simulator, like a fluid simulator or a soft body simulator.

In doubt, it may simply be necessary to actually store the positions of all vertices, frame for frame. You could achieve this as a sequence of OBJ files. But note that OBJ files contain ASCII data, and thus, are usually large and the geometry data has to be parsed before it can be used. Storing the data in a binary form that can be used by the application directly may be beneficial here.

There are some possible approaches for reducing the amount of data that has to be stored. For example, one could extend the concept of "key frames" here. For example, one could store the initial configuration of al vertices, and then only store the movement of vertices that acutally have changed between two frames. But of course, if (nearly) all vertices changed, the savings may be diminishing.

I'm not aware of any universally applicable, "standardized" solutions for this. Maybe some of the solutions that joojaa referred to in his answer could be helpful here.


Normalize $v_a$ and $v_b$ and your algorithm should no longer depend on the length of the edges, unless you consider numerical stability.

As for the "best way", I'm not sure there's a much better algorithm but I would try a better definition of the tolerance. Instead of a fixed value, I would use differentials to compute the change in the angle for a change of the inputs of a few ulps. If the resulting range overlaps 0, then I think the points are colinear within available precision.

Whether that complexity is justified depends on what you're doing with this, which you did not specify. A constant could be fine or even the best choice for some uses.


Controlling reflection and refraction with material properties in ray tracing

It sounds to me like you're slightly mixing up several concepts: how lighting due to point-light sources (Phong model) and lighting due to other objects in the scene (recursively raytraced reflection or transmission) combine with each other, as compared to how specular lighting combines with diffuse lighting and/or transmission (using the Fresnel equations).

A useful way to think about specular vs diffuse is that specular happens at the surface of the material while diffuse, transmissive, and refractive properties are associated with the bulk (the interior) of the material. So you can think of it as two layers composited together. Specular is on top: it reflects a portion of the incoming light, defined by the Fresnel equations, and any remaining light proceeds into the material to light the bulk.

The bulk will have some combination of scattering and absorption properties. Low scattering makes a transmissive material and high scattering makes a diffuse material (the diffuse color is due to scattering on a scale too small to see). Some materials have an intermediate degree of scattering too, but we can ignore that for simplicity and just treat it as a boolean: opaque (Lambert diffuse) or transmissive.

Meanwhile, each of these potential material components (specular, diffuse, and transmissive) has to deal with light arriving both from point light sources and from the rest of the scene. Physically, all light sources are just emissive objects that are part of the scene (area lights) but again for simplicity we often approximate by using point lights.

So, putting it all together, when lighting a surface:

For each point light

Calculate Fresnel factor based on light vector
Accumulate Phong specular weighted by Fresnel reflectance
[if opaque] Accumulate Lambert diffuse weighted by Fresnel transmittance
[if transparent] Nothing; light goes through and isn't scattered into the camera.

For reflecting the rest of the scene

Calculate Fresnel factor based on view vector
Cast a reflection ray and accumulate it weighted by Fresnel reflectance
[if opaque] Nothing. (Ideally, accumulate Lambert diffuse due to all objects visible in the normal hemisphere of the surface, i.e. indirect diffuse. However, for a basic recursive raytracer this isn't feasible; it requires path tracing, photon mapping, or suchlike.)
[if transparent] Cast a transmission ray and accumulate it weighted by Fresnel transmittance

Also, a note about your material properties: physically, not all of those properties are actually independent of each other. Ks, reflectivity, and refractive index are all the same thing. Ks is the same thing as reflectivity (literally the same value) because both control specular reflections; Ks is for point lights and reflectivity is for reflected environment rays, but it's the same physical mechanism for both. And Ks is determined in turn by the Fresnel equations based on refractive index (which can vary with wavelength).

So, strictly speaking, you should pick one of those variables and let it determine the others. For opaque materials where you don't really care about refractive index, you might let users set Ks (at normal incidence) and backsolve the Fresnel equations to get the appropriate refractive index from that. For transparent materials like glass you can let users set the refractive index and then let the Fresnel equations do their thing.

BTW, to get some physics background that will hopefully help clarify all this, I highly recommend Naty Hoffman's "Physics and Math of Shading" talk from the SIGGRAPH 2015 Physically-Based Shading Course. It gives a quick and not-too-mathematical intro to the physical basis of things like diffuse and specular lighting, BRDFs and so on.


Are these edges a result of display color accuracy, computational accuracy, or something else?

Updated answer, now that pow is getting valid arguments:

What you are seeing is called "banding" and it comes from the fact that the color channels are quantized into 8 bits - in other words, it comes from the fact that there are only 256 different shades of green, even though the math generating the colors is capable of much larger resolution.

The banding on my system isn't as bad as what you are describing but it's still present.  The reason it isn't as bad on my system is because you and I have different displays.

Some computer displays don't actually even have a full 256 different shades of green.  For instance, instead of having 8 bit color, they may only have 6 bit color, and they either dither it, temporally dither it, or neither and just let it look worse.

One way to get around this problem is by dithering.  There are various algorithms for dithering, with different computation and quality trade offs.

A good modern dithering algorithm is called "interleaved gradient noise" which is fairly close to blue noise - which is the ideal - and is pretty cheap to calculate.

Check out this link for more detailed information on different types of dithering, and how they compare:

https://bartwronski.com/2016/10/30/dithering-part-three-real-world-2d-quantization-dithering/

First answer, from when pow was getting a negative value parameter:

The issue is that you are passing a negative value to pow for the first argument, which has undefined behavior.

What shows up on my screen is completely different than what shows up on yours.

You shouldn't rely on the behavior of anything in the red pixel areas in this image:

More info:

https://www.opengl.org/sdk/docs/man4/html/pow.xhtml


How can I debug what is being rendered to a Frame Buffer Object in OpenGL?

Generally to see what is being rendered in the various steps of your pipeline I'd suggest the use of a tool for frame analysis. These usually provide you with a view on the content of each buffer for each API call and this can help you in your situation.

A very good one is Renderdoc which is both completely free and opensource. Also it is actively supported.

Another one is Intel GPA that unfortunately, according to its webpage, support up to OGL 3.3 Core.

Just for the sake of adding one more, I used to use gDEBugger, but it has passed long since the last update. This has evolved into AMD CodeXL that unfortunately I never used.


This question is a bit too broad to answer, but I'll try.

Right now you have a screenshot and you are not really sure whether it's correct or not, and if it turns out it is not, it will be difficult to know what is incorrect about it. I would recommend to split your problem into smaller, isolated, problems that can be validated separately.

For example, is the lighting information correct at each point? Are the coordinates (screen space vs world space) correct? To answer these questions you could maybe represent the distance from the light at each point, with something like the following, and check that the visual result is consistent with your scene.

color = fract(length(light.position - currentPosition));

Similarly, is the shadowing correct? This problem is unrelated to scattering so you can check that it  is correct by just having a constant light contribution (your commented code seems to indicate that you have tried so already). You will see quickly if a volumetric shadow is properly cast from occluding geometry.

if (!shadow)
{
    accumFog += light.diffuse;
}

Finally you can test the scattering without any shadow, and make sure it behaves as expected. For example, you should get the same fog color when changing the number of steps. Moreover, the scattering contribution should be less noticeable when there is geometry close to the camera (make sure you stop marching when you reach the scene depth).

Once you are confident the separate parts are working, you can combine them until you reach the final result. I have suggested a few points to check, but you might check other points as well, as you see fit.

  At the end should I just blend the volumetric light buffer with the deferred shading final buffer?

Since the lighting is additive, you just need to do an additive blend of the resulting buffer on your scene. However if you later decide to also handle distance fog (contrast fading in the distance), then a different blending will be necessary.

P.S.: I happen to have written about volumetric lighting in this article (see part "Volumetric lighting"). It doesn't go into details, but maybe it will be helpful nonetheless.


What is a good approach for handling uniforms in modern OpenGL?

Suballocating from a larger buffer is absolutely the way to go, with caveats. I'm coming more from a DirectX/Vulkan side of things, but this should apply equally to OpenGL (I just won't have direct API calls here in this answer). The things to consider are the following:

Do you need to index into the larger buffer, or are you OK with binding the resource to the offset each time?
Have you taken care of any/all alignment restrictions for your uniforms that are packed together (256 byte alignment is common)?

Newer graphics APIs have a "dynamic offset" you can specify with the draw command which is a pretty fast way to indirectly access a subregion of a buffer. However, assuming you are triple-buffering data that's mutable, there should be little to no contention in the driver for binding the data (just some fixed overhead).

In summary, yes, allocating larger regions of memory/buffers and subleasing those regions is considered best practice. This applies even to objects with different shaders (if your allocator can handle it).


Best technique to draw overlapping colored line segments that follow the same route

This question should most probably be asked on GD.SE or UX.SE. These sites specialize in how to design the graphics and how to choose the graphics for your purpose. But since you are here basic options are:

Routes dont overlap but are parallel like electronic layouts (image 1, A), Works well unless you have many overlapping things.
Space interleaved lines, differently colored line regions (image 1, B). Such as dashed lines, however many other interleaved shapes are possible.
Spatial solutions, lines on top of each other but laid over in 3d (image 1, C), Works well in isometric drawing.
Different line widths (image 1, D)
User suggested, overlay (image 1, E).
Same color but labels on segments (image 1, F), note not all segments need to be labeled as humans can interpolate values somewhat. This has benefit of being able to be paired with any other effect and works well with selection highlights with info intact.

Image 1: Different overlay methods.

Note: None, except the labeling works well when there are very many different lines as humans capability to see different colors in this kind of context as differing data is not very good (ask UX.SE).

As for what is best depends. None of these are particularity hard to do in code or vector drawing application (but that is out of scope ask GD.SE). Whichever is best for the user or the design is also not within the scope of this software (ask  GD.SE or UX.SE).


Solving a problem from *Foundations of Computer Graphics*:

Let me try to make my comment into a complete answer.

The general idea is to build a linear system using those 6 point pairs and solve for the desired 12 unknowns. You may find this paper[1] useful if you want to learn the theoretical background, and you can find some more details in my blog post[2].

Here is my solution using Mathematica.

In the first step I will build a symbolic array which contains all variables.

v = Table[Symbol["m" <> ToString[i]], {i, 0, 15}];
M = ArrayReshape[v, {4, 4}]

Next we write down the coordinates for all vertices.

k = f/ n;

X = {
    {l, t, -n, 1},
    {l, b, -n, 1},
    {r, b, -n, 1},
    {r, t, -n, 1},
    {k l, k t, -f, 1},
    {k l, k b, -f, 1},
    {k r, k b, -f, 1},
    {k r, k t, -f, 1}};

U = {
    {-1, 1, -1},
    {-1, -1, -1},
    {1, -1, -1},
    {1, 1, -1},
    {-1, 1, 1},
    {-1, -1, 1},
    {1, -1, 1},
    {1, 1, 1}};

where l,r,b,t,n,and f denote left, right, bottom, top, near, and far planes respectively,
and 'X' is the coordinates in world space, U is the coordinates in NDC.

Then next we implement a function Generate to generate an equation for each pair of coordinates:

Generate[{i_, j_}] := M[[j]].X[[i]] - U[[i, j]] M[[4]].X[[i]] == 0;

Because we do not care about the third row,
we only solve for the unknowns in partialv:

partialv = {m0, m1, m2, m3, m4, m5, m6, m7, m12, m13, m14, m15};

Now we can build our linear system and solve it:

solve[subset_] := Solve[
    Map[Generate, Table[{i, j}, {i, subset}, {j, 1, 2}], {2}] // Flatten,
    partialv] // Simplify

Let me explain it a little bit.

Because only 6 pairs of points are given, we select a subset of 6 elements from the 8 pairs of points.
Then we call Generate function on this subset
and the first two rows of matrix(we do not care about the third row. And the fourth row is the "homogeneous" part, which is also ignored).

Finally we enumerate all subsets from these 8 pairs of points, and try to find a solution:

solutions = Map[solve, Subsets[{1, 2, 3, 4, 5, 6, 7, 8}, {6}]];

We can view one solution by:

MatrixForm[ArrayReshape[v /. solutions[[1]], {4, 4}]]

$$
\left(
\begin{array}{cccc}
 \text{m0} & 0 & \frac{\text{m0} (l+r)}{2 n} & 0 \\
 0 & \frac{\text{m0} (l-r)}{b-t} & \frac{\text{m0} (l-r) (b+t)}{2 n (b-t)} & 0 \\
 \text{m8} & \text{m9} & \text{m10} & \text{m11} \\
 0 & 0 & \frac{\text{m0} (l-r)}{2 n} & 0 \\
\end{array}
\right)
$$

The m0 is the "unknown scalar factor" which is mentioned in problem statement.

In fact if we let m0 equals 2 n/(r - l),
we would get the conventional projection matrix used in OpenGL:

% /. m0 -> 2 n/(r - l)
% // Simplify // MatrixForm

$$
\left(
\begin{array}{cccc}
 \frac{2 n}{r-l} & 0 & \frac{l+r}{r-l} & 0 \\
 0 & -\frac{2 n}{b-t} & -\frac{b+t}{b-t} & 0 \\
 \text{m8} & \text{m9} & \text{m10} & \text{m11} \\
 0 & 0 & -1 & 0 \\
\end{array}
\right)
$$

[1] Heckbert, Paul S. Fundamentals of texture mapping and image warping. MS thesis. University of California, Berkeley, 1989.

[2] Inferring projective mappings Inferring projective mappings | linearconstraints.net


You shouldn’t need to use any trigonometry here at all. If you get the vector from v1 to v2 and divide it by the number of points you want along the line, each subsequent point is v1 + (that vector) × (the index of the point). Works in any number of dimensions, faster than the trigonometric approach.
v1 = Vector(100, 200, 300)
v2 = Vector(500, 400, 300)
steps = 25 # or however many subdivisions you want: int(length(Line(v1, v2))) if you want an equivalent to the original
stepAmount = Vector((v2.x - v1.x) / steps, (v2.y - v1.y) / steps, (v2.z - v1.z) / steps)
points = list()
for i in range(steps + 1)
    x = v1.x + stepAmount.x * i
    y = v1.y + stepAmount.y * i
    z = v1.z + stepAmount.z * i
    points.append(Vector(x, y, z))


Why is accessing textures much slower when calculating the texture coordinate in the fragment shader?

What you're talking about is commonly called "dependent texture reads" in the mobile development community. It's an implementation detail of certain hardware, and therefore it really depends on the GPU as to whether or not it has any performance implications. Typically it's something you see brought up for PowerVR GPU's in Apple hardware, since it was explicitly mentioned in both Imagination and Apple documentation. If I recall correctly, the issue basically came from hardware in the GPU that would start to pre-fetch textures before the fragment shader even began running, so that it could do a better job of hiding the latency. The docs I linked mention that it's no longer an issue on Series6 hardware, so at least on newer Apple hardware it's not something that you have to worry about. I'm honestly not sure about other mobile GPU's, since that's not my area of expertise. You should try and consult their documentation to find out for sure.

If you decide to do some Google searches on this issue, be aware that you'll probably find some older material that's talking about dependent texture fetches on older desktop hardware. Basic in the early days of pixel/fragment shaders, the term "dependent texture fetch" referred to using a UV address that relied on a previous texture fetch. The classic example was bump-mapped environment map rendering, where you wanted to use a reflection vector based on the normal map in order to sample the environment map. On this older hardware there was some major performance implications, and I think it wasn't even supported on some very old GPU's. With modern GPU's the hardware and the shader ISA is much more generalized, and so the performance situation is much more complicated.


Your point is correct when applied to ambient lighting as an approximation technique. This crude approach does indeed give the same lighting level to all surfaces regardless of their surroundings.

However, ambient occlusion is a different (less crude) approach, that models how the light levels in a real scene vary depending on how much reflected light can reach a surface. The confusion is due to two very different approaches having similar names.

Ambient occlusion is a much better approximation of what in real life is described as ambient light - light that arrives from the background rather than directly from a light source.

  If ambient lighting has no fixed position, how can points then be more or less exposed to it?

When used in the context of ambient occlusion, the term "ambient lighting" still means background light approximated as being the same in all directions, but the surface is lit based on how many of those directions are not occluded by objects in the scene.

This gives subtly varying light levels that give much more realistic images than with a single level of ambient lighting. This gives the effects seen in real life such as the corners of a room being slightly darker than elsewhere.


This issue looks like standard shadow map acne artifacts. Additionally your's lighting equation is incomplete or wrong. Light shouldn't influence faces with normals facing away from it. This also means that with a proper equation the "dark" side of the sphere shouldn't have any acne artifacts.

There are three sources of acne artifacts:

First acne source is shadow map precision. Make sure that near and far planes of shadow casting light are as tight as possible. All objects before the near plane can be pancaked, as their exact depth isn't important.
Second acne source is shadow map resolution. For a directional light you should be doing cascaded shadow maps with at least 3 1024x1012 cascades for ~100-200m of shadow distance. It's hard to cover similar shadow distance with one shadow map with uniform projection.
Third acne source is wide shadow map filter like PCF, as using a single depth depth comparison value across a wide kernel is insufficient. There are many methods to fix it, but none of them is robust.

To sum up, a tight frustum with a few cascades and some bias tweaking is enough to get the general case working (directional light). Start tweaking by disabling shadow filtering and tweak for filtering only when basic shadow maps are robust enough.

Additionally apart from the constant depth bias (which you currently use) you should also add slope depth bias and max slope depth bias. Both can be be implemented either as render state or as shader code during shadowmap rendering. Slope depth bias is simply a magic bias value scaled by dot( normal, lightDir ).

There are some additional interesting methods and most of them are implemented in the excellent demo: Matt Pettineo - "A sampling of shadow techniques".

Flip culling during shadow map rendering (trades acne for Peter panning).
Normal offset shadow mapping does wonders for bias issues, but requires to have vertex normals around during shading.
Variance based methods (ESM, VSM, EVSM) completely remove bias issues, but have other drawbacks (light leaking and/or performance issues).


Why do we multiply vertex from left side in vertex shader with matrices?

There is no 100% consensus on what order matrix multiplication should model things, worse the industry is split along this. Some sources use row major and some sources use column major matrices. Great care should be taken to verify which is being used in your source and what is being used by your API!

Some API's go as far as allowing one to multiply vectors on both sides of the matrix. The wisdom is thus in reading your API reference and verifying how it works. Old, fixed pipeline OpenGL and glm, does in fact use column major matrices. While (warning) DirectX had chosen to use row major matrices [source]. Obviously with programmable shaders any application can have any order, even mixed one

Due to fundamental properties of matrices this duality is the same as transposing the matrix. If you transpose the matrix you reverse the order of calculation, as you are now thinking outside in instead of inside out or vice versa.

The moral of the story is that you should always check your assumptions. And just printing out values does not always help as the printing function can screw things up, and memory not aligned the same way. Also simple testing may make you blind to having mixed conventions inside your program. Read what convention is in use! Because, as confusing as having two sources with different  conventions is, it pales in comparasion to a piece of software that mixes the two in the codebase.


Rendering Equation for photons carrying flux

It'd be hell a lot easier if this were on graphicsexchange, since I can't use latex here but anyways.

In the first pass of Photon Mapping you don't need to use the Flux form of rendering equation. You just divide the original flux coming from the light source among the N photons, then for each photon you use Russian Roulette to determine whether it reflects, transmits or gets absorbed and store the incoming flux.

In the second pass you need to find radiance but you have information in Flux form. The rendering equation can be simply transformed.

So all you need to do is divide the flux stored for each photon by the area. Note that the Area is often denoted as A= pi r^2 Since we are concerned with a hemisphere and the underlying area is just a circle.

About your 2nd point, I don't think the BRDF needs to be re-normalized or whatsoever since the equation is essentially the same. Note that the second term in the last equation is still a radiance written in flux form. (BRDF has units of inverse solid angle thus making the whole term radiance).


We had some related questions lately, so you might find them useful:

Paint pixels on the image around the cursor for painting application
How are textures projected onto 3d models in texture painting applications
Get normalized device coordinate of an image pixel

I can't tell you how the named applications do it, since I don't have access to the code, but painting onto a 3d model could be achieved as follows:
When you click the "paint" button, you cast a ray from the camera position into your 3d scene. For perspective projections, the direction depends on the mouse's screen position. You can look into my answer to this question to get an idea of how this is done.
Now perform a collision check of the ray and your model to find the affected triangle/triangles and vertices. Calculate the Barycentric coordinates where the ray hits the affected triangle. Use them to interpolate the texture coordinates of the intersection point. From the model data, you should know, which texture is used at the target location.
Now you know where you want to paint in your models' texture. How to proceed depends on what you want to paint. In case you just want to paint a pixel, you have to overwrite the texture's pixel at this location. You can do this either on the CPU or on the GPU using the rendering pipeline.
If you don't understand some things I mentioned yet, read the links I provided and continue learning OpenGL. Things will get clearer as soon as you know the rendering pipeline and get used to working with textures. If something is still not clear after that, feel free to ask.


An answer with details gotten elsewhere:

It's bad on mobile: ARM doesn't have a tessellation unit so is emulated in software. Tiled based architectures will write vertex data to main memory before pixel shading which is a bandwidth issue. qualcom switches to direct rendering mode if you enable tessellation which isn't cheap.
It's bad on AMD / Console
There is complexity in setting the tessellation factors to give good results
it sounds like it's difficult not to get too much resolution in some ways (aliasing), and not enough in others (low quality results). You essentially need to know screen space derivatives of vertices to calculate good tessellation factors, and this is not easy to get even decent approximations of.
A bummer about tessellation shaders is that they basically work by sampling discrete sample points (at vertices) so for a better result, you need a denser mesh to start with, which sort of defeats the purpose of dynamic tessellation.  Can be addressed a bit with mipmapping type strategies.
Tessellation is still used by major games in specific cases though, and terrain does seem to be one usage case


The open source implementation FreeGLUT has these functions.
The geometry code can be found in freeglut/freeglut/src/fg_geometry.c.
For the dodecahedron, these are the values:
#define DODECAHEDRON_NUM_VERT           20
#define DODECAHEDRON_NUM_FACES          12
#define DODECAHEDRON_NUM_EDGE_PER_FACE  5
#define DODECAHEDRON_VERT_PER_OBJ       (DODECAHEDRON_NUM_FACES*DODECAHEDRON_NUM_EDGE_PER_FACE)
#define DODECAHEDRON_VERT_ELEM_PER_OBJ  (DODECAHEDRON_VERT_PER_OBJ*3)
#define DODECAHEDRON_VERT_PER_OBJ_TRI   (DODECAHEDRON_VERT_PER_OBJ+DODECAHEDRON_NUM_FACES*4)    /* 4 extra edges per face when drawing pentagons as triangles */
/* Vertex Coordinates */
static GLfloat dodecahedron_v[DODECAHEDRON_NUM_VERT*3] =
{
               0.0f,  1.61803398875f,  0.61803398875f,
    -          1.0f,            1.0f,            1.0f,
    -0.61803398875f,            0.0f,  1.61803398875f,
     0.61803398875f,            0.0f,  1.61803398875f,
               1.0f,            1.0f,            1.0f,
               0.0f,  1.61803398875f, -0.61803398875f,
               1.0f,            1.0f, -          1.0f,
     0.61803398875f,            0.0f, -1.61803398875f,
    -0.61803398875f,            0.0f, -1.61803398875f,
    -          1.0f,            1.0f, -          1.0f,
               0.0f, -1.61803398875f,  0.61803398875f,
               1.0f, -          1.0f,            1.0f,
    -          1.0f, -          1.0f,            1.0f,
               0.0f, -1.61803398875f, -0.61803398875f,
    -          1.0f, -          1.0f, -          1.0f,
               1.0f, -          1.0f, -          1.0f,
     1.61803398875f, -0.61803398875f,            0.0f,
     1.61803398875f,  0.61803398875f,            0.0f,
    -1.61803398875f,  0.61803398875f,            0.0f,
    -1.61803398875f, -0.61803398875f,            0.0f
};
/* Normal Vectors */
static GLfloat dodecahedron_n[DODECAHEDRON_NUM_FACES*3] =
{
                0.0f,  0.525731112119f,  0.850650808354f,
                0.0f,  0.525731112119f, -0.850650808354f,
                0.0f, -0.525731112119f,  0.850650808354f,
                0.0f, -0.525731112119f, -0.850650808354f,

     0.850650808354f,             0.0f,  0.525731112119f,
    -0.850650808354f,             0.0f,  0.525731112119f,
     0.850650808354f,             0.0f, -0.525731112119f,
    -0.850650808354f,             0.0f, -0.525731112119f,

     0.525731112119f,  0.850650808354f,             0.0f,
     0.525731112119f, -0.850650808354f,             0.0f,
    -0.525731112119f,  0.850650808354f,             0.0f,
    -0.525731112119f, -0.850650808354f,             0.0f,
};

/* Vertex indices */
static GLubyte dodecahedron_vi[DODECAHEDRON_VERT_PER_OBJ] =
{
     0,  1,  2,  3,  4,
     5,  6,  7,  8,  9,
    10, 11,  3,  2, 12,
    13, 14,  8,  7, 15,

     3, 11, 16, 17,  4,
     2,  1, 18, 19, 12,
     7,  6, 17, 16, 15,
     8, 14, 19, 18,  9,

    17,  6,  5,  0,  4,
    16, 11, 10, 13, 15,
    18,  1,  0,  5,  9,
    19, 14, 13, 10, 12
};

As for the second part of your question, tools such as renderdoc allow you to export meshes which are used in a render. To do this, you must first take a capture, find the corresponding draw call in event browser, open the Mesh Viewer and click on the export icon. You can do this with almost any mesh input to a vertex shader drawn in a supported API.


Quick answers to your individual questions

What do you do if your image is not the same size on each axis though?

The paper uses square images with side lengths that are a power of 2. This is for ease of explanation, but is not necessary for the algorithm to work. See section 3.1:

  Without loss of generality, we can assume n is a power of 2.

That is, this assumption is not required in order for the algorithm to work.

Is it possible to run JFA on each axis separately and still get decent results?

Running on each axis separately is likely to give more incorrect pixel results, and take longer to run in most cases. In extreme cases where one of the image side lengths is less than 8 (the number of jump directions), it may be faster as the algorithm treats those 8 directions sequentially, but for any wider image, separating the axes loses the advantage of treating them in parallel.

In my situation ideally I'm looking to support both single point seeds, as well as arbitrary shaped seeds

The paper mentions arbitrary shaped seeds in section 6 under the subheading "Generalized Voronoi Diagram":

  ...our algorithms treat such generalized seeds as collections of point seeds and thus expect to inherit the good performance obtained for point seeds.

So provided it suits your purpose to model arbitrary shapes as collections of pixels, you don't need to make any adjustment to the algorithm. Simply feed in a texture that labels all pixels in an arbitrary shaped seed with the same seed number, but different locations.

Possibly even weighted seeds, where the distance to a seed is adjusted by a multiplier and/or an adder to give it more or less influence

For "weighting on seeds such as multiplicative and additive", the paper only mentions the possibility in passing in section 8, as potential future work. However, this should be straightforward to implement provided your desired weighting can be included in the seed data that is passed from pixel to pixel.

The current algorithm passes <s, position(s)> to specify a seed and its position, and only one seed is stored per pixel at any one time. Extending this to store <s, position(s), weight(s)> provides all the information required to weight the distance function and calculate whether the new seed being passed to a pixel is closer to it than the one it is currently storing.

You could even include two weights, one multiplicative and one additive, and just set the multiplicative one to 1 and the additive one to 0 when not required. Then your algorithm would include the possibility of being used for multiplicatively weighted seeds, additively weighted seeds, or even a combination of both at once or some of each. This would just need

<s, position(s), multiplicative(s), additive(s)>

and the current algorithm would be equivalent to this new algorithm using

<s, position(s), 1, 0>

Detailed explanation of why

As in the paper, all uses of $\log()$ refer to the base 2 logarithm.

The algorithm does not need to be adapted for different side lengths

If the side lengths are not equal, and are not powers of 2, there is no need to adapt the algorithm. It already deals with pixels on the edge of the image for which some of the jump directions lead outside the image. Since the algorithm already omits the seed information for directions that lead outside the image, a width or height that is not a power of 2 will not be a problem. For an image of width W and height H, the maximum jump size required will be

$$\large2^{\lceil\log(\max(W,H))\rceil-1}$$

For the case of equal width and height N, this reduces to

$$\large2^{\lceil\log(N)\rceil-1}$$

In the case of a side length N that is a power of 2, this reduces to

$$\large2^{\log(N)-1}=N/2$$

as used in the paper.

In more intuitive terms, round the maximum side length up to the next power of 2, and then halve that to get the maximum jump size.

This is always enough to cover every pixel in the image from every other pixel in the image, as the offset to any pixel will be in the range 0 to N-1 if the longest side length is N. Combinations of the powers of 2 from 0 to N/2 will cover every integer up to N-1 if N is a power of 2, and if N is not a power of 2 the range covered can only be larger than required, due to taking the ceiling of the logarithm (rounding up to the next power of 2).

Images with sides not a power of 2 will not be drastically less efficient

The number of jumps depends on the longest side length, say L. If L is a power of 2 then the number of jumps is $\log(L)$. If L is not a power of 2 then the number of jumps is $\lceil\log(L)\rceil$. For a reasonably large image this will not be a large difference.

For example, a 1024 by 1024 image will require 10 jump iterations. A 512 by 512 image will require 9 jump iterations. Anything between the two sizes will also require 10 iterations. Even in the worst case of an image only just over a power of 2, like a 513 by 513 image, it will only require 1 additional iteration, which at this scale is approximately 11% more (10 instead of 9).

Non-square images are less efficient per area

Since the number of iterations required is determined by the longest side length, the time taken for a 1024 by 1024 image will be the same as for a 1024 by 16 image. A square image allows a larger area to be covered in the same number of iterations.

Separating axes is likely to reduce quality

Section 5 of the paper describes possible errors. Every pixel is reachable by a path from every other pixel, but some paths do not bring the correct nearest seed, due to that seed not being the nearest to a previous pixel in the path. A pixel that does not allow a seed to propagate past it is said to have "killed" that seed. If the nearest seed to a pixel is killed on all paths to that pixel, then the pixel will record some other seed and there will be an incorrect colour in the final image.

Only one path needs to exist that does not kill the seed in order for the final result to be correct. Incorrect colours only occur when all paths from the correct seed to a given pixel are blocked.

If a path involves alternating horizontal and vertical jumps, separating axes will make this path impossible (all horizontal jumps will be taken before all vertical jumps, making alternating impossible). Diagonal jumps will not be possible at all. So any path that alternates or contains diagonal jumps will be excluded. Every pixel will still have a path to every other pixel, but since there are now fewer paths there is more chance of a given pixel being blocked from receiving the correct seed, so the final result will be more error prone.

Separating axes is likely to make the algorithm take longer

The efficiency would probably be reduced by separating axes, as the flooding would no longer be done in parallel, but would instead be repeated for each axis. For 2D this would likely take approximately twice as long, and for 3D approximately 3 times as long.

This may be somewhat mitigated by the lack of diagonal jumps, but I would still expect an overall reduction in efficiency.


I would think rendering to a downsampled FBO would be much faster than generating mipmaps.  glGenerateMipmap will generate the complete mipmap pyramid, which seems wasteful:

https://www.khronos.org/opengles/sdk/docs/man/xhtml/glGenerateMipmap.xml

Besides, only OpenGL ES 3.0 lets you manually select mip level when texturing:

https://www.khronos.org/opengles/sdk/docs/man3/html/textureLod.xhtml

Artistically, this might not give you the look you want.  In my experience, simply downsampling gives a grainy look.  I would recommend also doing a pair of one-dimensional blurs, as described in this article from GPU Gems:

http://developer.download.nvidia.com/books/HTML/gpugems/gpugems_ch21.html


I've figured a difference approach to projecting the point cloud inside Houdini.

I had assumed that the 'Scene Depth World Units' would provide the depth of each pixel along a vector parallel with the camera vector, rather than angled to point at the camera as a single point.

So rather than projecting the point cloud as if the depth was correct, and then trying to adjust for the distortion this creates, it makes more sense to just project the points in the same way they were created.

So, knowing the position and rotation of the camera, I can then use the depth pass to project the points outwards from this location.


Why is work-efficiency desired in GPU programming?

First of all, re: "GPUs are simply giant SIMD processors and should operate in lock-step", it's a bit more complicated than that. The entire GPU does not run in lockstep. Shader threads are organized into groups of 32 called "warps" (on NVIDIA; on AMD they're groups of 64 called "wavefronts", but same concept).  Within a warp, all the threads do run in lockstep as a SIMD array. However, different warps are not in lockstep with each other. In addition, some warps may be actively running while others may be suspended, much like CPU threads.  Warps can be suspended either because they're waiting for something (such as memory transactions to return or barriers to clear), or because there isn't a slot available for them (since the GPU can only actively run a certain number of warps at a time).

Now, back to your question. I can see two ways that the "work-efficient" algorithm from that paper looks like it would be more efficient than the "naive" algorithm.

The work-efficient version requires half as many threads to begin with. In the naive algorithm, they have one thread per array element; but in the work-efficient version, each thread operates on two adjacent elements of the array and so they need only half as many threads as array elements. Fewer threads means fewer warps, and so a larger fraction of the warps can be actively running.
Although the work-efficient version requires more steps, this is offset by the fact that the number of active threads decreases faster, and the total number of active threads over all the iterations is considerably smaller. If a warp has no active threads during an iteration, that warp will just skip to the following barrier and get suspended, allowing other warps to run. So, having fewer active warps can often pay off in execution time. (Implicit in this is that GPU code needs to be designed in such a way that active threads are packed together into as few warps as possible—you don't want them to be sparsely scattered, as even one active thread will force the whole warp to stay active.)

Consider the number of active threads in the naive algorithm.  Looking at Figure 2 in the article, you can see that all the threads are active except for the first 2k on the kth iteration.  So with N threads, the number of active threads goes like N − 2k. For example, with N = 1024, the number of active threads per iteration is:

1023, 1022, 1020, 1016, 1008, 992, 960, 896, 768, 512

If I convert this to number of active warps (by dividing by 32 and rounding up), I get:

32, 32, 32, 32, 32, 31, 30, 28, 24, 16

for a sum of 289.  On the other hand, the work-efficient algorithm starts with half as many threads, then it halves the number of active ones on each iteration until it gets down to 1, then starts doubling until it gets back up to half the array size again:

 512, 256, 128, 64, 32, 16, 8, 4, 2, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512

Converting this to active warps:

16, 8, 4, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 8, 16

The sum is 71, which is only a quarter as many. So you can see that over the course of the entire operation, the number of active warps is much smaller with the work-efficient algorithm. (In fact, for a lengthy run in the middle there are only a handful of active warps, which means most of the chip is not occupied. If there are additional compute tasks running, e.g. from other CUDA streams, they could expand to fill that unoccupied space.)

All that being said, it's unfortunate that the GPU Gems article does not clearly explain any of this, instead focusing on big-O "number of additions" analysis that, while not entirely irrelevant, misses a lot of the details about why this algorithm is faster.


What factors affect which shader language to learn?

The shader language is bound to the APIs/engines that support it (glsl to openGL & WebGL and hlsl to D3D). There are tools to translate from one to the other but they aren't perfect.

This means the main reason to pick one over the other is which platform you will be working on.

However by the end of the year a new binary shader representation will come out called Spir-V designed for Vulkan and an OpenGL extension for it is expected to follow shortly. This will allow you to write shaders in many languages as long as a translation tool to Spir-V exists (several are in the works including a Python and Haskell tool).


How to invert an affine matrix with small values?

I found a solution to my specific problem. Instead of computing the determinant and hitting the precision wall, I use the Gauss-Jordan method step by step.

In my specific case of affine transformation matrices and the range of values I use, I don't hit any precision problem this way.


What is the equivalent of rotation matrix in Euler angle space for Quaternions?

It's not true to say quaternions don't have a rotation matrix. They are both ways of completely describing a rotation in 3D, so they are freely convertible. Each quaternion corresponds to exactly one rotation matrix, and each rotation matrix corresponds to one quaternion and its negative. Quaternions are just more compact and easier to interpolate.

To work out the rotation matrix for a quaternion, you need to take the effect of the quaternion, $q = (w, x, y, z) = (w, \textbf{v}) $ on some point $\textbf{p}$, that is, $(w, \textbf{v}) * (0, \textbf{p}) * (w, \textbf{-v}) $, and express this formula as a matrix operating on $\textbf{p}$.

Multiplying out using the rules of quaternions gives us $(w^2 - \textbf{v.v})\textbf{p} + 2( w\textbf{v}\times \textbf{p} + \textbf{(v.p)v})  $

This can be rearranged into 3 matrix terms - the first part is just a constant value down the diagonal, the second part with the cross product has values everywhere except the diagonal, and the third one is simply $\textbf{v}_m *  \textbf{v}_n $ for rows m and columns n.

Adding these 3 matrices together gives you

$\begin{bmatrix} 1-2(y^2+z^2) &2(xy-zw) & 2(xz+yw)\\2(xy + zw)&1-2(x^2+z^2) & 2(yz-xw)\\2(xz-yw) & 2(yz + xw) & 1 - 2(x^2+y^2)\end{bmatrix}$

which is the conversion you need.

In your particular use case, though, you are just doing a simple rotate-x through 90 degrees, which simplifies things a lot. Remember a quaternion is $ ( cos(\theta/ 2) , sin(\theta/ 2)\textbf{a})$   where $\textbf{a}$  is your normalized axis of rotation. So the quaternion representing ±90 degrees rotation around x will be $( \frac{1}{\sqrt{2}},±\frac{1}{\sqrt{2}},0,0)$