523066680

这么多点，撸主搞图形的？

523066680

回复 7# tommytangtang


   
这里面88是什么意思？不是Z轴吧，数据类型不一样
100.477 85.911 88
96.747 84.002 88

523066680

回复 14# tommytangtang


　　按你的平行线的方法应该就可以的，他给的是书，可以慢慢看。
这个还是解决当下问题要紧吧？ 以前我还折腾一点儿，现在也记不起来咯。
楼主给个公式我可以用perl写写。（1楼那个公式太潦草……）
图我做出来了，还好不是凹多边形，（gl老版本，见笑）

523066680

其实这个问题应该看这本书：《实时碰撞检测算法技术》而且我手上有一本，
当时纯粹手贱买下的，从来没看过……  http://www.tup.com.cn/book/Showbook.asp?CPBH=033422-01&DJ=52


第五章就有  点与多边形之间的测试、点与三角形之间的测试（多边形可以分为多个三角形）

523066680

回复 19# CrLf


    你把筛选的顶点分两批晒一下，我可以脚本绘图 （纯粹找事做）

523066680

回复 21# CrLf

虽说Perl和python都有很多绘图接口，但是数据表达上的话还是python比较好用
不过我也忘了python怎么用了， =_= 正在用C语言+OPENGL绘图

523066680

回复 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

523066680

回复 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

第一张图是昨晚的。第二张图才是今天的

523066680

回复 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

#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 │
────────┘

523066680

回复 29# DAIC


    少将，你为何有两颗太阳！

523066680

回复 31# tommytangtang

　　这个是纯C，而且算法也不是碰撞法什么的，碰撞测试只是书名，专门介绍各种2D、3D图形、物理碰撞算法的书籍
这本书网上有电子档的：http://download.csdn.net/download/huangtao198805/5488335
该算法的具体解释我明天再弄了。

523066680

回复 40# tommytangtang

现在凸多边形和凹多边形的解决方案都有了， 你会VB的话就把函数照抄，转化为VB就可以了。
（话说就单单从文本读取数据，VB比gawk慢是真的吗？）

今天略忙，忘了分析算法了。太晚了，闪！

523066680

回复 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为原点，计算出另外两个点

的相对坐标，视为向量ａ、ｂ，计算这两个向量的向量积（cross product））

（高等数学第一章有介绍这个计算和性质，或者在网上搜索：叉积……）

计算向量积得到的是垂直于ａ和ｂ所在平面的向量V，遵循右手法则：右手握住V

向量（垂直于ａ和ｂ），手指从ａ向ｂ环绕，大拇指的方向和向量方向一致。

也就是如果ａ-> ｂ->原点p 是逆时针方向则V>0 ，是反方向则V<0，如果ａｂ平

行或者其中一方为0，则V=0

　　　　　　　　　　　　　　　　[attach]7613[/attach]

523066680

判断是否为逆时针的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                                        │
──────────────────────────┘