Segmentation based on normal directions

Bmesh script

Result of running script on grid above, each separated object is represented with a different colour Will notice some quad faces are "bent". Suggest Triangulating first

Find the average normal. Remember for objects like the cube this could be a zero vector. Could also consider using an arbitrary vector eg z axis rather than average normal

Since want to separate into surfaces have used face normals.

Then run thru the script and look at the tangent (grade) of the face normal to the average face normal.

Test script. Splits the current edit mesh into chunks with 0.05 grade intervals from the average normal.

When the script has run will still be in edit mode with a faceless mesh, all having been moved to separate objects.

import bpy
import bmesh
from mathutils import Vector
from math import tan
grade = 0.05

context = bpy.context
ob = context.object
me = ob.data
bm = bmesh.from_edit_mesh(me)

# use the average of all face normals
norm = sum((f.normal for f in bm.faces), Vector()) / len(bm.faces)

# use the local Z axis
#norm = (0, 0, 1)

# use the global z axis
#norm = ob.matrix_world.inverted() @ Vector((0, 0, 1))

# use the active element (will error if edge or none)
#norm = bm.select_history.active.normal

while bm.faces:
    for f in bm.faces:
        f.select = tan(f.normal.angle(norm)) < grade
    if any(f.select for f in bm.faces):
        bpy.ops.mesh.separate()
    grade += grade
    
# pop out of edit mode remove original