[From nobody Thu Feb 20 16:39:56 2003
X-Sieve: cmu-sieve 2.0
Return-Path: <ddnebert@fgdc.gov>
Received: from spider.refractions.net (spider [192.168.50.1])
	by lion.animals (8.10.1/8.10.1) with ESMTP id f6CKfNW01579
	for <dblasby@lion.animals>; Thu, 12 Jul 2001 13:41:23 -0700
Received: from prserv.net (out4.prserv.net [32.97.166.34])
	by spider.refractions.net (Postfix) with ESMTP id 34C7E3F13E
	for <dblasby@refractions.net>; Thu, 12 Jul 2001 13:41:01 -0700 (PDT)
Received: from fgdc.gov (slip-32-100-209-93.dc.us.prserv.net[32.100.209.93])
          by prserv.net (out4) with SMTP
          id <2001071220405220405ab0bge>; Thu, 12 Jul 2001 20:40:53 +0000
Message-ID: <3B4E0BA0.BB8BA9D8@fgdc.gov>
Date: Thu, 12 Jul 2001 16:42:08 -0400
From: Doug Nebert <fgdc2a@attglobal.net>
Reply-To: dn@mapcontext.com
Organization: FGDC
X-Mailer: Mozilla 4.77 [en] (Win98; U)
X-Accept-Language: en
MIME-Version: 1.0
To: Dave Blasby <dblasby@refractions.net>
Subject: Points in circle alrogithm
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
Sender: ddnebert@fgdc.gov

Here is an algorithm to extend your existing postgis bounding rectangle
search.

1. Get X,Y coordinate and radius in native cartesian system from user.
2. Calculate MBR for query circle to be:
    minx=X-radius
    maxx=X+radius
    miny=Y-radius
    maxy=Y+radius
3. Perform pre-emptive fetch from postgis of point records in MBR
based on above rectangle dimensions.
4. A speed-enhancing step now, especially useful for very large
areas or dense selections, would be to first calculate the outer MBR
into a selected set, then subselect for the points that fall within a
square that fits within the circle. Points that fall within the inner
square
do not need to be further tested. Only the point records falling
outside the inner square and inside the outer square need to be tested!
This set is defined as:

(((X+D) > A > (X-D))AND((Y+D) > B > (Y-D)) ) AND
(((X+R) > A > (X-R))AND((Y+R) > B > (Y-R)))

5. Move through each selected record to determine if the point
falls within the radius, basically if the following is true, where A,B
represents the returned coordinates found in MBR, and X,Y is
the test centroid and R=Radius posed by the user. If the following
statement is true, then the points fall within the circle:

((A-X)^2 + (B-Y)^2) le (R^2)

where ^2 means squared.

6. The resulting rows would be in or on the perimeter of the circle.
Within test would be a less-than not a less-than-equals operator.

Now, the final trick in my world is that the coordinates are in lat
long, so my circle is not even close to a circle except at the Equator
so I'd need to have my data projected into an equal area projection
before the test, or convert my circle query into an approximate
polygon (or ellipse?) to perform the overlay...

Hope this helps.

Doug.

]