I have it somewhat working. However, doing it this way is a lot slower.

            for each vertex v in vertex list
            {
                Vector3 n = Vector3.Zero 
                for each triangle t that shares v's vertex 
                {
                    n += Normalize(v+ t.normal)  
                }
                v.normal = Normalize(n) 
            }

That doesn't look right to me, though. I take it v is the vertex position. You can't mix that with normals. It's also evident from your screenshot: There's a color gradient going from start to end of your chunks.

You can take gradient calculation from that GPU Gems article (Example 1-1). The uvw is the vertex position here, which you feed to your noise/fractal generator.

And yes, it will get slower and slower, since you obviously don't cache the results.

Cool thanks.

My next question is this line

float d = 1.0/(float)voxels_per_block;

What's the voxels per block variable? And does the noise generator have to be anything like what I used for generator the scalar field?

It's the size of your grid cell. Or rather - if you have any scaling going on - the distance from point to point in one dimension you use to sample your scalar field (which is by the way what I meant with the noise generator). In short: d is meant to be in scalar field space. If you don't scale, forget what I just said.

You could also play with that d value a bit, depends on how your field was generated in the first place. Then again, your field shouldn't produce higher frequencies than this d value (d is sort of the wavelength, sampling theorem applies), so the grid cell dimension is probably fine.

Edit: You might even live with sampling at the very same spots you use for the grid cell corners. Though I guess it wouldn't be so far away from what you get using the face normal averaging. Anyway, an additional warning: Chunk seams will now be challenge.

After reading the source code. It seems like they are just sampling the normals from a raw file.

So I'm not sure how to translate that into perlin noise code.

        private Vector3 SmoothNormal(Vector3 position)
        {
            float dxz = 1.0f / 16.0f;
            float dy = 1.0f / 64.0f;
            Vector3 grad = Vector3.zero;

            grad.x = PerlinNoise.Noise(position.x + dxz, position.y, position.z) - PerlinNoise.Noise(position.x - dxz, position.y, position.z);
            grad.y = PerlinNoise.Noise(position.x, position.y + dy, position.z) - PerlinNoise.Noise(position.x, position.y - dy, position.z);
            grad.z = PerlinNoise.Noise(position.x, position.y, position.z + dxz) - PerlinNoise.Noise(position.x, position.y, position.z - dxz);

            grad.Normalize();
            return grad *= -1;
        }

position is the vertex position

I improved the averaging of vertex normals. I just need to improve the speed of it and chunk seams.