function [I,J] = InOutIndices(K,h)

i = round(sqrt(-1));

[x,y] = ndgrid(0:floor(1/h));
z = x + sqrt(-1)*y;
z = z(abs(z)<1/h);

I0 = [];
for k = 1:length(z)
  if sum(ismember(z,z(k)+[0 1 i 1+i])) == 4
    I0 = [I0 z(k)-2-2*i];
  end
end

I = [];
for p = -1:2
  for q = -1:2
    I = unique([I,I0+p+q*i]);
  end
end 
  
J = setdiff(K,I);