Without claiming much generality, I made the following. I'm using a slightly more complex image than your proposed one.
i = Binarize@Import@"https://i.stack.imgur.com/qDby8.png";

idi = ImageDimensions@i;
vertexI = SelectComponents[i, "Count", 5 < # < 100 &];
disk = 20 (*use some heuristics to ensure a proper vertex occlusion radii*);
p[disk_, fraction_] := IntegerPart[disk fraction](*proportionalty*)
g[x_, r___] := Graphics[x, PlotRange -> Transpose[{{0, 0}, idi}],
ImageSize -> idi, r]
vxRules = ComponentMeasurements[vertexI, "Centroid"];
vxPos = Range@Length@vxRules/. vxRules;
i1 = Binarize[Show[i, g[{White, Disk[#, disk] & /@ vxPos}]]];
i2 = ColorNegate@Erosion[i1, 1];
getMask[edges_, edge_] := SelectComponents[edges, {"Label", "Mask"}, #1 == edge &];
edges = MorphologicalComponents@DeleteSmallComponents[i2, 30];
(* masks "preserve" the mask number*)
masks = getMask[edges, #] & /@ Range@Max@Flatten@edges;
ImageAdd[#, g[{Red, (Disk[#1, disk] &) /@ vxPos}, Background -> Black]] & /@
(Image /@ masks)

(* tm may require Pruning[tm, nn] if the image is low quality *)
tm = Thinning@Image@Total@masks;
mbp = MorphologicalTransform[Binarize@tm, "SkeletonBranchPoints"];
(* get the "unique" branch points, like clustering by taking the mean of near points*)
mbpClustered = Union@MeanShift[ImageValuePositions[mbp, 1], p[disk, 1/2]];
(* Get the whole image of all multiples occluding branch points*)
segs = ImageMultiply[tm, Binarize@g[{Black, Disk[#1, p[disk, 1/4]] & /@ mbpClustered}]];
mcsegs = MorphologicalComponents[segs];
mcsegs // Colorize

I'm pretty sure the following function can be done better, for example by using @nikie's answer here
findContinuations[branchPoint_, i_, mcsegs_, disk_] :=
Module[{mm, coSegs, segmentsAtBranchPoint, tails, tgs, dests, a, b, x},
mm = Binarize@Image[g[{White, Disk[branchPoint, p[disk, 2/5]]},
Background -> Black], ImageSize -> idi];
coSegs = ImageMultiply[Image@mcsegs, mm] ;
segmentsAtBranchPoint = Select[Union@Flatten[ImageData@coSegs], # != 0 &];
tails = Position[ImageData@coSegs, #] & /@ segmentsAtBranchPoint;
tgs = a /. FindFit[#, a x + b, {a, b}, x] & /@ tails;
dests = Nearest[tgs -> segmentsAtBranchPoint, #, 2][[2]] & /@ tgs;
Sort /@ Transpose[{segmentsAtBranchPoint, dests}] // Union]
fc = findContinuations[#, i, mcsegs, disk] & /@ mbpClustered;
equiv = Flatten /@ Gather[Flatten[fc, 1], Intersection[#1, #2] =!= {} &];
rules = Reverse /@ Thread[First@# -> Rest@#] & /@ IntegerPart@equiv // Flatten;
unified = mcsegs //. rules;
f = Nearest[vxPos -> Automatic];
vxsForMask = Map[f[#, {Infinity, p[disk, 5/4]}] &,
ImageValuePositions[Image@unified, #] & /@ Range@Max@unified, {2}];
edgesFin = Rule @@@ DeleteCases[(Flatten /@ Union /@ vxsForMask /. {x_} :> {x, x}), {}];
GraphicsRow[{i,
Colorize[unified, ColorFunction -> ColorData@10,ColorFunctionScaling -> False],
GraphPlot[edgesFin, VertexCoordinateRules -> vxRules,
MultiedgeStyle -> 1/3, VertexLabeling -> True]}]

When running it on your image:

I'm using GraphPlot because in v9 multigraphs aren't supported
Finally, here you have the code "conveniently" packed into functions (usage example at the end)
p[disk_, fraction_] := IntegerPart[disk fraction](*proportionalty*)
g[x_, r___] := Graphics[x, PlotRange -> Transpose[{{0, 0}, idi}], ImageSize -> idi, r]
getProblemParms[i_Image] :=
Module[{idi, vertexI, disk, vxRules, vxPos},
idi = ImageDimensions@i;
vertexI = SelectComponents[i, "Count", 5 < # < 100 &];
disk = 20 (*find some heuristics to ensure a proper vertex occlusion radii*);
vxRules = ComponentMeasurements[vertexI, "Centroid"];
vxPos = Range@Length@vxRules /. vxRules;
{idi, disk, vxRules, vxPos}
]
getMasks[i_Image, disk_, vxPos_] := Module[{i1, i2, edges, getMask},
getMask[edges_, edge_] := SelectComponents[edges, {"Label", "Mask"}, #1 == edge &];
i1 = Binarize[Show[i, g[{White, Disk[#, disk] & /@ vxPos}]]];
i2 = ColorNegate@Erosion[i1, 1];
edges = MorphologicalComponents@DeleteSmallComponents[i2, 30];
(*masks "preserve" the mask number*)
getMask[edges, #] & /@ Range@Max@Flatten@edges
]
collectEdgesForests[masks_, disk_] :=
Module[{mIm, tm, mbp, posMbp, mbpClustered, segs},
mIm = Image@Total@masks;
tm = Thinning@mIm;
mbp = MorphologicalTransform[Binarize@tm, "SkeletonBranchPoints"];
posMbp = ImageValuePositions[mbp, 1];
(*get the "unique" branch points,
like clustering by taking the mean of near points*)
mbpClustered = Union@MeanShift[ImageValuePositions[mbp, 1], p[disk, 1/2]];
(*Get the whole image of all multiples occluding branch points*)
(*segs "preserve" the mask number*)
segs = ImageMultiply[tm, Binarize@g[{Black, Disk[#1, p[disk, 1/4]] & /@ mbpClustered}]];
{mbpClustered, MorphologicalComponents[segs]}
]
findContinuations[i_Image, disk_, idi_, branchPoint_, mcsegs_] :=
Module[{mm, coSegs, segmentsAtBranchPoint, tails, tgs, dests, a, b,
x},
mm = Binarize@ Image[g[{White, Disk[branchPoint, p[disk, 2/5]]},
Background -> Black], ImageSize -> idi];
coSegs = ImageMultiply[Image@mcsegs, mm];
segmentsAtBranchPoint = Select[Union@Flatten[ImageData@coSegs], # != 0 &];
tails = Position[ImageData@coSegs, #] & /@ segmentsAtBranchPoint;
tgs = a /. FindFit[#, a x + b, {a, b}, x] & /@ tails;
dests = Nearest[tgs -> segmentsAtBranchPoint, #, 2][[2]] & /@ tgs;
Sort /@ Transpose[{segmentsAtBranchPoint, dests}] // Union]
getEdges[i_Image, disk_, idi_, mbpClustered_, mcsegs_, vxPos_] :=
Module[{fc, equiv, rules, unified, f, vxsForMask},
fc = findContinuations[i, disk, idi, #, mcsegs] & /@ mbpClustered;
equiv = Flatten /@ Gather[Flatten[fc, 1], Intersection[#1, #2] =!= {} &];
rules = Reverse /@ Thread[First@# -> Rest@#] & /@ IntegerPart@equiv // Flatten;
unified = mcsegs //. rules;
f = Nearest[vxPos -> Automatic];
vxsForMask = Map[f[#, {Infinity, p[disk, 5/4]}] &,
ImageValuePositions[Image@unified, #] & /@ Range@Max@unified, {2}];
Rule @@@ DeleteCases[(Flatten /@ Union /@ vxsForMask /. {x_} :> {x, x}), {}]
]
(*Usage*)
i = Binarize@Import@"https://i.stack.imgur.com/58hg7.png";
i = Binarize@Import@"https://i.stack.imgur.com/qDby8.png";
{idi, disk, vxRules, vxPos} = getProblemParms[i];
masks = getMasks[i, disk, vxPos];
{branchPoints, allEdgeSegments} = collectEdgesForests[masks, disk];
edgesFin = getEdges[i, disk, idi, branchPoints, allEdgeSegments, vxPos];
GraphPlot[edgesFin, VertexCoordinateRules -> vxRules,
MultiedgeStyle -> 1/3, VertexLabeling -> True]