function [IJ,EJ] = NearestInner(I,J,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

r = -1:2;
r = [r-i r r+i r+2*i];

for j = 1:length(J)
  [m,f] = min(abs(J(j)-(I0+1/2+i/2)));
  f = f(1);
  IJ(j,:) = I0(f)+r;
end
 
for j=1:length(J)
  I = IJ(j,:);
  i1 = unique(real(I));
  i2 = unique(imag(I));
  for i=1:16
    v1 = i1==real(I(i));
    v2 = i2==imag(I(i));
    p1 = polyfit(i1,v1,3);
    p2 = polyfit(i2,v2,3);
    EJ(j,i) = round(polyval(p1,real(J(j)))*polyval(p2,imag(J(j))));
  end
end