You're not looking for edges (=borders between extended areas of high and low gray value), you're looking for ridges (thin lines darker or brighter than their neighborhood), so edge filters might not be ideal: An edge filter will give you two flanks (one on each side of the line) and a low response in the middle of the line:
ADD: If've been asked to explain the difference between an edge detector and a ridge detector more clearly. I apologize in advance if this answer is getting very long.
An edge detector is (usually) a first derivative operator: If you imagine the input image as a 3D landscape, an edge detector measures the steepness of the slope at each point of that landscape:
If you want to detect the border of an extended bright or dark region, this is just fine. But for the veins in the OP's image it will give you just the same: the outlines left and right of each vein:
That also explains the "double line pattern" in the Canny edge detector results:
So, how do you detect these thin lines (i.e. ridges), then? The idea is that the pixel values can be (locally) approximated by a 2nd order polynomial, i.e. if the image function is $g$, then for small values of $x$ and $y$:
$g(x,y)\approx \frac{1}{2} x^2 \frac{\partial ^2g}{\partial x^2}+x y \frac{\partial ^2g}{\partial x\, \partial y}+\frac{1}{2} y^2 \frac{\partial ^2g}{\partial y\, ^2}+x \frac{\partial g}{\partial x}+y \frac{\partial g}{\partial y}+g(0,0)$
or, in matrix form:
$g(x,y)\approx \frac{1}{2} \left(
\begin{array}{c}
x & y
\end{array}
\right).\left(
\begin{array}{cc}
\frac{\partial ^2g}{\partial x^2} & \frac{\partial ^2g}{\partial x\, \partial y} \\
\frac{\partial ^2g}{\partial x\, \partial y} & \frac{\partial ^2g}{\partial y\, ^2}
\end{array}
\right).\left(
\begin{array}{cc}
x \\
y
\end{array}
\right)+\left(
\begin{array}{cc}
\frac{\partial g}{\partial x} & \frac{\partial g}{\partial y}
\end{array}
\right).\left(
\begin{array}{c}
x \\
y
\end{array}
\right)+g(0,0)$
The second order derivative matrix $\left(
\begin{array}{cc}
\frac{\partial ^2g}{\partial x^2} & \frac{\partial ^2g}{\partial x\, \partial y} \\
\frac{\partial ^2g}{\partial x\, \partial y} & \frac{\partial ^2g}{\partial y\, ^2}
\end{array}
\right)$ is called the "Hessian matrix". It describes the 2nd order structure we're interested in.
The 2nd order part of this function can be transformed into the sum of two parabolas $\lambda _1 x^2 + \lambda _2 y^2$ rotated by some angle, by decomposing the Hessian matrix above to a rotation times a diagonal matrix of it's eigenvalues (Matrix decomposition). We don't care about the rotation (we want to detect ridges in any orientation), so we're only interested in $\lambda _1$ and $\lambda _2$
What kind of shapes can this function approximation have? Actually, not that many:
To detect ridges, we want to find areas in the image that look like the last of the plots above, so we're looking for areas where the major eigenvalue of the Hessian is large (compared to the minor eigenvalue). The simplest way to detect that is just to calculate the major eigenvalue at each pixel - and that's what the ridge filter below does.
A ridge filter will probably give better results. I've tried Mathematica's built in RidgeFilter (which calculates the major eigenvalue of the Hessian matrix at each pixel) on your image:
As you can see, there's only a single peak for every thin dark line. Binarizing and skeletonizing yields:
After pruning the skeleton and removing small components (noise) from the image, I get this final skeleton:
Full Mathematica code:
ridges = RidgeFilter[ColorNegate@src];
skeleton = SkeletonTransform[Binarize[ridges, 0.007]];
DeleteSmallComponents[Pruning[skeleton, 50], 50]
ADD:
I'm not a Matlab expert, I don't know if it has a built in ridge filter, but I can show you how to implement it "by hand" (again, using Matematica). As I said, the ridge filter is the major eigenvalue of the Hessian matrix. I can calculate that eigenvalue symbolically in Mathematica:
$\text{eigenvalue}=\text{Last}\left[\text{Eigenvalues}\left[\left(
\begin{array}{cc}
H_{\text{xx}} & H_{\text{xy}} \\
H_{\text{xy}} & H_{\text{yy}}
\end{array}
\right)\right]\right]$
=> $\frac{1}{2} \left(H_{\text{xx}}+H_{\text{yy}}+\sqrt{H_{\text{xx}}^2+4 H_{\text{xy}}^2-2 H_{\text{xx}} H_{\text{yy}}+H_{\text{yy}}^2}\right)$
So what you have to do is calculate the second derivatives $H_{\text{xx}}$, $H_{\text{xy}}$, $H_{\text{yy}}$ (using a sobel or derivative of gaussian filter) and insert them into the expression above, and you've got your ridge filter.
