Calculate set of rectangles covering pixel diffs? I'd like to develop a system for web animation similar to animated gifs, like this one. I capture a series of screenshots at fixed time intervals, and then I'd like to animate them. I can trivially calculate the pixels which have changed from one frame to another, but I'd then like to calculate a set of rectangles such that blitting those rectangles over the previous frame results in the new one. There's a tradeoff between the number of rectangles and the number of "wasted" pixels, i.e. the number of pixels stored and copied which have not changed between the two frames. Is there a good algorithm for this which would give me control over the tradeoff parameter in some form?

I would say it’s a weighted set covering problem.

In specific, the cost for a set of points is the number of extra pixels in the bounding rectangle of this set.

If the number of diff points is small, you can solve this by linear programming. Otherwise, try some approximated algorithms e.g. greedy.


Vulkan: using gl_PrimitiveID without geometryShader feature When compiling a glsl fragment shader which reads gl_PrimitiveID with glslc compiler I get the spirv-1.0 with Geometry capability enabled, which is technically correct, as it's stated in spirv spec.
Vulkan 1.0 spec says that Geometry capability must be supported if geometryShader feature is enabled.

So I assume that unless geometryShader is enabled, one can't use gl_PrimitiveID. Specifically, using gl_PrimitiveID is impossible on moltenVK, as it does not support this feature, isn't it?

In my pipeline I used gl_PrimitiveID to fetch per-face data from global ssbo and to render face ids to attachment for later reuse. Is it possible to achieve similar result without gl_PrimitiveID?

It seems likely that the mechanism by which the gl_PrimitiveID is generated is, in some way, based on the hardware that makes geometry shaders possible.

OpenGL ES 3.1 contains compute shaders but not geometry shaders (or tessellation). But it doesn't have gl_PrimitiveID. By contrast, ES 3.2 contains GS's and gl_PrimitiveID. The same goes for desktop OpenGL 3.1 vs. 3.2 and its GS functionality. Having a GS and having the ability to identify primitives in a sequence seem to be part of the same general functionality.

  In my pipeline I used gl_PrimitiveID to fetch per-face data from global ssbo and to render face ids to attachment for later reuse.

Without the hardware functionality, there is no general way to compute a count of primitives. There are only specific solutions for specific circumstances.

For example, if you are using non-index rendering of triangle lists, then you can use gl_VertexID / 3 to compute a primitive ID. Obviously that's a pretty tight restriction, but it would work in that circumstance.

You could also provide each vertex with a per-vertex face ID. The VS would output this value using flat interpolation. By carefully matching the face ID with the provoking vertex for the primitive topology, you can effectively assign a face ID to a face. This works best with indexed triangle lists, as you can play around with the order of indices in a triangle more than you can with strips. You basically have to decide which vertex maps to which face, and make sure that that face uses that vertex as the first one.


How to create a patterned object with variable pattern density in real-time? Assume there is some value p, calculated per-frame, that varies continuously over the surface of an object. The value of p determines the density of some pattern on the surface. For example, in a case with only two possible densities if p < 0.3 it is high density, otherwise it is low.

I've thought of a simple high-level solution: create two textures, each of different densities, and based on the value of p sample from the appropriate one. However, there's a problem with the boundary between high- and low-densities.

Here is an example to illustrate the problem (note my problem is NOT exclusive to this example pattern of dots. I describe the patterns I'm working with later on):

And here is the threshold between low- and high- (displayed on the high density texture but that's not relevant.) If under the line, it implies the high-density texture should be sampled.

And finally here is the comparison between what is desired and what would actually happen using this method:

The problem is that when a high-density-only circle crosses the line, it will be ignored when p indicates the low-density texture to be sampled, resulting in a truncated circle. I don't know how to solve this problem because p varies every frame, so I can't just 'bake' a boundary between the two densities. It's easy to prevent the reverse problem by creating the high-density texture from the lower one (i.e. if a circle is on the low-density texture ensure it is on the high-density texture.)

I'm interested if anyone has a way to improve my solution or even has another method entirely. The constraint here is that p is calculated per-frame in real-time. Another constraint is related to the pattern texture: the pattern is black and white, where black is the pattern and white is the background (like the circles in the example). The pattern may not just be identical shapes repeated, but any arrangement of arbitrary black shapes over a white background. (Maybe pattern is the wrong choice of word.)

I'm not familiar with the research in this field so I wasn't sure which keywords to search, so I'd even appreciate if anyone could point me in the right direction.

In order to ensure that the pattern shapes are always either wholly present or absent, never cut off, it's necessary to ensure that the same p value is used for all texels within the shape. In your example of circles, all the texels in a given circle need to agree on p.

I assume that you have some way of evaluating p at a given point on the surface (whether it is looked up from a texture or calculated from some function). Then one way to ensure a group of texels all get the same p value is to ensure they all look it up from the same point.

The UVs of this evaluation point could be stored in extra channels of the pattern texture. For instance, you could have the red and green channels store the UV coordinates at which to evaluate p, the blue channel store the threshold at which to turn on that pattern element, and the alpha store the antialiased gray level of the pattern to display. The UV+threshold data could also be in a separate, secondary texture if desired.

To generate this UV+threshold texture, starting from an input pattern texture, you could programmatically find connected components (e.g. by searching for black pixels and flood-filling). Set the evaluation point for all texels in each component to the UV of the component's center, and generate a random threshold for it. Then, when rendering, use a pixel shader that first samples this texture, then looks up p at the given evaluation point and compares it to the given threshold.

That way, each pattern shape will see a uniform p value and threshold, and will either turn on or turn off fully. As p increases, more of the shapes will pass their threshold and appear, giving the impression of a continuously varying density.


Creating vector shapes using only C++ I've been researching for quite a while and I couldn't find any explanation or at least a how to create a vector graphic application (specifically in windows) where I can create a 3D vector shape, say a box. I know there is OpenGL and DirectX for that and the new Vulkan API but I don't wan't to use those APIs and libraries. What I am asking is where can I find examples or at least tutorials for creating custom made API graphics library for the C++ language. If there is not a way in C++ language then any other way to do it. By example masm or nasm Assembly language.

If my question is too broad please specify any additional information you need to be able to answer my question.

What you are looking for is a way to get access to individual pixels in an image, in a way that you can modify those pixels with CPU code.

Once you have that you'll want to find software rendering / software rasterization methods.

A great place to start would be reading up in the various bresenham algorithms.

There is a lot of info out there on that. Here is one article for instance:
http://blog.demofox.org/2015/01/17/bresenhams-drawing-algorithms/


Is there any free software to visualize a texture mapped to a sphere? I am trying to study the texture deformations caused when mapping a rectangle to a sphere, if I had a software that took an image as an input and displayed how it looks like when texture mapped to a sphere, it would allow me to observe and study this deformations better.

Although this is borderline "off-topic" for the site, there is perhaps a simple answer:

Why not use a webgl approach such as this one?
http://learningwebgl.com/lessons/lesson11/index.html
(parent page here: http://learningwebgl.com/blog/?p=1253)

Just rotate the image of the moon around until you have one of the poles visible and you'll see the typical artefacts you get with mapping a rectangle to a sphere.
Though I've not tried, you could probably modify the code (use, say, "view page source" to read it, or save a local copy and edit) to load whatever texture your heart desired.


Compute the interpolation factor between the SLERP of 2 quaternions, given those quaternions and the SLERP result I am looking for a way to calculate the value of t of quaternion SLERP by having three quaternions.
Forward direction:
Eigen::Quaternionf q = q1.slerp(t, q2);

where q, q1 and q2 are quaternions from the Eigen library. And t is a float which interpolates between q1 and q2.
I am interested in the backward calculation, where q, q1 and q2 is given and the output is t.
To make it easier all quaternions are normalized and also lying on the same plane which crosses the origin. See the illustration below:

qdoes NOT need to be within q1 and q2.
I found this link here but I don't know how to use that.

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


Pipeline to create Voronoi Meshes I would like to implement a Maya plugin (this question is independent from Maya) to create 3D Voronoi patterns, Something like

I just know that I have to start from point sampling (I implemented the adaptive poisson sampling algorithm described in this paper).

I thought that, from those points, I should create the 3D wire of the mesh applying Voronoi (I tried to use (Python) scipy.spatial.Voronoi but the result was something different from what I expected).

I am missing something? Can anyone suggest the proper pipeline and algorithms I have to implement to create such patterns?

[EDIT] Here are a few example of what I get handling the result i get from scipy.spatial.Voronoi like this (as suggested here):

vor = Voronoi(points)
for vpair in vor.ridge_vertices:
    for i in range(len(vpair) - 1):
        if all(x >= 0 for x in vpair):
            v0 = vor.vertices[vpair[i]]
            v1 = vor.vertices[vpair[i+1]]
            create_line(v0.tolist(), v1.tolist())

The grey vertices are the sampled points (the original shape was a simple sphere):

Here is a more complex shape (an arm)

http://www.cs.sandia.gov/~samitch/papers/vor_final.pdf There isn't the computer algorithm in any programming languague, but u should be able to replicate it easily with some Plane Reflections and the algorithms from the above link.


Global illumination leading to image variation I have a 3D scene I use to generate some images which I then process in Matlab. I noticed that when global illumination (GI) is not enabled, rendering always generates exactly the same pixel-for-pixel image (as far as I can tell). However, when I enable GI, the image is not exactly the same. I have been told (I originally accidentally posted here) that GI is stochastic and so this is to be somewhat expected. However, I am curious if there's any way to avoid it (while still using GI).

Thanks

If you can make your samples deterministic (the same from run to run - don't base random sample points on time or any other non deterministic value), it will be a stable rendering that you get every time.


Normal Interpolation for Phong shading As the title says I have been trying to understand the normal interpolation for the Phong shading/lighting model. I am unsure of the equation to calculate the normal but I have come across this equation:

N=Na*α+Nb(1−α)

Where N is the interpolated vector, Na and Nb are unit vectors to interpolate between and I'm assuming α is the percentage between the two. However in using this model as shown here:

public static Vector3 lerp(Vector3 a, Vector3 b, double alpha){
     return (a.multiply(alpha)).add(b.multiply(1-alpha)).normalize();
}

I am not obtaining the correct results, so any advice on the topic, such as what I may be doing wrong etc, would be appreciated.

EDIT: I found my mistake in the interpolation, here is the updated method for those who are interested:

public static Vector3 lerp(Vector3 a, Vector3 b, double gradient){
    return a.add((b.minus(a)).multiply(gradient));
}

Your first version was correct, except that alpha and 1-alpha should be swapped (the result should be a when alpha == 0 and b when alpha == 1).

For Phong shading, you do want to keep the normalize on the end; that's the key thing that makes Phong different from Gouraud shading. But also, normalize isn't part of the usual definition of "lerp" (linear interpolation), so I would keep it out of the lerp function, and move the normalize to the Phong shading function where you call lerp.


Why is the approximation valid, in the formula provided by Brian Karis? I'm sorry I'm bringing this topic up again, but I need to expound some of the topics.

With reference to this question, I was wondering if someone can help me out in expanding some of the details.

In the answer to that question we have:

  However, since $L_i(l_k)$ is usually not a constant, this integral is
  not separable. Yet, in some cases (like when dealing with diffuse
  lighting or, more specifically, diffuse ambient lighting), the
  information that $f(l_k,v)$ carries is low-frequency enough that it
  can reasonably be approximated as a constant with respect to the
  distribution of $L_i(l_k)$. Experimentally, it can be shown that the
  final result is not too different from what you would expect.

I'm not sure I understand the specific bit "the information that $f(l_k,v)$ carries is low-frequency enough" and also the bit "with repsect to the distribution of $L_i(l_k)$".

I've been reading something about radiometry lately. The only clue I have is that since $f(l_k,v)$ is also function of the wavelength, which is inverse w.r.t. the frequency. Therefore low frequency means high wavelength, which in turn translates some specific assumption about the colour of both incoming and reflected light I guess. Experimental results are also mentioned and I'd like to find out those if they've been published.

Is my reasoning correct? If yes, I still have the feeling I'm missing something. If not what am I missing yet? Experimental results are also mentioned, can you point out any paper so I can read about those?

Thank you.

"the information that $f(l_k,v)$ carries is low-frequency enough": As IneQuation explains, low-frequency was used to refer to the detail of the brdf function. I did actually mean that $f_r$ was low frequency though (which is the case with diffuse lighting), not $L_i$.

"with respect to $L_i$": what this means it that, since there aren't any large peaks in $f_r$, no part of $L_i$ will be wheighted much more heavily than any other in the integral.

If each function was, say, reciprocal of the other (meaning in this case that the peaks in each one would cancel out when multiplied by the other) the correct integral would be very different to the separated integral where this cancelation does not occurr. but since the two functions are not correlated (or they shouldn't be) this difference is likely to be fairly small.

Also, if you use a really high frequency BRDF that takes into acount only incident light around the normal and you separate it, you will get something that's way off. in this example the estimate using separation is 4x larger that the correct result: https://www.desmos.com/calculator/c4zgypthrc

Experimental results can be generated easily. here is an example: https://www.desmos.com/calculator/poyswojjsa it shows a low frequency brdf $f_r$ and a high frequency incident light $L_i$ function being multiplied before and after integration. the ratio between the correct result and the one obtained using separation is of about 1.6

Addendum regarding your latest comment on my other answer:

$$\frac{1}{N}\sum^{N}_{k=1}f(L_k) = \sum^{N}_{k=1}(f(L_k)\cdot \frac{1}{N})$$

Here we have "moved" or distributed $\frac{1}{N}$ into the sum. This formulation is equivalent but looks a bit more similar to an integral. if we set $\Delta L = \frac{1}{N}$ then, we can get

$$\sum^{N}_{k=1}f(L_k)\Delta L\approx \int_{L}f(L)dL$$

keeping in mind that the separation that occurs is, conceptually, separation of integrals, and that every integral "needs" its own $dL$. it follows that every montecarlo estimator needs its own $1/N$

I hope this helps you understand where the extra$\frac{1}{N}$ comes from


OpenGL: Why glTextureSubImage2D crashes? I'm trying to create single-channel (no mipmaps) texture initialized with my custom pixel data, but it crashes on glTextureSubImage2D. Any idea why?

const std::vector<GLubyte> myData(1025 * 1025, 0);
GLuint texture;
glCreateTextures(GL_TEXTURE_2D, 1, &texture);
glTextureStorage2D(texture, 1, GL_R8, 1025, 1025);
glTextureSubImage2D(texture, 0, 0, 0, 1025, 1025, GL_RED, GL_UNSIGNED_BYTE, &myData[0]);

Pixel transfer operations are governed by more than just your parameters to a pixel transfer function. They are also governed by parameters passed to the glPixelStore functions.

Of particular interest to you is GL_PACK_ALIGNMENT. This tells OpenGL what the byte alignment for each row of data is in your array. By default, this value is four, not one.

This means that OpenGL will believe that the second row of your image starts at this address: &myData[1028], since 1025 is not aligned to 4 bytes. And that means it will run off the end of your allocated buffer eventually.

If you are using data with 1-byte alignment for rows, you must inform OpenGL of this by calling glPixelStorei(GL_PACK_ALIGNMENT, 1);. Note that you should not do this by default; do it only when dealing with byte aligned data.


Confusion between usages of linear RGB and sRGB Suppose you have generated an image using linear values for RGB channels. E.g. you linearly interpolated it when doing blending, etc.. When you're going to present this on the screen, should the values be transformed as $C_{\mathrm{linear}}\to C_{\mathrm{sRGB}}$ as given here or backwards as described here?

It would seem that the image should be encoded in sRGB, thus the forward transformation is necessary. But I've measured response of my monitor, and it appears that to have linearly changing brightness with input RGB, the color data have to be transformed as $C_{\mathrm{sRGB}}\to C_{\mathrm{linear}}$ as if my image were already in sRGB color space (the behavior is consistent with that described in another question).

So is it correct that to render and present an image one has to do the following:

Draw in linear RGB to have correct blending
Treating the data as if they were sRGB, convert to linear RGB
Present on screen, so that monitor now converts this back

Then the result will be linear as expected. Or am I missing something?

The monitor interprets its input data as encoded in sRGB color space. This means that to actually display them, i.e. to get brightness of the subpixels it does the transformation $C_{\mathrm{sRGB}}\to C_{\mathrm{linear}}$. This transformation is basically application of gamma of $\approx2.2$ to the value, so that output brightness will be superlinearly dependent on input data. Since non-extremal input values are $<1$, they'll all be darker than with linear response.

Now if you have rendered an image in linear RGB color space, you need to compensate the transformation the monitor does. Thus you convert the image to sRGB, i.e. apply the $C_{\mathrm{linear}}\to C_{\mathrm{sRGB}}$ transformation. This transformation is inverse of what the monitor does: it applies a gamma of $\approx1/2.2$ to the values, so all non-extremal values will be larger than original.

In the OP I was coming to the wrong conclusion because of a mistake when trying to match the output of monitor with my input. To match it correctly one needs to either apply $C_{\mathrm{sRGB}}\to C_{\mathrm{linear}}$ transformation to input data to simulate what the monitor does with them, or apply backward transformation to the measured output — to undo monitor's processing.


Area Lights in Path Tracing I have been following the scratchapixel ray tracing tutorials online which have been great at explaining the theory behind everything & how to implement it. However, I reached the point where I want to implement Area Lights but I am clueless on how to do it. The tutorial on the website hasn't been uploaded yet (all though there is a placeholder link with a bunch of others).

There's not much online regarding implementation so this has lead me here. Could someone please explain how I can implement Area Lights with Ray Tracing as opposed to Point Lights?

I'm going to provide the simple/naive/brute force answer to this question, which does work, and gives accurate results.

There are better answers however, which make the rendering converge faster with fewer samples, but uses more advanced math.  I'm hoping someone else will provide an answer along those lines.

The simple answer is this:
Area lights are just geometry which have an emissive value.

To add support for emissive values to path tracing is pretty easy luckily.

As you know, the amount of light heading towards a pixel is the amount of light leaving the surface that a ray down that pixel would intersect.

That is calculated recursively by seeing how much light is hitting that point, integrating the BRDF, multiplying by the cosine and all that jazz.

Emissive lighting is just added to that result.  So, whatever you have now as "outgoing light" from a surface, just add to it the emissive color of that surface.  (You could also optionally do something similar to texture mapping and have emissive values vary over the object, possibly based on a texture).

That's all there is to it, it's very simple and easy.

In the rendering equation below, $L_e(\omega_o)$ represents the emissive term (that's what the e stands for, emissive), and the rest of the equation is the light integrated over the hemisphere.  You can see that you just add the emissive to the result you return as your outgoing light from a point ($L_o( \omega_o)$) to apply that emissive lighting.

$L_o( \omega_o)= L_e(\omega_o)+\int_{\Omega}{f(\omega_i, \omega_o)L_i(\omega_i)(\omega_i \cdot n)\mathrm{d}\omega_i}$

You can read more about a simple path tracer with area lights at this link, which also has commented C++ source code.

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


Points of a line given in global coordinates to calculate the pixels that are set using the midpoint algorithm without anti-aliasing I've the following problem:

Let the camera opening angle be $\frac{3}{4}\pi$ and the window be $15 \times 15$ pixels large.
Which pixels does the midpoint algorithm (without anti-aliasing) set for the line from $p_1 = (1, 1, 4)$ to $p_2 = (2, 1, 1)$, given in global coordinates?

My question are:

Why is the detail that the points are given in global coordinates important?

Could the points be given in local coordinates? If yes, how would they would look like? How would that affect the answer?

What would that mean local coordinates in this case?

I think I understood what are local coordinates, i.e. coordinates of an object given with respect to a local system to the object. I think that global coordinates at the end are just local coordinates with respect to the camera or the viewer. Am I right?
Note: I'm not asking you to solve the problem, i.e. telling me which pixels are set, maybe I'm going to ask another question regarding that.

When drawing objects that exist in 3D space to a 2D plane (like your monitor or an image), there are a number of spaces that are useful to work in:

Global Coordinates - Somewhere in your 3D world there is an origin, and objects will be placed relative to that origin. The camera will also be placed somewhere relative to that origin and will point in some direction in the world.
Object Coordinates - If you have 3D objects in your scene, they are generally created outside of the world and then transformed in some way to be placed in the world. You might create, say, a chair, with the origin of the chair where the front right leg touches the ground. In Object space, that's (0,0,0). But you'll transform that chair to where it needs to go in the scene and put it at some location in world (global) coordinates.
View Coordinates - Your virtual camera, as mentioned above, is placed somewhere and is pointing at something. Whatever your camera is pointing at will generate a space relative to the camera. You can think of the camera as always being at the origin of camera space.
Clip Space - eventually during rendering, you will use a projection matrix of some sort (perspective, orthographic, etc.) to render the final image. When running your geometry through this matrix, any coordinates that end up outside of [-1..1] will be clipped out. (In some cases [0..1] is used instead, but this detail doesn't change the point.)
Screen Space - finally, after you've clipped everything and gotten it into normalized coordinates, you then have to flatten it to 2D and display it to the user. This is screen space, and is often the same as the pixels (or some portion of them) on the user's physical display (or in an image if your rendering to an image).

So in the original problem, you've been given the requirement of having a field of view for your camera of 135° (3/4π radians), and a window that's 15x15 pixels large. So you need to construct a projection matrix with a 135° field of view and an aspect ratio of 1:1.

The view parameters aren't specified in your problem, so I guess the camera is at the world's origin, so camera/view space and world space are the same for now.

Normally, you'd take the world coordinates, (p1 and p2) and multiply them by the view matrix (identity here, so no work to do), then multiply that result by the projection matrix to get to clip space. Once you have coordinates in the normalized range, you can multiply them by the window's size (15x15 in this case) to get them into screen space and determine which pixels to set.


How to reduce the number of draw calls when rendering one object multiple times? Especially when rendering particle effects, the same object needs to be rendered several times with slightly modified properties. But these changes are often limited to properties like pose or textures, and are not in the geometry of the object itself.

A different pose (translation, rotation) is merely a matrix multiplication with the vertices, and a different texture with an applied texture atlas are merely different texture coordinates.

Nevertheless each of the particles needs to be drawn and I currently do that by calling glDrawElements for every particle, while setting uniforms appropriately. Until now this is sufficient for rendering the scene in real-time, but for more complex scenes or simply for more particles emitted by the particle system this could easily lead to dropping frame rates.

So, is there a way to reduce the number of draw calls for very similar objects?

There are many, many ways to draw things in OpenGL, so this is naturally confusing sometimes. The first method you describe, setting the shader parameters and issuing one draw call per object is usually the most inefficient, due to the high API overhead. The second one, using instanced drawing is a much smarter approach for objects with the same parameters.

When it comes to particles, specifically, there are two other methods which I'm aware of and have tested:

The first one, more traditional and easier to implement, is to generate a unique quadrilateral for each particle in the application code. Then use one of the several optimized buffer streaming paths of OpenGL to upload this data and issue a single draw call. This is the most straightforward method and provides good results. It will involve very few API calls if you can map the vertex/index buffers (glMapBuffer/MapBufferRange).
Another method is to move the whole thing to a shader program, using Transform Feedback. This method is a little more complicated to get up and running, but you can find a lot of references on the subject, such as this tutorial. This one should be the optimal path, since it moves the whole simulation to the GPU.

Those are some of the optimized ways of rendering particle effects. But OpenGL provides several other rendering paths that are better for different cases, one of such is indirect draw, which isn't available on ES at the moment, but is probably one of the fastest drawing paths available on modern PC OpenGL. Transform Feedback also requires a geometry shader, so it is not available for current OpenGL-ES.

For more on the optimized rendering paths of OpenGL, I recommend watching this very good presentation by NVidia: Approaching Zero Driver Overhead in OpenGL. At the end of the talk, they show a very interesting benchmark of several methods and how they compare with the Direct3D equivalents.


Is there a way to perform calculations of Mandelbrot set using only integer numbers? I would like to create program in JavaScript (JS) which draws Mandelbrot set with arbitrary precision (zoom). In JS there is build in integer type BigInt which support simple operations like +,*,/,power on arbitrary precision integer numbers. JS not support calculations on arbitrary precision floating point numbers.

Question: Is there a way to perform calculations of Mandelbrot set using calculations based only on integers (with arbitrary precision) - if yes how to do it?

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


Why would a 4d texture be used for color conversion? You can perform color conversion as a post process step to convert a rendered frame's colors to a different set of colors.

You might make darker colors more blue for instance, to give it a night time feel, or just modify colors to give it a different mood.

One way of doing this is using a volume texture.  You use the RGB value of a source pixel as texture coordinates into a volume texture, where you read out the RGB value that you should replace that pixel with.

This page on the 4d texture GL extension says this though:

  Four-dimensional textures are used primarily as color lookup tables
  for color conversion.

Can anyone explain why you would use a 4d texture lookup for color grading?

Alpha (transparency) lookup doesn't seem like a plausible explanation, but perhaps something like normalizing an HDR color and using the intensity as the fourth component could be?  Or perhaps time, to be able to interpolate between two different color gradings, when moving from inside to outside and exposure is changing or something?

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


Converting Vertex and Index buffers to a File Format I have many sets of vertex and index buffers that I’d like to render on to the screen. Ideally each model in different locations with-in world space. I’m actually doing a statistical study on spatial relationships between models in world space.

I have done a lot of reading on the DirectX graphics pipeline to try and get a basic understanding of how it works.

Instead of programming my own model viewer from the ground up, I thought I’d first take a look at some pre-built solutions first to save me the time. I found this: https://github.com/walbourn/directxtkmodelviewer

From what I can tell, it only supports loading models from DirectX SDK .SDKMESH files and .VBO files. So 2 questions:

1)  How do I convert my vertex and index buffers to one of those file formats?

2)  For changing the models world space in the viewer, isn’t that done by changing the view matrix before uploading the model?

Also, if this can be done in standard 3D modeling application, that would be ideal. All I want to do is view the models in relation to each other. I haven’t been able to find any applications that accept vertex and index buffers in raw form though.

Thank you.

No idea about DirectX, but in OpenGL the buffers are just a block of binary data so you can easily load and save them.

In any 3D framework, changing the object position to achieve a different relative arrangement of a scene is not done through the view matrix (that represents the camera position and viewing direction) but through the model matrix of each model. Note that the vertex and index buffers stay unaffected by matrix transformations, that is the main benefit of such pipelines.

Perhaps you could do your study with less effort using a 3D game engine, to list just a few:

https://godotengine.org/
https://www.unrealengine.com/
https://unity.com/

I am not aware that any of them would support loading the raw buffers straight away, but they can load OBJ files. OBJ is one of the oldest text-based 3D model formats. It has a simple structure, so you can easily convert your buffer data to an OBJ file.


Javascript color extraction library I have a small project that needs to extract the perceptually-salient colors of an image, in javascript. There are a few libraries out there (color-thief.js, vibrant.js). The trouble is, I'm finding that adobe color's image extraction does a noticeably better job of recognizing the little high-contrast highlight colors which make a big visual aesthetic difference, even though they represent a small number of pixels. (This is the "color mood" option). For example, the color of lipstick a person wears, or small decorations on a big cake.

I think (?) Adobe doesn't provide an API for their Color CC any more, but in any case I'd rather not rely on it. Is there a javascript library out there with similar functionality that provides access to some sort of functionality that maximizes dissimilarity between hues?

In the end I decided to use Leon Sorokin's RGBQuant.js for quantization, because it offered better flexibility over the color histogram/clustering method. I'm probably going to convert the quantized palette RGB output to HSV and implement custom clustering that picks out the "highlight" colors, favoring diversity of hue for colors with high value and saturation, and low occurrence frequency.


SVG opengl implementation - strokes overlap on transparency I am trying to implement svg rendering by opengl using meshes.
In paths with opacity less than 1, the color of overlapping strokes add up.

But it should be like this

the opengl calls are:

    glDepthMask(false);
    glEnable(GL20.GL_BLEND);
    glBlendFunc(GL20.GL_SRC_ALPHA, GL20.GL_ONE_MINUS_SRC_ALPHA);
    ....
    glDrawArrays(....);
    ....
    glDepthMask(true);

is there any opengl calls, shader or other method to have same transparent value around the whole stroke without modifying the geometry of the shape?

You can use the depth buffer to ensure that the second time you get on a point the pixels aren't written.

You start with the maximum depth value and each new path (that you want to add to the existing drawn image) you decrement the depth.


Using transformations to get from A to B image I try to solve this,

First, I used scaling: $Scale(S_x,S_y)=S(\sqrt 5,\sqrt 5)$

So, I enlarge the bottom of the house and the sides, to the same size of B image.

Later, I used rotation: $R(\theta)=R(cos^{-1}(2/\sqrt 5))$
the vertex that was $[\sqrt 5, 0, 1]$ moved $[2, 1, 1]$ I do that because I thought to translate it to $[5,3,1]$

But I stuck with that.

What is the idea of this question?

Thanks.

Just use linear algebra.

Since parallel lines have remained parallel (i.e. meaning we can probably assume the w component of homogeneous coordinates is always 1), we can just pick 3 starting points and their destinations, e.g.

(0,0,1)->(3,2,1),  (0,1,1)->(1,3,1), and  (1,0,1)->(5,3,1),

define M to be your (homogeneous) transformation matrix as

$$ M =\begin{bmatrix}
a & b & c \\
d & e & f \\
0 & 0 & 1\\
\end{bmatrix}  $$
and then we have
$$ M.\begin{bmatrix}
0 & 0 & 1 \\
0 & 1 & 0 \\
1 & 1 & 1\\
\end{bmatrix} =
\begin{bmatrix}
3 & 1 & 5 \\
2 & 3 & 3 \\
1 & 1 & 1\\
\end{bmatrix}
 $$

You just need to find the inverse of the 'starting point' matrix (which is left as an exercise for the reader ) and rearrange to compute...
$$ M=
\begin{bmatrix}
3 & 1 & 5 \\
2 & 3 & 3 \\
1 & 1 & 1\\
\end{bmatrix}
.\begin{bmatrix}
0 & 0 & 1 \\
0 & 1 & 0 \\
1 & 1 & 1\\
\end{bmatrix}^{-1}$$

To be really thorough, you could also then check that all the remaining points are transformed correctly by your result.

If you don't want to work through the maths yourself, the resulting matrix is

   $ M= \begin{bmatrix} 2 & -2 & 3 \\ 1 & 1 & 2 \\ 0 & 0 & 1\\ \end{bmatrix}$


Repeating alpha over/under operation multiple times I have two RGBA colors in linear [0..1] space, src and dst and I want to alpha blender src over/under dst n times. Of course, I could call the respective operation n times using the result from the previous one (which is what I am currently doing), but I want to optimize it into a single operation. I think just calculating the correct alpha on src should be sufficient, but I may be wrong. Here's what I currently use:

public static final Color alphaOver(Color dst, Color src) {
    float src1A = 1 - src.a;
    float outA = src.a + dst.a * src1A;

    if (outA == 0)
        return Color.TRANSPARENT;
    return new Color(
            outA,
            (src.r * src.a + dst.r * dst.a * src1A) / outA,
            (src.g * src.a + dst.g * dst.a * src1A) / outA,
            (src.b * src.a + dst.b * dst.a * src1A) / outA);
}

public static final Color alphaUnder(Color dst, Color src) {
    return alphaOver(src, dst);
}

public static final Color alphaOver(Color dst, Color src, int times) {
    Color ret = dst;
    for (int i = 0; i < times; i++)
        ret = alphaOver(ret, src);
    return ret;
}

public static final Color alphaUnder(Color dst, Color src, int times) {
    Color ret = dst;
    for (int i = 0; i < times; i++)
        ret = alphaUnder(ret, src);
    return ret;
}

How can I optimize the last two functions?

So taking source over as an example, the math works out to:

  dst = a * src + (1 - a) * dst

Taking Nathan Reed's suggestion of replacing (1 - a) with b, that gives us:

  dst = a * src + b * dst

If we perform the next iteration, it becomes:

  dst = a * src + ab * src + b^2 * dst

Doing it again, we get:

  dst = a * src + ab * src + ab^2 * src + b^3 * dst

And then:

  dst = a * src + ab * src + ab^2 * src + ab^3 * src + b^4 * dst

The pattern appears to be:

int exponent;
for (exponent = 0; exponent < numIterations; exponent++)
{
    newDst += a * pow(b, exponent) * src;
}
newDst += pow(b, exponent) * dst;


Image warping using pixel normals I want to insert a logo into an RGB image with proper 3D transformations. I can get the estimated pixel normals of this image using the code here: https://github.com/yindaz/surface_normal. Pixel normals are encoded into RGB channels as described as follows:

  The estimated normal will be saved in 16-bit PNG format, where 0-65535 in R,G,B channel correspond to [-1, 1] for the X,Y,Z component of the normal vector. We use the camera coordinates defined as - X points to the camera right, Y points to the camera forward, and Z points to the camera upward. For example, right facing wall are very red, floor are very blue, and you rarely see green as it's parallel to the camera viewing direction.

How can I transform the 2D logo to put it onto a planar surface in the image? I have tried to calculate the rotation angles of the normal to rotate it so that it has the same direction as the camera normal, and rotate the logo according to these angles, but I have no success so far.

As far as I understand from the quote above, normal of the untransformed logo should be $(x:0, y:-1, z:0)$ and let the pixel normal be $(x:a, y:b, z:c)$ with magnitude $M$. I rotate the logo around $X$ by $arccos(a/M)*180/pi$, around $Y$ by $arccos(b/M)*180/pi$, around $Z$ by $arccos(c/M)*180/pi$, and if $b < 0$ rotate around $Y$ by an additional 180 degrees.

Is there anything reasonable in these? Is it possible that I can do this using the estimated normals?

Maybe you should try to find Virtual Tangent Plane like in Surface Normal In the Wild


How to design something like this? I am very eager to know what technology might been used to implement this piece of wonder:

Any idea?

I think, it's ray tracing. What do you think?

It does look like raytracing to me as well.  It could also possibly be ray marching, which would be easier to make the tubes with since ray marching is all based on distances from things (center of tube in this case).

Refraction is what's used to bend the light that goes through the tubes.

You are then just left with how to model the tubes. One way could be to use 3d bezier curves (:

This may not have been made in real time, but may have been rendered frame by frame and then assembled into the smooth "movie" you are showing.


Does UV-mapping generate for some points of the 3D model several points in the 2D projection? Assume my 3D model is a pyramid with a rectangular base and I numerate four corners of the base with 1, 2, 3, 4 and give 5 to the top point. If I project this pyramid into a 2D star-shaped object, would not all four outer points of the star be associated with the point number 5 in the 3D model?

The way it is usually defined, every vertex in a mesh corresponds to a UV-coordinate, and there must be exactly one such coordinate per vertex, because that is how a renderer would fetch texels. So it is impossible to have multiple UV coordinates for a single vertex unless you have multiple parameterizations for your mesh.

However, meshes often get split up during parameterizations and the same vertices can end up in multiple locations, just like in your example. Most often, those vertices are simply replicated while loading/rendering so the vertex-to-uv mapping remains 1-to-1 or many-to-1. So a renderer will likely interpret your example as the following (v5 is split into 4):


qt - view geomview (.off ) file I want to open a .off file and view the content 3d model using libQGLViewer, qt-creator.

I don't know how to display the 3d model after opening its .off file.

QString fileName = QFileDialog::getOpenFileName(this,tr("Open .off model"), "/home", tr("*.off"));

From a quick gander through the examples it looks like QGLViewer is a thin wrapper around a opengl window using the legacy api.

So what you really need is a way to load a .off file to render in legacy opengl. it looks like a very simple format, even simpler than .obj. So simple that IMO it is worth it for you to figure out on your own.


Stylized math like rendering? How can I render plots and graphics that mimic how mathematicians draw diagrams?
For example Look at this shirt:

The diagrams represent 2D and 3D shapes, however they have characteristics of human drawings. The curves are not perfect and exhibit lot's of noise and waviness that a computer rendering would normally not have.
If you know any papers regarding this topic, please add them as reference to your answer.
As a side note and perhaps somewhat off topic: Is there a C++ library that can do tex rendering onto textures?

You might want to check matplotlib-cpp. You'll find a "funny-looking xkcd-styled example" in the README.
I also saw xkcd related entries on matplotlib's Python documentation. So, it should be available in Python too. However, I did not use these libraries myself. I just stumbled upon matplotlib-cpp while I was looking for a plotting library for C++.


Euler Angles - Gimbal lock, why non-orthogonal axes Already asked this question on mathexchange but didn't got any promising answers so posting it here Since I'm from a CG background.

I really can't understand how and why gimbal lock occurs using euler angles.
First of all let me clear this "I know that gimbal lock occurs when 2 gimbals/axes coincide thus losing 1 degree of freedom"

However what I'm interested in why the axes are coinciding in the first place when they are supposed to remain perpendicular when we move anyone.

According to wikipedia here :

  The general problem of decomposing a rotation into three composed
  movements about intrinsic axes was studied by P. Davenport, under the
  name "generalized Euler angles", but later these angles were named
  "Davenport angles" by M. Shuster and L. Markley.

  Davenport proved that any orientation can be achieved by composing
  three elemental rotations using non-orthogonal axes.

That's what I want to hear about more. Why "non-orthogonal" ? I can't see any problems with making the axes orthogonal? What does making the axes non-orthogonal give us?

2) After the above one has been answered. How is this all tied up to rotational matrices? How can we achieve gimbal lock through the matrices when they only rotate a given point/vector around global axis. For example If i multiply a column vector with a Rotation matrix around X-axis Rx  and then with Ry It will rotate around the global X-axis first, then global Y-axis second. So how can I achieve the lock situation using matrices?

EDIT:-
To make it more clear I've also heard about rotation orders like in the order Y-Z-X when Y axis rotates then Z and X rotate with it but when Z rotates only X rotates with it. thus making the axes non-orthogonal. I am assuming that's what it means by non-orthogonal mentioned in the wiki article. Here is a picture.

As you can see in the Y-Z-X order, the Y axis remains there as in 3rd picture causing the axes to coincide...

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

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

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

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

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

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

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


Objects in depth buffer are always black (0.0) My problem is exactly what titles and screenshot below says - rendered object depth value is 0.0 in depth buffer. I'm using DirectX 12.

Matrices:

// VS code
PixelInputType main(VertexInputType input)
{
    PixelInputType output;

    float4 position = float4(input.position, 1.0f);
    position = mul(position, g_worldMatrix);
    position = mul(position, g_viewMatrix);
    position = mul(position, g_projMatrix);
    output.position = position;

    return output;
}
// PS code
float4 main(PixelInputType input) : SV_TARGET
{
    return float4(g_texture.Sample(g_sampler, input.uv).rgb, 1.0f);
}

I am rendering single object for now, so world matrix is just identity matrix.
Regarding view/proj matrix - they're updated each frame and saved to the CB:
void Renderer::CreateViewAndPerspective()
{
    const DirectX::XMVECTOR up = DirectX::XMVectorSet(0.0f, 1.0f, 0.0f, 0.f);
    constexpr float conv{ 0.0174532925f };

    // Create the rotation matrix from the yaw, pitch, and roll values.
    const XMMATRIX rotationMatrix = XMMatrixRotationRollPitchYaw(m_cameraRotation.x * conv, m_cameraRotation.y * conv, m_cameraRotation.z * conv);

    const DirectX::XMVECTOR eye = DirectX::XMVectorSet(m_cameraPosition.x, m_cameraPosition.y, m_cameraPosition.z, 0.0f);

    //Setup target (look at object position)
    XMVECTOR target = XMVector3TransformCoord(DirectX::XMVECTOR{ 0, 0, 1, 0 }, rotationMatrix);
    target = XMVector3Normalize(target);
    target = { m_cameraPosition.x + target.m128_f32[0], m_cameraPosition.y + target.m128_f32[1], m_cameraPosition.z + target.m128_f32[2], 0.0f };

    //Create view matrix
    m_constantBuffer.value.view = DirectX::XMMatrixTranspose(DirectX::XMMatrixLookAtLH(eye, target, up));

    //Create perspective matrix
    constexpr float FOV = 3.14f / 4.0f;
    float aspectRatio = m_windowSize.x / m_windowSize.y /* 1280 x 720 */;
    m_constantBuffer.value.projection = DirectX::XMMatrixTranspose(DirectX::XMMatrixPerspectiveFovLH(FOV, aspectRatio, Z_NEAR /* 0.5f */, Z_FAR /* 100.0f */));
}

Edit: Depth buffer implementation details below (full code - https://github.com/komilll/YARE/blob/main/YARE/Renderer.cpp#L131 ):
// Create DSV
D3D12_DESCRIPTOR_HEAP_DESC desc{};
desc.NumDescriptors = 1;
desc.Type = D3D12_DESCRIPTOR_HEAP_TYPE_DSV;
desc.Flags = D3D12_DESCRIPTOR_HEAP_FLAG_NONE;
device->CreateDescriptorHeap(&desc, IID_PPV_ARGS(&heap));

// Create Depth Stencil
D3D12_DEPTH_STENCIL_VIEW_DESC depthStencilDesc = {};
depthStencilDesc.Format = DXGI_FORMAT_D32_FLOAT;
depthStencilDesc.ViewDimension = D3D12_DSV_DIMENSION_TEXTURE2D;
depthStencilDesc.Flags = D3D12_DSV_FLAG_NONE;

D3D12_CLEAR_VALUE depthOptimizedClearValue = {};
depthOptimizedClearValue.Format = DXGI_FORMAT_D32_FLOAT;
depthOptimizedClearValue.DepthStencil.Depth = 1.0f;
depthOptimizedClearValue.DepthStencil.Stencil = 0;

ThrowIfFailed(device->CreateCommittedResource(
    &CD3DX12_HEAP_PROPERTIES(D3D12_HEAP_TYPE_DEFAULT),
    D3D12_HEAP_FLAG_NONE,
    &CD3DX12_RESOURCE_DESC::Tex2D(format, width, height, 1, 0, 1, 0, D3D12_RESOURCE_FLAG_ALLOW_DEPTH_STENCIL),
    D3D12_RESOURCE_STATE_DEPTH_WRITE,
    &depthOptimizedClearValue,
    IID_PPV_ARGS(&depthStencil)
));

device->CreateDepthStencilView(depthStencil.Get(), &depthStencilDesc, dsvHeap->GetCPUDescriptorHandleForHeapStart());

/* ... */
// Usage of DSV
const float clearColor[] = { 0.0f, 0.2f, 0.4f, 1.0f };
m_commandList->ClearRenderTargetView(rtvHandle, clearColor, 0, nullptr);
m_commandList->ClearDepthStencilView(m_dsvHeap->GetCPUDescriptorHandleForHeapStart(), D3D12_CLEAR_FLAG_DEPTH, 1.0f, 0, 0, nullptr);
m_commandList->IASetPrimitiveTopology(D3D_PRIMITIVE_TOPOLOGY_TRIANGLELIST);
m_commandList->IASetVertexBuffers(0, 1, &m_modelSphere->GetVertexBufferView());
m_commandList->DrawInstanced(m_modelSphere->GetIndicesCount(), 1, 0, 0);

That's how I used to initialize viewport (CD3DX12_VIEWPORT). But I didn't realize that there are two additional fields minDepth and maxDepth. Therefore, I have min/max depth set to 0.0 and objects rendered to depth buffer are always black:
m_viewport.TopLeftY = 0.0f;
m_viewport.TopLeftX = 0.0f;
m_viewport.Width = static_cast<float>(m_windowSize.x);
m_viewport.Height = static_cast<float>(m_windowSize.y);

Fixed that by simply using constructor which default minDepth to 0.0 and maxDepth to 1.0:
m_viewport = CD3DX12_VIEWPORT{ 0.0f, 0.0f, static_cast<float>(m_windowSize.x), static_cast<float>(m_windowSize.y) };


How to get the 3D position for the point with (0,0) UV coordinates? I have a 3D triangular mesh(vertices, indices, uv coords) that I'm rendering to the screen. Let's assume that the UV mapping is one-to-one. I'm trying to find a way to find the 3D position of the point with UV coordinates equal to (0,0). I searched the internet but I only find answers that I don't find convincing.
The solution that I found:

Find, in UV space, the triangle that contains (0,0).. let's call it T
Calculate barycentric coordinates for (0,0) with respect to T
interpolate the 3D positions of T's vertices using barycentric coords to get the result.

this seems wrong to me. Here's why:
Let M be the mapping between 3D space and UV space that associates UV coords for every vertex.
Let A,B and C be the vertices of T.
Let P be the origin of UV space ( P = (0,0) ).

We have P = alpha*A + beta * B + gamma*C (alpha,beta and gamma are the barycentric coords of P with respect to T).
We assumed the UV mapping to be one-to-one, so let M° be the inverse of M.
We have :

M°(P) = M°(alpha*A + beta * B + gamma*C)

The solution in question assumes that M° is linear.. if that was the case you can have:

M°(P) = alpha*M°(A) + beta * M°(B) + gamma*M°(C)

But that is not the case (correct me if I'm wrong).

So is there a way to find the 3D position of a point with specific UV coords?

Thanks in advance.

I realize that the solution I'm giving further below is basically the one you mentioned in your question. So let's begin by addressing this:

  The solution in question assumes that M° is linear.. if that was the case you can have:

  M°(P) = alpha*M°(A) + beta * M°(B) + gamma*M°(C)

  But that is not the case (correct me if I'm wrong).

The solution in question works directly in UV space. I don't see where it would make any assumptions concerning the nature or even existence of a "mapping from 3D to UV space"!? All it does assume is that UV coordinates vary linearly across each triangle (note: each triangle, not the entire mesh)…

$$
\def\mvec#1{\begin{pmatrix}#1\end{pmatrix}}
\def\ivec#1{\left(\begin{smallmatrix}#1\end{smallmatrix}\right)}
\def\mmat#1{\begin{bmatrix}#1\end{bmatrix}}
\def\vec#1{\mathrm{\mathbf{#1}}}
\def\mat#1{\mathrm{\mathbf{#1}}}
$$

In general, the $u$ and $v$ coordinates at a point on a triangle's plane given by the affine coordinates $\lambda_1$ and $\lambda_2$ are
\begin{align}
  u &= u_1 + \lambda_1 (u_2 - u_1) + \lambda_2 (u_3 - u_1) \\
  v &= v_1 + \lambda_1 (v_2 - v_1) + \lambda_2 (v_3 - v_1)
\end{align}
where $u_i$, $v_i$ are the $u$, $v$ coordinates at vertex $i$. We can rewrite this as
\begin{equation}
  \mmat{
    u_2 - u_1 & u_3 - u_1 \\
    v_2 - v_1 & v_3 - v_1
  } \cdot \mvec{
    \lambda_1 \\
    \lambda_2
  } = \mvec{
    u - u_1 \\
    v - v_1
  }.
\end{equation}
Knowing $u_i$ and $v_i$ at each vertex of a given triangle, we can compute
\begin{equation}
  \mat M = \mmat{
    u_2 - u_1 & u_3 - u_1 \\
    v_2 - v_1 & v_3 - v_1
  }.
\end{equation}
We can then find the affine coordinates corresponding to the point with whatever $u$, $v$ coordinates we want:
\begin{equation}
  \mvec{
    \lambda_1 \\
    \lambda_2
  } =
  \mat M^{-1} \cdot \mvec{
    u - u_1 \\
    v - v_1
  }.
\end{equation}
Note that these affine coordinates always exist as long as $\mat M^{-1}$ exists (which makes sense; there will always be exactly one point somewhere in the plane of the triangle where $u$ and $v$ reach whatever value you may be looking for as long as $u$ and $v$ vary independently). However, they will not necessarily fall inside the triangle.

So one way to find the point you're looking for would be to iterate over all triangles of the mesh and, for each triangle, compute $\mat M^{-1}$, and $\lambda_1$ and $\lambda_2$ for $u = 0$ and $v = 0$, and check whether
\begin{align}
  \lambda_1 &\geq 0 & &\land & \lambda_2 &\geq 0 & &\land & \lambda_1 + \lambda_2 &\leq 1.
\end{align}
If this is the case, you have found a triangle that contains a point at which $u = 0$ and $v = 0$. You can then compute the coordinates of that point as
\begin{equation}
  \vec p = (1 - \lambda_1 - \lambda_2) \vec v_1 + \lambda_1 \vec v_2 + \lambda_2 \vec v_3
\end{equation}
where $\vec v_1$, $\vec v_2$, and $\vec v_3$ are the coordinates of the respective triangle vertices.

It may make sense to, for example, use an initial check against a given triangle's bounding rectangle in UV space and/or other methods to quickly rule out triangles that cannot contain the coordinates of interest to speed up the search…


3ds max: Z-fighting when working with Large Objects I had a z-fighting/z-clipping issue with a scene containing some very large objects yesterday in 3ds max. It can be seen in the image below:

For some reason downscaling the scene significantly (a factor of 100 approx) resolves the z-fighting issues... Please help me understand how reducing scene scale could resolve this issue and why this issue only occurs with very large scenes/objects.

From what i understand z-fighting occurs in this case because the near and far planes are very far apart compared to the minimum and maximum distance of the rendered object(s) relative to the viewport.

Please correct me if I'm wrong, my basic understanding of the z-buffer is as follows:

1) The distance of each rendered pixel of the object(s) is compared to a standard value stored in the z-buffer.

2) if the depth of the current pixel is smaller than the currently stored value for this pixel location in this image...the value for this pixel in the z-buffer is overwritten with the new and lower value.

3) before the z-depth value is stored in the z-buffer it is divided by a certain value (floating point type value), not sure why the number is devided. perhaps to reduce memory usage by storing smaller values...

I managed to fix the z-fighting in the viewports by increasing the near plane (in viewport settings) from 0.1 ish to a larger value.

Please take a look at the image below:

The red planes are the camera's near and far planes. The green planes I have drawn myself and they indicate the minimum distance of the object relative to the camera and the maximum distance of the object relative to the camera. When scaling the scene by a factor x the ratio between the variables A,B,C,D (please see the dimensions drawn in the image) should stay approximately the same (as in the ratios A/B , A/C, A/D, B/A etc. should be approximately equal before and after the scaling process takes place).
Therefore the max and min depth relative to the near and far viewport planes should be approximately equal. And as a result the depth precision of the object should not change when scaling (I think).
My understanding is clearly incorrect since the z-fighting disappears when downscaling. Why is it that scaling has any effect on the z-fighting? Does 3ds max automatically adjust the far and near plane distances relative to the viewport when scaling the object/scene?

I've read about something called depth precision (https://developer.nvidia.com/content/depth-precision-visualized). Assuming this  (objects(their pixels) that are closer to the camera being stored more accurately in the z-buffer) is true in 3ds max it makes sense that when zooming in really close to the object in the viewport...the z-fighting disappears. however, when I zoom out quite a lot in the viewport, the z-fighting also isn't visible anymore even though it should in my mind get worse since the precision of the depth info stored in the z-buffer is lower for objects further away from the camera. Am I missing something ?

PS: the viewport clipping distances are set to the default in the image below and the first image of this post.

If I have been unclear or incomplete about anything please let me know and I will try my best to explain/elaborate. Any help will be greatly appreciated. Thank you very much!

Unfortunately, I don't know how 3ds max handles its Z-buffer or its near and far planes so I can't give a definite answer. There are many ways the software could handle this.

However, you are on the right track. Z-fighting occurs when there is not enough numerical precision to distinguish two "objects" from the scene rendered at the same pixel. In such cases, the two pixels "fight" one another over which one is closest to the camera and flickering occurs.

These factors influence how a scene is susceptible to Z-fighting:

Number of bits used by the depth buffer (typically 24 or 32).
Whether the depth buffer uses float values or integers.
Whether the Z value 0 is used for the near or far plane (inverted depth buffer).
How close objects are to each other.
The size of the range between the near and far planes.

Of these, you really only have control over the near and far planes from what you say.  It looks like your model is too big for the default way 3ds max handles this so I suggest making the near plane value bigger. An increase of the near plane from 0.1 to 1.0 (*10) will give the same precision increase as dividing the far plane by 10.

For your other question, you usually get less Z-fighting when an object is closer because there is more precision in the depth buffer closer to the near plane.

Also, for your other question: Why do GPUs divide clip space Z by W, for position?

A nice Z precision visualizer: http://dev.theomader.com/depth-precision/


Getting texture value of cartesian coordinate from a Obj file I have an obj file together with an associated mtl file and texture jpg file. I have managed to find some code that sets a specific camera view and raytraces the image to output a depth map plus a point cloud. The code is using assimp and CGAL to do this. I have access to the specific face the ray intersects but I don't know how to find the pixel color value at the intersection point. Can this be done? Or is there are better way to do this using some software or library?

My final goal is to create depth and color images of different viewpoints from a set of viewpoints for given camera parameters.

Your question is very broad because it is close again to "how to I render 3D objects". So some people might give you a more precise answer than mine but at least I will try to give you some pointers.

First look on the web. There are plenty of tutorials that explain to you exactly what you to do (which is to raytrace texture mapped objects). While it's great to ask questions, there is really no excuse for posting such a generic question here, considering the amount of information that is available on the topic on the Web.
Though I will try to give you some pointers nonetheless (because we are kind and not here to judge each other;0)

What you want to do is called texture mapping. The way it's done is simple. Each triangle of your geometry holds what we call vertex coordinates, basically it explains how a triangle is actually mapped onto an image (or texture). Imagine your triangle laying flat on the image. The part that it lays on is that part of the texture that should appear on the triangle once rendered. Note that texture coordinates are 2D coordinates and that the shape of the triangle in the texture space (when it lays flat on the texture) doesn't have to match the shape of the triangle in 3D space (in which case the texture will be stretched).

So your 3D object when you render it should have a set of 3D points (the vertex position) 3D normals (the normals of the object at the vertex position) and 2D coordinates (read from the OBJ file).

When the ray hits a triangle your compute the texture coordinates at the intersection point using the barycentric coordinates of the intersection point and the texture coordinates set at the 3 vertices of the hit triangle. This is very common technique you will find many tutorials about that.

Then once you have your interpolated texture coordinates you use it to do a lookup into the texture. The texture coordinates are in the range [0:1] (in general), where (0,0) might be the top-left corner if your texture and (1,1) the bottom right. You need to remap that texture coordinate to pixel space (multiply it by image width and height and cast to integers eventually).

Now just read the pixel value at that position (remapped to pixel coordinates), and that's the color of the mapped intersection point.

Easy, breezy!


Anti-aliasing / Filtering in Ray Tracing In ray tracing / path tracing, one of the simplest way to anti-alias the image is to supersample the pixel values and average the results. IE. instead of shooting every sample through the center of the pixel, you offset the samples by some amount.

In searching around the internet, I've found two somewhat different methods to do this:

Generate samples however you want and weigh the result with a filter

One example is PBRT

Generate the samples with a distribution equal to the shape of a filter

Two examples are smallpt and Benedikt Bitterli's Tungsten Renderer

Generate and Weigh

The basic process is:

Create samples however you want (randomly, stratified, low-discrepancy sequences, etc.)
Offset the camera ray using two samples (x and y)
Render the scene with the ray
Calculate a weight using a filter function and the distance of the sample in reference to the pixel center. For example, Box Filter, Tent Filter, Gaussian Filter, etc.)

Apply the weight to the color from the render

Generate in the shape of a filter

The basic premise is to use Inverse Transform Sampling to create samples that are distributed according to the shape of a filter. For example a histogram of a samples distributed in the shape of a Gaussian would be:

This can either be done exactly, or by binning the function into a discrete pdf/cdf. smallpt uses the exact inverse cdf of a tent filter. Examples of binning method can be found here

Questions

What are the pros and cons of each method? And why would you use one over the other? I can think of a few things:

Generate and Weigh seems to be the most robust, allowing any combination of any sampling method with any filter. However, it requires you to track the weights in the ImageBuffer and then do a final resolve.

Generate in the Shape of a Filter can only support positive filter shapes (ie. no Mitchell, Catmull Rom, or Lanczos), since you can not have a negative pdf. But, as mentioned above, it's easier to implement, since you don't need to track any weights.

Though, in the end, I guess you can think of method 2 as a simplification of method 1, since it's essentially using an implicit Box Filter weight.

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

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

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

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


Are lookup textures still used for anything? I know in the not so long ago (5-10 years?) that it was popular / efficient to bake data out into textures and then read the data from the textures, often using the built in texture interpolation to get linear interpolation of the baked out data.

Now that computing time is cheaper compared to texture lookup time, this practice has definitely lessened if not all together disappeared.

My question is, are baked out textures still used for anything?  Does anyone have any usage cases for them in modern architecture?  Does it seem likely they will ever make a come back? (say, if memory technology or basic GPU architecture changes)

Yes, lookup textures are still used. For example, pre-integrated BRDFs (for ambient lighting, say), or arbitrarily complicated curves baked down to a 1D texture, or a 3D lookup texture for color grading, or a noise texture instead of a PRNG in the shader.

ALU is generally cheaper than a texture sample, true, but you still have a limited amount of ALU per frame. GPUs are good at latency hiding and small lookup textures are likely to be in the cache. If your function is complicated enough, it may still be worth using a lookup texture.


Is there any reason to prefer Direct3D over OpenGL? So I was reading this, I sort of got the reason why there are a lot more games on Microsoft windows than on any other OS. The main issue presented was that Direct3D is preferred over OpenGL.

What I don't understand is why would any developer sacrifice compatibility? That is simply a financial loss to the company. I understand that OpenGL is kind of a mess, but that should hardly be a issue for experts. Even if it is, I think that people would go a extra mile than to incur a financial loss.

Also if I'm not wrong then many cross platform applications use both Direct3D and OpenGL. I think they switch between the APIs.

This is weird as they can just use OpenGL, why even care about Direct3D?

So the question is, are there any technical issues with OpenGL or is there any support that Direct3D provides that OpenGL lacks?

I am aware that this question might be closed as being off-topic or too broad, I tried my best to narrow it down.

About the fact that there are more games for Windows, some reasons are

Windows has the majority of the market and in the past to develop cross platform games was more complicated than it is today.
DirectX comes with way better tools for developing (e.g. debugging)
Big innovations are generally first created/implemented in DirectX, and then ported to/implemented in OpenGL.
As for Windows vs Linux, you have to consider that when there is an actual standard due to marketing and historical reasons (for that see Why do game developers prefer Windows? | Software Engineering, as I said in the comments), it has its inertia.

The inertia thing is very important. If your team develop for DirectX, targeting 90% of the market (well... if you play games on pc you probaly have windows, so... 99% of the market?), why would you want to invest in OpenGL? If you already develop in OpenGL, again targeting 99% of the market, you will stick to it as long as you can. For exampe Id Tech by Id Software is an excellent game engine (powering the DOOM series) that uses OpenGL.

About the topic of your discussion, a comment.

As of today, there are many many APIs, and a common practice is to use a Game Engine that abstracts over them. For example, consider that

On most mobile platform you have to use OpenGl ES.
On PC you can use both DirectX and OpenGL
I think on XBOX you have to use DirectX.
I think that on PS you use their own API.
With old hardware you use DirectX9 or OpenGL 3, or OpenGL ES 2.
With more recent hardware you can (and want to) use DirectX 11, OpenGL 4, OpenGL ES 3.

Recently, with the advent of the new low overhead APIs - whch is a major turning point for graphical programming - there are DirectX 12 for Windows and XBOX, Metal for iOS and Vulkan (the new OpenGL) for Windows and Linux (including Android and Tizen).

There still are games that only target Windows and XBOX, but IMHO today that can only be marketing choice.


Progressive Path Tracing with Explicit Light Sampling I understood the logic behind the importance sampling for the BRDF part. However, when it comes to sampling light sources explicitly, all becomes confusing. For example, if I have one point light source in my scene and if I directly sample it at each frame constantly, should I count it as one more sample for the monte carlo integration? That is, I take one sample from cosine-weighted distribution and other from the point light. Is it two samples in total or just one? Also, should I divide the radiance coming from the direct sample to any term?

There are multiple areas in path tracing that can be importance sampled. In addition, each of those areas can also use Multiple Importance Sampling, first proposed in Veach and Guibas's 1995 paper. To better explain, let's look at a backwards path tracer:

void RenderPixel(uint x, uint y, UniformSampler *sampler) {
    Ray ray = m_scene->Camera->CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);
    SurfaceInteraction interaction;

    // Bounce the ray around the scene
    const uint maxBounces = 15;
    for (uint bounces = 0; bounces < maxBounces; ++bounces) {
        m_scene->Intersect(ray);

        // The ray missed. Return the background color
        if (ray.GeomID == INVALID_GEOMETRY_ID) {
            color += throughput * m_scene->BackgroundColor;
            break;
        }

        // Fetch the material
        Material *material = m_scene->GetMaterial(ray.GeomID);
        // The object might be emissive. If so, it will have a corresponding light
        // Otherwise, GetLight will return nullptr
        Light *light = m_scene->GetLight(ray.GeomID);

        // If we hit a light, add the emission
        if (light != nullptr) {
            color += throughput * light->Le();
        }

        interaction.Position = ray.Origin + ray.Direction * ray.TFar;
        interaction.Normal = normalize(m_scene->InterpolateNormal(ray.GeomID, ray.PrimID, ray.U, ray.V));
        interaction.OutputDirection = normalize(-ray.Direction);

        // Get the new ray direction
        // Choose the direction based on the bsdf
        material->bsdf->Sample(interaction, sampler);
        float pdf = material->bsdf->Pdf(interaction);

        // Accumulate the weight
        throughput = throughput * material->bsdf->Eval(interaction) / pdf;

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.Origin = interaction.Position;

        // Reset the other ray properties
        ray.Direction = interaction.InputDirection;
        ray.TNear = 0.001f;
        ray.TFar = infinity;

        // Russian Roulette
        if (bounces > 3) {
            float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
            if (sampler->NextFloat() > p) {
                break;
            }

            throughput *= 1 / p;
        }
    }

    m_scene->Camera->FrameBufferData.SplatPixel(x, y, color);
}

In English:

Shoot a ray through the scene
Check if we hit anything. If not we return the skybox color and break.
Check if we hit a light. If so, we add the light emission to our color accumulation
Choose a new direction for the next ray. We can do this uniformly, or importance sample based on the BRDF
Evaluate the BRDF and accumulate it. Here we have to divide by the pdf of our chosen direction, in order to follow the Monte Carlo Algorithm.
Create a new ray based on our chosen direction and where we just came from
[Optional] Use Russian Roulette to choose if we should terminate the ray
Goto 1

With this code, we only get color if the ray eventually hits a light. In addition, it doesn't support punctual light sources, since they have no area.

To fix this, we sample the lights directly at every bounce. We have to do a few small changes:

void RenderPixel(uint x, uint y, UniformSampler *sampler) {
    Ray ray = m_scene->Camera->CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);
    SurfaceInteraction interaction;

    // Bounce the ray around the scene
    const uint maxBounces = 15;
    for (uint bounces = 0; bounces < maxBounces; ++bounces) {
        m_scene->Intersect(ray);

        // The ray missed. Return the background color
        if (ray.GeomID == INVALID_GEOMETRY_ID) {
            color += throughput * m_scene->BackgroundColor;
            break;
        }

        // Fetch the material
        Material *material = m_scene->GetMaterial(ray.GeomID);
        // The object might be emissive. If so, it will have a corresponding light
        // Otherwise, GetLight will return nullptr
        Light *light = m_scene->GetLight(ray.GeomID);

        // If this is the first bounce or if we just had a specular bounce,
        // we need to add the emmisive light
        if ((bounces == 0 || (interaction.SampledLobe & BSDFLobe::Specular) != 0) && light != nullptr) {
            color += throughput * light->Le();
        }

        interaction.Position = ray.Origin + ray.Direction * ray.TFar;
        interaction.Normal = normalize(m_scene->InterpolateNormal(ray.GeomID, ray.PrimID, ray.U, ray.V));
        interaction.OutputDirection = normalize(-ray.Direction);

        // Calculate the direct lighting
        color += throughput * SampleLights(sampler, interaction, material->bsdf, light);

        // Get the new ray direction
        // Choose the direction based on the bsdf
        material->bsdf->Sample(interaction, sampler);
        float pdf = material->bsdf->Pdf(interaction);

        // Accumulate the weight
        throughput = throughput * material->bsdf->Eval(interaction) / pdf;

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.Origin = interaction.Position;

        // Reset the other ray properties
        ray.Direction = interaction.InputDirection;
        ray.TNear = 0.001f;
        ray.TFar = infinity;

        // Russian Roulette
        if (bounces > 3) {
            float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
            if (sampler->NextFloat() > p) {
                break;
            }

            throughput *= 1 / p;
        }
    }

    m_scene->Camera->FrameBufferData.SplatPixel(x, y, color);
}

First, we add "color += throughput * SampleLights(...)". I'll go into detail about SampleLights() in a bit. But, essentially, it loops through all the lights, and returns their contribution to the color, attenuated by the BSDF.

This is great, but we need to make one more change in order to make it correct; specifically, what happens when we hit a light. In the old code, we added the light's emission to the color accumulation. But now we directly sample the light every bounce, so if we added the light's emission, we would "double dip". Therefore, the correct thing to do is... nothing; we skip accumulating the light's emission.

However, there are two corner cases:

The first ray
Perfectly specular bounces (aka mirrors)

If the first ray hits the light, you should see the light's emission directly. So if we skip it, all the lights will show up as black, even though the surfaces around them are lit.

When you hit a perfectly specular surfaces you can't directly sample a light, because an input ray has only one output. Well, technically, we could check if the input ray is going to hit a light, but there's no point; the main Path Tracing loop is going to do that anyway. Therefore, if we hit a light just after we hit a specular surface, we need to accumulate the color. If we don't, lights will be black in mirrors.

Now, let's delve into SampleLights():

float3 SampleLights(UniformSampler *sampler, SurfaceInteraction interaction, BSDF *bsdf, Light *hitLight) const {
    std::size_t numLights = m_scene->NumLights();

    float3 L(0.0f);
    for (uint i = 0; i < numLights; ++i) {
        Light *light = &m_scene->Lights[i];

        // Don't let a light contribute light to itself
        if (light == hitLight) {
            continue;
        }

        L = L + EstimateDirect(light, sampler, interaction, bsdf);
    }

    return L;
}

In English:

Loop through all the lights
Skip the light if we hit it

Don't double dip

Accumulate the direct lighting from all the lights
Return the direct lighting

Finally, EstimateDirect() is just evaluating $BSDF(p, \omega_{\text{i}}, \omega_{\text{o}}) L_{\text{i}}(p, \omega_{\text{i}})$

For punctual light sources, this is simple as:

float3 EstimateDirect(Light *light, UniformSampler *sampler, SurfaceInteraction &interaction, BSDF *bsdf) const {
    // Only sample if the BRDF is non-specular
    if ((bsdf->SupportedLobes & ~BSDFLobe::Specular) != 0) {
        return float3(0.0f);
    }

    interaction.InputDirection = normalize(light->Origin - interaction.Position);
    return bsdf->Eval(interaction) * light->Li;
}

However, if we want lights to have area, we first need to sample a point on the light. Therefore, the full definition is:

float3 EstimateDirect(Light *light, UniformSampler *sampler, SurfaceInteraction &interaction, BSDF *bsdf) const {
    float3 directLighting = float3(0.0f);

    // Only sample if the BRDF is non-specular
    if ((bsdf->SupportedLobes & ~BSDFLobe::Specular) != 0) {
        float pdf;
        float3 Li = light->SampleLi(sampler, m_scene, interaction, &pdf);

        // Make sure the pdf isn't zero and the radiance isn't black
        if (pdf != 0.0f && !all(Li)) {
            directLighting += bsdf->Eval(interaction) * Li / pdf;
        }
    }

    return directLighting;
}

We can implement light->SampleLi however we want; we can choose the point uniformly, or importance sample. In either case, we divide the radiosity by the pdf of choosing the point. Again, to satisfy the requirements of Monte Carlo.

If the BRDF is highly view dependent, it may be better to choose a point based on the BRDF, instead of a random point on the light. But how do we choose? Sample based on the light, or based on the BRDF?

Why not both? Enter Multiple Importance Sampling. In short, we evaluate $BSDF(p, \omega_{\text{i}}, \omega_{\text{o}}) L_{\text{i}}(p, \omega_{\text{i}})$ multiple times, using different sampling techniques, then average them together using weights based on their pdfs. In code this is:

float3 EstimateDirect(Light *light, UniformSampler *sampler, SurfaceInteraction &interaction, BSDF *bsdf) const {
    float3 directLighting = float3(0.0f);
    float3 f;
    float lightPdf, scatteringPdf;

    // Sample lighting with multiple importance sampling
    // Only sample if the BRDF is non-specular
    if ((bsdf->SupportedLobes & ~BSDFLobe::Specular) != 0) {
        float3 Li = light->SampleLi(sampler, m_scene, interaction, &lightPdf);

        // Make sure the pdf isn't zero and the radiance isn't black
        if (lightPdf != 0.0f && !all(Li)) {
            // Calculate the brdf value
            f = bsdf->Eval(interaction);
            scatteringPdf = bsdf->Pdf(interaction);

            if (scatteringPdf != 0.0f && !all(f)) {
                float weight = PowerHeuristic(1, lightPdf, 1, scatteringPdf);
                directLighting += f * Li * weight / lightPdf;
            }
        }
    }

    // Sample brdf with multiple importance sampling
    bsdf->Sample(interaction, sampler);
    f = bsdf->Eval(interaction);
    scatteringPdf = bsdf->Pdf(interaction);
    if (scatteringPdf != 0.0f && !all(f)) {
        lightPdf = light->PdfLi(m_scene, interaction);
        if (lightPdf == 0.0f) {
            // We didn't hit anything, so ignore the brdf sample
            return directLighting;
        }

        float weight = PowerHeuristic(1, scatteringPdf, 1, lightPdf);
        float3 Li = light->Le();
        directLighting += f * Li * weight / scatteringPdf;
    }

    return directLighting;
}

In English:

First, we sample the light

This updates interaction.InputDirection
Gives us the Li for the light
And the pdf of choosing that point on the light

Check that the pdf is valid and the radiance is non-zero
Evaluate the BSDF using the sampled InputDirection
Calculate the pdf for the BSDF given the sampled InputDirection

Essentially, how likely is this sample, if we were to sample using the BSDF, instead of the light

Calculate the weight, using the light pdf and the BSDF pdf

Veach and Guibas define a couple different ways to calculate the weight. Experimentally, they found the power heuristic with a power of 2 to work the best for most cases. I refer you to the paper for more details. The implementation is below

Multiply the weight with the direct lighting calculation and divide by the light pdf. (For Monte Carlo) And add to the direct light accumulation.
Then, we sample the BRDF

This updates interaction.InputDirection

Evaluate the BRDF
Get the pdf for choosing this direction based on the BRDF
Calculate the light pdf, given the sampled InputDirection

This is the mirror of before. How likely is this direction, if we were to sample the light

If lightPdf == 0.0f, then the ray missed the light, so just return the direct lighting from the light sample.
Otherwise, calculate the weight, and add the BSDF direct lighting to the accumulation
Finally, return the accumulated direct lighting

.

inline float PowerHeuristic(uint numf, float fPdf, uint numg, float gPdf) {
    float f = numf * fPdf;
    float g = numg * gPdf;

    return (f * f) / (f * f + g * g);
}

There are a number of optimizations / improvements you can do in these functions, but I've pared them down to try to make them easier to comprehend. If you would like, I can share some of these improvements.

Only Sampling One Light

In SampleLights() we loop through all the lights, and get their contribution. For a small number of lights, this is fine, but for hundreds or thousands of lights, this gets expensive. Fortunately, we can exploit the fact that Monte Carlo Integration is a giant average. Example:

Let's define

$$h(x) = f(x) + g(x)$$

Currently, we're estimating $h(x)$ by:

$$h(x) = \frac{1}{N} \sum_{i=1}^N f(x_i) + g(x_i)$$

But, calculating both $f(x)$ and $g(x)$ is expensive, so instead we do:

$$h(x) = \frac{1}{N} \sum_{i=1}^N \frac{r(\zeta, x)}{pdf}$$

Where $\zeta$ is a uniform random variable, and $r(\zeta, x)$ is defined as:

$$r(\zeta, x) =
\begin{cases}
  f(x), & 0.0 \leq \zeta \lt 0.5 \\
  g(x), & 0.5 \leq \zeta \lt 1.0
\end{cases}$$

In this case $pdf = \frac{1}{2}$ because the pdf must integrate to 1, and there are 2 functions to choose from.

In English:

Randomly choose either $f(x)$ or $g(x)$ to evaluate.
Divide the result by $\frac{1}{2}$ (since there are two items)
Average

As N gets large, the estimate will converge to the correct solution.

We can apply this same principle to light sampling. Instead of sampling every light, we randomly pick one, and multiply the result by the number of lights (This is the same as dividing by the fractional pdf):

float3 SampleOneLight(UniformSampler *sampler, SurfaceInteraction interaction, BSDF *bsdf, Light *hitLight) const {
    std::size_t numLights = m_scene->NumLights();

    // Return black if there are no lights
    // And don't let a light contribute light to itself
    // Aka, if we hit a light
    // This is the special case where there is only 1 light
    if (numLights == 0 || numLights == 1 && hitLight != nullptr) {
        return float3(0.0f);
    }

    // Don't let a light contribute light to itself
    // Choose another one
    Light *light;
    do {
        light = m_scene->RandomOneLight(sampler);
    } while (light == hitLight);

    return numLights * EstimateDirect(light, sampler, interaction, bsdf);
}

In this code, all the lights have an equal chance of being picked. However, we can importance sample, if we like. For example, we can give larger lights a higher chance of being picked, or lights closer to the hit surface. You just have to divide the result by the pdf, which would no longer be $\frac{1}{\text{numLights}}$.

Multiple Importance Sampling the "New Ray" Direction

The current code only importance samples the "New Ray" direction based on the BSDF. What if we want to also importance sample based on the location of lights?

Taking from what we learned above, one method would be to shoot two "new" rays and weight each based on their pdfs. However, this is both computationally expensive, and hard to implement without recursion.

To overcome this, we can apply the same principles we learned by sampling only one light. That is, randomly choose one to sample, and divide by the pdf of choosing it.

// Get the new ray direction

// Randomly (uniform) choose whether to sample based on the BSDF or the Lights
float p = sampler->NextFloat();

Light *light = m_scene->RandomLight();

if (p < 0.5f) {
    // Choose the direction based on the bsdf
    material->bsdf->Sample(interaction, sampler);
    float bsdfPdf = material->bsdf->Pdf(interaction);

    float lightPdf = light->PdfLi(m_scene, interaction);
    float weight = PowerHeuristic(1, bsdfPdf, 1, lightPdf);

    // Accumulate the throughput
    throughput = throughput * weight * material->bsdf->Eval(interaction) / bsdfPdf;

} else {
    // Choose the direction based on a light
    float lightPdf;
    light->SampleLi(sampler, m_scene, interaction, &lightPdf);

    float bsdfPdf = material->bsdf->Pdf(interaction);
    float weight = PowerHeuristic(1, lightPdf, 1, bsdfPdf);

    // Accumulate the throughput
    throughput = throughput * weight * material->bsdf->Eval(interaction) / lightPdf;
}

That all said, do we really want to importance sample the "New Ray" direction based on the light? For direct lighting, the radiosity is affected by both the BSDF of the surface, and the direction of the light. But for indirect lighting, the radiosity is almost exclusively defined by the BSDF of the surface hit before. So, adding light importance sampling doesn't give us anything.

Therefore, it is common to only importance sample the "New Direction" with the BSDF, but apply Multiple Importance Sampling to the direct lighting.


Drawing Directed Lines (with arrow endings) I am working on a simple graph viewer with OpenGL. My problem with is directed edges. I tried to draw arrow endings with GL_TRIANGLES but somehow my calculations are off. I've used direction cosine vector to calculate triangle position.

Actually triangle positions are almost correct but angles are a bit off for some triangles. Here is a screenshot:

and here is my code to draw arrow endings:

foreach (Edge e in m_graph.EdgeList.Values)
{
    if (e.Source == e.Target)
    {
        //TO DO: loops
    }
    if (e.Directed)
    {
        GL.Begin(PrimitiveType.Triangles);

        Vector2 p1, p2;
        e.Source.GetAttribute("pos", out p1);
        e.Target.GetAttribute("pos", out p2);

        Vector2 direction = p2 - p1;
        direction.Normalize();

        Vector2 perpendicular = (Vector2)direction.Clone();
        perpendicular.Perpendicu();

        Vector2 c_pt = p2 - direction * VERTICLE_SIZE * 3;
        Vector2 pt1 = c_pt - perpendicular * VERTICLE_SIZE;
        Vector2 pt2 = c_pt + perpendicular * VERTICLE_SIZE;

        GL.Vertex2(p2.X, p2.Y);
        GL.Vertex2(pt1.X, pt1.Y);
        GL.Vertex2(pt2.X, pt2.Y);

        GL.End();
    }
}

also Normalize and Perpendicu methods from Vector2 class:

    /// <summary>
    /// Normalizes the vector
    /// </summary>
    public void Normalize()
    {
        m_x /= this.Length;
        m_y /= this.Length;
    }

    /// <summary>
    /// Perpedicus the vector
    /// </summary>
    public void Perpendicu()
    {
        float _x = m_x;
        m_x = -m_y;
        m_y = _x;
    }

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

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

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