About the digital processing part, you should avoid using the general purpose sharpening methods like unsharp mask, as these methods will only increase the local contrast — making details more visible but you won't get details back that have become invisible. It is better to use methods that are based on actually reversing the blurring due to imperfect focus, and for that you need to know the point spread function (so, in the part of the image that is not in focus, a single point will appear to be a small disk with some brightness profile, this is the so-called point spread function).
You can calculate the point spread function by zooming in to high contrast areas of the image. If you know that there is a sharp edge across which the brightness changes by some amount, you can easily calculate what the point spread function is that would yield the profile you see in the image.
A rough approximation is obtained by assuming that the point spread function is a uniform disk of radius R. If you zoom in to some sharp edge then the curvature will become very small in the magnified picture so you can assume that you have a straight line, on one side the brightness is v1 on the other side of the line it will be v2. The brightness a distance d from the line v(d) in the picture will be a smooth function due to the unsharpness that tends to v1 on one side and v2 on the other side. Near d = 0, the function g(d) = [v(d) - v1]/[2(v2 - v1)] behaves as follows:
g(d) = 1/4 - d/(pi R) + d^3/(6 pi R^3) +...
So, by making a linear fit of the function g(d) - 1/4 in terms of the x and y coordinates of the picture near the line, you will get a result of the form:
g(x,y) = A + b x + c y
and then it follows that:
1/(pi R) = sqrt[b^2 + c^2]
So, the point spread function can be computed with some effort from the picture, if you are not too bad at math.
Then if you have computed the point spread function, then inverting the defocus blur is a piece of cake using algorithms such as Wiener deconvolution or Richardson–Lucy deconvolution. Such algorithms are usually included in image processing software, but you should run them by using the actual point spread function that applies to your image, not some standard Gaussian blur. E.g. this ImageJ plugin has various deconvolution algorithms that require you to specify the point spread function. And ImageJ can be obtained here.
I should add here that this should be done in a linear color space. So, you need to transform to linear RGB or to XYZ color space first, do the sharpening operation there and then transform back to sRGB.