trsh said:
when you zoom into something it's just lots of tiny cubes.

Yeah, but that's simple to do, and just deforming the vertices to get smooth results seems an option.

Code below is how i do this, density Gradient can be used as normal for lighting.

It reminds about something similar from this post, see Surface Nets: https://0fps.net/2012/07/12/smooth-voxel-terrain-part-2/

Maybe his method is better, with smoother results. Need to try that at some time. He used the average of density / edge intersection.

He also has source code, but because i have requirement on manifolds i decided to start from scratch - given code examples always contain look up tables which makes them hard to understand.

		static void VertexPositionAndNormalFromDensityField (
			vec &displacedPosition, vec &densityGradient, 
			const vec gridPosition, 
			const int x, const int y, const int z, // grid indices 	
			const std::vector<float> &density, const int regDim, // density volume and dimension
			const float cellSize, const float densCentDisp = 2.0f, const float gradDisp = 1.0f) 
		{
			vec p = gridPosition;


			vec cd(0);
			float dSum = 0;


			float d[8];
			int n = 0;
			for (int Z=-1; Z<1; Z++)
			for (int Y=-1; Y<1; Y++)
			for (int X=-1; X<1; X++)
			{
				float D = density[Index(x+X,y+Y,z+Z, regDim)];
				cd += D * vec(float(X)+0.5f, float(Y)+0.5f, float(Z)+0.5f);
				dSum += D;
				d[n++] = D;
			}
			cd *= cellSize / dSum;
			p += cd * densCentDisp;


			vec g;
			g[0]  = d[0b001] - d[0b000];
			g[0] += d[0b011] - d[0b010];
			g[0] += d[0b101] - d[0b100];
			g[0] += d[0b111] - d[0b110];
			g[1]  = d[0b010] - d[0b000];
			g[1] += d[0b110] - d[0b100];
			g[1] += d[0b011] - d[0b001];
			g[1] += d[0b111] - d[0b101];
			g[2]  = d[0b100] - d[0b000];
			g[2] += d[0b101] - d[0b001];
			g[2] += d[0b110] - d[0b010];
			g[2] += d[0b111] - d[0b011];
			g *= cellSize / 4;
			p -= g * gradDisp;


			displacedPosition = p;
			densityGradient = g;
		} 

Oh, notice such vertex deformation can be done for each quad vertex without knowing neighbours and connectivity.

This means redunant work (precessing the same vertex position multiple times), but may outweight the complexity to detect shared vertices.

So you would have more vertices and bandwidth, but very low complexity - seems very GPU friendly to me.

Also much less programming work.

Need to update my code snippet from above, which turned out pretty wrong.

This improved version snaps to the surface exactly, no more blockyness or tuning parameters. Nice smooth surface :D

		static void VertexPositionAndNormalFromDensityField (
			vec &displacedPosition, vec &surfaceNormal, 
			const vec gridPosition, 
			const int x, const int y, const int z, // grid indices 	
			const std::vector<float> &density, const int regDim, // density volume and dimension
			const float cellSize) 
		{
			vec p = gridPosition;

			float dSum = 0;
			float d[8];
			int i = 0;
			for (int Z=-1; Z<1; Z++)
			for (int Y=-1; Y<1; Y++)
			for (int X=-1; X<1; X++)
			{
				float D = density[Index(x+X,y+Y,z+Z, regDim)] - 0.5f;
				dSum += D;
				d[i++] = D;
			}
			float avD = dSum/8;

			vec n;
			n[0]  = d[0b001] - d[0b000];
			n[0] += d[0b011] - d[0b010];
			n[0] += d[0b101] - d[0b100];
			n[0] += d[0b111] - d[0b110];
			
			n[1]  = d[0b010] - d[0b000];
			n[1] += d[0b110] - d[0b100];
			n[1] += d[0b011] - d[0b001];
			n[1] += d[0b111] - d[0b101];
			
			n[2]  = d[0b100] - d[0b000];
			n[2] += d[0b101] - d[0b001];
			n[2] += d[0b110] - d[0b010];
			n[2] += d[0b111] - d[0b011];

			float mag = n.Length();
			if (mag > FP_EPSILON2)
			{
				n /= mag;
			
				float numer = 0, denom = 0;
				int i = 0;
				for (int Z=-1; Z<1; Z++)
				for (int Y=-1; Y<1; Y++)
				for (int X=-1; X<1; X++)
				{
					float t = n.Dot( vec(float(X)+0.5f, float(Y)+0.5f, float(Z)+0.5f) );
					float r = d[i] - avD;
					if (t<0) {t = -t; r = -r;}
					numer += t;
					denom += r;
					i++;
				}
				if (denom > FP_EPSILON2) 
				{
					float fieldScale = numer / denom;
					p += n * -avD * fieldScale * cellSize;
				}
			}
			displacedPosition = p;
			surfaceNormal = n;
		}

@JoeJ What are u building?

@trsh Long story. Primary goal is seamless UV parametrization for my realtime GI, but i aim for automatic LOD generation as well so all this work is worth it.

Actually the pipeline is soild voxelization → isosurface → quadrangulation → seamless UVs → do all this per LOD → generate surfel hierarchy for lighting and mesh hierarchy for rendering.

I'll also try geomorph / progressive mesh ideas for continuous LOD. And stochastic point splatting would be very easy, must try, and displacement mapping on every surface would work too because of seamless UVs.

I did not really expect to get lost in geometry processing so deep, but in the end i hope it works out…

@JoeJ Yeah but what the end product? Game? Whats the gameplay?

Middleware for lighting. Game is plan B.

@JoeJ do you have a public repo for this project?