标题: [文本处理] 批处理用怎么实现快速判断点与多边形关系呢 [打印本页]
作者: tommytangtang 时间: 2014-8-17 19:38 标题: 批处理用怎么实现快速判断点与多边形关系呢
多边形顶点数据如下:
A
96.802,87.23
89.094,78.623
92.268,64.808
115.391,61.863
124.913,77.265
115.845,90.174
102.696,92.439
96.802,87.23
待判断的点
如果点在多边形外就保存下来>>c.txt
如果点在多边形内或者多边形边上就剔除它>>d.txt
基本思路是用射线法,计算以待判断点(yt,xt)为起点,
平行向右或者向上发射射线,依次与每条多边形的边发射,判断有没有交点,有偶数个交点包括0,就说明在多边形内,那么这一行数据>>d.txt。奇数个就>>c.txt
当某条边的两个顶点x坐标 x(1),x(2)满足如下条件时
If xt >= x1 And xt < x2 or xt <= x1 And xt > x2 then ::就是说要依次找到相交的可能性,再去判断下面的
set yt=yi+(y2-y1)(xt-x1)/(x2-x1)
if yt > y0 Then Ncross = Ncross + 1 ::如果有交点,交点数Ncross就+1
每个点遍历每条边,if Ncross%2=0,>>d.txt
else if Ncross%2=1,>>c.txt
因为需要判断的点数据量非常大(有几百万个),用vb还可以做出来,但是速度太慢了,最近感受到了gawk的速度,想用gawk或者sed去解决,但是又涉及到二维数组,头大,哪位大神能帮帮忙
作者: CrLf 时间: 2014-8-17 19:47
sed 不用考虑,gawk 可能可以
如果数量大,还是用 perl 或者其他什么更合适的语言吧
另,Ncross+1 改用 Ncross= not Ncross 更好,这样就不需要计算 Ncross%2 了
作者: tommytangtang 时间: 2014-8-17 20:06
回复 2# CrLf
批处理应该要比用vb编出来的速度快多了吧,我想先用批处理试试,批处理的数组像vb一样吗,论坛里一搜出来的都是vbs数组,vbs又是什么
作者: CrLf 时间: 2014-8-17 20:10
唔...其实批处理做这个是极其慢的
vbs 效率远比 vb 低,还是用 c 写个好了,最快
作者: 523066680 时间: 2014-8-17 20:23
本帖最后由 523066680 于 2014-8-17 20:29 编辑
这么多点,撸主搞图形的?
作者: tommytangtang 时间: 2014-8-17 20:25
回复 4# CrLf
额。。。前几天处理固定格式文本数据的时候你让我扔了个gawk到win32里
让我感受到了它处理文本速度上的快感,我还以为这种处理也很速度呢,好好的研究了几天,全网搜索了好多教程,但是一写起来,它的判断啊,数组让我崩溃了。。。
原来他不擅长这种啊。。。
作者: tommytangtang 时间: 2014-8-17 20:28
回复 5# 523066680
算是吧,经常用autocad的,工作需要自己弄点小程序
作者: neorobin 时间: 2014-8-17 20:29
回复 6# tommytangtang
推荐书籍
计算机图形学的算法基础
计算机图形学几何工具算法详解
维基百科: http://en.wikipedia.org/wiki/Point_in_polygon
作者: 523066680 时间: 2014-8-17 20:29
回复 7# tommytangtang
这里面88是什么意思?不是Z轴吧,数据类型不一样
100.477 85.911 88
96.747 84.002 88
作者: tommytangtang 时间: 2014-8-17 20:33
回复 9# 523066680
哈哈,被你发现了,最后一列是我乱加上去的,本来应该是不同的数据表示z轴的,88吉利
作者: tommytangtang 时间: 2014-8-17 20:43
回复 8# neorobin
谢谢推荐,我说的思路和维基百科的差不多(全英文的只看懂了大概。。。)
作者: CrLf 时间: 2014-8-17 20:46
回复 6# tommytangtang
纯批很慢,第三方工具很快
作者: neorobin 时间: 2014-8-17 21:20
找到些代码
http://rosettacode.org/wiki/Ray-casting_algorithm
作者: tommytangtang 时间: 2014-8-17 21:35
回复 13# neorobin
大神给的都太高端了,全E文经常看不懂,他处理的多边形只是矩形吗,那么多语言写出来了,好像就freebasic好理解一点
作者: 523066680 时间: 2014-8-17 21:37
回复 14# tommytangtang
按你的平行线的方法应该就可以的,他给的是书,可以慢慢看。
这个还是解决当下问题要紧吧? 以前我还折腾一点儿,现在也记不起来咯。
楼主给个公式我可以用perl写写。(1楼那个公式太潦草……)
图我做出来了,还好不是凹多边形,(gl老版本,见笑)
作者: tommytangtang 时间: 2014-8-17 21:52
回复 15# 523066680
谢谢热心的版版,我这个也不急的,这个是我的举例,如果数据量只有这么几个数的话,我用autocad也可以解决的,主要是大批量的判断,就非常难搞了
作者: neorobin 时间: 2014-8-17 22:10
回复 14# tommytangtang
处理的是一般多边形, 在C代码中有测试用的多边形, 顶点坐标可以看出是一个 正4边形, 和一个 凹 8边形- vec vsq[] = { {0,0}, {10,0}, {10,10}, {0,10},
- {2.5,2.5}, {7.5,0.1}, {7.5,7.5}, {2.5,7.5}};
-
- polygon_t sq = { 4, vsq }, /* outer square */
- sq_hole = { 8, vsq }; /* outer and inner square, ie hole */
复制代码
作者: 523066680 时间: 2014-8-17 22:45
其实这个问题应该看这本书:《实时碰撞检测算法技术》而且我手上有一本,
当时纯粹手贱买下的,从来没看过…… http://www.tup.com.cn/book/Showbook.asp?CPBH=033422-01&DJ=52
第五章就有 点与多边形之间的测试、点与三角形之间的测试(多边形可以分为多个三角形)
作者: CrLf 时间: 2014-8-18 00:37
按顶楼算法实现了基本功能,也不懂得到的结果是不是正确的,撸主搞个样本来比对下吧:
http://bbs.bathome.net/viewthread.php?tid=31486
作者: 523066680 时间: 2014-8-18 00:51
回复 19# CrLf
你把筛选的顶点分两批晒一下,我可以脚本绘图 (纯粹找事做)
作者: CrLf 时间: 2014-8-18 00:55
回复 20# 523066680
那敢情好,有劳您老人家了!
顺便请教一下用什么脚本绘图比较方便
作者: 523066680 时间: 2014-8-18 01:12
回复 21# CrLf
虽说Perl和python都有很多绘图接口,但是数据表达上的话还是python比较好用
不过我也忘了python怎么用了, =_= 正在用C语言+OPENGL绘图
作者: 523066680 时间: 2014-8-18 01:28
回复 21# CrLf
C.txt 有一个点在里面, 是哪个还不清楚[attach]7585[/attach]
B.txt 的坐标全在外面
我自己抄了个算法,结果也不对,难道是绘图有问题?晕。
2014-08-18, 按碰撞检测算法书上的算法做了一个
C/C++的代码以及链接参考附件
Perl代码:- use IO::Handle;
- STDOUT->autoflush(1);
-
- my $coordfile="a.txt";
- my @poly=(
- 96.802,87.23,
- 89.094,78.623,
- 92.268,64.808,
- 115.391,61.863,
- 124.913,77.265,
- 115.845,90.174,
- 102.696,92.439,
- #96.802,87.23,
- );
-
- our @V;
- my $i=0;
- while (scalar(@poly)>0) {
- $V[$i]{x}=shift @poly;
- $V[$i]{y}=shift @poly;
- $i++;
- }
-
- my ($x, $y);
- my %dot;
- open READ,"<",$coordfile or die "$!";
- open INPOLY, ">", "inside.txt" or die;
- open OUTPOLY, ">", "outside.txt" or die ;
- foreach (<READ>) {
- ($x, $y, undef) = split(/ |,/, $_);
- #print $x," , ",$y,"\n";
- $dot{x} = $x;
- $dot{y} = $y;
-
- if (&pointInPolygon(\%dot) eq "true") {
- print INPOLY "$x, $y\r\n";
- print "($x, $y) in polygon\n";
- } else {
- print OUTPOLY "$x, $y\r\n";
- }
- }
- close READ;
- close INPOLY;
- close OUTPOLY;
-
- sub pointInPolygon {
- our @V;
- my $p = shift;
- my $npoints = $#V+1;
- my ($low, $high) = (0, $npoints);
- do {
- $mid = ($low + $high) / 2;
- $mid=int($mid)+1 if ($mid > int($mid));
- if ( &TriangleIsCCW($V[0], $V[$mid], $p) eq "true" ) {
- $low = $mid;
- } else {
- $high = $mid;
- }
- } while ( ($low + 1) < $high);
-
- if( ($low == 0) or ($high == $npoints) ) {
- #print "$npoints $low $mid $high\n";
- return "false";
- }
- return &TriangleIsCCW($V[$low], $V[$high], $p);
- }
-
- sub TriangleIsCCW {
- my ($p0, $p1, $p2) = @_;
- my $cross;
- my ($d0, $d1, $d2, $d3);
- $d0 = ($p1->{x} - $p0->{x});
- $d1 = ($p1->{y} - $p0->{y});
- $d2 = ($p2->{x} - $p0->{x});
- $d3 = ($p2->{y} - $p0->{y});
-
- $cross = &Cross($d0, $d1, $d2, $d3);
- if ($cross > 0) {
- return "true";
- } else {
- return "false";
- }
- }
-
- sub Cross {
- #T Cross(T x0, T y0, T x1, T y1) {
- #//|x0 y0| = |a b|
- #//|x1 y1| |c d|
- #//ad * bc
- #//x0 * y1 - x1 * y0
- my ($x0, $y0, $x1, $y1) = @_;
- my $result;
- $result = ($x0 * $y1) - ($y0 * $x1);
- return $result;
- }
复制代码
在区域范围内的点:
100.477, 85.911
96.747, 84.002
96.817, 75.030
105.692, 67.907
97.253, 70.233
111.950, 80.117
105.931, 76.600
111.138, 85.624
114.157, 74.060
作者: tommytangtang 时间: 2014-8-18 15:55
回复 23# 523066680
结果不对啊,挑出来的坐标和之前我给的a.txt都对不上了,之前给的x没有小于10的,基本都是50以上
版主绘图没问题,我打印到图上也和你一样的
图中,红色的点是a.txt。蓝色的点是得到的结果。
作者: 523066680 时间: 2014-8-18 16:01
本帖最后由 523066680 于 2014-8-18 16:09 编辑
回复 24# tommytangtang
这里面分明没有小于10的
100.477, 85.911
96.747, 84.002
96.817, 75.030
105.692, 67.907
97.253, 70.233
111.950, 80.117
105.931, 76.600
111.138, 85.624
114.157, 74.060
第一张图是昨晚的。第二张图才是今天的
作者: CrLf 时间: 2014-8-18 16:04
回复 23# 523066680
昨天刚好看到书上说可以把:- my $result;
- $result = ($x0 * $y1) - ($y0 * $x1);
- return $result;
复制代码
简化成:- ($x0 * $y1) - ($y0 * $x1);
复制代码
作者: 523066680 时间: 2014-8-18 16:06
本帖最后由 523066680 于 2014-8-18 18:41 编辑
回复 26# CrLf
谢谢提醒,昨天手工将C转Perl后BUG一堆,所以很谨慎地把句子分开了 =_=
就是这里有个取反…… 然后return cross > 0 真精简……- template<class T>
- bool triangleIsCCW(
- T x0, T y0,
- T x1, T y1,
- T x2, T y2) {
- auto cross = -(Cross<T>(x1-x0,y1-y0,x2-x0,y2-y0));
- return cross > 0;
- }
复制代码
作者: 523066680 时间: 2014-8-18 21:02 标题: 测试点与凸多边形的关系,纯C版本
本帖最后由 523066680 于 2014-8-18 21:33 编辑
- #include <stdio.h>
- #define bool int
- #define false 0
- #define true 1
-
- typedef struct {
- float x;
- float y;
- } Point;
-
- Point poly[]={
- 96.802,87.23,
- 89.094,78.623,
- 92.268,64.808,
- 115.391,61.863,
- 124.913,77.265,
- 115.845,90.174,
- 102.696,92.439,
- };
-
- float Cross(
- float x0, float y0,
- float x1, float y1
- ) {
- return (x0 * y1) - (x1 * y0);
- }
-
- int triangleIsCCW( //判断三个点的顺序是否为逆时针
- Point P0, Point P1, Point P2
- ) {
- auto cross = (
- Cross(
- P1.x-P0.x, P1.y-P0.y,
- P2.x-P0.x, P2.y-P0.y
- )
- );
- return cross >= 0;
- }
-
- bool pointInPolygon(
- Point vip, Point poly[], int n_poly
- ) {
- // points represented by [(x0,y0), ... (xn,yn)]
- int low = 0, high = n_poly;
- int mid;
- do {
- mid = (low + high) / 2;
- if (
- triangleIsCCW(poly[0], poly[mid], vip)
- ) {
- low = mid;
- } else {
- high = mid;
- }
- } while (low + 1 < high);
- if (low == 0 || high == n_poly) return false;
- return triangleIsCCW(poly[low], poly[high], vip);
- }
-
- int main(int argc, char *argv[])
- {
- FILE *fp = fopen("a.txt","r");
- FILE *wrt = fopen("inside.txt","w");
- Point vip;
- char s[1024];
- int n_poly = (sizeof(poly) / sizeof(poly[0].x)) / 2;
- int in_polygon;
- while ( ! feof(fp) ) {
- fgets(s, 1024, fp);
- if ( sscanf(s, "%f %f", &vip.x, &vip.y) < 2 ) {
- sscanf(s, "%f,%f", &vip.x, &vip.y);
- }
- in_polygon = pointInPolygon(vip, poly, n_poly);
- if ( in_polygon ) {
- printf("%.3f, %.3f\n", vip.x, vip.y);
- fprintf(wrt, "%.3f, %.3f\n", vip.x, vip.y);
- }
- }
- fclose(fp);
- fclose(wrt);
- return 0;
- }
复制代码
算法来自《实时碰撞检测算法技术》5.4.1 点与多边形之间的测试 P139
这段C的处理结果和Perl又那么一两个点的误差,主要是有些点几乎就在边缘上,
零点零几的区别 和 0 判断产生不同的结果。
────────┐
100.477, 85.911 │
96.747, 84.002 │
96.817, 75.030 │
105.692, 67.907 │
97.253, 70.233 │
115.392, 61.864 │
111.950, 80.117 │
105.931, 76.600 │
111.138, 85.624 │
114.157, 74.060 │
120.702, 70.397 │
92.268, 64.808 │
115.391, 61.863 │
────────┘
作者: DAIC 时间: 2014-8-18 21:14
回复 3# tommytangtang
vbs 就是 visual basic script
作者: 523066680 时间: 2014-8-18 21:19
回复 29# DAIC
少将,你为何有两颗太阳!
作者: tommytangtang 时间: 2014-8-18 21:26
本帖最后由 tommytangtang 于 2014-8-18 21:33 编辑
回复 28# 523066680
我之前的算法还有个问题,当有个点在多边形左边,坐标x刚好和某个多边形顶点一样,当它往右发射时,交了两次,其中刚好有个点过多边形一个顶点,算法里就统计成三个点了,判断点在内,实际上在外…
碰撞法真心没看懂…这个是c还是c#还是c++?
作者: 523066680 时间: 2014-8-18 21:37
本帖最后由 523066680 于 2014-8-18 21:44 编辑
回复 31# tommytangtang
这个是纯C,而且算法也不是碰撞法什么的,碰撞测试只是书名,专门介绍各种2D、3D图形、物理碰撞算法的书籍
这本书网上有电子档的:http://download.csdn.net/download/huangtao198805/5488335
该算法的具体解释我明天再弄了。
作者: tommytangtang 时间: 2014-8-18 21:47
回复 29# DAIC
哈,原来和vb有关,刚好有点vb基础,看了下vbs教程,好亲切,感觉是bat和vb的结合体呢
作者: tommytangtang 时间: 2014-8-18 21:50
回复 32# 523066680
原本是感觉到了gawk在处理文本的速度上的优势才发这个贴的,c完全没有自己编过,汗…算了,多学点吧。
作者: CrLf 时间: 2014-8-18 22:05
回复 31# tommytangtang
这可以通过判断交叉点 x 是否在 [x1,x2) 内来区分,但仔细想了一下,射线法还有两个问题比较郁闷:- 1、出现与射线平行的边线时,在线上的点应视为不在区域内才能避免被多次计算
- 2、有两个角点重合时,会计算两次而导致误判
复制代码
我原来的代码确实有问题,想太简单了。为了不误导人,那帖还是删了为好
作者: tommytangtang 时间: 2014-8-18 22:24
回复 35# CrLf
仔细想了想,这两种情况可能要先用算法计算后区分,要单独去判断,感觉越来越复杂了。。。
作者: neorobin 时间: 2014-8-18 22:39
回复 31# tommytangtang
我测试了一个C代码
作者: CrLf 时间: 2014-8-18 22:53
卧槽,你们这群技术宅...Holy high
作者: CrLf 时间: 2014-8-19 23:11
回复 30# 523066680
突然发现我也是两个太阳...话说这里为什么没有列出版主用户:http://bbs.bathome.net/memcp.php?action=usergroups
你看看在线时间下面那一排太阳吧
作者: tommytangtang 时间: 2014-8-19 23:57
C不懂啊,这个算是解决了吗?
作者: 523066680 时间: 2014-8-20 00:52
本帖最后由 523066680 于 2014-8-20 01:03 编辑
回复 40# tommytangtang
现在凸多边形和凹多边形的解决方案都有了, 你会VB的话就把函数照抄,转化为VB就可以了。
(话说就单单从文本读取数据,VB比gawk慢是真的吗?)
今天略忙,忘了分析算法了。太晚了,闪!
作者: tommytangtang 时间: 2014-8-20 21:30
回复 41# 523066680
真的没看懂,要翻译好难,求助版版
作者: 523066680 时间: 2014-8-20 22:06 标题: 点与凸多边形之间的测试
回复 40# tommytangtang
一下内容是原文照抄(有小小修改;我终于也沦为搬运工了啊……)
点与凸多边形之间的测试
假设n边凸多边形顶点V0, V1, …,Vn-1 呈逆时针排列。通过测试顶点P与方向直线V0、Vk(k=n/2)
的位置关系,可将多边形的测试范围缩至一半。该过程重复执行并相应地调整K值,直至发现顶点P位于多
边形的外部或位于有向直线V0Vk ~ V0Vk+1之间。对于后者,若顶点P位于有向直线VkVk+1的左侧(呈
逆时针排列),则其位于多边形内部。
以一个8条边的凸多边形作为例。首先,顶点P针对直线V0V4进行测试,并将多边形顶点范围降至一半。
由于顶点P位于A的左侧,所以丢弃多边形的右半部分,并继续对左半部分执行相同处。下一步,在顶点P与
直线B之间执行测试,最后对直线C进行测试,此时,顶点P位于邻接顶点构成的两条直线之间。最后,顶点
P与该邻接顶点构成的直线进行测试。
测试结果为:顶点P位于多边形内部。
[attach]7612[/attach]
图为:凸多边形顶点的二分搜索使得点P的包含测试的时间复杂度为O(log n)(这里测试了4条边:A~D)
另外,预计算函数TriangleIsCCW()将计算三角形顶点是否为逆时针方向排列。该测试的高效实现
方法如下列代码所示:
──────────────────────────────────────────┐
// Test if point p lies inside ccw-specified convex n-gon given by vertices v[] │
int PointInConvexPolygon(Point p, int n, Point v[]) │
{ │
// Do binary search over polygon vertices to find the fan triangle │
// (v[0], v[low], v[high]) the point p lies within the near sides of │
int low = 0, high = n; │
do { │
int mid = (low + high) / 2; │
if (TriangleIsCCW(v[0], v[mid], p)) │
low = mid; │
else │
high = mid; │
} while (low + 1 < high); │
│
// If point outside last (or first) edge, then it is not inside the n-gon │
if (low == 0 || high == n) return 0; │
│
// p is inside the polygon if it is left of │
// the directed edge from v[low] to v[high] │
return TriangleIsCCW(v[low], v[high], p); │
} │
──────────────────────────────────────────┘
该点包含测试方法类似[Preparata85]中所描述的相关算法结构,其中n边多边形内部也存在一个与
此处的顶点V0对应。
---------以下内容可能造成不适,请[跳过]或者[对上面内容进行消化]后再访问--------
该算法通过多边形的一个原点和其他顶点,对多边形进行划分(二分),通过TriangleIsCCW()函数
逐渐找到最靠近点P的两个边界,在示例中最后找到为:V0->V5 和 V0->V6 。这个过程手工演算一
次就知道原理了。
最后一步 return TriangleIsCCW(v[low], v[high], p); 即为判断:如果V5->V6->P 呈逆时针顺序,
则点P在多边形内部。
TriangleIsCCW() 函数判断传入的三个顶点数据是否呈逆时针排列,书中没有给出,但是在网上找到:
─────────────────────────────┐
template<class T> │
T Cross(T x0, T y0, T x1, T y1) { │
//|x0 y0| = |a b| │
//|x1 y1| |c d| │
│
//ad * bc │
//x0 * y1 - x1 * y0 │
│
return (x0 * y1) - (y0 * x1); │
} │
│
template<class T> │
bool triangleIsCCW( │
T x0, T y0, │
T x1, T y1, │
T x2, T y2) { │
auto cross = (Cross<T>(x1-x0,y1-y0,x2-x0,y2-y0)); │
return cross > 0; │
} │
─────────────────────────────┘
原理就是对传入的三个顶点坐标进行向量运算(以P为原点,计算出另外两个点
的相对坐标,视为向量a、b,计算这两个向量的向量积(cross product))
(高等数学第一章有介绍这个计算和性质,或者在网上搜索:叉积……)
计算向量积得到的是垂直于a和b所在平面的向量V,遵循右手法则:右手握住V
向量(垂直于a和b),手指从a向b环绕,大拇指的方向和向量方向一致。
也就是如果a-> b->原点p 是逆时针方向则V>0 ,是反方向则V<0,如果ab平
行或者其中一方为0,则V=0
[attach]7613[/attach]
作者: 523066680 时间: 2014-8-20 23:02
判断是否为逆时针的VBS实例贴出来啦,我也算图文都弄上来了,其他的楼主自己折腾。
──────────────────────────┐
│
msgbox triangleIsCCW(1.0, 0.0, 10.0, 3.0, 0.0, 0.0) │
│
msgbox triangleIsCCW(10.0, 3.0, 1.0, 0.0, 0.0, 0.0) │
│
function Cross(x0,y0, x1,y1) │
Cross = (x0 * y1) - (y0 * x1) │
end function │
│
function triangleIsCCW(x0, y0, x1, y1, x2, y2) │
result = Cross(x1-x0, y1-y0, x2-x0, y2-y0) │
if (result > 0) then │
triangleIsCCW=True │
else │
triangleIsCCW=False │
end if │
end function │
──────────────────────────┘
作者: 523066680 时间: 2014-8-20 23:09
本帖最后由 523066680 于 2014-8-20 23:10 编辑
Crlf确定加的是技术分吗……
这几段码出来我就可以作为原创内容发到各个文档网站赚点积分方便以后用吧,大概这么想的。
作者: tommytangtang 时间: 2014-8-22 00:32
回复 44# 523066680
谢谢啦,我慢慢研究一下
作者: 523066680 时间: 2014-9-5 10:33
回复 37# neorobin
发现一个文章量不算多,但是图文并茂的学术博客:
http://alienryderflex.com/polygon/
http://alienryderflex.com/tutorials.shtml
Geometry/Math
欢迎光临 批处理之家 (http://www.bathome.net/) |
Powered by Discuz! 7.2 |