The idea: take Simplexy's output as the initialization for a better fit.

The setup:

• take a 31x31 pixel cutout around each Simplexy peak.
• initialize a Gaussian + flat background model with Simplexy's output
• use a standard nonlinear optimization routine with various error functions and strategies.

The data:

• thor2.png : 400x300, 131 peaks. Includes diffraction spikes, near-saturation (but no pixels have value 1.0)

Simplexy's output:

Simplexy's error:

These images consider only the cutout areas around each star. The background estimates are averaged in regions where the cutouts overlap, and then the peaks are added.

Another way of looking at it:

Notice that Simplexy's standard deviation estimate is wrong (it's just defined to be 1.0): the error image shows halos around each star. Also notice that the centers and flux estimates for the brightest stars are completely wrong - large black holes are left. Also observe that the background estimate isn't great - the background grey colour is zero error, and many of the cutouts stand out against that background.

In this document, sigma refers to the standard deviation of the Gaussian, not sky noise.

Step A

• fit each peak independently
• standard sum-of-squares of the pixel error values

The fit is better - note that the bright stars now have very little error. However, notice a fews pairs of white and black blobs. These are places where a peak moved from its original position in the center of the cutout to a different star within the cutout, or where it tried to move to compromise between two nearby stars.

Here are the fit sigmas:

Step B

One obvious problem with the above setup is that the model is heavily penalized for having "unexplained flux", so it will tend to migrate to the largest peak within the cutout, or try to merge two or more nearby stars. We want to allow the model to explain only its central star but not allow it to overestimate the flux - we don't want to allow it to invent photons. To do this, we can skew the error curve, penalizing overestimates more than underestimates.

The results, however, are strange. In some cases, the model uses the Gaussian to create an background gradient instead of fitting a peak. Pretty resourceful, actually, but not an improved model for our purposes.

Step C

How can we fix this behaviour? We sort of want to regularize the sigma estimates. For astronomical images, we expect sigma to be pretty much constant across the image - it is determined by the PSF of the telescope. One way of estimating the PSF while also regularizing is to jointly fit the sigma values.

In these results, we alternate between fitting the center, height, and background of each peak independently, and fitting sigma for all peaks simultaneously.

This is a slow method - it might be faster to fit the whole model simultaneously, but for ease of implementation I'm alternating. The result is that sigma grows slowly at each alternation. The results below are after 20 alternations.

(The jumps are the alternations)

The resulting fit is nice:

Except note that many stars have halos - the brightest stars are bright in the middle, while the dimmer ones are dark in the middle. I think this happens because the model wants to increase the peak width for the brightest stars, but the dimmer stars prevent this. It can increase the height of the brightest peaks a little, but not too much since it gets penalized by the high cost of overestimating.

Step D

An important aspect of astronomical images, and one place where Simplexy goes quite wrong, is saturation.

Although this particular image has no actual saturation - no pixel has brightness 1.0 - it looks like the curves have been adjusted to flatten the tops of the brightest stars. This is soft saturation, and a strategy for dealing with real saturation should do fairly well.

A very simple saturation strategy is to simply clamp the model's predictions at 1.0. Here are the results:

The fit is much better, especially for the brightest stars. The saturated cores are almost perfectly fit. There are still slight rings just outside the core, possibly showing the softness of the saturation in this particular image.

Observations

A couple of things to think about:

• a smooth gradient background model might work better. What about having a grid of sample points with bilinear interpolation between them.
• using cutouts isn't great because whenever the cutouts overlap, each model is trying to ignore the other cutout's star - it would be nice to do a more EM-style optimization
• Jon Say: On that note, shouldn't the "flat" backgrounds all be fitted together? Realistically, all single pseudo-flat backgrounds should be fit to the entire image at once. Cutouts are the most problematic here, because either all flat backgrounds are fitted separately, which is wrong, or all flat backgrounds are fitted together, which is also wrong because then each cutout isn't aware of its location on the image. What if all things are fitted together, but just for speed, only pixels in a stars cutout "region" are used in computing the Gaussian. And what if a big bendy surface is fit to all of the rest of the image, like a big bendy Bezier surface.
  • Dstn says: yes. I actually think a piecewise (bi)linear background would work fine - have a grid of sample points and let the background vary linearly between them.
• run simplexy on the residuals!