# Gas station without pumps

## 2012 April 21

### Fitting a sphere

Filed under: Accelerometer,magnetometer,Robotics — gasstationwithoutpumps @ 21:19
Tags: , , , , ,

Today I had to write a program to fit a sphere to a bunch of points that were supposedly near the surface of a sphere, but were noisy and sampled in a very biased way.  Since this is obviously not a new problem, I started out doing web research.  but I didn’t look for fitting a sphere, but for fitting a circle, since that is a simpler related problem.

I found a lot of papers, including several review papers, on how to fit a circle to a bunch of points.  The “obvious” method is to  do a least-squares fit to minimize the distance between the points and the circle, minimizing $\sum_i (r- \sqrt{ (x_i-a)^2+(y_i-b)^2})^2$, where $r$ is the radius and $(a,b)$ is the center of the circle.  Unfortunately, that is a difficult problem to solve, and even numerical methods require a lot of iterations to get decent solutions.  What most people do is to change to a slightly different problem that optimizes a different fitness function.  For example, Kåsa’s method minimizes $\sum_i (r^2 - (x_i-a)^2 - (y_i-b)^2)^2$.

There is a very nice, but very formal, presentation of the methods in a paper by Vaughn Pratt from 1987: Direct Least-Squares Fitting of Algebraic Surfaces.  This paper introduced Pratt’s method, which was later slightly improved to make Taubin’s method. I did not read these original papers (other than skimming Pratt’s paper).  Kåsa’s paper (A curve fitting procedure and its error analysis. IEEE Trans. Inst. Meas. 25: 8–14) does not seem to be available on-line.  The IEEE digital library is missing the whole 1976 year.

I did find a recent paper that does careful error analysis of both the geometric approach and several of the algebraic approaches (including the most popular ones: Kåsa, Pratt, and Taubin):

Error analysis for circle fitting algorithms
Electronic Journal of Statistics
Vol. 3 (2009) 886–911 ISSN: 1935-7524 DOI: 10.1214/09-EJS419

This paper shows that Taubin’s method is theoretically superior to Pratt’s which is theoretically superior to Kåsa’s (having less essential bias), and gives a very weak example showing it is also tru empirically.  More interestingly, it also gives a “hyperaccurate” algorithm that has less bias even than Taubin’s method.  I did not read the error analysis, but I did read the description of their Hyper algorithm and the implementations of it that Chernov has on his website.

Since I needed Python code, not Matlab code, and I needed spheres rather than circles, I spent a few hours today reimplementing Chernov’s Hyperfit algorithm.  I noticed that the basis suggested by Pratt for spheres, $(x^2+y^2+z^2,x,y,z,1)$, was a simple modification of the one used in both Pratt’s paper and Chernov’s paper for circles, $(x^2+y^2,x,y,1)$.  I decided to generalize to $n$ dimensions, and use the Numpy package in Python for all the matrix stuff.  I hope I got the generalization right!

From starting to look for papers until getting the code working was about 6 hours, but I had lunch in there as well, so this felt like pretty speedy development.  I’ve released the code with a Creative Commons Attribution-ShareAlike 3.0 Unported License, and would welcome corrections and improvments to it.

Of course, after all this buildup, you are probably wondering why I needed to fit a sphere to points—that is not a common problem for a bioinformatician to have.  Well, it is for the robotics club, of course.  They’ve been having a lot of trouble with the magnetometer calibration and heading code, so we decided to try doing an external calibration of the magnetometer, which has an enormous arbitrary 3D offset.  By waving the magnetometer around in different orientations (which means tumbling the ROV once the magnetometer is installed), we can sample the magnetic field in many orientations, though far from uniformly.  The center of  the sphere fitted to the readings gives us the 3D offset for the magnetometer.

My son and I tested it out with Python code and Arduino code that he had written to get the data from the magnetometer to the laptop, and the magnetometer readings do seem to be nicely centered around (0,0,0) after we do the correction.  We’re still having trouble using the accelerometer to get a tilt correction to give us clean compass headings, but that is a problem for tomorrow morning, I think.

## 2011 March 20

### Real bioinformaticians write code

Filed under: Uncategorized — gasstationwithoutpumps @ 22:32
Tags: , , , ,

The What You’re Doing Is Rather Desperate blog has a post Real bioinformaticians write code, which makes the point that bioinformaticians have to program.  It will never be the case that every interesting question can be answered by applying a pre-existing tool. I find myself writing code a lot, even when there are tools or packages that supposedly do the job for me.  Quite often, I find that the existing tools only do 80% of the job I need, and adding the extra 20% of functionality is more work than starting over from scratch.  That is not always the case, though, and I do appreciate well-written tools and packages when I can find them.

I used to do most of my programming in large systems I wrote in C++, but lately I find myself doing more rapid prototyping in Python, as it is easier for doing small programs in, and there are some nice general-purpose packages that are easy to use (like NumPy for numerical work).

The post was a response to the frequent form of request on BioStar: “I’m looking for a resource … “, which usually has one of two answers: “Let me Google that for you” or “you’ll have to write a program to do that yourself”.  If you do have a good bioinformatics question, the BioStar website can be a good place to get it answered, though.

One of the comments on the original blog post points to BioMart, an open-source data management system for biologists, for those looking for some of the standard data handling and query stuff.  I’ve never used the tool, as I almost always have questions that require scanning entire data sets, not ones that can be answered by single-point queries, as BioMart is designed for.