Sobol sampler (xor scrambled)

HummaWhite · HummaWhite · commit 4789c1325334 · 2022-10-03T15:48:47.000-04:00
diff --git a/sobol_10k_200.bin b/sobol_10k_200.bin
diff --git a/src/common.h b/src/common.h
@@ -1,10 +1,7 @@
 #pragma once
 #include <iostream>
 
-#define MESH_DATA_STRUCT_OF_ARRAY false
-#define MESH_DATA_INDEXED false
-
-#define DEV_SCENE_PASS_BY_CONST_MEM false
+#define SAMPLER_USE_SOBOL true
 
 #define SCENE_LIGHT_SINGLE_SIDED true
 
diff --git a/src/material.h b/src/material.h
@@ -67,7 +67,7 @@ __device__ inline float smithG(float cosWo, float cosWi, float alpha) {
     return schlickG(glm::abs(cosWo), alpha) * schlickG(glm::abs(cosWi), alpha);
 }
 
-__device__ static float ggxDistrib(float cosTheta, float alpha) {
+__device__ static float GTR2Distrib(float cosTheta, float alpha) {
     if (cosTheta < 1e-6f) {
         return 0.f;
     }
@@ -78,18 +78,19 @@ __device__ static float ggxDistrib(float cosTheta, float alpha) {
     return nom / denom;
 }
 
-__device__ static float ggxPdf(glm::vec3 n, glm::vec3 m, glm::vec3 wo, float alpha) {
-    return ggxDistrib(glm::dot(n, m), alpha) * schlickG(glm::dot(n, wo), alpha) *
+__device__ static float GTR2Pdf(glm::vec3 n, glm::vec3 m, glm::vec3 wo, float alpha) {
+    return GTR2Distrib(glm::dot(n, m), alpha) * schlickG(glm::dot(n, wo), alpha) *
         Math::absDot(m, wo) / Math::absDot(n, wo);
 }
 
 /**
-* Sample GGX microfacet distribution, but only visible normals.
+* Sampling (Trowbridge-Reitz/GTR2/GGX) microfacet distribution, but only visible normals.
 * This reduces invalid samples and make pdf values at grazing angles more stable
 * See [Sampling the GGX Distribution of Visible Normals, Eric Heitz, JCGT 2018]:
 * https://jcgt.org/published/0007/04/01/
+* Note: 
 */
-__device__ static glm::vec3 ggxSample(glm::vec3 n, glm::vec3 wo, float alpha, glm::vec2 r) {
+__device__ static glm::vec3 GTR2Sample(glm::vec3 n, glm::vec3 wo, float alpha, glm::vec2 r) {
     glm::mat3 transMat = Math::localRefMatrix(n);
     glm::mat3 transInv = glm::inverse(transMat);
 
@@ -175,9 +176,9 @@ struct Material {
             return glm::vec3(0.f);
         }
 
-        glm::vec3 f = fresnelSchlick(glm::dot(h, wo), baseColor * metallic);
+        glm::vec3 f = fresnelSchlick(glm::dot(h, wo), glm::mix(glm::vec3(.08f), baseColor, metallic));
         float g = smithG(cosO, cosI, alpha);
-        float d = ggxDistrib(glm::dot(n, h), alpha);
+        float d = GTR2Distrib(glm::dot(n, h), alpha);
 
         return glm::mix(baseColor * PiInv * (1.f - metallic), glm::vec3(g * d / (4.f * cosI * cosO)), f);
     }
@@ -186,7 +187,7 @@ struct Material {
         glm::vec3 h = glm::normalize(wo + wi);
         return glm::mix(
             Math::satDot(n, wi) * PiInv,
-            ggxPdf(n, h, wo, roughness * roughness) / (4.f * Math::absDot(h, wo)),
+            GTR2Pdf(n, h, wo, roughness * roughness) / (4.f * Math::absDot(h, wo)),
             1.f / (2.f - metallic)
         );
     }
@@ -198,7 +199,7 @@ struct Material {
             sample.dir = Math::sampleHemisphereCosine(n, r.x, r.y);
         }
         else {
-            glm::vec3 h = ggxSample(n, wo, alpha, glm::vec2(r));
+            glm::vec3 h = GTR2Sample(n, wo, alpha, glm::vec2(r));
             sample.dir = -glm::reflect(wo, h);
         }
 
diff --git a/src/pathtrace.cu b/src/pathtrace.cu
@@ -141,15 +141,17 @@ __device__ Ray sampleCamera(const Camera& cam, int x, int y, glm::vec4 r) {
 * motion blur - jitter rays "in time"
 * lens effect - jitter ray origin positions based on a lens
 */
-__global__ void generateRayFromCamera(Camera cam, int iter, int traceDepth, PathSegment* pathSegments) {
+__global__ void generateRayFromCamera(
+	DevScene* scene, Camera cam, int iter, int traceDepth, PathSegment* pathSegments
+) {
 
 	int x = (blockIdx.x * blockDim.x) + threadIdx.x;
 	int y = (blockIdx.y * blockDim.y) + threadIdx.y;
 
 	if (x < cam.resolution.x && y < cam.resolution.y) {
 		int index = x + (y * cam.resolution.x);
 		PathSegment& segment = pathSegments[index];
-		thrust::default_random_engine rng = makeSeededRandomEngine(iter, index, 0);
+		Sampler rng = makeSeededRandomEngine(iter, index, traceDepth, scene->sampleSequence);
 
 		segment.ray = sampleCamera(cam, x, y, sample4D(rng));
 		segment.throughput = glm::vec3(1.f);
@@ -171,7 +173,7 @@ __global__ void previewGBuffer(
 		return;
 	}
 	int index = y * width + x;
-	thrust::default_random_engine rng = makeSeededRandomEngine(iter, index, 0);
+	Sampler rng = makeSeededRandomEngine(iter, index, 0, scene->sampleSequence);
 
 	Ray ray = sampleCamera(cam, x, y, sample4D(rng));
 	Intersection intersec;
@@ -256,7 +258,7 @@ __global__ void pathIntegSampleSurface(
 	int numPaths,
 	bool sortMaterial
 ) {
-	const int SamplesConsumedOneIter = 10;
+	const int SamplesOneIter = 7;
 
 	int idx = blockDim.x * blockIdx.x + threadIdx.x;
 	if (idx >= numPaths) {
@@ -289,7 +291,7 @@ __global__ void pathIntegSampleSurface(
 		return;
 	}
 
-	thrust::default_random_engine rng = makeSeededRandomEngine(iter, idx, 4 + depth * SamplesConsumedOneIter);
+	Sampler rng = makeSeededRandomEngine(iter, idx, 4 + depth * SamplesOneIter, scene->sampleSequence);
 
 	Material material = scene->getTexturedMaterialAndSurface(intersec);
 
@@ -380,7 +382,7 @@ __global__ void singleKernelPT(
 	glm::vec3 accRadiance(0.f);
 
 	int index = y * width + x;
-	thrust::default_random_engine rng = makeSeededRandomEngine(iter, index, 0);
+	Sampler rng = makeSeededRandomEngine(iter, index, 0, scene->sampleSequence);
 
 	Ray ray = sampleCamera(cam, x, y, sample4D(rng));
 	Intersection intersec;
@@ -486,7 +488,7 @@ __global__ void BVHVisualize(int iter, DevScene* scene, Camera cam, glm::vec3* i
 	}
 	int index = y * width + x;
 
-	thrust::default_random_engine rng = makeSeededRandomEngine(iter, index, 0);
+	Sampler rng = makeSeededRandomEngine(iter, index, 0, scene->sampleSequence);
 	Ray ray = sampleCamera(cam, x, y, sample4D(rng));
 
 	Intersection intersec;
@@ -527,7 +529,7 @@ void pathTrace(uchar4* pbo, int frame, int iter) {
 		(cam.resolution.x + blockSize2D.x - 1) / blockSize2D.x,
 		(cam.resolution.y + blockSize2D.y - 1) / blockSize2D.y);
 
-	generateRayFromCamera<<<blocksPerGrid2D, blockSize2D>>>(cam, iter, Settings::traceDepth, devPaths);
+	generateRayFromCamera<<<blocksPerGrid2D, blockSize2D>>>(hstScene->devScene, cam, iter, Settings::traceDepth, devPaths);
 	checkCUDAError("PT::generateRayFromCamera");
 	cudaDeviceSynchronize();
 
diff --git a/src/sampler.h b/src/sampler.h
@@ -4,23 +4,49 @@
 #include <thrust/random.h>
 #include "mathUtil.h"
 #include "cudaUtil.h"
-
-#define SAMPLER_USE_SOBOL false
+#include "common.h"
 
 #if SAMPLER_USE_SOBOL
+#define SobolSampleNum 10000
+#define SobolSampleDim 200
+
+struct Sampler {
+	__device__ Sampler() = default;
+
+	__device__ Sampler(int ptr, uint32_t scramble, uint32_t* data) :
+		ptr(ptr), scramble(scramble), data(data) {}
+
+	__device__ float sample() {
+		uint32_t r = data[ptr++] ^ scramble;
+		scramble = Math::utilhash(scramble);
+		return r * 0x1p-32f;
+	}
+
+	uint32_t* data;
+	uint32_t scramble;
+	int ptr;
+};
+
+__device__ static Sampler makeSeededRandomEngine(int iter, int index, int dim, uint32_t* data) {
+	return Sampler(iter * SobolSampleDim + dim, Math::utilhash(index), data);
+}
+
+__device__ inline float sample1D(Sampler& sampler) {
+	return sampler.sample();
+}
+
 #else
 using Sampler = thrust::default_random_engine;
-#endif
 
-__host__ __device__
-static thrust::default_random_engine makeSeededRandomEngine(int iter, int index, int depth) {
-	int h = Math::utilhash((1 << 31) | (depth << 22) | iter) ^ Math::utilhash(index);
-	return thrust::default_random_engine(h);
+__device__ static Sampler makeSeededRandomEngine(int iter, int index, int dim, uint32_t* data) {
+	int h = Math::utilhash((1 << 31) | (dim << 22) | iter) ^ Math::utilhash(index);
+	return Sampler(h);
 }
 
-__device__ inline float sample1D(thrust::default_random_engine& rng) {
-	return thrust::uniform_real_distribution<float>(0.f, 1.f)(rng);
+__device__ inline float sample1D(Sampler& sampler) {
+	return thrust::uniform_real_distribution<float>(0.f, 1.f)(sampler);
 }
+#endif
 
 __device__ inline glm::vec2 sample2D(Sampler& sampler) {
 	return glm::vec2(sample1D(sampler), sample1D(sampler));
diff --git a/src/scene.cpp b/src/scene.cpp
@@ -38,26 +38,6 @@ MeshData* Resource::loadOBJMesh(const std::string& filename) {
     }
     bool hasTexcoord = !attrib.texcoords.empty();
 
-#if MESH_DATA_INDEXED
-    model->vertices.resize(attrib.vertices.size() / 3);
-    model->normals.resize(attrib.normals.size() / 3);
-    memcpy(model->vertices.data(), attrib.vertices.data(), attrib.vertices.size() * sizeof(float));
-    memcpy(model->normals.data(), attrib.normals.data(), attrib.normals.size() * sizeof(float));
-    if (hasTexcoord) {
-        model->texcoord.resize(attrib.texcoords.size() / 2);
-        memcpy(model->texcoords.data(), attrib.texcoords.data(), attrib.texcoords.size() * sizeof(float));
-    }
-    else {
-        model->texcoord.resize(attrib.vertices.size() / 3);
-    }
-
-    for (const auto& shape : shapes) {
-        for (auto idx : shape.mesh.indices) {
-            model->indices.push_back({ idx.vertex_index, idx.normal_index,
-                hasTexcoord ? idx.texcoord_index : idx.vertex_index });
-        }
-    }
-#else
     for (const auto& shape : shapes) {
         for (auto idx : shape.mesh.indices) {
             model->vertices.push_back(*((glm::vec3*)attrib.vertices.data() + idx.vertex_index));
@@ -68,7 +48,7 @@ MeshData* Resource::loadOBJMesh(const std::string& filename) {
             );
         }
     }
-#endif
+
     std::cout << "\t\t[Vertex count = " << model->vertices.size() << "]" << std::endl;
     meshDataPool[filename] = model;
     return model;
@@ -163,8 +143,6 @@ void Scene::createLightSampler() {
 }
 
 void Scene::buildDevData() {
-#if MESH_DATA_INDEXED
-#else
     int primId = 0;
     for (const auto& inst : modelInstances) {
         const auto& material = materials[inst.materialId];
@@ -196,7 +174,7 @@ void Scene::buildDevData() {
             }
         }
     }
-#endif
+
     if (primId == 0) {
         std::cout << "[No mesh data loaded, quit]" << std::endl;
         exit(-1);
@@ -503,6 +481,14 @@ void DevScene::create(const Scene& scene) {
     }
 
     checkCUDAError("DevScene::meshData");
+
+#if SAMPLER_USE_SOBOL
+    std::ifstream sobolFile("sobol_10k_200.bin", std::ios::in | std::ios::binary);
+    std::vector<char> sobolData(SobolSampleNum * SobolSampleDim * sizeof(uint32_t));
+    sobolFile.read(sobolData.data(), byteSizeOf(sobolData));
+    cudaMalloc(&sampleSequence, byteSizeOf(sobolData));
+    cudaMemcpyHostToDev(sampleSequence, sobolData.data(), byteSizeOf(sobolData));
+#endif
 }
 
 void DevScene::destroy() {
@@ -524,4 +510,6 @@ void DevScene::destroy() {
     cudaSafeFree(lightUnitRadiance);
     lightSampler.destroy();
     envMapSampler.destroy();
+
+    cudaSafeFree(sampleSequence);
 }
diff --git a/src/scene.h b/src/scene.h
@@ -26,9 +26,6 @@ struct MeshData {
     std::vector<glm::vec3> vertices;
     std::vector<glm::vec3> normals;
     std::vector<glm::vec2> texcoords;
-#if MESH_DATA_INDEXED
-    std::vector<glm::ivec3> indices;
-#endif
 };
 
 struct ModelInstance {
@@ -423,6 +420,8 @@ struct DevScene {
     float sumLightPowerInv;
     DevTextureObj* envMap = nullptr;
     DevDiscreteSampler1D<float> envMapSampler;
+
+    uint32_t* sampleSequence = nullptr;
 };
 
 class Scene {
