trajectory.cpp

一个可执行的nbc算法设计与实现 有助大家互相探讨学习

#include "stdafx.h"
#include "trajectory.h"


double caculate_perpen(point *si,point *ei,point *sj,point *ej)
{
	
	double x1=sj->x-si->x;
	double y1=sj->y-si->y;
	
	double x2=ei->x-si->x;
	double y2=ei->y-si->y;
	
	double x3=ej->x-si->x;
	double y3=ej->y-si->y;
	
	double x=pow(x2,2)+pow(y2,2);
	
	double u1=(x1*x2+y1*y2)/x;
	double u2=(x3*x2+y3*y2)/x;
    
    double psx=si->x+u1*x2;
	double psy=si->y+u1*y2;
	
	double pex=si->x+u2*x2;
	double pey=si->y+u2*y2;
	
	double Lper1=sqrt(pow(psx-sj->x,2)+pow(psy-sj->y,2));
	double Lper2=sqrt(pow(pex-ej->x,2)+pow(pey-ej->y,2));
	double d_perpen;
	
	if(Lper1==0&&Lper2==0)
	{
		d_perpen=0;
	}
	else
	{
		d_perpen=(pow(Lper1,2)+pow(Lper2,2))/(Lper1+Lper2);
	}
	
	
	return d_perpen;
	
	
}


double caculate_angle(point *si,point *ei,point *sj,point *ej)
{
	double x1=ej->x-sj->x;
	double y1=ej->y-sj->y;
	
	double x2=ei->x-si->x;
	double y2=ei->y-si->y;
	
	double xi=sqrt(pow(x1,2)+pow(y1,2));
	double yi=sqrt(pow(x2,2)+pow(y2,2));
	
	double cos=(x1*x2+y1*y2)/(xi*yi);
	float m=cos*cos;
	float s=1-m;
	float t=sqrt(s);

	double angle=sqrt(pow(x1,2)+pow(y1,2))*t;
	
	return angle;
	
}

double distance(point *p1,point *p2)
{
	double x=p1->x-p2->x;
	double y=p1->y-p2->y;
	return sqrt(pow(x,2)+pow(y,2));
}


void Trajectory::insert_point(double px,double py,int num)
{//向轨道链表中插入读入的点
	point *p=new point;	
	p->x=px;
	p->y=py;
    p->num=num;
	p->next=NULL;

	if(first==NULL)
	{
	  first=last=p;		  
	}
	else
	{
		last->next=p;	
		last=p;	
	}
}


void Trajectory::add_character(point *p)
{
  point *np=new point;
  np->x=p->x;
  np->y=p->y;
  np->num=p->num;
  np->next=NULL;
  last_character_point->next=np;  
  last_character_point=np;
}


double Trajectory::caculate_mdl_par(point *pstart,point *pcurrent)
{//
  	double LH;
	double LDH;
	
	LH=log10(sqrt(pow(pstart->x-pcurrent->x,2)+pow(pstart->y-pcurrent->y,2)))/log10(2);
	
	double temp_perpen=0;
	double temp_angle=0;

	point *p=pstart;
	do {
		if(distance(pstart,pcurrent)>=distance(p,p->next))
		{
			temp_perpen+=caculate_perpen(pstart,pcurrent,p,p->next);
		    temp_angle+=caculate_angle(pstart,pcurrent,p,p->next);
		}
		else
		{
			temp_perpen+=caculate_perpen(p,p->next,pstart,pcurrent);
			temp_angle+=caculate_angle(p,p->next,pstart,pcurrent);
		}		
		p=p->next;
	} while(p!=pcurrent);

	if(temp_perpen==0||temp_angle==0)
	{
		LDH=0;
	}
	else
	{
		LDH=(log10(temp_perpen)+log10(temp_angle))/log10(2);
	}
	
	return (LH+LDH);

}


double Trajectory::caculate_mdl_nopar(point *pstart,point *pcurrent)
{//
 	point *p=pstart;
	double tempdis=0;

	do {
		tempdis+=distance(p,p->next);
		p=p->next;
	} while(p!=pcurrent);

	return log10(tempdis)/log10(2);
}


void Trajectory::Trajectory_patition(double s)
{//轨道划分算法核心函数
   point *cp=new point;
   cp->x=first->x;
   cp->y=first->y;//轨道的第一个点就是特征点
   cp->num=first->num;
   first_character_point=last_character_point=cp;
   cp->next=NULL;
   
   point *start_p;
   start_p=first;

   point *current_p;//当前处理的点
   point *prev_p;   //前一个点
   
   prev_p=start_p;
   current_p=start_p->next;

   while(current_p->next!=NULL)
   {
      
	  double cost_par=caculate_mdl_par(start_p,current_p);
	  double cost_nopar=caculate_mdl_nopar(start_p,current_p);
	  if(cost_par>(cost_nopar+s))
	  {
		  add_character(prev_p);
          add_character(prev_p);
		  start_p=prev_p;		  
	  }
	  else
	  {
		  prev_p=current_p;
		  current_p=current_p->next;
	  }
   }
   add_character(current_p);//最后一个点
   
}


line *Trajectory::get_first_line()
{
	return first_ls;
}

line *Trajectory::get_last_line()
{
	return last_ls;
}


void Trajectory::del_point()
{
	point *p1;
	p1=first;
	while (p1!=NULL)
	{
	 point *p2;
	 p2=p1;
	 p1=p1->next;
	 delete p2;
	}
	first=last=NULL;
}

void Trajectory::del_chara_point()
{
	first_character_point=last_character_point=NULL;
}

void Trajectory::get_linelist()
{
	line *lp=new line;
	point *tp=first_character_point;
	

	lp->clusterId=0;
	lp->line_mark=0;
	lp->tra_id=temp_tra_id;
	lp->kNB=0;
	lp->R_kNB=0;
	lp->p1=tp;
	lp->p2=tp->next;
	lp->neigh_next=NULL;
	lp->next=NULL;
	if(first_ls==NULL)
	{
		first_ls=last_ls=lp;
	}
	else
	{
	  last_ls->next=lp;
      last_ls=lp;
	}


	tp=tp->next->next;	
	
	while(tp!=NULL)
	{          
	      line *ls=new line;
		  ls->clusterId=0;
		  ls->line_mark=0;
		  ls->tra_id=temp_tra_id;
		  ls->kNB=0;
		  ls->R_kNB=0;
		  ls->p1=tp;
		  ls->p2=tp->next;
		  ls->neigh_next=NULL;
		  ls->next=NULL;
		  last_ls->next=ls;
		  last_ls=ls;
		 
		  tp=tp->next->next;
	}	
	temp_tra_id++;
}

void Trajectory::ATP_Algorithm(char *filename,double s)
{
	ifstream infile;
	infile.open(filename,ios::in);
	if(!infile)
	{
		cout<<"打开文件失败！"<<endl;
		exit(0);
	}
	int num;
	int tra_id;
	int tempid=0;
	double x;
	double y;
	
	while(!infile.eof())
	{
		infile>>num;
		infile>>tra_id;
		infile>>x;
		infile>>y;
		if(tra_id==tempid)
		{
			insert_point(x,y,num);			
		}
		else
		{
			Trajectory_patition(s);
            get_linelist();
			del_point();
            del_chara_point();

			tempid=tra_id;
			insert_point(x,y,num);
		}
	}

	Trajectory_patition(s);
	get_linelist();
	del_point();
	del_chara_point();
}


void Trajectory::out_ls()
{
	line *ls=first_ls;
	int num=0;
	while(ls!=NULL)
	{
		num++;
		ls=ls->next;
	}
	cout<<num<<endl;
}


void Trajectory::output_ls()
{
    ofstream outfile;
	outfile.open("ls.txt",ios::out);
	if(!outfile)
	{
		exit(0);
	}
	
	line *ls=first_ls;
	int num=0;
	
	while(ls!=NULL)
	{
		outfile<<num<<"     "<<ls->tra_id<<"       "<<ls->p1->x<<"      "<<ls->p1->y<<"      "<<ls->p2->x<<"      "<<ls->p2->y<<endl;
        num++;
		ls=ls->next;
	}
	outfile.close();
}