What's a parent vector space or coordinate system?

You have to step outside the world of mathematics for a while and look at what we are trying to achieve. Mathematics in its purest form only tells us what kinds of properties certain constructs have, it does not tell us why and how that would be useful.

The modeler's perspective on things

So we have to look at the problem from modeling perspective. In this perspective it is mostly irrelevant how the mathematics work. As long as it provides the modeler with a system that yields with a understandable toolset to do his work everything is as it should be.

It is frequently beneficial to describe things such as physical objects (creatures, mechanisms etc.), forces or any number of other things independent of each other. Either because it's easier to think in one orientation or just simply because you want to move/rotate/scale objects out of convenience.

In practice you want at minimum that the camera is separated from the rest of the objects. While in more complex cases you want to have a more rich system. For example a bug is walking straight ahead on a straw. For practical reasons it would be neat if the bug only move forward in relation to the straw while the straw sways in the wind in relation to the world and the camera is independent of the straw. Similar modeling cases happen a lot, in human hands, robotics, car simulations, Finite object modeling etc.

Computational perspective

The rendering and other functions need to be able to transform elements from these modeling hierarchies to each other. For example you need to know the point location in camera space to render. Preferably all this space transform happens automatically. For this purpose the modeler tells your application how the modeling elements are structured in relation to each other, a hierarchy, for modeler's benefit.

Although this could be done in any number of ways it is commonly done with a tree structure, much like a filesystem with much the same effect. Each level stores in addition to tree information a affine matrix that describes the transform to its predecessor in the hierarchy. Mainly because this takes care of two core issues. The tree structure is easy to traverse and find the path to any other part. The matrices solve the problem of needing to do a lot of work transforming form a transform from one hierarchy space to another.

So the user provides a tree where each object describes the transform required to transform the object to the coordinates of the object it is described in relation to, the parent. This parent in turn in turn describes the transform to its parent all the way to the world. Since all objects share one common ancestor, world, we can easily traverse the tree to find a transform path from any one coordinate system to another coordinate system.

Since linear transforms are invertible, we do not need to store a description separately for moving in the opposite direction, away from world towards the children/leaves of the tree structure. All we need to do is to invert the matrix. This is almost always needed as it's practical to describe camera as a separate object under the world, there is no reason to describe the cameras transform from world to camera just treat it like any other object.

To compute a point from one space to another you do as follows:

Use the tree to find a relationship between spaces.
For each step on the way pre multiply your matrix/ or its inverse if you're traversing away from world. **

Example:  You have the following hierarchy because it makes sense (this way you do not need to model each petal separately)

World
  |
  +---- Flower
  |      |
  |      +--- Petal 1
  |      |
  |      +--- Petal 2
  |      |
  |      +--- Petal 3
  |
  +---- Track
         |
         +--- Boom
               |
               +--- Camera

To render Petal 1 you need to transform the vertices of petal one into the camera object space and from there to perspective space etc. Your calculation would be:

$$
T_{\text{Perspective}} \cdot T_{\text{Camera}}^{-1} \cdot T_{\text{Boom}}^{-1} \cdot T_{\text{Track}}^{-1} \cdot T_{\text{Flower}} \cdot T_{\text{Petal 1}} \cdot \vec{v}_{\text{PetalPoint}}
$$

This way you can model a flower in a flowerpot totally independently from  from the robot shaking the camera about. At the same time freeing the modeler of the need to do complex trigonometric calculations, and in fact a big load of mathematics. Pretty neat.

Now that you understand why, you should understand also that we do not have any limit to complexity of the hierarchy to
facilitate a broad range of things you could do.

This kind of hierarchy is called a scene graph and such things can be found in web pages, vector and CAD applications, many scientific modeling applications, most games.

This should answer why as well as the less exciting how.

Do we absolutely need these?

No, we do not. We can dispense all the matrices and describe all of our data directly in the camera's final coordinates. It's just that this is horribly complicated in 98% of the cases.

In practice we need this hierarchy if we try to do even something more complex than a stationary box. Most of this deign is to make scene modelers' life easier. So that simple operations make more sense.

* why do we call it a parent child relationship? Because it's a tree, and earliest tree graphs are genealogical studies, aka family trees. Where the previous level is simply the parents and the next one is their children.

** or the other way around if you want to transpose your calculations.


Rendering same vertex-array in different modes

Your shaders are fine, the problem is topological, that is the way the vertices are connected to create triangles.

For cases 1 & 2: You cannot really draw a grid of points in this mode (GL_TRIANGLE_FAN) as vertex 0 (the top left one) will be used in every triangle, as well as other issues. You should be using Triangle/Line strips, although these aren't so straightforward as some vertices are shared between each row.

For example, imagine you have vertices arranged like so:

  0 1 2

  3 4 5

  6 7 8

You would draw strips by sending the above vertices in the following order:

  0,3,1 (outputs triangle 0,1,3)

  4,    (outputs triangle 1,3,4)

  2,    (outputs triangle 1,2,4)

etc...

This you can do in code by scanning across 1 row of vertices and including the vertex on the next row for every other vertex.

e.g.

int numColumns = 5;   // As per your example..
int numRows = 5;
int vertexOffset = 0;
int vertexSize = 3;          // 3 floats per vertex
int nextRow = numColumns * vertexSize;

// This loop is for each row-strip, which requires 2 vertex rows..
for(int j = 0; j < numRows - 1; j++) {

    // Inner loop sends vertices for each row-strip
    for(int n = 0; n < numColumns; n++, vertexOffset += vertexSize) {
       glVertex3fv(vertices + vertexOffset);            // Current row
       glVertex3fv(vertices + vertexOffset + nextRow);  // Row below
    }
    // This is special, we send a duplicate vertex to restart a new strip for the next row (This can be omitted if j == numRows - 1)
    glVertex3fv(vertices + vertexOffset - vertexSize + nextRow);
}

If you are concerned about sending duplicate vertices you could use an additional index buffer, as that will allow you to use the exact same vertex buffer you used for the GL_POINTS version.


Triangulation of vertices of an ellipsoid

If you sample the two parameters $\eta$ and $\omega$ with steps $d\eta$ and $d\omega$, then you'll get a grid of points $v_{ij} = f(i\;d\eta,j\;d\omega)$.  Any four adjacent points will define a quadrilateral.  To get triangles, you just have to split each quad in two by a diagonal.

So in the example, you'd split the quadrilateral $\{v_{00},v_{01},v_{11},v_{10}\}$ into two triangles $\{v_{00},v_{01},v_{10}\}$ and $\{v_{01},v_{11},v_{10}\}$.  You'll want to make sure the points come in the right order so the triangle orientations match (or you'll end up backface-culling half of your triangles.)

Your second equation, for $n_{xyz}$, is finding the normal of each vertex, which you can use for shading.


How to capture rendered frames from an application and stream over network

I think you're confused on how rendering and animation work.

A traditional model is nothing more than a bunch of triangles. So a model file coming from say Blender, etc. is just a list of vertices for the triangles. (with some added stuff if you want)

Rendering takes the scene definition and transforms a 3D scene (aka the models oriented somewhere in space) and creates a 2D picture.

Animation creates the illusion of movement by showing static pictures very fast (24 frames per second or more).

Therefore, to have a spinning mug, your application would render each frame, then present them, at, say 30 fps.

Now, to stream the data, you can take two paths:

Render each frame on the server and stream them to the client with something like H.264
Stream the graphics commands to the client, and let them render them.

The first option is what something like OBS or Fraps and (maybe) Nvidia GameStream is doing.
There's a paper that covers the 2nd option here


What is the difference between world coordinate, viewing coordinate and device coordinate in computer graphics?

The 2D pipeline involves ... [coordinate transformation terms]
Can someone give me the detail differentiation among these?

This is something I very recently learned while trying to understand how software handles 3D graphics behind the scenes. I've never encountered these terms in reference to 2D graphics before, but the ideas still apply. 2D graphics are just 3D graphics without a few fancy tricks.
Very Generally
Most of these different coordinate spaces don't technically exist. They come about by considering the objects you wish to draw from some reference point. When you change that reference point, you change the coordinate space. Math gets involved, but I won't do any here.
The Drawing Device and Device Coordinates

... converting normalized coordinate to device coordinate.
Isn't viewing coordinate and device coordinate the same things?

Short answer, no, but they can be in certain 2D cases. In 3D, pretty much never.
There seems to be a confusion not just between view and device coordinates, but also with the distinction between normalized and device coordinates. Lets start with normalized and device, and we'll get to view later.
In 3D systems (and 2D systems setup to work like them) you are presented with a "drawing device" consisting of two parts: drawing tools (usually functions accessed from application code) and a drawing area. Ignoring the tools, the important thing is that the drawing area has a coordinate system independent of pixels which adds a lot of flexibility and generality to the device.
Any time something is setup to work numerically in a range from 0 to 1 it is "normal", and often the range 0 to -1 is included too. With a drawing device the coordinates for both the X and Y directions of your drawing area range from -1 to 1, with -1 meaning "all the way left" for X and "all the way down" for Y, and 1 meaning "all the way right" or "all the way up" for X and Y respectively.
You can easily add a Z direction too (also ranging from -1 to 1) for "3D" drawing, but all this does is give you a layered 2D image. There is no perspective (or "depth") in your final drawing which is great for 2D games and graphics. The reason for this is that the drawing device will simply divide the space up into pixel size steps along X and Y then draw each pixel back to front (-Z to +Z). If you want real 3D though, you'll need some extra mathematical steps.
It is possible (and common) for the device to understand and work with coordinates outside the -1 to 1 range. However, if we want an object to appear on screen than it must be "normalized" so that any coordinates it refers to are within -1 to 1. After this is done to an object the terms "normalized coordinates" and "device coordinates" are practically synonymous, just one refers to the object you altered so that it would fit within the drawing area and the other refers to the drawing area itself.
Object and World Spaces

The 2D [or 3D] pipeline involves with the construction of world coordinate scene ...

Yes, but you start in "object space".
When someone creates a model of an object (or draws a sprite) all that matters to them, coordinate-wise, is where zero is within the space they're working. This consideration is called "object space" because the artist isn't concerned with the -1 to 1 drawing area, only accurately representing the object they'd like to draw. A model of a pen might be 6 units long and a model of a person might be 5.5 units tall. The important thing is that all of the points that make up some object are positioned with regard to the coordinate (0, 0, 0) and nothing else.
When later loaded into memory they will be too big for the drawing area (which is only 2 units tall and wide) and the pen will appear larger than the person. It's up to us humans to know what the units mean and adjust them properly. Resizing and positioning all objects in relation to each other instead of each individually being sized with relation to zero transforms the objects into "world space".
From the shortness of your question I wasn't sure if you knew that yet or not, but I thought I'd include that for completeness, as a good foundation for answering the next part of your question.
View Space

... followed by converting world coordinate to viewing coordinate ...
Isn't viewing coordinate and device coordinate the same things?

Note that in my explanation above nothing ever happened to the drawing device, or the coordinates of the drawing area. The only thing that ever changed was the objects we wanted to draw. The drawing area is fixed and unchanging, only ever drawing any object that fits within the coordinates -1 to 1.
Each object loaded appears at the center of the drawing area at whatever size it was designed to have (object space) which is practically always too big for the drawing device. After rearranging and resizing every object around each other (object space → world space) we still have a final scene that is too big for the drawing area, and we have no way to move around within the scene. This is where we transform everything into "view space" and work with "eye coordinates". Again, there is no such space, but what we do to work within it is to consider an imaginary camera and, treating the entire world as one giant object, we resize and reposition everything with respect to the position of the camera. That is "view space", the movement of everything around a camera, and the terms "view coordinates" (or "eye coordinates") refer to the positions of everything in relation to the camera.
View to Viewport
At this point everything that you want to see should be within the device space, but there are still a few last tasks to perform to get everything to draw.
If you were working on a 3D scene instead of a 2D one, you would perform a "perspective transformation" to create the illusion of depth. This is where you squish and warp things to make them look like they're shrinking into the distance.
Anything outside the drawing area needs to be thrown out. This saves on work for the computer, and it also is considered another transformation from view space to "clip space" because you are "clipping away" the extra data.
And finally, after clipping, you'll do a "viewport" transformation that turns the -1 to 1 coordinates into pixel coordinates.
Summary
Well... I didn't intend for this answer to get so long. My apologies if I covered anything you already knew. It's hard to guess what someone does or does not already understand. Please leave a comment if you have any further questions.
Object Coordinates: The positions of points relative to a zero point. The points referred to are the building blocks of some single object.
World Coordinates: The positions of collections of points (objects) relative to a single shared standard zero point.
View Coordinates: The position of all points (the whole scene) relative to a camera.
Device Coordinates: Any position within -1 to 1 on both the X and Y axis (and Z too for 3D). Any point within the drawing area of the device.
Normalized Coordinates: The positions of points of any object (or the whole scene) that have been transformed to fit within the -1 to 1 space of the drawing device.


Those are moire patterns. They are  an aliasing artifact that usually occurs when sampling on a regular grid. Did you jitter the positions of your samples? If you just sampled an evenly spaced 10x10 grid within each pixel, that could explain it. Also, numerical errors or inaccuracy could it.


Method for rendering heightmap from a 3D scene

Using OpenGL, or any such library, you could transform your scene with an orthogonal matrix and move your camera view point to be looking above (and parallel to) the scene (in a vertex shader) and then use an empty fragment shader which is going to write to the depth buffer of the current attached framebuffer. Finally retrieving this depth buffer gives you the heightmap of the current view.

If your scene is too big, you could move the camera higher and also move it around taking "chunks" of the map at a time and assembling your heightmap manually by "stitching" these "chunk" heightmaps together.

Hope this helps.


Jos Stam's "Stable Fluids" boundary conditions

It's Dirichlet boundary condition. The quantities at domain boundaries will take either the same value as its inner neighbor or the negated value, depending on the condition $\mathbf{b}$. Take velocity field for instance, velocity of fluid will either gets reflected or not change at domain boundary. The condition $\mathbf{b}$ is a user option, $\mathbf{b}==1$ is used to set reflected values at the horizonal boundaries, and $\mathbf{b}==2$ is to set reflected values at vertical boundaries. To be more clear:

$\mathbf{b}==1$: field values at horizonal boundaries get reflected, field values at vertical boundaries do not change.
$\mathbf{b}==2$: field values at vertical boundaries get reflected, field values at horizonal boundaries do not change.
$\mathbf{b}== other\ value$: field values at neither vertical nor horizonal boundaries do not change.

set_bnd could be called multiple times to set desired boundary values.


It stands for texture coordinate.

It's a 2D value between 0 and 1 to indicate where in the texture to sample from for the color.

For sphere that means creating a map projection to map the curved surface to the flat plane, a simple option is to use the lat/long remapped to 0-1 though that gives awful distortion at the poles.


You're seeing channels and samplers there, not individual keyframes.  In this case, one animation is simultaneously targeting the translation, rotation, and scale channels of the same node, with three separate "samplers" (sets of keyframes per channel).

Not shown above, your glTF file has a section called accessors that tells viewers how to slice out portions of binary data from an associated .bin file.  The .bin file contains mostly mesh position/normal/UV data, but also contains the keyframe times and values.

So in the above snippet, that first channel targets the translation of node 3 with sampler 0 (the first sampler).  The sampler there asks for LINEAR interpolation of keyframes from animation times (input: 8) for translation keyframe values (output: 9).  In the accessors section, accessors 8 and 9 (using zero-based indexing) are references to the keyframe times (in 8) and values (in 9).  They are separate because the data type of the output depends on what's being targeted (rotations are vec4 quaternions for example, while translations and scales are vec3, and the input times are just floats).

The accessors in question (8 and 9 for example) will contain a value "count" that effectively tells you the number of keyframes.  If you have 5 keyframes for example, accessor 8 should have "count": 5 keyframe times, and accessor 9 will also have "count": 5 keyframe values.  The glTF Validator will flag an error if these counts don't match.


Think about it this way. To calculate YCbCr of (1, 0, 0), you'll end up with Y = 0.299, Cb = -0.16873, Cr = 0.5. If you take the magnitude of the chroma channels (sqrt(Cb*Cb + Cr*Cr)), you get ~0.53 for the "saturation". The YCbCr color space has a different shape than RGB, so it may represent values differently enough that you don't get the exact same answer in both cases. Furthermore, our perception of luminance and saturation is not consistent among colors.

If all you want is to calculate the saturation of an RGB component, it can be calculated by:

float saturation = max(r, g, b) - min(r, g, b);

Note that if you attempt to calculate perceptually uniform amounts of color saturation you end up with systems like the Munsell Color System where the shape of the color space is highly irregular compared with the cubes and cylinders of other color spaces.

If you tell us what you're trying to achieve, we might be able to offer an actual solution to whatever the problem is.


Unwinding an image on a spiral to make it long and flat

A quickly formulated method, read first one that popped in my brain (not best), could be. Find the closest points on a parametric spiral for each sample (read A Pixel Is Not A Little Square3). Then place the samples on a line by placing the pints in one axis by how far they are from your spiral line and the other by what the closest point is.

You can then use something like poisson disk sampling to reconstruct the image (maybe the didsk should be longer ovals the further away they are form the center line and less long as you progress further into the spiral).

Now that I think of the problem a bit more. A simpler approach could be to map the image on mesh, using UV mapping, then unwrap the mesh. A triangular mesh might exhibit some artifacts so using a B-spline interpolation could be better. Several image mesh algorithms exist and you even have these in Photoshop.

Image 2: A proof of concept of unwound spiral with the mesh method. Note i scaled the unwound spirals length to make it fit the site requirements better.

From initial test in my 3D app of choice seems that a dense enough triangular mesh works fine.


Methods for grid traversal in a glsl pixel shader?

why not building a bounding box (or spheres) hierarchy ?

(but for a shadertoy implementation, the lack of dynamic loop length might spoil the gain ).


So you have a series of points and, at each point, a supplied derivative?
Is a piecewise cubic sufficient or does it need higher derivative continuity?  If the former is ok, then Cubic Hermite Splines will do the job. (If you need, they can be trivially mapped into cubic Beziers)


How's this laser cutting effect with backlight achieved?

This question may be off-topic here, but this is almost certainly a 2D effect, and it's achievable without much difficulty in After Effects or any similar package once you've broken the effect down into individual elements.
This isn't exactly the same effect, but it gives most of the ideas: https://www.youtube.com/watch?v=uHFmWu__Q4k
The key point here is that you need to look not at the effect as a whole, but break it down into the individual elements that make it up. The parts make the whole.


Computing a rotation: complex numbers vs rotation matrix

Both methods end up doing the same calculations when you break it down.

Rotating a vector $u$ with a matrix:
$$\begin{bmatrix}\cos\theta & -\sin\theta \\ \sin\theta & \cos\theta\end{bmatrix} \begin{bmatrix}u_x\\u_y\end{bmatrix} =
\begin{bmatrix}u_x \cos\theta - u_y\sin\theta \\ u_x \sin\theta + u_y \cos\theta \end{bmatrix}$$

Rotating a vector $u$ using complex numbers:
$$\begin{aligned}
(\cos\theta + i\sin\theta)(u_x + iu_y) &= u_x\cos\theta + iu_y\cos\theta + iu_x\sin\theta - u_y\sin\theta \\
&=(u_x\cos\theta - u_y\sin\theta) + i(u_x\sin\theta + u_y\cos\theta)
\end{aligned}$$

I wouldn't expect either one to be appreciably faster or slower than the other, since they all end up doing the same set of basic operations, i.e. 4 multiplies and 2 adds.

Conceivably, complex numbers could be faster when you have a large number of rotations stored in an array, because they have only two components instead of four and therefore more of them fit into each cache line.

  On a related note, is there some spatial transformation that complex numbers can do but matrices cannot?

No. Matrices are more general than complex numbers. Any complex number $z$ can be represented by a matrix as:
$$\begin{bmatrix}\text{Re}(z) & -\text{Im}(z) \\ \text{Im}(z) & \text{Re}(z) \end{bmatrix}$$
This corresponds to rotation by the phase of $z$ combined with scaling by the magnitude of $z$. Complex numbers can only represent rotation and uniform scaling. Matrices can represent those, but also nonuniform scaling and shearing. Another way to see it is that complex numbers have only two degrees of freedom, while 2×2 matrices have four degrees of freedom.


The rasterizer determines, which pixel is rendered depending on vertex positions that come out of the vertex shader. For each rasterized pixel, the rasterizer knows where it lies inside the triangle formed by its three positions. With the known position inside the triangle, the rasterizer can also interpolate the vertex values for each pixel. So it also calculates the interpolated color, in case your vertex shader passes a color down the pipeline. The interpolated value is then provided to the fragment shader. However, the fragment shader can still modify the color before it is finally rendered to the screen.

  Does that mean rasterizer has access to color attribute

It has full access to everything that comes out of the vertex shader. A vertex shader must always have a position as output (gl_Position). In some special cases also other outputs are required. However, if your initial vertices had a color, but the color is not passed down the pipeline by the vertex shader, the rasterizer won't produce an interpolated color. It only interpolates each variable/attribute, that comes out of the vertex shader.

  how is the color attribute passed to the rasterizer, as the inputs for rasterizer are vertices?

You can define arbitrary outputs in your vertex shader. Even if the initially submitted vertex, did not have a color attribute, you can add a color output in the vertex shader. The data the rasterizer gets per vertex is just the vertex shader outputs. A vertex is a data point containing a position and optional additional values/attributes. So if a vertex has a color attribute, the rasterizer can see it, since it is part of the vertex.

Note that the vertices you pass into the vertex shader and the vertices that come out of it are not the same and might differ.

  Or is it happening in the fragment shader?

No, the fragment shader receives already interpolated values as input. But it can decide to pass it directly, ignore it or modify it as you wish. If you just pass it through, you will see the interpolated vertex colors.

Further information can be found on this site


What you're ultimately trying to do is, between draw calls, change some state that allows one to select which set of read-only data to use in the rendering process. There are many techniques for this.

Push constants can be used to provide an index into an array of arbitrary size stored in an SSBO. Indeed, a single push constant can be used to provide multiple indices into multiple SSBOs, for different frequencies of state (per-frame vs. per-object). The upper 16-bits can be for one frequency, with the lower 16-bits for another.

That may sound slower than other methods, since you're accessing global data rather than an actual UBO. But on some hardware, it won't change anything. AMD hardware for example doesn't actually have any specialized hardware for UBOs, so they implement UBOs as just read-only SSBOs. Even on hardware with dedicated UBOs, the state data you're reading will likely be quickly cached, since all invocations will be reading the same values from the array. So unless your thrashing the SSBO caches, this shouldn't impact performance much.

Dynamic UBO descriptors also exist, given the VK_DESCRIPTOR_TYPE_UNIFORM_BUFFER_DYNAMIC type. Basically, the idea is that, when you call vkCmdBindDescriptorSets, you provide a dynamic offset for any dynamic descriptors. The offset represents a byte offset from the buffer attached to that descriptor to where the UBO data you want to fetch lives. So you would have all of your data in a large buffer, and you would select which particular data you're using for a particular render call by invoking vkCmdBindDescriptorSets with a different dynamic offset.

Now vkCmdBindDescriptorSets may seem like a more heavyweight call than setting push constant. But the implementation has all the information needed to be able to see that you are using the same descriptor as you used in the previous call in that command buffer, and thus not do the actual heavy part of that operation, simply applying the dynamic offsets to any dynamic descriptors. So that part should not be a concern.

Push descriptors are mainly intended for texture/sampler/image bindings, which don't have a convenient mechanism like dynamic offsets that would allow you to quickly switch which image/sampler you're using. That being said, push descriptors are a useful concept even for UBO/SSBOs. Whether they are faster than dynamic descriptors is a matter of profiling.

The principle downside is the obvious: push descriptors is an extension. Dynamic UBO offsets and push constants are mandatory, required by every Vulkan 1.0 implementation. Push descriptors aren't even optional features of Vulkan. They're important enough to get a KHR extension, but they didn't go into Vulkan 1.1 as an optional feature. This is not a theoretical issue: NVIDIA and Intel seem to support it on desktop, but no official AMD drivers do.

So I would go with dynamic UBOs; this is exactly what they were made for, after all.


Multiple buffers and calling glBufferSubData

It depends. There are a few competing factors at work here.

First, consolidating terrain data into fewer buffers (or one) may allow you to combine multiple terrain patches together in a single draw call—assuming there aren't other state changes between patches that would prevent this. So, you may be able to reduce draw-call overhead this way, or at least reduce the state-change overhead between draws. That being said, 50–100 draw calls is not that many, so this may not be a real performance concern.

In addition, do you anticipate needing to use the terrain data for other shader effects? For example, maybe a shader procedurally generates grass growing on the terrain, or you use height above the terrain to influence lighting or AO or some such. In that case having all the terrain data in one big buffer that you can pass to those shaders may simplify things.

On the other hand, if you're having the terrain deformed/edited in real-time (which it sounds like you are), then it may be faster to update the entirety of a small buffer than to update a small portion of a large buffer. The problem is when the CPU and GPU are working asynchronously, and you update a buffer from the CPU, the GPU may still be using the previous version of that buffer. To fix this, the driver internally keeps 2–3 instances of the buffer (much like double/triple buffering the screen). If you always update the entire buffer at once, the driver can just put the updated data into the "next" instance of that buffer. But if you update only part of the buffer, the driver must copy all the rest of the data from the current instance to the next instance. This extra copy will consume some performance either on the CPU or on the GPU.

So, if the terrain data was entirely static I'd say it's better to put it all in one big buffer. However, if you're going to be updating portions of the data in real-time, it may be better to keep it in smaller buffers so that you can update exactly the ones that have changed while leaving others alone.

In the end, you'll have to do some profiling to see which approach is better.


Which 3D applications currently uses Vulkan API?

I'm a newbie in terms of CG world. I'm curious to know about the applications that use Vulkan API and their performances compared to other APIs

There are some popular magazine websites, such as PCGamer and AnandTech that has some comparison charts between other rendering APIs.

PCGamer compares OpenGL vs Vulkan for Doom (2016)

AnandTech compares Vulkan vs DirectX12 API Overhead

I would argue that comparisons make sense on the application-level - one API wouldn't necessarily be better on all occasions than the others. It totally depends on how the graphics developers of the software have architected their system to unlock the APIs full potential.

  As I read, Vulkan has more access to the graphic cards which helps in performance. If so, do the current rendering engines use it?
  <
  Solidworks, Revit, Sketchup, Unity3D, Unreal Engine, After Effects, etc. ?

Wikipedia: Software that supports Vulkan page has a list of software that supports Vulkan API. As for your question, quoting from Wiki page:

  Unreal Engine 4 – In February 2016, Epic Games announced Unreal Engine 4 support for Vulkan at Samsung's Galaxy S7 Unpacked
  event.[64][65]
  Unity – The engine has support for Vulkan since version 5.6.

There's plenty of others which are also popular, video games and game engines listed in the Wikipedia page, such as:

Dota 2
Doom
vkQuake
Star Citizen
Wolfenstein II
CryEngine
Source 2


Why not calibrate all display devices to sRGB?

Simple answer: because it costs. :) Some higher-end monitor models are indeed calibrated at the factory (and come with a calibration report sheet) but it would not be cost-effective to do so for all models in the product line.

Color reproduction varies even between batches of the exact same panel model, so it isn't possible to just use a set calibration for a given monitor model, it need to be done on a unit-by-unit basis. So to save costs, they instead determine the color reproduction curves of the "average" panel and use that for an ok but inconsistent result.


Actually what it meant by Affine Covariant regions is that covariant regions in two images which are related by some affine transformation. So the regions found in one image and those are exactly same regions in other image which have been transformed through affine transformation.


GLSL has built-in fragment shader inputs, and one of them is gl_FragCoord. You can get the $x$ and $y$ value from this and use that as the lookup coordinate of the texture map in a sequence of two passes, but that could be a good or bad thing depending on the application (such as caching the result). This would both be rotation invariant while allowing you to scale and move the mesh.

First Pass

Transform the mesh without rotation. Bind a texture to the fragment shader and write the interpolated UV coordinates to it.

Second Pass

Transform the mesh correctly this time. Bind the texture from the last pass and read from it this time. You preserved the old UV coordinates, and now you can use these to look up whatever texture you have before. Obviously, this works best with a mesh that's a sphere (such as an eye).

vec2 texCoords = gl_FragCoord.xy; // do some sort of scaling
vec2 oldTexCoords = texture(samplerOfUVTexture, texCoords); // sample UV

Now you have both the old UV coordinates and the new UV coordinates. You can use the new UV coordinates for you diffuse map and the old UV coordinates for the occlusion map.


Converting from RGB Color Picker to Surface Reflectance

Yes, I think you're right. The values from a typical color picker would be in a gamma encoding or sRGB encoding and would need to be converted to linear to use as reflectance. This is probably just a mistake in the book.
(BTW, in case you're not aware, sRGB isn't actually gamma 2.2 but a more complicated transfer function.)


You can use a stencil instead, This allows the gpu to do an early out test before the fragment shader even starts.

You only have to build the stencil once (possibly with a discard in a simple shader) and you can reuse the stencil in future draws.


Tessellation Speed in terms of interpolated points per second

If you’re only tessellating the model once and never changing the density of the mesh then tessellation is probably a waste of time here—any measurements you make would include all the extra work the GPU has to do to generate new points / triangles, which is orthogonal to what it sounds like you’re trying to measure. Tessellating the mesh ahead of time (in your tool of choice) and doing the interpolation / evaluation in a vertex shader would probably give you a better benchmark against a CPU implementation working on the same mesh.


When depth testing is enabled, the testing occurs after the fragment shader executes. The procedure is relatively simple (and spelled out pretty well here):

Get depth of fragment
Test against depth buffer
Discard if it fails

Early-z discard is when a rasterizer discards a fragment before a fragment shader runs because it's clear early on that it won't pass the depth test. However, the rasterizer can't safely do this if certain things happen (such as calling discard() or modifying the fragment depth in a fragment shader). This is a hardware optimization and isn't required to occur unless you force it with OpenGL 4.2 or with an extension.

There are things that can prevent early-z discard from occurring. To achieve it:

Don't use discard() in a fragment shader
Don't modify the depth of a fragment in a fragment shader
Don't use alpha testing (deprecated)
Make sure you have depth testing turned on
Make sure you have a depth buffer


Generally edge detection boils down to detect areas of the image with high gradient value.

In our case we can crudely see the gradient as the derivative of the image function, therefore the magnitude of the gradient gives you an information on how much your image changes locally (in regards of neighbouring pixels/texels). Now, an edge is as you say an indication of discontinuity, so now that we defined the gradient is clear that this info is all we need. Once we find the gradient of an image, it's just a matter of applying a threshold to it to obtain a binary value edge/non-edge.

How do you find this gradient is really what you are asking and I am yet to answer :)

Lots of ways! Here a couple :)

Built in shader functions

Both hlsl and glsl offer derivative functions. In GLSL you have dFdx and dFdy that give you respectively gradient information in x and y direction. Typically these functions are evaluated in a block of 2x2 fragments.
Unless you are interested in a single direction, a good way to have a compact result that indicates how strong is the gradient in the region is fwidth that gives you nothing else but the sum of the absolute value of dFdy and dFdy. You are likely to be interested in an edge on the overall image rather than on a specific channel, so you might want to transform your image function to luma. With this in mind,when it comes to edge detection your shader could include something to the like of:

  float luminance = dot(yourFinalColour,vec3(0.2126, 0.7152, 0.0722));
  float gradient = fwidth(luminance );
  float isEdge = gradient > threshold;

With an high threshold you will find coarser edges and you might miss some, conversely, with a low threshold you might detect false edges. You have to experiment to find the threshold that better suits your needs.

The reason why these functions work is worth mentioning but I don't have time for it now, I am likely to update this answer later on :)

Screen space post-process

You could go fancier than this, now the field of Edge detection in image processing is immense. I could cite you tens of good ways to detect edge detection according to your needs, but let's keep it simple for now, if you are interested I can cite you more options!

So the idea would be similar to the one above, with the difference that you could look at a wider neighbourhood and use a set of weights on sorrounding samples if you want. Typically, you run a convolution over your image with a kernel that gives you as a result a good gradient info.  A very common choice is the Sobel kernel

Which respectively give you gradients in x and y directions:

You can get the single value out of the gradient as $ GradientMagnitude = \sqrt{ (Gradient_x) ^ 2 + (Gradient_y) ^ 2 } $

Then you can threshold as the same way I mentioned above.

This kernel as you can see give more weight to the central pixel, so effectively is computing the gradient + a bit of smoothing which traditionally helps (often the image is gaussian blurred to eliminate small edges).

The above works quite well, but if you don't like the smoothing you can use the Prewitt kernels:

(Note I am in a rush, will write proper formatted text instead of images soon! )

Really there are plenty more kernels and techniques to find edge detection in an image process-y way rather than real time graphics, so I have excluded more convoluted (pun not intended) methods as probably you'd be just fine with dFdx/y functions.


To answer your question we just need to write it as linear algebra equations and solve them. Although your question doesn't state it, I assume that $v$ and $d$ are unit vectors. Let's call the projected point $x$.

First, because the projected point is in the direction $d$, we can write:
$$\vec{vx} = \lVert\vec{vx}\rVert d $$

Second, because $p$ and $x$ are on the plane and $n$ is the plane normal, we can write:
$$\vec{vx}\cdot n = \vec{vp}\cdot n$$

We can now work on the equation until we get a definition for $\vec{vx}$:
$$\begin{align}
\vec{vx}\cdot n & = \vec{vp}\cdot n \\
\lVert\vec{vx}\rVert d\cdot n & = \vec{vp}\cdot n \\
\lVert\vec{vx}\rVert & = \frac{\vec{vp}\cdot n}{d\cdot n} \\
\vec{vx} & = \frac{\vec{vp}\cdot n}{d\cdot n}d \\
\end{align}$$

Or written differently:
$$x = v + \frac{(p-v)\cdot n}{d\cdot n}d$$

The result is undefined when the scalar product $d \cdot v$ is $0$, which happens when $d$ and $n$ are orthogonal, when $d=0$, or when $n=0$. The first case means $d$ is parallel to the plane and the projection doesn't have a point solution, and the two other cases mean that either the plane or the projection are not defined.


Figured it out, I did skip a step, you are supposed to multiply by the homogeneous coordinate.
Basically the actually projected point is (courtesy of GLM):
        vec<4, T, Q> tmp = vec<4, T, Q>(obj, static_cast<T>(1));
        tmp = model * tmp;
        tmp = proj * tmp;

        tmp /= tmp.w;

Kinda obvious in hindsight.


GLSL. Can someone explain why gl_FragCoord.xy / screenSize is performed and for what reason?

First, gl_FragCoord.xy are screen space coordinates of current pixel based on viewport size. So if viewport size is width=5, height=4 then each fragment contains:

Why are uvs needed? For example I rendered geometry to screen quad and then I need to apply some postprocessing on this quad in another rendering pass. To sample from that quad I need texture coordinates in range [0, 1]. So to calculate them I perform division gl_FragCoord.xy / viewPortSize.


Do we use 3x3 matrices in computer graphics?

That quote is a very strange way to phrase things. We definitely use 3x3 matrices in computer graphics. They tend to be most useful for doing affine transformations of 2D objects. It allows you to have scale, rotation, shearing, and translation (in 2D), but not perspective transformations. I believe that's what the quote is trying to say.

To get a perspective projection of a 3D scene onto a 2D plane (such as your computer screen), at some point you will need to divide the x and y components of each point in your geometry by some factor that foreshortens lines in the z dimension. That's the divide that they are referring to. See this article for more details. In it, they describe:

  The simplest perspective projection uses the origin as the center of projection, and z = 1 as the image plane. The functional form of this transformation is then x' = x / z; y' = y / z.


How to correctly implement Lambertian BRDF with point light

It doesn’t matter that it’s a point light; the 1/π normalization factor is what’s causing your surface to be 1/π as bright as you’d like it to be. You can either up your light intensity so the maximum brightness is π or just get rid of the normalization factor. There’s a good article with some notes about this and how it interacts with more complex lighting techniques here.


Given two arbitrary surfaces, what is the accurate and quite fast algo to find collision within X distance?

A key search phrase is "collision detection." You may find some references at this MO posting, "Intersection of Polyhedra."
And here is a more recent reference:

  Jiménez, Juan J., and Rafael J. Segura. "Collision detection between complex polyhedra." Computers & Graphics 32.4 (2008): 402-411. (Journal link.)

Google Scholar shows this has been cited by 33 papers since.


You can detect a matrix that can't be decomposed in TRS form by taking its 3×3 upper-left submatrix, interpreting its columns as vectors, and dotting them together in all combinations (1 with 2, 2 with 3, and 1 with 3). For a TRS matrix, the three dot products should all be zero, i.e. the column vectors should be orthogonal to each other. If any of the dot products are nonzero (beyond some epsilon to account for roundoff error) then the matrix is not decomposable as TRS.

This derives from the property that a rotation matrix is an orthogonal matrix, i.e. its column vectors form an orthonormal basis. When you multiply the rotation by a scaling matrix, the column vectors may no longer be unit vectors, but they should still be orthogonal. And the translation part doesn't affect the rotation-scaling part.


How can I combine two filter kernels for image processing to create a single kernel?

Why not just do them one after the other?

Img = (Img * median) * sharpen;

I think you should look into how these filter operations actually work before attempting to use them. It's pretty simple, use wikipedia or even any image processing textbook. Even the Matlab docs have good explanations.

Median filtering is a non linear filtering operation which means that there is no fixed kernel function for it. You may define a neighborhood size but other than that, nope, no kernel. Have you looked into how it works? I dont know how you could generate a "kernel" for a non-linear filtering operation.

Maybe a basic refresher in Image Processing is in order? Perhaps even going over basics of Signal Processing / DSP? I would suggest looking over 2D filtering if you can.


Blurring image with RGB values without convolving it with a kernel

Blurring an image depends on being able to read from the image. If you don’t have access to the pixel data, or external information about the image contents (e.g. average hue / luminance), there’s no way to find an overlay color that will look like you’ve blurred what it’s on top of.


There is no general algorithm for packing problems. Only some of the special cases have known, and optimal, solutions. If you are packing one shape then finding a reasonable solution is possible. Like the known cases of hexagonal packing etc.

However, if you have multiple diffenently sized objects then easy just flew out of the door. Some heurestics have been developed, since there is a very big demand for such algorithms in manufacturing. But they are not usually very fast (performant?) or simple to implement. Many of them use a genetic algorithm. A free one that you may test can be found atsvgnest.com but faster ones and ones working on shaking and collision also exist.

Anyway your question is sufficiently vague, so i dont actually know if
 you want to do packing. I mean it is perfectly possible you mean tiling, which is pretty trivial.


smallpt: ray scattering and Importance Sampling

Let me first address some misconceptions you have:

  the theory states that we need to shoot x number of rays for each
  intersection

No, the "theory" doesn't state such a thing. Note also that the paper @gallickgunner is referring to is inapplicable in this case since smallpt is based on the rendering equation and not the limited variant present in Cook's paper. Thus you actually need to refer to Kajiya's paper or Veach's path formulation. What @gallickgunner is right about is the canceling of the cosine term (but there's a missing $\pi$ factor).

Now onto the bigger misconception you have:

  Multiple rays should be scattered for each hitting point, and their contribution needs to be calculated to return the color value of the specific pixel.

There's no such requirement, in fact Kajiya proved that scattering multiple rays at each bounce is bad since you get an exponential explosion of the ray count later down the bounce tree where the actual contribution is actually very small. In MC terms you perform splitting where the contribution is small thus wasting computational power which would be a lot more useful in reducing variance in other places. To understand why this is the case, one has to dig deeper into the theory of how raytracing works and how you estimate the rendering equation.

You start with a Fredholm integral equation of the second kind:

$$L_o(x,\omega_o) = L_e(x, \omega_o) + \int_{\Omega}{f(\omega_o,x,\omega_i)L_i(x,\omega_i)\cos\theta\,d\omega_i}$$

As long as the integral converges (which happens when the operator is a contraction, or when it is not but enough energy can escape the scene for there to be a steady state), you can find the solution as the Liouville-Neumann expansion. To do this however I'll first rewrite the equation in operator form:
$$L = L_e + TL$$
$$(I-T)L = L_e$$
$$L = (I-T)^{-1}L_e = \sum_{i=0}^{\infty}T^iL_e$$
Where the last equality holds as long as $T$ is a contraction ($||T||<1$). Note that $T^0$ represents $0$ bounces, $T^1$ represents $1$ bounce, and $T^n$ represents $n$ bounces. Thus you can look at the solution as accumulating contributions from all light paths of length $0$, $1$, $2$ etc. till infinity. Now to actually write this out explicitly it's more convenient to rewrite the rendering equation in terms of integration over all scene surface points $S$ rather than the upper hemisphere (this was actually how Kajiya derived it initially, not precisely the form below but in an equivalent form):
$$L_o(x_1\rightarrow x_0) = L_e(x_1 \rightarrow x_0) + \int_{S}{f(x_2 \rightarrow x_1 \rightarrow x_0)L_o(x_2 \rightarrow x_1) G(x_1 \leftrightarrow x_2)\,dA(x_2)}$$
Where $G$ is defined as:
$$G(x \leftrightarrow y) = \frac{|\cos\theta_x\cos\theta_y|}{||x-y||}V(x,y)$$
Where $\theta_x$ is the angle between the vector $y-x$ and the normal at $x$, and $\theta_y$ is the angle betwen the vector $x-y$ and the normal at $y$ (the cosines can be written as dot products provided you have the normals). $V(x,y)$ is the visibility function and it is 1 only if there's no surface point between $x$ and $y$ (if you shoot a ray from $x$ to $y$ the first thing you hit will be $y$), and $0$ otherwise. The visibility function is implicit in the solid angle formulation inside $L_i$ (through the ray-casting function). Now let us rewrite the operator expansion in terms of what we have above, we get:
$$L_o(x_1\rightarrow x_0) = L_e(x_1 \rightarrow x_0) + \int_{S}{f(x_2 \rightarrow x_1 \rightarrow x_0)L_e(x_2 \rightarrow x_1) G(x_1 \leftrightarrow x_2)\,dA(x_2)} + $$
$$\int_{S}\int_{S}{{L_e(x_3 \rightarrow x_2)f(x_3 \rightarrow x_2 \rightarrow x_1) G(x_3 \leftrightarrow x_2)f(x_2 \rightarrow x_1 \rightarrow x_0) G(x_2 \leftrightarrow x_1)\,dA(x_3)\,dA(x_2)}}$$
$$+ ...$$
This sum continues till infinity with iterated integrals corresponding to $T^i$, the part expanded above is $T^0L_e + T^1L_e + T^2L_e$. Note that usually you have an extra integration over the area of the pixel (then there's a factor $W$ to account for the sensor sensitivity). When you do pathtracing (what smallpt does) you're actually estimating this infinite sum of integrals through Monte Carlo (this is also the reason why Monte Carlo is employed, since the convergence doesn't depend on the dimension, which in this case is infinite). You can estimate each term in the sum by Monte Carlo integration, that is for a $n$-iterated integral your samples are really paths starting from the camera with length $n$: $x_0x_1x_2...x_{n+1}$ for a path of length $n$. Note that you do not require that you split the path for each subsequent integral.
What path tracing implementations do, is create this path $x_0x_1...x_{n+1}$ incrementally, and reuse the path for earlier terms in the sum. That means that you start by picking some point $x_0$ inside a pixel, pick some point $x_1$ in the scene compute the value of the integrand of the term $T^0L_e$ for the path $x_0x_1$, pick a point $x_2$ in the scene and compute the value of the integrand of the term $T^1L_e$ for the path $x_0x_1x_2$, and so on. In an implementation you usually use the solid angle formulation (at least for bsdf sampling), but there's not major difference, you still end up with paths of the form $x_0x_1...x_{n+1}$.

The takeaway from this is that pathtracing works with paths where, a path starts at the camera and ends at a lightsource. Using different paths that do not split (in general) to estimate the sum of integrals is usually more efficient than splitting each path into more paths at each bounce. I am not saying there's no place for splitting (actually smallpt uses it for the first few bounces for materials that refract), you just have to know when to apply it (when the contribution is large and is inadequately sampled by only 1 or a few samples). Note that the opposite of splitting is early path termination in the form of the Russian roullete - basically kill off paths with small contributions since it's not worth it wasting computational resources on those (it also lets you make your estimator unbiased, but in practice there's usually a ceiling to the number of bounces allowed, since for most scenes your estimator being "fully unbiased" is not worth it in terms of performance).


Is there a way to interpolate color across the line with only integer calculation ?`

Yes, it is possible to use only integer calculations. I will describe how, but bear in mind that the difference in speed between integer arithmetic and floating point arithmetic is not as great as it was historically. If you want your code to run faster, it is best to profile and identify which parts of the code are taking up most time, before considering whether to convert a particular aspect to use integers.

A straight line through colour space

Looking at the code, it appears you first convert from RGB to HSB then create a linear gradient, and convert it back to RGB.

A linear gradient is a straight line through the colour space, so you can use an integer based straight line algorithm to generate it without the need for floating point arithmetic. For example, the integer based version of Bresenham's algorithm.

H, S, and B are being treated independently, so you can calculate each one separately. For example, the H component can be calculated by treating the line as a straight line from (x0, H0) to (x1, H1), instead of the usual (x0, y0) to (x1, y1).

If the range of y values along the line is greater than the range of x values along the line then Bresenham's algorithm will step along integer y values instead of x values. The algorithm can still be applied in exactly the same way, just using (y0, H0) to (y1, H1).

You can also combine all 4 calculations into a single loop, keeping track of a parameter each for determining y value, H value, S value and B value.

The only difference from plotting a line on the screen (using x and y) is that not every value of H, S, and B needs to be taken. If the x range is greater than the y range, the algorithm steps through each of the x values checking whether to increment y. If the y range is greater then it instead steps through the values of y checking whether to increment x. However, if the range of H is greater than the range of x, it makes no sense to step through H values, as many of them will correspond to the same x value, wasting calculations. If the range of H is smaller than the range of x then the algorithm can proceed as usual. If the range of H is often going to be larger than the range of x then it is worth considering a different approach. Similarly for S and B.

Fitting the colour space to your line

One way to eliminate this problem is to tailor your colour space to the x (or y if greater) range required for this specific line. Since your input colours and final output colours are in RGB, you are not constrained to a particular range for the HSB colours. You are free to convert from 256 levels of R, G, and B to 1000 levels of H, S, and B (or 1000000, or whatever you choose). If you choose to use a colour space where the range of H values used is the same integer as the range of x values used, then you won't even need to use a line drawing algorithm for the colours. Simply store the required values for converting back to RGB, and each of the H, S, and B values will be x plus a constant.

This will require modifying your conversion function, but this will probably be necessary anyway, as it currently uses radians, which are not going to help with keeping things integer.

More accurate colour spaces

Note that using a colour space like HSB will not give gradients which look as even as if using a more perceptually uniform color space such as Lab. However, since such spaces require calculating roots and you are looking to keep to integer arithmetic as much as possible, I'm assuming you are happy with HSB. If you did want to use something like Lab, you could still potentially keep to integers by using look up tables.

What inputs does this need to work for?

You also mention having to make decisions in the code which also takes up time. I would question how much difference this makes since the decisions are outside the loops and therefore only performed once. If you really want to eliminate them, ask yourself which inputs the code needs to work for.

Currently the code checks the start and end values of the hue and decides whether to change the hue clockwise or anticlockwise. If you eliminate this check the code will still work, but it may not always give the colour direction the user intended. Deciding this automatically also may not give the direction the user intended (since different users have different expectations, and the requirement may change from one use to the next). You could simply take the direction as an additional parameter, which would avoid having to make the decision.

Code review

Meaningful variable names

The use of "hue" in the naming of variables that store all three of hue, saturation and brightness makes the code less intuitive to read, potentially causing the reader to wonder whether saturation and brightness have been overlooked. It also makes the code difficult to reason about, since startHueInRadians[0] is a hue, while startHueInRadians[1] is a saturation, which is not only not a hue, but also presumably not in radians.

In addition to making the code more accessible to others, you may also find it easier to work with it yourself if the naming is consistent, especially if you look back at the code after a few months.

A potential bug

If one end of the colour gradient is at hue zero (red), the code checks whether the other end is greater than or equal to two thirds of the way around the hue circle. This corresponds to blue. If it is less than blue, then the gradient proceeds from zero upwards (anticlockwise), otherwise it proceeds from 2 pi downwards (clockwise). This means that a gradient from red to a colour with hue just above blue will go from red to magenta to blue. A gradient from red to a colour with hue just below blue will go from red to yellow to green to cyan to blue, which may not be what is expected. Intuitively, I would expect the cut off to be half way around the circle rather than two thirds, so that the gradient always includes the fewest colours (taking the shorter, more direct route).

If you did intend to have this asymmetry, note that choosing a similar red with hue just above or just below zero will override this behaviour, forcing the long route or the short route respectively.

If you want the gradient to consistently take the shorter route, then the code will need to check which route to take for all values, not just zero. For example, currently choosing the start hue just above zero and the finish hue just below zero will give start and finish colours that both look red, but the gradient will pass through red, yellow, green, cyan, blue, magenta and back to red, rather than just fading subtly from one hue of red to another.


Is it possible to access adjacent vertex data in tessellation evaluation shader?

If you want to prevent cracks between the tessellated lines in a single curve, what you need to know is how to miter the quads correctly. Your GS's extrusion is presumably based on taking the 2D perpendicular to the line direction and extending the quad in that direction. Well, don't do that.

Instead, have the TES compute the perpendicular direction for each vertex in the tessellation. You can compute this by taking the first derivative of the Bezier curve (its tangent), which can be computed exactly, and then finding the perpendicular of that. Each vertex will have its own perpendicular value.

It should be this perpendicular which you feed to the GS to do the quad extrusion. And you should be guaranteed to have invariant results. The TES's computation of the position and perpendicular will be based solely on the positions of the control points and the tessellation coordinates. The control points are fixed for all invocations of a TES over a patch, and the tessellator guarantees the invariance of gl_TessCoord for the "same vertex".

Invariance between patches will be something of a problem. That is, invariance between the end of one set of 4 control points and the beginning of an attached one. The tessellation system only guarantees invariance at edges between patches if the computations that lead to the positions of the vertices are all based solely on data that is binary identical between the two edge vertices. And that's not generally the case for bezier curves; you achieve C1 continuity by ensuring that the slope between boundary control points is the same. And that slope is not passed as a parameter; it is computed. And thus, two edge vertices on different patches could have slightly different slope computations, which leads to different vertices, which leads to invariance.

The only way to handle this is to have your Bezier spline defined by six control points rather than four. You pass as the first and last control points the previous/next control points for the Bezier spline (you'll have to make something up at the edges). This allows your computations to be invariant, since at the edges, you will be doing computations based on the exact same data.


OpenGL - drawing icosahedrons and cylinders and combining them

When drawing in OpenGL and related APIs you usually have one or more models. In your case, this would be the cylinder and the icosahedron. The models are generally generated at some default size which may or may not suit your needs for a particular task. You move them, scale them, and rotate them by setting up a model matrix which describes how you want to transform them.

So how would this work in your case? I recommend having 2 functions. One which draws a unit cylinder and another which draws a unit icosahedron. Say you have a data structure for your coordinates like this (I'm assuming a C-like language):

typedef struct Position3D {
    float x;
    float y;
    float z;
} Position3D;

You can then construct a triangle by making a structure for one:

typedef struct Triangle3D {
    Position3D p0;
    Position3D p1;
    Position3D p2;
} Triangle3D;

Then you'd have functions for generating the vertices you need. Something like this:

void generateCylinder (Triangle3D* triangles) // <-assumes you allocated triangles before calling and that you have enough space for all of them!
{
    int i = 0;
    const float kRadius = 1.0;
    const float kDeltaH = 0.1; // <- change this to make the mesh more or less dense
    const float kDeltaAngle = M_PI / 10.0; // <- change this to make the mesh more or less curved
    for (float h = -1.0; h < 1.0; h += kDeltah)
    {
        for (float angle = 0; angle < 2.0 * M_PI; angle += kDeltaAngle)
        {
            float x = kRadius * cos(angle);
            float y = h;
            float z = kRadius * sin(angle);
            //... add points to triangles array here...
            // p0 = <x,y,z>
            // p1 = <x, y + kDeltaH, z>
            // p1 = <x + radius * cos(angle + kDeltaAngle), y, z  + radius * sin(angle + kDeltaAngle)>
        }
    }
}

And a similar function for the icosahedron.

Next, you'll create a matrix for every object you draw and send that to your shader before drawing it. So you'll have the translations for the vertices of the (big) icosahedron, and you want to draw a (small) icosahedron at each one. You'd do the following:

float translationMatrix [ 16 ] = {
    1, 0, 0, 0,
    0, 1, 0, 0,
    0, 0, 1, 0,
    0, 0, 0, 1
};
translationMatrix [ 12 ] = xPos;
translationMatrix [ 13 ] = yPos;
translationMatrix [ 14 ] = zPos;
GLuint modelMatrixLoc = glGetUniformLocation(program, "modelMatrix");
glUniformMatrix4fv(modelMatrixLoc, 1, GL_FALSE, translationMatrix);

Then you'd draw your icosahedron. For the cylinders you would need to scale vertically, rotate to the appropriate angle and translate it to the position between 2 icosahedrons.


How is screen space ambient occlusion implemented?

As Alan and trichoplax mention in the comments, the effect that ambient occlusion simulates is not the occlusion of a surface from the camera but the surface’s occlusion from its surroundings.

Think of it this way: say you have even illumination coming from every direction, so that the total incoming light at any point sums up to a value of 1. If you place a flat plane in that environment and look at one side of it, that side is going to receive 50% of that illumination, or 0.5, because the other half is blocked by the plane itself. In other words, any point on the plane’s surface can only “see” the light coming from half of the environment, so it’s half as brightly illuminated. If you fold that plane towards your viewpoint (a “valley” fold), then you decrease the incoming illumination to that side of the plane further, to some value below 0.5, because, again, each point on the plane “sees” a bit less of the light coming from the surroundings.

Screen-space ambient occlusion works more or less by looking for these “folds”—areas where the depth changes sharply, as defined by comparing the depths of neighboring pixels—and darkening them to simulate the decreased illumination from the points’ environment.


It looks to me like the vertical dimension is the one to use for scaling:

$$\frac{50\text{ px}}{720\text{ px}} \approx \frac{30\text{ px}}{480\text{ px}}$$

There's a slight error there, but it's only a few pixels. I'm not sure if the object sizes you mentioned were intended to be precise to the nearest pixel (30 and 50 are roundish numbers).

Typically when switching between aspect ratios, the vertical FOV would remain constant while the horizontal FOV widens or narrows to match the aspect ratio.


I fixed it, It was just a coding mistake, there was nothing wrong with the implementation. So the error was in fact in the 2nd for-loop when I was parsing the image buffer. I was reaching for the last pixel in the vector and then popping it off from the back of the stack.

This means that in the for loop condition, the vector size will be changing but k will still be incrementing the same, so it will effectively only do half of the elements per vector. Because as K is incrementing stack.size() is decrementing so they meet in the middle and the loop is done.

To fix this, I just used index notation, so:

rgbpacket &rgbpixel = stack[index];

and then I deleted the line where I was popping off from the stack.

Even though this is a stupid bug,
I hope this helps other people.


What program can bake transparent light maps?

Not sure what sort of lightmaps are produced by blender. Perhaps if you can get it to render a non-transparent lightmap (for example one that has all the colours you want but on a black background), you could then use a program to remove all of the black. I believe GIMP and Photoshop both have a colour to alpha tool. Be sure to pick a background colour that will not be the same colour as any of the lights.


As you said you want to create something like this
As you know these are not graphic intense rendering.
Thus, (i guess but not 100%) Windows do these rendering on CPU not on GPU.
GPU is for intense rendering with lighting, textures, 3D models, etc. Where you can take the benefit of algorithm built in GPU circuit, Which is what OpenGL does.
So i would recommend to use built-in libraries that don't use GPU and use CPU for rendering, Although you can use OpenGL or Direct3D as well.
Some libraries that are built in different programming languages are :-

Java - Jframe, Graphics, Buffered image,  Canvas. (Note : all these libraries have different functionality and you have to use them together for software rendering See this for more detail)
C# - System.drawing,  (i don't have much experience of C#, So i can't name any more libraries), also C# form application.
JavaScript - Canvas check this out
C++ - you need third party libraries like QT, GTX etc. Qt learning resource, related resource
Python - Pygame (No personnal eperience) See this for detail


What is performance better - add four materials to mesh or duplicate every vertex and use packed texture?

Hello and welcome Jonek2208!

Unity specific questions probably are better suited for the Unity Graphics Forum. However, I'll try to give you an answer as I do have some experience with the engine.

In my opinion, the only way of knowing for sure which method is the best would be to implement both and profile them on your target device. Do not fool yourself in thinking that some basic algorithmic analysis is all you need to do. Always profile and compare your solutions, no matter the circumstances.

One thing I can do is to tell you the trade-offs between the two methods, and also provide a hint at another possible solution. Let's go!

When using the built-in rendering pipeline, Unity batches objects by materials. A batch is considered to be a change in GPU state followed by one or more draw calls. In other words, if you have four materials, you increase the number of batches when compared to having just one.

If you are targeting mobile devices, reducing the batch count is one of the keys for achieving good performance, and most mobile developers would choose to use a single material. On the other hand, duplicating vertices incurs a cost in memory, and hence you get better batching at the expense of using more memory. If you are targeting consoles or PCs, you can have a much higher batch count and memory budget, so the solution that you think is simpler to implement, expand, and maintain should be your best bet.

Another solution would be to use textures instead of vertex colors. For this, you would need to compute the uv-map procedurally when generating the vertices. You could place all textures into an atlas and try to use a single material. This is the usual solution for mobile. If you are targeting PCs/consoles, you could even afford procedural textures in your fragment shader.

If you are using the latest version of Unity, you can use the SRP Batcher, which batches per shader instead of per material. Hence, you can use 4 materials, but have a single batch if all of them use the same shader. Note some devices may not have this capability, especially on mobile.

All in all, there are a lot of things to consider in any kind of application. That is why I started talking about profiling on your target device. Honestly, it is the only way of knowing for sure which method is the best. You should build a test scene with thousands of generated meshes and then profile each method. Select the one that fits best in your budget (FPS/memory/etc) and the long-term goals of your project.

Good luck! :)


Implementing the example of Walter et al. (2007) paper

Initially and still I'm not so sure how to get the color of the object behind the glass into the formulas.

Normal, this was written with path / ray tracers in mind where that part is easy. Your idea of a cube map might be ok for some use cases (eg. reasonably rough surface for the distance to the surrounding objects). It is not a general solution however. That would likely require GPU ray tracing.

  Going by my first thought again I just multiplied the result of the BTFD with the sample of the o_t vector but apparently the BTDF was really big so the whole slate turned white (possible bug?).

I'm writing this from memory but I think you need to divide by the PDF. Which means that your weights will be close to 1. But not quite 1 as the sampling distributions don't include the shadowing terms.

  Also I think I read somewhere that the reflection and transmission should be handled seperatly depending on the result of the Fresnel formula. Or should it be handled at the same time (possibly somethink like : BRDF * sample_of_reflection + BTDF * sample_of_refraction)?

A classic path tracer would randomly pick either reflection or refraction based on the fresnel coefficient. That does efficient sampling. But if you can easily get both from a cube map, you can also blend them with the fresnel coefficient.

  After some playing around with this function I noticed that is basically something akin to a gaussian distribution pointing towards the positive z-axis. But shouldn't the distribution give samples which are distributed around the surface normal?

Yes, you need to transform the generated direction. It is mentionned in section 5.2, after equation (24), that $\theta_m$ is the angle between $m$ (microfacet normal) and $n$ (surface normal). They don't go into great detail as this way of defining microfacet models goes back to much older articles.

The short version is that you build an orthonormal basis using the surface normal as z and two other perpendicular vectors as x and y. Here it does not matter which as the distribution is isotropic, as long as they are perpendicular to each other and to the normal. Then you transform using that basis.


About the use of Russian Roulette in Smallpt

The reflected and transmitted values need to be scaled because they are not sampled proportionally to their contribution. You'll notice how the decision whether to reflect or transmit is made based on P, not Re.

About why that scaling has been done in the first place, only the author can give you the definitive answer. My guess is that this is done to prevent extreme values that would come from tiny probabilities. Imagine that Re is 0.001 and is used without scaling. Then only one in 1000 rays would reflect, and the contribution of that ray would be multiplied by 1000. At a render with, say, 64 rays per pixel, this would leave very visible artefacts, as not every pixel would get a reflected ray. The scaling from this code snipped on the other hand makes sure that the ray budget is more evenly distributed between reflection and transmission and that each potential path taken is guaranteed at least 25% or the rays.

For further reading about path tracing, there is no better and comprehensive resource than www.pbrt.org - both the book and the source code are worth reading.


Triangulated model rasterization: light-weight Python library

I've found VTK library and used bindings to Python  3 available from different Anaconda organizations: anaconda itself, conda-forge, clinicalgraphics and menpo.
Also it can be compiled from sources.

Works fast (I've rendered a mesh with half a million vertices without any problems) and easier in use than low-level OpenGL shaders with GLUT.

menpo has its own wrapper for VTK and mayavi called menpo3d.
I have not tried it yet, but examples look very simple.


What methods/technologies to reduce required performance for virtual reality are there?

Avoid stereo when possible
There was some research published recently measuring when users can or cannot tell whether the specular contribution, which is view point dependent, is different between eyes.
Perception of Highlight Disparity at a Distance in Consumer Head-Mounted Displays
In the future, this could be use to share some computation between the two eye renders. Unfortunately this sounds hard to take advantage of in a typical rasterization pipeline.
Some games already avoid stereo when rendering far objects and background.

Reduce shading when possible
Lens distortion
Part of the NVidia VRWorks SDK is a feature, Multi-Res Shading, that allow to reduce the amount of shading, to take into account the fact that some pixels contribute less to the final image due to lens distortion.
In Alex Vlachos' GDC 2015 presentation, he also mentions the idea of rendering the image in two parts of different resolution, to achieve the same goal.
Foveated rendering
Finally, like your question mentions, foveated rendering aims at scaling the resolution to take into account the fact that eyes have more precision in the fovea and less in the periphery. There is some research on the topic, but this requires eye tracking like in the Fove HMD.


Can one VAO store multiple calls to glVertexAttribPointer?

In version 4.3+ the vao layout can be set independently from the VBO buffer which allows you to rebind the vbo without having to call vertexAttribPointer every time:

//in init
glBindVertexArray(vao);

glVertexAttribFormat(posAttrLoc, 3, GL_FLOAT, false, offsetof(Vertex, pos));
// set the details of a single attribute
glVertexAttribBinding(posAttrLoc, 0);
// which buffer binding point it is attached to
glEnableVertexAttribArray(posAttrLoc);

glVertexAttribFormat(normalAttrLoc, 3, GL_FLOAT, false, offsetof(Vertex, normal));
glVertexAttribBinding(normalAttrLoc, 0);
glEnableVertexAttribArray(normalAttrLoc);

glVertexAttribFormat(texAttrLoc, 2, GL_FLOAT, false, offsetof(Texture, tex));
glVertexAttribBinding(texAttrLoc, 1);
glEnableVertexAttribArray(texAttrLoc);

Note that there is no mention of any VBO during init, just the binding points. That binding point is then used during draw to specify which vbo to pull the vertex data from:

//draw

glBindVertexArray(vao);
foreach(mesh in meshes){
    glBindVertexBuffer(0, mesh.vbo, mesh.vboOffset, sizeof(Vertex));
    glBindVertexBuffer(1, mesh.texVbo, mesh.texVboOffset, sizeof(Texture));
    // bind the buffer to the binding point
    glDrawArrays(GL_TRIANGLES, 0, mesh.vertexCount);
    //draw
}

In pre 4.3 you would create a VAO per mesh and per mesh call glVertexAttribPointer for each attribute:

glBindVertexArray(mesh.vao);
glBindBuffer(GL_ARRAY_BUFFER, mesh.vbo);

glVertexAttribPointer (posAttrLoc, 3, GL_FLOAT, false, sizeof(Vertex), mesh.vboOffset + offsetof(Vertex, pos));
glEnableVertexAttribArray(posAttrLoc);
//...

Then during draw bind the per mesh vao before drawing it:

foreach(mesh in meshes){
    glBindVertexArray(mesh.vao);
    glDrawArrays(GL_TRIANGLES, 0, mesh.vertexCount);
}


When your GPU displays a new frame on the screen, it transfers the image over the HDMI cable (or whatever kind) in a process called "scanout".  The pixels are sent out in linear order, usually left-to-right and top-to-bottom.  The process is timed so that it takes most of the duration of a refresh interval to do this.  For instance, at 60Hz, one frame is ~17 ms.  Each scanout will take probably around 15-16 ms, with 1-2 ms of vblank in between (the exact values vary according to the display and video mode).

Traditionally, rendering is double-buffered, which means there are two buffers stored in GPU memory: one that is currently being scanned out ("front buffer"), and one that is being rendered to ("back buffer").  Each frame, the two are swapped.  The GPU never renders to the same buffer that's being scanned out, which prevents artifacts due to potentially seeing parts of an incomplete frame.  However, a side effect of this is increased latency, since each frame may sit around in the buffer for several ms before it starts being scanned out.

VR is very latency-sensitive, so this isn't desirable. An alternative approach is to render directly to the front buffer, but time things out very carefully so that you have rendered each line of the image shortly before the scanout gets there. That's called "scanline racing" or "racing the beam" (the "beam" harkening back to the CRT days of yore). This more or less requires that you render the image in scanline order, i.e. the same order that the pixels get scanned out.  It doesn't literally have to be rendered one line at a time—it could be rendered in thin strips a few pixels high, but it does have to be done in order, as you can't go back and edit pixels that have already been scanned out.

There are a lot of disadvantages to this approach; it has very stringent performance requirements, has to be timed very carefully against vsync, and it greatly complicates the rendering process. But in principle it can shave milliseconds off your latency, which why VR folks are interested in it.


PBR : Understanding the right part of Split sum integration of specular IBL

1)What is the V vector? How did they compute that

It's the vector toward the camera (view vector), i.e. direction of the reflected ray. The lookup table they're building is parameterized in terms of NdotV, so they are working backward from NdotV to obtain a corresponding V vector. For purposes of doing this integration it doesn't matter what N and V are exactly, just that they have the correct NdotV. So, just to make things concrete, they assume that $N = (0,0,1)$ and V is aligned with the positive x-axis (i.e. in the positive xz half-plane). Then you have $V = (\sin θ, 0, \cos θ)$ where θ is the angle from N to V (this is just the equation for a circular arc from z down to x). Now $N \cdot V = \cos θ$, so by a standard trig identity $\sin θ = \sqrt{1 - \cos^2 θ} = \sqrt{1 - (N \cdot V)^2}$.

2)Why have they assumed N=(0,0,1)?

This whole thing is a pre-integration of the BRDF alone, it doesn't depend on any scene or lighting environment. So it doesn't matter where the BRDF is or how it is oriented in space, only the relative orientation between the N, L, and V vectors is important. So just for the sake of convenience they have chosen to fix N = (0, 0, 1), let V vary only in the xz half-plane as mentioned above, and let L vary over the whole hemisphere. To put it another way, they have chosen to set the coordinate system to be aligned with N and V.

3)why is NDotL=max(L.z,0)? How did they simplifying a dot product to an min/max Operation same question I have for NDotH

Since N = (0, 0, 1), dotting N with any vector just takes the z component of that vector. $N \cdot L = 0L_x + 0L_y + 1L_z = L_z$. (This is why it's convenient to assume a coordinate system in which N = (0, 0, 1).) Taking the max with 0 clamps it to the upper hemisphere, otherwise the dot product would go negative if the vector happens to be pointing below the hemisphere. This clamp is a standard thing to do when dealing with opaque BRDFs to avoid accidentally getting unwanted negative values in your calculations— although it is not really necessary here since they check if (NdotL > 0.0) immediately afterward.

4)G_Vis=GVDoth/(NDotHNDotL) why the division by (NDOTH*NDOTL)? that was nowhere in the equation.

It is a standard part of the microfacet BRDF formulation:
$$
f_r(p, w_i, w_o) = \frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}
$$
Note that $\omega_o = V$, and $\omega_i = L$, these are just different notations/names for the incoming and outgoing light vectors.
However I think there is a typo in the code, the denominator should contain NdotV rather than NdotH. I am also not sure why there is a multiplication by VdotH in the numerator, or where the factor of 4 went, but perhaps they have accounted for it elsewhere. I would probably double-check the BRDF equations here against some other sources as I don't think it's entirely correct.

5)what is vector L? it looks like they have reflected vector V about H to get this but why?

As noted, L is the vector toward the incoming light. H is the "half-angle" vector exactly halfway between L and V, so reflecting V around H will give you L, or vice versa. The significance of H is that in microfacet BRDFs, it is the normal vector of those microfacets that are "active" for a given L and V. So, to importance-sample a BRDF, we sample a random normal vector from the NDF, and use that as H. Since H and V are specified, we then calculate L from those.


It seems that sampler2D_float doesn't allow to interpolate shadow lookup  linearly. So I had to do it by hand. Here's an example of interpolated shadowing.

float texture2DCompare(sampler2D depths, vec2 uv, float compare){
    float depth = texture2D(depths, uv).r;
    return step(compare, depth);
}

float texture2DShadowLerp(sampler2D depths, vec2 size, vec2 uv, float compare){
    vec2 texelSize = vec2(1.0)/size;
    vec2 f = fract(uv*size+0.5);
    vec2 centroidUV = floor(uv*size+0.5)/size;

    float lb = texture2DCompare(depths, centroidUV+texelSize*vec2(0.0, 0.0), compare);
    float lt = texture2DCompare(depths, centroidUV+texelSize*vec2(0.0, 1.0), compare);
    float rb = texture2DCompare(depths, centroidUV+texelSize*vec2(1.0, 0.0), compare);
    float rt = texture2DCompare(depths, centroidUV+texelSize*vec2(1.0, 1.0), compare);
    float a = mix(lb, lt, f.y);
    float b = mix(rb, rt, f.y);
    float c = mix(a, b, f.x);
    return c;
}


Ask for detailed derivation of a formula in "Advanced Global Illumination"

The only thing you really need to know in order to derive (2.22) from (2.21) is that the $\delta$-distribution satisfies

$$ \int \delta(x) f(x)\,\text{d}x = f(0). $$

(The $\delta$-distribution is not a function in the common sense, so the above integral is merely a useful notation and not to be read as e.g. a Lebesgue or Riemann integral.)

Unfortunately, the notation that the authors use for the substitution of the incident radiance distribution is misleading. The $\Theta$ in $L(x\leftarrow\Psi) = L_{in}\delta(\Psi-\Theta)$ is completely unrelated to the integration variable $\omega_\Theta$, so I chose to denote it by $\tilde\Psi$ instead.

Here is the derivation of (2.22) from (2.21):

Substitute $L(x\leftarrow\Psi) = L_{in}\delta(\Psi-\tilde\Psi)$:

$$\frac{\int_{\Omega_x} \int_{\Omega_x} f_r(x, \Psi \rightarrow \Theta) L_{in}\delta(\Psi-\tilde\Psi)\cos(N_x,\Theta)\cos(N_x,\Psi)\,\text{d}\omega_\Psi\,\text{d}\omega_\Theta}{\int_{\Omega_x}L_{in}\delta(\Psi-\tilde\Psi)\cos(N_x,\Psi)\,\text{d}\omega_\Psi}$$

Substitute $\Psi$ with $\Psi+\tilde\Psi$:

$$\frac{\int_{\Omega_x} \int_{\Omega_x-\tilde\Psi} f_r(x, \Psi+\tilde\Psi \rightarrow \Theta) L_{in}\delta(\Psi)\cos(N_x,\Theta)\cos(N_x,\Psi+\tilde\Psi)\,\text{d}\omega_\Psi\,\text{d}\omega_\Theta}{\int_{\Omega_x-\tilde\Psi}L_{in}\delta(\Psi)\cos(N_x,\Psi+\tilde\Psi)\,\text{d}\omega_\Psi}$$

Use the above property of the $\delta$-distribution:

$$\frac{\int_{\Omega_x} f_r(x, \tilde\Psi \rightarrow \Theta) L_{in}\cos(N_x,\Theta)\cos(N_x,\tilde\Psi)\,\text{d}\omega_\Theta}{L_{in}\cos(N_x,\tilde\Psi)}$$

Simplify:

$$\int_{\Omega_x} f_r(x, \tilde\Psi \rightarrow \Theta) \cos(N_x,\Theta)\,\text{d}\omega_\Theta$$

Now you only need to notice that we want equation (2.21) to be true for all such "test functions" $L_{in}\delta(\Psi-\tilde\Psi)$ for the incident radiance, i.e. for all possible values of $\tilde\Psi$.


Why Diffuse Light use max(N · H, 0) instead of just letting it be negative?

If lighting with multiple lights, they add together to make the final lighting value. A negative light would darken other lights, which is incorrect.  This is also true if you add in specular reflections, emissive lighting, or other sources of lighting. The negative lighting from one source of lighting would darken some other source of non negative lighting.


What to reproject when doing temporal filtering in stochastic screen space reflection?

Usual reprojection finds UV of current visible point into previous frame. For diffuse and rough objects it simply adds motion vector: vec2 UV_reprojected = UV + motion;
But for mirror-like surfaces it doesn't work, because we see not object (mirror) itself but what reflects in mirror (virtual object). The general idea to reproject reflected point is:

Find world space position for reflected point visible in reflector (Here we need to know reflected path length).
Substract world space delta between frames, so we found previous world position.
Project this point into reflector using previous camera position.
Project world space reflector point back to the screen.

The most hardest part is p.3, and is not easy analytically find this point, because reflector can change normal. Some interesting links about this topic:

Reflection reprojection in Uncharted. Here we assume that reflector is plane that changes transform, both camera and receiver are also dynamic.
Ray Tracing Gems chapter 32. (32.4.3) Theoretical approach that finds reflector position in terms of least squares. And interesting geometrical optic approach (here we assume that reflector point is locally sphere-like).


Rather than fit bilinear patches to your sample points, which are only C0 continuous$^\spadesuit$ may I instead suggest fitting the surface using Catmull-Rom splines? These will give you C2 continuity.
Catmull-Rom is an interpolating spline which means it will pass through your sample points, as did the bilinear surface, but will do so smoothly from patch to patch.

$^\spadesuit$ because the first derivatives won't generally match at the boundaries of your patches

Background/Assumption

You have a function $$F:\Re^2\to\Re^3$$ of the form  $$F(s,t)=[s,t,f(s,t)]$$ where $f:\Re^2\to\Re$. The function, $f$, can be considered a 'black box' in the sense that we might not be given it explicitly.

$F$ has been regularly sampled as $$S_{i,j}=F(x_{scale}.i, y_{scale}.j, f(x_{scale}.i,y_{scale}.j))$$ where $i,j\in\mathbb{Z}$.  You then want interpolate between these sample points using piecewise bicubic Catmull-Rom patches.

Constructing patches

To form a particular bicubic Catmull-Rom patch, $C_{i,j}$, which has a given 2x2 set of samples $\{S_{i,j},  S_{i+1,j}, S_{i,j+1},  S_{i+1,j+1}\}$ at the corners, you need the set of 4x4 samples formed from the 2x2 at the centre and the surrounding 12 samples, i.e.$\{S_{i-1,j-1}$ through $S_{i+2,j+2}\}$,  as shown below.
(My apologies that this is 'hand placed' so is far from ideal)

If we just consider the behaviour of the surface at sample, $S_{i,j}$, then the partial derivative in the $i$ direction is given by
$$\frac{\partial C_{ij}}{\partial i}=\frac{S_{i+1,j}-S_{i-1,j}}{2}$$
and similarly, in the $j$ direction by
$$\frac{\partial C_{ij}}{\partial j}=\frac{S_{i,j+1}-S_{i,j-1}}{2}$$

Remembering these are vectors, you can use the cross product to generate a surface normal at any of the $S_{i,j}$ points$^\clubsuit$. (Though you may need an extra set of boundary samples)

  Ali: You mentioned that both partial derivatives are vectors(in i,j directions). could you write each vector components plz

Certainly: If we say $S_{i,j}=[x_{i,j},y_{i,j},z_{i,j}]^T$ then $$\frac{\partial C_{ij}}{\partial i}=\frac{1}{2}.[(x_{i+1,j}-x_{i-1,j}),(y_{i+1,j}-y_{i-1,j}),(z_{i+1,j}-z_{i-1,j})]^T$$
$$\frac{\partial C_{ij}}{\partial j}=\frac{1}{2}.[(x_{i,j+1}-x_{i,j-1}),(y_{i,j+1}-y_{i,j-1}),(z_{i,j+1}-z_{i,j-1})]^T$$

If you've sampled the surface 'finely' enough you could probably choose to render the surface by just using 2 or 4 triangles per original patch.

Alternatively, you can recursively subdivide your quad mesh to get a finer set of triangles. To double the number of points in the 'i' direction, simply insert a new point between each $S_{i,j}$ and $S_{i+1,j}$ that equals ...
$$\frac{9(S_{i,j}+S_{i+1,j}) - (S_{i-1,j}+S_{i+2,j})}{16}$$
 Once you have the finer grid in 'i' simply repeat with 'j' with the appropriate equivalent set of points. You can repeat this ad nauseam to get whatever triangle resolution you desire ;-)

$^\clubsuit$ (Copying some text from the comments below): I've made a couple of valid assumptions here.  Noting that if a given first (partial) derivative is non-zero (which is guaranteed in your case because of the way you construct the surface), then it also is locally tangential to that surface. Further, we know your two partial derivates aren't parallel, so the cross product will construct a valid normal.
FWIW For some arrangements of control points, i.e. in more general use cases, you can get so-called 'degenerate' surfaces where you can't naively apply the above rules, but you're safe here.


Ray Tracing Shadows: The Shadow Line Artifact

The shadow terminator problem (as it's also known) is still something of an open problem in ray-tracing. As you say, disabling self-shadowing is not a very good workaround, because non-convex meshes often need to self-shadow. There's a possible refinement to this method. Instead of throwing away all shadow ray intersections with the same object, you can only throw away the intersection if it's with the same object and it's through the reverse face of the triangle (i.e. discard the intersection if $d\dot\ N > 0$, where $N$ is the geometric normal of the intersected polygon and $d$ is the shadow ray direction).

One alternative that I've used successfully is to offset the origin of the shadow ray a short distance along the surface normal. This way, the shadow ray originates not at the geometric surface of the mesh (which is only an approximation) but at the hypothetical surface of the real object.

In theory, you should choose the distance along the normal according to the point's distance from that polygon's vertices, so that the origin is on the approximated surface, but in practice, a small fixed offset works for almost all situations.

You should also be sure to use the normal rather than the geometric surface to decide whether the light is "behind" the surface. That is, if $N\dot\ L < 0$, the shading point is in shadow (without needing to cast a ray), where $N$ is the interpolated and shaded normal.


Does the blending matrix change between calculating various curve segments in a uniform cubic B-splines approximation?

No, the $B$ matrix (basis coefficient matrix) does not change from one segment to the next. It is a property of the type of spline you're using, in this case cubic B-splines. If you used Bézier splines or Hermite splines instead, you'd have a different $B$ matrix.


Errors on Compiling my first program on OpenGL on Ubuntu

You need to link with the OpenGL library. Try adding -lgl to the compilation command.


Why random monte carlo sampling instead of uniform sampling?

Sample locations with a uniform pattern will create aliasing in the output, whenever there are geometric features of size comparable to or smaller than the sampling grid. That's the reason why "jaggies" exist: because images are made of a uniform square pixel grid, and when you render (for example) an angled line without antialiasing, it crosses rows/columns of pixels at regular intervals, creating a regular pattern of stairstep artifacts in the resulting image.

Supersampling at a finer uniform grid will improve things, but the image will still have similar artifacts—just less badly. You can see this with MSAA, as in this comparison image from an NVIDIA presentation on temporal antialiasing:

The 8x MSAA image (which is not quite a grid but still has a repeating pattern) still plainly has jaggies in it, although they're antialiased jaggies. Compare to the TXAA result, which both has a higher effective sample count (due to temporal reuse) and uses a Gaussian rather than box filter for accumulating the samples.

On the other hand, random sampling produces noise instead of aliasing. There is no pattern to the sample locations, so no pattern to the resulting errors. Both aliasing and noise are errors due to not having enough samples to form a clean image, but arguably noise is the less visually-objectionable artifact.

On the other other hand, perfectly random sampling (in the sense of i.i.d. random variables) tends to exhibit a certain level of clumping. Purely by chance, some areas in the domain will have denser-than-average clumps of samples, and other areas will be lacking in samples; those areas will be, respectively, over-represented and under-represented in the resulting estimate.

The convergence rate of the Monte Carlo process can often by improved by using things like stratified sampling, low-discrepancy sequences, or blue noise. These are all strategies to generate "de-clumped" samples that are spaced out a bit more evenly than i.i.d. samples, but without creating any regular patterns that might lead to aliasing.


What is the name of these software on 2 screen?

Judging by the interface that would be blender:

https://www.blender.org/

Please note that here is a bender-specific stackexchange site as well:

https://blender.stackexchange.com/


How to get polygon coordinates in screen space

Edit: This answer is only helpful in 3D

If you want to do it geometrically...

The inverse of the view-projection matrix, $K^{-1}$ is the matrix you want. Where $K = View * Projection$

If $\vec{v}$ is a point in homogeneous screen coordinates, $[x, y, 1]$ where $x$ and $y$ are whatever coordinates your $K$ matrix projects onto.

$$\vec{d} = K^{-1} * \vec{v}$$

Where $\vec{d}$ now represents a line in homogeneous coordinates, coming out of your camera origin. Now you'd iterate through the objects in view, and choose the closest one in front of the camera that this line intersects with.

But there are other ways!

Color Picking

Here, we render (into a buffer) each object in view with a slightly different color, which corresponds to the object ID. We then read it using glReadPixels(mouse_coordinates...) and decide which object we clicked. This is pretty easy to do in the new pipeline - you just write a color passthrough shader. Would not recommend pre OpenGL 3.0 (You'll have to do shenanigans to render the scene without interference)

Selection Buffer (Old pipeline)

If you're on the fixed pipeline, you get a nice little toy called the selection buffer to do the work for you.