Реализация простого Z-буфера в CUDA

У меня есть трехмерное облако точек, и я проецирую пиксели на плоскость изображения. Поскольку некоторые 3D-точки отображаются на один и тот же пиксель, я хочу, чтобы для моей камеры использовался только пиксель с самым низким значением Z. Я использую Z-буфер — массив с плавающей точкой — чтобы отслеживать значения глубины. Вот некоторый псевдокод:

// Initialize z-Buffer with max depth (99999.9f)
// Go through every point of point cloud...
// Project Point (x,y,z) to image plane (u,v)
int newIndex = v*imgWidth+u;
float oldDepth = zbuffer[newIndex];

if (z < oldDepth){
zbuffer[newIndex] = z; // put z value in buffer
outputImg[newIndex] = pointColor[i]; // put pixel in resulting image
}

У меня есть отлично работающая одноядерная версия процессора.

Версия CUDA выглядит хорошо и очень быстро, но только области, где z-тест играет роль я думаю, что они очень «полосатые», что означает, что некоторые точки фона перезаписывают пиксели переднего плана. Кроме того, когда я смотрю на цветное изображение, я вижу случайные цветные полосы с цветами, которых нет на изображении.

Версия CUDA выглядит примерно так:

//Initialize, kernel, project, new coordinates...

const float oldDepth = outputDepth[v * outputMaxWidth + u];

if (z < oldDepth){
outputDepth[v * outputMaxWidth + u] = z;
const int inputColorIndex = yIndex * inputImageStep + 3*xIndex;
const int outputColorIndex = yIndex * outputImageStep + 3*xIndex;
outputImage[outputColorIndex] = inputImage[inputColorIndex]; //B
outputImage[outputColorIndex + 1] = inputImage[inputColorIndex + 1]; //G
outputImage[outputColorIndex + 2] = inputImage[inputColorIndex + 2]; //R
}

Я думаю, что параллелизм является проблемой здесь. Один поток может записать ближайшее Z-значение этого пикселя в Z-буфер, но в то же время другой поток считывает старое значение и перезаписывает правильное значение.

Как я могу предотвратить это в CUDA?

Edit1:
Уменьшение размера блока с (16,16) до (1,1) приведет к появлению менее полосатых рисунков, но это будет выглядеть как дыры в 1 пиксель.

Edit2:
Вот минимальный пример:

#include "cuda_runtime.h"#include "device_launch_parameters.h"
#include <stdio.h>

cudaError_t insertToZBuffer(int *z, const int *a, unsigned int size);

__global__ void zbufferKernel(int *z, const int *a)
{
int i = threadIdx.x;
if (a[i] < z[0]){
z[0] = a[i]; //  all mapped to pixel index 0
}
}

int main(){
for (int i = 0; i < 20; ++i){
const int arraySize = 5;
const int a[arraySize] = { 1, 7, 3, 40, 5 }; // some depth values which get mapped all to index 0
int z[arraySize] = { 999 }; // large depth value

insertToZBuffer(z, a, arraySize);

printf("{%d,%d,%d,%d,%d}\n", z[0], z[1], z[2], z[3], z[4]);
cudaDeviceReset();

}
return 0;
}

cudaError_t insertToZBuffer(int *z, const int *a, unsigned int size){
int *dev_a = 0;
int *dev_z = 0;
cudaError_t cudaStatus;
cudaStatus = cudaSetDevice(0);
cudaStatus = cudaMalloc((void**)&dev_z, size * sizeof(int));
cudaStatus = cudaMalloc((void**)&dev_a, size * sizeof(int));
cudaStatus = cudaMemcpy(dev_a, a, size * sizeof(int), cudaMemcpyHostToDevice);
cudaStatus = cudaMemcpy(dev_z, z, size * sizeof(int), cudaMemcpyHostToDevice);
zbufferKernel<<<1, size >>>(dev_z, dev_a);
cudaStatus = cudaGetLastError();
cudaStatus = cudaDeviceSynchronize();
cudaStatus = cudaMemcpy(z, dev_z, size * sizeof(int), cudaMemcpyDeviceToHost);

cudaFree(dev_z);
cudaFree(dev_a);

return cudaStatus;
}

Значение от z в индексе 0 должно быть 1, потому что это самое низкое значение, но это 5, которое является последним значением a.