Geomentric Intersection:O( (n+k) * log(n))算法 Bentiey & Ottmann Algorithm
标签: 计算几何
课程:计算几何:Geometric Intersection
书籍:计算几何:算法与应用
BO算法用于求解线段组中所有的交点。在实现中,我忽略了课程中邓老师所描述的退化情况。
实现:
POINT_H
//Point.h
#pragma once
#ifndef POINT_H
#define POINT_H
#include<iostream>
using namespace std;
class Point{
private:
float x, y;
int sign; //1---线段的起点 2---线段的终点 3---两条线段的交点
public:
Point() :x(0.0f), y(0.0f) {}
Point(float xx, float yy) :x(xx), y(yy), sign(0){}
Point(float xx, float yy, int s) :x(xx), y(yy), sign(s){}
Point(const Point& p) :x(p.x),y(p.y), sign(p.sign){}
~Point() {}
float getX() { return x; }
float getY() { return y; }
void setSign(int s) { sign = s; }
int getSign() { return sign; };
friend ostream& operator<<(ostream& os,const Point& p) {
if (p.sign == 1)
os << "begin_point";
else if (p.sign == 2)
os << "end_point";
else if(p.sign == 3)
os << "intersect_point";
os << "( " << p.x << "," << p.y << " )";
return os;
}
bool operator == (const Point& p)const {
if (this->x == p.x && this->y == p.y)
return true;
return false;
}
//重定义Priority Queue排序
bool operator < (const Point& p)const {
if (this->x < p.x || (this->x == p.x && this->y < p.y))
return false;
return true;
}
};
#endif // !POINT_H
SEGMENT_H
#pragma once
#ifndef SEGMENT_H
#define SEGMENT_H
#include "Point.h"
class Segment {
private:
Point begin, end; //线段的两个端点
public:
Segment() {}
Segment(Point a, Point b) {
if (a.getX() < b.getX() || (a.getX() == b.getX() && a.getY() < b.getY())) {
begin = a;
end = b;
}
else {
begin = b;
end = a;
}
begin.setSign(1);
end.setSign(2);
}
Segment(Segment& temp) :begin(temp.begin), end(temp.end){}
~Segment() {}
Point getBeginPoint() {
return begin;
}
Point getEndPoint() {
return end;
}
bool ifIntersect(Segment s) { //判断两个线段是否相交---四次ToLeft测试,如果线段两端点分别在另一条线段的左右两侧,则相交
bool judge1, judge2;
judge1 = judge2 = false;
bool tlt1 = ToLeftTest(this->begin, this->end, s.begin);
bool tlt2 = ToLeftTest(this->begin, this->end, s.end);
bool tlt3 = ToLeftTest(s.begin, s.end, this->begin);
bool tlt4 = ToLeftTest(s.begin, s.end, this->end);
if (tlt1 != tlt2)
judge1 = true;
if (tlt3 != tlt4)
judge2 = true;
if (judge1 && judge2)
return true;
return false;
}
bool ToLeftTest(Point b, Point e, Point t) {
//vector(b,e) X vector(b,t)
float result = e.getX() * t.getY() + b.getX() * e.getY() + t.getX() * b.getY()
- (e.getX() * b.getY() + b.getX() * t.getY() + t.getX() * e.getY());
if (result > 0.0f) //left
return true;
else if (result < 0.0f) //right
return false;
}
friend ostream& operator<<(ostream& os, Segment& temp) {
os << "Segment:" << " begin_point" << temp.getBeginPoint() << "end_point" << temp.getEndPoint() << endl;
return os;
}
};
#endif // !SEGMENT_H
AVL_H
#pragma once
#ifndef AVL_H
#define AVL_H
#include<vector>
#include "Segment.h"
using namespace std;
class AVLTreeNode{
public:
Segment segment;
Point current_point;
int height;
AVLTreeNode *lchild, *rchild;
AVLTreeNode *prenode, *nextnode;
AVLTreeNode() {}
AVLTreeNode(Segment s,Point p, AVLTreeNode *l, AVLTreeNode *r) :
segment(s), current_point(p), lchild(l), rchild(r),height(0){ }
AVLTreeNode(AVLTreeNode &t) :segment(t.segment),current_point(t.current_point),height(t.height),lchild(t.lchild),
rchild(t.rchild), prenode(t.prenode), nextnode(t.nextnode){}
friend ostream& operator<< (ostream& os, AVLTreeNode& temp) {
os << "segment" << temp.segment << endl;
return os;
}
};
class AVLTree {
private:
AVLTreeNode *mRoot;
vector<AVLTreeNode*>arr; //用于查找某一线段节点前驱和后继的索引数组
public:
AVLTree() :mRoot(NULL){}
~AVLTree() {}
//获取树的高度
int getHeight() {
return height(mRoot);
}
//获取树的高度
int max(int a, int b) {
return a > b ? a : b;
}
//插入线段
AVLTreeNode* insertElement(Segment s,float current_x) {
return insert(mRoot, s, current_x);
}
//删除线段
AVLTreeNode* removeElement(Point s) {
AVLTreeNode* z;
if ( (z = search(mRoot, s)) != NULL) {
mRoot = remove(mRoot, z);
}
return z;
}
//获得AVL树中最大节点
AVLTreeNode* getMaximum() {
return maximum(mRoot);
}
//获得AVL中最小节点
AVLTreeNode* getMinimum() {
return minimum(mRoot);
}
//打印AVL
void print() {
printAVLTree(mRoot);
}
//更新所有线段的当前节点
void update(float cx) {
updateCP(mRoot, cx);
}
//重置索引数组
int resetvector() {
arr.clear();
resetvector(mRoot);
return arr.size();
}
//输出queue
void printvector() {
printvector(arr);
}
//返回该线段在索引数组中的下标,以此来得到该线段的前驱与后继
int search_pre_next(Point p, int type) {
return search_pre_next(arr, p, type);
}
//根据索引返回线段
Segment getArrSegment(int index) {
if (arr.size())
return arr[index]->segment;
}
//交换AVL树中的两个线段节点
void exchange(int n1, int n2) {
AVLTreeNode* t1 = arr[n1];
AVLTreeNode* t2 = arr[n2];
AVLTreeNode temp = *t1;
t1->segment = t2->segment;
t2->segment = temp.segment;
}
private:
//获取树的高度
int height(AVLTreeNode* tree) {
if (tree != NULL)
return tree->height;
return 0;
}
//调整树的平衡
//LL
AVLTreeNode* leftLeftRotation(AVLTreeNode* k2) {
AVLTreeNode* k1 = k2->lchild;
k2->lchild = k1->rchild;
k1->rchild = k2;
k2->height = max(height(k2->lchild), height(k2->rchild)) + 1;
k1->height = max(height(k1->rchild), k2->height) + 1;
return k1;
}
//RR
AVLTreeNode* rightRightRotation(AVLTreeNode* k2) {
AVLTreeNode* k1 = k2->rchild;
k2->rchild = k1->lchild;
k1->lchild = k2;
k2->height = max(height(k2->lchild), height(k2->rchild)) + 1;
k1->height = max(height(k1->rchild), k2->height) + 1;
return k1;
}
//RL == LL + RR
AVLTreeNode* rightLeftRotation(AVLTreeNode* k2) {
k2->rchild = leftLeftRotation(k2->rchild);
return rightRightRotation(k2);
}
//LR == RR + LL
AVLTreeNode* leftRightRotation(AVLTreeNode* k2) {
k2->lchild = rightRightRotation(k2->lchild);
return leftLeftRotation(k2);
}
//将节点插入到AVL树中
AVLTreeNode* insert(AVLTreeNode* &tree, Segment s,float cx) {
if (tree == NULL) {
tree = new AVLTreeNode(s, s.getBeginPoint(), NULL, NULL);
if (tree == NULL) {
cout << "ERROR: create avltree node failed!" << endl;
}
}
else if(s.getBeginPoint() < tree->current_point){ //插入到左子树中
tree->lchild = insert(tree->lchild, s,cx);
if (height(tree->lchild) - height(tree->rchild) == 2) {
if (s.getBeginPoint() < tree->current_point)
tree = leftLeftRotation(tree);
else
tree = leftRightRotation(tree);
}
}
else if (tree->current_point < s.getBeginPoint()) { //插入到右子树中
tree->rchild = insert(tree->rchild, s,cx);
if (height(tree->rchild) - height(tree->lchild) == 2){
if (tree->current_point < s.getBeginPoint())
tree = rightRightRotation(tree);
else
tree = rightLeftRotation(tree);
}
}
else {
cout << "the element has been existed" << endl;
}
tree->height = max(height(tree->lchild), height(tree->rchild)) + 1;
return tree;
}
//将节点从AVL树中删除,返回被删除的节点
AVLTreeNode* remove(AVLTreeNode* &tree, AVLTreeNode* z) {
if (tree == NULL || z == NULL)
return NULL;
if (z->current_point == tree->current_point) { //找到删除节点,进行删除操作
if (tree->lchild != NULL && tree->rchild != NULL) {
if (height(tree->lchild) > height(tree->rchild)) {
AVLTreeNode* max = maximum(tree->lchild);
tree->segment = max->segment;
tree->current_point = max->current_point;
tree->lchild = remove(tree->lchild, max);
}
else {
AVLTreeNode* min = minimum(tree->rchild);
tree->segment = min->segment;
tree->current_point = min->current_point;
tree->rchild = remove(tree->rchild, min);
}
}
else {
AVLTreeNode* temp = tree;
tree = (tree->lchild != NULL) ? tree->lchild : tree->rchild;
delete temp;
}
}
else if (z->current_point < tree->current_point) { //继续从左子树中查找
tree->lchild = remove(tree->lchild, z);
if (height(tree->rchild) - height(tree->lchild) == 2) {
AVLTreeNode* r = tree->rchild;
if (height(r->lchild) > height(r->rchild))
tree = rightLeftRotation(tree);
else
tree = rightRightRotation(tree);
}
}
else if (tree->current_point < z->current_point ) { //继续从右子树中查找
tree->rchild = remove(tree->rchild, z);
if (height(tree->lchild) - height(tree->rchild) == 2) {
AVLTreeNode* l = tree->lchild;
if (height(l->rchild) > height(l->lchild))
tree = leftRightRotation(tree);
else
tree = leftLeftRotation(tree);
}
}
return tree;
}
//从AVL树中查找特定权值的节点
AVLTreeNode* search(AVLTreeNode* tree, Point p) {
if (tree == NULL || tree->current_point == p)
return tree;
if (p < tree->current_point)
return search(tree->lchild, p);
else
return search(tree->rchild, p);
}
//从AVL树中找到权值最大的节点
AVLTreeNode* maximum(AVLTreeNode* tree) {
if (tree == NULL)
return NULL;
while (tree->rchild != NULL)
tree = tree->rchild;
return tree;
}
//从AVL树中找到权值最小的节点
AVLTreeNode* minimum(AVLTreeNode* tree) {
if (tree == NULL)
return NULL;
while (tree->lchild != NULL)
tree = tree->lchild;
return tree;
}
//打印AVL树
void printAVLTree(AVLTreeNode* tree) {
if (tree == NULL)
return;
printAVLTree(tree->lchild);
cout << " " << tree->segment;
printAVLTree(tree->rchild);
}
//更新所有线段的当前节点
void updateCP(AVLTreeNode* mRoot, float cx) {
if (mRoot == NULL)
return;
updateCP(mRoot->lchild, cx);
float cy, bx, by, ex, ey;
bx = mRoot->segment.getBeginPoint().getX();
by = mRoot->segment.getBeginPoint().getY();
ex = mRoot->segment.getEndPoint().getX();
ey = mRoot->segment.getEndPoint().getY();
cy = by + (ey - by) * (cx - bx) / (ex - bx);
mRoot->current_point = Point(cx, cy);
updateCP(mRoot->rchild, cx);
}
//设置索引数组
void resetvector(AVLTreeNode* tree) {
if (tree == NULL)
return;
resetvector(tree->lchild);
arr.push_back(tree);
resetvector(tree->rchild);
}
//输出索引数组
void printvector(vector<AVLTreeNode*> arr) {
for (int i = 0; i < arr.size(); i++) {
cout << arr[i]->segment;
}
cout << "******" << endl;
}
//返回该线段在索引数组中的下标,以此来得到该线段的前驱与后继
int search_pre_next(vector<AVLTreeNode*>arr, Point p, int type) {
if (type == 0) { //从线段的起始点中寻找
for (int i = 0; i < arr.size(); i++) {
if (p == arr[i]->segment.getBeginPoint()) {
return i;
}
}
}
else if(type == 1){ //从线段的终点中寻找
for (int i = 0; i < arr.size(); i++) {
if (p == arr[i]->segment.getEndPoint())
return i;
}
}
else { //从当前节点中寻找
for (int i = 0; i < arr.size(); i++) {
if (p == arr[i]->current_point)
return i;
}
}
return -1;
}
};
#endif // !AVL_H
OB ALGORITHM_CPP
/* 输入:平面线段集s
** 输出:s中各线段之间的所有交点
** 算法:OB Algorithm
*/
#include "Segment.h"
#include "Avl.h"
#include<queue>
void InitializeSegment(Segment* s, int num); //初始化线段
void PrintPriorityQueue(priority_queue<Point> eq); //输出
void CalculateIntersectionPoint(Segment* S, int segmentNum, AVLTree* tree, priority_queue<Point> eq); //计算所有交点
Segment findSegment(Segment* S, int segmentNum, Point begin); //根据起始点寻找线段
Point GetIntersectPoint(Segment s1, Segment s2); //求两线段的交点
int main() {
//初始化线段信息
int segmentNum;
cin >> segmentNum;
Segment* S = new Segment[segmentNum];
InitializeSegment(S, segmentNum);
//输出线段信息
for (int i = 0; i < segmentNum; i++) {
cout << S[i];
}
//OB算法的两个数据结构:SS(Status Structure)---BBST、 EQ(Event Queue)---Priority Queue
//EQ
priority_queue<Point> eq;
for (int i = 0; i < segmentNum; i++) {
eq.push(S[i].getBeginPoint());
eq.push(S[i].getEndPoint());
}
PrintPriorityQueue(eq);
//SS
AVLTree *tree = new AVLTree();
//开始计算交点
CalculateIntersectionPoint(S, segmentNum, tree, eq);
return 0;
}
void InitializeSegment(Segment* S, int segmentNum) {
float a, b, c, d;
for (int i = 0; i < segmentNum; i++) {
cin >> a >> b >> c >> d;
S[i] = Segment(Point(a, b), Point(c, d));
}
}
void PrintPriorityQueue(priority_queue<Point> eq) {
while (!eq.empty()) {
Point *temp = new Point(eq.top());
cout << *temp << endl;
delete temp;
temp = NULL;
eq.pop();
}
}
void CalculateIntersectionPoint(Segment* S, int segmentNum, AVLTree* tree, priority_queue<Point> eq) {
while (!eq.empty()) {
Point temp = Point(eq.top());
float current_x = temp.getX();
//更新线段树节点中,线段与当前扫描线交点
tree->update(current_x);
if (temp.getSign() == 1) { //线段起始点:增加新的线段、判断新线段与其相邻线段是否有交点
Segment cur_s = findSegment(S, segmentNum, temp);
AVLTreeNode* newnode = tree->insertElement(cur_s, current_x);
//tree->printvector();
//将新加入的线段与它的前驱和后继线段进行比较
//判断是否会有新交点出现
int nodesize = tree->resetvector();
int index = tree->search_pre_next(temp, 0);
int pre, next;
Segment pre_s, next_s;
if (index == 0)
pre = -1;
else
pre = index - 1;
if (index == nodesize - 1)
next = -1;
else
next = index + 1;
if (pre != -1) {
pre_s = tree->getArrSegment(pre);
if (cur_s.ifIntersect(pre_s)) { //存在交点
eq.push(GetIntersectPoint(cur_s, pre_s));
}
}
if (next != -1) {
next_s = tree->getArrSegment(next);
if (cur_s.ifIntersect(next_s)) { //存在交点
eq.push(GetIntersectPoint(cur_s, next_s));
}
}
}
else if (temp.getSign() == 2) { //线段终点:删除指定线段、判断是否有两个不相邻线段变成相邻线段,并求相邻线段交点
int nodesize = tree->resetvector();
int index = tree->search_pre_next(temp, 1);
AVLTreeNode* deletenode = tree->removeElement(temp);
int pre, next;
Segment pre_s, next_s;
if (index == 0)
pre = -1;
else
pre = index - 1;
if (index == nodesize - 1)
next = -1;
else
next = index + 1;
if (pre != -1 && next != -1) {
pre_s = tree->getArrSegment(pre);
next_s = tree->getArrSegment(next);
if (next_s.ifIntersect(pre_s)) { //存在交点
Point ip = GetIntersectPoint(next_s, pre_s);
if(ip.getX() >= current_x)
eq.push(ip);
}
}
}
else if (temp.getSign() == 3) { //线段交点:交换两个线段的次序,判断新的相邻关系,求新的交点
cout << "the intersection point :" << temp << endl;
/*1、根据交点找到两条交线 以及前驱、后继
**2、在AVL树中交换两个线段
**3、做两次相交判断
*/
//1、找交线以及前驱、后继
int nodesize = tree->resetvector();
int index,pre,next;
Segment intersect_s1, intersect_s2, pre_s, next_s;
index = tree->search_pre_next(temp, 2);
pre = index - 1;
next = index + 2;
intersect_s1 = tree->getArrSegment(index);
intersect_s2 = tree->getArrSegment(index + 1);
//2、交换AVL树中的两个线段节点
tree->exchange(index, index + 1);
//3、由于优先级发生变化,要做相交判断
if (pre >= 0) {
pre_s = tree->getArrSegment(pre);
if (intersect_s2.ifIntersect(pre_s)) {
Point ip = GetIntersectPoint(intersect_s2, pre_s);
if (ip.getX() >= current_x)
eq.push(ip);
}
}
if (next < nodesize) {
next_s = tree->getArrSegment(next);
if (intersect_s1.ifIntersect(next_s)) {
Point ip = GetIntersectPoint(intersect_s1, next_s);
if (ip.getX() >= current_x)
eq.push(ip);
}
}
}
eq.pop();
}
}
Segment findSegment(Segment* S, int segmentNum, Point begin) {
for (int i = 0; i < segmentNum; i++) {
if (S[i].getBeginPoint() == begin)
return S[i];
}
exit(-1);
}
Point GetIntersectPoint(Segment s1, Segment s2) {
Point a1 = s1.getBeginPoint();
Point a2 = s1.getEndPoint();
Point b1 = s2.getBeginPoint();
Point b2 = s2.getEndPoint();
float t = ((b1.getX() - a1.getX()) * (b2.getY() - b1.getY()) - (b1.getY() - a1.getY()) * (b2.getX() - b1.getX()))
/ ((a2.getX() - a1.getX()) * (b2.getY() - b1.getY()) - (a2.getY() - a1.getY()) * (b2.getX() - b1.getX()));
float newx = a1.getX() + t * (a2.getX() - a1.getX());
float newy = a1.getY() + t * (a2.getY() - a1.getY());
return Point(newx, newy, 3);
}
用例:
6
6 18 49 39
80 29 67 34
90 76 10 4
0 1 -5 -9
0 -3 -6 -9
0 -4 -4 0
测试: