# Describe overlapping parts of an image with the fewest points

Lets say I have an image composed of 19x19 pixel red squares, each of which represents a single point (i.e. a plotmarker). I wish to obtain an estimate for the actual point each square represents, assuming that there are as little overlapping points as possible. This can be done by hand, to get the corresponding dots (second image, with a line indicating one spot where a single square could be in several different locations).

I've been wondering if there's a way to do this automatically, but using something like:

HitMissTransform[Binarize[img], ConstantArray[1, {19, 19}]]
(*or*)
Erosion[Binarize[img], 9]


ends up giving wide overlapping regions (seen in the third image). There should be a better way of doing it using corners/edges to more effect, but I can't seem to figure it out.

For square markers of known size, this is actually straightforward with morphological transforms. Using:

img = Import["https://i.stack.imgur.com/nLlB5.png"];
bin = Binarize[img];


We can use HitMissTransform with a 2x2 corner template to find the corners of each square, e.g.

HighlightImage[bin, HitMissTransform[bin, {{-1, 1}, {1, 1}}]]


for the bottom right corner:

If we find each of the 4 corners and translate them to the box centers, we get most of the points:

centers = ImageApply[Max, {
ImageTransformation[HitMissTransform[bin, {{-1, 1}, {1, 1}}],
ImageTransformation[HitMissTransform[bin, {{1, -1}, {1, 1}}],
ImageTransformation[HitMissTransform[bin, {{1, 1}, {-1, 1}}],
ImageTransformation[HitMissTransform[bin, {{1, 1}, {1, -1}}],
}]


But for some squares, all 4 corners are occluded. These are missing now:

missing = ColorNegate[bin] - Dilation[centers, radius]


To find these centers, I'll use HitMissTransform to find the top, left, right, bottom edges of each box, choose only the "missing" ones and translate them to the box centers:

centers2 = ImageApply[Max, {
ImageTransformation[
HitMissTransform[bin, {{-1, -1}, {1, 1}}]*missing,