Reconstructing axis of pot with one shard via vertex normals
Proof of Concept

Result on single selected edge of convex hull
Further to comment, have added this as going some way to proof of concept,
Fistly here is a script that copies your object and mesh to another, goes into edit mode and converts it to a convex hull.
Select the shard and run
import bpy
bpy.ops.object.mode_set()
bpy.ops.object.duplicate(linked=False)
dupe = bpy.context.object
dupe.display_type = 'WIRE'
bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.select_all(action='SELECT')
bpy.ops.mesh.convex_hull()

after which new wire frame of the convex hull of original in edit mode with all geometry selected.
The next script goes thru the edges of the hull, finds the closest point on the mesh to its middle point, uses these to create a circle from chord as described here How can I create a mathematically correct arc/circular segment?
To visualize have added a vert at circle center and the two joining edges.  As data would save as radius, centre coordinate, and normal (axis of rotation the normalized cross product of two edge vectors)
Test script, creates predicted circle &quot;wedges&quot; for each selected edge.  Run with convex hull mesh in edit mode, edges of interest selected.

Result on all edges of convex hull
import bpy
import bmesh
from math import asin, degrees
context = bpy.context
scene = context.scene
ob = context.object
me = ob.data
bm = bmesh.from_edit_mesh(me)
shard = scene.objects.get(&quot;3D_Scherbe_Model_50K&quot;)
#edges = bm.edges[:]  # all edges 
edges = [e for e in bm.edges if e.select]
#edges = [e for e in bm.select_history if isinstance(e, bmesh.types.BMEdge)]  
for edge in edges:

    o = (edge.verts[1].co + edge.verts[0].co) / 2

    hit, loc, _, _ = shard.closest_point_on_mesh(o)

    if hit:
        h = (loc - o).length
        if h &lt; 0.1:
            print(&quot;On surface&quot;)
            continue
        a = edge.calc_length() / 2
        r = (a * a + h * h) / (2 * h)
        if abs(a / r) &gt; 1:
            # math domain error on arcsin
            print(&quot;N/A&quot;)
        else:
            angle = 2 * asin(a / r)    
            print(f&quot;{r} {degrees(angle)}&quot;)
            vc = bm.verts.new(o + r * (o - loc).normalized())
            for v in edge.verts:
                bm.edges.new((v, vc))  
bmesh.update_edit_mesh(me)
me.update()

Notes.

Instead of predicting a circle from the closest point to edge centre, could walk the edge and make  sample points to crunch into   https://meshlogic.github.io/posts/jupyter/curve-fitting/fitting-a-circle-to-cluster-of-3d-points/ and https://github.com/ndvanforeest/fit_ellipse as suggested by @RobinBetts.

Similarly could use the generated circle estimate to test against actual mesh surface.

Look at the normal returned from closest point on mesh.

Narrow down the selection, is there historic data that suggests radii or wedge angle within a certain range.

Project (closest point on mesh) equal length &quot;subchords&quot; of edge onto mesh, if each has same radius and angle would be a perfect match.  Minimize for best fit.

Look at the bounding box dimensions. If an edge is shorter than some fraction of minimum bbox dimension it is probably not major pot axis. Consider shaving off some percentage.

Decimating cleaning or smoothing shard mesh in some way.