線分の交差判定

[みけCAT]

Windows Vista Home Premium SP2 32ビット
Dev-C++5.2.0.3
gcc4.7.0

線分の交差判定のプログラムを作っています。
参考までに、方程式を使った古いプログラムを載せます。

► スポイラーを表示

コード:

int cross(double x1,double y1,double x2,double y2,
	double x3,double y3,double x4,double y4) {
	double kaA,seA,kaB,seB,x,y,uex,uey,sitax,sitay;
	if(x1==x3 && y1==y3)return 1;
	if(x1==x4 && y1==y4)return 1;
	if(x2==x3 && y2==y3)return 1;
	if(x2==x4 && y2==y4)return 1;
	if(x1!=x2 && x3!=x4) {
		kaA=(y2-y1)/(x2-x1);
		seA=y1-x1*kaA;
		kaB=(y4-y3)/(x4-x3);
		seB=y3-x3*kaB;
		if(kaA==kaB) {
			if(seA==seB) {
				if(x1>x2){uex=x1;sitax=x2;} else {uex=x2;sitax=x1;}
				if((x3>=sitax && x3<=uex) || (x4>=sitax && x4<=uex))return 1;
				return 0;
			} else return 0;
		}
		x=(seB-seA)/(kaA-kaB);
		y=kaA*x+seA;
		if(x1>x2){uex=x1;sitax=x2;} else {uex=x2;sitax=x1;}
		if(y1>y2){uey=y1;sitay=y2;} else {uey=y2;sitay=y1;}
		if(y1!=y2) {
			if(x<sitax || x>uex || y<sitay || y>uey)return 0;
		} else {
			if(x<sitax || x>uex)return 0;
		}
		if(x3>x4){uex=x3;sitax=x4;} else {uex=x4;sitax=x3;}
		if(y3>y4){uey=y3;sitay=y4;} else {uey=y4;sitay=y3;}
		if(y3!=y4) {
			if(x<sitax || x>uex || y<sitay || y>uey)return 0;
		} else {
			if(x<sitax || x>uex)return 0;
		}
	}
	if(x1==x2 && x3!=x4) {
		kaB=(y4-y3)/(x4-x3);
		seB=y3-x3*kaB;
		y=kaB*x1+seB;
		if(y1>y2){uey=y1;sitay=y2;} else {uey=y2;sitay=y1;}
		if(x3>x4){uex=x3;sitax=x4;} else {uex=x4;sitax=x3;}
		if(y<sitay || y>uey || x1<sitax || x1>uex)return 0;
	}
	if(x1!=x2 && x3==x4) {
		kaA=(y2-y1)/(x2-x1);
		seA=y1-x1*kaA;
		y=kaA*x3+seA;
		if(y3>y4){uey=y3;sitay=y4;} else {uey=y4;sitay=y3;}
		if(x1>x2){uex=x1;sitax=x2;} else {uex=x2;sitax=x1;}
		if(y<sitay || y>uey || x3<sitax || x3>uex)return 0;
	}
	if(x1==x2 && x3==x4) {
		if(x1!=x3)return 0;
		if(y1>y2){uey=y1;sitay=y2;} else {uey=y2;sitay=y1;}
		if(y3>y4){uex=y3;sitax=y4;} else {uex=y4;sitax=y3;}
		if(sitax>uey)return 0;
		if(uex>sitay)return 0;
	}
	return 1;
}

さて、問題はこの古いコードではなく、この新しいコードです。

コード:

int cross(int x1,int y1,int x2,int y2,
		int x3,int y3,int x4,int y4) {
	/* ベクトル */ 
	long long b12x,b12y,b34x,b34y;
	long long b24x,b24y;
	long long b13x,b13y,b14x,b14y,b32x,b32y;
	/* ベクトルの外積のz成分 */ 
	long long g1213,g1214,g3431,g3432;
	/* ベクトルの内積 */ 
	long long n1213,n1214,n3431,n3432;
	long long n2123,n2124,n4341,n4342;
	/* 座標を引き算してベクトルを計算する */ 
	b12x=x2-x1;b12y=y2-y1;
	b34x=x4-x3;b34y=y4-y3;
	b13x=x3-x1;b13y=y3-y1;
	b14x=x4-x1;b14y=y4-y1;
	b32x=x2-x3;b32y=y2-y3;
	b24x=x4-x2;b24y=y4-y2;
	/* 平面ベクトルをz成分が0の空間ベクトルとみなし、
	   外積のz成分を計算する */ 
	g1213=b12x*b13y-b13x*b12y;
	g1214=b12x*b14y-b14x*b12y;
	g3431=b34x*(-b13y)-(-b13x)*b34y;
	g3432=b34x*b32y-b32x*b34y;
	/* ベクトルの内積を計算する */ 
	n1213=b12x*b13x+b12y*b13y;
	n1214=b12x*b14x+b12y*b14y;
	n3431=b34x*(-b13x)+b34y*(-b13y);
	n3432=b34x*b32x+b34y*b32y;
	n2123=(-b12x)*(-b32x)+(-b12y)*(-b32y);
	n2124=(-b12x)*b24x+(-b12y)*b24y;
	n4341=(-b34x)*(-b14x)+(-b34y)*(-b14y);
	n4342=(-b34x)*(-b24x)+(-b34y)*(b24y);

	/* ベクトルの外積の符号が違っていれば、交わっている */
	if(
		((g1213<0 && g1214>0) || (g1213>0 && g1214<0)) &&
		((g3431<0 && g3432>0) || (g3431>0 && g3432<0))
	) return 1;
	/* 線分の端が他方の線分上にあるので、交わっている */
	if(
		(g1213==0 && n1213>=0 && n2123>=0) ||
		(g1214==0 && n1214>=0 && n2124>=0) ||
		(g3431==0 && n3431>=0 && n4341>=0) ||
		(g3432==0 && n3432>=0 && n4342>=0)
	) return 1;
	/* 交わっていない */ 
	return 0;
}

※引数の型が古い関数ではdoubleで、新しい関数ではintですが、仕様です

「ベクトルの外積の符号が違っていれば、交わっている」の部分では、
点3と点4が直線12（点1と点2を結ぶ直線）に関して反対側にあり、
かつ点1と点2が直線34に関して反対側にある場合に、「交わっている」と判定しています。

この考え方と実装で本当に正しく動作するでしょうか？

DXライブラリを用いて視覚化したプログラムを添付します。
動かしたい点の数字のキー（テンキーはだめ）を押しながらカーソルキーで点を動かせます。
このとき、シフトキーを押していると1ドットづつ動かせます。

よろしくお願いします。

[みけCAT]

なかなか返信が来ないので、マルチポストさせていただきました。
http://sasuke.main.jp/bbs/index.cgi?mod ... &rev=&no=0

オフトピック

禁止でしたっけ？

[みけCAT]

POJのJack Strawsでこの判定コードを使用したところ、Wrong Answerになりました。
しかし、バグを修正してもう一度Submitすると、きちんとAcceptedになりました。
したがって、これで正しいものとみなします。

バグを修正したコード

コード:

int cross(int x1,int y1,int x2,int y2,
		int x3,int y3,int x4,int y4) {
	/* ベクトル */ 
	long long b12x,b12y,b34x,b34y;
	long long b24x,b24y;
	long long b13x,b13y,b14x,b14y,b32x,b32y;
	/* ベクトルの外積のz成分 */ 
	long long g1213,g1214,g3431,g3432;
	/* ベクトルの内積 */ 
	long long n1213,n1214,n3431,n3432;
	long long n2123,n2124,n4341,n4342;
	/* 座標を引き算してベクトルを計算する */ 
	b12x=x2-x1;b12y=y2-y1;
	b34x=x4-x3;b34y=y4-y3;
	b13x=x3-x1;b13y=y3-y1;
	b14x=x4-x1;b14y=y4-y1;
	b32x=x2-x3;b32y=y2-y3;
	b24x=x4-x2;b24y=y4-y2;
	/* 平面ベクトルをz成分が0の空間ベクトルとみなし、
	   外積のz成分を計算する */ 
	g1213=b12x*b13y-b13x*b12y;
	g1214=b12x*b14y-b14x*b12y;
	g3431=b34x*(-b13y)-(-b13x)*b34y;
	g3432=b34x*b32y-b32x*b34y;
	/* ベクトルの内積を計算する */ 
	n1213=b12x*b13x+b12y*b13y;
	n1214=b12x*b14x+b12y*b14y;
	n3431=b34x*(-b13x)+b34y*(-b13y);
	n3432=b34x*b32x+b34y*b32y;
	n2123=(-b12x)*(-b32x)+(-b12y)*(-b32y);
	n2124=(-b12x)*b24x+(-b12y)*b24y;
	n4341=(-b34x)*(-b14x)+(-b34y)*(-b14y);
	n4342=(-b34x)*(-b24x)+(-b34y)*(-b24y);

	/* ベクトルの外積の符号が違っていれば、交わっている */
	if(
		((g1213<0 && g1214>0) || (g1213>0 && g1214<0)) &&
		((g3431<0 && g3432>0) || (g3431>0 && g3432<0))
	) return 1;
	/* 線分の端が他方の線分上にあるので、交わっている */
	if(
		(g1213==0 && n1213>=0 && n2123>=0) ||
		(g1214==0 && n1214>=0 && n2124>=0) ||
		(g3431==0 && n3431>=0 && n4341>=0) ||
		(g3432==0 && n3432>=0 && n4342>=0)
	) return 1;
	/* 交わっていない */ 
	return 0;
}

検証用コード

コード:

#include <stdio.h>

#define UNION_FIND_MAX 12

int uf_table[UNION_FIND_MAX];

void uf_init(void) {
	int i;
	for(i=0;i<UNION_FIND_MAX;i++)uf_table[i]=0;
}

void uf_merge(int a,int b) {
	int next;
	while(uf_table[a]>0) {
		next=uf_table[a]-1;
		if(uf_table[next]>0)uf_table[a]=uf_table[next];
		a=next;
	}
	while(uf_table[b]>0) {
		next=uf_table[b]-1;
		if(uf_table[next]>0)uf_table[b]=uf_table[next];
		b=next;
	}
	if(a!=b)uf_table[a]=b+1;
}

int uf_issame(int a,int b) {
	int next;
	while(uf_table[a]>0) {
		next=uf_table[a]-1;
		if(uf_table[next]>0)uf_table[a]=uf_table[next];
		a=next;
	}
	while(uf_table[b]>0) {
		next=uf_table[b]-1;
		if(uf_table[next]>0)uf_table[b]=uf_table[next];
		b=next;
	}
	return a==b;
}

int cross(int x1,int y1,int x2,int y2,
		int x3,int y3,int x4,int y4) {
	long long b12x,b12y,b34x,b34y;
	long long b24x,b24y;
	long long b13x,b13y,b14x,b14y,b32x,b32y;
	long long g1213,g1214,g3431,g3432;
	long long n1213,n1214,n3431,n3432;
	long long n2123,n2124,n4341,n4342;

	b12x=x2-x1;b12y=y2-y1;
	b34x=x4-x3;b34y=y4-y3;
	b13x=x3-x1;b13y=y3-y1;
	b14x=x4-x1;b14y=y4-y1;
	b32x=x2-x3;b32y=y2-y3;
	b24x=x4-x2;b24y=y4-y2;

	g1213=b12x*b13y-b13x*b12y;
	g1214=b12x*b14y-b14x*b12y;
	g3431=b34x*(-b13y)-(-b13x)*b34y;
	g3432=b34x*b32y-b32x*b34y;

	n1213=b12x*b13x+b12y*b13y;
	n1214=b12x*b14x+b12y*b14y;
	n3431=b34x*(-b13x)+b34y*(-b13y);
	n3432=b34x*b32x+b34y*b32y;
	n2123=(-b12x)*(-b32x)+(-b12y)*(-b32y);
	n2124=(-b12x)*b24x+(-b12y)*b24y;
	n4341=(-b34x)*(-b14x)+(-b34y)*(-b14y);
	n4342=(-b34x)*(-b24x)+(-b34y)*(-b24y);

	if(
		((g1213<0 && g1214>0) || (g1213>0 && g1214<0)) &&
		((g3431<0 && g3432>0) || (g3431>0 && g3432<0))
	) return 1;
	if(
		(g1213==0 && n1213>=0 && n2123>=0) ||
		(g1214==0 && n1214>=0 && n2124>=0) ||
		(g3431==0 && n3431>=0 && n4341>=0) ||
		(g3432==0 && n3432>=0 && n4342>=0)
	) return 1;

	return 0;
}

int n;
int x1[12],y1[12],x2[12],y2[12];

int main(void) {
	int i,j;
	while(1) {
		scanf("%d",&n);
		if(n==0)break;
		for(i=0;i<n;i++)scanf("%d%d%d%d",&x1[i],&y1[i],&x2[i],&y2[i]);
		uf_init();
		for(i=0;i<n;i++) {
			for(j=0;j<n;j++) {
				if(cross(x1[i],y1[i],x2[i],y2[i],
					x1[j],y1[j],x2[j],y2[j]))uf_merge(i,j);
			}
		}
		while(1) {
			scanf("%d%d",&i,&j);
			if(i==0 && j==0)break;
			if(!uf_issame(i-1,j-1))printf("NOT ");
			puts("CONNECTED");
		}
	}
	return 0;
}