function [Ix,Iy] = gradientKeepDirectionalMax(Ix,Iy)

G = (Ix.*Ix+Iy.*Iy).^0.5;
M = G*0;
[h,w] = size(Ix);
theta = atan2(Iy, Ix);
tx = (2:w-1);
ty = (2:h-1);
[X,Y] = meshgrid(tx,ty);
X1 = X+sqrt(2)*cos(theta(ty,tx));
Y1 = Y+sqrt(2)*sin(theta(ty,tx));
G1 = interp2(X,Y,G(ty,tx),X1,Y1);
X2 = X-sqrt(2)*cos(theta(ty,tx));
Y2 = Y-sqrt(2)*sin(theta(ty,tx));
G2 = interp2(X,Y,G(ty,tx),X2,Y2);
M(ty,tx) = (G(ty,tx)>G1).*(G(ty,tx)>G2);
Ix = Ix.*M;
Iy = Iy.*M;
end

