Define an ellipse with a given angle and find if a point lies inside its boundaries

I have a number of points represented by an X a Y and an angle component, and I need to define an elliptical region around them, I currently do that like this:

x1 = 5;
y1 = 8;
radius = 8;
eccentricity = .7;
angle = pi/4;

x2 = x1 + radius * cos(angle);
y2 = y1 + radius * sin(angle);
a = 1/2*sqrt((x2-x1)^2+(y2-y1)^2);
b = a*sqrt(1-eccentricity^2);

and I plot it by finding some elements of that ellipse

t = linspace(0,2*pi);
X = a*cos(t);
Y = b*sin(t);
w = atan2(y2-y1,x2-x1);
x = (x1+x2)/2 + X*cos(w) - Y*sin(w);
y = (y1+y2)/2 + X*sin(w) + Y*cos(w);

here x1 and y1 represent the position of the red dot, and angle represents its angle. Then I define x2 and y2 to be another point at the other side of the ellipse to the desired direction in order to form the ellipse rotated as I want it to.

Is this a good way of doing this or is there a potentially better way of defining an ellipse around a point with a specific direction (note that the eccentricity and other variables are not all that important and we can set them as we wish)?

Another way that worked for me was finding the distance of each point in the ellipse from its center given an angle

a = 5;
e = .7;
theta = -pi : 0.05 : pi;

r = (a*(1-e^2)) ./ (1-e.*(cos(theta)));

x=r.*cos(theta);
y=r.*sin(theta);

so if we plot x and y here we will get an ellipse but it will not be rotated.

Secondly I need to find if some other point with its own x and y position is inside this particular ellipse and I would like to do that without having to use inpolygon that needs a big number of points in the ellipse to work. That is why I would prefer having a formula for the ellipse that passing it another point would return <1 for points inside that ellipse and >1 otherwise.

Answer 1

[Jan]

The most easy approach is defining the ellipse by the two foci F1 and F2. Then an ellipse is the set of points, which sum of distances to F1 and F2 is 2*a, where a is the semi-major axis length. Convert this to the polar form to draw the ellipse, but use this to determine points inside:

P      = [x, y];  % A point in 2D
dist   = norm(P - F1) + norm(P - F2);
inside = (dist <= 2*a);

[Raldi]

I did a few tests but I am not sure I managed to do it correctly.

el_radius = 8;
eccentricity = .7;
angle = pi/4;

x1 = 5;
y1 = 8;
x2 = x1 + el_radius * cos(angle);
y2 = y1 + el_radius * sin(angle);

%center point
xc = (x1 + x2)/2;
yc = (y1 + y2)/2;

fociposxy = [xc yc];

a = 1/2*sqrt((x2-x1)^2+(y2-y1)^2);
b = a*sqrt(1-eccentricity^2);

foci = sqrt(a^2-b^2);
fx1 = xc + foci * cos(angle);
fy1 = yc + foci * sin(angle);
fx2 = xc - foci * cos(angle);
fy2 = yc - foci * sin(angle);

foci1 = [fx1 fy1];
foci2 = [fx2 fy2];

%test points
xq=rand(300, 1)*20;
yq=rand(300, 1)*20;
P=[xq, yq];

norm1 = sqrt(sum((P-foci1).^2, 2));
norm2 = sqrt(sum((P-foci2).^2, 2));
dist = norm1 + norm2;
in = (dist <= 2*a);

%theta = 0:.01:2*pi;
%x=a*sin(theta)+fociposxy(1)-foci;
%y=b*cos(theta)+fociposxy(2);

hold on
plot([x1 x2],[y1 y2],'mo');
plot([fx1 fx2],[fy1 fy2],'bo');
plot(xq(~in), yq(~in), 'b*');
plot(xq(in), yq(in), 'r*');

I failed to correctly draw the ellipse, but I think I am computing the inside points correctly. I feel I am not following your suggestions as I should though since it ended up more complicated I think.