Relative coordinates -- cumulative sum

Question

I have a vertex shader happily producing all the vertices I want for a 2D plot. Now I want to also render a plot of the same data, but each point relative to the previous one, like SVG Paths using l instead of L.

I have tried a geometry shader, but the primitive would have to be a line strip with all the vertices -- 500,000 or more.

I'm pretty sure a vertex shader can't use the output of the previous instance, and computing the sum of 1..n is difficult to parallelise.

The only method I can think of is;

1. Use a compute shader to calculate the data
2. Calculate cumulative sums with CPU
3. Send data to vertex shader that does very little

Is there a way of doing this with a vertex, geometry or mesh shader? Or with an extension that might help?

It seems to be something basic enough, and useful enough to be included in SVG.

What is the best method to tackle this?

Answer 1

Can a GPU do this? Sure. Can it do it fast? Not really, no.

The nature of the algorithm is antagonistic to how a GPU wants to do work. GPUs gain performance through massive parallelism. Your algorithm is inherently single-threaded: no answer can be computed until the prior answer is computed.

Will the processing time of such a single-threaded operation be greater than the cost of a full GPU->CPU->GPU synchronization? That depends on how much data you're talking about. That being said, unless the data is being computed by the GPU, the CPU could just hold onto a copy of it and compute the values itself.

Answer 2

I am rewriting this answer because of misunderstandings and a better explanation of the algorithm.

The render pipeline is not well suited for this kind of task. It is better to use compute shaders instead.

The following algorithm is designed to be run as a compute shader.

First, let's start with a simple example:

Suppose we have 50000 values stored in a shader storage buffer object (SSBO1). And let us further assume that all values are 1. So figure 1 shows the input data.

Figure 1: 50000 data points, each has the value 1.

This algorithm calculates the cumulative sum. For our simple example, the result looks like Figure 2.

Figure 2: The final result of the simple example (cumulative sum / discrete integral).

Compute shaders can have different workgroup sizes. In our case, we try to keep the workgroup size as small as possible, with respect to the minimum warp/wave length of 32 (Nvidia) 64 (AMD). Let's assume we are using an AMD graphics card, then our workgroup size is (x = 64, y = 1, z = 1). The algorithm needs additional SSBOs to store intermediate values. The additional SSBOs will contain the following number of values:

unsigned int workgroupSize = 64;
unsigned int dataCount = 50000;

unsigned int SSBO2_count = static_cast<unsigned int>(ceil(static_cast<float>(inputDataCount) / workgroupSize)); // 782
unsigned int SSBO3_count = static_cast<unsigned int>(ceil(static_cast<float>(SSBO2_count) / workgroupSize)); //13

The first compute shader will work as follows: Each invocation will load the global value from SSBO1 (original data) and store the value in shared memory so that other invocations within that workgroup can quickly access that value. Then a barrier is needed to ensure that all data is loaded correctly. Then the first invocation within each workgroup accumulates these values while updating the original SSBO1, and stores the accumulated workgroup result to SSBO2.

#version 430 core
layout (local_size_x = 64, local_size_y = 1, local_size_z = 1) in;

layout(std430, binding = 1) buffer SSBOInput
{
    float values[]; //first execution: 50000 values. Second execution: 782 values
}ssboInput;

layout(std430, binding = 2) buffer SSBOOutput
{
    float values[]; //first execution: 782 values. Second execution 13 values
}ssboOutput;    

shared float loadedValues[64];
uniform int countData;
uniform int writeToOutput;

void main()
{
    loadedValues[gl_LocalInvocationIndex] = 0.0;
    if(gl_GlobalInvocationID.x < countData) //check if inside ssbo1
    {
        loadedValues[gl_LocalInvocationIndex] = ssboInput.values[gl_GlobalInvocationID.x];
    }
    barrier(); //ensure, that all values are loaded to shared memory
    if(gl_LocalInvocationIndex == 0)
    {
        float sum = 0;
        for(int i = 0; i < 64; ++i)
        {
            sum += loadedValues[i];
            if(gl_GlobalInvocationID.x + i < countData)
            {
                ssboInput.values[gl_GlobalInvocationID.x + i] = sum;
            }
        }
        if(writeToOutput == 1)
            ssboOutput.values[gl_WorkGroupID.x] = sum;
    }
}

Bind SSBO1 to binding location 1 (ssboInput) and SSBO2 to binding location 2 (ssboOutput). The uniforms countData = 50000 and writeToOutput = 1. After running 782 workgroups, the values will look like Figure 3.

Figure 3: These are the values stored in SSBO1 after the compute shader is executed.

In our example, the SSBO2 (782 data points) will all store the value 64. This is the case because the input values are all 1. If we ask for the cumulative value of index (e.g. 200), we only need to calculate the following:

$$indexMax = floor(200 / 64)-1$$ $$CumulativeSum = SSBO1[200] + \sum_{i = 0}^{indexMax}(SSBO2[i])$$

If we look at the equation, we will see that the problem is still the same, except that the difference is between 50000 and 782 data points.

So we do the same thing again. First, bind SSBO2 to binding slot 1. And bind SSBO3 to binding slot 2. Change the variable countData to 782. writeToOutput = 1. And run the compute shader with 13 workgroups.

Let's take a look at the SSBO value:. SSBO1 is not touched and has the values from figure 3. SSBO2 has changed! The values are accumulated depending on the workgroup size (see Figure 4). SSBO3 is now filled, the values are all 4096.

[]

Figure 4: These are the values stored in SSBO2 after the second compute shader call.

So we do the same thing again for the last time: This time it is the last execution of this compute shader. We bind SSBO3 to binding slot 1 and disable writeToOutput = 0. countData = 13. We only need to execute one workgroup.

Final step!

The precalculations are finished. Now we can compute the cumulative sum of the individual values within $log_{64}(50000)$ steps. The following compute shader stores the cumulative sum in SSBO1.

#version 430 core
layout (local_size_x = 64, local_size_y = 1, local_size_z = 1) in;

layout(std430, binding = 1) buffer SSBO1
{
    float values[]; //50000 values
}ssbo1;

layout(std430, binding = 2) buffer SSBO2
{
    float values[]; //782 values
}ssbo2;   

layout(std430, binding = 3) buffer SSBO3
{
    float values[]; //13 values
}ssbo3;    

uniform int countData;

void main()
{
    if(gl_GlobalInvocationID.x < countData) //check if inside ssbo1
    {
        float sum1 = ssbo1.values[gl_GlobalInvocationID.x];
        int index2 = gl_GlobalInvocationID.x / 64 - 1;
        float sum2 = index2 >= 0 ? ssbo2.values[index2] : 0;
        int index3 = index2 / 64 - 1;
        float sum3 = index3 >= 0 ? ssbo3.values[index3] : 0;

        //store cululative sum
        ssbo1.values[gl_GlobalInvocationID.x] = sum1 + sum2 + sum3;
    }
}

The cumulative sum is stored in SSBO1 (see Figure 2).

This code has not been tested yet. So there might be some sintax errors in it.... Sorry for that.