16
$\begingroup$

How can I find the simplest curve that gives me a scroll like this?

enter image description here

This is where I am at the moment and doing it by hand is becoming more and more cumbersome.

First, take a drawing of the scroll since the outline is already processed:

enter image description here

Then, define a spiral in a parametric plot with several parameters that move things around:

Manipulate[
 Show[ParametricPlot[
   1/d (1 + a Sin[ϕ/2 + f]/4 + a Sin[ϕ/4 + f]/4) Norm[
     2 + ϕ^1.1] {Sin[-ϕ], 0.075 + Cos[ϕ]}, {ϕ, 0, 
    5.25 π}, PlotRange -> {-1, 1}],
  Epilog -> Disk[{0, 0}, .075],
  Prolog -> {Texture[img], 
    Polygon[{Scaled[{0, 0}], Scaled[{1, 0}], Scaled[{1, 1}], 
      Scaled[{0, 1}]}, 
     VertexTextureCoordinates -> {{0 + b, c}, {1 + b, c}, {1 + b, 
        1 + c}, {0 + b, 1 + c}}]}], {a, .5, 
  2}, {b, .091}, {c, .148}, {d, 30, 40}, {f, 0, π}]

(b and c are there to center the image only, could not find a nicer way)

What I tried to do was to increase the radius of the spiral as a function of $\phi$ but I have the impression a computer would be smarter than me in finding this. This is the best I can do atm.

enter image description here

I know there is a new FindFormula but cannot figure out how to make it work for parametric functions.

I would very much like something that does not involve BSplineCurve. A piecewise function on $\phi$ is welcome. Violin shapes were made with compass and rulers so I imagine that there must be a simple function for the curvature and how it changes as the scroll turns. With simple function I mean something like this: enter image description here

In this case the radius of each segment on the scroll is given by the projection from the lower triangle. This can be taken to the continuous case as a simple function.

Using the method provided by nikie I can obtain different scrolls varying the functions I use for the fit. This is the one I most like fitFn = a + b theta + c theta^2 + d /theta + e/theta^2 + h Sin[ f (-theta)^.5 + g];

Maybe you can find out to which one it corresponds

enter image description here

$\endgroup$
  • 3
    $\begingroup$ This may be helpful. $\endgroup$ – xslittlegrass Apr 22 '16 at 1:29
  • $\begingroup$ Okay I'm confused: splines are still piecewise functions after all, so you seem to have conflicting desires. Also, "with compass and rulers" can be taken to mean that the arcs could very well be made up of circle arcs and line segments suitably pasted together. $\endgroup$ – J. M. will be back soon Apr 22 '16 at 1:46
  • $\begingroup$ Maybe I should rephrase it as "don't use BSplineCurve" which would be easy enough to fit with enough points. Let me try to edit the question $\endgroup$ – tsuresuregusa Apr 22 '16 at 1:50
  • $\begingroup$ With respect to "simple": the lituus seems to be a good starting point as any for a volute. $\endgroup$ – J. M. will be back soon Apr 22 '16 at 2:08
  • $\begingroup$ @J.M. I tried a few Archimedes' spirals but it doesn't really work, you need to add some kind of "asymmetry" (for lack of a better word), that's what my 1 + a Sin[\[Phi]/2 + f]/4 + a Sin[\[Phi]/4 + f]/4 is 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 with FindFormula but don't really know yet how to yet. $\endgroup$ – tsuresuregusa Apr 22 '16 at 2:38
20
$\begingroup$

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]]]}]]]

enter image description here

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]}]

enter image description here

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]

enter image description here

And that there are no path cuts:

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

enter image description here

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]
    ]}]

enter image description here

(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]]

enter image description here

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

enter image description here

$\endgroup$
  • $\begingroup$ Beautiful. At the moment I;m using this to fit: a + b theta + c Cos[theta + e] + f theta^2 + g Sqrt[Abs[theta]] I don't know either what function was, but I guess it should be easy to obtain with compass and ruler via different projections. Now, which kind of projections are those... that's another question. Thanks a lot! $\endgroup$ – tsuresuregusa Apr 22 '16 at 19:58
  • $\begingroup$ I've added a linear fit. No idea if this can be constructed using ruler and compass. $\endgroup$ – Niki Estner Apr 22 '16 at 20:17
  • 1
    $\begingroup$ This might be useful for taking care of the branch cuts in the angles. $\endgroup$ – J. M. will be back soon Apr 22 '16 at 22:30
  • $\begingroup$ One of the most beautiful things of the scroll is that has a certain "movement". If you plot the difference between your last fit and the curve you obtain a very marked oscillation with a defined frequency during at least 2 cycles. I think that's what make the scrolls alive somehow. I added some figures obtained with your code. Thanks again. $\endgroup$ – tsuresuregusa Apr 22 '16 at 23:37
  • $\begingroup$ I thought the "oscillation" was because I missed the actual center point... $\endgroup$ – Niki Estner Apr 23 '16 at 7:06
13
$\begingroup$

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 manuallyenter image description here

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]

enter image description here

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}]]

enter image description here

$\endgroup$
  • $\begingroup$ thanks a lot, I didn't know about the coordinate tool, quite useful. Is there a way to automatically fix the branches? Will leave the question open for a bit to see if another answer appears but this is kinda exactly what I was looking for, now I will check the coefficients for different scrolls and see how the different "schools" behave. $\endgroup$ – tsuresuregusa Apr 22 '16 at 18:56
  • $\begingroup$ Yes. You could write a function that goes through a list and when it sees a value drop by 2Pi adds 2Pi to all subsequent values. $\endgroup$ – Searke Apr 22 '16 at 19:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.