Looking at the images above, it seems to me that the boundary of the articular surface is not marked everywhere by sudden change in the surface normal direction (see green arrows). No matter how you smooth the normals or constrain region growing, you would not get the actual articular surface in these areas.
Instead, you could determine what part of the mesh is the articular surface by rotating the mesh around a bit around an axis and keep that mesh regions that is close to the original (non-rotated) surface. The center of rotation on the surface (marked with an X in the image) and the axis orientation can be easilty set approximately. You only need 4 parameters (shift of the center of rotation on the surface in two directions; tilting of the rotation axis in two directions) to get the true, optimal center of rotation axis position from the initial guess, which is a simple optimization problem (small number of unknowns, smooth objective function, fast computation of the maximized value). It should be very easy to write a small Python function that computes the size of the articular surface for the current axis position and angle (by rotating the surface a bit, computing distance using polydata distance filter, applying threshold filter, getting largest connected component) and use it in scipy.optimize. The axis position and orientation could be very useful for computing many other metrics for the joint.