FWHM Focusing. Help!

[Yoddha]

Hello,

I'm working on FWHM focusing implementation, but facing with my lack of mathematical knowledge I'll be very thankful for any hint, idea, link, book name or thought!!!

My research showed that FWHM focusing problem is divided in these steps:

1. Locating a start in the image - easy to solve

2. Calculating the light values of the pixels in the star diameter - easy to solve

3. Finding a Gaussian curve that fits best the light values - not solved yet

4. Calculating the FWHM value - simple calculation knowing the Gaussian curve parameters.

The problematic step is number 3 There is need to find a good curve fitting algorithm (or method) to determine the Gaussian curve parameters that describes best the light values of the half tones that builds the star image. I tried various searches in Google, but without much success. There are few curve fitting methods described in pure mathematical aspects... The common in them is that is needed as starting function.

In short all the information that I found rises two questions:

1. Which is the most suitable curve fitting method to be used in the analysis of astronomical image?

2. What is the starting function that have to be used?

It looks like the answers of these two questions is shrouded quite good (maybe by of my lack of knowledge).

"Brainstorming" the problems, leads me to the though that there should be a special case of the methods because the number of the light values are quite limited 0-255, they are whole numbers and also the count of these values near the the focus should be small (limited by the diameter of the star image).

Can anybody help? Many, many thanks in advance!

[SteveP]

Ivo

I don't know if you're already aware that this technique is also used by a number of other focusing regimes including FocusMax, Maxim's own embedded focuser and @Focus (I think). I guess that Robofocus will also use something based on FWHM as well. This might help when you're tracking down the Gaussian formulae

Steve

[FraserClarke]

What language?

Here's an implementation of fitting a 2d Gaussian courtesy of scipy.org;

Cookbook/FittingData -

[Yoddha]

Steve, thanks for the hints! FocusMax is implementing other techniques which can help when you are quite away from the focus point. But in fact is using the same mathematical instruments that I try to understand MaximDL do not provide details how they implemented the FWHM...

TeaDwarf, thanks for the link! Any programming language will be fine (except Mathlab and similar scripts) if it shows the algorithm which will allow me to implement it in C++. For bad luck the CookBook is using a library named matplotlib to do the main job.

For good luck there is mentioned this link:

gaussfitter.py - agpy - Project Hosting on Google Code

Which I'm inspecting in the moment

[SteveP]

Ivo

Alternatively, have you thought about using an alternative focusing method - Half Flux Diameter? If you're struggling to find a suitable Gaussian distribution to fit a star, the HFD method might be a touch simpler. You calculate the diameter of a circle in pixels such that the flux of the particular star outside that circle is the same as its flux inside the circle. Smaller is better! Apparently it works well even if the star is grossly out of focus

HTH

Steve

[Yoddha]

Thanks Steve! HFD looks good alternative!

[MartinB]

HFD is used by FocusMax.

Don't know if this is any use PixInsight — Deconvolution and HDRWaveletTransform: NGC 5189 from Gemini Observatory South Have a look at the section on deconvolution, PSF

[Yoddha]

Thanks Martin! It is interesting reading!

[Yoddha]

Reading and evaluating the problem again and again leads me to the thought that maybe I'm trying to solve the task by the hard way. It is like proving a theorem rather than just use it. By the theory is known that there is Gaussian curve that describes best the brightness values of the pixels that makes the star image. And again by the theory it is known that Half Maximum will give the focus point.

So what if we just use these facts, rather than prove them and find the mathematical solution of the task in its common case?

What do you think about this simple implementation:

1. Find the brightest pixel(s) in the star image - M

2. Compute the half of this brightness - H = M/2

3. Find the diameter D in the image that is built of pixels with brightness H

4. D is the FWHM value that we are looking for

There will be some complications in the programming, but nothing so hard to do...

More I'm reading about non-linear least squares fitting (regression) methods more I understand that they are complex to calculate and will not be suitable for real-time image processing like FWHM of LiveView stream...

Do you think that there is a mistake in these thoughts?

[roundycat]

earlier you mentioned that the number of light values is quite limited, ie, to 0-255. All the cameras I use are outputting 16bit so the numbers go up to 65,000.

Dennis

[FraserClarke]

Do you think that there is a mistake in these thoughts?

Depends what you are trying to measure...

Your approximation is probably OK if you have relatively well sampled data (FWHM of several pixels at least); below that you'll be limited by the sampling of the data, and need to fit a function to fit.

For a very quick measure, you could also try just plotting the maximum pixel brightness against focus position; that should be correlated with focus (brightest==smallest PSF==best focus). You'll also be sensitive to flux changes (clouds/scintillation) too though.

[Yoddha]

Thanks Dennis!

[Yoddha]

TeaDwarf,

But even we have fitting function, when there is no enough data (not well sampled - the star is saturated, right?) we will have a table-land of the peak value which is not Gaussian curve. The other case is when the start is too faint so the curve will be very small and will cover too few pixels to measure something.

Is it possible the FWHM to be part of the saturated area? If yes, then how to measure it? My understanding is that with or without fitting function the algorithm ends with half-tone (defined from the Gaussian curve) of the maximum and then measures the diameter of intersection of the 2D Gaussian at this level. If at this level we still have pixels with the maximum brightness then we can't measure. Of course it is quite possible to be in big mistake

I'm not sure that understand your idea about "maximum pixel brightness against focus position". You suggest to measure the diameter of the area with the brightest value?

A lot of thanks!

[FraserClarke]

By 'not well sampled' I mean that the Gaussian FWHM is less than 2-3 pixels wide.

Yes, if the data are saturated things get very tricky. It is still possible to fit a Gaussian to the 'wings' (the bits which aren't saturated), but you can't do anything with the saturated area. The 'half-flux' method based on the raw data (without a fit) obviously fails completely if you have a saturated star, because as you say, you don't know what the real maximum value is.

The brightness of the star (assuming it is not saturated) does not affect the size of the FWHM (this is a relatively common misconception because of the way images are displayed). The same Gaussian will fit bright and faint stars. The width will be the same, just the height scaled up for brighter stars. For faint stars you may of course not have enough signal to reliably measure the shape, but the intrinsic shape is the same.

For the "maximum pixel brightness" method, consider what happens when you have an out of focus star. All the light from the star is spread out over more and more pixels; so each individual pixel has less light in it. As you get closer to focus, the star is concentrated into fewer pixels, so the middle pixels get brighter (the total amount of light stays the same of course). When you are in the best focus, the star is most concentrated, so the middle pixels (or pixel) are the brightest they can be. So your best focus is when the middle pixel has the most light in it...

You do have to assume the total amount of light stays the same, which is a problem if you have clouds etc, so it's not a great method -- but it is very quick because all you need to do is find the value of the brightest pixel in the star.

[Yoddha]

TeaDwarf,

Sorry for the delayed answer, but had to think on it on clear mind... I spend two nights in reading tons of formulas and articles.

The "maximum pixel brightness" idea is very clever!

About FWHM... I'm in big hesitation. In one hand we have the definitions

f(x) = exp(- x^2 / 2*(sigma^2) )

and

FWHM = 2.35 * sigma

So in the common case FWHM is not the half of the maximum brightens, right?

In the other hand we have the real world... Here is a star taken with ZS 66SD

It is less than 18 pixels wide. In the best case we will have 18 different pixels with values from 0-65000. MaximDL gives FWHM ~ 10.70 in the three channels for this image. The "half-flux" is 11 pixels by vertical and 12 pixels by horizontal, but it is limited to whole number value.

What is interesting is that MaximDL calculates FWHM very fast. Which makes me think that there should be easy to calculate way... Maybe you are right that there is "fixed" Gaussian that fits the stars model and which is scaled somehow?!?!

[FraserClarke]

About FWHM... I'm in big hesitation. In one hand we have the definitions

f(x) = exp(- x^2 / 2*(sigma^2) )

and

FWHM = 2.35 * sigma

So in the common case FWHM is not the half of the maximum brightens, right?

Sigma represents the 'standard deviation' of the Gaussian function (important in statistics; less so in this case). FWHM stands for 'Full Width at Half Maximum', and is what it says; the full width of the function at the point where the flux (brightness) is half of the maximum. For a Gaussian, there is a relation between the two parameters as you wrote down.

What is interesting is that MaximDL calculates FWHM very fast. Which makes me think that there should be easy to calculate way... Maybe you are right that there is "fixed" Gaussian that fits the stars model and which is scaled somehow?!?!

The Gaussian fitting algorithm can be calculated very quickly on modern computers. It's a pain to code (the usual implementation has some big matrix transformation subroutines), but when it works it will run quickly (<0.01s for a standard image I'd expect). There are examples of gaussian fitting algorithms on the web.

[Yoddha]

Thanks TeaDwarf!

I spent most of the day in reading and calculations. Now I'm sure that Half Maximum is f(max)/2

Half flux has potential, but definitely is not the goal. So more reading and searching is needed to figure out how to make the fitting. Hope that will find explanation how to make the calculations