Fitting Stradivari's scroll

Step 1: find the curve in the black and white image. For this, I will use a shortest path search, and to make sure it finds the path I'm looking for (instead of one of the short cuts through the labels), I will give it a few "path markers" along the way that the path should visit in order:

First, load the image, get the black pixels:

img = Import["http://i.stack.imgur.com/XBEkk.png"];    
blackPts = PixelValuePositions[Binarize[img], 0];    
nearestFn = Nearest[blackPts -> Range[Length[blackPts]]];    

Next, enter rough path markers using a LocatorPane

pathMarkers = {{373, 5}, {290, 595}, {475, 588}, {495, 371}, {382, 476}};    
LocatorPane[Dynamic[pathMarkers], 
 Dynamic@Show[img, 
   Graphics[{Red, 
     Line[blackPts[[First /@ nearestFn /@ pathMarkers]]]}]]]

Then, turn the black pixels into a graph, with point distances as edge weight:

vertices = Range[Length[blackPts]];
edges = Select[
   Flatten[Thread /@ Thread[vertices <-> nearestFn[blackPts, 20]]],
   #[[1]] < #[[2]] &];
edgeWeights = 
  edges /. {i_ <-> j_ :> Norm[N@blackPts[[i]] - blackPts[[j]]] - .1};
g = Graph[vertices, edges, EdgeWeight -> edgeWeights];

(Don't try to evaluate g in a notebook. It's a very large graph, displaying it crashed the kernel on my machine. But we don't have to display it, we just want to find paths in it, which is very fast.)

Then find the shortest path that touches out path markers:

pathPoints = 
  Flatten[FindShortestPath[g, ##] & @@@ 
    Partition[First /@ nearestFn /@ pathMarkers, 2, 1]];
path = blackPts[[pathPoints]];
HighlightImage[img, {Thick, Line[path]}]

That looks promising.

Step 2: transform to polar coordinates.

I find handling polar coordinates easier working in the complex plane:

center = {402, 465};
pathComplex = N[(# - center & /@ path).{1, I}];

The radii are obvious:

radii = Abs[pathComplex];

The branch cut takes a bit more effort:

first, multiply each point with the conjugate of the next point:

pointPairs = 
  MapThread[#2*Conjugate[#1] &, {pathComplex[[;; -2]], 
    pathComplex[[2 ;;]]}];

Now, the phase of pointPairs is the angle between each point pairs. And since the points are close to each other, this will not contain any discontinuities due to the branch cut:

anglesBetweenPointPairs = Im@Log[pointPairs];

I just have to sum these angles up to get the original phase again:

angles = Prepend[Accumulate[anglesBetweenPointPairs], 0] + 
   Im[Log[pathComplex[[1]]]];

We can check that it's the same path:

ListPolarPlot[{angles, radii}\[Transpose], Joined -> True]

And that there are no path cuts:

ListLinePlot[{angles, radii}\[Transpose]]

You can even move the center around and watch as the polar curve changes dynamically to see where the best center point might be:

center = {402, 465};
Row[{
  LocatorPane[Dynamic[center], Image[img, ImageSize -> All]],
  Dynamic@
   Module[{pathComplex, radii, pointPairs, anglesBetweenPointPairs, 
     angles},
    pathComplex = N[(# - center & /@ path).{1, I}];
    radii = Abs[pathComplex];
    pointPairs = 
     MapThread[#2*Conjugate[#1] &, {pathComplex[[;; -2]], 
       pathComplex[[2 ;;]]}];
    anglesBetweenPointPairs = Im@Log[pointPairs];
    angles = 
     Prepend[Accumulate[anglesBetweenPointPairs], 0] + 
      Im[Log[pathComplex[[1]]]];
    ListLinePlot[{angles, radii}\[Transpose], ImageSize -> 500]
    ]}]

(I have no idea what kind of function might have been used for the radius/angle relationship, so I'm not going to attempt the actual curve fitting. IMHO, if you don't know anything about a functional relationship, a spline usually is the best guess.)

FWIW, I get a somewhat decent approximation using these terms:

Clear[g, eq]
fitFn = a + b theta + c theta^2 + d/theta + e/theta^2;
eq[theta_] = 
  fitFn /. FindFit[{angles, radii}\[Transpose], 
    fitFn, {a, b, c, d, e, f, g, h}, theta];

Show[
 ListLinePlot[{angles, radii}\[Transpose], ImageSize -> 500],
 Plot[eq[theta], {theta, Min[angles], Max[angles]}, 
  PlotStyle -> Directive[Dashed, Red], PlotRange -> All]]

HighlightImage[img,
 Line[Array[center + ReIm[Exp[#*I + Log[eq[#]]]] &, 500, 
   MinMax[angles]]]]

I didn't get a good answer, but this is how I would do it and probably how I would evaluate the quality of an answer. Much of what I do below involves manual work, but it could be automated:

First, we need some points along the curve. I selected them manually

If I hadn't manually selected them, I would have used some image processing including a thinning transform and then I would have used FindCurvePath.

We also want to select a point that seems like the "center" of the spiral.

center = {{165.92578125`, 147.453125`}};
pts = {{183.65234375`, 292.625`}, {201.4765625`, 
   287.3671875`}, {215.21484375`, 282.234375`}, {226.21484375`, 
   277.7578125`},.....};

And we want to normalize the pts so the "center" is the origin:

normalizedPts = (# - First@center) & /@ pts;

We want to work with this in polar coordinates:

ToPolarCoordinates@normalizedPts

But there's a problem. There's a branch cut. I've manually edited the results so that there isn't a branch cut. This involves subtracting or adding 2Pi to the results so that they form a continuous curve:

arc1 = {{146.25014293696273`, 
   1.4492904217728535`}, .... {82.16481360563434`, -3.113682087252578`}}

arc2 = {{79.36095240482918`, 2.9819303345670916`},... {39.305506825650895`, -2.883525354070561`}}

arc3 = {{36.32153376025393`, 3.0092472204025533`}, {19.075687437056864`, -2.8414075099042826`}}

branchCutsRemoved = Reverse/@ 
 Join[arc1, (# - 2 Pi) & /@ arc2, (# - 4 Pi) & /@ arc3];

I reversed all the points so that they're in {arg, abs} form. This means we are looking for the absolute value of the point as a function of arg. I use theta for arg or the angle.

We now have a fairly reasonable function of theta to try and approximate:

ListPlot[branchCutsRemoved, PlotRange -> All]

You can try to fit to this:

eq = Normal@
  NonlinearModelFit[branchCutsRemoved, 
   a + b theta + c  Cos[d theta + e] + f  Cos[d theta/2 + e], {a, b, 
    c, d, e, f, g, h}, theta]

The fit seems... okayish:

Show[ListPlot[sorted, PlotStyle -> Red], Plot[eq, {theta, -16, 5}]]