diff --git a/Analysis.cpp b/Analysis.cpp
index 17e4b904425e99eb5e5cbbb1868067b8d312b149..60ec457f1509c2d75f758135e864a53aec72a79e 100644
--- a/Analysis.cpp
+++ b/Analysis.cpp
@@ -161,6 +161,7 @@ void Analysis::InitArgs(ArgumentParser* args)
_StopFramesMethodD = args->GetStopFramesMethodD();
_IndividualFDFlags = args->GetIndividualFDFlags();
_plotTimeseriesD=args->GetIsPlotTimeSeriesD();
+ _geoPoly = ReadGeometry(args->GetGeometryFilename(), _areaForMethod_D);
}
if(args->GetIsMethodI()) {
_DoesUseMethodI = true;
@@ -173,6 +174,7 @@ void Analysis::InitArgs(ArgumentParser* args)
_StopFramesMethodI = args->GetStopFramesMethodI();
_IndividualFDFlags = args->GetIndividualFDFlags();
_plotTimeseriesI=args->GetIsPlotTimeSeriesI();
+ _geoPoly = ReadGeometry(args->GetGeometryFilename(), _areaForMethod_I);
}
_deltaF = args->GetDelatT_Vins();
@@ -186,7 +188,6 @@ void Analysis::InitArgs(ArgumentParser* args)
_IgnoreBackwardMovement =args->GetIgnoreBackwardMovement();
_grid_size_X = int(args->GetGridSizeX());
_grid_size_Y = int(args->GetGridSizeY());
- _geoPoly = ReadGeometry(args->GetGeometryFilename(), _areaForMethod_D);
_geometryFileName=args->GetGeometryFilename();
_projectRootDir=args->GetProjectRootDir();
_trajFormat=args->GetFileFormat();
@@ -199,11 +200,11 @@ void Analysis::InitArgs(ArgumentParser* args)
std::map<int, polygon_2d> Analysis::ReadGeometry(const fs::path& geometryFile, const std::vector<MeasurementArea_B*>& areas)
{
-
_building = new Building();
_building->LoadGeometry(geometryFile.string());
// create the polygons
_building->InitGeometry();
+ // _building->AddSurroundingRoom();
double geo_minX = FLT_MAX;
double geo_minY = FLT_MAX;
@@ -221,10 +222,8 @@ std::map<int, polygon_2d> Analysis::ReadGeometry(const fs::path& geometryFile, c
for (auto&& it_sub : it_room.second->GetAllSubRooms())
{
SubRoom* subroom = it_sub.second.get();
-
point_2d point(0,0);
boost::geometry::centroid(area->_poly,point);
-
//check if the area is contained in the obstacle
if(subroom->IsInSubRoom(Point(point.x()/M2CM,point.y()/M2CM)))
{
@@ -245,7 +244,6 @@ std::map<int, polygon_2d> Analysis::ReadGeometry(const fs::path& geometryFile, c
geoPoly[area->_id].inners().resize(k++);
geoPoly[area->_id].inners().back();
model::ring<point_2d>& inner = geoPoly[area->_id].inners().back();
-
for(auto&& tmp_point:obst->GetPolygon())
{
append(inner, make<point_2d>(tmp_point._x*M2CM, tmp_point._y*M2CM));
@@ -267,9 +265,8 @@ std::map<int, polygon_2d> Analysis::ReadGeometry(const fs::path& geometryFile, c
_highVertexY = geo_maxY;
_lowVertexX = geo_minX;
_lowVertexY = geo_minY;
- // using boost::geometry::dsv;
+ using boost::geometry::dsv;
// cout<<"INFO: \tGeometry polygon is:\t" << dsv(geoPoly[1])<<endl;
-
return geoPoly;
}
@@ -488,7 +485,7 @@ int Analysis::RunAnalysis(const fs::path& filename, const fs::path& path)
bool result_I = method_I.Process(data,_scriptsLocation, _areaForMethod_I[i]->_zPos);
if(result_I)
{
- Log->Write("INFO:\tSuccess with Method I uing measurement area id %d!\n",_areaForMethod_I[i]->_id);
+ Log->Write("INFO:\tSuccess with Method I using measurement area id %d!\n",_areaForMethod_I[i]->_id);
std::cout << "INFO:\tSuccess with Method I using measurement area id "<< _areaForMethod_I[i]->_id << "\n";
if(_plotTimeseriesI[i])
{
diff --git a/methods/Method_I.cpp b/methods/Method_I.cpp
index b1b14959c5ace5f02ac189d4dcb3c948fda3ae11..e12bff72412a0f5c7ddddb4769ccd00b4de1b585 100644
--- a/methods/Method_I.cpp
+++ b/methods/Method_I.cpp
@@ -152,665 +152,656 @@ bool Method_I::Process (const PedData& peddata,const fs::path& scriptsLocation,
XInFrame.erase(XInFrame.begin() + i);
YInFrame.erase(YInFrame.begin() + i);
VInFrame.erase(VInFrame.begin() + i);
+ Log->Write("Warning:\t Pedestrian removed");
i--;
}
}
int NumPeds = IdInFrame.size();
//---------------------------------------------------------------------------------------------------------------
- if(NumPeds>3)
+ if(_isOneDimensional)
{
- if(_isOneDimensional)
- {
- CalcVoronoiResults1D(XInFrame, VInFrame, IdInFrame, _areaForMethod_I->_poly,str_frid);
- }
- else
+ CalcVoronoiResults1D(XInFrame, VInFrame, IdInFrame, _areaForMethod_I->_poly,str_frid);
+ }
+ else
+ {
+ if(IsPointsOnOneLine(XInFrame, YInFrame))
{
- if(IsPointsOnOneLine(XInFrame, YInFrame))
+ if(fabs(XInFrame[1]-XInFrame[0])<dmin)
{
- if(fabs(XInFrame[1]-XInFrame[0])<dmin)
- {
- XInFrame[1]+= offset;
- }
- else
- {
- YInFrame[1]+= offset;
- }
+ XInFrame[1]+= offset;
}
- std::vector<std::pair<polygon_2d, int> > polygons_id = GetPolygons(XInFrame, YInFrame, VInFrame, IdInFrame);
- // std::cout << ">> polygons_id " << polygons_id.size() << "\n";
- vector<polygon_2d> polygons;
- for (auto p: polygons_id)
- polygons.push_back(p.first);
+ else
+ {
+ YInFrame[1]+= offset;
+ }
+ }
+ std::vector<std::pair<polygon_2d, int> > polygons_id = GetPolygons(XInFrame, YInFrame, VInFrame, IdInFrame);
+ // std::cout << ">> polygons_id " << polygons_id.size() << "\n";
+ vector<polygon_2d> polygons;
+ for (auto p: polygons_id)
+ polygons.push_back(p.first);
- if(!polygons.empty())
+ if(!polygons.empty())
+ {
+ OutputVoronoiResults(polygons, str_frid, VInFrame); // TODO polygons_id
+ if(_calcIndividualFD)
{
- OutputVoronoiResults(polygons, str_frid, VInFrame); // TODO polygons_id
- if(_calcIndividualFD)
+ if(!_isOneDimensional)
{
- if(!_isOneDimensional)
- {
- GetIndividualFD(polygons,VInFrame, IdInFrame, _areaForMethod_I->_poly, str_frid); // TODO polygons_id
- }
- }
- if(_getProfile)
- { // field analysis
- GetProfiles(str_frid, polygons, VInFrame); // TODO polygons_id
- }
- if(_outputVoronoiCellData)
- { // output the Voronoi polygons of a frame
- OutputVoroGraph(str_frid, polygons_id, NumPeds, XInFrame, YInFrame,VInFrame);
+ GetIndividualFD(polygons,VInFrame, IdInFrame, _areaForMethod_I->_poly, str_frid); // TODO polygons_id
}
}
- else
+ if(_getProfile)
+ { // field analysis
+ GetProfiles(str_frid, polygons, VInFrame); // TODO polygons_id
+ }
+ if(_outputVoronoiCellData)
+ { // output the Voronoi polygons of a frame
+ OutputVoroGraph(str_frid, polygons_id, NumPeds, XInFrame, YInFrame,VInFrame);
+ }
+ }
+ else
+ {
+ for(int i=0;i<(int)IdInFrame.size();i++)
{
- for(int i=0;i<(int)IdInFrame.size();i++)
- {
- std::cout << XInFrame[i]*CMtoM << " " << YInFrame[i]*CMtoM << " " << IdInFrame[i] << "\n";
- }
- Log->Write("WARNING: \tVoronoi Diagrams are not obtained!. Frame: %d (minFrame = %d)\n", frid, minFrame);
+ std::cout << XInFrame[i]*CMtoM << " " << YInFrame[i]*CMtoM << " " << IdInFrame[i] << "\n";
}
+ Log->Write("WARNING: \tVoronoi Diagrams are not obtained!. Frame: %d (minFrame = %d)\n", frid, minFrame);
}
- } // if N >3
- else
+ }
+ fclose(_fVoronoiRhoV);
+ if(_calcIndividualFD)
{
- Log->Write("INFO: \tThe number of the pedestrians is small (%d). Frame = %d (minFrame = %d)\n", NumPeds, frid, minFrame);
+ fclose(_fIndividualFD);
}
+ return return_value;
}
- fclose(_fVoronoiRhoV);
- if(_calcIndividualFD)
- {
- fclose(_fIndividualFD);
- }
- return return_value;
-}
-bool Method_I::OpenFileMethodD()
-{
+ bool Method_I::OpenFileMethodD()
+ {
- std::string voroLocation(VORO_LOCATION);
- fs::path tmp("_id_"+_measureAreaId+".dat");
- tmp = _outputLocation / voroLocation / ("rho_v_Voronoi_" + _trajName.string() + tmp.string());
+ std::string voroLocation(VORO_LOCATION);
+ fs::path tmp("_id_"+_measureAreaId+".dat");
+ tmp = _outputLocation / voroLocation / ("rho_v_Voronoi_" + _trajName.string() + tmp.string());
// _outputLocation.string() + voroLocation+"rho_v_Voronoi_"+_trajName+"_id_"+_measureAreaId+".dat";
- string results_V= tmp.string();
+ string results_V= tmp.string();
- if((_fVoronoiRhoV=Analysis::CreateFile(results_V))==nullptr)
- {
- Log->Write("ERROR: \tcannot open the file to write Voronoi density and velocity\n");
- return false;
- }
- else
- {
- if(_isOneDimensional)
+ if((_fVoronoiRhoV=Analysis::CreateFile(results_V))==nullptr)
{
- fprintf(_fVoronoiRhoV,"#framerate:\t%.2f\n\n#Frame \t Voronoi density(m^(-1))\t Voronoi velocity(m/s)\n",_fps);
+ Log->Write("ERROR: \tcannot open the file to write Voronoi density and velocity\n");
+ return false;
}
else
{
- fprintf(_fVoronoiRhoV,"#framerate:\t%.2f\n\n#Frame \t Voronoi density(m^(-2))\t Voronoi velocity(m/s)\n",_fps);
+ if(_isOneDimensional)
+ {
+ fprintf(_fVoronoiRhoV,"#framerate:\t%.2f\n\n#Frame \t Voronoi density(m^(-1))\t Voronoi velocity(m/s)\n",_fps);
+ }
+ else
+ {
+ fprintf(_fVoronoiRhoV,"#framerate:\t%.2f\n\n#Frame \t Voronoi density(m^(-2))\t Voronoi velocity(m/s)\n",_fps);
+ }
+ return true;
}
- return true;
}
-}
-bool Method_I::OpenFileIndividualFD()
-{
- fs::path trajFileName("_id_"+_measureAreaId+".dat");
- fs::path indFDPath("Fundamental_Diagram");
- indFDPath = _outputLocation / indFDPath / "IndividualFD" / (_trajName.string() + trajFileName.string());
- string Individualfundment=indFDPath.string();
- if((_fIndividualFD=Analysis::CreateFile(Individualfundment))==nullptr)
- {
- Log->Write("ERROR:\tcannot open the file individual\n");
- return false;
- }
- else
+ bool Method_I::OpenFileIndividualFD()
{
- if(_isOneDimensional)
+ fs::path trajFileName("_id_"+_measureAreaId+".dat");
+ fs::path indFDPath("Fundamental_Diagram");
+ indFDPath = _outputLocation / indFDPath / "IndividualFD" / (_trajName.string() + trajFileName.string());
+ string Individualfundment=indFDPath.string();
+ if((_fIndividualFD=Analysis::CreateFile(Individualfundment))==nullptr)
{
- fprintf(_fIndividualFD,"#framerate (fps):\t%.2f\n\n#Frame \t PedId \t Individual density(m^(-1)) \t Individual velocity(m/s) \t Headway(m)\n",_fps);
+ Log->Write("ERROR:\tcannot open the file individual\n");
+ return false;
}
else
{
- fprintf(_fIndividualFD,"#framerate (fps):\t%.2f\n\n#Frame \t PedId \t Individual density(m^(-2)) \t Individual velocity(m/s) \t Voronoi Polygon \t Intersection Polygon\n",_fps);
+ if(_isOneDimensional)
+ {
+ fprintf(_fIndividualFD,"#framerate (fps):\t%.2f\n\n#Frame \t PedId \t Individual density(m^(-1)) \t Individual velocity(m/s) \t Headway(m)\n",_fps);
+ }
+ else
+ {
+ fprintf(_fIndividualFD,"#framerate (fps):\t%.2f\n\n#Frame \t PedId \t Individual density(m^(-2)) \t Individual velocity(m/s) \t Voronoi Polygon \t Intersection Polygon\n",_fps);
+ }
+ return true;
}
- return true;
}
-}
-std::vector<std::pair<polygon_2d, int> > Method_I::GetPolygons(vector<double>& XInFrame, vector<double>& YInFrame, vector<double>& VInFrame, vector<int>& IdInFrame)
-{
- VoronoiDiagram vd;
- //int NrInFrm = ids.size();
- double boundpoint =10*max(max(fabs(_geoMinX),fabs(_geoMinY)), max(fabs(_geoMaxX), fabs(_geoMaxY)));
- std::vector<std::pair<polygon_2d, int> > polygons_id;
- polygons_id = vd.getVoronoiPolygons(XInFrame, YInFrame, VInFrame,IdInFrame, boundpoint);
- // std:: cout << " GetPolygons " << polygons_id.size() << "\n";
-
- polygon_2d poly ;
- if(_cutByCircle)
- {
- polygons_id = vd.cutPolygonsWithCircle(polygons_id, XInFrame, YInFrame, _cutRadius,_circleEdges);
- }
- // std:: cout << " GetPolygons cirlces " << polygons_id.size() << "\n";
- polygons_id = vd.cutPolygonsWithGeometry(polygons_id, _geoPoly, XInFrame, YInFrame);
- // std:: cout << " GetPolygons geometry " << polygons_id.size() << "\n";
- for(auto && p:polygons_id)
+ std::vector<std::pair<polygon_2d, int> > Method_I::GetPolygons(vector<double>& XInFrame, vector<double>& YInFrame, vector<double>& VInFrame, vector<int>& IdInFrame)
{
- poly = p.first;
- ReducePrecision(poly);
- // TODO update polygon_id?
+ VoronoiDiagram vd;
+ //int NrInFrm = ids.size();
+ double boundpoint =10*max(max(fabs(_geoMinX),fabs(_geoMinY)), max(fabs(_geoMaxX), fabs(_geoMaxY)));
+ std::vector<std::pair<polygon_2d, int> > polygons_id;
+ polygons_id = vd.getVoronoiPolygons(XInFrame, YInFrame, VInFrame,IdInFrame, boundpoint);
+ polygon_2d poly ;
+ if(_cutByCircle)
+ {
+ polygons_id = vd.cutPolygonsWithCircle(polygons_id, XInFrame, YInFrame, _cutRadius,_circleEdges);
+ }
+ // std:: cout << " GetPolygons cirlces " << polygons_id.size() << "\n";
+ polygons_id = vd.cutPolygonsWithGeometry(polygons_id, _geoPoly, XInFrame, YInFrame);
+ // std:: cout << " GetPolygons geometry " << polygons_id.size() << "\n";
+ for(auto && p:polygons_id)
+ {
+ poly = p.first;
+ ReducePrecision(poly);
+ // TODO update polygon_id?
+ }
+ // std:: cout << " GetPolygons leave " << polygons_id.size() << "\n";
+ return polygons_id;
}
- // std:: cout << " GetPolygons leave " << polygons_id.size() << "\n";
- return polygons_id;
-}
/**
* Output the Voronoi density and velocity in the corresponding file
*/
-void Method_I::OutputVoronoiResults(const vector<polygon_2d>& polygons, const string& frid, const vector<double>& VInFrame)
-{
- double VoronoiVelocity=1;
- double VoronoiDensity=-1;
- std::tie(VoronoiDensity, VoronoiVelocity) = GetVoronoiDensityVelocity(polygons,VInFrame,_areaForMethod_I->_poly);
- fprintf(_fVoronoiRhoV,"%s\t%.3f\t%.3f\n",frid.c_str(),VoronoiDensity, VoronoiVelocity);
-}
+ void Method_I::OutputVoronoiResults(const vector<polygon_2d>& polygons, const string& frid, const vector<double>& VInFrame)
+ {
+ double VoronoiVelocity=1;
+ double VoronoiDensity=-1;
+ std::tie(VoronoiDensity, VoronoiVelocity) = GetVoronoiDensityVelocity(polygons,VInFrame,_areaForMethod_I->_poly);
+ fprintf(_fVoronoiRhoV,"%s\t%.3f\t%.3f\n",frid.c_str(),VoronoiDensity, VoronoiVelocity);
+ }
/*
* calculate the voronoi density and velocity according to voronoi cell of each pedestrian and their instantaneous velocity "Velocity".
* input: voronoi cell and velocity of each pedestrian and the measurement area
* output: the voronoi density and velocity in the measurement area (tuple)
*/
-std::tuple<double,double> Method_I::GetVoronoiDensityVelocity(const vector<polygon_2d>& polygon, const vector<double>& Velocity, const polygon_2d & measureArea)
-{
- double meanV=0;
- double density=0;
- int temp=0;
- for (auto && polygon_iterator:polygon)
- {
- polygon_list v;
- intersection(measureArea, polygon_iterator, v);
- if(!v.empty())
- {
- meanV+=Velocity[temp]*area(v[0]);
- density+=area(v[0])/area(polygon_iterator);
- if((area(v[0]) - area(polygon_iterator))>J_EPS)
+ std::tuple<double,double> Method_I::GetVoronoiDensityVelocity(const vector<polygon_2d>& polygon, const vector<double>& Velocity, const polygon_2d & measureArea)
+ {
+ double meanV=0;
+ double density=0;
+ int temp=0;
+ for (auto && polygon_iterator:polygon)
+ {
+ polygon_list v;
+ intersection(measureArea, polygon_iterator, v);
+ if(!v.empty())
{
- std::cout<<"----------------------Now calculating density-velocity!!!-----------------\n ";
- std::cout<<"measure area: \t"<<std::setprecision(16)<<dsv(measureArea)<<"\n";
- std::cout<<"Original polygon:\t"<<std::setprecision(16)<<dsv(polygon_iterator)<<"\n";
- std::cout<<"intersected polygon: \t"<<std::setprecision(16)<<dsv(v[0])<<"\n";
- std::cout<<"this is a wrong result in density calculation\t "<<area(v[0])<<'\t'<<area(polygon_iterator)<< " (diff=" << (area(v[0]) - area(polygon_iterator)) << ")" << "\n";
+ meanV+=Velocity[temp]*area(v[0]);
+ density+=area(v[0])/area(polygon_iterator);
+ if((area(v[0]) - area(polygon_iterator))>J_EPS)
+ {
+ std::cout<<"----------------------Now calculating density-velocity!!!-----------------\n ";
+ std::cout<<"measure area: \t"<<std::setprecision(16)<<dsv(measureArea)<<"\n";
+ std::cout<<"Original polygon:\t"<<std::setprecision(16)<<dsv(polygon_iterator)<<"\n";
+ std::cout<<"intersected polygon: \t"<<std::setprecision(16)<<dsv(v[0])<<"\n";
+ std::cout<<"this is a wrong result in density calculation\t "<<area(v[0])<<'\t'<<area(polygon_iterator)<< " (diff=" << (area(v[0]) - area(polygon_iterator)) << ")" << "\n";
+ }
}
+ temp++;
}
- temp++;
+ meanV = meanV/area(measureArea);
+ density = density/(area(measureArea)*CMtoM*CMtoM);
+ return std::make_tuple(density, meanV);
}
- meanV = meanV/area(measureArea);
- density = density/(area(measureArea)*CMtoM*CMtoM);
- return std::make_tuple(density, meanV);
-}
// and velocity is calculated for every frame
-void Method_I::GetProfiles(const string& frameId, const vector<polygon_2d>& polygons, const vector<double>& velocity)
-{
- std::string voroLocation(VORO_LOCATION);
- fs::path tmp("field");
- fs::path vtmp ("velocity");
- fs::path dtmp("density");
- tmp = _outputLocation / voroLocation / tmp;
- vtmp = tmp / vtmp / ("Prf_v_"+_trajName.string()+"_id_"+_measureAreaId+"_"+frameId+".dat");
- dtmp = tmp / dtmp / ("Prf_d_"+_trajName.string()+"_id_"+_measureAreaId+"_"+frameId+".dat");
- //string voronoiLocation=_outputLocation.string() + voroLocation+"field/";
-
- // string Prfvelocity=voronoiLocation+"/velocity/Prf_v_"+_trajName.string()+"_id_"+_measureAreaId+"_"+frameId+".dat";
- // string Prfdensity=voronoiLocation+"/density/Prf_d_"+_trajName.string()+"_id_"+_measureAreaId+"_"+frameId+".dat";
- string Prfvelocity = vtmp.string();
- string Prfdensity = dtmp.string();
-
- FILE *Prf_velocity;
- if((Prf_velocity=Analysis::CreateFile(Prfvelocity))==nullptr) {
- Log->Write("cannot open the file <%s> to write the field data\n",Prfvelocity.c_str());
- exit(EXIT_FAILURE);
- }
- FILE *Prf_density;
- if((Prf_density=Analysis::CreateFile(Prfdensity))==nullptr) {
- Log->Write("cannot open the file to write the field density\n");
- exit(EXIT_FAILURE);
- }
-
- int NRow = (int)ceil((_geoMaxY-_geoMinY)/_grid_size_Y); // the number of rows that the geometry will be discretized for field analysis
- int NColumn = (int)ceil((_geoMaxX-_geoMinX)/_grid_size_X); //the number of columns that the geometry will be discretized for field analysis
- for(int row_i=0; row_i<NRow; row_i++) { //
- for(int colum_j=0; colum_j<NColumn; colum_j++) {
- polygon_2d measurezoneXY;
- {
- const double coor[][2] = {
- {_geoMinX+colum_j*_grid_size_X,_geoMaxY-row_i*_grid_size_Y}, {_geoMinX+colum_j*_grid_size_X+_grid_size_X,_geoMaxY-row_i*_grid_size_Y}, {_geoMinX+colum_j*_grid_size_X+_grid_size_X, _geoMaxY-row_i*_grid_size_Y-_grid_size_Y},
- {_geoMinX+colum_j*_grid_size_X, _geoMaxY-row_i*_grid_size_Y-_grid_size_Y},
- {_geoMinX+colum_j*_grid_size_X,_geoMaxY-row_i*_grid_size_Y} // closing point is opening point
- };
- assign_points(measurezoneXY, coor);
- }
- correct(measurezoneXY); // Polygons should be closed, and directed clockwise. If you're not sure if that is the case, call this function
-
- double densityXY;
- double velocityXY;
- std::tie(densityXY, velocityXY) = GetVoronoiDensityVelocity(polygons,velocity,measurezoneXY);
- fprintf(Prf_density,"%.3f\t",densityXY);
- fprintf(Prf_velocity,"%.3f\t",velocityXY);
+ void Method_I::GetProfiles(const string& frameId, const vector<polygon_2d>& polygons, const vector<double>& velocity)
+ {
+ std::string voroLocation(VORO_LOCATION);
+ fs::path tmp("field");
+ fs::path vtmp ("velocity");
+ fs::path dtmp("density");
+ tmp = _outputLocation / voroLocation / tmp;
+ vtmp = tmp / vtmp / ("Prf_v_"+_trajName.string()+"_id_"+_measureAreaId+"_"+frameId+".dat");
+ dtmp = tmp / dtmp / ("Prf_d_"+_trajName.string()+"_id_"+_measureAreaId+"_"+frameId+".dat");
+ //string voronoiLocation=_outputLocation.string() + voroLocation+"field/";
+
+ // string Prfvelocity=voronoiLocation+"/velocity/Prf_v_"+_trajName.string()+"_id_"+_measureAreaId+"_"+frameId+".dat";
+ // string Prfdensity=voronoiLocation+"/density/Prf_d_"+_trajName.string()+"_id_"+_measureAreaId+"_"+frameId+".dat";
+ string Prfvelocity = vtmp.string();
+ string Prfdensity = dtmp.string();
+
+ FILE *Prf_velocity;
+ if((Prf_velocity=Analysis::CreateFile(Prfvelocity))==nullptr) {
+ Log->Write("cannot open the file <%s> to write the field data\n",Prfvelocity.c_str());
+ exit(EXIT_FAILURE);
+ }
+ FILE *Prf_density;
+ if((Prf_density=Analysis::CreateFile(Prfdensity))==nullptr) {
+ Log->Write("cannot open the file to write the field density\n");
+ exit(EXIT_FAILURE);
}
- fprintf(Prf_density,"\n");
- fprintf(Prf_velocity,"\n");
- }
- fclose(Prf_velocity);
- fclose(Prf_density);
-}
-void Method_I::OutputVoroGraph(const string & frameId, std::vector<std::pair<polygon_2d, int> >& polygons_id, int numPedsInFrame, vector<double>& XInFrame, vector<double>& YInFrame,const vector<double>& VInFrame)
-{
- //string voronoiLocation=_projectRootDir+"./Output/Fundamental_Diagram/Classical_Voronoi/VoronoiCell/id_"+_measureAreaId;
+ int NRow = (int)ceil((_geoMaxY-_geoMinY)/_grid_size_Y); // the number of rows that the geometry will be discretized for field analysis
+ int NColumn = (int)ceil((_geoMaxX-_geoMinX)/_grid_size_X); //the number of columns that the geometry will be discretized for field analysis
+ for(int row_i=0; row_i<NRow; row_i++) { //
+ for(int colum_j=0; colum_j<NColumn; colum_j++) {
+ polygon_2d measurezoneXY;
+ {
+ const double coor[][2] = {
+ {_geoMinX+colum_j*_grid_size_X,_geoMaxY-row_i*_grid_size_Y}, {_geoMinX+colum_j*_grid_size_X+_grid_size_X,_geoMaxY-row_i*_grid_size_Y}, {_geoMinX+colum_j*_grid_size_X+_grid_size_X, _geoMaxY-row_i*_grid_size_Y-_grid_size_Y},
+ {_geoMinX+colum_j*_grid_size_X, _geoMaxY-row_i*_grid_size_Y-_grid_size_Y},
+ {_geoMinX+colum_j*_grid_size_X,_geoMaxY-row_i*_grid_size_Y} // closing point is opening point
+ };
+ assign_points(measurezoneXY, coor);
+ }
+ correct(measurezoneXY); // Polygons should be closed, and directed clockwise. If you're not sure if that is the case, call this function
- fs::path voroLocPath(_outputLocation);
- fs::path voro_location_path (VORO_LOCATION); // TODO: convert
- // this MACRO to
- // path. Maybe
- // remove the MACRO?
- voroLocPath = voroLocPath / voro_location_path / "VoronoiCell";
- polygon_2d poly;
- if(!fs::exists(voroLocPath))
- {
- if(!fs::create_directories(voroLocPath))
- {
- Log->Write("ERROR:\tcan not create directory <%s>", voroLocPath.string().c_str());
- std::cout << "can not create directory "<< voroLocPath.string().c_str() << "\n";
- exit(EXIT_FAILURE);
- }
- else
- std::cout << "create directory "<< voroLocPath.string().c_str() << "\n";
+ double densityXY;
+ double velocityXY;
+ std::tie(densityXY, velocityXY) = GetVoronoiDensityVelocity(polygons,velocity,measurezoneXY);
+ fprintf(Prf_density,"%.3f\t",densityXY);
+ fprintf(Prf_velocity,"%.3f\t",velocityXY);
+ }
+ fprintf(Prf_density,"\n");
+ fprintf(Prf_velocity,"\n");
+ }
+ fclose(Prf_velocity);
+ fclose(Prf_density);
}
- fs::path polygonPath=voroLocPath / "polygon";
- if(!fs::exists(polygonPath))
+ void Method_I::OutputVoroGraph(const string & frameId, std::vector<std::pair<polygon_2d, int> >& polygons_id, int numPedsInFrame, vector<double>& XInFrame, vector<double>& YInFrame,const vector<double>& VInFrame)
{
- if(!fs::create_directory(polygonPath))
+ //string voronoiLocation=_projectRootDir+"./Output/Fundamental_Diagram/Classical_Voronoi/VoronoiCell/id_"+_measureAreaId;
+
+ fs::path voroLocPath(_outputLocation);
+ fs::path voro_location_path (VORO_LOCATION); // TODO: convert
+ // this MACRO to
+ // path. Maybe
+ // remove the MACRO?
+ voroLocPath = voroLocPath / voro_location_path / "VoronoiCell";
+ polygon_2d poly;
+ if(!fs::exists(voroLocPath))
{
- Log->Write("ERROR:\tcan not create directory <%s>", polygonPath.string().c_str());
- exit(EXIT_FAILURE);
+ if(!fs::create_directories(voroLocPath))
+ {
+ Log->Write("ERROR:\tcan not create directory <%s>", voroLocPath.string().c_str());
+ std::cout << "can not create directory "<< voroLocPath.string().c_str() << "\n";
+ exit(EXIT_FAILURE);
+ }
+ else
+ std::cout << "create directory "<< voroLocPath.string().c_str() << "\n";
}
- }
- fs::path trajFileName(_trajName.string()+"_id_"+_measureAreaId+"_"+frameId+".dat");
- fs::path p = polygonPath / trajFileName;
- string polygon = p.string();
- ofstream polys (polygon.c_str());
- if(polys.is_open())
- {
- //for(vector<polygon_2d> polygon_iterator=polygons.begin(); polygon_iterator!=polygons.end(); polygon_iterator++)
- for(auto && p:polygons_id)
+ fs::path polygonPath=voroLocPath / "polygon";
+ if(!fs::exists(polygonPath))
{
- poly = p.first;
- for(auto&& point:poly.outer())
+ if(!fs::create_directory(polygonPath))
{
- point.x(point.x()*CMtoM);
- point.y(point.y()*CMtoM);
+ Log->Write("ERROR:\tcan not create directory <%s>", polygonPath.string().c_str());
+ exit(EXIT_FAILURE);
}
- for(auto&& innerpoly:poly.inners())
+ }
+ fs::path trajFileName(_trajName.string()+"_id_"+_measureAreaId+"_"+frameId+".dat");
+ fs::path p = polygonPath / trajFileName;
+ string polygon = p.string();
+ ofstream polys (polygon.c_str());
+
+ if(polys.is_open())
+ {
+ //for(vector<polygon_2d> polygon_iterator=polygons.begin(); polygon_iterator!=polygons.end(); polygon_iterator++)
+ for(auto && p:polygons_id)
{
- for(auto&& point:innerpoly)
+ poly = p.first;
+ for(auto&& point:poly.outer())
{
point.x(point.x()*CMtoM);
point.y(point.y()*CMtoM);
}
+ for(auto&& innerpoly:poly.inners())
+ {
+ for(auto&& point:innerpoly)
+ {
+ point.x(point.x()*CMtoM);
+ point.y(point.y()*CMtoM);
+ }
+ }
+ polys << p.second << " | " << dsv(poly) << endl;
+ //polys <<dsv(poly)<< endl;
}
- polys << p.second << " | " << dsv(poly) << endl;
- //polys <<dsv(poly)<< endl;
}
- }
- else
- {
- Log->Write("ERROR:\tcannot create the file <%s>",polygon.c_str());
- exit(EXIT_FAILURE);
- }
- fs::path speedPath=voroLocPath / "speed";
- if(!fs::exists(speedPath))
- if(!fs::create_directory(speedPath))
+ else
{
- Log->Write("ERROR:\tcan not create directory <%s>", speedPath.string().c_str());
+ Log->Write("ERROR:\tcannot create the file <%s>",polygon.c_str());
exit(EXIT_FAILURE);
}
- fs::path pv = speedPath /trajFileName;
- string v_individual= pv.string();
- ofstream velo (v_individual.c_str());
- if(velo.is_open())
- {
- for(int pts=0; pts<numPedsInFrame; pts++)
+ fs::path speedPath=voroLocPath / "speed";
+ if(!fs::exists(speedPath))
+ if(!fs::create_directory(speedPath))
+ {
+ Log->Write("ERROR:\tcan not create directory <%s>", speedPath.string().c_str());
+ exit(EXIT_FAILURE);
+ }
+ fs::path pv = speedPath /trajFileName;
+ string v_individual= pv.string();
+ ofstream velo (v_individual.c_str());
+ if(velo.is_open())
{
- velo << fabs(VInFrame[pts]) << endl;
+ for(int pts=0; pts<numPedsInFrame; pts++)
+ {
+ velo << fabs(VInFrame[pts]) << endl;
+ }
}
+ else
+ {
+ Log->Write("ERROR:\tcannot create the file <%s>",pv.string().c_str());
+ exit(EXIT_FAILURE);
+ }
+
+ /*string point=voronoiLocation+"/points"+_trajName+"_id_"+_measureAreaId+"_"+frameId+".dat";
+ ofstream points (point.c_str());
+ if( points.is_open())
+ {
+ for(int pts=0; pts<numPedsInFrame; pts++)
+ {
+ points << XInFrame[pts]*CMtoM << "\t" << YInFrame[pts]*CMtoM << endl;
+ }
+ }
+ else
+ {
+ Log->Write("ERROR:\tcannot create the file <%s>",point.c_str());
+ exit(EXIT_FAILURE);
+ }*/
+
+ if(_plotVoronoiCellData)
+ {
+ //@todo: use fs::path
+ string parameters_rho=" " + _scriptsLocation.string()+"/_Plot_cell_rho.py -f \""+ voroLocPath.string() + "\" -n "+ _trajName.string()+"_id_"+_measureAreaId+"_"+frameId+
+ " -g "+_geometryFileName.string()+" -p "+_trajectoryPath.string();
+ string parameters_v=" " + _scriptsLocation.string()+"/_Plot_cell_v.py -f \""+ voroLocPath.string() + "\" -n "+ _trajName.string() + "_id_"+_measureAreaId+"_"+frameId+
+ " -g "+_geometryFileName.string()+" -p "+_trajectoryPath.string();
+
+ if(_plotVoronoiIndex)
+ parameters_rho += " -i";
+
+ Log->Write("INFO:\t%s",parameters_rho.c_str());
+ Log->Write("INFO:\tPlotting Voronoi Cell at the frame <%s>",frameId.c_str());
+ parameters_rho = PYTHON + parameters_rho;
+ parameters_v = PYTHON + parameters_v;
+ system(parameters_rho.c_str());
+ system(parameters_v.c_str());
+ }
+ //points.close();
+ polys.close();
+ velo.close();
}
- else
+
+/*std::string polygon_to_string2(const polygon_2d & polygon)
+ {
+ string polygon_str = "((";
+ for(auto point: boost::geometry::exterior_ring(polygon) )
+ {
+ double x = boost::geometry::get<0>(point);
+ double y = boost::geometry::get<1>(point);
+ polygon_str.append("(");
+ polygon_str.append(std::to_string(x));
+ polygon_str.append(", ");
+ polygon_str.append(std::to_string(y));
+ polygon_str.append("), ");
+ }
+ for(auto ring: boost::geometry::interior_rings(polygon) )
+ {
+ for(auto point: ring )
+ {
+ double x = boost::geometry::get<0>(point);
+ double y = boost::geometry::get<1>(point);
+ polygon_str.append("(");
+ polygon_str.append(std::to_string(x));
+ polygon_str.append(", ");
+ polygon_str.append(std::to_string(y));
+ polygon_str.append("), ");
+ }
+ }
+ polygon_str.pop_back(); polygon_str.pop_back(); //remove last komma
+ polygon_str.append("))");
+ return polygon_str;
+ }*/
+
+ void Method_I::GetIndividualFD(const vector<polygon_2d>& polygon, const vector<double>& Velocity, const vector<int>& Id, const polygon_2d& measureArea, const string& frid)
{
- Log->Write("ERROR:\tcannot create the file <%s>",pv.string().c_str());
- exit(EXIT_FAILURE);
+ double uniquedensity=0;
+ double uniquevelocity=0;
+ int uniqueId=0;
+ int temp=0;
+ for (const auto & polygon_iterator:polygon)
+ {
+ polygon_list v;
+ intersection(measureArea, polygon_iterator, v);
+ if(!v.empty()) {
+
+ string polygon_str = polygon_to_string(polygon_iterator);
+ // string measureArea_str =
+ // polygon_to_string(measureArea); // maybe used for debugging
+ string v_str = polygon_to_string(v[0]);
+
+ uniquedensity=1.0/(area(polygon_iterator)*CMtoM*CMtoM);
+ uniquevelocity=Velocity[temp];
+ uniqueId=Id[temp];
+ fprintf(_fIndividualFD,"%s\t%d\t%.3f\t%.3f\t%s\t%s\n",
+ frid.c_str(),
+ uniqueId,
+ uniquedensity,
+ uniquevelocity,
+ polygon_str.c_str(),
+ v_str.c_str());
+ }
+ temp++;
+ }
}
- /*string point=voronoiLocation+"/points"+_trajName+"_id_"+_measureAreaId+"_"+frameId+".dat";
- ofstream points (point.c_str());
- if( points.is_open())
- {
- for(int pts=0; pts<numPedsInFrame; pts++)
- {
- points << XInFrame[pts]*CMtoM << "\t" << YInFrame[pts]*CMtoM << endl;
- }
- }
- else
- {
- Log->Write("ERROR:\tcannot create the file <%s>",point.c_str());
- exit(EXIT_FAILURE);
- }*/
-
- if(_plotVoronoiCellData)
+ void Method_I::SetCalculateIndividualFD(bool individualFD)
{
- //@todo: use fs::path
- string parameters_rho=" " + _scriptsLocation.string()+"/_Plot_cell_rho.py -f \""+ voroLocPath.string() + "\" -n "+ _trajName.string()+"_id_"+_measureAreaId+"_"+frameId+
- " -g "+_geometryFileName.string()+" -p "+_trajectoryPath.string();
- string parameters_v=" " + _scriptsLocation.string()+"/_Plot_cell_v.py -f \""+ voroLocPath.string() + "\" -n "+ _trajName.string() + "_id_"+_measureAreaId+"_"+frameId+
- " -g "+_geometryFileName.string()+" -p "+_trajectoryPath.string();
-
- if(_plotVoronoiIndex)
- parameters_rho += " -i";
-
- Log->Write("INFO:\t%s",parameters_rho.c_str());
- Log->Write("INFO:\tPlotting Voronoi Cell at the frame <%s>",frameId.c_str());
- parameters_rho = PYTHON + parameters_rho;
- parameters_v = PYTHON + parameters_v;
- system(parameters_rho.c_str());
- system(parameters_v.c_str());
+ _calcIndividualFD = individualFD;
}
- //points.close();
- polys.close();
- velo.close();
-}
-/*std::string polygon_to_string2(const polygon_2d & polygon)
-{
- string polygon_str = "((";
- for(auto point: boost::geometry::exterior_ring(polygon) )
- {
- double x = boost::geometry::get<0>(point);
- double y = boost::geometry::get<1>(point);
- polygon_str.append("(");
- polygon_str.append(std::to_string(x));
- polygon_str.append(", ");
- polygon_str.append(std::to_string(y));
- polygon_str.append("), ");
- }
- for(auto ring: boost::geometry::interior_rings(polygon) )
- {
- for(auto point: ring )
- {
- double x = boost::geometry::get<0>(point);
- double y = boost::geometry::get<1>(point);
- polygon_str.append("(");
- polygon_str.append(std::to_string(x));
- polygon_str.append(", ");
- polygon_str.append(std::to_string(y));
- polygon_str.append("), ");
- }
- }
- polygon_str.pop_back(); polygon_str.pop_back(); //remove last komma
- polygon_str.append("))");
- return polygon_str;
-}*/
-
-void Method_I::GetIndividualFD(const vector<polygon_2d>& polygon, const vector<double>& Velocity, const vector<int>& Id, const polygon_2d& measureArea, const string& frid)
-{
- double uniquedensity=0;
- double uniquevelocity=0;
- int uniqueId=0;
- int temp=0;
- for (const auto & polygon_iterator:polygon)
- {
- polygon_list v;
- intersection(measureArea, polygon_iterator, v);
- if(!v.empty()) {
-
- string polygon_str = polygon_to_string(polygon_iterator);
- // string measureArea_str =
- // polygon_to_string(measureArea); // maybe used for debugging
- string v_str = polygon_to_string(v[0]);
-
- uniquedensity=1.0/(area(polygon_iterator)*CMtoM*CMtoM);
- uniquevelocity=Velocity[temp];
- uniqueId=Id[temp];
- fprintf(_fIndividualFD,"%s\t%d\t%.3f\t%.3f\t%s\t%s\n",
- frid.c_str(),
- uniqueId,
- uniquedensity,
- uniquevelocity,
- polygon_str.c_str(),
- v_str.c_str());
- }
- temp++;
+ void Method_I::SetStartFrame(int startFrame)
+ {
+ _startFrame=startFrame;
}
-}
-void Method_I::SetCalculateIndividualFD(bool individualFD)
-{
- _calcIndividualFD = individualFD;
-}
-
-void Method_I::SetStartFrame(int startFrame)
-{
- _startFrame=startFrame;
-}
+ void Method_I::SetStopFrame(int stopFrame)
+ {
+ _stopFrame=stopFrame;
+ }
-void Method_I::SetStopFrame(int stopFrame)
-{
- _stopFrame=stopFrame;
-}
+ void Method_I::Setcutbycircle(double radius,int edges)
+ {
+ _cutByCircle=true;
+ _cutRadius = radius;
+ _circleEdges = edges;
+ }
-void Method_I::Setcutbycircle(double radius,int edges)
-{
- _cutByCircle=true;
- _cutRadius = radius;
- _circleEdges = edges;
-}
+ void Method_I::SetGeometryPolygon(polygon_2d geometryPolygon)
+ {
+ _geoPoly = geometryPolygon;
+ }
-void Method_I::SetGeometryPolygon(polygon_2d geometryPolygon)
-{
- _geoPoly = geometryPolygon;
-}
+ void Method_I::SetGeometryBoundaries(double minX, double minY, double maxX, double maxY)
+ {
+ _geoMinX = minX;
+ _geoMinY = minY;
+ _geoMaxX = maxX;
+ _geoMaxY = maxY;
+ }
-void Method_I::SetGeometryBoundaries(double minX, double minY, double maxX, double maxY)
-{
- _geoMinX = minX;
- _geoMinY = minY;
- _geoMaxX = maxX;
- _geoMaxY = maxY;
-}
+ void Method_I::SetGeometryFileName(const fs::path& geometryFile)
+ {
+ _geometryFileName=geometryFile;
+ }
-void Method_I::SetGeometryFileName(const fs::path& geometryFile)
-{
- _geometryFileName=geometryFile;
-}
+ void Method_I::SetTrajectoriesLocation(const fs::path& trajectoryPath)
+ {
+ _trajectoryPath=trajectoryPath;
+ }
-void Method_I::SetTrajectoriesLocation(const fs::path& trajectoryPath)
-{
- _trajectoryPath=trajectoryPath;
-}
+ void Method_I::SetGridSize(double x, double y)
+ {
+ _grid_size_X = x;
+ _grid_size_Y = y;
+ }
-void Method_I::SetGridSize(double x, double y)
-{
- _grid_size_X = x;
- _grid_size_Y = y;
-}
+ void Method_I::SetCalculateProfiles(bool calcProfile)
+ {
+ _getProfile =calcProfile;
+ }
-void Method_I::SetCalculateProfiles(bool calcProfile)
-{
- _getProfile =calcProfile;
-}
+ void Method_I::SetOutputVoronoiCellData(bool outputCellData)
+ {
+ _outputVoronoiCellData = outputCellData;
+ }
-void Method_I::SetOutputVoronoiCellData(bool outputCellData)
-{
- _outputVoronoiCellData = outputCellData;
-}
+ void Method_I::SetPlotVoronoiGraph(bool plotVoronoiGraph)
+ {
+ _plotVoronoiCellData = plotVoronoiGraph;
+ }
+ void Method_I::SetPlotVoronoiIndex(bool plotVoronoiIndex)
+ {
+ _plotVoronoiIndex = plotVoronoiIndex;
+ }
-void Method_I::SetPlotVoronoiGraph(bool plotVoronoiGraph)
-{
- _plotVoronoiCellData = plotVoronoiGraph;
-}
-void Method_I::SetPlotVoronoiIndex(bool plotVoronoiIndex)
-{
- _plotVoronoiIndex = plotVoronoiIndex;
-}
+ void Method_I::SetMeasurementArea (MeasurementArea_B* area)
+ {
+ _areaForMethod_I = area;
+ }
-void Method_I::SetMeasurementArea (MeasurementArea_B* area)
-{
- _areaForMethod_I = area;
-}
+ void Method_I::SetDimensional (bool dimension)
+ {
+ _isOneDimensional = dimension;
+ }
-void Method_I::SetDimensional (bool dimension)
-{
- _isOneDimensional = dimension;
-}
+ void Method_I::ReducePrecision(polygon_2d& polygon)
+ {
+ for(auto&& point:polygon.outer())
+ {
+ point.x(round(point.x() * 100000000000.0) / 100000000000.0);
+ point.y(round(point.y() * 100000000000.0) / 100000000000.0);
+ }
+ }
-void Method_I::ReducePrecision(polygon_2d& polygon)
-{
- for(auto&& point:polygon.outer())
+ bool Method_I::IsPedInGeometry(int frames, int peds, double **Xcor, double **Ycor, int *firstFrame, int *lastFrame)
{
- point.x(round(point.x() * 100000000000.0) / 100000000000.0);
- point.y(round(point.y() * 100000000000.0) / 100000000000.0);
+ for(int i=0; i<peds; i++)
+ for(int j =0; j<frames; j++)
+ {
+ if (j>firstFrame[i] && j< lastFrame[i] && (false==within(point_2d(round(Xcor[i][j]), round(Ycor[i][j])), _geoPoly)))
+ {
+ Log->Write("Error:\tPedestrian at the position <x=%.4f, y=%.4f> is outside geometry. Please check the geometry or trajectory file!", Xcor[i][j]*CMtoM, Ycor[i][j]*CMtoM );
+ return false;
+ }
+ }
+ return true;
}
-}
-bool Method_I::IsPedInGeometry(int frames, int peds, double **Xcor, double **Ycor, int *firstFrame, int *lastFrame)
-{
- for(int i=0; i<peds; i++)
- for(int j =0; j<frames; j++)
+ void Method_I::CalcVoronoiResults1D(vector<double>& XInFrame, vector<double>& VInFrame, vector<int>& IdInFrame, const polygon_2d & measureArea,const string& frid)
+ {
+ vector<double> measurearea_x;
+ for(unsigned int i=0;i<measureArea.outer().size();i++)
{
- if (j>firstFrame[i] && j< lastFrame[i] && (false==within(point_2d(round(Xcor[i][j]), round(Ycor[i][j])), _geoPoly)))
+ measurearea_x.push_back(measureArea.outer()[i].x());
+ }
+ double left_boundary=*min_element(measurearea_x.begin(),measurearea_x.end());
+ double right_boundary=*max_element(measurearea_x.begin(),measurearea_x.end());
+
+ vector<double> voronoi_distance;
+ vector<double> Xtemp=XInFrame;
+ vector<double> dist;
+ vector<double> XLeftNeighbor;
+ vector<double> XLeftTemp;
+ vector<double> XRightNeighbor;
+ vector<double> XRightTemp;
+ sort(Xtemp.begin(),Xtemp.end());
+ dist.push_back(Xtemp[1]-Xtemp[0]);
+ XLeftTemp.push_back(2*Xtemp[0]-Xtemp[1]);
+ XRightTemp.push_back(Xtemp[1]);
+ for(unsigned int i=1;i<Xtemp.size()-1;i++)
+ {
+ dist.push_back((Xtemp[i+1]-Xtemp[i-1])/2.0);
+ XLeftTemp.push_back(Xtemp[i-1]);
+ XRightTemp.push_back(Xtemp[i+1]);
+ }
+ dist.push_back(Xtemp[Xtemp.size()-1]-Xtemp[Xtemp.size()-2]);
+ XLeftTemp.push_back(Xtemp[Xtemp.size()-2]);
+ XRightTemp.push_back(2*Xtemp[Xtemp.size()-1]-Xtemp[Xtemp.size()-2]);
+ for(unsigned int i=0;i<XInFrame.size();i++)
+ {
+ for(unsigned int j=0;j<Xtemp.size();j++)
{
- Log->Write("Error:\tPedestrian at the position <x=%.4f, y=%.4f> is outside geometry. Please check the geometry or trajectory file!", Xcor[i][j]*CMtoM, Ycor[i][j]*CMtoM );
- return false;
+ if(fabs(XInFrame[i]-Xtemp[j])<1.0e-5)
+ {
+ voronoi_distance.push_back(dist[j]);
+ XLeftNeighbor.push_back(XLeftTemp[j]);
+ XRightNeighbor.push_back(XRightTemp[j]);
+ break;
+ }
}
}
- return true;
-}
-void Method_I::CalcVoronoiResults1D(vector<double>& XInFrame, vector<double>& VInFrame, vector<int>& IdInFrame, const polygon_2d & measureArea,const string& frid)
-{
- vector<double> measurearea_x;
- for(unsigned int i=0;i<measureArea.outer().size();i++)
- {
- measurearea_x.push_back(measureArea.outer()[i].x());
- }
- double left_boundary=*min_element(measurearea_x.begin(),measurearea_x.end());
- double right_boundary=*max_element(measurearea_x.begin(),measurearea_x.end());
-
- vector<double> voronoi_distance;
- vector<double> Xtemp=XInFrame;
- vector<double> dist;
- vector<double> XLeftNeighbor;
- vector<double> XLeftTemp;
- vector<double> XRightNeighbor;
- vector<double> XRightTemp;
- sort(Xtemp.begin(),Xtemp.end());
- dist.push_back(Xtemp[1]-Xtemp[0]);
- XLeftTemp.push_back(2*Xtemp[0]-Xtemp[1]);
- XRightTemp.push_back(Xtemp[1]);
- for(unsigned int i=1;i<Xtemp.size()-1;i++)
- {
- dist.push_back((Xtemp[i+1]-Xtemp[i-1])/2.0);
- XLeftTemp.push_back(Xtemp[i-1]);
- XRightTemp.push_back(Xtemp[i+1]);
- }
- dist.push_back(Xtemp[Xtemp.size()-1]-Xtemp[Xtemp.size()-2]);
- XLeftTemp.push_back(Xtemp[Xtemp.size()-2]);
- XRightTemp.push_back(2*Xtemp[Xtemp.size()-1]-Xtemp[Xtemp.size()-2]);
- for(unsigned int i=0;i<XInFrame.size();i++)
- {
- for(unsigned int j=0;j<Xtemp.size();j++)
- {
- if(fabs(XInFrame[i]-Xtemp[j])<1.0e-5)
+ double VoronoiDensity=0;
+ double VoronoiVelocity=0;
+ for(unsigned int i=0; i<XInFrame.size(); i++)
+ {
+ double ratio=getOverlapRatio((XInFrame[i]+XLeftNeighbor[i])/2.0, (XRightNeighbor[i]+XInFrame[i])/2.0,left_boundary,right_boundary);
+ VoronoiDensity+=ratio;
+ VoronoiVelocity+=(VInFrame[i]*voronoi_distance[i]*ratio*CMtoM);
+ if(_calcIndividualFD)
{
- voronoi_distance.push_back(dist[j]);
- XLeftNeighbor.push_back(XLeftTemp[j]);
- XRightNeighbor.push_back(XRightTemp[j]);
- break;
+ double headway=(XRightNeighbor[i] - XInFrame[i])*CMtoM;
+ double individualDensity = 2.0/((XRightNeighbor[i] - XLeftNeighbor[i])*CMtoM);
+ fprintf(_fIndividualFD,"%s\t%d\t%.3f\t%.3f\t%.3f\n",frid.c_str(), IdInFrame[i], individualDensity,VInFrame[i], headway);
}
}
+ VoronoiDensity/=((right_boundary-left_boundary)*CMtoM);
+ VoronoiVelocity/=((right_boundary-left_boundary)*CMtoM);
+ fprintf(_fVoronoiRhoV,"%s\t%.3f\t%.3f\n",frid.c_str(),VoronoiDensity, VoronoiVelocity);
+
}
- double VoronoiDensity=0;
- double VoronoiVelocity=0;
- for(unsigned int i=0; i<XInFrame.size(); i++)
+ double Method_I::getOverlapRatio(const double& left, const double& right, const double& measurearea_left, const double& measurearea_right)
{
- double ratio=getOverlapRatio((XInFrame[i]+XLeftNeighbor[i])/2.0, (XRightNeighbor[i]+XInFrame[i])/2.0,left_boundary,right_boundary);
- VoronoiDensity+=ratio;
- VoronoiVelocity+=(VInFrame[i]*voronoi_distance[i]*ratio*CMtoM);
- if(_calcIndividualFD)
+ double OverlapRatio=0;
+ double PersonalSpace=right-left;
+ if(left > measurearea_left && right < measurearea_right) //case1
{
- double headway=(XRightNeighbor[i] - XInFrame[i])*CMtoM;
- double individualDensity = 2.0/((XRightNeighbor[i] - XLeftNeighbor[i])*CMtoM);
- fprintf(_fIndividualFD,"%s\t%d\t%.3f\t%.3f\t%.3f\n",frid.c_str(), IdInFrame[i], individualDensity,VInFrame[i], headway);
+ OverlapRatio=1;
}
+ else if(right > measurearea_left && right < measurearea_right && left < measurearea_left)
+ {
+ OverlapRatio=(right-measurearea_left)/PersonalSpace;
+ }
+ else if(left < measurearea_left && right > measurearea_right)
+ {
+ OverlapRatio=(measurearea_right - measurearea_left)/PersonalSpace;
+ }
+ else if(left > measurearea_left && left < measurearea_right && right > measurearea_right)
+ {
+ OverlapRatio=(measurearea_right-left)/PersonalSpace;
+ }
+ return OverlapRatio;
}
- VoronoiDensity/=((right_boundary-left_boundary)*CMtoM);
- VoronoiVelocity/=((right_boundary-left_boundary)*CMtoM);
- fprintf(_fVoronoiRhoV,"%s\t%.3f\t%.3f\n",frid.c_str(),VoronoiDensity, VoronoiVelocity);
-
-}
-
-double Method_I::getOverlapRatio(const double& left, const double& right, const double& measurearea_left, const double& measurearea_right)
-{
- double OverlapRatio=0;
- double PersonalSpace=right-left;
- if(left > measurearea_left && right < measurearea_right) //case1
- {
- OverlapRatio=1;
- }
- else if(right > measurearea_left && right < measurearea_right && left < measurearea_left)
- {
- OverlapRatio=(right-measurearea_left)/PersonalSpace;
- }
- else if(left < measurearea_left && right > measurearea_right)
- {
- OverlapRatio=(measurearea_right - measurearea_left)/PersonalSpace;
- }
- else if(left > measurearea_left && left < measurearea_right && right > measurearea_right)
- {
- OverlapRatio=(measurearea_right-left)/PersonalSpace;
- }
- return OverlapRatio;
-}
-bool Method_I::IsPointsOnOneLine(vector<double>& XInFrame, vector<double>& YInFrame)
-{
- double deltaX=XInFrame[1] - XInFrame[0];
- bool isOnOneLine=true;
- if(fabs(deltaX)<dmin)
+ bool Method_I::IsPointsOnOneLine(vector<double>& XInFrame, vector<double>& YInFrame)
{
- for(unsigned int i=2; i<XInFrame.size();i++)
+ double deltaX=XInFrame[1] - XInFrame[0];
+ bool isOnOneLine=true;
+ if(fabs(deltaX)<dmin)
{
- if(fabs(XInFrame[i] - XInFrame[0])> dmin)
+ for(unsigned int i=2; i<XInFrame.size();i++)
{
- isOnOneLine=false;
- break;
+ if(fabs(XInFrame[i] - XInFrame[0])> dmin)
+ {
+ isOnOneLine=false;
+ break;
+ }
}
}
- }
- else
- {
- double slope=(YInFrame[1] - YInFrame[0])/deltaX;
- double intercept=YInFrame[0] - slope*XInFrame[0];
- for(unsigned int i=2; i<XInFrame.size();i++)
+ else
{
- double dist=fabs(slope*XInFrame[i] - YInFrame[i] + intercept)/sqrt(slope*slope +1);
- if(dist > dmin)
+ double slope=(YInFrame[1] - YInFrame[0])/deltaX;
+ double intercept=YInFrame[0] - slope*XInFrame[0];
+ for(unsigned int i=2; i<XInFrame.size();i++)
{
- isOnOneLine=false;
- break;
+ double dist=fabs(slope*XInFrame[i] - YInFrame[i] + intercept)/sqrt(slope*slope +1);
+ if(dist > dmin)
+ {
+ isOnOneLine=false;
+ break;
+ }
}
}
+ return isOnOneLine;
}
- return isOnOneLine;
-}
diff --git a/methods/VoronoiDiagram.cpp b/methods/VoronoiDiagram.cpp
index a31985b6710a86d2311dd0dc432e6b8ca9ef1834..6babd9be93816b7a69de9a810b2dfeceffe79b15 100644
--- a/methods/VoronoiDiagram.cpp
+++ b/methods/VoronoiDiagram.cpp
@@ -49,11 +49,39 @@ VoronoiDiagram::~VoronoiDiagram()
std::vector<std::pair<polygon_2d, int> > VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame, vector<double>& YInFrame,
vector<double>& VInFrame, vector<int>& IdInFrame, const double Bound_Max)
{
- const int numPedsInFrame = IdInFrame.size();
+ int numPedsInFrame = IdInFrame.size();
vector<int> XInFrame_temp;
vector<int> YInFrame_temp;
vector<double> VInFrame_temp;
vector<int> IdInFrame_temp;
+ // in case 1 or 2 pedestrians are in the geometry
+ // add dummy pedestrians around to enable voronoi calculations
+ // @todo: maybe use negative ids for these dummy pedestrians to exclude
+ // them from any analysis.
+ if(numPedsInFrame == 1 || numPedsInFrame == 2)
+ {
+ numPedsInFrame += 4;
+ // up right
+ XInFrame.push_back(XInFrame[0]+10);
+ YInFrame.push_back(YInFrame[0]+10);
+ VInFrame.push_back(VInFrame[0]);
+ IdInFrame.push_back(IdInFrame[0]+1);
+ // up left
+ XInFrame.push_back(XInFrame[0]-10);
+ YInFrame.push_back(YInFrame[0]+10);
+ VInFrame.push_back(VInFrame[0]);
+ IdInFrame.push_back(IdInFrame[0]+2);
+ // down right
+ XInFrame.push_back(XInFrame[0]+10);
+ YInFrame.push_back(YInFrame[0]-10);
+ VInFrame.push_back(VInFrame[0]);
+ IdInFrame.push_back(IdInFrame[0]+3);
+ // down left
+ XInFrame.push_back(XInFrame[0]-10);
+ YInFrame.push_back(YInFrame[0]-10);
+ VInFrame.push_back(VInFrame[0]);
+ IdInFrame.push_back(IdInFrame[0]+4);
+ }
for (int i = 0; i < numPedsInFrame; i++)
{
@@ -64,6 +92,7 @@ std::vector<std::pair<polygon_2d, int> > VoronoiDiagram::getVoronoiPolygons(vect
IdInFrame_temp.push_back(IdInFrame[i]);
}
+
VD voronoidiagram;
construct_voronoi(points.begin(), points.end(), &voronoidiagram);
int Ncell = 0;