From: "Rob Comer" <>
Newsgroups: comp.soft-sys.matlab
Subject: Re: how to find pixels inside ellipse ?
Date: Fri, 19 Aug 2011 10:55:11 +0000 (UTC)
Organization: The MathWorks Inc
Lines: 47
Message-ID: <j2lfef$coo$>
References: <g1bsct$ggh$> <g1c6l2$b04$> <j2ellm$gt5$> <j2eomo$q1e$> <j2gopq$78p$> <j2hcup$fh9$> <j2ja2d$ebk$> <j2l4p3$cus$>
Reply-To: "Rob Comer" <>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: 1313751311 13080 (19 Aug 2011 10:55:11 GMT)
NNTP-Posting-Date: Fri, 19 Aug 2011 10:55:11 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1184234
Xref: comp.soft-sys.matlab:740544

"Roger Stafford" wrote in message <j2l4p3$cus$>...
> "Rob Comer" <> wrote in message <j2ja2d$ebk$>...
> > This sounds reasonable, Roger, but I think Cody's problem can be solved more easily (and just as easily on an oblate spheroid as on a sphere) by using the distance function in Mapping Toolbox.
> > 
> > The Mapping Toolbox function ellipse1 defines an ellipse on the surface of a sphere or spheroid by generalizing the parametric form of a planar ellipse:
> > 
> > x = a*cos(theta)
> > y = b*sin(theta)
> > 
> > where the center is a (0,0), the semimajor and semiminor axes are a and b, and the major axis is aligned with the X-axis in the plane. Taking the same approach to this problem, consider an ellipse on a sphere or spheroid defined by its center -- at (lat0,lon0), semimajor and semiminor axes a and b, and a rotation angle between the axes and the local meridian on which the center falls (the meridian for which lon = lon0).
> > 
> > 1. Use the distance function to compute the great circle (geodesic) distance d from the center (lat0, lon0) to a point in question, along with the azimuth angle (the second output of distance -- no need for a separate call to the azimuth function).
> > 
> > 2. Determine theta from the computed azimuth and the rotation angle for the ellipse.
> > 
> > 3. The distance to the elliptical curve, following a geodesic that departs from the center at the computed azimuth, is simply defined as e = hypot(a*cos(theta), b*sin(theta)).
> > 
> > 4. If d < e, the point is inside the ellipse. If d == e the point is on the ellipse. If d > e the point is outside the ellipse.
> > 
> > Hope this helps,
> > 
> > Rob Comer
> > Mapping Toolbox Development
> > MathWorks
> - - - - - - - - -
> Hello Rob,
>   If I understand your parametric method correctly, it is equivalent to generalizing the polar coordinate form of a planar ellipse.  The cartesian coordinate equation of an ellipse
>  x^2/a^2 + y^2/b^2 = 1
> has the polar coordinate form
>  r = 1/sqrt((cos(phi)/a)^2+(sin(phi)/b)^2)
> where x = r*cos(phi), y = r*sin(phi).
>   As applied to a spherical or ellipsoidal surface this latter equation can be used to generalize r to a geodesic distance from the ellipse center to a point on the ellipse and phi to the angle between that geodesic and the major (or minor) axis at the center of the ellipse.  As you say, the Mapping Toolbox can readily be used for this purpose.

Roger, I agree that this is the case at least figuratively. One does need take into account that the polar angle and parametric angle are not equivalent except when the ellipse is a circle (just as the geocentric and parametric latitudes on an oblate spheroid are not the same). -- Rob

>   I agree with Bruno that this definition of an ellipse may well differ from one involving a constant sum of distances from two fixed foci as applied to a spherical or ellipsoidal surface.

I also agree with Bruno on that. Of course, for small ellipses (relative to the size of the planet) they are similar. And the definition used in ellipse1 is pragmatic and relatively easy to compute. The main application of which I'm aware is to indicate the degree of uncertainty with respect to the location of a point on the surface -- an earthquake epicenter, for example. Here a and b relate to the eigenvalues of a covariance matrix, and everything is somewhat approximate to begin with. (Cody, what's your application?)

Intersecting a cone with the surface, as Bruno suggested, is a legitimate and interesting problem also. I'd have to study it before declaring the answer to be an ellipse. Certainly the cone can miss the planet entirely, resulting in the empty set, or partially, resulting in a curve that has sharp cusps where it needs to "shortcut" back along a great circle.
-- Rob