In cases where there are several equally good solutions, algorithm's output is not always the most elegant one. The simplest case is a rectangle: the set of solutions lies on a line and algorithm returns a point that coincides with one of the line's endpoints. Current implementation covers such case when calculating initial best cell, however these lines can occur also in more complex polygons containing parallel edges.

I propose using a cost function when estimating cell quality, which would consider both distance from polygon and distance from centroid in such a way it would prefer points closer to centroid, e.g.:
0.5 * (distFromPoly + distFromPoly / (distFromCtr/distFromPoly+1); distFromPoly >= 0
distFromPoly; distFromPoly < 0

So for example, if pia point is very far away from centroid and its distance from polygon is not much greater than centroid's distance, we would get centroid as optimal point. Weights on both values can be manipulated by cost function. What are optimal set of weights is of course subjective.

Cost function would be then used when comparing two cells in while loop and max attribute of a cell can be then calculated in a following way (java code):
// for distance from polygon take max possible distance
// for distance from mass centre take minimal possible distance
double maxDist = dist + h*SQRT2;
if (maxDist < 0) {
return maxDist;
}
double distFromCtr = Math.max(distance(this, massCentreCell) - h*SQRT2, 0);
return costFunction.compute(distFromCtr, maxDist);

I've implemented this in Java and I'm very satisfied with the results.