Non least squares formulation to fit catmull rom spline I have a set of unordered points that I'm getting from an image attached. I'd like to simply fit a parametrized curve such as a catmull-rom curve to with n control points (n = 4 to 10, and can be assumed to be known).

I've explored ways of doing this and it seems there doesn't seem to be some library function that simply does this. I explore b-spline fitting (as I read that it's possible to convert b-spline to catmull-rom afterwards) with scipy that simply seems to unable to fit this properly.

import scipy.interpolate as interpolate
tck, u = interpolate.splprep([x, y], s=6)

N = 100
xmin, xmax = 0, 1.01
unew = np.linspace(xmin, xmax, N)
spline = interpolate.splev(unew, tck)

plt.plot(x, y, 'bo', label='Original points')
plt.plot(spline[0], spline[1], 'r', label='BSpline')
plt.show()

I'm now thinking of doing a robust regression approach now (https://scipy-cookbook.readthedocs.io/items/robust_regression.html). But, I'm lost on how to convert my set of [x, y] points to a parametric formulation for catmull-rom spline to write the loss function needed for least squares fit.

Any help is appreciated with a relevant code sample.

Result:
.

Original:

Disclaimer: I've done some work in the past with wavelet decomposition of image data but it occurred to me that it may be applicable to your problem.  Admittedly, I don't actually know if the following will work well in practice but, given how simple the approach is, and, as no other answer has been put forward, it seems to warrant investigation

Let's assume you have a sequence of 2D (or nD) vertices, $\small[P_0,P_1,P_2,...P_{N-1}]$, and you want to reduce those to a sequence of just $M$ which form your Catmull-Rom control points. For the moment let's assume $\small N=2^kM$, for some small value of $\small k$. We will apply k wavelet decompositions to reduce N to M points.

For our purposes, a single wavelet decomposition step takes a sequence of 2R points, $\small P=[P_0,P_1,...P_{2R-1}]$, and replaces it with two sequences of R points each: a low frequency set, $\small L=[L_0,...L_{R-1}]$, and a high frequency set, $\small H=[H_0,...H_{R-1}]$. I suspect in your use case you might only be interested in the low frequency sequence(s).

There are multiple different wavelets that can be used - the simplest being the Haar wavelet - but I think a linear (see Sweldens & Schröder for an easy introduction)  or perhaps a 'Catmull-Rom-based' cubic wavelet might be suitable (though I don't know if this appears in any literature).

For each wavelet decomposition step we perform two passes; a 'lifting' operation to create all the high frequency terms, i.e. either...

  Linear:  $$H_i := P_{2i+1} - \frac{(P_{2i} +P_{2i+2})}{2}$$
  CR-cubic:$$H_i := P_{2i+1} - \frac{9(P_{2i} +P_{2i+2})-(P_{2i-2} +P_{2i+4})}{16}$$

...followed by an 'update'  which filters some of the high frequency terms back into the 'even'/low frequency terms as a correction. (We'll just assume, perhaps sub-optimally, this is the same for both linear and cubic)

  $$L_i := P_{2i} + \frac{(H_{i-1} + H_{i})}{4}$$

Notes: Actually, in practice, with clever indexing, these calculations can be done 'in place' on the original array.
This has also glossed over what happens when exceeding the bounds of the arrays, e.g. at $P_{2N}$ and $L_{-1}$, but one approach is to simply repeat values and/or pad with zeros.

To do the next wavelet decomposition, you repeat the process with the low frequency data, e.g. "P:=L" and repeat.

I imagine that the final set of low frequency points, L, will form a reasonable set of control points for a piecewise Catmull-Rom curve.

More notes: In this I've assumed that all the original points were 'equally important' and probably 'reasonably evenly spread'. It could be that you may want to look at the magnitudes of some of the {H_i} terms and if some are relatively large (compared to their neighbours) perhaps take action, e.g. say, re-insert them into the lower frequency set?

As requested in the comments, I've put together a little bit of C code to do the "cubic" decomposition. For simplicity, this just uses separate arrays for each of the lower frequency results rather than, say, doing the calculation "in place" in the original array.

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

#define IS_EVEN(X) (((X)&1)==0)
#define MAX(X,Y) ((X)>=(Y)?(X):(Y))
#define MIN(X,Y) ((X)<=(Y)?(X):(Y))

typedef unsigned int UINT;

#define NUM_DIMS (2)

typedef struct
{
    float pos[NUM_DIMS];
}vec_t;

/*
// Lift
//
// Applies a wavelet "lift" step, i.e. computes high frequency terms, by
// predicting using Catmull-Rom spline through even points, and computing
// differences from odd
*/
static void Lift(const vec_t  PointsIn[],
                 const UINT NumPtsIn,
                 vec_t HFOut[])
{
    //this *example* assumes an even number of input values
    assert(IS_EVEN(NumPtsIn) &&  NumPtsIn>=2);

    /*
    // step through odd values doing a "lifting" step, i.e. predicting
    // a position from neighbourhood of even points and computing the difference
    */
    for(UINT i = 1; i < NumPtsIn; i+= 2)
    {
        /*
        // get 2 neighbours on either side of this odd point. IF we step out of range, then
        // clamp to nearest valid point
        //
        // If you want more efficient code, move these decisions out of the loop
        */
        const vec_t * pv_m3 = &PointsIn[MAX(i, 3) - 3];
        const vec_t * pv_m1 = &PointsIn[i-1];
        const vec_t * pv_p1 = &PointsIn[MIN(i+1, NumPtsIn - 1)];
        const vec_t * pv_p3 = &PointsIn[MIN(i+3, NumPtsIn - 1)];
        const vec_t * pv    = &PointsIn[i];

        vec_t *pResult = &HFOut[i/2];

        for(UINT j = 0; j < NUM_DIMS; j++)
        {
            //Predict using Catmull-Rom interpolation -
            //i.e. mid point of CR spline through the even points {p[i-3], p[i-1], p[i+1], p[i+3]}
            const float PredictedPos = (9.0/16.0) * (pv_m1->pos[j] + pv_p1->pos[j]) - (1.0/16.0) * (pv_m3->pos[j] + pv_p3->pos[j]);

            //Compute the difference from actual pos
            pResult->pos[j] =pv->pos[j] - PredictedPos;
        }
    }
}

/*
// Update
//
// Applies a wavelet "Update" step, i.e. computes the low frequency terms,
// by filtering in the 'errors' in the prediction (i.e the HF terms) into the
// 'even' vals
*/
static void Update(const vec_t  PointsIn[],
                   const UINT NumPtsIn,
                   const vec_t HFIn[],
                   vec_t LFOut[])
{

    static const vec_t TheZeroVec = {0};

    //this *example* assumes an even number of input values
    assert(IS_EVEN(NumPtsIn));

    //init as if for the "i = -1" iteration
    const vec_t * pHFToRight = &TheZeroVec;

    /*
    // step through even values doing an "update" step,
    */
    for(UINT i = 0; i < NumPtsIn; i+= 2)
    {
        /*
        // move one spot over in HF data
        */
        const vec_t *pHFToLeft = pHFToRight;
        pHFToRight = &HFIn[i/2];

        const vec_t * pv    = &PointsIn[i];

        vec_t *pResult = &LFOut[i/2];

        for(UINT j = 0; j < NUM_DIMS; j++)
        {
            //Add in 1/2 of the average of the adjacent HF/Odd terms
            pResult->pos[j] = pv->pos[j] + (pHFToLeft->pos[j] + pHFToRight->pos[j]) * 0.25;
        }
    }
}

//Some test data ...
#define NUM_ORIG_POINTS (512u)
static const vec_t PointsIn[NUM_ORIG_POINTS]=
{
#include "fjord.h"
};

#define NUM_ADDITIONAL_PASSES (4)
extern int main(int argc, const char *argv[])
{

    vec_t HighFreq[NUM_ORIG_POINTS/2];
    vec_t LowFreqs[NUM_ADDITIONAL_PASSES+1][NUM_ORIG_POINTS/2];

    /*
    // Do initial wavelet pass
    */
    Lift(PointsIn, NUM_ORIG_POINTS, HighFreq);
    Update(PointsIn, NUM_ORIG_POINTS, HighFreq, LowFreqs[0]);

    /*
    // Do a few more wavelet passes
    */
    for(UINT Pass=0; Pass < NUM_ADDITIONAL_PASSES; Pass++)
    {
        const UINT CurrentNumSrcPoints = NUM_ORIG_POINTS >> (Pass+1);

        Lift(LowFreqs[Pass],   CurrentNumSrcPoints, HighFreq);
        Update(LowFreqs[Pass], CurrentNumSrcPoints, HighFreq, LowFreqs[Pass+1]);
    }
}

The above has been run on the outline on a section of a map (a section of a fjord apparently) which I've plotted using Excel's scatter plot tools (which I hope will give a 'reasonable' approximation of Catmull-Rom splines without me having to write (or find) my own rendering code)

This shows an original set of 512 points:

..the set of 256 low-frequency points obtained from one wavelet step...

...through to the 128 points obtained by two steps:


GLSL - increasing line spacing with y I'm trying to force perspective into the horizon of a Shadertoy scene by drawing horizontal lines at decreasing distances as y grows to create an outrun-esque scene, but can't figure out which functions to use. I know the spaces need to be asymptotic up to uv.y/2, and thought I could use a combination of smoothstep and modulo like:

if (mod(smoothstep(-0.07, 2.0, uv.y/1.5), 0.03) < 0.001) {
    fragColor = vec4(1.0);
}

This does work, but the lines start wide and quickly become constant width rather than getting smaller. I'd like to offset these lines with iTime so will need to modulo that as well (eg below), but in other tests I get an effect which is nearly correct but the lines bug and either don't draw sometimes, or draw too thick.

if (mod(round(fragCoord.y+mod(iTime*150.0, 500.0)), iResolution.y/2.0-fragCoord.y) < 10.0) {
        fragColor = 1.0;
    }

Am I on the right lines?

We tackle this problem by computing, like shaders do, the color of a certain pixel with position $y$ at timestamp $t$, namely $I_{x, y}(t)$. All we want to know if that pixel must draw a bar on timestamp $t$. So simple as that, we start with

float horizon = 0.5;

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
   // Normalized pixel coordinates (from 0 to 1)
   vec2 uv = fragCoord / iResolution.xy;

   float position = uv.y;
   float time = iTime;

   // everything above our horizon we make white
   if (position > horizon - 0.01) {
      fragColor = vec4(1.0);
      return;
   }

   // check if pixel is on the bar, and if so color it white.
   float color = isOnBar(time, position);

   fragColor = vec4(color, color, color, 1.0);
}

Next step is to determine whether that pixel actually is on the bar.
This problem can be solved with the rules of gravity. We visualize the bars falling down from the horizon. For any pixel $I_{x, y}$, we can say the pixel is on the first bar if

$h_{bar}(t) = I_{y}(t)$.

We know however that bars fall with gravity from the horizon, so:

$h_{bar}(t) = h_{horizon} - gt^{2}$.

Hence, the falling time $ t_{\text{fall}} $, when the bar hits the pixel, can be calculated with

$ t_{\text{fall}} = \sqrt{ \frac{h_{horizon} - I_{y}(t_{\text{fall}})}g } $.

However, this was just our first bar. Any other time a bar hitting that pixel, is when $
(t - t_{\text{fall}}) \bmod a = 0 $, where $ a $ is our falling frequency.

e.g. first bar hits a pixel after 3 seconds, and every 2 seconds another bar falls, then for instance on the 7th second a bar again hits our pixel: $ (7 - 3) \bmod 2 $ is indeed $0$.

How much room we give when computing that modulo close to zero, is the width of our bar. Right now, we can incorporate the widening of the bars when they drop.

float gravity = 0.1;
float timeBetweenBars = 0.5;

float isOnBar(float time, float position) {
   float timeOfHit = sqrt((horizon - position) / gravity);
   float moddedTime = mod(time - timeOfHit, timeBetweenBars);

   // the lower the pixel, the more room we give to catching the bar with this pixel.
   if (moddedTime < 0.1 * (1.0 - position)) {
      return 1.0;
   };
   return 0.0;
}

Here is the code on ShaderToy


Do we really need a fence event? One of the hardest things for me with the modern closer to the metal graphics APIs is understanding synchronization. I usually try to simplify things in order to understand the concept better.

In Directx12, we can use a set of "x" fences (x being the number of back buffers) to allow the CPU to work ahead of the GPU. This could be the code we use to wait for a previous fence (https://www.braynzarsoft.net):

if (fence[frameIndex]->GetCompletedValue() < fenceValue[frameIndex])
{
    hr = fence[frameIndex]->SetEventOnCompletion(fenceValue[frameIndex], fenceEvent);
    if (FAILED(hr))
    {
        Running = false;
    }
    WaitForSingleObject(fenceEvent, INFINITE);
}

Is there a reason why we can´t just:

while(fence[frameIndex]->GetCompletedValue() < fenceValue[frameIndex])
{
}

Thanks!

In your first example, the call to WaitForSingleObject blocks the thread. This puts it into a waiting state, which takes it off the run queue and allows other threads to run on this CPU core. If no threads can run, the CPU will idle. When the fenceEvent completes, it wakes up this thread.

Your proposed replacement is called a spinlock. The thread keeps running, using 100% of a core until the fence's completed flag changes. No other code can run on this core, so there's no interleaving of CPU and GPU work. Worse yet, this code has to read from memory to check the flag every time around the loop, so it'll slow down any threads on other cores that are also accessing memory. Finally, because it's keeping the core and the memory bus active, it'll burn more power and generate more heat.

This is why fence events are needed: to enable the CPU to wait efficiently.


How to Sample 3D Points Outside and Inside the Mesh Surface I have a mesh of an object enclosed in unit volume cube and I would like to sample points inside and outside the mesh surface. What are the different ways of doing it? Is there any sample code available?

What you want is something that can perform the inside-outside test for meshes. The simplest solution would be to use Trimesh's implementation which can be ran with Embree to accelerate ray queries.
I recommend you check out recent literature on neural implicit representations if you want more alternatives. Your problem arises when one wants to learn 3D shapes using small neural networks; it is then necessary to sample points and their groundtruth attributes with respect to the mesh (signed distances or occupancies of points). Here's one example for SDFs written in CUDA. NVIDIA's Kaolin also has a CUDA kernel for this IIRC.
For more troublesome meshes, you can also use the generalized winding number, which has an implementation in libigl but I think their code is CPU-only.


How to improve on bilateral upsampling in real time scenarios It is a well known "standard" to use bilateral upscaling when it comes to comes to combine a low resolution target and an higher res one.

I have personally noticed that using the basic algorithm (with weights based on depth differences between high and low res depth values) is far from perfect in situations where high res and low res are blended, say for example an high res object inside a low res effect.

Jagged edges are still noticeable
"haloes" around edges start to appear when the objects/effects in low res occupy a smaller area on the target. By this I mean holes around the silouhette.

I am well aware that there is no way to eliminate completely these problems, but is there a way to improve bilateral upsampling to partially improve on one of the two points above?

And as side question, is there an equivalently good algorithm or variation that is cheaper?

Nearest-depth filtering is an alternative to bilaterial filtering that was specifically developed for upsampling low-resolution transparent rendering. It's a bit simpler than bilateral sampling, in that it only requires one sample from your low-resolution texture. However it can still have issues, particularly with high-frequency geometry that isn't well-represented in your low-resolution depth buffer.

At SIGGRAPH 2013, Bungie presented a new approach for handling very-low-resolution particles using what they called a "variance depth map". You can find the slides here in PDF and PPTX format.

One interesting approach that hasn't been used yet (at least, to my knowledge) would be to leverage MSAA. The idea is that you take your full-resolution depth buffer, and then alias it (or convert it, if the API doesn't allow you to do so) as a 4xMSAA depth buffer with half the width and height. You could then render your low-resolution transparents to a half-resolution render target with 4xMSAA, using your "4xMSAA" depth buffer for depth testing. Then to upscale, for each high-resolution pixel you would access the corresponding subsample from your low-resolution MSAA render target, and blend it on top. This would give you full-resolution depth testing, but half-resolution shading (with framebuffer bandwidth somewhere in between full-res and half-res due to compression). The catch here is that you need programmable MSAA sample points to do this, since you'll want to use an ordered grid sampling pattern for your low-resolution rendering instead of the typical rotated grid pattern. Newer GPU's (Maxwell from Nvidia, GCN from AMD) support doing this through vendor extensions in D3D11 and OpenGL, and it can be done natively on consoles.


Get vector length with GLM I am quite confused over how GLM library is behaving or I am using it improperly.

glm::vec2 testVec(6,-4);

float len = testVec.length();

I get the value 2 with the above code snippet. I believe I am trying to get the length of the vector defined by testVec. You know very well that it is not the correct length of the vector. What am I missing here?

Sorry folks for posting such a trivial issue! The issue is solved. I was using the wrong function. Here goes the correct one:

glm::vec2 testVec(6,-4);
float len  = glm::length(testVec);

The member function of the same name returns the number of components instead (i.e. vec2::length will always yield 2, vec3::length will always yield 3, etc.).


Is there a standard specification for 3D graphic representation? I'm reading a book on computer graphics, and at some point, it shows a 3D model, created from a .dat file.

Here are a few rules for the .dat file creation:

Made of a list of vertex positions / faces
Written counter-clockwise
Faces end with a dot (.)
Fractions are also acceptable

I want to know if this is the standard, and if not, what is the standard way of defining a 3D object.

Also, if there is a specification, who defined it and where to find it?

When talking about file formats, we are talking about persisting some data related to a 3D model/geometry. There is no universal standard on file formats for persisting 3D geometry. There are only a few formats more dominant than others.

Just as it is with image file formats, the PNGs and JPEGs are the most common formats out there today, but there is no universal agreement between applications on using one or the other. Each app uses the best fit for its purposes.

The same happens with file formats storing 3D data. Each 3D modeling software will generally have one or a set of preferred formats. Most actually define custom formats that only work with the specific version of the tool. This can happen for many reasons, from simplifying the inner workings of the application, or making loading of files faster, to binding the user to a given tool on purpose.

This .dat format you describe is a custom format created by the authors of the book, which was probably designed with simplicity in mind. It seems to be a text file similar to the Wavefront OBJ format, which in turn is a very popular format for storing static geometry, though a bit outdated by now. One could say that the .OBJ format is the .BMP of 3D model formats.

Other popular 3D model formats include:

The Collada format.
The 3DS format.
The PLY format.

And many others. More general explanation here.

Also note that I said at the beginning that these are formats for storing/persisting a 3D model or geometry on an offline storage. It doesn't mean (and it usually isn't the case) that the applications that use them will store the data internally in memory using the same layout as the file. Normally, this kind of data will undergo a lot of processing after being loaded from a file until it is, for instance, displayed on the screen.


How can I model distant rain? Close up, rain can be modeled as transparent balls of water with appropriate motion blur. This seems impractical for large volumes, which would be necessary for scenes of rain in the distance.

At distances for which the human eye cannot resolve individual raindrops, how can the effect of rain-filled atmosphere be modeled? This does not need to be efficient enough for real time use, but it does need to allow multiple frames to be produced offline for animation without jarring discontinuities or flickering.

What approaches are available for this?

When you think about rain in the distance (generally), you will have rain at several different depths from the camera, some close up, some very far away, all of which will look slightly different as you won't be able to focus on all of them. But the effect of them all layering on top of each other as they go off into the distance is what helps give the look you're aiming for, you need to recreate that depth and layering to really give the look of rain further in the distance.

Particle emitters are the general method of creating rain, as usually they use very small meshes, or even just a texture sprite of a drop. This keeps them pretty efficient - to a degree as the more particles you use and have on-screen at a time will change how heavy or resource dependant the effect becomes.

Generally you would only use large detailed meshes for rain, when you have a close up of a puddle and wanted to animate the splash, or even a close up of the single drop for some reason.

This can be combined with distance fog to further increase the sense of depth in your scene by occluding the objects in the far distance.

Again you could use another particle emitter to create dark, heavy clouds that go off into the distance, as having clouds in the scene will help to sell the look.

A lot of this will depend on what software you're using, as well as its intended platform, as different software have different solutions but particles and fog are pretty commonplace.

tldr: Particle Effects & Distance Fog


I need advice and good reference in geometric modeling? I am a student of computer science and this semester is my first semester of MSc. Next semester I must choose an advisor and the subject of my dissertation. I somewhat know what geometric modeling is and I want to research in this area, but I don't know where should I start. Can anyone who is a professional researcher in geometric modeling help me with where I should start?

Does this field have problems to work on?

As a graphics professional, not a researcher, it's a bit hard to know what exactly is or is not part of geometric modeling, but looking at this page (http://graphics.cs.ucdavis.edu/~joy/GeometricModelingLectures/) it looks like it's just all about different ways to model surfaces, whether that's with triangles, mathematical descriptions, or other things.

If that is the case, yeah, there are still open problems.

Here is a journal with some papers on the topic:
http://jcgt.org/

I would also look at the SIGGRAPH conference topics over the years, there is a lot of that sort of thing there.

IMO you are going to find that the problems are broken up into two main groups:

Real Time Graphics
Non Real Time Graphics

Video games and real time simulation software has the primary need of needing representations that can be rendered in real time, and are also fairly efficient as far as size, so that they don't take up too much memory.

Modern day movies, and also "offline processes" (such as baking out lighting into textures to use in real time applications) use non real time rendering techniques to get high quality results.  Speed and memory still matter to them of course, but the restrictions are much more relaxed as it isn't uncommon for a single frame of a movie to take a day to render, using a cluster of machines in a rendering farm.

Both real time and non real time graphics also have needs beyond rigid body geometry though, and also are interested in finding good ways to model liquids, gasses like smoke, and other things.

Besides the obvious representations of geometry like triangle and quad meshes, or bezier and other analytical shapes, there is also a lesser known concept of the signed distance field which can be used for both 3d modeling as well as 2d vector graphics.  In both cases, when you zoom into these shapes they have inifinite detail - such as being completely smooth, or having a perfectly sharp corner or whatever else.

Here's a video about one way subdivision surfaces can be achieved in real time graphics, shown at SIGGRAPH 2016, regarding the call of duty game:
https://www.youtube.com/watch?v=PTGTviYwolE

If I were in your shoes, part of what I would research would be what the "demo scene" is up to.  By demo scene I mean the people who make 64KB executables that when you run them, display photo realistic movies and have music and audio as well.  Those folks have some very creative ways of doing things that often ties into cutting edge research topics.

Shadertoy.com is a place where you can write webgl pixel shaders and share them with others.  There is a lot of really interesting and exotic graphics stuff there, including signed distance field rendering, ray marching, ray tracing, path tracing, and more.
https://www.shadertoy.com/

Here's a demoscene website:
http://www.pouet.net/

Lastly, if I may humbly point you towards a paper I wrote, I came up with a technique for getting the GPU texture sampler to calculate Bezier curve points of arbitrary dimension for you as part of a texture lookup, using the non programmable N-dimensional texture sample linear interpolator.

I believe it can be expanded to support Bezier surfaces as well, but haven't yet looked into the details of that.

http://blog.demofox.org/2016/02/22/gpu-texture-sampler-bezier-curve-evaluation/


Why the Bresenham's algorithm does not work for those two points? I am using the Bresenham algorithm from wikipedia:

 plotLine(x0,y0, x1,y1)
      dx = x1 - x0
      dy = y1 - y0
      D = 2*dy - dx
      y = y0

     for x from x0 to x1
        plot(x,y)
        if D > 0
           y = y + 1
           D = D - 2*dx
        end if
        D = D + 2*dy

I have the following two points:

P1 = (10, 120) and P2 = (18, 117)

I wrote a C program to run the algorithm and I got as output:

ePlot the following points:
10, 120
11, 120
12, 120
13, 120
14, 120
15, 120
16, 120
17, 120

Which is totally wrong, is there any other version of this algorithm which works for this scenario or it does not exist?

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

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

There are 8 cases:

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

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


Modify camera matrix so that it renders a "chunk" of the frustum Given a perspective projection matrix. How can one from it derive a set of matrices representing "chunks" of the viewport? So that the rendered chunks can be "stiched"
to form an image "chunk times higher" resolution than the viewport.

Here's an image (drawn in ms paint, ik) describing what a chunk of the frustum would look like:

I'm planning on using this for a camera mod library to take images in "super-resolution", in reality the chunks should overlap a little, to allow for exposure compensation and the blending of effects like lens flare. Most linear algebra in this project is done with a header-only library called linalg , so it'd be nice if any examples shown used that library. Thanks!

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

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

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

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

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


Getting the color palette from an SVG folder I am a discord server owner, and I want to add a few new emojis to the server. But because the color palette should match, I downloaded Twemoji as SVGs. Now I need the color palette from these emojis, and for that, I need help. I already have ways to get a palette from a folder filled with PNGs. Now I just need a way to bulk-convert SVGs to PNGs, or a way to directly extract the palette from a SVG-folder. Either way, it can't have any color blending. Neither can it have mipmapping, antialiasing or anything else that changes the color.
I appreciate any attempt to help, no matter if it is helpful or not, because this is the first time I've done something like this, and this is a learning experience for me.

You can convert SVGs to PNGs using Inkscape
inkscape in.svg --export-type=png --export-filename=out.png
You can batch convert using
inkscape *.svg --export-type=png


Approximating Geodesics in a half edge DS, how can I refine my mesh to get good approximations I implemented Djikstra's shortest path algorithm to approximate Geodesics on arbitrary meshes. Djikstra's works, but I noticed a problem inherent to the discretization of my meshes.
Consider the following figure squence:

...

This is my current refinement algorithm which is the easiest/standard face subdivision. Now consider the approximation of a geodesic in 2 points:

The blue point is where I think the actual geodesic intersects that edge, which is quite far from where the approximated geodesic passes. However that path ISN'T wrong.
Consider a square grid. The distance between any 2 points in the grid is the manhattan distance  |x| + |y|.
So as far as Djikstra's is concerned, a path that goes all the way down and then to the left has the same length as a path that goes diagonally in a staircase pattern. Refining the mesh won't change the distance either. In other words, the limit of the shortest path found by Djikstra in a regular square grid as the size of the squares goes to 0 is NOT the straigtht line connecting 2 points.
Now the actual question, does anyone know of a way to subdivide my surface that is fairly straightforward but will actually converge to the geodesic?

As you pointed out, the problem here is the discretization/subdivision of the mesh. If your mesh was made of quadrilaterals instead of triangles, the obvious subdivision strategy would be to split each quadliteral into four equally sized smaller quadliterals:
$\hspace{2cm}$

$\hspace{2cm}$

For any two points $P_1$ and $P_2$, Dijkstra's algorithm would yield sets of shortest paths between these points $P_1$ and $P_2$. The more you refine the discretization with this subdivision strategy, the more shortest paths you'll find between the two points. However, intuitively it is clear that for every subdivision level $l\in\mathbb{N}$ you could pick one of these shortest paths $p_l$ such that the sequence $(p_l)_{l\in\mathbb{N}}$ converges to the actual geodesic between $P_1$ and $P_2$ (with respect to some norm that has to be specified, for example the supremum norm).
Unfortunately, the same is not true for the standard subdivision strategy of splitting a triangle into four smaller triangles, which you have proven with your example. I believe that, at its core, the problem is that there is no way to reach the center of the triangle with a straight line from each of its edges. This can be achieved by splitting a triangle in each subdivision step into 6 smaller triangles like this:
$\hspace{2cm}$

$\hspace{2cm}$

I do not have a proof that this subdivision is more useful to compute geodesics with Dijkstra's algorithm, but it seems quite likely to me. I would be very interested to see what your results look like with this subdivision strategy! However, no matter what you do, in the end you may or may not end up with sets of shortest paths instead of a single shortest path. In this case you will need some kind of heuristic or additional algorithm to decide which path resembles the true geodesic most closely.


When unsetting a VAO, should you also unbind the associated VBOs? I'm working with a js/wasm/c++ program that uses object oriented programming when updating the OpenGL state.  There are corresponding c++ classes for VBOs (vertex buffer) and VAOs (vertex array).  When we call vao.unbind();, it simply calls glBindVertexArray(0);.  I'm seeing that unbinding the VAO does not unbind the corresponding VBOs.  One of the VAO's VBOs was still set as ARRAY_BUFFER_BINDING target.
My question is: is it best practice to automatically unset the associated VBOs when you unset the VAO?  Or should I write a separate function that unsets the VAO's VBOs?
With WebGL I'm seeing that it does not unbind the VBO when you unset the VAO, despite what this reddit post is saying.

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

glEnableVertexAttribArray
glDisableVertexAttribArray
glVertexAttribPointer
glBindBuffer(GL_ELEMENT_ARRAY_BUFFER)

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


Does a bumpmap need to be linearized from sRGB when being read? By a bumpmap, I mean the black and white height map, not the rgb normal map style one. If I were to convert a bumpmap to a normal map then I don't think I should linearize the values because if an artist creates two pixels, one looking twice as bright to them, and therefore high, let's say at 0.2 and 0.4 values in sRGB, this is how they should be read and what the normals should be based on, right?

Converting them from sRGB to linear would mean that what looked liked two pixels with a height difference of 2 times (0.2 and 0.4) would get read as (0.03 and 0.13), if I got the calculation right.

We can not answer this question. This is important, and you need to get this info from the artist (or vice versa you need to give this info to the artist before he begins) It entirely depends on what kind of assumptions were made when the artist was working with the image.

Was the artist working by numbers
Was the artist working by eye
Was the artist working with a sculpting application
Was the data built from a photograph

You would be right if the artist worked by number. But if the artist was working by eye, which is incredibly hard, then you would be a bit off. sRGB does not entirely match human brightness sense so you should really fiddle a bit with the exponent if this was the case. But you might as well skip this. For sculpting and baking the image is most likely to be linear already to begin with, although the profile might not reflect this. And for the photograph case depends heavily on the method.

Safest bet is to do no conversion at all. But it is not really defined so strictly. Kind of the entire world of image manipulation has skipped how to tag images that arent you know display images but pure data. Who knows how it is encoded.


perspective matrix derivation Now there has already been an answer here in this stackexchange, but I don't understand certain aspects of it.  The author has MORE DETAILED DESCRIPTION HERE(#1) but the problem is he doesn't describe what the linear mapping is done or how is the frustum transformed into the bounding box, there's NOWHERE IN THE WORLD mentioned how these two subspaces overlap or one is transformed to the other, for example:

Is this how these two spaces exist and after projecting the point from inside the view frustum to the near plane a linear mapping is done where the author says:

"Our input to the linear mapping is x=x_near, and our output is the x coordinate in NDC space y=x_ndc, whose range is [-1.0, 1.0]" from the above page

It looks as if he somehow connects a line from the frustum to the NDC space but not clear on how it is done e.g. does he connect the top left vertex to the top left front vertex of the NDC cube?   Are the four corners of the frustum mapped to four corners of the cube?

As much I'd like to make this answer self contained, I don't know if I can copy his materials in here, please be kind enough to check the link out and help me understand this. Thanks in advance.

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

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

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


3d Math Primer book equation derivation - Projecting One vector onto another I'm reading the book 3d Math Primer for Graphics and Game Development 1st edition by Fletcher Dunn and Ian Parberry. On page 61 there's this:

How does it derive $ v_{||} $ there? It doesn't follow from what preceded this point in the book, unless I'm missing something. Also I'm not sure, is $n$ supposed to be unit vector here? This should probably be in the errata but I didn't find anything. Anyway, how did he deduce this equation? Thanks in advance.

Since v∥ is parallel to n, it can be expressed as some multiple of n. The multiple is the ratio of the length of v∥ to the length of n. The n vector doesn’t have to be normalized, but if it is, then |n| is 1, so the equation becomes v∥ = n * |v∥|.


Do occlusion queries respect the stencil test at all? In a project where I want to select visible objects using using a screen-aligned polygon and various conditions for what "selected" means with regards to that polygon, I'm basically drawing the polygon into the stencil buffer and then drawing the individual objects with stencil testing and using occlusion queries for measuring how much the object and selection polygon overlap in screen space.

This works pretty well. But imagine my surprise when stumbling over the statement that, according to the OpenGL Wiki, occlusion queries supposedly don't actually care about the stencil test:

  They detect whether the scoped drawing commands pass the depth test and if so, how many samples in those triangles passed the Depth Test. Note that the sample count and pass/fail are based only on the depth test. The Stencil Test, or fragment shader discard is irrelevant with queries, even though they technically happen after the fragment shader...

This clearly contradicts my practical experience and made me a little worried/confused. I took a look at the specification (both the newest 4.6 core as well as good old 2.0, which we're actually targetting in our project for compatibility) and it also says that occlusion queries count the number of samples passing the depth test. It does however not make such a definite statement as the wiki pertaining to the stencil test. And since the depth test technically comes after the stencil test, I'm not entirely sure if that maybe already includes the stencil test anyway and the wording of the specification is simply using the last stage for convenience.

So, my question is, do occlusion queries respect the stencil test (or any other earlier per-fragment tests, I'm actually using scissor testing with occlusion queries in another unrelated algorithm)? I see a few possiblities here:

They do definitely respect all tests and the wording in the spec is just using the last test for convenience reasons. This would mean that fragments failing earlier tests and thus not even reaching the depth test would technically count as failing that. I'm not sure this comes too clearly out of the specification's wording, though, but maybe I'm missing something and people better-versed in the intricacies of specification laywering can shed more light on this. This would also mean the above Wiki article is slightly misleading (or I just misread it).
They definitely don't respect the stencil test and really only care about the depth test. This would mean that a fragment failing the stencil test would have to still go through depth testing. It would also mean that what I experience on my specific hardware (it's an nVidia GeForce GTX 1xxx) is clearly a bug.
Or this might simply not be a well-specified part of the API and it leaves room for interpretation. This might be rather odd, though, as this isn't a particularly rarely used feature I'd suppose.

(In practice, I could definitely restrucure my algorithm to use depth testing instead, but it is a problem that makes me curious about the theoretic implications and underpinnings nevertheless.)

The specification has been clarified on this point, as has the Wiki article. So let's explain that.

The new wording in the spec says:

  the samples-passed count is incremented for each fragment still being processed after the depth test

This is better wording, as it makes clear that the order in which things are processed matters. OpenGL defines the depth test as happening after the stencil test, scissor test, etc. So this wording clearly states that the depth test, and everything before it, matter.

So it does respect the stencil test. And fragment shader discard, since the FS is specified to come before all of these tests.

Well, unless you force these tests to happen before the fragment shader. If you do that, then discarding the fragment in the FS will not prevent that fragment from being counted by the occlusion query. So if you enforce early tests, you lose the right to use discard to affect the occlusion query count.


Understanding of Camera Up Vector Assume eye position e, gaze direction g, view-up vector t, right-handed base. As the graphs shows:

$w=-\frac{g}{||g||}$
$u=\frac{t\times w}{||t\times w||}$
$v=w\times u$

I am really confused about the 2nd equation $u=\frac{t\times w}{||t\times w||}$.
What if the camera rotate around its gaze direction? The v rotates while the t stays the same.
Then the vector t will not in the plane that w and v lies in which means the second equation is wrong.

I know I have something misunderstood.
Could someone please give me some clues?

What if the camera rotate around its gaze direction? The v rotates while the t stays the same.

The camera's orientation is defined by two directions: g and t. If the camera "rotates around its gaze direction", g will not be changing; that is what it means to "rotate around its gaze direction". If g isn't changing, and the camera's orientation is defined by two directions, then you can only get a different orientation if t is changing.

Therefore, your second statement is a contradiction of the first. Changing t is precisely how you cause the camera to "rotate around its gaze direction".

To put it a different way, there are changes to t which will not affect the values of u, v, or w. However, none of those changes will cause the camera to "rotate around its gaze direction".


Threads lockstep and conditions in compute shader I'm working with DirectCompute, but this question can be applied to general gpu programming I suppose.

As far as I know the threads in a group works in a lockstep. It means that every instruction for every thread executes at the same time, right? But what if one thread out of 1024 entered if/else condition? All other 1023 will just wait or lockstep condition will be violated?

Not all threads will execute in lockstep but they are split into groups whose threads are locked to each other.

This means that if only 1 thread out of all threads enters a branch then only 1 group will need to enter that branch while all the others will skip it.

In that group that has to execute both branches it will actually execute both branches but threads will throw away the result of the branch that it didn't need to go into.

For more information see this question.


using importance sampling to reduce the error of outscattered radiance http://renderman.pixar.com/view/implementing-a-skin-bssrdf

In this paper on subsurface scattering, I'm trying to understand how importance sampling is used to compute single scattering.

It says that the outscattered radiance for the single term, $L_o$, can be rewritten as the product of the exponential falloff in $s'_o$ with another function:

  $L_o(x_o,\vec{\omega}_o) = \bigg[
 L_i(x_i,\vec{\omega}_i)\cfrac{\sigma_s(x_o)F_p(\vec{\omega}'_i \cdotp
 \vec{\omega}'_o )}{\sigma_{tc}}e^{-s'_i \sigma_t(x_i)} \bigg]
 e^{-s_o\sigma(x_o)}$

  By picking samples according to $x \sim \sigma_t e^{-\sigma_t x}$, the
  integration in $L_o$ can be approximated by a summation by applying
  importance sampling. Also by choosing $s'_o =
 \cfrac{-\mathbb{log(random())}}{\sigma_{tr}}$, you can sum up all
  contributions without computing the falloff term in $s_o$.

I understand you can approximate an integral using importance sampling, but I didn't follow how it applies to the rendering equation. Please could you describe how this works in greater detail?

$E[f(x)]=\int f(x)p(x) = \int f(x)\cfrac{p(x)}{q(x)}q(x) \approx \cfrac{1}{n}\sum \limits_{i=1}^n  f(x_i)\cfrac{p(x_i)}{q(x_i)}$

Also, where does single scattering fit into the general rendering equation?

$L_{\text{o}}(\mathbf x,\, \omega_{\text{o}}) \,=\, L_e(\mathbf x,\, \omega_{\text{o}}) \ +\, \int_\Omega f_r(\mathbf x,\, \omega_{\text{i}},\, \omega_{\text{o}})\, L_{\text{i}}(\mathbf x,\, \omega_{\text{i}})\, (\omega_{\text{i}}\,\cdot\,\mathbf n)\, \operatorname d \omega_{\text{i}}$

In general, the idea of importance sampling is to distribute samples in a way that matches the function being integrated (or more practically, matching some factor in it). This reduces the variance of the output values, allowing the Monte Carlo integration to converge faster.

In volume scattering, when a medium is homogeneous, we have the Beer–Lambert law, which states that the transmittance falls off exponentially with depth. That's the source of factors like $e^{-\sigma x}$ in those equations. So when raymarching through a medium, the regions closer to the start of the ray generally contribute more light to the output, while the farther regions have more of their light absorbed along the way.

So, if you sample points along the ray from an exponential distribution, the samples cluster toward the closer, more important regions. Then you just have to evaluate the scattering at those points and average over them.

BTW, importance sampling in the context of surface scattering is explained pretty well by this answer to another question.

On your second question, "where does single scattering fit into the general rendering equation?", well, it doesn't! What people call the "general rendering equation" actually only includes surface scattering, not volume scattering. So it's not really that "general". A general rendering equation for volume scattering might look something like:

$$ \nabla_\omega \, L(x, \omega) = \sigma_e(x) L_e(x, \omega) - \sigma_t(x) L(x, \omega) + \sigma_s(x) \int p(\omega, \omega') \, L(x, \omega') \, d\omega' $$

This states that the derivative of radiance as you move along a ray is given by the sum of emitted light from the volume, minus extinction, plus in-scattering.

To combine surface and volume scattering, we usually just say that the surface version of the equation applies at points on surfaces, and the volume version applies everywhere else.


Was Evans & Sutherland CT5 really created in 1981? A few days ago I found a video on Twitter showing the Evans & Sutherland CT5 Flight Simulator system. It was labeled as 1981 footage and looked really amazing. The first time I saw similar 3D grahics in computer games was in the late 90s and in ealy 2000s. I wasn't suspisous because it was described as a product made for the military and I thought it's normal the technology to be more advanced.

Yesterday I showed the video to my colleagues, all of them experienced developers and computer gurus, and they said it's impossible such graphics to be achieved in the 80's. I tried to search for a proof but I found none except a short documentary on YouTube and some forum topics where people were discussing how it was achieved without mentioning anything about the exact year. I found no definitive proof that the year of creation of CT5 system was 1981. No product specifications, no publications in wellknown newspapers and not even a Wikipedia article. I'm amazed how little information is out there about something so astonishing. And this is one of the things that fuel the doubt of my colleagues.

Was Evans & Sutherland CT5 really created in 1981 or at least from the same decade?

My colleagues are with a high self esteem and it's hard to change their opinion on something. I'll appreciate any resource that proves or disproves the existance of such advanced graphics in the early 80s as well as some details how it was made.

In the link below you can see the two parts of the documentary I mentioned earlier and a short forum discussion about the CT5 and similar technologies. My colleagues don't consider it a sufficient proof.

CT5 Evans Sutherland Simulator - How did it work?

The CT5 was designed before 1981. I started working at E&S in May 1981. The CT5 was fully designed by then and 3 had been built. All 3 were still at E&S and hadn't shipped yet when I started. Two of them were on the test floor and the third was the engineering prototype that was later sold. The prototype was in the "snake room", so called because of all the black power cords hanging from the ceiling. I also remember seeing the prototype in the snake room when I made a recruiting visit to E&S in early January 1981. So, the prototype, at least, was up and running before 1981.

The first two CT5 systems were CH-46 helicopter simulators and the first one went to Cherry Point NAS. I don't think it shipped in 1981 as I remember it being there on the test floor for about a year after I joined the company. So, it probably shipped in 1982. As I recall, the two systems built after the CH-46 systems were AV8B Harrier jet simulators.


Any way to fake textureLod on GLES2? I've been having nightmares at work this week trying to get Unity's Standard shader, and my own custom car paint shader modded from it, to look good on Android.

The most obvious problem was fractured specular highlights, which I  identified as a float precision issue, and solved by changing all the directional math to use float3 instead of half3.

It now looks fine on GLES3, but on GLES2, which we want to support, the glossy reflections coming from reflection probes look terrible, with noisy matte areas all over the place.

I finally tracked this down to the cube map texture lookup itself - different blur levels are stored in different MIP levels, and a textureCubeLod lookup samples the reflection at the right smoothness level. Trouble is, GLES2 doesn't natively support lod texture lookups, so it relies on an extension. If the extension is not found, and it seems most Android devices don't have it, it falls back to a regular texture sample, with the lod parameter becoming a bias. This explains the matte areas - sharp features are causing the tex lookup to be kicked down the mip chain to the rougher maps, and the bias applied from there. From certain angles, a single bit of noise in the normal map is enough to cause the effect.

What I want to know is, does any workaround exist for this? My first thought was to 'bias the bias' using dfdx, dfdy, but unfortunately these are also GLES3 only. I may just have to kill glossy reflections altogether on older droids otherwise.

One way would be to break out the mip levels of the cube into a bunch of separate textures, then implement the mip selection and trilinear filtering yourself, by branching based on the LOD value.

That doesn't sound fast, though, given you're targeting older devices. You could speed it up by forgoing trilinear filtering and just using the nearest-neighbor mip level. Gloss values would be effectively quantized to the number of mip levels, but maybe it won't look too bad.


What is a "light source" in global illumination? In the context of global illumination and the scattering equation, do we need to understand "light sources" as surfaces with a "purely emissive" material on which no scattering is present?

In contrast, any surface in the scene geometry could be equipped with a self-emissive material and hence being somehow a "light source" (it's kind of a philosophical question, I guess).

However, since we often talk about light-carrying paths of "length" $k$, where one endpoint is located on the camera and the other on a light source, I guess that's how we need to see it. Am I correct or is there something I'm missing?

Like @lightxbulb said, a light source (i.e. an emissive surface whose L_e term in the rendering equation is greater than zero) can also have a BRDF.
People usually do not model light sources really accurately down to the Tungsten filament geometry inside of a glass-enclosed gas chamber for incandescent light bulbs. Such light sources usually are modelled as simple geometries, like spheres, capsules, lines, rectangles or more generic polygons, etc.
So, most often, light sources are usually not exactly physically accurate anyways, because it would be too computationally expensive (imagine computing the intersection of a ray with the Tungsten filament in addition to account for a very little bit of volumetric scattering through the gas).
Therefore, in the simplest cases, such an incandescent light source is down to a simple sphere anyways. In that case, it would make most sense to model this light source with a translucent BTDF (Bidirectional Transmittance Distribution Function), because most of the rays will probably not hit the filament but would go straight through the glass.
What we would have to consider, however, is Fresnel and retro-reflection on the glass surface of the light bulb.
But it all depends on how accurately you want to model your light sources.
If, on the other hand, your light source is a milk glass light bulb, then you would probably need a BRDF + BTDF (or more generally a BSDF - bidirectional scattering distribution function) to account for the light that is reflected off the glass surface and the amount of light that is trasmitting through the glass.

In any case, stopping a ray at the light source (as you pointed out) would only be actually physically correct, if the surface of the light source is a black body not reflecting any light but only emitting light due to its temperature.


What is the mathematics behind inkscape's power stroke path effect interpolator_type and interpolator_beta attributes? Inkscape has a feature called Power Stroke which is a form of path effect.  When you apply it to a stroke you gain 3 control points that allow you to adjust the thickness of the stroke.  (I have not yet figured out how to add more control points through the UI, but it can be done in the XML editor).

Inspecting the SVG reveals that the stroke gets turned into a fill with a path@d that specifies the outline, allowing the SVG file to be rendered by applications that know nothing of Power Stroke.

Applications that do understand Power Stroke will be looking for these attributes on the <path>:

   inkscape:path-effect="#path-effect881"
   inkscape:original-d="m 109.59491,132.23789 c -19.148152,28.51854 -37.80734,35.62606 -46.969458,22.13393 -9.383889,-13.81871 -39.69277,-3.57313 -46.276359,5.34215 -16.01795169,21.691 17.695947,71.79071 39.993907,62.61593 17.000232,-6.99496 23.845432,-20.30176 23.845432,-20.30176"

The definition of the PowerStroke control points occurs inside the <defs>

<inkscape:path-effect
   effect="powerstroke"
   id="path-effect881"
   is_visible="true"
   offset_points="0.11606368,4.5282059 | 0.61857805,2.0941625 | 1.992627,6.538108 | 2.7945971,4.4492524"
   sort_points="true"
   interpolator_type="CubicBezierFit"
   interpolator_beta="0.2"
   start_linecap_type="round"
   linejoin_type="bevel"
   miter_limit="4"
   end_linecap_type="butt" />

I have concluded that the @offset_points list can be interpreted as a list of t,r tuples, where t=0 is the first "major" control point in the path@d, and t=1 is the next.  A segment that is a cubic bezier curve is traditionally interpreted as having 4 control points, but for the purpose of path effects, the first and last control points are what I have termed "major" control points, and are shared with the segments on either side.

What I do not yet understand is the @interpolator_type, likewise the @interpolator_beta.  I assume this controls how you interpolate the radius you calculate between the different offset_points, but I'm not sure what exact formula to use.  When I think Bezier, I think of Bezier curves, and how they only match up with their control points at the end of the segment.  I think the _beta corresponds to the "smoothness" field in the UI.

If I wanted to calculate the Power Stroke radius at t=1.32, what formula would I use?

Not a direct answer, but I can shed light on adding and removing pink control points using the GUI:

initially PowerStrok gives you 3 pink control points.
CTRL-Click on one of them duplicates it. Drag it along the path to a new position.
CTRL-ALT-Click removes the clicked one (don't remove all)
ALT-Click on a pink control point opens a text entry for exact coordinates.


(Ray tracing) Trouble converting to screen coordinates, objects being stretched I followed along with Ray Tracing in One Weekend and managed to get the final output but I wanted to understand more about creating a camera and "painting" the screen since he didn't go over it much.
When I tried using a different way of creating a camera by spheres actually get elongated, making them look like more like ellipses. I've tried modifying the x and y assignments in screenCords but I've only managed to make more errors (such as objects wrapping around to the opposite side)
Camera.h:
#pragma once

#include "../Matrix.h"
#include "../Defs.h"
#include "Defs.h"

template<typename O>
using Point3 = Vec3<O>;

template<typename O>
using Color = Vec3<O>;

template <typename O>
class Camera{
  O Height;
  O Width;
  Vec3<O> Forward, Right, Up;
  Point3<O> Origin;

public:
  Camera(O fov, O aspect_ratio, Point3<O> origin, Point3<O> target, Vec3<O> upguide) {
    Height = atan(degrees_to_radians(fov));
    Width = Height * aspect_ratio;

    Origin = origin;

    Forward = target - origin;
    Forward.normalize();
    Right = Forward.cross(upguide);
    Right.normalize();
    Up = Right.cross(Forward);

    }

    Ray<O> get_raydir(O right, O up){
      Vec3<O> result(Forward + right * Width * Right + up * Height * Up); result.normalize();

      return Ray<O>(Origin, result);
    }

    void screenCords(O &x, O &y, O width, O height){
      x = ((2.0f * x) / width) -1.0f;
      y = ((2.0f * y) / height);
    }
};

Main.cpp
#include <iostream>
#include <cmath>
#include "../Matrix.h"
#include "Camera.h"
#include <vector>
#include "Image.h"
#include "Shapes.h"
#include "Tracer.h"
#include "../Defs.h"

template<typename O>
using Point3 = Vec3<O>;

template<typename O>
using Color = Vec3<O>;

int main(){
  const int img_ratio = 2;
  const int img_width = 640;
  const int img_height = 480;
  const int depth = 50; float t_Max = infinity; float t_Min = 0.001;

  float inv_width = 1 / float(img_width);
  float inv_height = 1 / float(img_height);

  std::vector<Sphere<float>> shapes;

  Camera<float> cam1(20.0f, img_ratio, Point3<float>(0.0f, 0.0f, 0.0f), Point3<float>(0.0f, 0.0f, -1.0f), Vec3<float>(0.0f, 1.0f, 0.0f));

  Sphere<float> cir1(0.2f, Point3<float>(0.2f, 0.0f, -1.0f));
  Sphere<float> cir2(7.0f, Point3<float>(0.0f, -7.0f, -1.0f));
  Sphere<float> cir3(0.5f, Point3<float>(1.0f, 0.0f, -1.0f));
  shapes.push_back(cir1);
  //shapes.push_back(cir2);
  //shapes.push_back(cir3);

  Tracer<float> tracer(shapes);

  std::cout << "P3\n" << img_width << ' ' << img_height << "\n255" << std::endl;

  Ray<float> ray(Point3<float>(0.0f), Vec3<float>(0.0f));

  for (int j = 0; j < img_height; j++)
  {
    std::cerr << "\rScanlines remaining: " << j << ' ' << std::flush;
    for (int i = 0; i < img_width; i++){

        float x = i;
        float y = j;

        cam1.screenCords(x, y, img_width, img_height);

        ray = cam1.get_raydir(x, y);
        //ray = Ray<float>(Vec3<float>(x1, y1, 1), Point3<float>(0.0f, 0.0f, 0.0f));
        tracer.iterator(ray, depth, t_Max, t_Min);
    }
  }
  std::cerr << "\n done " << std::endl;
}

I suspect the error is in one of these files since the spheres are actually being drawn with the colors based on normals (with the top and bottom normal colors unsurprisingly being bugged)
Here are a few examples of the output:

It looks like the parameters to your camera class are not correct
 const int img_ratio = 2;
 ...
 Camera<float> cam1(20.0f, img_ratio, Point3<float>(0.0f, 0.0f, 0.0f), Point3<float>(0.0f, 0.0f, -1.0f), Vec3<float>(0.0f, 1.0f, 0.0f));

In your code you set img_ratio to 2, I would expect it should be something like img_width/img_height (or the reverse img_height/img_width). img_ratio should also be also be a floating point to match the  type the camera class is being instantiated with.


What are Affine Transformations? What are Affine Tranformations? Do they apply just to points or to other shapes as well? What does it mean that they can be "composed"?

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

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

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

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

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

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

Which yields the same 2 equations above.

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

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

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

Identity transform

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

Scaling

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

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

Translation

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

Skew x by y

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

Skew y by x

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

Rotation

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

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

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


Is glDispatchCompute(a,b,c) any faster than glDispatchCompute(a*b*c,1,1)? In my application, there is no inherent 2D or 3D structure to the computer shader instances. So having multiple dimensions doesn't benefit me semantically in any way. I am wondering, whether it still might be beneficial to split up my workload over three dimensions to increase performance.

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


Is there a performance implication if the frame buffer resolution and the physical screen resolution are different? I was under the impression that the frame buffer resolution always has to be the same as the physical screen resolution, so changing the software resolution settings would not affect rasterization. But after I made a demo and ran a profiler against it, the profiler told me that the frame buffer resolution changed according to my settings and fewer pixels were rasterized.

This brings up other questions I am unable to figure out. Which part of the rendering pipeline does re-sample the frame buffer to match the screen resolution? And the re-sampling stage does not show up in the profiler (Intel GPA for Android), is it so fast that I can safely ignore it or I am missing something that would cost me here?

Thank you.

For PCs with discrete displays, the resampling to the screen resolution may happen either on the computer prior to the signal being sent to the display, as part of window compositing, or on the display side. I don't think it's anything to worry about in terms of FPS, in either case, but it possibly could affect latency.

It depends whether you're in "true" fullscreen (exclusive) mode or not. If the app renders to a window, including a borderless fullscreen window, then the window compositor (part of the OS) will be running each frame, and any upsampling needed can happen "for free" as part of that operation. In other words, I wouldn't expect an upsample during compositing to cost noticeably more than it does to do compositing in general (which is very cheap).

On the other hand, if you're in true fullscreen and your app has exclusive control of the display, it's possible to set the video mode to something non-native and send your rendered image directly to the display, in which case the display's scaler chip will handle it, again without performance cost to you, since it's happening in parallel to anything going on with the GPU.

You mentioned Android, so I assume you're working on a mobile device. In that case, I'm not sure if there is a scaler chip, or if the mobile GPU is just expected to always output at native resolution, in which case scaling would always be done on the GPU by the compositor software. (By the way, the reason the scaling doesn't show up in profiling is likely because it's being done by the OS rather than by your app.)

That being said, though upscaling is very cheap, it does imply an extra pass that could increase input latency. With a borderless fullscreen window, if no other windows are visible on top of your app, and no scaling is required, the compositor could be short-circuited so that your app writes directly to the display output rather than passing through the compositor. (This is an optimization that the OS would have to perform when it detects that all the conditions are right.) Similarly, on the display, if no extra image processing is required including scaling, color correction, etc. then some stages of the display hardware might be bypassed. In either case, you could see lower latency when you render directly to native resolution than when you require scaling. This is a matter of concern for VR in particular.


Is there some kind of Bresenham algorithm or equivalent for scanline rendering a rotated ellipse? Back in the day when you often had to write your own low level rendering algorithms we all used to learn the Bresenham algorithms for lines and circles.

It was almost trivially easy to extend the Bresenham circle algorithm to cover ellipses whose axes were parallel to the X and Y axes.

But it always eluded me if there was some Bresenham-style way to render an arbitrary ellipse whose axes were at some arbitrary angle.

Does such a scanline method for rotated ellipses exist?

The classic book Computer Graphics: Principles and Practice (second edition) by Foley, van Dam, et al. describes such an algorithm in section 19.2.

The explanation in the book seems to come from an MSc thesis, Raster Algorithms for 2D Primitives by Dilip Da Silva.

See also these papers:

Curve-drawing algorithms for raster displays by Van Aken and Novak (ACM TOG, 1985).
A tracking algorithm for implicitly defined curves by Chandler (IEEE Computer Graphics and Applications, 1988).


Which techniques do I need to animate a 3D model in my homemade engine Suppose I am making my own C++ 3D engine; I can already load a 3D model using assimp, display it with phong lighting, and I now want to animate it.

I'd like features in my engine to, say, load a run animation with assimp, and apply it to my human character, preferably on the GPU side. Which techniques do I need to learn to develop this?

This is a "grocery list" kind of question, where I don't necessarily need a step-by-step explanation, but more of an overview of the techniques involved to know what I need to research. (stuff like "read this OpenGL tutorial" or "learn how X works")

First decide how you are going to animate the model.

The simplest way is to just have a separate mesh per frame of animation. This requires more memory to be dedicated to mesh data more artist time to create them and it locks you into only those poses.

The next step upwards is interpolating between full meshes (keyword morph targets). This works great for facial expressions especially if you can give different weights for each vertex. It's not that great for full body animation because it's difficult to get proper circular motion into it.

The most common technique for full body animation is skeletal animation. Instead of full meshes for each frame of animation you instead have a set of model transform per "bone" and each vertex has a set of weights to decide which bones it follows. In the vertex shader you then interpolate the transforms as appropriate and then apply the interpolated transform to the vertex. This has the upside that different meshes can use the same animation data assuming the bones are set up the same.


How to match Blenders polynomial lens distortion model with GLSL? I'm trying to re-implement Blenders (libmv) polynomial lens distortion model using shadertoy, in order to apply the lens distortion in another application ideally by generating a STMap for that. I think I'm close, but when passing the same values k1=0.1 k2=0.0 k3=0.0 the resulting distortion is still too heavy in comparison to the distorted image out of Blender.

// Based on:
// https://www.shadertoy.com/view/lslGRN
// https://www.shadertoy.com/view/XdfGzH

vec3 color = vec3(0.4,0.4,0.4);
float lens = 24.0;

//note: period of 2, [-1;1]=/\
float sawtooth( float t )
{
    return abs(mod(abs(t), 2.0)-1.0);
}

vec2 BarrelDistort(vec2 p)
{
    float theta  = atan(p.y, p.x);
    float radius = length(p);
    radius = pow(radius, 1.2);
    p.x = radius * cos(theta);
    p.y = radius * sin(theta);
    return (1.0 * (p));
}

vec2 BlenderDistort(vec2 p)
{
    // https://en.wikipedia.org/wiki/Distortion_(optics)
    // https://github.com/sobotka/blender/blob/master/intern/libmv/libmv/simple_pipeline/distortion_models.h#L66
    float x = p.x;
    float y = p.y;

    float k1 = .1/2.0;
    float k2 = .0;
    float k3 = .0;

    float p1 = .0;
    float p2 = .0;

    float r2 = x*x + y*y;
    float r4 = r2 * r2;
    float r6 = r4 * r2;
    //float r_coeff = (T(1) + k1*r2 + k2*r4 + k3*r6);
    float r_coeff = (1.0 + k1*r2 + k2*r4 + k3*r6);
    //float xd = x * r_coeff + T(2)*p1*x*y + p2*(r2 + T(2)*x*x);
    float xd = x * r_coeff + 2.0*p1*x*y + p2*(r2 + 2.0*x*x);
    float yd = y * r_coeff + 2.0*p2*x*y + p1*(r2 + 2.0*y*y);

    return vec2(xd, yd);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = fragCoord.xy / iResolution.xy;
    vec2 normalizedUV = (fragCoord.xy/(iResolution.xy/2.0))-1.0;

    vec2 dist = BlenderDistort(normalizedUV);
    vec2 uv_transform = vec2(.5)+vec2(.5)*dist;

    vec3 outcol = vec3(0.0);
    //outcol += texture( iChannel0, vec2(0,1)+vec2(1,-1)*uv ).rgb;
    outcol += texture( iChannel0, uv_transform ).rgb;

    // Draw Lines
    vec2 div = 30.0 * vec2(1.0, iResolution.y / iResolution.x );
    float lines = 0.0;
    lines += smoothstep( 0.2, 0.0, sawtooth( dist.x*2.0*div.x ) );
    lines += smoothstep( 0.2, 0.0, sawtooth( dist.y*2.0*div.y ) );
    lines = clamp( lines, 0.0, 1.0 );

    //outcol *= vec3(1.0-lines); //black
    //outcol += vec3(lines); //white

    //note: force black outside valid range
    vec2 valid = step( vec2(-1.0), dist ) * step( dist, vec2(1.0) );
    outcol *= valid.x*valid.y;

    fragColor = vec4( outcol,1.0 );

}

Test image undistorted

Test image distorted in Blender k1=0.1 k2=0.0 k3=0.0

Q: What I'm missing in order to match Blender's polynomial lens distortion model?

Not sure what's wrong with your implementation, but it looks like the blender's distortion model is not proportional to the image aspect ratio, producing equal curvature in both axis..

Can you compare the input coordinates between the two implementations?


Thorough, accessible material about Fourier basis functions and Spherical Harmonics? Spherical harmonics appear in several computer graphics techniques.
I feel that in order to be a better computer graphics developer, I need to have a deep understanding of what they are and how they are used.

It seems that the reference most often recommended to understand Spherical Harmonics is "Stupid Spherical Harmonics Tricks" by Peter-Pike Sloan.

I started reading it but did not find a "satisfying" definition of SH, it seems like the document mostly relies on other references for the "basics".
Other references introduce the Fourier Basis functions as a "simpler version" of SH, but once again it seems hard to find good material explaining them.

What are thorough, accessible references to understand Fourier basis functions and Spherical Harmonics ?

wil, you largely have enough scholar background, you must have done Fourier and Laplace transforms in second year, and maybe in your engineering school again as part of signal processing classes.

If you read "stupid tricks" there is not much more you can do to find a condensed course at this point.

The second most famous paper that goes with SH for graphics is by Robin Green called "the gritty details":
http://www.research.scea.com/gdc2003/spherical-harmonic-lighting.pdf

And the third most important is the one by Ramamoorthi (the original paper preceding "An Efficient Representation for Irradiance Environment Maps"), which was called On the relationship between radiance and
irradiance: determining
the illumination from
images of a convex Lambertian object

And I think they mention somewhere that SH were previously most used by another science field, forgot which one, physics maybe, and that most of their base material came from these papers. So if you want to dig in the roots, you've got to pull out these mid last century references.


How is Iridescence handled in a Path Tracer? I have little knowledge of this topic but from what I gather in a Ray Tracer each ray is given its own wavelength as opposed to an RGB value. This automatically simulates Dispersion and Fluorescence. Is this called Spectral Rendering?

What are the approaches usually taken to render Iridescence? Could someone give a brief overview of the method and a reference to a resource I can use to learn more in-depth knowledge I can use for implementation?

Iridesence can be seen in several variations: diffraction as seen in CDs and DVDs (see Stam's SIGGRAPH 1999 paper on this), phase shifting and modulation in multi-layered surfaces as you see in oilspills on water, or structural coloration as you see in butterfly wings. (Also there is the related pearlescence). However, it's not trivial in any way, probably worth a paper or two.


How can I model ink dropped into still water? I want to be able to model an opaque liquid being dropped into a transparent liquid with sufficient velocity to cause turbulence and the resulting chaotic mixing. Assuming the two liquids are of the same density and viscosity, how could I model this?

This doesn't need to be rendered in real time. I would like the turbulent flow to gradually settle down as the two fluids mix, but it will not be required to interact with other objects.

I've tried modelling an expanding interface between the two fluids, and accelerating the vertices of the interface based on the local curvature, adding new vertices as they move apart. This seemed to work in the early stages, giving a constant exaggeration of small surface imperfections to create at least the appearance of turbulence, but it ran into problems very early which seemed to result from the interface self-intersecting. I've been trying to think of ways to redirect vertices that are near to causing a surface self-intersection, but I wonder if I would be better off starting from scratch with a more amenable method.

Simulating drops of ink is probably one of the best applications of the Vorton Method.

Unlike fluid solvers based on the velocity form of the Navier-Stokes equations, which track the density and velocity of a fluid, the Vorton Method tracks vorticity instead. It does so by simulating a large number of vortons, each carrying a small amount of vorticity; you can think of a vorton as being a tiny particle that induces a small spinning vortex around itself. The sum of velocities caused by each vortex defines a velocity field in the fluid, which is used to in turn advect the vortons themselves.

Vorton methods are difficult to combine with boundary conditions and have a few efficiency caveats, but for fluid flow where intricate turbulance is the key feature and boundaries can be avoided, the Vorton method far outshines e.g. Eulerian grid solvers, which would need a large resolution to simulate turbulence at a comparable scale.

Michael Gourlay has a very accessible article series on the vorton method. Although his articles are targeted at video games, they hold up similarly for off-line simulation. Here's also an example video of a drop of ink simulated using the vorton method.


Rendering a particle system as a smooth blob What is a good algorithm for rendering a particle system as a coherent mesh, in realtime? I am running a smoothed-particle hydrodynamics fluid simulation, and I would like to render the particles not as balls but as a body of water.
I understand that the basis is to find all places where the combined density of the particle system is equal to a certain boundary value, which will yield a surface which is smooth. However, how do you actually do so quickly and translate it into vertices and triangles to render?
Another option is raytracing, raymarching specifically. Just march until the density reaches the boundary value. However, I'm still stuck on how to get the surface normal in this case, as well as how to make it performant. Surely you don't need to iterate through all particles, the far away particles won't have any effect, but I don't know how to spatially cull particles.

What are some algorithms for spatially culling far away particles with no local contribution?
What are some algorithms for getting the surface normal of a body of water made of particles?
What are some algorithms for getting the surface of a body of water made of particles, in the form of a mesh?

For culling faraway particles: as lightxbulb said, the keyword to search for is "acceleration structure". There are a variety of options: BVH, octree, kd-tree, uniform grid, or some combination of these. Which one will be most suitable is going to depend on the details of your situation—how large/dense is the simulation, how much you care about performance vs memory vs implementation complexity, and so on.
They all work by subdividing space into volumes in some fashion and keeping track of which particles are in each volume. Then, when you're marching you only need to keep track of which volume you're in, and that will tell you which particles can potentially contribute to the density at that point. For instance if it's a uniform grid, you only need to look at the current cell and possibly some of its immediate neighbors if they're within the maximum particle radius of your current point.
Getting a smoothed mesh with surface normals etc. is usually done by defining a density field for the particle system. This is a function on 3D space that consists of summing up a falloff function for each particle over a defined radius:
$$
D(\vec r) = \sum_{i \in \text{particles}} f \left(\frac{|\vec r - \vec r_i|}{R_i} \right)
$$
where $f(u)$ is some function that falls off smoothly from $f(0) = 1$ to $f(1) = 0$ (for example, a simple cubic falloff: $f(u) = 1 - 3u^2 + 2u^3$), and $\vec r_i, R_i$ are the position and radius of the $i$th particle. Then you set some threshold density value and let the surface be defined by all the points where $D(\vec r) = \text{threshold}$.
For getting the surface normal: the surface normal is just the (normalized) gradient vector of the density field. This can be calculated by finite differences between nearby points in the field:
$$
\begin{aligned}
\nabla D(\vec r) &= \left( \frac{\partial D}{\partial x}, \frac{\partial D}{\partial y}, \frac{\partial D}{\partial z} \right) \\
&\approx \frac{1}{\epsilon} \left( D(\vec r + \epsilon \vec x) - D(\vec r), D(\vec r + \epsilon \vec y) - D(\vec r), D(\vec r + \epsilon \vec z) - D(\vec r) \right)
\end{aligned}
$$
or analytically by adding up gradient contributions from each particle.
$$
\begin{aligned}
\nabla D(\vec r) &= \sum_i \nabla \left( f \left(\frac{|\vec r - \vec r_i|}{R_i} \right) \right) \\
&= \sum_i \text{normalize}(\vec r - \vec r_i) \frac{1}{R_i} \frac{df}{du}
\end{aligned}
$$
For getting the surface in the form of a mesh: the keyword to search for is "isosurface extraction". Marching cubes is one example of such an algorithm, as lightxbulb mentioned. It evaluates the density field at evenly spaced grid points, uses those values to estimate where the intersection of the surface with the cubes' edges will be, and constructs from that a tiny mesh within each cube. Stitching all those meshes together yields a representation of the full surface. Naturally, the accuracy of it depends on the size of the grid used.
There is also a variant that uses tetahedra instead of cubes, and another similar algorithm called "dual contouring". Here's a nice article showing how these work interactively (in 2D). There's also a great deal of literature on how to improve accuracy, performance, etc. with these kinds of methods.


Texture coordinates are not updating according to the vertex divisor specified I want to draw a string of text with a single draw call (Instancing). I am reading a bitmap for each character texture. I want to pass the texture coordinates for each vertex of a quad (a single character). Now these texture coordinates are supposed to be updated for each vertex processed by the vertex shader because I set the glVertexAttribDivisor() to 0 which means it should be. But it's not happening.

Here's the source code where I am sending the data to the GPU...

//TEXTURE COORDS
glBindBuffer(GL_ARRAY_BUFFER, mVBOTexCoords);
glBufferData(GL_ARRAY_BUFFER, sizeof(glm::vec2)*mTexCoords.size(), mTexCoords.data(), GL_STATIC_DRAW);

glEnableVertexAttribArray(3);
glVertexAttribPointer(3, 2, GL_FLOAT, GL_FALSE, 0, 0);

glVertexAttribDivisor(3,0);

and the vertex shader...

#version 400

out vec2 tex_coords;

layout(location = 0) in vec2 pos;
layout(location = 1) in vec2 offset;
layout(location = 2) in vec2 scale;
layout(location = 3) in vec2 tex;

uniform mat4 projection;

void main()
{
    tex_coords = tex;

 vec4 position = vec4(pos,0.0f,1.0f);
     position.xy*=scale;
     position.xy+=offset;

     gl_Position = projection*position;

}

Sorry for the mess. I hope you get it. I am not getting why its not updating the texture coordninates for each vertex.

Now this texture coordinates supposed to be updated for each vertex processed by vertex shader because i set the glVertexAttribDivisor() to 0 which means it should be.

No, it doesn't.

A divisor of 0 means that the texture coordinate work exactly as if you hadn't called glVertexAttribDivisor for that attribute at all. 0 is not a special case; it's the default.

That means that the texture coordinate will be updated at the same rate as the position. So it will get the same value for each instance, exactly like it does for the positions.

There is no way to use instancing to get an attribute to always increase. It either is a per-instance attribute (changing more-or-less based on gl_InstanceID) or it is a per-vertex attribute (changing based on gl_VertexID). To get the per-instance behavior, you use a divisor 1 or greater. To get the per-vertex behavior, you use a divisor of 0.


Vulkan: how does host coherence work? I have a memory that is persistently mapped. It is both HOST_VISIBLE and HOST_COHERENT but not DEVICE_LOCAL. (My physical device offers no memory that is exclusively HOST_VISIBLE without also being HOST_COHERENT.)

This memory contains vertex buffers. I send commands every frame to draw those buffers. Some of those, I change the content every frame, and others more rarely.

I believe buffers somehow need to be transfered by the driver into the GPU before it can draw them, yes? Here, I never need to invalidate/flush anything because the memory is HOST_COHERENT. How does the driver know what to update/transfer and when? Does it have to a way to detect parts where I change stuff, or does it compulsively retransfer data every frame?

How does the driver know what to update/transfer and when?

It doesn't. There is no transfer or update in that fashion.

All memory allocated through a Vulkan device represents memory that some device operations can read directly. If a Vulkan device advertises that it can use visible/coherent memory as source vertex data for rendering operations, that does not mean that some kind of upload is going on to make this happen. The device is reading that data from the memory you write to. The vkFlushMappedRange function merely expunges writes from CPU cache lines, so that they actually reach the memory.

Memory that is host-coherent and host-visible would represent memory that either:

Is not cached by the CPU, and therefore the GPU sees the memory as it is. This also means that CPU writes and reads tend to be slow.
The GPU can directly access the CPU's caches, to one degree or another.

Usually, #2 happens when the GPU is on-die with the CPU. Which would also mean that this would be device-local memory. Since you've said that it isn't device-local, odds are good that it is simply not cached. It may use write-combining (which improves the efficiency of sequential writes), but CPU writes otherwise are delivered as they are. Given that the memory does not seem to include "cached" in its description, I'd say #1 is the most likely.

Note that just because you don't need to explicitly flush memory does not change the fact that you still need a dependency between the CPU operation writing that data and the GPU command(s) that read it. This would generally be some from of event with a memory dependency.


How to understand Z-Fighting? I searched for Z Buffer on Google and I discovered that's a 2 dimensional array used by the graphical card and whose aim is to enable the handling of two pixels' depth in order to display one of these pixels in front of the other.

To illustrate my words, let an image which has only 1 pixel.

What is a Z-Buffer? What is its utility?

Introduction

Let $P_1$ be the pixel we want to draw in front of $P_2$ (thus, $P_1$ will hide $P_2$ at drawing). So $P_1$ and $P_2$ have the same $X$ and $Y$ coordinates.

We decide that the $Z$ coordinate of $P_1$ (ie. its depth) is less than $P_2$'s one (it means $P_1$ is closer to the computer's user than $P_2$): $Z(P_1) = 1$, and $Z(P_2) = 2$.

Our Z-Buffer has, by definition, the same size as our image: 1 case. Each case of the Z-Buffer contains one pixel's depth and each case is associated with one image's pixel.

First, we initialize the (only one) case of our Z Buffer at the maximum value we can store in its cases (depending on the Z Buffer's type, on the OS, etc.: double/float/integer).

$P_1$

Now, we want to first draw $P_1$. The Z Buffer is used: the depth of the $(X(P_1); Y(P_1))$ case's value is compared to the $P_1$'s depth. Of course, $Z(P_1)$ is lesser than $Z( (X(P_1) ; Y(P_1)) )$. The results of this comparison are: first, we replace the old case value with the new $Z$ value (which is: $Z(P_1)$ ); then, $P_1$ is drawn.

$P_2$

Finally, we want to draw $P_2$. Again the Z Buffer is used. $(X(P_2); Y(P_2))$ case's value is compared to $Z(P_2)$ (note that: $X(P_2) = X(P_1)$ and $Y(P_2) = Y(P_1)$ - cf. "So $P_1$ and $P_2$ have the same $X$ and $Y$ coordinates.").
The results of this comparison are: $(Z(P_2) = 2)$ is superior to $( Z( (X(P_2); Y(P_2)) ) = 1)$; first, we don't replace any case; then, $P_2$ is not drawn.

That's how a Z Buffer is used to display a pixel in front of another.

Example of application

http://raphaello.univ-fcomte.fr/IG/ZBuffer/Images/ZBuffer-03.gif

Z-Fighting

What is Z-Fighting?

If two polygons are very closed AND (note the "AND") if the camera is moving, Z-Fighting phenomenon appears: https://icculus.org/neverball/mapping/zfighting.png

How to reproduce Z-Fighting (ie. how does it work? why does it appear?)

I know Z Fighting is caused by the resolution in depth of the Z Buffer, due to some conversion/casts mathematical and precision errors. But I don't understand it...

In particular, I don't understand why the camera must move.

My question is: could you explain this please?

The camera does not need to move for this problem to exist. You can see the mixed polygons as in your linked image even with a static camera.

Things are worse with a moving camera because it makes the position of the polygons change, which can lead to different rounding. This means it may not be the same polygon which comes up in front at a given location for every frame. Even a very small change in the camera could swap which polygon is visible. So you get flickering which is visually a lot worse because it draws attention.

One possible solution is to progressively blend surfaces together when their Z values become too similar. However, I'm not certain how easily that can be implemented in traditional Z-buffer rendering.


Calculating Skeleton-Deformed Mesh Normals I'm implementing mesh skinning in a project of mine. I can calculate vertex positions based on bone transformations, but I've run into a problem with calculating their normals.

From the resources I've read, it appears that a given vertex's skinned normal is the weighted sum of each influencing bone's transform applied to the vertex's normal.

What I don't understand, is how this would produce equivalent normals to what you would get if you calculated each vertex's normal directly. (I.E. from the sum of the normals of its adjacent faces)

For example, consider the case of three connected vertices with the normals shown below. Note the bone in the center, which has the same position as the center vertex.

Next, suppose this bone is turned 45 degrees CCW. If this bone only influenced the vertex to its right, my understanding is that it would produce the following result:

This is because there is only one influencing bone, and the normal is rotated by the same rotation as the bone, 45 degrees CCW.

Next, consider that the bone influenced the middle and right vertices. Each vertex's normal would be rotated exactly as in the previous image.

But, neither of these would produce the "correct" normal for the middle vertex, shown here:

So, my question is, is this the intended behavior for this algorithm? Is it only meant to be an approximation, or have I missed something critical? How are these normals meant to be calculated?

For what it's worth, the resources I've been following are presented here and here.

You're right, standard skinning pipelines don't produce correct normals. You get similar problems when using translation on joints, which is very common in facial animation. The only way to get a correct result is to take the average of the surrounding face normals - this is not information you will be able to get in the vertex shader, so you'll have to use a geometry shader or run your skinning on CPU or a compute shader.