```Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: how to find pixels inside ellipse ?
Date: Wed, 17 Aug 2011 21:48:09 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 30
Message-ID: <j2hcup\$fh9\$1@newscl01ah.mathworks.com>
References: <g1bsct\$ggh\$1@fred.mathworks.com> <g1c6l2\$b04\$1@fred.mathworks.com> <j2ellm\$gt5\$1@newscl01ah.mathworks.com> <j2eomo\$q1e\$1@newscl01ah.mathworks.com> <j2gopq\$78p\$1@newscl01ah.mathworks.com>
NNTP-Posting-Host: www-06-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1313617689 15913 172.30.248.38 (17 Aug 2011 21:48:09 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Wed, 17 Aug 2011 21:48:09 +0000 (UTC)
Xref: news.mathworks.com comp.soft-sys.matlab:740367

"Cody Oppermann" wrote in message <j2gopq\$78p\$1@newscl01ah.mathworks.com>...
> Thanks Roger, I was planning on using the flat plane assumption as my ellipses are relatively small. Just to clarify though, say I have a point, for example, at 30°,-30° and the center of my ellipse is at, say, 29°,-29°. Would it be most accurate if x-x0 is the distance from 29°,-29° to 29°,-30° (as opposed to 30°,-30° to 30°,-29°) and y-y0 is the distance from 29°,-29° to 30°,-29° (as opposed to 30°,-30° to 29°,-30°)?
> Thanks again.
> 30,-30 •     • 30,-29
> 29,-30 •     • 29,-29
- - - - - - - - -
That's a good question, Cody, but I'm afraid I don't know the answer.  However, in any case I am not satisfied with the method I described to you and am now suggesting an approach less prone to error.

Suppose we define an "ellipse" on the earth's (supposedly) spherical surface as the locus of all points on that surface for which the sum of their great circle arc length distances from two fixed surface focal points is a given constant.  If the location of such foci for your ellipse can be found, it is easy to determine whether points lie within it by comparing such a sum of arc lengths with the given constant.

Transform all latitude, longitude values for points to earth-centered cartesian coordinates using the equations:

x = R*cosd(theta)*cosd(phi)
y = R*cosd(theta)*sind(phi)
z = R*sind(theta)

where R is the earth's radius, theta is the (north-positive) latitude, and phi is the (east-positive) longitude and where theta and phi are in degrees.

Let v1 = [x1,y1,z1] and v2 = [x2,y2,z2] be vectors for the two given foci and let v = [x,y,z] be the vector for an arbitrary point to be tested.  (These are vectors pointing from the earth's center to the surface points.)  Then the sum of the two great circle arc lengths referred to above is:

t = R*(atan2(norm(cross(v,v1)),dot(v,v1)) + ...
atan2(norm(cross(v,v2)),dot(v,v2)));

The point (x,y,z) will lie within the ellipse whenever t is less than or equal to the above given constant.

(Note: You can actually assume that R is equal to one in all the above provided you make the corresponding correction to the assumed arc length constant sum for the ellipse.)

I would venture to guess that, provided you use double precision, the major source of inaccuracy in the above would be the assumption that the earth is a sphere, as opposed to an oblate spheroid.  However the computations for the latter would be much more complicated.

Roger Stafford
```