This hardly feels like a complete solution, but I haven't got the reputation to contribute it as a comment or in chat. Here's what I've got:
You can use the undocumented ReturnMeshObject method like @Simon Woods used here to get ListSurfacePlot3D to do the smoothing for you. With this option added, it returns a GraphicsComplex ready to be used by DiscretizeGraphics.
Graphics`Mesh`MeshInit[];
mc = MeshCoordinates[surface];
extractedmesh = DiscretizeGraphics[
First@ListSurfacePlot3D[mc, Method -> {"ReturnMeshObject" -> True},
Mesh -> {15, 15, 15}, MaxPlotPoints -> 10]]
This returns a new mesh with many more interpolated sample points. One can play with the plotting parameters or directly operate on the new region to get the surface to acceptable smoothness.
As for BSplineSurface et al., the approach taken here looks good, where @Belisarius chopped the datapoints into z-slices, interpolated, and used regular points from the interpolated curves to construct BSplineFunctions. I'm presently trying to replicate this and will add that if I succeed.