Fitting Vogel's formula for phyllotaxis to the image of a flower

After having read this question I was inspired to investigate how we could take an image of a plant displaying phyllotaxis and fit Vogel's model of phyllotaxis to the image. The theory is very simple.

The model is this:

$r = c\sqrt{n}$

$\theta = 137.508 n$

And courtesy of shrx we have the original image of the sunflower in that other question, located here. To begin with I wanted to retrieve the positions of the florets in the image, which I did using LocatorPane:

These are the points:

florets = {{182.4,288.},{150.,286.},{166.6,274.5},{181.4,261.},{200.6,274.},{214.6,256.},{234.2,256.},{220.6,279.5},{196.2,247.5},{220.,237.},{238.,235.},{259.2,224.5},{257.8,252.},{261.8,202.5},{284.,193.},{256.8,183.},{274.2,169.5},{257.6,148.5},{247.2,166.5},{233.4,179.},{239.2,196.},{239.4,215.},{224.4,218.5},{223.8,198.},{212.,179.},{222.,163.5},{236.4,152.},{249.4,131.},{250.2,110.5},{268.6,127.},{227.4,82.5},{233.2,102.},{231.,125.},{221.6,142.5},{213.2,123.5},{213.2,102.},{207.4,82.},{196.2,64.},{183.2,84.},{169.4,72.},{194.6,105.},{194.4,122.},{204.,139.5},{203.6,161.},{194.8,178.},{201.4,195.5},{207.2,208.5},{204.2,228.5},{186.2,232.5},{172.,244.},{155.8,258.},{136.,264.},{116.8,259.5},{96.4,252.5},{74.8,246.5},{84.,226.},{102.4,235.},{105.6,215.5},{122.4,221.5},{118.6,239.5},{138.,242.5},{155.4,235.5},{167.6,223.5},{187.8,212.5},{186.4,194.5},{171.6,208.5},{153.2,212.5},{139.8,223.5},{134.,207.},{119.,201.},{120.8,182.},{137.2,185.5},{149.6,193.5},{165.2,193.5},{177.6,178.5},{161.,180.},{165.2,164.5},{184.4,163.5},{188.6,146.5},{178.,136.5},{177.6,114.5},{169.8,97.},{150.2,85.5},{130.4,74.5},{134.2,96.5},{151.8,105.5},{162.,125.},{156.4,140.},{172.4,150.},{149.4,155.5},{145.8,173.5},{130.8,166.5},{124.8,151.},{140.6,140.5},{143.,123.},{121.8,130.},{125.8,111.5},{106.4,110.5},{113.8,89.5},{102.8,129.},{102.4,149.},{110.6,168.},{102.,190.},{91.2,172.5},{80.8,158.},{84.6,136.},{84.8,112.5},{57.4,148.},{64.4,171.5},{48.2,193.5},{65.,210.},{76.2,188.5},{89.8,203.},{123.8,284.}};


Later I will use the fact that $n = (\frac{r}{c})^2$ such that to order the list according to index n, we can use SortBy[florets,Norm].

In order to fit Vogel's model to this data I will try to use the following degrees of freedom:

• The constant c
• The origin {x0, y0}
• Angular offset

And I will use NMinimize. My Mathematica code for this is as follows:

phi[n_, ang_] := ang n
r[n_, c_] := c Sqrt[n]
pos[phi_, r_, {x0_, y0_}, phi0_] := {x0, y0} + {r Cos[(phi + phi0) Degree], r Sin[(phi + phi0) Degree]}

NMinimize[{
Total[
Norm[#[[2]] - #[[1]]] & /@ Partition[
Riffle[
SortBy[Array[pos[phi[#, ang], r[#, c], {x0, y0}, phi0] &, Length@florets], Norm],
SortBy[florets, Norm]
], 2]
],
c > 0, 137 < ang < 138
}, {ang, c, x0, y0, phi0}]


Riffle and the SortBy trick puts the corresponding points with the same index n next to each other. Partition finds these pairs and then I measure the difference between them. When the points are almost on top of each other the difference should be very small, instead I get something quite large:

{8318.44, {ang -> 137.055, c -> 1.21577, x0 -> 169.497, y0 -> 176.741,
phi0 -> 41.2544}}


Now Mean[florets] gives {168.811, 178.544} so the {x0, y0} is probably correct. But why is the fit as bad as it is?

I've worked through this problem using only Mathematica. If the problem could be for example precision in the way that I find the positions, please suggest a different way to do it using Mathematica. What you see is the best I could do...

-
I haven't looked at your code in detail, but using Riffle and SortBy inside NMinimize sounds like a bad idea. It makes your function hard to differentiate, and possible even discontinuous, too. –  Rahul Jan 4 '14 at 1:20
In any case, shouldn't you be sorting by Norm[# - {x0,y0}]& instead of just Norm? –  Rahul Jan 4 '14 at 1:21
@RahulNarain As for discontinuity etc. that's OK as long as it doesn't give me the wrong result. Because NMinimize is very fast, it doesn't appear to have a problem with it. As for sorting, I tried that now and it doesn't change the result. –  Pickett Jan 4 '14 at 1:27
What is pos?? –  s0rce Jan 4 '14 at 3:04
@s0rce Sorry... It's included now. –  Pickett Jan 4 '14 at 3:17

I tried working in polar coordinates and fitting r and theta independently, I fixed x0 and y0 as the mean value. Not a great fit, but not terrible. I'm not sure how much including a free x0 and y0 would improve the fit or if the data simply doesn't match the model perfectly. We can improve the r(n) by adjusting the exponent during the fit instead of fixing it at 1/2 (sqrt), it ends up being close to 0.6.

florets = Rescale[Rescale@florets, {0, 1}, {-1, 1}];
florets = # - Mean@florets & /@ florets;

polar = Transpose@
Reverse@Transpose@
CoordinateTransform["Cartesian" -> "Polar", florets];
polar = SortBy[polar, Last];
const = First@
LinearModelFit[polar[[All, 2]], Sqrt[n], n,
IncludeConstantBasis -> False]["BestFitParameters"]
fit1 = NonlinearModelFit[polar[[All, 2]], a n^b, {{a, 0.5}, b}, n]
ListPlot[{polar[[All, 2]], const Sqrt[Range[Length@florets]]}]
ListPlot[{polar[[All, 2]], fit1 /@ Range[Length@florets]}]


much better.

(*I think something is wrong with the angle model*)

fit = NonlinearModelFit[polar[[All, 1]], {Mod[ang n - offset, 2 Pi, -Pi], -Pi <= offset <= Pi}, {{ang, 1}, {offset, 1}}, n]


angles = Differences[polar[[All, 1]]]
angles = 2 Pi*UnitStep[angles] - angles
Histogram[angles, {0, 2 Pi, 0.2}]


Doesn't look constant to me :(

fitr = fit1 /@ Range[Length@florets];
fittheta = fit /@ Range[Length@florets];
fitpolar = Transpose[{fittheta, fitr}];

ListPolarPlot[{polar, fitpolar}]


Still not great, looks like the angle between each floret is not a constant.

Moving the center position doesn't seem to dramatically improve the fit of r(n)

rdata[{x0_, y0_}] :=
Module[{},
fit = LinearModelFit[
SortBy[(Transpose@
Reverse@Transpose@
CoordinateTransform[
"Cartesian" -> "Polar", # - {x0, y0} & /@ florets])[[All,
2]], Last], Sqrt[n], n, IncludeConstantBasis -> False];
{x0, y0, First@fit["BestFitParameters"], fit["AdjustedRSquared"]}
]
coords = Flatten[
Table[{x, y}, {x, -0.01, 0.01, 0.001}, {y, -0.01, 0.01, 0.001}], 1];
output = rdata /@ coords;
output[[First@First@Position[output[[All, 4]], Max[output[[All, 4]]]]]]
Histogram[output[[All, 4]]]


-
I think it's probably unrealistic to expect the structure of a particular actual flower to match an idealized mathematical model. In fact the linked question already shows clearly enough that some arms of the spiral are of different lengths and/or have different origins (offsets relative to the centre) than the others. –  Oleksandr R. Jan 4 '14 at 6:24
I wasn't confident enough in my code to blame the model, and rightly so since you found errors in it, but now I guess we know. Thank you very much. –  Pickett Jan 4 '14 at 23:31
+1 for a deserved Enlightened badge. –  Mr.Wizard Aug 7 '14 at 6:17