I am closing this bug as fixed. The example in comment #1 behaves correctly in the current version of Octave (4.2.0).

The original example still has problems, but they appear to be solely due to bug #37501, because the input data contains NaN values. If I replace all NaNs in the ZI matrix with 0, then the plot appears correct.

Another example is below

N = 601;
D = round (N/3);
M = round ((N-1) / 10) + 1;
z1 = rand (M, M);
x1 = 10 * ((1:M) - 1) + 1;
[x1, y1] = meshgrid (x1, x1);

[x, y] = meshgrid (1:N, 1:N);
z = 0.5 * (round (2 * interp2 (x1, y1, z1, x, y)) == 2);

r = min (abs (x - D * round (x / D)), abs (y - D * round (y / D)));
z(r < ceil (D/20)) = 1;

contourf (z)
caxis ([0 1])
colorbar