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["https://i.sstatic.net/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]]]]

BSplineCurve" which would be easy enough to fit with enough points. Let me try to edit the question $\endgroup$1 + a Sin[\[Phi]/2 + f]/4 + a Sin[\[Phi]/4 + f]/4is doing, modulating the radius of the spiral to try to fit it the real scroll. I thought about subtracting a simple spiral to the edge and fitting the difference to a function withFindFormulabut don't really know yet how to yet. $\endgroup$