1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
|
/**********************************************************************
* $Id: NonRobustCGAlgorithms.cpp,v 1.14 2004/07/08 19:34:49 strk Exp $
*
* GEOS - Geometry Engine Open Source
* http://geos.refractions.net
*
* Copyright (C) 2001-2002 Vivid Solutions Inc.
*
* This is free software; you can redistribute and/or modify it under
* the terms of the GNU Lesser General Public Licence as published
* by the Free Software Foundation.
* See the COPYING file for more information.
*
**********************************************************************
* $Log: NonRobustCGAlgorithms.cpp,v $
* Revision 1.14 2004/07/08 19:34:49 strk
* Mirrored JTS interface of CoordinateSequence, factory and
* default implementations.
* Added DefaultCoordinateSequenceFactory::instance() function.
*
* Revision 1.13 2004/07/02 13:28:26 strk
* Fixed all #include lines to reflect headers layout change.
* Added client application build tips in README.
*
* Revision 1.12 2004/03/17 02:00:33 ybychkov
* "Algorithm" upgraded to JTS 1.4
*
* Revision 1.11 2003/11/07 01:23:42 pramsey
* Add standard CVS headers licence notices and copyrights to all cpp and h
* files.
*
*
**********************************************************************/
#include <geos/geosAlgorithm.h>
#include <stdio.h>
namespace geos {
NonRobustCGAlgorithms::NonRobustCGAlgorithms(){
}
NonRobustCGAlgorithms::~NonRobustCGAlgorithms(){
}
/**
* ring is expected to contain a closing point;
* i.e. ring[0]=ring[length-1]
*/
bool
NonRobustCGAlgorithms::isPointInRing(const Coordinate& p,
const CoordinateSequence* ring){
int i,i1; // point index;i1=i-1 mod n
double xInt; // x intersection of e with ray
int crossings=0; // number of edge/ray crossings
double x1,y1,x2,y2;
int nPts=ring->getSize();
/* For each line edge l=(i-1,i),see if it crosses ray from test point in positive x direction. */
for(i=1;i<nPts;i++){
i1=i-1;
const Coordinate& p1=ring->getAt(i);
const Coordinate& p2=ring->getAt(i1);
x1=p1.x-p.x;
y1=p1.y-p.y;
x2=p2.x-p.x;
y2=p2.y-p.y;
if(((y1>0) && (y2<=0)) || ((y2>0)&&(y1<=0))) {
/* e straddles x axis,so compute intersection. */
xInt=(x1*y2-x2*y1)/(y2-y1);
//xsave=xInt;
/* crosses ray if strictly positive intersection. */
if(0.0<xInt)
crossings++;
}
}
/* p is inside if an odd number of crossings. */
if((crossings%2)==1)
return true;
else
return false;
}
//bool NonRobustCGAlgorithms::isOnLine(const Coordinate& p,const CoordinateSequence* pt) const {
// for(int i=1;i<pt->getSize();i++){
// const Coordinate& p0=pt->getAt(i-1);
// const Coordinate& p1=pt->getAt(i);
// li->computeIntersection(p,p0,p1);
// if(li->hasIntersection())
// return true;
// }
// return false;
//}
/**
* Computes whether a ring defined by an array of {@link Coordinate} is
* oriented counter-clockwise.
* <p>
* This will handle coordinate lists which contain repeated points.
*
* @param ring an array of coordinates forming a ring
* @return <code>true</code> if the ring is oriented counter-clockwise.
* @throws IllegalArgumentException if the ring is degenerate (does not contain 3 different points)
*/
bool NonRobustCGAlgorithms::isCCW(const CoordinateSequence* ring){
// # of points without closing endpoint
int nPts=ring->getSize()-1;
// check that this is a valid ring - if not, simply return a dummy value
if (nPts<4) return false;
// algorithm to check if a Ring is stored in CCW order
// find highest point
Coordinate hip=ring->getAt(0);
int hii=0;
for(int i=1;i<=nPts;i++) {
Coordinate p=ring->getAt(i);
if(p.y>hip.y){
hip.setCoordinate(p);
hii=i;
}
}
// find different point before highest point
int iPrev=hii;
do {
iPrev = (iPrev - 1) % nPts;
} while (ring->getAt(iPrev)==hip && iPrev != hii);
// find different point after highest point
int iNext = hii;
do {
iNext = (iNext + 1) % nPts;
} while (ring->getAt(iNext)==hip && iNext != hii);
Coordinate prev=ring->getAt(iPrev);
Coordinate next=ring->getAt(iNext);
if (prev==hip || next==hip || prev==next)
throw new IllegalArgumentException("degenerate ring (does not contain 3 different points)");
// translate so that hip is at the origin.
// This will not affect the area calculation,and will avoid
// finite-accuracy errors(i.e very small vectors with very large coordinates)
// This also simplifies the discriminant calculation.
double prev2x=prev.x-hip.x;
double prev2y=prev.y-hip.y;
double next2x=next.x-hip.x;
double next2y=next.y-hip.y;
// compute cross-product of vectors hip->next and hip->prev
//(e.g. area of parallelogram they enclose)
double disc=next2x*prev2y-next2y*prev2x;
/* If disc is exactly 0,lines are collinear. There are two possible cases:
(1)the lines lie along the x axis in opposite directions
(2)the line lie on top of one another
(2)should never happen,so we're going to ignore it!
(Might want to assert this)
(1)is handled by checking if next is left of prev==>CCW
*/
if(disc==0.0){
// poly is CCW if prev x is right of next x
return(prev.x>next.x);
} else {
// if area is positive,points are ordered CCW
return(disc>0.0);
}
}
int
NonRobustCGAlgorithms::computeOrientation(const Coordinate& p1,
const Coordinate& p2,const Coordinate& q){
double dx1=p2.x-p1.x;
double dy1=p2.y-p1.y;
double dx2=q.x-p2.x;
double dy2=q.y-p2.y;
double det=dx1*dy2-dx2*dy1;
if(det>0.0)return 1;
if(det<0.0)return-1;
return 0;
}
}
|
