Fixes most common artifacts between bricks (issue #99, see also issue #…

…75) in rayDeepBrick and raySurfaceTrilinearBrick.

Also hides other artifacts in gResample and gSprayDeposit by using a larger value of epsilon. Partly
fixes gvdbRaytrace by using trilinear interpolation instead of tricubic (since tricubic requires
an apron size of 2 for correctness, which isn't supported yet.) Finally, adds some more
comments to HDDAState.

Short version: If one still sees artifacts, they should try using VolumeGVDB::SetEpsilon(0.01f, 256);
this seems to fix the remaining inter-brick volume artifact, while introducing some banding artifacts
that are hard to see. I have some ideas for an HDDA approach that might avoid having to use an
epsilon in the algorithm at all.

For the artifacts in issue #99, it turned out that there were two issues in raySurfaceTrilinearBrick.
One was that a line that quantized samples to discrete increments of t was commented out. Since the
value of t when entering a brick otherwise depends on where within a voxel the ray entered the brick
(I think), this led to visible bands the size of voxels. rayDeepBrick also had this problem, except
there the issue was that quantization was performed after determining the voxel position, instead of
before.

The other issue in raySurfaceTrilinearBrick was that the loop to skip empty voxels inside bricks
(which appears to have been intended as an optimization, but I'm not sure how effective it was)
was incorrect - it would always step by SCN_DIRECTSTEP at least once, even if the initial sample
was inside a brick.

Together, these two issues appear to have been responsible for visible banding along the planes
separating bricks, and for visible gaps between bricks in the volume view of g3DPrint and several
other samples.

However, volume rendering isn't totally artifact-free on all samples: in some samples, positioning
one's camera just right and at a particular range of distances - such that one of the planes between
bricks is greatly foreshortened - will show floating-point glitches between bricks. (This only seems
to happen for some planes and on some samples - e.g. I haven't seen it on g3DPrint yet - but is
possible to reliably reproduce on gDepthMap and gInteractiveGL). This appears to be due to the
value of epsilon, which is used in rayCast to move samples slightly into nodes. In these cases,
setting epsilon to a larger value - e.g. 0.01f - using VolumeGVDB::SetEpsilon(0.01f, 256) - appears
to fix these artifacts, although it adds some subtle banding artifacts that look like circles in gDepthMap.
As such, I've only applied this fix to gResample and gSprayDeposit, where it doesn't appear to add new
visible artifacts.

I'll probably take a look at replacing HDDAState with a new approach that expresses index-space position
differently; hopefully this should avoid having to use an epsilon, although functions that access HDDAState's
class members directly will probably have to be modified.

This also fixes some problems with gSprayDeposit! It turns out that some artifacts were because not
all aprons were being updated. Another issue is that gvdbRaytrace ignored the shading method passed
to it and uses tricubic intersection instead. Because tricubic intersection technically requires an
apron size of 2 (which isn't supported at the moment), this lead to rays sometimes intersecting the
boundaries between bricks, which presented as adding density to spaces in the air. Changing gvdbRaytrace
to use trilinear intersection instead fixes this - but the correct solution would be to respect the
shading parameter passed to the function.

source/gResample/main_resample.cpp

@@ -214,6 +214,7 @@ bool Sample::init()
 	gvdb.getScene()->LinearTransferFunc ( 0.50f, 0.75f, Vector4DF(1,0,0,0.02f), Vector4DF(.2f,.2f,0.2f,0.02f) );
 	gvdb.getScene()->LinearTransferFunc ( 0.75f, 1.00f, Vector4DF(.2f,.2f,0.2f,0.02f), Vector4DF(.1f,.1f,.1f,.1f) );
 	gvdb.CommitTransferFunc ();
+	gvdb.SetEpsilon(0.01f, 256); // Use a larger value of epsilon than the default (0.001) to avoid artifacts between bricks
 
 	// Create Camera 
 	Camera3D* cam = new Camera3D;						

source/gSprayDeposit/main_spray_deposit.cpp

@@ -148,6 +148,7 @@ bool Sample::init()
 	gvdb.getScene()->LinearTransferFunc ( 0.2f, 0.3f,   Vector4DF(1,1,1,0.05f), Vector4DF(1,1,1,0.05f) );	
 	gvdb.getScene()->LinearTransferFunc ( 0.3f, 1.0f,   Vector4DF(1,1,1,0.05f), Vector4DF(0,0,0,0.0) );	
 	gvdb.CommitTransferFunc (); 
+	gvdb.SetEpsilon(0.01f, 256); // Use a larger epsilon than default to avoid artifacts between bricks
 
 	// Configure a new GVDB volume
 	gvdb.Configure ( 3, 3, 3, 3, 5 );
@@ -291,7 +292,7 @@ void Sample::simulate()
 	// Returns the hit position and normal of each ray
 	gvdb.Raytrace ( m_rays, 0, SHADE_TRILINEAR, 0, -0.0001f);
 
-	// Insert Points into the GVDBb grid
+	// Insert Points into the GVDB grid
 	// This identifies a grid cell for each point. The SetPointStruct function accepts an arbitrary 
 	// structure, which must contain a vec3f position input, and uint node offset and index outputs.
 	// We can use the ray hit points as input directly from the ScnRay data structure.
@@ -307,7 +308,7 @@ void Sample::simulate()
 	gvdb.InsertPointsSubcell (subcell_size, m_numrays, radius, Vector3DF(0,0,0), scPntLen);
 	gvdb.GatherDensity (subcell_size, m_numrays, radius, Vector3DF(0,0,0), scPntLen, 0, 1, true ); // true = accumulate
 
-	gvdb.UpdateApron ( 0 );
+	gvdb.UpdateApron ();
 
 	// Smooth the volume
 	// A smoothing effect simulates gradual erosion

source/gvdb_library/kernels/cuda_gvdb_dda.cuh

@@ -40,6 +40,13 @@ inline __device__ float4 transfer(VDBInfo* gvdb, float v)
 // Essentially, we're always marching within a parent node, over the children of that node.
 // We can always prepare this state given a point and a node.
 // At each step, we advance to the next child node in the X, Y, or Z directions.
+//
+// This mainly gets used in cuda_gvdb_raycast.cuh.
+// Typically, a function will call SetFromRay once per ray, then call Prepare initially
+// and whenever the traversal level changes (this essentially performs some precomputation).
+// It'll then usually alternate between calling Next (determines next cell intersection)
+// and Step (moves to that next intersection).
+//
 // Reference: http://ramakarl.com/pdfs/2016_Hoetzlein_GVDB.pdf
 struct HDDAState {
 	// Constant per ray:
@@ -54,6 +61,8 @@ struct HDDAState {
 	float3 tSide; // Value of t.x that intersects the next plane in the x, y, and z direction.
 	int3 mask; // Which coordinates to move along next, each 1 or 0.
 
+	// Initializes the HDDA given an index-space point (startPos), direction (startDir), and
+	// a starting value of t (startT).
 	__device__ void SetFromRay(float3 startPos, float3 startDir, float3 startT) {
 		pos = startPos;
 		dir = startDir;

source/gvdb_library/kernels/cuda_gvdb_module.cu

@@ -231,7 +231,7 @@ extern "C" __global__ void gvdbRaytrace ( VDBInfo* gvdb, uchar chan, int num_ray
 	// raytrace
 	rays[x].hit = make_float3(NOHIT, NOHIT, NOHIT);
 	float4 hclr = make_float4(1,1,1,1);
-	rayCast ( gvdb, chan,  rays[x].orig, rays[x].dir, rays[x].hit, rays[x].normal, hclr, raySurfaceTricubicBrick );
+	rayCast ( gvdb, chan,  rays[x].orig, rays[x].dir, rays[x].hit, rays[x].normal, hclr, raySurfaceTrilinearBrick );
 
 	if ( rays[x].hit.z != NOHIT ) rays[x].hit -= rays[x].dir * bias;
 }

source/gvdb_library/kernels/cuda_gvdb_raycast.cuh

@@ -294,10 +294,10 @@ __device__ void raySurfaceVoxelBrick ( VDBInfo* gvdb, uchar chan, int nodeid, fl
 __device__ void raySurfaceTrilinearBrick ( VDBInfo* gvdb, uchar chan, int nodeid, float3 t, float3 pos, float3 dir, float3& hit, float3& norm, float4& hclr )
 {
 	float3 vmin;
-	VDBNode* node	= getNode ( gvdb, 0, nodeid, &vmin );			// Get the VDB leaf node	
-	float3  o = make_float3( node->mValue ) ;				// Atlas sub-volume to trace	
-	float3	p = pos + t.x*dir - vmin;					// sample point in index coords			
-	t.x = SCN_DIRECTSTEP * ceil ( t.x / SCN_DIRECTSTEP );
+	VDBNode* node	= getNode ( gvdb, 0, nodeid, &vmin );	// Get the VDB leaf node	
+	float3  o = make_float3( node->mValue ) ;				// Atlas sub-volume to trace
+	t.x = SCN_DIRECTSTEP * ceil(t.x / SCN_DIRECTSTEP);		// Start on sampling wavefront (avoids subvoxel banding artifacts)
+	float3	p = pos + t.x*dir - vmin;						// sample point in index coords			
 
 	for (int iter=0; iter < MAX_ITER && p.x >=0 && p.y >=0 && p.z >=0 && p.x < gvdb->res[0] && p.y < gvdb->res[0] && p.z < gvdb->res[0]; iter++) 
 	{	
@@ -498,34 +498,22 @@ __device__ void rayShadowBrick ( VDBInfo* gvdb, uchar chan, int nodeid, float3 t
 __device__ void rayDeepBrick ( VDBInfo* gvdb, uchar chan, int nodeid, float3 t, float3 pos, float3 dir, float3& hit, float3& norm, float4& clr )
 {
 	float3 vmin;
-	VDBNode* node = getNode ( gvdb, 0, nodeid, &vmin );			// Get the VDB leaf node		
+	VDBNode* node = getNode ( gvdb, 0, nodeid, &vmin );		// Get the VDB leaf node		
 
-	//t.x = SCN_DIRECTSTEP * ceil( t.x / SCN_DIRECTSTEP );						// start on sampling wavefront	
-
-	float3 o = make_float3( node->mValue );					// atlas sub-volume to trace
-	float3 wp = pos + t.x*dir;	
-	float3 p = wp-vmin;					// sample point in index coords	
-	float3 wpt = SCN_DIRECTSTEP*dir;					// world increment
-	float4 val = make_float4(0,0,0,0);
-	float4 hclr;
-	int iter = 0;
-	float dt = length(SCN_DIRECTSTEP*dir);
+	t.x = SCN_DIRECTSTEP * ceil( t.x / SCN_DIRECTSTEP );	// Start on sampling wavefront (avoids subvoxel banding artifacts)
+
+	float3 o = make_float3( node->mValue );	// Atlas sub-volume to trace
+	float3 wp = pos + t.x*dir;				// Sample position in index space
+	float3 p = wp - vmin;					// Sample point in sub-volume (in units of voxels)
+	const float3 wpt = SCN_DIRECTSTEP*dir;	// Increment in units of voxels
+	const float dt = length(wpt);			// Change in t-value per step
 	const float tDepthIntersection = getRayDepthBufferMax(dir); // The t.x at which the ray intersects the depth buffer
 
-	// record front hit point at first significant voxel
+	// Record front hit point at first significant voxel
 	if (hit.x == 0) hit.x = t.x; // length(wp - pos);
 
-	// skip empty voxels
-	for (iter=0; val.w < SCN_MINVAL && iter < MAX_ITER && p.x >= 0 && p.y >=0 && p.z >=0 && p.x < gvdb->res[0] && p.y < gvdb->res[0] && p.z < gvdb->res[0]; iter++) {		
-		val.w = transfer ( gvdb, tex3D<float> ( gvdb->volIn[chan], p.x+o.x, p.y+o.y, p.z+o.z ) ).w;
-		p += SCN_DIRECTSTEP*dir;
-		wp += wpt;
-		t.x += dt;
-	}	
-
-	// accumulate remaining voxels
-	for (; clr.w > SCN_ALPHACUT && iter < MAX_ITER && p.x >=0 && p.y >=0 && p.z >=0 && p.x < gvdb->res[0] && p.y < gvdb->res[0] && p.z < gvdb->res[0]; iter++) {			
-
+	// Accumulate remaining voxels
+	for (int iter = 0; clr.w > SCN_ALPHACUT && iter < MAX_ITER && p.x >=0 && p.y >=0 && p.z >=0 && p.x < gvdb->res[0] && p.y < gvdb->res[0] && p.z < gvdb->res[0]; iter++) {			
 		// Test to see if we've intersected the depth buffer (if there is no depth buffer, then this will never happen):
 		if (t.x > tDepthIntersection) {
 			hit.y = length(wp - pos);
@@ -534,15 +522,22 @@ __device__ void rayDeepBrick ( VDBInfo* gvdb, uchar chan, int nodeid, float3 t,
 			return;
 		}
 
-		val = transfer ( gvdb, tex3D<float> ( gvdb->volIn[chan], p.x+o.x, p.y+o.y, p.z+o.z ) );
-		val.w = exp ( SCN_EXTINCT * val.w * SCN_DIRECTSTEP );
-		hclr = (gvdb->clr_chan==CHAN_UNDEF) ? make_float4(1,1,1,1) : getColorF (gvdb, gvdb->clr_chan, p+o );
-		clr.x += val.x * clr.w * (1 - val.w) * SCN_ALBEDO * hclr.x;
-		clr.y += val.y * clr.w * (1 - val.w) * SCN_ALBEDO * hclr.y;
-		clr.z += val.z * clr.w * (1 - val.w) * SCN_ALBEDO * hclr.z;
-		clr.w *= val.w;		
+		// Get the value of the volume at this point. Only consider it if it's greater than SCN_MINVAL.
+		const float rawSample = tex3D<float>(gvdb->volIn[chan], p.x + o.x, p.y + o.y, p.z + o.z);
+		if (rawSample >= SCN_MINVAL) {
+			// Apply transfer function; integrate val.w to get transmittance according to the Beer-Lambert law:
+			float4 val = transfer(gvdb, rawSample);
+			val.w = exp(SCN_EXTINCT * val.w * SCN_DIRECTSTEP);
+			// RGB color from color channel (alpha component is unused):
+			const float4 hclr = (gvdb->clr_chan == CHAN_UNDEF) ? make_float4(1, 1, 1, 1) : getColorF(gvdb, gvdb->clr_chan, p + o);
+			clr.x += val.x * clr.w * (1 - val.w) * SCN_ALBEDO * hclr.x;
+			clr.y += val.y * clr.w * (1 - val.w) * SCN_ALBEDO * hclr.y;
+			clr.z += val.z * clr.w * (1 - val.w) * SCN_ALBEDO * hclr.z;
+			clr.w *= val.w;
+		}
 
-		p += SCN_DIRECTSTEP*dir;		
+		// Step forwards.
+		p += wpt;
 		wp += wpt;		
 		t.x += dt;
 	}