Express the elliptical region in equality form
(x-c).'*Q*(x-c)<=1
where c is the centroid, Q is the 2x2 matrix,
Q=R.'*diag(1./[a,b]).^2*R;
and R is a 2x2 rotation matrix.
Then test all pixels to see which satisfy the inequality. This can be vectorized as follows,
[x,y]=ndgrid((1:Nx)-cx,(1:Ny)-cy);
delta=[x(:),y(:)].';
mask= all( sum( (Q*delta).*delta )<=1 );
mask=reshape(mask, Nx,Ny);
