diff --git a/IO/OutputHandler.cpp b/IO/OutputHandler.cpp index f8dc91ade8c76bfc6f97d0316996c4315c845eff..2a0641cfebb5e88567b74e95cf7161813ba47f61 100644 --- a/IO/OutputHandler.cpp +++ b/IO/OutputHandler.cpp @@ -89,53 +89,53 @@ void OutputHandler::ProgressBar(double TotalPeds, double NowPeds) } void OutputHandler::Write(const char* message,...) - { - char msg[CLENGTH]=""; - va_list ap; - va_start(ap, message); - vsprintf(msg, message, ap); - va_end(ap); - - string str(msg); - - if (str.find("ERROR") != string::npos) - { - cerr << msg << endl; - cerr.flush(); - incrementErrors(); - } - else if (str.find("WARNING") != string::npos) - { - cerr << msg << endl; - cerr.flush(); - incrementWarnings(); - } - else - { // infos - cout << msg << endl; - cout.flush(); - } +{ + char msg[CLENGTH]=""; + va_list ap; + va_start(ap, message); + vsprintf(msg, message, ap); + va_end(ap); + + string str(msg); + + if (str.find("ERROR") != string::npos) + { + cerr << msg << endl; + cerr.flush(); + incrementErrors(); + } + else if (str.find("WARNING") != string::npos) + { + cerr << msg << endl; + cerr.flush(); + incrementWarnings(); + } + else + { // infos + cout << msg << endl; + cout.flush(); + } } void STDIOHandler::Write(const string& str) { - if (str.find("ERROR") != string::npos) - { - cerr << str << endl; - cerr.flush(); - incrementErrors(); - } - else if (str.find("WARNING") != string::npos) - { - cerr << str << endl; - cerr.flush(); - incrementWarnings(); - } - else - { // infos - cout << str << endl; - cout.flush(); - } + if (str.find("ERROR") != string::npos) + { + cerr << str << endl; + cerr.flush(); + incrementErrors(); + } + else if (str.find("WARNING") != string::npos) + { + cerr << str << endl; + cerr.flush(); + incrementWarnings(); + } + else + { // infos + cout << str << endl; + cout.flush(); + } } FileHandler::FileHandler(const char *fn) @@ -156,19 +156,19 @@ FileHandler::~FileHandler() void FileHandler::Write(const string& str) { - if (this != NULL) { - _pfp << str << endl; - _pfp.flush(); - } - - if (str.find("ERROR") != string::npos) - { - incrementErrors(); - } - else if (str.find("WARNING") != string::npos) - { - incrementWarnings(); - } + if (this != NULL) { + _pfp << str << endl; + _pfp.flush(); + } + + if (str.find("ERROR") != string::npos) + { + incrementErrors(); + } + else if (str.find("WARNING") != string::npos) + { + incrementWarnings(); + } } void FileHandler::Write(const char* str_msg,...) @@ -184,11 +184,11 @@ void FileHandler::Write(const char* str_msg,...) string str(msg); if (str.find("ERROR") != string::npos) { - incrementErrors(); + incrementErrors(); } else if (str.find("WARNING") != string::npos) { - incrementWarnings(); + incrementWarnings(); } } diff --git a/IO/OutputHandler.h b/IO/OutputHandler.h index acc3bb512ecffe7abb88340ff7861e702431e751..1f9d98c9981ac8f2a64a428170aac70d578c8974 100644 --- a/IO/OutputHandler.h +++ b/IO/OutputHandler.h @@ -24,7 +24,7 @@ * * **/ - + #ifndef OUTPUT_HANDLER_H_ #define OUTPUT_HANDLER_H_ diff --git a/general/ArgumentParser.cpp b/general/ArgumentParser.cpp index b3eec6805f0687f053a9b9ac6d24560723e2e3d7..4259a965b9d338b19cc874cd44798800c051f979 100644 --- a/general/ArgumentParser.cpp +++ b/general/ArgumentParser.cpp @@ -129,7 +129,7 @@ ArgumentParser::ArgumentParser() - bool ArgumentParser::ParseArgs(int argc, char **argv) +bool ArgumentParser::ParseArgs(int argc, char **argv) { //special case of the default configuration ini.xml if (argc == 1) @@ -241,8 +241,8 @@ bool ArgumentParser::ParseIniFile(const string& inifile) } else { - Log->Write("Error: \tthe given trajectory format is not supported. Supply '.xml' or '.txt' format!"); - return false; + Log->Write("Error: \tthe given trajectory format is not supported. Supply '.xml' or '.txt' format!"); + return false; } string unit = xmltoa(xMainNode->FirstChildElement("trajectories")->Attribute("unit"), "m"); @@ -283,9 +283,9 @@ bool ArgumentParser::ParseIniFile(const string& inifile) // ignore the project root in this case if ( (boost::algorithm::contains(_trajectoriesLocation,":")==false) && //windows (boost::algorithm::starts_with(_trajectoriesLocation,"/") ==false)) //linux - // &&() osx + // &&() osx { - _trajectoriesLocation=_projectRootDir+_trajectoriesLocation; + _trajectoriesLocation=_projectRootDir+_trajectoriesLocation; } // in the case no file was specified, collect all xml files in the specified directory @@ -302,11 +302,11 @@ bool ArgumentParser::ParseIniFile(const string& inifile) if (boost::algorithm::ends_with(filename, fmt)) //if (filename.find(fmt)!=std::string::npos) - { - //_trajectoriesFiles.push_back(_projectRootDir+filename); - _trajectoriesFiles.push_back(filename); - Log->Write("INFO: \tInput trajectory file is\t<"+ (filename)+">"); - } + { + //_trajectoriesFiles.push_back(_projectRootDir+filename); + _trajectoriesFiles.push_back(filename); + Log->Write("INFO: \tInput trajectory file is\t<"+ (filename)+">"); + } } closedir (dir); } @@ -417,10 +417,10 @@ bool ArgumentParser::ParseIniFile(const string& inifile) _timeIntervalA = xmltoi(xMethod_A->FirstChildElement("timeInterval")->GetText()); Log->Write("INFO: \ttime interval used in Method A is <%d> frame",_timeIntervalA); for(TiXmlElement* xMeasurementArea=xMainNode->FirstChildElement("method_A")->FirstChildElement("measurementArea"); - xMeasurementArea; xMeasurementArea = xMeasurementArea->NextSiblingElement("measurementArea")) + xMeasurementArea; xMeasurementArea = xMeasurementArea->NextSiblingElement("measurementArea")) { - _areaIDforMethodA.push_back(xmltoi(xMeasurementArea->Attribute("id"))); - Log->Write("INFO: \tMeasurement area id <%d> will be used for analysis", xmltoi(xMeasurementArea->Attribute("id"))); + _areaIDforMethodA.push_back(xmltoi(xMeasurementArea->Attribute("id"))); + Log->Write("INFO: \tMeasurement area id <%d> will be used for analysis", xmltoi(xMeasurementArea->Attribute("id"))); } } } diff --git a/general/ArgumentParser.h b/general/ArgumentParser.h index 54395c3ee0cd135963b81560f329089276e60042..176732850e39c6ba6418e7a17c6f77785f2411d4 100644 --- a/general/ArgumentParser.h +++ b/general/ArgumentParser.h @@ -48,90 +48,90 @@ extern OutputHandler* Log; class ArgumentParser { private: - std::string _geometryFileName; - std::string _errorLogFile; - std::string _trajectoriesLocation; - std::string _trajectoriesFilename; - std::string _projectRootDir; - FileFormat _fileFormat; - std::vector<std::string> _trajectoriesFiles; + std::string _geometryFileName; + std::string _errorLogFile; + std::string _trajectoriesLocation; + std::string _trajectoriesFilename; + std::string _projectRootDir; + FileFormat _fileFormat; + std::vector<std::string> _trajectoriesFiles; - char _vComponent; - bool _isMethodA; - bool _isMethodB; - bool _isMethodC; - bool _isMethodD; - bool _isCutByCircle; - double _cutRadius; - int _circleEdges; - bool _isOutputGraph; - bool _isIndividualFD; - bool _isGetProfile; - double _steadyStart; - double _steadyEnd; - int _delatTVInst; - int _timeIntervalA; - std::vector<int> _areaIDforMethodA; - std::vector<int> _areaIDforMethodB; - std::vector<int> _areaIDforMethodC; - std::vector<int> _areaIDforMethodD; - float _scaleX; - float _scaleY; - int _log; + char _vComponent; + bool _isMethodA; + bool _isMethodB; + bool _isMethodC; + bool _isMethodD; + bool _isCutByCircle; + double _cutRadius; + int _circleEdges; + bool _isOutputGraph; + bool _isIndividualFD; + bool _isGetProfile; + double _steadyStart; + double _steadyEnd; + int _delatTVInst; + int _timeIntervalA; + std::vector<int> _areaIDforMethodA; + std::vector<int> _areaIDforMethodB; + std::vector<int> _areaIDforMethodC; + std::vector<int> _areaIDforMethodD; + float _scaleX; + float _scaleY; + int _log; - std::map <int, MeasurementArea*> _measurementAreas; - void Usage(const std::string file); + std::map <int, MeasurementArea*> _measurementAreas; + void Usage(const std::string file); public: - // Konstruktor - ArgumentParser(); + // Konstruktor + ArgumentParser(); - const std::string& GetTrajectoriesFilename() const; - const std::vector<std::string>& GetTrajectoriesFiles() const; - const std::string& GetTrajectoriesLocation() const; - const FileFormat& GetFileFormat() const; - const std::string& GetGeometryFilename() const; - const std::string& GetErrorLogFile() const; - const std::string& GetProjectRootDir() const; + const std::string& GetTrajectoriesFilename() const; + const std::vector<std::string>& GetTrajectoriesFiles() const; + const std::string& GetTrajectoriesLocation() const; + const FileFormat& GetFileFormat() const; + const std::string& GetGeometryFilename() const; + const std::string& GetErrorLogFile() const; + const std::string& GetProjectRootDir() const; - double GetLengthMeasurementArea() const; - polygon_2d GetMeasureArea() const; - double GetLineStartX() const; - double GetLineStartY() const; - double GetLineEndX() const; - double GetLineEndY() const; + double GetLengthMeasurementArea() const; + polygon_2d GetMeasureArea() const; + double GetLineStartX() const; + double GetLineStartY() const; + double GetLineEndX() const; + double GetLineEndY() const; - char GetVComponent() const; - int GetDelatT_Vins() const; - int GetTimeIntervalA() const; - bool GetIsMethodA() const; - bool GetIsMethodB() const; - bool GetIsMethodC() const; - bool GetIsMethodD() const; - std::vector<int> GetAreaIDforMethodA() const; - std::vector<int> GetAreaIDforMethodB() const; - std::vector<int> GetAreaIDforMethodC() const; - std::vector<int> GetAreaIDforMethodD() const; - bool GetIsCutByCircle() const; - double GetCutRadius() const; - int GetCircleEdges() const; - bool GetIsOutputGraph() const; - bool GetIsIndividualFD() const; - double GetSteadyStart() const; - double GetSteadyEnd() const; - bool GetIsGetProfile() const; - float GetScaleX() const; - float GetScaleY() const; - int GetLog() const; - bool ParseArgs(int argc, char **argv); + char GetVComponent() const; + int GetDelatT_Vins() const; + int GetTimeIntervalA() const; + bool GetIsMethodA() const; + bool GetIsMethodB() const; + bool GetIsMethodC() const; + bool GetIsMethodD() const; + std::vector<int> GetAreaIDforMethodA() const; + std::vector<int> GetAreaIDforMethodB() const; + std::vector<int> GetAreaIDforMethodC() const; + std::vector<int> GetAreaIDforMethodD() const; + bool GetIsCutByCircle() const; + double GetCutRadius() const; + int GetCircleEdges() const; + bool GetIsOutputGraph() const; + bool GetIsIndividualFD() const; + double GetSteadyStart() const; + double GetSteadyEnd() const; + bool GetIsGetProfile() const; + float GetScaleX() const; + float GetScaleY() const; + int GetLog() const; + bool ParseArgs(int argc, char **argv); - MeasurementArea* GetMeasurementArea(int id); + MeasurementArea* GetMeasurementArea(int id); - /** - * parse the initialization file - * @param inifile - */ - bool ParseIniFile(const std::string& inifile); + /** + * parse the initialization file + * @param inifile + */ + bool ParseIniFile(const std::string& inifile); }; #endif /*ARGPARSER_H_*/ diff --git a/general/Macros.h b/general/Macros.h index 044163b68cebe6a27846dd029a46c8c57adaea5d..d178e718894ed51f1d959e05a2b7cdb608f319e6 100644 --- a/general/Macros.h +++ b/general/Macros.h @@ -25,7 +25,7 @@ * * **/ - + #ifndef _MACROS_H #define _MACROS_H @@ -117,12 +117,12 @@ enum RoutingStrategy { }; enum OperativModels { - MODEL_GFCM=1, - MODEL_GOMPERTZ, -// MODEL_ORCA, -// MODEL_CFM, -// MODEL_VELO -// MODEL_GNM + MODEL_GFCM=1, + MODEL_GOMPERTZ, + // MODEL_ORCA, + // MODEL_CFM, + // MODEL_VELO + // MODEL_GNM }; enum AgentColorMode { @@ -247,26 +247,26 @@ inline std::string concatenate(std::string const& name, int i) { inline void _printDebugLine(const std::string& fileName, int lineNumber) { -unsigned found = fileName.find_last_of("/\\"); -std::cerr << "["<< lineNumber << "]: ---"<< fileName.substr(found+1)<< " ---"<<std::endl; + unsigned found = fileName.find_last_of("/\\"); + std::cerr << "["<< lineNumber << "]: ---"<< fileName.substr(found+1)<< " ---"<<std::endl; } #define dtrace(...) \ - (_printDebugLine(__FILE__, __LINE__), \ - fprintf(stderr, __VA_ARGS__), \ - (void) fprintf(stderr, "\n")) + (_printDebugLine(__FILE__, __LINE__), \ + fprintf(stderr, __VA_ARGS__), \ + (void) fprintf(stderr, "\n")) #define derror(...) \ - (_printDebugLine(__FILE__, __LINE__), \ - fprintf(stderr, "ERROR: "), \ - fprintf(stderr, __VA_ARGS__) \ - ) + (_printDebugLine(__FILE__, __LINE__), \ + fprintf(stderr, "ERROR: "), \ + fprintf(stderr, __VA_ARGS__) \ + ) #else #define dtrace(...) ((void) 0) #define derror(...) \ - (fprintf(stderr, __VA_ARGS__) \ - ) + (fprintf(stderr, __VA_ARGS__) \ + ) #endif /* TRACE_LOGGING */ #endif /* _MACROS_H */ diff --git a/methods/Method_B.cpp b/methods/Method_B.cpp index cc97b7b150bb46d65b1d467446ca00272c8c2a59..13a8c83025946fa9b9f11048fa771bc26d93ef7c 100644 --- a/methods/Method_B.cpp +++ b/methods/Method_B.cpp @@ -34,14 +34,14 @@ using std::vector; Method_B::Method_B() { - _xCor = NULL; - _yCor = NULL; - _tIn = NULL; - _tOut = NULL; - _DensityPerFrame = NULL; - _fps = 10; - _NumPeds =0; - _areaForMethod_B = NULL; + _xCor = NULL; + _yCor = NULL; + _tIn = NULL; + _tOut = NULL; + _DensityPerFrame = NULL; + _fps = 10; + _NumPeds =0; + _areaForMethod_B = NULL; } Method_B::~Method_B() @@ -51,100 +51,100 @@ Method_B::~Method_B() bool Method_B::Process (const PedData& peddata) { - _trajName = peddata.GetTrajName(); - _projectRootDir = peddata.GetProjectRootDir(); - _fps =peddata.GetFps(); - _peds_t = peddata.GetPedsFrame(); - _NumPeds = peddata.GetNumPeds(); - _xCor = peddata.GetXCor(); - _yCor = peddata.GetYCor(); - _measureAreaId = boost::lexical_cast<string>(_areaForMethod_B->_id); - _tIn = new int[_NumPeds]; // Record the time of each pedestrian entering measurement area - _tOut = new int[_NumPeds]; - for (int i=0; i<_NumPeds; i++) - { - _tIn[i] = 0; - _tOut[i] = 0; - } - GetTinTout(peddata.GetNumFrames()); - Log->Write("------------------------Analyzing with Method B-----------------------------"); - if(_areaForMethod_B->_length<0) - { - Log->Write("Error:\tThe measurement area length for method B is not assigned!"); - exit(0); - } - else - { - GetFundamentalTinTout(_DensityPerFrame,_areaForMethod_B->_length); - } - delete []_tIn; - delete []_tOut; - return true; + _trajName = peddata.GetTrajName(); + _projectRootDir = peddata.GetProjectRootDir(); + _fps =peddata.GetFps(); + _peds_t = peddata.GetPedsFrame(); + _NumPeds = peddata.GetNumPeds(); + _xCor = peddata.GetXCor(); + _yCor = peddata.GetYCor(); + _measureAreaId = boost::lexical_cast<string>(_areaForMethod_B->_id); + _tIn = new int[_NumPeds]; // Record the time of each pedestrian entering measurement area + _tOut = new int[_NumPeds]; + for (int i=0; i<_NumPeds; i++) + { + _tIn[i] = 0; + _tOut[i] = 0; + } + GetTinTout(peddata.GetNumFrames()); + Log->Write("------------------------Analyzing with Method B-----------------------------"); + if(_areaForMethod_B->_length<0) + { + Log->Write("Error:\tThe measurement area length for method B is not assigned!"); + exit(0); + } + else + { + GetFundamentalTinTout(_DensityPerFrame,_areaForMethod_B->_length); + } + delete []_tIn; + delete []_tOut; + return true; } void Method_B::GetTinTout(int numFrames) { - bool* IsinMeasurezone = new bool[_NumPeds]; - for (int i=0; i<_NumPeds; i++) - { - IsinMeasurezone[i] = false; - } - _DensityPerFrame = new double[numFrames]; - Method_C method_C; - for(int frameNr = 0; frameNr < numFrames; frameNr++ ) - { - vector<int> ids=_peds_t[frameNr]; - int pedsinMeasureArea=0; - for(unsigned int i=0; i< ids.size(); i++) - { - int ID = ids[i]; - int x = _xCor[ID][frameNr]; - int y = _yCor[ID][frameNr]; - if(within(make<point_2d>( (x), (y)), _areaForMethod_B->_poly)) - { - pedsinMeasureArea++; - } - if(within(make<point_2d>( (x), (y)), _areaForMethod_B->_poly)&&!(IsinMeasurezone[ID])) { - _tIn[ID]=frameNr; - IsinMeasurezone[ID] = true; - } - if((!within(make<point_2d>( (x), (y)), _areaForMethod_B->_poly))&&IsinMeasurezone[ID]) { - _tOut[ID]=frameNr; - IsinMeasurezone[ID] = false; - } - } - _DensityPerFrame[frameNr] = pedsinMeasureArea/(area(_areaForMethod_B->_poly)*CMtoM*CMtoM); - } - delete []IsinMeasurezone; + bool* IsinMeasurezone = new bool[_NumPeds]; + for (int i=0; i<_NumPeds; i++) + { + IsinMeasurezone[i] = false; + } + _DensityPerFrame = new double[numFrames]; + Method_C method_C; + for(int frameNr = 0; frameNr < numFrames; frameNr++ ) + { + vector<int> ids=_peds_t[frameNr]; + int pedsinMeasureArea=0; + for(unsigned int i=0; i< ids.size(); i++) + { + int ID = ids[i]; + int x = _xCor[ID][frameNr]; + int y = _yCor[ID][frameNr]; + if(within(make<point_2d>( (x), (y)), _areaForMethod_B->_poly)) + { + pedsinMeasureArea++; + } + if(within(make<point_2d>( (x), (y)), _areaForMethod_B->_poly)&&!(IsinMeasurezone[ID])) { + _tIn[ID]=frameNr; + IsinMeasurezone[ID] = true; + } + if((!within(make<point_2d>( (x), (y)), _areaForMethod_B->_poly))&&IsinMeasurezone[ID]) { + _tOut[ID]=frameNr; + IsinMeasurezone[ID] = false; + } + } + _DensityPerFrame[frameNr] = pedsinMeasureArea/(area(_areaForMethod_B->_poly)*CMtoM*CMtoM); + } + delete []IsinMeasurezone; } void Method_B::GetFundamentalTinTout(double *DensityPerFrame,double LengthMeasurementarea) { - FILE *fFD_TinTout; - Log->Write("---------Fundamental diagram from Method B will be calculated!------------------"); - string fdTinTout=_projectRootDir+"./Output/Fundamental_Diagram/TinTout/FDTinTout_"+_trajName+"_id_"+_measureAreaId+".dat";; - if((fFD_TinTout=Analysis::CreateFile(fdTinTout))==NULL) - { - Log->Write("cannot open the file to write the TinTout data\n"); - exit(0); - } - fprintf(fFD_TinTout,"#person Index\t density_i(m^(-2))\t velocity_i(m/s)\n"); - for(int i=0; i<_NumPeds; i++) - { - double velocity_temp=_fps*CMtoM*LengthMeasurementarea/(_tOut[i]-_tIn[i]); - double density_temp=0; - for(int j=_tIn[i]; j<_tOut[i]; j++) - { - density_temp+=DensityPerFrame[j]; - } - density_temp/=(_tOut[i]-_tIn[i]); - fprintf(fFD_TinTout,"%d\t%f\t%f\n",i+1,density_temp,velocity_temp); - } - fclose(fFD_TinTout); + FILE *fFD_TinTout; + Log->Write("---------Fundamental diagram from Method B will be calculated!------------------"); + string fdTinTout=_projectRootDir+"./Output/Fundamental_Diagram/TinTout/FDTinTout_"+_trajName+"_id_"+_measureAreaId+".dat";; + if((fFD_TinTout=Analysis::CreateFile(fdTinTout))==NULL) + { + Log->Write("cannot open the file to write the TinTout data\n"); + exit(0); + } + fprintf(fFD_TinTout,"#person Index\t density_i(m^(-2))\t velocity_i(m/s)\n"); + for(int i=0; i<_NumPeds; i++) + { + double velocity_temp=_fps*CMtoM*LengthMeasurementarea/(_tOut[i]-_tIn[i]); + double density_temp=0; + for(int j=_tIn[i]; j<_tOut[i]; j++) + { + density_temp+=DensityPerFrame[j]; + } + density_temp/=(_tOut[i]-_tIn[i]); + fprintf(fFD_TinTout,"%d\t%f\t%f\n",i+1,density_temp,velocity_temp); + } + fclose(fFD_TinTout); } void Method_B::SetMeasurementArea (MeasurementArea_B* area) { - _areaForMethod_B = area; + _areaForMethod_B = area; } diff --git a/methods/Method_D.h b/methods/Method_D.h index 7c1db75a1aa17b7e12436b21d724db8300a086db..a53d1f15163a5ff763d46c3bdce9fe9f8b39c0e2 100644 --- a/methods/Method_D.h +++ b/methods/Method_D.h @@ -85,14 +85,14 @@ private: bool OpenFileIndividualFD(); std::vector<polygon_2d> GetPolygons(std::vector<int> ids, std::vector<double>& XInFrame, std::vector<double>& YInFrame, - std::vector<double>& VInFrame, std::vector<int>& IdInFrame); + std::vector<double>& VInFrame, std::vector<int>& IdInFrame); void OutputVoronoiResults(const std::vector<polygon_2d>& polygons, int frid, const std::vector<double>& VInFrame); double GetVoronoiDensity(const std::vector<polygon_2d>& polygon, const polygon_2d & measureArea); double GetVoronoiDensity2(const std::vector<polygon_2d>& polygon, double* XInFrame, double* YInFrame, const polygon_2d& measureArea); double GetVoronoiVelocity(const std::vector<polygon_2d>& polygon, const std::vector<double>& Velocity, const polygon_2d & measureArea); void GetProfiles(const std::string& frameId, const std::vector<polygon_2d>& polygons, const std::vector<double>& velocity); void OutputVoroGraph(const std::string & frameId, std::vector<polygon_2d>& polygons, int numPedsInFrame,std::vector<double>& XInFrame, - std::vector<double>& YInFrame,const std::vector<double>& VInFrame); + std::vector<double>& YInFrame,const std::vector<double>& VInFrame); void GetIndividualFD(const std::vector<polygon_2d>& polygon, const std::vector<double>& Velocity, const std::vector<int>& Id, const polygon_2d& measureArea, int frid); /** diff --git a/methods/PedData.cpp b/methods/PedData.cpp index 46377286c3b8f8202d5491f4e58b3a5a0fb3ec0b..20dbb125f81b952aee72e77febba1864492e8db7 100644 --- a/methods/PedData.cpp +++ b/methods/PedData.cpp @@ -144,8 +144,8 @@ bool PedData::InitializeVariables(const string& filename) Ids_temp.erase(unique(Ids_temp.begin(), Ids_temp.end()), Ids_temp.end()); if((unsigned)_numPeds!=Ids_temp.size()) { - Log->Write("Error:\tThe index of ped ID is not continuous. Please modify the trajectory file!"); - return false; + Log->Write("Error:\tThe index of ped ID is not continuous. Please modify the trajectory file!"); + return false; } CreateGlobalVariables(_numPeds, _numFrames); @@ -287,7 +287,7 @@ bool PedData::InitializeVariables(TiXmlElement* xRootNode) vector<double> PedData::GetVInFrame(int frame, const vector<int>& ids) const { - vector<double> VInFrame; + vector<double> VInFrame; for(unsigned int i=0; i<ids.size();i++) { int id = ids[i]; @@ -301,17 +301,17 @@ vector<double> PedData::GetVInFrame(int frame, const vector<int>& ids) const vector<double> PedData::GetXInFrame(int frame, const vector<int>& ids) const { - vector<double> XInFrame; - for(int id:ids) - { - XInFrame.push_back(_xCor[id][frame]); - } + vector<double> XInFrame; + for(int id:ids) + { + XInFrame.push_back(_xCor[id][frame]); + } return XInFrame; } vector<double> PedData::GetYInFrame(int frame, const vector<int>& ids) const { - vector<double> YInFrame; + vector<double> YInFrame; for(unsigned int i=0; i<ids.size();i++) { int id = ids[i]; @@ -322,11 +322,11 @@ vector<double> PedData::GetYInFrame(int frame, const vector<int>& ids) const vector<int> PedData::GetIdInFrame(const vector<int>& ids) const { - vector<int> IdInFrame; + vector<int> IdInFrame; for(int id:ids) { - id = id +_minID; - IdInFrame.push_back(id); + id = id +_minID; + IdInFrame.push_back(id); } return IdInFrame; } @@ -425,7 +425,7 @@ int PedData::GetNumFrames() const int PedData::GetNumPeds() const { - return _numPeds; + return _numPeds; } int PedData::GetFps() const @@ -445,23 +445,23 @@ map<int , vector<int> > PedData::GetPedsFrame() const double** PedData::GetXCor() const { - return _xCor; + return _xCor; } double** PedData::GetYCor() const { - return _yCor; + return _yCor; } int* PedData::GetFirstFrame() const { - return _firstFrame; + return _firstFrame; } int* PedData::GetLastFrame() const { - return _lastFrame; + return _lastFrame; } string PedData::GetProjectRootDir() const { - return _projectRootDir; + return _projectRootDir; } diff --git a/methods/VoronoiDiagram.cpp b/methods/VoronoiDiagram.cpp index dec9ff1558aa6d836f62ba45ab895d5cf6aa5d4f..c5786bd36a3e91a3a2c7da07f3265dd9e80204f5 100644 --- a/methods/VoronoiDiagram.cpp +++ b/methods/VoronoiDiagram.cpp @@ -47,453 +47,453 @@ VoronoiDiagram::~VoronoiDiagram() // Traversing Voronoi edges using cell iterator. vector<polygon_2d> VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame, vector<double>& YInFrame, - vector<double>& VInFrame, vector<int>& IdInFrame, const int numPedsInFrame, const double Bound_Max) + vector<double>& VInFrame, vector<int>& IdInFrame, const int numPedsInFrame, const double Bound_Max) { - int XInFrame_temp[numPedsInFrame]; - int YInFrame_temp[numPedsInFrame]; - double VInFrame_temp[numPedsInFrame]; - int IdInFrame_temp[numPedsInFrame]; - - for (int i = 0; i < numPedsInFrame; i++) - { - points.push_back(point_type2(round(XInFrame[i]), round(YInFrame[i]))); - XInFrame_temp[i] = round(XInFrame[i]); - YInFrame_temp[i] = round(YInFrame[i]); - VInFrame_temp[i] = VInFrame[i]; - IdInFrame_temp[i] = IdInFrame[i]; - } - - VD voronoidiagram; - construct_voronoi(points.begin(), points.end(), &voronoidiagram); - int Ncell = 0; - std::vector<polygon_2d> polygons; - - double Bd_Box_minX = -Bound_Max; - double Bd_Box_minY = -Bound_Max; - double Bd_Box_maxX = Bound_Max; - double Bd_Box_maxY = Bound_Max; - polygon_2d boundingbox; - boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(-Bound_Max, -Bound_Max)); - boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(-Bound_Max, Bound_Max)); - boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(Bound_Max, Bound_Max)); - boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(Bound_Max, -Bound_Max)); - boost::geometry::correct(boundingbox); - - for (voronoi_diagram<double>::const_cell_iterator it = voronoidiagram.cells().begin(); - it != voronoidiagram.cells().end(); ++it) - { - polygon_2d poly; - vector<point_type2> polypts; - point_type2 pt_s; - point_type2 pt_e; - const voronoi_diagram<double>::cell_type& cell = *it; - const voronoi_diagram<double>::edge_type* edge = cell.incident_edge(); - point_type2 pt_temp; - point_type2 thispoint = retrieve_point(*it); - for (int k = 0; k < numPedsInFrame; k++) - { - if (XInFrame_temp[k] == thispoint.x() && YInFrame_temp[k] == thispoint.y()) - { - VInFrame[Ncell] = VInFrame_temp[k]; - IdInFrame[Ncell] = IdInFrame_temp[k]; - break; - } - } - XInFrame[Ncell] = thispoint.x(); - YInFrame[Ncell] = thispoint.y(); - int NumVertex = 0; - int index_end = 0; - bool infinite_s = false; - bool infinite_e = false; - // This is convenient way to iterate edges around Voronoi cell. - - do - { - if (edge->vertex0()) - { - if(edge->vertex1()) - { - if(fabs(edge->vertex0()->x()) < Bound_Max && (fabs(edge->vertex0()->y()) < Bound_Max) - &&fabs(edge->vertex1()->x()) < Bound_Max && (fabs(edge->vertex1()->y()) < Bound_Max)) - { - polypts.push_back(point_type2(edge->vertex0()->x(), edge->vertex0()->y())); - pt_temp = point_type2(edge->vertex0()->x(), edge->vertex0()->y()); - NumVertex++; - } - else if((fabs(edge->vertex0()->x()) > Bound_Max || (fabs(edge->vertex0()->y()) > Bound_Max)) - &&fabs(edge->vertex1()->x()) < Bound_Max && (fabs(edge->vertex1()->y()) < Bound_Max)) - { - pt_s = getIntersectionPoint(point_2d(edge->vertex0()->x(), edge->vertex0()->y()), point_2d(edge->vertex1()->x(), edge->vertex1()->y()), boundingbox); - polypts.push_back(point_type2(pt_s.x(), pt_s.y())); - NumVertex++; - infinite_s = true; - } - else if((fabs(edge->vertex0()->x()) < Bound_Max && (fabs(edge->vertex0()->y()) < Bound_Max)) - &&(fabs(edge->vertex1()->x()) > Bound_Max || (fabs(edge->vertex1()->y()) > Bound_Max))) - { - polypts.push_back(point_type2(edge->vertex0()->x(), edge->vertex0()->y())); - pt_e = getIntersectionPoint(point_2d(edge->vertex0()->x(), edge->vertex0()->y()), point_2d(edge->vertex1()->x(), edge->vertex1()->y()), boundingbox); - polypts.push_back(point_type2(pt_e.x(), pt_e.y())); - NumVertex+=2; - index_end = NumVertex; - infinite_e = true; - } - } - else - { - if(fabs(edge->vertex0()->x()) < Bound_Max && (fabs(edge->vertex0()->y()) < Bound_Max)) - { - polypts.push_back(point_type2(edge->vertex0()->x(), edge->vertex0()->y())); - pt_e = clip_infinite_edge(*edge, Bd_Box_minX, Bd_Box_minY, Bd_Box_maxX, - Bd_Box_maxY); - polypts.push_back(point_type2(pt_e.x(), pt_e.y())); - NumVertex+=2; - index_end = NumVertex; - infinite_e = true; - } - } - } - else - { - if(edge->vertex1()&&fabs(edge->vertex1()->x()) < Bound_Max && (fabs(edge->vertex1()->y()) < Bound_Max)) - { - pt_s = clip_infinite_edge(*edge, Bd_Box_minX, Bd_Box_minY, Bd_Box_maxX, - Bd_Box_maxY); - polypts.push_back(point_type2(pt_s.x(), pt_s.y())); - NumVertex++; - infinite_s = true; - } - } - edge = edge->next(); - - } while (edge != cell.incident_edge()); - Ncell++; - - if (infinite_s && infinite_e) - { - vector<point_type2> vertexes = add_bounding_points(pt_s, pt_e, pt_temp, - Bd_Box_minX, Bd_Box_minY, Bd_Box_maxX, Bd_Box_maxY); - - if (!vertexes.empty()) - { - polypts.reserve(polypts.size() + vertexes.size()); - polypts.insert(polypts.begin() + index_end, vertexes.begin(), - vertexes.end()); - } - } - for (int i = 0; i < (int) polypts.size(); i++) - { - boost::geometry::append(poly, point_2d(polypts[i].x(), polypts[i].y())); - } - - boost::geometry::correct(poly); - polygons.push_back(poly); - } - return polygons; + int XInFrame_temp[numPedsInFrame]; + int YInFrame_temp[numPedsInFrame]; + double VInFrame_temp[numPedsInFrame]; + int IdInFrame_temp[numPedsInFrame]; + + for (int i = 0; i < numPedsInFrame; i++) + { + points.push_back(point_type2(round(XInFrame[i]), round(YInFrame[i]))); + XInFrame_temp[i] = round(XInFrame[i]); + YInFrame_temp[i] = round(YInFrame[i]); + VInFrame_temp[i] = VInFrame[i]; + IdInFrame_temp[i] = IdInFrame[i]; + } + + VD voronoidiagram; + construct_voronoi(points.begin(), points.end(), &voronoidiagram); + int Ncell = 0; + std::vector<polygon_2d> polygons; + + double Bd_Box_minX = -Bound_Max; + double Bd_Box_minY = -Bound_Max; + double Bd_Box_maxX = Bound_Max; + double Bd_Box_maxY = Bound_Max; + polygon_2d boundingbox; + boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(-Bound_Max, -Bound_Max)); + boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(-Bound_Max, Bound_Max)); + boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(Bound_Max, Bound_Max)); + boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(Bound_Max, -Bound_Max)); + boost::geometry::correct(boundingbox); + + for (voronoi_diagram<double>::const_cell_iterator it = voronoidiagram.cells().begin(); + it != voronoidiagram.cells().end(); ++it) + { + polygon_2d poly; + vector<point_type2> polypts; + point_type2 pt_s; + point_type2 pt_e; + const voronoi_diagram<double>::cell_type& cell = *it; + const voronoi_diagram<double>::edge_type* edge = cell.incident_edge(); + point_type2 pt_temp; + point_type2 thispoint = retrieve_point(*it); + for (int k = 0; k < numPedsInFrame; k++) + { + if (XInFrame_temp[k] == thispoint.x() && YInFrame_temp[k] == thispoint.y()) + { + VInFrame[Ncell] = VInFrame_temp[k]; + IdInFrame[Ncell] = IdInFrame_temp[k]; + break; + } + } + XInFrame[Ncell] = thispoint.x(); + YInFrame[Ncell] = thispoint.y(); + int NumVertex = 0; + int index_end = 0; + bool infinite_s = false; + bool infinite_e = false; + // This is convenient way to iterate edges around Voronoi cell. + + do + { + if (edge->vertex0()) + { + if(edge->vertex1()) + { + if(fabs(edge->vertex0()->x()) < Bound_Max && (fabs(edge->vertex0()->y()) < Bound_Max) + &&fabs(edge->vertex1()->x()) < Bound_Max && (fabs(edge->vertex1()->y()) < Bound_Max)) + { + polypts.push_back(point_type2(edge->vertex0()->x(), edge->vertex0()->y())); + pt_temp = point_type2(edge->vertex0()->x(), edge->vertex0()->y()); + NumVertex++; + } + else if((fabs(edge->vertex0()->x()) > Bound_Max || (fabs(edge->vertex0()->y()) > Bound_Max)) + &&fabs(edge->vertex1()->x()) < Bound_Max && (fabs(edge->vertex1()->y()) < Bound_Max)) + { + pt_s = getIntersectionPoint(point_2d(edge->vertex0()->x(), edge->vertex0()->y()), point_2d(edge->vertex1()->x(), edge->vertex1()->y()), boundingbox); + polypts.push_back(point_type2(pt_s.x(), pt_s.y())); + NumVertex++; + infinite_s = true; + } + else if((fabs(edge->vertex0()->x()) < Bound_Max && (fabs(edge->vertex0()->y()) < Bound_Max)) + &&(fabs(edge->vertex1()->x()) > Bound_Max || (fabs(edge->vertex1()->y()) > Bound_Max))) + { + polypts.push_back(point_type2(edge->vertex0()->x(), edge->vertex0()->y())); + pt_e = getIntersectionPoint(point_2d(edge->vertex0()->x(), edge->vertex0()->y()), point_2d(edge->vertex1()->x(), edge->vertex1()->y()), boundingbox); + polypts.push_back(point_type2(pt_e.x(), pt_e.y())); + NumVertex+=2; + index_end = NumVertex; + infinite_e = true; + } + } + else + { + if(fabs(edge->vertex0()->x()) < Bound_Max && (fabs(edge->vertex0()->y()) < Bound_Max)) + { + polypts.push_back(point_type2(edge->vertex0()->x(), edge->vertex0()->y())); + pt_e = clip_infinite_edge(*edge, Bd_Box_minX, Bd_Box_minY, Bd_Box_maxX, + Bd_Box_maxY); + polypts.push_back(point_type2(pt_e.x(), pt_e.y())); + NumVertex+=2; + index_end = NumVertex; + infinite_e = true; + } + } + } + else + { + if(edge->vertex1()&&fabs(edge->vertex1()->x()) < Bound_Max && (fabs(edge->vertex1()->y()) < Bound_Max)) + { + pt_s = clip_infinite_edge(*edge, Bd_Box_minX, Bd_Box_minY, Bd_Box_maxX, + Bd_Box_maxY); + polypts.push_back(point_type2(pt_s.x(), pt_s.y())); + NumVertex++; + infinite_s = true; + } + } + edge = edge->next(); + + } while (edge != cell.incident_edge()); + Ncell++; + + if (infinite_s && infinite_e) + { + vector<point_type2> vertexes = add_bounding_points(pt_s, pt_e, pt_temp, + Bd_Box_minX, Bd_Box_minY, Bd_Box_maxX, Bd_Box_maxY); + + if (!vertexes.empty()) + { + polypts.reserve(polypts.size() + vertexes.size()); + polypts.insert(polypts.begin() + index_end, vertexes.begin(), + vertexes.end()); + } + } + for (int i = 0; i < (int) polypts.size(); i++) + { + boost::geometry::append(poly, point_2d(polypts[i].x(), polypts[i].y())); + } + + boost::geometry::correct(poly); + polygons.push_back(poly); + } + return polygons; } point_type2 VoronoiDiagram::retrieve_point(const cell_type& cell) { - source_index_type index = cell.source_index(); - return points[index]; + source_index_type index = cell.source_index(); + return points[index]; } point_type2 VoronoiDiagram::clip_infinite_edge(const edge_type& edge, double minX, double minY, - double maxX, double maxY) + double maxX, double maxY) { - const cell_type& cell1 = *edge.cell(); - const cell_type& cell2 = *edge.twin()->cell(); - point_type2 origin, direction, pt; - // Infinite edges could not be created by two segment sites. - - if (cell1.contains_point() && cell2.contains_point()) - { - point_type2 p1 = retrieve_point(cell1); - point_type2 p2 = retrieve_point(cell2); - origin.x((p1.x() + p2.x()) * 0.5); - origin.y((p1.y() + p2.y()) * 0.5); - direction.x(p1.y() - p2.y()); - direction.y(p2.x() - p1.x()); - } - else - { - cout<<"A cell cannot find its center point!"<<endl; - exit(0); - } - - double side = maxX - minX; - double koef = side / (std::max)(fabs(direction.x()), fabs(direction.y())); - if (edge.vertex0() == NULL) - { - pt.x(origin.x() - direction.x() * koef); - pt.y(origin.y() - direction.y() * koef); - } - else if (edge.vertex1() == NULL) - { - pt.x(origin.x() + direction.x() * koef); - pt.y(origin.y() + direction.y() * koef); - } - - polygon_2d boundingbox; - boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(minX, minY)); - boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(minX, maxY)); - boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(maxX, maxY)); - boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(maxX, minY)); - boost::geometry::correct(boundingbox); - - pt = getIntersectionPoint(point_2d(pt.x(), pt.y()), point_2d(origin.x(), origin.y()), boundingbox); - return pt; + const cell_type& cell1 = *edge.cell(); + const cell_type& cell2 = *edge.twin()->cell(); + point_type2 origin, direction, pt; + // Infinite edges could not be created by two segment sites. + + if (cell1.contains_point() && cell2.contains_point()) + { + point_type2 p1 = retrieve_point(cell1); + point_type2 p2 = retrieve_point(cell2); + origin.x((p1.x() + p2.x()) * 0.5); + origin.y((p1.y() + p2.y()) * 0.5); + direction.x(p1.y() - p2.y()); + direction.y(p2.x() - p1.x()); + } + else + { + cout<<"A cell cannot find its center point!"<<endl; + exit(0); + } + + double side = maxX - minX; + double koef = side / (std::max)(fabs(direction.x()), fabs(direction.y())); + if (edge.vertex0() == NULL) + { + pt.x(origin.x() - direction.x() * koef); + pt.y(origin.y() - direction.y() * koef); + } + else if (edge.vertex1() == NULL) + { + pt.x(origin.x() + direction.x() * koef); + pt.y(origin.y() + direction.y() * koef); + } + + polygon_2d boundingbox; + boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(minX, minY)); + boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(minX, maxY)); + boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(maxX, maxY)); + boost::geometry::append(boundingbox, boost::geometry::make<point_2d>(maxX, minY)); + boost::geometry::correct(boundingbox); + + pt = getIntersectionPoint(point_2d(pt.x(), pt.y()), point_2d(origin.x(), origin.y()), boundingbox); + return pt; } double VoronoiDiagram::area_triangle(const point_type2& tri_p1, const point_type2& tri_p2, const point_type2& tri_p3) { - double x = tri_p2.x() - tri_p1.x(); - double y = tri_p2.y() - tri_p1.y(); + double x = tri_p2.x() - tri_p1.x(); + double y = tri_p2.y() - tri_p1.y(); - double x2 = tri_p3.x() - tri_p1.x(); - double y2 = tri_p3.y() - tri_p1.y(); + double x2 = tri_p3.x() - tri_p1.x(); + double y2 = tri_p3.y() - tri_p1.y(); - return abs(x * y2 - x2 * y) / 2.; + return abs(x * y2 - x2 * y) / 2.; } bool VoronoiDiagram::point_inside_triangle(const point_type2& pt, const point_type2& tri_p1, const point_type2& tri_p2, const point_type2& tri_p3) { - double s = area_triangle(tri_p1, tri_p2, tri_p3); - double ss = area_triangle(pt, tri_p1, tri_p2); - double ss2 = area_triangle(pt, tri_p2, tri_p3); - double ss3 = area_triangle(pt, tri_p1, tri_p3); - - if (ss + ss2 + ss3 - s > 1e-6) - { - return false; - } - return true; + double s = area_triangle(tri_p1, tri_p2, tri_p3); + double ss = area_triangle(pt, tri_p1, tri_p2); + double ss2 = area_triangle(pt, tri_p2, tri_p3); + double ss3 = area_triangle(pt, tri_p1, tri_p3); + + if (ss + ss2 + ss3 - s > 1e-6) + { + return false; + } + return true; } vector<point_type2> VoronoiDiagram::add_bounding_points(const point_type2& pt1, const point_type2& pt2, const point_type2& pt, - double minX, double minY, double maxX, double maxY) + double minX, double minY, double maxX, double maxY) { - vector<point_type2> bounding_vertex; - if (fabs(pt1.x() - pt2.x()) > J_EPS && fabs(pt1.y() - pt2.y()) > J_EPS) - { - if ((fabs(pt1.y() - maxY) < J_EPS && fabs(pt2.x() - maxX) < J_EPS) - || (fabs(pt2.y() - maxY) < J_EPS && fabs(pt1.x() - maxX) < J_EPS)) - { - point_type2 vertex(maxX, maxY); - if (point_inside_triangle(pt, vertex, pt1, pt2)) - { - bounding_vertex.push_back(point_type2(minX, maxY)); - bounding_vertex.push_back(point_type2(minX, minY)); - bounding_vertex.push_back(point_type2(maxX, minY)); - } - else - { - bounding_vertex.push_back(point_type2(maxX, maxY)); - } - } - else if ((fabs(pt1.y() - maxY) < J_EPS && fabs(pt2.x() - minX) < J_EPS) - || (fabs(pt2.y() - maxY) < J_EPS && fabs(pt1.x() - minX) < J_EPS)) - { - point_type2 vertex(minX, maxY); - if (point_inside_triangle(pt, vertex, pt1, pt2)) - { - bounding_vertex.push_back(point_type2(minX, minY)); - bounding_vertex.push_back(point_type2(maxX, minY)); - bounding_vertex.push_back(point_type2(maxX, maxY)); - } - else - { - bounding_vertex.push_back(point_type2(minX, maxY)); - } - } - else if ((fabs(pt1.y() - minY) < J_EPS && fabs(pt2.x() - minX) < J_EPS) - || (fabs(pt2.y() - minY) < J_EPS && fabs(pt1.x() - minX) < J_EPS)) - { - point_type2 vertex(minX, minY); - if (point_inside_triangle(pt, vertex, pt1, pt2)) - { - bounding_vertex.push_back(point_type2(maxX, minY)); - bounding_vertex.push_back(point_type2(maxX, maxY)); - bounding_vertex.push_back(point_type2(minX, maxY)); - } - else - { - bounding_vertex.push_back(point_type2(minX, minY)); - } - } - else if ((fabs(pt1.y() - minY) < J_EPS && fabs(pt2.x() - maxX) < J_EPS) - || (fabs(pt2.y() - minY) < J_EPS && fabs(pt1.x() - maxX) < J_EPS)) - { - point_type2 vertex(maxX, minY); - if (point_inside_triangle(pt, vertex, pt1, pt2)) - { - bounding_vertex.push_back(point_type2(maxX, maxY)); - bounding_vertex.push_back(point_type2(minX, maxY)); - bounding_vertex.push_back(point_type2(minX, minY)); - } - else - { - bounding_vertex.push_back(point_type2(maxX, minY)); - } - } - else if ((fabs(pt1.y() - minY) < J_EPS && fabs(pt2.y() - maxY) < J_EPS) - || (fabs(pt2.y() - minY) < J_EPS && fabs(pt1.y() - maxY) < J_EPS)) - { - if (fabs(pt1.x() - pt2.x()) < J_EPS) - { - if (pt1.x() < pt.x()) - { - bounding_vertex.push_back(point_type2(minX, maxY)); - bounding_vertex.push_back(point_type2(minX, minY)); - } - else - { - bounding_vertex.push_back(point_type2(maxX, minY)); - bounding_vertex.push_back(point_type2(maxX, maxY)); - } - } - else - { - double tempx = pt1.x() - + (pt2.x() - pt1.x()) * (pt.y() - pt1.y()) / (pt2.y() - pt1.y()); - if (tempx < pt.x()) - { - bounding_vertex.push_back(point_type2(minX, maxY)); - bounding_vertex.push_back(point_type2(minX, minY)); - } - else - { - bounding_vertex.push_back(point_type2(maxX, minY)); - bounding_vertex.push_back(point_type2(maxX, maxY)); - } - } - } - else if ((fabs(pt1.x() - minX) < J_EPS && fabs(pt2.x() - maxX) < J_EPS) - || (fabs(pt2.x() - minX) < J_EPS && fabs(pt1.x() - maxX) < J_EPS)) - { - if (fabs(pt1.y() - pt2.y()) < J_EPS) - { - if (pt1.y() < pt.y()) - { - bounding_vertex.push_back(point_type2(minX, minY)); - bounding_vertex.push_back(point_type2(maxX, minY)); - } - else - { - bounding_vertex.push_back(point_type2(maxX, maxY)); - bounding_vertex.push_back(point_type2(minX, maxY)); - } - } - else - { - double tempy = pt1.y() - + (pt2.y() - pt1.y()) * (pt.x() - pt1.x()) / (pt2.x() - pt1.x()); - if (tempy < pt.y()) - { - bounding_vertex.push_back(point_type2(minX, minY)); - bounding_vertex.push_back(point_type2(maxX, minY)); - } - else - { - bounding_vertex.push_back(point_type2(maxX, maxY)); - bounding_vertex.push_back(point_type2(minX, maxY)); - } - } - } - } - return bounding_vertex; + vector<point_type2> bounding_vertex; + if (fabs(pt1.x() - pt2.x()) > J_EPS && fabs(pt1.y() - pt2.y()) > J_EPS) + { + if ((fabs(pt1.y() - maxY) < J_EPS && fabs(pt2.x() - maxX) < J_EPS) + || (fabs(pt2.y() - maxY) < J_EPS && fabs(pt1.x() - maxX) < J_EPS)) + { + point_type2 vertex(maxX, maxY); + if (point_inside_triangle(pt, vertex, pt1, pt2)) + { + bounding_vertex.push_back(point_type2(minX, maxY)); + bounding_vertex.push_back(point_type2(minX, minY)); + bounding_vertex.push_back(point_type2(maxX, minY)); + } + else + { + bounding_vertex.push_back(point_type2(maxX, maxY)); + } + } + else if ((fabs(pt1.y() - maxY) < J_EPS && fabs(pt2.x() - minX) < J_EPS) + || (fabs(pt2.y() - maxY) < J_EPS && fabs(pt1.x() - minX) < J_EPS)) + { + point_type2 vertex(minX, maxY); + if (point_inside_triangle(pt, vertex, pt1, pt2)) + { + bounding_vertex.push_back(point_type2(minX, minY)); + bounding_vertex.push_back(point_type2(maxX, minY)); + bounding_vertex.push_back(point_type2(maxX, maxY)); + } + else + { + bounding_vertex.push_back(point_type2(minX, maxY)); + } + } + else if ((fabs(pt1.y() - minY) < J_EPS && fabs(pt2.x() - minX) < J_EPS) + || (fabs(pt2.y() - minY) < J_EPS && fabs(pt1.x() - minX) < J_EPS)) + { + point_type2 vertex(minX, minY); + if (point_inside_triangle(pt, vertex, pt1, pt2)) + { + bounding_vertex.push_back(point_type2(maxX, minY)); + bounding_vertex.push_back(point_type2(maxX, maxY)); + bounding_vertex.push_back(point_type2(minX, maxY)); + } + else + { + bounding_vertex.push_back(point_type2(minX, minY)); + } + } + else if ((fabs(pt1.y() - minY) < J_EPS && fabs(pt2.x() - maxX) < J_EPS) + || (fabs(pt2.y() - minY) < J_EPS && fabs(pt1.x() - maxX) < J_EPS)) + { + point_type2 vertex(maxX, minY); + if (point_inside_triangle(pt, vertex, pt1, pt2)) + { + bounding_vertex.push_back(point_type2(maxX, maxY)); + bounding_vertex.push_back(point_type2(minX, maxY)); + bounding_vertex.push_back(point_type2(minX, minY)); + } + else + { + bounding_vertex.push_back(point_type2(maxX, minY)); + } + } + else if ((fabs(pt1.y() - minY) < J_EPS && fabs(pt2.y() - maxY) < J_EPS) + || (fabs(pt2.y() - minY) < J_EPS && fabs(pt1.y() - maxY) < J_EPS)) + { + if (fabs(pt1.x() - pt2.x()) < J_EPS) + { + if (pt1.x() < pt.x()) + { + bounding_vertex.push_back(point_type2(minX, maxY)); + bounding_vertex.push_back(point_type2(minX, minY)); + } + else + { + bounding_vertex.push_back(point_type2(maxX, minY)); + bounding_vertex.push_back(point_type2(maxX, maxY)); + } + } + else + { + double tempx = pt1.x() + + (pt2.x() - pt1.x()) * (pt.y() - pt1.y()) / (pt2.y() - pt1.y()); + if (tempx < pt.x()) + { + bounding_vertex.push_back(point_type2(minX, maxY)); + bounding_vertex.push_back(point_type2(minX, minY)); + } + else + { + bounding_vertex.push_back(point_type2(maxX, minY)); + bounding_vertex.push_back(point_type2(maxX, maxY)); + } + } + } + else if ((fabs(pt1.x() - minX) < J_EPS && fabs(pt2.x() - maxX) < J_EPS) + || (fabs(pt2.x() - minX) < J_EPS && fabs(pt1.x() - maxX) < J_EPS)) + { + if (fabs(pt1.y() - pt2.y()) < J_EPS) + { + if (pt1.y() < pt.y()) + { + bounding_vertex.push_back(point_type2(minX, minY)); + bounding_vertex.push_back(point_type2(maxX, minY)); + } + else + { + bounding_vertex.push_back(point_type2(maxX, maxY)); + bounding_vertex.push_back(point_type2(minX, maxY)); + } + } + else + { + double tempy = pt1.y() + + (pt2.y() - pt1.y()) * (pt.x() - pt1.x()) / (pt2.x() - pt1.x()); + if (tempy < pt.y()) + { + bounding_vertex.push_back(point_type2(minX, minY)); + bounding_vertex.push_back(point_type2(maxX, minY)); + } + else + { + bounding_vertex.push_back(point_type2(maxX, maxY)); + bounding_vertex.push_back(point_type2(minX, maxY)); + } + } + } + } + return bounding_vertex; } //-----------In getIntersectionPoint() the edges of the square is vertical or horizontal segment-------------- - point_type2 VoronoiDiagram::getIntersectionPoint(const point_2d& pt0, const point_2d& pt1, const polygon_2d& square) +point_type2 VoronoiDiagram::getIntersectionPoint(const point_2d& pt0, const point_2d& pt1, const polygon_2d& square) { - vector<point_2d> pt; - segment edge0(pt0, pt1); - vector<point_2d> const& points = square.outer(); - for (vector<point_2d>::size_type i = 1; i < points.size(); ++i) - { - segment edge1(points[i], points[i-1]); - if(intersects(edge0, edge1)) - { - intersection(edge0, edge1, pt); - break; - } - } - if(pt.empty()) - { - segment edge1(points[3], points[0]); - intersection(edge0, edge1, pt); - } - point_type2 interpts(pt[0].x(), pt[0].y()); - return interpts; + vector<point_2d> pt; + segment edge0(pt0, pt1); + vector<point_2d> const& points = square.outer(); + for (vector<point_2d>::size_type i = 1; i < points.size(); ++i) + { + segment edge1(points[i], points[i-1]); + if(intersects(edge0, edge1)) + { + intersection(edge0, edge1, pt); + break; + } + } + if(pt.empty()) + { + segment edge1(points[3], points[0]); + intersection(edge0, edge1, pt); + } + point_type2 interpts(pt[0].x(), pt[0].y()); + return interpts; } std::vector<polygon_2d> VoronoiDiagram::cutPolygonsWithGeometry(const std::vector<polygon_2d>& polygon, - const polygon_2d& Geometry, const vector<double>& xs, const vector<double>& ys) + const polygon_2d& Geometry, const vector<double>& xs, const vector<double>& ys) { - vector<polygon_2d> intersetionpolygons; - int temp = 0; - for (const auto & polygon_iterator:polygon) - { - polygon_list v; - intersection(Geometry, polygon_iterator, v); - - BOOST_FOREACH(auto & it, v) - { - if (within(point_2d(xs[temp], ys[temp]), it)) - { - //check and remove duplicates - //dispatch::unique (it); - polygon_2d simplified; - simplify(it,simplified,J_EPS); - - correct(simplified); - intersetionpolygons.push_back(simplified); - } - } - - temp++; - } - return intersetionpolygons; + vector<polygon_2d> intersetionpolygons; + int temp = 0; + for (const auto & polygon_iterator:polygon) + { + polygon_list v; + intersection(Geometry, polygon_iterator, v); + + BOOST_FOREACH(auto & it, v) + { + if (within(point_2d(xs[temp], ys[temp]), it)) + { + //check and remove duplicates + //dispatch::unique (it); + polygon_2d simplified; + simplify(it,simplified,J_EPS); + + correct(simplified); + intersetionpolygons.push_back(simplified); + } + } + + temp++; + } + return intersetionpolygons; } std::vector<polygon_2d> VoronoiDiagram::cutPolygonsWithCircle(const std::vector<polygon_2d>& polygon, - const vector<double>& xs, const vector<double>& ys, double radius, int edges) + const vector<double>& xs, const vector<double>& ys, double radius, int edges) { - std::vector<polygon_2d> intersetionpolygons; - int temp = 0; - for (const auto & polygon_iterator:polygon) - { - polygon_2d circle; - { - for (int angle = 0; angle <=edges; angle++) - { - double ptx= xs[temp] + radius * cos(angle * 2*PI / edges); - double pty= ys[temp] + radius * sin(angle * 2*PI / edges); - append(circle, make<point_2d>(ptx, pty)); - } - } - correct(circle); - - polygon_list v; - intersection(circle, polygon_iterator, v); - BOOST_FOREACH(auto & it, v) - { - if (within(point_2d(xs[temp], ys[temp]), it)) - { - //check and remove duplicates - //dispatch::unique (it); - polygon_2d simplified; - simplify(it,simplified,J_EPS); - correct(simplified); - intersetionpolygons.push_back(simplified); - } - } - temp++; - } - return intersetionpolygons; + std::vector<polygon_2d> intersetionpolygons; + int temp = 0; + for (const auto & polygon_iterator:polygon) + { + polygon_2d circle; + { + for (int angle = 0; angle <=edges; angle++) + { + double ptx= xs[temp] + radius * cos(angle * 2*PI / edges); + double pty= ys[temp] + radius * sin(angle * 2*PI / edges); + append(circle, make<point_2d>(ptx, pty)); + } + } + correct(circle); + + polygon_list v; + intersection(circle, polygon_iterator, v); + BOOST_FOREACH(auto & it, v) + { + if (within(point_2d(xs[temp], ys[temp]), it)) + { + //check and remove duplicates + //dispatch::unique (it); + polygon_2d simplified; + simplify(it,simplified,J_EPS); + correct(simplified); + intersetionpolygons.push_back(simplified); + } + } + temp++; + } + return intersetionpolygons; } diff --git a/methods/VoronoiDiagram.h b/methods/VoronoiDiagram.h index eef482a858c832aac1b7be793695a1f5dbcac2fb..d7006d29f87e8bfeda3918950a2c6909f44d9ddf 100644 --- a/methods/VoronoiDiagram.h +++ b/methods/VoronoiDiagram.h @@ -66,23 +66,23 @@ typedef VD::cell_type::source_index_type source_index_type; class VoronoiDiagram { private: - std::vector<point_type2> points; - point_type2 retrieve_point(const cell_type& cell); - point_type2 clip_infinite_edge( const edge_type& edge, double minX, double minY, double maxX, double maxY); - double area_triangle(const point_type2& tri_p1, const point_type2& tri_p2, const point_type2& tri_p3); - bool point_inside_triangle(const point_type2& pt, const point_type2& tri_p1, const point_type2& tri_p2, const point_type2& tri_p3); - std::vector<point_type2> add_bounding_points(const point_type2& pt1, const point_type2& pt2, const point_type2& pt, double minX, - double minY, double maxX, double maxY); - point_type2 getIntersectionPoint(const point_2d& pt0, const point_2d& pt1, const polygon_2d& square); + std::vector<point_type2> points; + point_type2 retrieve_point(const cell_type& cell); + point_type2 clip_infinite_edge( const edge_type& edge, double minX, double minY, double maxX, double maxY); + double area_triangle(const point_type2& tri_p1, const point_type2& tri_p2, const point_type2& tri_p3); + bool point_inside_triangle(const point_type2& pt, const point_type2& tri_p1, const point_type2& tri_p2, const point_type2& tri_p3); + std::vector<point_type2> add_bounding_points(const point_type2& pt1, const point_type2& pt2, const point_type2& pt, double minX, + double minY, double maxX, double maxY); + point_type2 getIntersectionPoint(const point_2d& pt0, const point_2d& pt1, const polygon_2d& square); public: - VoronoiDiagram(); - virtual ~VoronoiDiagram(); + VoronoiDiagram(); + virtual ~VoronoiDiagram(); - std::vector<polygon_2d> getVoronoiPolygons(std::vector<double>& XInFrame, std::vector<double>& YInFrame, std::vector<double>& VInFrame, - std::vector<int>& IdInFrame, const int numPedsInFrame, const double Bound_Max); - std::vector<polygon_2d> cutPolygonsWithGeometry(const std::vector<polygon_2d>& polygon, const polygon_2d& Geometry, const std::vector<double>& xs, const std::vector<double>& ys); - std::vector<polygon_2d> cutPolygonsWithCircle(const std::vector<polygon_2d>& polygon, const std::vector<double>& xs, const std::vector<double>& ys, double radius, int edges); + std::vector<polygon_2d> getVoronoiPolygons(std::vector<double>& XInFrame, std::vector<double>& YInFrame, std::vector<double>& VInFrame, + std::vector<int>& IdInFrame, const int numPedsInFrame, const double Bound_Max); + std::vector<polygon_2d> cutPolygonsWithGeometry(const std::vector<polygon_2d>& polygon, const polygon_2d& Geometry, const std::vector<double>& xs, const std::vector<double>& ys); + std::vector<polygon_2d> cutPolygonsWithCircle(const std::vector<polygon_2d>& polygon, const std::vector<double>& xs, const std::vector<double>& ys, double radius, int edges); }; #endif /* VORONOIDIAGRAM_H_ */