function solution = pimg_fit(pimg, psf, StartPoint)
    % pimg_fit(pimg, psf, startpoint)
    %
    % Input:
    %   pimg        N x N   particle image
    %   psf         psf model structure, with the following fields
    %   .kernel     K x K   PSF model kernel
    %   .aa_factor  1 x 1   anti-aliasing factor
    %   .irel_aa    AA*N x AA*N coordinate axis for upsampled image
    %   .jrel_aa    AA*N x AA*N coordinate axis for upsampled image
    %   .n_support  1 x 1   size of support for particle image
    %   startpoint  starting point for optimisation, with fields
    %   .centroid   sub-pixel coordinate of image center
    %   .radius     image radius
    %   .bw         black and white values

    %% create initial guess
    X0          = [StartPoint.centroid(1);
                   StartPoint.centroid(2);
                   StartPoint.radius(1);
                   StartPoint.bw(1);
                   StartPoint.bw(2)];

    X0_lb       = [-2 -2 1 0 0];
    X0_ub       = [+2 +2 psf.n_support inf inf];
    %% solve
    opts        = optimoptions('lsqcurvefit', ...
        'SpecifyObjectiveGradient',true, ...
        'Display','off');
    [X,pimg_resnorm,pimg_res,flag] = ...
        lsqcurvefit(@(X, xdata) pimg_fit_modelfun(X, psf), X0, [], pimg, X0_lb, X0_ub, opts);

    %% format
    solution    = struct('centroid', X(1:2), 'radius', X(3), 'bw', X(4:5), 'flag', flag);
end
