See links at the end for documentation links of Mathematica functions used in this answer. One possibility might be to copy paste the current page to a separate tab and scroll to the bottom there to check documentation when needed.
I will focus here on the clumps of points that make it hard to obtain the markers.
Edit
The previous version was slow and did not provide a tight packing of disks in the clumps.
The idea here will be to choose a circle that is close to the boundary of the clump using Erosion[img,DiskMatrix@radius]
where radius is the typical radius of disk in the image. Then we will obtain the curve from the eroded contour with ImageMeasurements
and then choose a random point on this contour which will have a radius distance with the original boundary. Then we define the new region to be the previous region minus the disk found. We repeat the process above recursively.
First we upscale the image to have less pixeled circles later:
Note : Image@ImageGraphics changes the size of the images in certain versions such as mathematica 13.1. Hence the scales given in the other answer would have to be re computed or adjusted to the image sizes here.
img3 = Image@ImageGraphics[img2, Method -> "Exact"];
Second, we remove the components that are completely attached to clumps. We need to separate the little links they have with one another. This is done using Erosion.
HighlightImage[img3, {Magenta, ColorNegate@img4, Blue,
Erosion[ColorNegate@img4, 3]}]
The image below shows the original image, the binarized image in red and the eroded binarized image in blue. Notice that the eroded regions are better separated than in the binarized image.

We can then extract just the clumps using:
img5 = img4 // ColorNegate // Erosion[#, 3] & //
SelectComponents[#, Large] & // Dilation[#, 3] & // ColorNegate;

Now we erode the clump and center a disk at a random part of the boundary of the eroded region. I chose a radius of 5.6 a bit arbitrarily you might want to adjust that.
contours =
img5 // ColorNegate // Erosion[#, DiskMatrix[5.6]] & //
ImageMeasurements[#, "Contours"] &;
point = contours // RegionUnion // RandomPoint;
disk = Disk[point, 5.6];
Next we make a mask from disk. As I understand, recent versions of Mathematica have a bug that makes it difficult to convert regions to image while preserving coordinates. This work around worked for me. You should check if it works for you.
mask = Binarize@ImageResize[#, ImageDimensions[img5]] &@
ImageCrop@
Image@Graphics@
RegionDifference[
BoundaryDiscretizeRegion@
Rectangle[{0, 0}, ImageDimensions@img5],
BoundaryDiscretizeRegion@disk]
HighlightImage[img5, {contours, Blue, disk, Green, mask}]
The image below shows the contour of the eroded region and the disk
centered at a point of the boundary of the eroded region. Notice that
the circle is very close to the boundary of the original image. You should check whether the blue and green overlap which implies that the mask overlaps well with the disk position.

We now remove the disk from the original img5:
img5 = ImageDifference[img5, mask]

We then repeat the process. The function below repeats the process above:
getDisks[img_, radius_, iterations_] :=
Module[{contours, point, disk,
diskList, imgaux, scale, mask},
imgaux = Binarize@img;
diskList = {};
Do[
contours =
imgaux // ColorNegate // Erosion[#, DiskMatrix[radius]] & //
ImageMeasurements[#, "Contours"] &;
point = contours // RegionUnion // RandomPoint;
disk = Disk[point, radius];
diskList = {\[FormalW] @@ point, diskList};
mask =
Binarize@ImageResize[#, ImageDimensions[img]] &@
ImageCrop@
Image@Graphics@
RegionDifference[
BoundaryDiscretizeRegion@
Rectangle[{0, 0}, ImageDimensions@img],
BoundaryDiscretizeRegion@disk];
imgaux = ImageDifference[imgaux, mask];
, iterations];
Flatten[diskList, Infinity] /. \[FormalW] -> List //
Map[Disk[#, radius] &]
]
Test and example
(NI decreased the disk radius because I was trying to fit more disks. You might need to adjust the radius to your needs)
HighlightImage[img5, getDisks[img5, 5.4, 51]]

If time is not much of an issue you can also consider using getDisks on the entire image (other than the labels) (img5 will be defined in the preprocessing step below):
HighlightImage[img5, getDisks[img5, 4, 100]]
100 disks (you can probably ask for more but it takes a while if time is important it might be best to use the code above only on the clumps )

The preprocessing:
The cropped image which is set to img2

vectorize the image:
img3 = Image@ImageGraphics[img2, Method -> "Exact"]
remove noise :
img4 = ColorNegate@DeleteSmallComponents[#, 40] &@
ColorNegate@Binarize@img3
remove horizontal lines :
img5 = Erosion[#, 3] &@Dilation[img4, 3]
Previous version
The idea is to randomly pack markers within that region. Getting markers to fit within the region is non trivial given the shape of the clumps. Here is an attempt:
First we can focus on the first clump by cropping the image:

That image is called img2
in the following. To extract the clump we could maybe use :
mesh = ImageMesh@ColorNegate@Binarize@img2
Without any parameters, Binarize
will remove the light gray markers. You can use something like Binarize[img2, 0.9]
to get some or all of the gray markers (I did not check if all are kept). In the following, the gray markers are discarded.
Then one may obtain disks that are roughly/approximately contained within the region:
r = 2 (* Radius of disk as an example. I did not check what the radius should be *);
circlearea = π*r^2
disks = Select[(1/circlearea)
RegionMeasure@
RegionIntersection[BoundaryDiscretizeRegion@#, mesh] > 0.9 &]@
Thread@Disk[RandomPoint[mesh, 10000], 2];
nooverlaps =
DeleteDuplicates[disks,
Not[RegionIntersection[#1, #2] === EmptyRegion[2]] &];
visualization of the markers:
Graphics[{Red, Opacity[0.2], mesh}~Join~{Blue}~Join~nooverlaps,
Background -> White]

With 10000 markers it takes a lot of time and the region was not filled with the random points I got.
Links below are generated automatically using Mathematica on the text of this answer . May contain errors .
{Erosion,DiskMatrix,ImageMeasurements,First,Image,ImageGraphics,Method,HighlightImage,Magenta,ColorNegate,Blue,SelectComponents,Large,Dilation,Now,Contours,RegionUnion,RandomPoint,Disk,Next,Binarize,ImageResize,ImageDimensions,ImageCrop,Graphics,RegionDifference,BoundaryDiscretizeRegion,Rectangle,Green,ImageDifference,Module,Do,[[FormalW]](https://reference.wolfram.com/language/ref/\[FormalW].html),Flatten,Infinity,List,Map,If,DeleteSmallComponents,Previous,Here,ImageMesh,In,Select,RegionMeasure,RegionIntersection,Thread,DeleteDuplicates,Not,EmptyRegion,Red,Opacity,Join,Background,White,With,Links}
ImageGraphics[img, Method -> "DualMarchingSquares"]
but the online options might look nicer. You can choose other methods if you prefer. $\endgroup$