added functionality to pytorch Module

added octohedral mapping
added feature maps

.gitignore

@@ -78,6 +78,7 @@ cmake-build-*/
 *.vtk
 *.mesh
 
+Data/EnvironmentMaps/*
 Data/CloudDataSets/*
 Data/CameraPaths/*
 Data/Checkpoints/*

Data/Shaders/VPT/Clouds.glsl

@@ -72,8 +72,8 @@ void main() {
 #endif
     );
 #elif defined(USE_SPECTRAL_DELTA_TRACKING)
-    ScatterEvent firstEvent = ScatterEvent(false, x, 0.0, w, 0.0);
-    vec3 result = deltaTrackingSpectral(x, w);
+    ScatterEvent firstEvent = ScatterEvent(false, x, 0.0, w, 0.0, 0.0, 0.0);
+    vec3 result = deltaTrackingSpectral(x, w, firstEvent);
 #elif defined(USE_RATIO_TRACKING)
     ScatterEvent firstEvent;
     vec3 result = ratioTracking(x, w, firstEvent);
@@ -94,12 +94,37 @@ void main() {
     }
 #endif
 
+    // Accumulate cloudOnly
+    vec4 cloudOnlyOld = frame == 0 ? vec4(0) : imageLoad(cloudOnlyImage, imageCoord);
+    vec4 cloudOnly = firstEvent.hasValue ? vec4(result, 1) : vec4(0);
+    cloudOnly = mix(cloudOnlyOld, cloudOnly, 1.0 / float(frame + 1));
+    imageStore(cloudOnlyImage, imageCoord, cloudOnly);
+
     // Accumulate result
     vec3 resultOld = frame == 0 ? vec3(0) : imageLoad(accImage, imageCoord).xyz;
     result = mix(resultOld, result, 1.0 / float(frame + 1));
     imageStore(accImage, imageCoord, vec4(result, 1));
     imageStore(resultImage, imageCoord, vec4(result,1));
 
+    vec4 positionOld = frame == 0 ? vec4(0) : imageLoad(positionImage, imageCoord);
+    vec4 position = firstEvent.hasValue ? vec4(firstEvent.x, 1) : vec4(0);
+    position = mix(positionOld, position, 1.0 / float(frame + 1));
+    imageStore(positionImage, imageCoord, position);
+
+    vec4 depthOld = frame == 0 ? vec4(0) : imageLoad(depthDensityImage, imageCoord);
+    vec4 depth = firstEvent.hasValue ? vec4(firstEvent.depth, firstEvent.depth * firstEvent.depth, firstEvent.density * .001, firstEvent.density * firstEvent.density * .001 * .001) : vec4(0);
+    depth = mix(depthOld, depth, 1.0 / float(frame + 1));
+    imageStore(depthDensityImage, imageCoord, depth);
+
+
+    vec2 octoUV = worldToOctohedralUV(w);
+    vec3 octoCol = texture(environmentMapOctohedralTexture, octoUV).rgb;
+    octoCol = octoCol * (1.-cloudOnly.a) + cloudOnly.rgb;
+    //octoCol = octohedralUVToWorld(octoUV);
+    //octoCol.r = octoUV.x > .5?1.:0.;
+    //octoCol.g = octoUV.y > .5?1.:0.;
+    imageStore(depthDensityImage, imageCoord, vec4(octoCol,1.));
+
     //vec3 resultOld = frame == 0 ? vec3(0) : imageLoad(accImage, imageCoord).xyz;
     //result += resultOld;
     //imageStore(accImage, imageCoord, vec4(result, 1));

Data/Shaders/VPT/DecompositionTracking.glsl

@@ -33,7 +33,7 @@
 #if !defined(USE_NANOVDB) || defined(USE_SUPER_VOXEL_GRID)
 
 vec3 analogDecompositionTracking(vec3 x, vec3 w, out ScatterEvent firstEvent) {
-    firstEvent = ScatterEvent(false, x, 0.0, w, 0.0);
+    firstEvent = ScatterEvent(false, x, 0.0, w, 0.0, 0.0, 0.0);
 
 #ifdef USE_NANOVDB
     pnanovdb_readaccessor_t accessor = createAccessor();
@@ -134,7 +134,7 @@ vec3 analogDecompositionTracking(vec3 x, vec3 w, out ScatterEvent firstEvent) {
  * Analog decomposition tracking implemented using NanoVDB's HDDA algorithm for traversing the sparse grid.
  */
 vec3 analogDecompositionTracking(vec3 x, vec3 w, out ScatterEvent firstEvent) {
-    firstEvent = ScatterEvent(false, x, 0.0, w, 0.0);
+    firstEvent = ScatterEvent(false, x, 0.0, w, 0.0, 0.0, 0.0);
 
     pnanovdb_readaccessor_t accessor = createAccessor();
     pnanovdb_buf_t buf = pnanovdb_buf_t(0);

Data/Shaders/VPT/DeltaTracking.glsl

@@ -30,11 +30,13 @@
  * P. Kutz, R. Habel, Y. K. Li, and J. Novák. Spectral and decomposition tracking for rendering heterogeneous volumes.
  * ACM Trans. Graph., 36(4), Jul. 2017.
  */
-vec3 deltaTrackingSpectral(vec3 x, vec3 w) {
+vec3 deltaTrackingSpectral(vec3 x, vec3 w, out ScatterEvent firstEvent) {
 #ifdef USE_NANOVDB
     pnanovdb_readaccessor_t accessor = createAccessor();
 #endif
 
+    firstEvent = ScatterEvent(false, x, 0.0, w, 0.0, 0.0, 0.0);
+
     float majorant = maxComponent(parameters.extinction);
 
     vec3 weights = vec3(1, 1, 1);
@@ -88,12 +90,33 @@ vec3 deltaTrackingSpectral(vec3 x, vec3 w) {
             float xi = random();
 
             if (xi < Pa) {
+                if (!firstEvent.hasValue) {
+                    firstEvent.x = x;
+                    firstEvent.pdf_x = 0; // TODO
+                    firstEvent.w = vec3(0.);
+                    firstEvent.pdf_w = 0;
+                    firstEvent.hasValue = true;
+                    firstEvent.density = density * maxComponent(parameters.extinction);
+                    firstEvent.depth = tMax - d + t;
+                }
+
                 return vec3(0); // weights * sigma_a / (majorant * Pa) * L_e; // 0 - No emission
             }
 
             if (xi < 1 - Pn) { // scattering event
                 float pdf_w;
                 w = importanceSamplePhase(parameters.phaseG, w, pdf_w);
+
+                if (!firstEvent.hasValue) {
+                    firstEvent.x = x;
+                    firstEvent.pdf_x = 0; // TODO
+                    firstEvent.w = w;
+                    firstEvent.pdf_w = pdf_w;
+                    firstEvent.hasValue = true;
+                    firstEvent.density = density * maxComponent(parameters.extinction);
+                    firstEvent.depth = tMax - d + t;
+                }
+
                 if (rayBoxIntersect(parameters.boxMin, parameters.boxMax, x, w, tMin, tMax)) {
                     x += w*tMin;
                     d = tMax - tMin;
@@ -128,7 +151,7 @@ vec3 deltaTracking(
     float depth = 0.0;
 #endif
 
-    firstEvent = ScatterEvent(false, x, 0.0, w, 0.0);
+    firstEvent = ScatterEvent(false, x, 0.0, w, 0.0, 0.0, 0.0);
 
 #ifdef USE_NANOVDB
     pnanovdb_readaccessor_t accessor = createAccessor();
@@ -205,6 +228,15 @@ vec3 deltaTracking(
             float xi = random();
 
             if (xi < Pa) {
+                if (!firstEvent.hasValue) {
+                    firstEvent.x = x;
+                    firstEvent.pdf_x = 0;
+                    firstEvent.w = vec3(0.);
+                    firstEvent.pdf_w = 0;
+                    firstEvent.hasValue = true;
+                    firstEvent.density = density * parameters.extinction.x;
+                    firstEvent.depth = tMax - d + t;
+                }
                 return vec3(0); // weights * sigma_a / (majorant * Pa) * L_e; // 0 - No emission
             }
 
@@ -221,6 +253,8 @@ vec3 deltaTracking(
                     firstEvent.w = w;
                     firstEvent.pdf_w = pdf_w;
                     firstEvent.hasValue = true;
+                    firstEvent.density = density * parameters.extinction.x;
+                    firstEvent.depth = tMax - d + t;
                 }
 
                 if (rayBoxIntersect(parameters.boxMin, parameters.boxMax, x, w, tMin, tMax)) {

Data/Shaders/VPT/OctohedralMapper.glsl

@@ -0,0 +1,52 @@
+-- Compute
+
+#version 460
+
+#extension GL_EXT_scalar_block_layout : enable
+
+layout(local_size_x = BLOCK_SIZE, local_size_y = BLOCK_SIZE) in;
+
+layout(binding = 0) uniform sampler2D environmentMapTexture;
+
+layout(binding = 1, rgba32f) uniform writeonly image2D outputImage;
+
+const float PI = 3.14159265359;
+
+vec3 octohedralUVToWorld(vec2 uv) {
+    uv = uv * 2. - 1.;
+    vec2 sgn = sign(uv);
+    uv = abs(uv);
+
+    float r = 1. - abs(1. - uv.x - uv.y);
+    float phi = .25 * PI * ( (uv.y - uv.x) / r + 1. );
+    if (r == 0.){
+        phi = 0.;
+    }
+
+    float x = sgn.x * cos(phi) * r * sqrt(2 - r * r);
+    float y = sgn.y * sin(phi) * r * sqrt(2 - r * r);
+    float z = sign(1. - uv.x - uv.y) * (1. - r * r);
+    return vec3(x, y, z);
+}
+
+
+vec3 sampleSkybox(in vec3 dir) {
+    // Sample from equirectangular projection.
+    vec2 texcoord = vec2(atan(dir.z, dir.x) / (2. * PI) + 0.5, -asin(dir.y) / PI + 0.5);
+    vec3 textureColor = texture(environmentMapTexture, texcoord).rgb;
+    //textureColor = vec3(texcoord, 0.);
+    return textureColor;
+}
+
+void main() {
+    ivec2 writePos = ivec2(gl_GlobalInvocationID.xy);
+    ivec2 outputImageSize = imageSize(outputImage);
+
+    vec2 uv = vec2(gl_GlobalInvocationID.xy) / vec2(outputImageSize);
+
+    if (writePos.x < outputImageSize.x && writePos.y < outputImageSize.y) {
+        vec3 dir = octohedralUVToWorld(uv);
+        vec3 skybox = sampleSkybox(dir);
+        imageStore(outputImage, writePos, vec4(skybox, 1));
+    }
+}

Data/Shaders/VPT/RatioTracking.glsl

@@ -24,7 +24,7 @@
 
 #ifdef USE_RATIO_TRACKING
 vec3 ratioTracking(vec3 x, vec3 w, out ScatterEvent firstEvent) {
-    firstEvent = ScatterEvent(false, x, 0.0, w, 0.0);
+    firstEvent = ScatterEvent(false, x, 0.0, w, 0.0, 0.0, 0.0);
 
 #ifdef USE_NANOVDB
     pnanovdb_readaccessor_t accessor = createAccessor();

Data/Shaders/VPT/ResidualRatioTracking.glsl

@@ -83,7 +83,7 @@ float residualRatioTrackingEstimator(
 }
 
 vec3 residualRatioTracking(vec3 x, vec3 w, out ScatterEvent firstEvent) {
-    firstEvent = ScatterEvent(false, x, 0.0, w, 0.0);
+    firstEvent = ScatterEvent(false, x, 0.0, w, 0.0, 0.0, 0.0);
 
 #ifdef USE_NANOVDB
     pnanovdb_readaccessor_t accessor = createAccessor();

Data/Shaders/VPT/VptHeader.glsl

@@ -89,6 +89,13 @@ layout (binding = 10, r32f) uniform image2DArray primaryRayAbsorptionMomentsImag
 layout (binding = 11, r32f) uniform image2DArray scatterRayAbsorptionMomentsImage;
 #endif
 
+layout (binding = 12, rgba32f) uniform image2D cloudOnlyImage;
+layout (binding = 13, rgba32f) uniform image2D depthDensityImage;
+layout (binding = 14, rgba32f) uniform image2D positionImage;
+
+layout (binding = 15) uniform sampler2D environmentMapOctohedralTexture;
+
+
 /**
  * This code is part of an GLSL port of the HLSL code accompanying the paper "Moment-Based Order-Independent
  * Transparency" by Münstermann, Krumpen, Klein, and Peters (http://momentsingraphics.de/?page_id=210).
@@ -97,7 +104,7 @@ layout (binding = 11, r32f) uniform image2DArray scatterRayAbsorptionMomentsImag
  * This port is released under the terms of the MIT License.
  */
 /*! This function implements complex multiplication.*/
-layout(std140, binding = 12) uniform MomentUniformData {
+layout(std140, binding = 16) uniform MomentUniformData {
     vec4 wrapping_zone_parameters;
     //float overestimation;
     //float moment_bias;

Data/Shaders/VPT/VptUtils.glsl

@@ -34,6 +34,8 @@ struct ScatterEvent {
     bool hasValue;
     vec3 x; float pdf_x;
     vec3 w; float pdf_w;
+    float depth;
+    float density;
 };
 
 
@@ -129,12 +131,52 @@ vec3 sRGBToLinearRGB(in vec3 linearRGB) {
     return mix(pow((linearRGB + 0.055) / 1.055, vec3(2.4)), linearRGB / 12.92, lessThanEqual(linearRGB, vec3(0.04045)));
 }
 
+vec3 octohedralUVToWorld(vec2 uv) {
+    uv = uv * 2. - 1.;
+    vec2 sgn = sign(uv);
+    uv = abs(uv);
+
+    float r = 1. - abs(1. - uv.x - uv.y);
+    float phi = .25 * PI * ( (uv.y - uv.x) / r + 1. );
+    if (r == 0.){
+        phi = 0.;
+    }
+
+    float x = sgn.x * cos(phi) * r * sqrt(2 - r * r);
+    float y = sgn.y * sin(phi) * r * sqrt(2 - r * r);
+    float z = sign(1. - uv.x - uv.y) * (1. - r * r);
+    return vec3(x, y, z);
+}
+
+vec2 worldToOctohedralUV(vec3 dir) {
+    dir = normalize(dir);
+    vec3 sgn = sign(dir);
+    dir = abs(dir);
+
+    float phi = atan(dir.y , dir.x);
+    float r = sqrt(1. - dir.z);
+
+    float v = r * 2. / PI * phi;
+    float u = r - v;
+
+    vec2 uv = vec2(u,v);
+    if (sgn.z < 0.){
+        uv = vec2(1.-v, 1.-u);
+    }
+    uv *= sgn.xy;// + 1.- abs(sgn.xy);
+
+    return uv * .5 + .5;
+}
 
 #ifdef USE_ENVIRONMENT_MAP_IMAGE
 vec3 sampleSkybox(in vec3 dir) {
     // Sample from equirectangular projection.
     vec2 texcoord = vec2(atan(dir.z, dir.x) / TWO_PI + 0.5, -asin(dir.y) / PI + 0.5);
     vec3 textureColor = texture(environmentMapTexture, texcoord).rgb;
+
+    // Make sure there is no infinity in the skybox
+    textureColor = min(textureColor , vec3(10000, 10000, 10000));
+
 #ifdef ENV_MAP_IMAGE_USES_LINEAR_RGB
 #ifdef USE_LINEAR_RGB
     return parameters.environmentMapIntensityFactor * textureColor;