I tried do use 1000 samples on the 2-sphere surrounding the centre vertex. The mesh is manifold, but it looks like the results from Marching Cubes.

def quat_normal_from_function(f, d=0.000001):
    """Given a sufficiently smooth 3d function, f, returns a function approximating of the gradient of f.
    d controls the scale, smaller values are a more accurate approximation."""
    def norm(x, y, z):

        num_samples = 1000
        sample_points = []

        for i in range(0, num_samples):
            sample_point = V3(0, 0, 0)

            u = random.randint(0, 1000)/1000.0
            v = random.randint(0, 1000)/1000.0

            pi = 4.0*math.atan(1.0);

            theta = 2.0*pi*u
            phi = math.acos(2.0*v - 1.0)

            # generate pseudorandom location on 2-sphere of radius d
            sample_point.x = d*math.cos(theta)*math.sin(phi)
            sample_point.y = d*math.sin(theta)*math.sin(phi)
            sample_point.z = d*math.cos(phi)

            sample_points.append(sample_point)

        normal = V3(0, 0, 0)

        for sample_counter in range(0, len(sample_points)):
            if f(sample_points[sample_counter].x, sample_points[sample_counter].y, sample_points[sample_counter].z) < 0:
                normal.x += -sample_points[sample_counter].x
                normal.y += -sample_points[sample_counter].y
                normal.z += -sample_points[sample_counter].z

        if 0 == normal.x*normal.x + normal.y*normal.y + normal.z*normal.z:
            normal.x = 1;
            
        return normal.normalize()
    return norm

I found the curvature once, marked in dark gray. I'm not sure if this is related?

I don't believe normals are your problem, actually. The adaptivity (i.e. choice of vertex location) is also dependent on the function being smooth, which your is not. It expects height of the function to roughly correspond to "depth from the surface". Your function is basically zero inside the unit sphere, and infinity (or 4) outside, with little to no transition, violates this.

The julia set is really just an yes/no measurement, it's not meant as function returning a float.

The easiest thing to do is to turn off adaptivity, and use marching cubes, as in that case a yes/no measurement is all you need. Then just crank up the resolution.

It was worth trying, anyway. Thanks again for the code!

I wonder if there's a way to use the Marching Cubes surface normals to bootstrap the Dual Contouring algorithm. Yes, there is... I can generate a portion of the surface near the vertex in question, using Marching Cubes, and then average the triangle normals to get a single normal. I'll let you know when it's done.

Anything has got to be better than Marching Cubes by itself.

Any chance that you can rewrite the Marching Cubes code to have a function like

bool get_fractal_surface_triangles(vector<triangle> &triangles, float x_grid_min, float x_grid_max, float y_grid_min, float y_grid_max, float z_grid_min, float z_grid_max, size_t x_res, size_t y_res, size_t z_res)?

I tried to write such a function, but I can't figure out how to alter the marching_cubes_3d_single_cell(f, x, y, z): function.

def marching_cubes_3d_2(f, x_grid_min = -1.5, x_grid_max = 1.5, y_grid_min = -1.5, y_grid_max = 1.5, z_grid_min = -1.5, z_grid_max = 1.5, x_res = 100, y_res = 100, z_res = 100):
    """Iterates over a cells of size one between the specified range, and evaluates f to produce
        a boundary by Marching Cubes. Returns a Mesh object."""
    # For each cube, evaluate independently.
    # If this wasn't demonstration code, you might actually evaluate them together for efficiency

    x_step_size = (x_grid_max - x_grid_min) / (x_res - 1)
    y_step_size = (y_grid_max - y_grid_min) / (y_res - 1)
    z_step_size = (z_grid_max - z_grid_min) / (z_res - 1)

    x_loc = x_grid_min
    y_loc = y_grid_min
    z_loc = z_grid_min

    mesh = Mesh()
    for x in range(0, x_res):

        y_loc = y_grid_min
        
        for y in range(0, y_res):

            z_loc = z_grid_min

            for z in range(0, z_res):
                
                cell_mesh = marching_cubes_3d_single_cell(f, x_loc, y_loc, z_loc)
                mesh.extend(cell_mesh)

                z_loc += z_step_size

            y_loc += y_step_size

        x_loc += x_step_size
    return mesh

If you want surface nets out of my code, it's pretty easy. If you turn off adaptivity and run the dual contouring, you get unsmoothed surface nets. If you want smoothed surface nets, replace solve_qef_3d with numpy.mean(positions, axis=0).

I put up the C++ code for a DLL for Windows on GitHub. The DLL exports a function that takes in all of the quaternion Julia set parameters, as well as parameters related to the 3D grid extent, resolution, etc. I also included a .py file to show how to call the function from the DLL using Python (the function returns the surface normal gotten by the Marching Cubes algorithm).

https://github.com/sjhalayka/marching_cubes_dll

The relevant C++ file is mc_dll.cpp and the relevant Python file is mc_dll.py

It still suffers from a jittery surface.

Can't say that I didn't try.