I think your idea of fitting a quadratic makes sense but in your example you are choosing three arbitrary points around the minimum thus not taking advantage of the other data. I suggest using a least-squares fit such as given here http://mathforum.org/library/drmath/view/72047.html That would use all the data points. Since the data curve is quite flat, curve fitting has the potential of giving a better result than trying to do a linear fit of the sides of the curve. In my own experience, I often find that there are one or more clear outlier points. Perhaps an approach that does two curve fits would work better. The first attempt would use all the data and the second attempt would remove data points that are off the curve by some amount. If there are too many divergent points that would be a clue that the process is failing.
Also, the idea of selecting a subset of stars makes sense to reduce the computation. I assume that giving preference to stars near the center of the image would give a better result. Perhaps using some sort of weighting approach where stars near the center contribute more to the result.
Lastly, as was pointed out by another poster, when you start getting donuts the star diameters are clearly increasing. So perhaps a hybrid algorithm that looks at HFR and star diameter would be the way to go. The current algorithm starts off by moving the focuser to one end of the curve as defined by the user. Perhaps a better approach would be to take one image it the current focuser position to establish a baseline. If the star diameters are getting much larger than the baseline you know there is a problem. The current algorithm is predicated on the idea that the focus is near the ideal so why not take advantage of that fact. Moving the focuser in steps away from the current position (rather than moving to the curve edge) would give feedback about whether you are getting into large diameters. Of course, backlash has to be taken into account.