Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

   Norm[#[[2]] - #[[1]]] & /@ Partition[
      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...

share|improve this question
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 Narain Jan 4 at 1:20
In any case, shouldn't you be sorting by Norm[# - {x0,y0}]& instead of just Norm? –  Rahul Narain Jan 4 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 at 1:27
What is pos?? –  s0rce Jan 4 at 3:04
@s0rce Sorry... It's included now. –  Pickett Jan 4 at 3:17

1 Answer 1

up vote 10 down vote accepted

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

Mathematica graphics Mathematica graphics

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]

Mathematica graphics

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

Mathematica graphics

Doesn't look constant to me :(

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

ListPolarPlot[{polar, fitpolar}]

Mathematica graphics

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_}] := 
  fit = LinearModelFit[
           "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]]]

Mathematica graphics

share|improve this answer
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 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 at 23:31
+1 for a deserved Enlightened badge. –  Mr.Wizard Aug 7 at 6:17

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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