this is nothing trivial. Currently there is no good solution to do this with native Grasshopper. Surface approximation is difficult, because you need to fit a surface as close as possible to a mesh without making it a complete bad surface. This means you also need some sort of smoothing algorithm.
Furthermore approximating a curve and creating a surface patch out of it does not work either, because if you match a boundary, that doesn’t match the inner part of your surface at all.
In Rhino 5 there was are great plugin called VSR, which is obsolete. However there are some commerical “reverse engineering” software solutions out there like Resurf for Rhino, Mesh2Surface and many others.
But of course you can try to script it. Usually a great way to start is using Linear regression to quickly match to a raw surface as close as possible. Since it is bit difficult to map this to a Nurbs polynom, it might be better to move the controlpoints just recursivly “up” and “down” until you fit well enough. In theory, and due to its regression-alike character, it could be solved using Galapagos, which however is not a good way of doing things. Smoothing is another deal. I could thing of something like a soap effect, which reduces the surface tension or you just directly reduce the cp count as much as possible, with the intention to reduce the curvature delta. This might require to move the cps also in other directions as in projection direction (other than “up” and “down”)