Solver to fit a 2D gaussian distribution

See more:

optimization

C#6.0

Hi,

I am using C# and the Solver to fit a 2D Gaussian. Using simulated data (no noise) with various sigmas, intensities and center, it was working perfectly. However, with real data, it works only about 50% of time even though all the test data are images of about the same quality and approx. same sigmas for both x & y. Can I submit the code so someone who has more experience can see what I did wrong.

Images are about 100x100 with 5 mm resolution and sigmas of about 7mm for both x & y.

When failed, it failed with values of sigmas and intensity both close to 0.0.

It is such a simple problem, I don't think the Solver can't handle it and the data are not that noisy at all. It must be something subtle that I didn't construct the model properly.

I would greatly appreciate if someone can help.

Thanks,

Daniel

Here is the main code

var solver = SolverContext.GetContext();
            solver.ClearModel();
            var model = solver.CreateModel();
                 
            var decisionsA = Spot_Params.Spot_ParamsA.Select(
                it => new Decision(Domain.RealNonnegative, it.Key));
            model.AddDecisions(decisionsA.ToArray());
            var decisionsB = Spot_Params.Spot_ParamsB.Select(
                it => new Decision(Domain.Real, it.Key));
            model.AddDecisions(decisionsB.ToArray());
          
            var objective = new SumTermBuilder(gaussMatrix.Data_1D.Count());
            var sx = model.Decisions.ElementAt(0);           
            var sy = model.Decisions.ElementAt(1);           
            var sf = model.Decisions.ElementAt(2);            
            var cx = model.Decisions.ElementAt(3);            
            var cy = model.Decisions.ElementAt(4);           

            var sx_sq_2_inv = 1.0 / (sx * sx * 2.0);
            var sy_sq_2_inv = 1.0 / (sy * sy * 2.0);
                      
            int n = 0;
            for (int j = 0; j < gaussMatrix.Nrow; j++)
            {
                var dy = gaussMatrix.Ypos[j] - cy;

                for (int i = 0; i < gaussMatrix.Ncol; i++)
                {                              
                    var dx = gaussMatrix.Xpos[i] - cx;
                    var v = (dx * dx) * sx_sq_2_inv + (dy * dy) * sy_sq_2_inv;                  
                    var value = sf * Model.Exp(-v);
                    var diff = gaussMatrix.Data_1D[n++] - value;

                    objective.Add(diff * diff);
                }
            }

            model.AddGoal("OptimizedSpotParams", GoalKind.Minimize, objective.ToTerm());
                      
            // Extract fitted parameters
          
            var spot_params = fittedParams.Spot_Params;
            Solution solution;
            solution = await Task.Run(() => solver.Solve());

What I have tried:

Reviewed code and input data. Couldn't see anything obviously wrong