Return to Answer

1. find the subset of pixel coordinates for which both the left and right images have foreground (are non-zero);
2. reduce to a map from left value to right value at those coordinates;
3. demote to background anything not in the map; and
4. apply the map.

def main() -> None:
    images = make_samples()
    show(images, 'Before matching')

    both_foreground = np.logical_and.reduce(images, axis=0)
    mapping = np.unique(images[:, both_foreground], axis=1)

    for image, value_set in zip(images, mapping):
        image[~np.isin(image, value_set)] = 0

    # Substitution routine loosely inspired by https://stackoverflow.com/questions/3403973
    k, v = mapping
    sparse_map = np.zeros(v.max() + 1, dtype=np.int32)
    sparse_map[v] = k
    images[1, ...] = sparse_map[images[1, ...]]

    show(images, 'After matching')
    plt.show()

1. find the subset of pixel coordinates for which both the left and right images have foreground (are non-zero);
2. reduce to a map from left value to right value at those coordinates;
3. apply the map; and
4. demote to background anything not in the map.

def main() -> None:
    images = make_samples()
    show(images, 'Before matching')

    both_foreground = np.logical_and.reduce(images, axis=0)
    mapping = np.unique(images[:, both_foreground], axis=1)

    # Substitution routine loosely inspired by https://stackoverflow.com/questions/3403973
    v, k = mapping
    sparse_map = np.zeros(k.max() + 1, dtype=np.int32)
    sparse_map[k] = v
    images[1, ...] = sparse_map[images[1, ...]]
    images[~np.isin(images, mapping[0, :])] = 0

    show(images, 'After matching')
    plt.show()

def main() -> None:
    images = make_samples()
    show(images, 'Before matching')

    both_foreground = np.logical_and.reduce(images, axis=0)
    mapping = np.unique(images[:, both_foreground], axis=1)

    for image, value_set in zip(images, mapping):
        image[~np.isin(image, value_set)] = 0
    for right, left in mapping:
        images[1, images[1,:] == right] = left

    show(images, 'After matching')
    plt.show()

def main() -> None:
    images = make_samples()
    show(images, 'Before matching')

    both_foreground = np.logical_and.reduce(images, axis=0)
    mapping = np.unique(images[:, both_foreground], axis=1)

    for image, value_set in zip(images, mapping):
        image[~np.isin(image, value_set)] = 0 

    # Substitution routine loosely inspired by https://stackoverflow.com/questions/3403973
    k, v = mapping
    sparse_map = np.zeros(v.max() + 1, dtype=np.int32)
    sparse_map[v] = k
    images[1, ...] = sparse_map[images[1, ...]]

    show(images, 'After matching')
    plt.show()

This:

summing (equivalent to and in bool)

is untrue. "And" is equivalent to multiplication.

As for your code:

Simplification

is it possible to reduce the complexity of the matching part of my code ?

Yes, extremely. I think the algorithm may have been overthought. For example, you go through the trouble of building a zero collection but then never use it. This really boils down to

1. find the subset of pixel coordinates for which both the left and right images have foreground (are non-zero);
2. reduce to a map from left value to right value at those coordinates;
3. demote to background anything not in the map; and
4. apply the map.

The results are the same:

def main() -> None:
    images = make_samples()
    show(images, 'Before matching')

    both_foreground = np.logical_and.reduce(images, axis=0)
    mapping = np.unique(images[:, both_foreground], axis=1)

    for image, value_set in zip(images, mapping):
        image[~np.isin(image, value_set)] = 0
    for right, left in mapping:
        images[1, images[1,:] == right] = left

    show(images, 'After matching')
    plt.show()