Union using polyshape and polybool

I used the voronoin function and obtained the voronoi cell and the vertices. I then found the voronoi cell areas using the polyshape and polyarea functions.

[vor_vert,vor_cell] = voronoin(centroids);

I know certain voronoi cells intersect as I check, using ismember, whether they call the same vertex. This is done within a for loop but I'll only post the conditional if statement:

sum(ismember(vor_cell{jj},vor_cell{ii})) > 0

I then use:

polyout = union(polyshape1,polyshape2)

where polyshape1 and polyshape2 are polygons of two known connected voronoi cells. However, sometimes I do not get a polyout with one region. Instead, I get two separate regions.

If I use polybool however:

[xi,yi] = polybool('union'...)
polyout2 = polyshape(xi,yi)

polyout2 is now a polyshape with one region. This is what I would expect.

Can anyone explain why these two different scenario's have occurred? I am aware that in a certain minority of cases, there may be one and only one point where the two polygons intersect. But there are many cases where this cannot be and yet "union" will return two region polygon and polybool will behave as expected...

I am curious because there now exists a recommendation to use polyshape so I would like to avoid using polybool if possible.