A comment from a referee led me to find a method to describe curvature of membranes. This is a quick write-up of our solution.
I couldn’t find a solution readily available in Fiji, so I created one using a combination of Fiji for tracing the curvature and IgorPro to do the fitting. If there is any interest I can convert this code into an all-Fiji routine. The code is simple and a gist is available here.
The method
We use a least-squares circle fit, and here’s a short technical description. Click here to skip if you just want to see the code in action.
A curve that represents the membrane is a finite set of points (xy coordinates), how can we find the circle that best fits this set of points.
In mathematical terms, we can define the centre of our set, by taking the average of x and y coordinates.
\(\bar{x} = \frac{1}{n}\sum_{i} x_i\) and \(\bar{y} = \frac{1}{n}\sum_{i} y_i\)
Now, if we use this centre to offset the coordinates,
Now, if we use this centre to offset the coordinates, \(u_{i} = x_{i} – \bar{x}\), \(v_{i} = y_{i} – \bar{y}\)
we can solve the problem in \((u,v)\) coordinate space and then transform back to \((x,y)\).
We can say that the circle we want to find has a centre \((u_{c},v_{c})\) and radius \(R\).
Now, we want to minimise \(S = \sum_{i} (g(u_{i},v_{i}))^2\), this is the “least-squares” bit. Here \(g(u,v) = (u – u_{c})^2 + (v – v_{c})^2 – \alpha\), and \(\alpha = R^2\)
We can minimise \(S\) by differentiating \(S(\alpha,u_{c},v_{c})\).
\(\frac{\delta S}{\delta \alpha} = -2 \sum_{i} g(u_{i},v_{i})\)and this would be zero if \(\sum_{i} g(u_{i},v_{i}) = 0\)
We can continue… \(\frac{\delta S}{\delta u_{c}} = -4 \sum_{i} u_{i} g(u_{i},v_{i}) + 4 u_{c} \sum_{i} g(u_{i},v_{i})\)
note that this final sum term will equal 0 as described above.
So, when \(\sum_{i} g(u_{i},v_{i}) = 0\) then \(\frac{\delta S}{\delta u_{c}} = 0\) if \(\sum_{i} u_{i} g(u_{i},v_{i}) = 0\)
and it follows that if we need \(\frac{\delta S}{\delta v_{c}} = 0\) then \(\sum_{i} v_{i} g(u_{i},v_{i}) = 0\)
So far, so hopefully not too confusing.
We can expand these equations (the u side is shown here)
\(\sum_{i} u_{i} (u_{i}^2 – 2 u_{i} u_{c} + u_{c}^2 + v_{i}^2 – 2 v_{i} v_{c} + v_{c}^2 – \alpha) = 0 \)and using this kind of notation for simplicity… \(S_{uu} = \sum_{i} u_{i}^2\) etc.
we can rewrite the above equation as
\(S_{uuu} – 2 u_{c} S_{uu} + u_{c}^2 S_{u} + S_{vvv} – 2 v_{c} S_{uv} + v_{c}^2S_{u} – \alpha S_{u} = 0\)and since \(S_{u} = 0\), we can simplify this to:
\(u_{c}S_{uu} + v_{c}S_{uv} = \frac{1}{2}(S_{uuu} + S{uvv})\)The v side can be done similarly so that we can use
\(S_{v} = 0\) to give
\(u_{c}S_{uv} + v_{c}S_{vv} = \frac{1}{2}(S_{vvv} + S{vuu})\)We can solve these two equations simultaneously to give
\((u_{c},v_{c})\). Then we can find the centre of the fitted circle in the original coordinate system \((x_{c},y_{c})\) by \((u_{c},v_{c}) + (\bar{x},\bar{y})\).
Nearly there! We just need to find the radius of the fitted circle. Making use of the expanded equation above and using
\(S_{u} = S_{v} = 0\) we get \(\alpha = u_{c}^2 + v_{c}^2 + \frac{S_{uu} + S_{vv}}{N}\)
since we know that \(R = \sqrt{\alpha}\) we can get the radius out.
While that looked kind of complicated, it is actually quite straightforward to write a function to compute the centre and radius for the fitted circle. If you follow the gist you can see that once all the sums (the Suu… and so on) are calculated, we simply construct two matrices and use a linear solving function to find the centre. From there the radius follows.
Some results
It’s always a good idea to test out your solution (and indeed your implementation of it) with some examples.
The code performs well with the simple case of three points. Here at 1,3 2.5,4.5 and 4,3
The points do not have to be uniformly distributed along the circumference, for the fitting to work.
Nor do they need to be a complete circle. Here just an arc of 0.3 pi radians is used to generate some random points on the circumference of a circle.
This method requires points to be on (or approximating to) the circumference of a circle. Here random points within a circle are used and the method fits a smaller circle to the points rather than the “correct” solution.
This is worth doing to know the limitations of our function.
So what about a real example?
Using the freehand line tool in Fiji, a contour on an EM image was drawn and imported into Igor for fitting. The line on the left is a series of points and the resulting circle is fitted to this curve (grey, right).
In this example (in pixel coordinates) of the original EM image we find the radius of curvature of this profile of interest. This result can be scaled back to real world units.
One final snag…
The fitting will produce a circle and, as long as the points approximate a curve on some scale, the fit should be reasonable. However, we may be measuring positive or negative curvature, with respect to some other object. Our simple method does not “know” which side of the membrane is which. I dealt with this by manually annotated the curves and reading that information in Igor. A future solution would be to place an ROI on one side of the curve and compute whether the curvature was positive or negative with respect to this ROI.
—
The post title is taken from “In Circles” by T.2. from their album “It’ll All Work Out In Boomland”. I’m amazed I haven’t used this track before as a title. I love this record. Reissued in the early 90s, it’s probably one of my favourite prog rock albums. I bought a copy on vinyl around this time and sold it like an idiot in the 00s. I made the mistake of looking it up on discogs just now and realised that the copy I had was likely an original Decca pressing…