How can I estimate molecule positions based on a density map (similar to a minesweeper map)?

Question

What I want: My goal is to locate molecules based on some measurements. With some maths I managed to extract some kind of density map from the data, that indicates how many molecules contributed to the data measured on each pixel. So where the density map says 3, there are 3 molecules within a certain radius. I want to use this map to find an initial estimate of the molecule locations, so that I can use that as a starting point for further optimalization.

Some things to consider:

• I can estimate the total number of molecules in the image reasonably well, so that variable can be assumed to be known.
• All molecules leave a circular 'footprint' of approximately the same diameter (+- 4 pixels). However, this diameter can vary a bit.
• The density map is noisy.
• A molecule can be everywhere in a pixel.

In this example, there are 16 molecules on a 15x15 grid with a pixel size of 166.25, so that the image is 2493.75 squared.

The density map:

np.array([
    [0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 1, 2, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 1, 1, 2, 2, 1, 2, 0, 0, 0, 2, 3, 3, 0],
    [0, 0, 1, 1, 1, 2, 3, 3, 3, 0, 3, 2, 3, 3, 3],
    [0, 0, 0, 3, 3, 3, 2, 3, 3, 3, 3, 3, 4, 3, 4],
    [0, 0, 3, 3, 3, 3, 2, 3, 4, 4, 5, 4, 4, 4, 3],
    [0, 0, 2, 2, 2, 3, 3, 3, 4, 3, 4, 4, 5, 4, 4],
    [0, 0, 2, 2, 2, 3, 4, 3, 3, 2, 3, 4, 4, 4, 3],
    [0, 0, 1, 2, 2, 3, 3, 3, 3, 2, 3, 4, 3, 5, 4],
    [0, 0, 0, 0, 2, 2, 2, 2, 2, 3, 0, 2, 2, 3, 0],
    [0, 0, 0, 0, 0, 1, 1, 2, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0]
])

The true positions:

np.array([
    [1413.316, 2129.608],
    [1609.36,  1025.322],
    [1167.155, 2048.27 ],
    [ 779.355, 1960.638],
    [1086.153,  617.528],
    [1102.6,   2236.773],
    [1389.079, 1318.307],
    [ 688.23,  2181.688],
    [ 646.378, 2009.764],
    [ 301.399, 1075.606],
    [ 872.145, 1286.817],
    [ 920.3,   1177.433],
    [ 482.789,  622.468],
    [1233.384,  676.559],
    [1216.433, 1594.56 ],
    [2339.047, 2068.726]
])

What I tried: I already tried to do a deconvolution with the impulse response function of the system, but that gave me pure noise. I also tried to do the inverse operation: I drew circles on every pixel with the value given by the density map and added everything together, hoping that that would give some maximums on the molecule locations, but that also didn't result in anything.

The best method that I found so far is using K-means clustering. I spawn n points on random positions in a pixel that has value n on the density map. Next, I do K-means clustering with the number of centroids equalling the number of molecules. As K-means depends a lot on the initialization, I repeat this a number of times and select the best output, however, this still does not give satisfactory results.

This could possibly be a problem that is very suitable to solve with machine learning, however, in this project it's preferable to be able to calculate it with our own mathematics.

Can anyone help me out? Thanks in advance!

Side note: this is my first question here, so please let me know if the question needs any clarifications:)