GIS点在多边形内算法

(vl

：

np)

内的程 序：

以下是引用片段：

Copy code

int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */

const vertex_t* v)

{

int i, j, k1, k2, c;

rect_t rc;

vertex_t w;

if (np < 3)

return 0;

vertices_get_extent(vl, np, &rc);

if (v->x < _x || v->x > _x || v->y < _y || v->y > _y)

return 0;

/* Set a horizontal beam l(*v, w) from v to the ultra right */

w.x = _x + DBL_EPSILON;

w.y = v->y;

c = 0; /* Intersection points counter */

for(i=0; i

{

j = (i+1) % np;

if(is_intersect(vl+i, vl+j, v, &w))

{

c++;

}

else if(vl[i].y==w.y)

{

k1 = (np+i-1)%np;

while(k1!=i && vl[k1].y==w.y)

k1 = (np+k1-1)%np;

k2 = (i+1)%np;

while(k2!=i && vl[k2].y==w.y)

k2 = (k2+1)%np;

if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0)

c++;

if(k2 <= i)

break;

i = k2;

}

}

return c%2;

}

本想配些插图说明问题，但是，

CS DN

的文章里放图片我还没用过。以后再试吧

!

实践

证明，

本程序算法的适应性极强。

但是，

对于点正好落在多边形边上的极端情形，

有可能 得

出

2

种不同的结果。

3.

判断一个几何点在多边形内的例子：

#define

X

0

#define

Y

1

typedef

enum

{

FALSE,

TRUE

}

bool;

#define

DIM

2

/*

Dimension

of

points

*/

typedef

int

tPointi[DIM];

/*

type

integer

point

*/

typedef

double

tPointd[DIM];

/*

type

double

point

*/

#define

PMAX

10000

/*

Max

#

of

pts

in

polygon

*/

typedef

tPointi

tPolygoni[PMAX];/*

type

integer

polygon

*/

char

InPoly(

tPointi

q,

tPolygoni

P

,

int

n

);

void

PrintPoly(

int

n,

tPolygoni

P

);

void

PrintPoint(

tPointi

p

);

void

Copy

(

tPolygoni

a,

tPolygoni

b

,

int

n

);

main()

{

int

n;

tPolygoni

P

,

Porig;

tPointi

q;

n

=

ReadPoly(

P

);

Copy(

P

,

Porig,

n

);

while(

scanf(

%d

&q[X],

&q[Y])

!=

EOF

)

{

printf(

=

%cn

InPoly(

q,

P

,

n

)

);

/*

Refill

the

destroyed

polygon

with

original.

*/

Copy(

Porig,

P

,

n

);

}

}

/*

InPoly

returns

a

char

in

{i,o,v,e}.

See

above

for

definitions.

*/

char

InPoly(

tPointi

q,

tPolygoni

P

,

int

n

)

{

int

i,

i1;

/*

point

index;

i1

=

i-1

mod

n

*/

int

d;

/*

dimension

index

*/

double

x;

/*

x

intersection

of

e

with

ray

*/

int

Rcross

=

0;

/*

number

of

right

edge/ray

crossings

*/

int

Lcross

=

0;

/*

number

of

left

edge/ray

crossings

*/

printf(

q

=

PrintPoint(q);

putchar('n');

/*

Shift

so

that

q

is

the

origin.

Note

this

destroys

the

polygon.

This

is

done

for

pedogical

clarity.

*/

for(

i

=

0;

i

<

n;

i++

)

{

for(

d

=

0;

d

<

DIM;

d++

)

P[i][d]

=

P[i][d]

-

q[d];

}

/*

For

each

edge

e=(i-1,i),

see

if

crosses

ray.

*/

for(

i

=

0;

i

<

n;

i++

)

{

/*

First

see

if

q=(0,0)

is

a

vertex.

*/

if

(

P[i][X]==0

&&

P[i][Y]==0

)

return

'v';

i1

=

(

i

+

n

-

1

)

%

n;

/*

printf(

i1,

i);

*/

/*

if

e

the

x-axis...

*/

/*

The

commented-out

statement

is

logically

equivalent

to

the

one

following.

*/

/*

if(

(

(

P[i][Y]

>

0

)

&&

(

P[i1][Y]

<=

0

)

)

||

(

(

P[i1][Y]

>

0

)

&&

(

P[i]

[Y]

<=

0

)

)

)

{

*/

if(

(

P[i][Y]

>

0

)

!=

(

P[i1][Y]

>

0

)

)

{

/*

e

straddles

ray

,

so

compute

intersection

with

ray

.

*/

x

=

(P[i][X]

*

(double)P[i1][Y]

-

P[i1][X]

*

(double)P[i][Y])

/

(double)(P[i1][Y]

-

P[i][Y]);

/*

printf(

x

=

%gt

x);

*/

/*

crosses

ray

if

strictly

positive

intersection.

*/

if

(x

>

0)

Rcross++;

}

/*

printf(

cross=%dt

Rcross);

*/

/*

if

e

straddles

the

x-axis

when

reversed...

*/

/*

if(

(

(

P[i]

[Y]

<

0

)

&&

(

P[i1][Y]

>=

0

)

)

||

(

(

P[i1][Y]

<

0

)

&&

(

P[i]

[Y]

>=

0

)

)

)

{

*/

if

(

(

P[i][Y]

<

0

)

!=

(

P[i1][Y]

<

0

)

)

{

/*

e

straddles

ray

,

so

compute

intersection

with

ray

.

*/

x

=

(P[i][X]

*

P[i1][Y]

-

P[i1][X]

*

P[i][Y])

/

(double)(P[i1][Y]

-

P[i][Y]);

/*

printf(

x

=

%gt

x);

*/

/*

crosses

ray

if

strictly

positive

intersection.

*/

if

(x

<

0)

Lcross++;

}

/*

printf(

cross=%dn

Lcross);

*/

}

/*

q

on

the

edge

if

left

and

right

cross

are

not

the

same

parity.

*/

if(

(

Rcross

%

2

)

!=

(Lcross

%

2

)

)

return

'e';

/*

q

inside

iff

an

odd

number

of

crossings.

*/

if(

(Rcross

%

2)

==

1

)

return

'i';

else

return

'o';

}

void

PrintPoint(

tPointi

p

)

{

int

i;

putchar('(');

for

(

i

=

0;

i

<

DIM;

i++

)

{

printf(

p[i]);

if

(

i

!=

DIM-1

)

putchar(',');

}

putchar(')');

}

/*

Reads

in

the

coordinates

of

the

vertices

of

a

polygon

from

puts

them

into

P

,

and

returns

n,

the

number

of

vertices.

Formatting

conventions:

etc.

*/

int

ReadPoly(

tPolygoni

P

)

{

int

i,

n;

do

{

printf(

the

number

of

vertices:n

scanf(

&n

);

if

(

n

<=

PMAX

)

break;

printf(

in

read_poly:

too

many

points;

max

is

%dn

}

while

(

1

);

printf(

);

printf(

i

x

yn

for

(

i

=

0;

i

<

n;

i++

)

{

stdin,

PMAX);

scanf(

%d

&P[i][0],

&P[i][1]

);

printf(

i,

P[i][0],

P[i][1]);

}

printf(

=

%3d

vertices

readn

putchar('n');

return

n;

}

void

PrintPoly(

int

n,

tPolygoni

P

)

{

int

i;

printf(

printf(

i

x

yn

for(

i

=

0;

i

<

n;

i++

)

printf(

i,

P[i][0],

P[i][1]);

}

/*

Copy

polygon

a

to

b

(overwriting

b).

*/

void

Copy(

tPolygoni

a,

tPolygoni

b,

int

n

)

{

int

i,

j;

for

(

i=0;

i

<

n;

i++)

for

(

j

=

0;

j

<

DIM;

j++

)

b[i][j]

=

a[i][j];

}

bool

EqPoint(

tPointi

a,

tPointi

b

)

{

int

i;

for

(

i

=

0;

i

<

DIM;

i++

)

if

(

a[i]

!=

b[i])

return

FALSE;

return

TRUE;

}

1.

已知点

point(x,y)

和多边形

Polyg on

（

x1,y1;x2,y2;….xn,yn;

2.

< br>以

point

为起点，以无穷远为终点作平行于

X

轴的直线

line(

x,y; -

∞

,y)

；

3.

循环取得

(f or(i=0;i

多边形的每一条边

side

(xi,yi;xi+1,yi+1)

，

且判断是否平行于

X

轴，

如果平行< /p>

continue

，

否则，

i++

；

4.

同时判断

point(x,y)

是否在

side

上，如果是，则返回

1(

点在多边形

上

)

，否则继续下面的判断；

5.

判断线

side

与

line

是否有交点，如果有则

count++