diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..564f67366f2ec7f7d26db2df33a0ff074479cc50
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,87 @@
+# Created by .ignore support plugin (hsz.mobi)
+### C template
+# Object files
+*.o
+*.ko
+*.obj
+*.elf
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Libraries
+*.lib
+*.a
+*.la
+*.lo
+
+# Shared objects (inc. Windows DLLs)
+*.dll
+*.so
+*.so.*
+*.dylib
+
+# Executables
+*.exe
+*.out
+*.app
+*.i*86
+*.x86_64
+*.hex
+
+# Debug files
+*.dSYM/
+*.su
+### C++ template
+# Compiled Object files
+*.slo
+*.lo
+*.o
+*.obj
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Compiled Dynamic libraries
+*.so
+*.dylib
+*.dll
+
+# Fortran module files
+*.mod
+*.smod
+
+# Compiled Static libraries
+*.lai
+*.la
+*.a
+*.lib
+
+# Executables
+*.exe
+*.out
+*.app
+### CMake template
+CMakeCache.txt
+CMakeFiles
+CMakeScripts
+Makefile
+cmake_install.cmake
+install_manifest.txt
+CTestTestfile.cmake
+
+.gitignore
+.idea/
+Utest/sample_jpsreport/
+bin/
+build/
+lib/
+scripts/overlap.py
+scripts/results_rho_v_Voronoi_CROSSING_90_A_07_f_dir/
+scripts/rho_v_Voronoi_CROSSING_90_A_07_f_dir.txt_id_1.dat
+scripts/stefan copy/
+scripts/stefan/
+test.py
+/Release/
diff --git a/Analysis.cpp b/Analysis.cpp
index 06b6d57538a56293b5d4b8cba6f7a3c625a41951..5290b5de0c67fcda6e051b326b280ca330d3bc85 100644
--- a/Analysis.cpp
+++ b/Analysis.cpp
@@ -2,7 +2,7 @@
* \file Analysis.cpp
* \date Oct 10, 2014
* \version v0.7
- * \copyright <2009-2015> Forschungszentrum J��lich GmbH. All rights reserved.
+ * \copyright <2009-2016> Forschungszentrum Juelich GmbH. All rights reserved.
*
* \section License
* This file is part of JuPedSim.
@@ -316,12 +316,12 @@ int Analysis::RunAnalysis(const string& filename, const string& path)
}
if(false==IsInBuilding)
{
- Log->Write("Warning:\tAt %dth frame pedestrian at <x=%.4f, y=%.4f> is not in geometry!", frameNr+data.GetMinFrame(), XInFrame[i]*CMtoM, YInFrame[i]*CMtoM );
+ Log->Write("Warning:\tAt frame=%d pedestrian at <x=%.4f, y=%.4f> is outside the building!", frameNr+data.GetMinFrame(), XInFrame[i]*CMtoM, YInFrame[i]*CMtoM );
}
}
}
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------
+ //-----------------------------------------------------------------------------------------------------------------------------------------------------------------
if(_DoesUseMethodA) //Method A
{
#pragma omp parallel for
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 36fbd91ad6f41ae6b9143a3c5efb8f86786acd68..0b58cc4fef1e1423b0a34013e39325b2394b7f73 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,16 @@
# Change Log
All notable changes to `jpsreport` will be documented in this file.
+## v0.8.2 [unreleased]
+
+### Added
+
+- Added an error warning when the number of agents in the trajectory is not corresponding to total ids or the ped ids are not continuous in the first frame.
+
+### Changed
+
+### Fixed
+- Fixed SegFault due to reading files from different OS. (9a42c9dd)
+
## v0.8.1 [11.10.2016]
@@ -126,4 +137,4 @@ removing these pedestrians from the list.
### Fixed
-- Fixed bug for dealing with obstacles inside geometry.
+- Fixed bug for dealing with obstacles inside geometry.
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b7ae408c73b535bde4bec976f07e6f0a332d34c5..38a02bb044f19736fb1fbd3eddbcf5bbf09894a2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -36,7 +36,7 @@ message( STATUS "EXECUTABLE_OUTPUT_PATH: " ${EXECUTABLE_OUTPUT_PATH} )
message( STATUS "CMAKE_VERBOSE_MAKEFILE: " ${CMAKE_VERBOSE_MAKEFILE} )
if(NOT BUILD_TESTING)
- set(BUILD_TESTING ON) # test units are generated.
+ set(BUILD_TESTING OFF) # test units are generated.
endif(NOT BUILD_TESTING)
message( STATUS "BUILD_TESTING: " ${BUILD_TESTING} )
diff --git a/IO/OutputHandler.cpp b/IO/OutputHandler.cpp
index 4c9d19ee230206cb02639de4e294c35b05b89d98..b0e7f5c9f99e5c2d60e980c0ef515b7d89a64955 100644
--- a/IO/OutputHandler.cpp
+++ b/IO/OutputHandler.cpp
@@ -1,7 +1,7 @@
/**
* \file OutputHandler.cpp
* \date Nov 20, 2010
- * \version v0.5
+ * \version v0.8
* \copyright <2009-2015> Forschungszentrum Jülich GmbH. All rights reserved.
*
* \section License
diff --git a/geometry/Building.cpp b/geometry/Building.cpp
index 2a9b44080db7d772d8b3f663bf8c7c8b5512ad27..6722ea018fa07328d5e16a2e3d8ffc1cf9af496d 100644
--- a/geometry/Building.cpp
+++ b/geometry/Building.cpp
@@ -1,7 +1,7 @@
/**
* \file Building.cpp
* \date Oct 1, 2014
- * \version v0.7
+ * \version v0.8
* \copyright <2009-2015> Forschungszentrum Jülich GmbH. All rights reserved.
*
* \section License
diff --git a/main.cpp b/main.cpp
index 817b9d279e677402c5b12c0dd72b27793b7d6278..65a8328b36a2de56bd7fb37d42abea66a3eeab6a 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,8 +1,8 @@
/**
* \file main.cpp
* \date Oct 10, 2014
- * \version v0.7
- * \copyright <2009-2015> Forschungszentrum J��lich GmbH. All rights reserved.
+ * \version v0.8.1
+ * \copyright <2009-2015> Forschungszentrum Juelich GmbH. All rights reserved.
*
* \section License
* This file is part of JuPedSim.
diff --git a/methods/Method_A.cpp b/methods/Method_A.cpp
index 475a5113158e70a7d429c3d680adb07c43068dce..12ecf012649a78b809c492450ed913bbafbe957c 100644
--- a/methods/Method_A.cpp
+++ b/methods/Method_A.cpp
@@ -86,8 +86,8 @@ bool Method_A::Process (const PedData& peddata,const string& scriptsLocation, co
}
vector<int> ids=_peds_t[frameNr];
vector<int> IdInFrame = peddata.GetIdInFrame(frameNr, ids, zPos_measureArea);
- const vector<double> XInFrame = peddata.GetXInFrame(frameNr, ids, zPos_measureArea);
- const vector<double> YInFrame = peddata.GetYInFrame(frameNr, ids, zPos_measureArea);
+ // const vector<double> XInFrame = peddata.GetXInFrame(frameNr, ids, zPos_measureArea);
+ // const vector<double> YInFrame = peddata.GetYInFrame(frameNr, ids, zPos_measureArea);
const vector<double> VInFrame = peddata.GetVInFrame(frameNr, ids, zPos_measureArea);
if(IdInFrame.size()>0)
{
@@ -105,7 +105,7 @@ bool Method_A::Process (const PedData& peddata,const string& scriptsLocation, co
}
else
{
- Log->Write("Warning: No any pedestrian exists on the plane of the selected Measurement area!!");
+ Log->Write("Warning: No pedestrian exists on the plane of the selected Measurement area!!");
}
delete []_passLine;
return true;
diff --git a/methods/PedData.cpp b/methods/PedData.cpp
index 41e989382ff5fb26340d07ee21b747ff0858e7d3..fe222140c0787a7e92bf3bdbfef740943765a49c 100644
--- a/methods/PedData.cpp
+++ b/methods/PedData.cpp
@@ -106,7 +106,7 @@ bool PedData::InitializeVariables(const string& filename)
int pos_vd=0; //velocity direction
while ( getline(fdata,line) )
{
- //looking for the framerate which is suppposed to be at the second position
+ //looking for the framerate which is supposed to be at the second position
if(line[0] == '#')
{
std::vector<std::string> strs;
@@ -122,7 +122,7 @@ bool PedData::InitializeVariables(const string& filename)
{
std::vector<std::string> strs1;
line.erase(0,1);
- boost::split(strs1, line , boost::is_any_of("\t"),boost::token_compress_on);
+ boost::split(strs1, line , boost::is_any_of("\t "),boost::token_compress_on);
vector<string>::iterator it_id;
it_id=find(strs1.begin(),strs1.end(),"ID");
pos_id = std::distance(strs1.begin(), it_id);
@@ -143,6 +143,7 @@ bool PedData::InitializeVariables(const string& filename)
{
std::vector<std::string> strs;
+ boost::algorithm::trim_right(line);
boost::split(strs, line , boost::is_any_of("\t "),boost::token_compress_on);
if(strs.size() < 5)
{
@@ -252,7 +253,7 @@ bool PedData::InitializeVariables(const string& filename)
return true;
}
-// initialize the global variables variables
+// initialize the global variables
bool PedData::InitializeVariables(TiXmlElement* xRootNode)
{
if( ! xRootNode ) {
@@ -324,23 +325,29 @@ bool PedData::InitializeVariables(TiXmlElement* xRootNode)
}
//int frameNr=0;
for(TiXmlElement* xFrame = xRootNode->FirstChildElement("frame"); xFrame;
- xFrame = xFrame->NextSiblingElement("frame")) {
+ xFrame = xFrame->NextSiblingElement("frame"))
+ {
int frameNr = atoi(xFrame->Attribute("ID")) - _minFrame;
//todo: can be parallelized with OpenMP
for(TiXmlElement* xAgent = xFrame->FirstChildElement("agent"); xAgent;
- xAgent = xAgent->NextSiblingElement("agent")) {
-
+ xAgent = xAgent->NextSiblingElement("agent"))
+ {
//get agent id, x, y
double x= atof(xAgent->Attribute("x"));
double y= atof(xAgent->Attribute("y"));
double z= atof(xAgent->Attribute("z"));
int ID= atoi(xAgent->Attribute("ID"))-_minID;
-
+ if(ID>=_numPeds)
+ {
+ Log->Write("ERROR:\t The number of agents are not corresponding to IDs. Maybe Ped IDs are not continuous in the first frame, please check!!");
+ return false;
+ }
_peds_t[frameNr].push_back(ID);
_xCor[ID][frameNr] = x*M2CM;
_yCor[ID][frameNr] = y*M2CM;
_zCor[ID][frameNr] = z*M2CM;
+ // std::cout << "_zcor[ " << ID << "][" << frameNr <<"]=" << z << std::endl;
if(_vComponent == "F")
{
if(xAgent->Attribute("VD"))
@@ -370,7 +377,6 @@ bool PedData::InitializeVariables(TiXmlElement* xRootNode)
}
//frameNr++;
}
-
for(int id = 0; id<_numPeds; id++)
{
int actual_totalframe= totalframes[id];
@@ -396,7 +402,7 @@ vector<double> PedData::GetVInFrame(int frame, const vector<int>& ids, double zP
double v = GetInstantaneousVelocity1(frame, Tpast, Tfuture, id, _firstFrame, _lastFrame, _xCor, _yCor);
if(zPos<1000000.0)
{
- if(fabs(_zCor[id][frame]-zPos*M2CM)<J_EPS_EVENT)
+ if(fabs(_zCor[id][frame]-zPos*M2CM)<0.5*M2CM)
{
VInFrame.push_back(v);
}
@@ -416,7 +422,7 @@ vector<double> PedData::GetXInFrame(int frame, const vector<int>& ids, double zP
{
if(zPos<1000000.0)
{
- if(fabs(_zCor[id][frame]-zPos*M2CM)<J_EPS_EVENT)
+ if(fabs(_zCor[id][frame]-zPos*M2CM)<0.5*M2CM)
{
XInFrame.push_back(_xCor[id][frame]);
}
@@ -441,23 +447,23 @@ vector<double> PedData::GetXInFrame(int frame, const vector<int>& ids) const
vector<double> PedData::GetYInFrame(int frame, const vector<int>& ids, double zPos) const
{
- vector<double> YInFrame;
- for(unsigned int i=0; i<ids.size();i++)
- {
- int id = ids[i];
- if(zPos<1000000.0)
- {
- if(fabs(_zCor[id][frame]-zPos*M2CM)<J_EPS_EVENT)
- {
- YInFrame.push_back(_yCor[id][frame]);
- }
- }
- else
- {
- YInFrame.push_back(_yCor[id][frame]);
- }
- }
- return YInFrame;
+ vector<double> YInFrame;
+ for(unsigned int i=0; i<ids.size();i++)
+ {
+ int id = ids[i];
+ if(zPos<1000000.0)
+ {
+ if(fabs(_zCor[id][frame]-zPos*M2CM)<0.5*M2CM)
+ {
+ YInFrame.push_back(_yCor[id][frame]);
+ }
+ }
+ else
+ {
+ YInFrame.push_back(_yCor[id][frame]);
+ }
+ }
+ return YInFrame;
}
vector<double> PedData::GetYInFrame(int frame, const vector<int>& ids) const
@@ -500,10 +506,10 @@ vector<int> PedData::GetIdInFrame(int frame, const vector<int>& ids, double zPos
{
if(zPos<1000000.0)
{
- if(fabs(_zCor[id][frame]-zPos*M2CM)<J_EPS_EVENT)
- {
- IdInFrame.push_back(id +_minID);
- }
+ if(fabs(_zCor[id][frame]-zPos*M2CM)<0.5*M2CM)
+ {
+ IdInFrame.push_back(id +_minID);
+ }
}
else
{
diff --git a/scripts/2Dto1D_ID.py b/scripts/2Dto1D_ID.py
new file mode 100644
index 0000000000000000000000000000000000000000..6248b85ecec26f26b1f1abfe45b4a0950bcddcfb
--- /dev/null
+++ b/scripts/2Dto1D_ID.py
@@ -0,0 +1,161 @@
+from numpy import *
+import matplotlib
+import matplotlib.pyplot as plt
+import argparse
+from numpy import *
+import glob,os
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import matplotlib.cm as cm
+import datetime
+import shutil
+
+def getParserArgs():
+ parser = argparse.ArgumentParser(description='Combine French data to one file')
+ parser.add_argument("-N", "--pedNum", type=int , default=16, help='The number of pedestrians in the run (default 16)')
+ parser.add_argument("-r", "--runName", default="XP01", help='give the name of each run')
+ args = parser.parse_args()
+ return args
+
+def get_dateString():
+ date = str(datetime.datetime.today()) #
+ date = date.replace('-', ' ')
+ date = date.replace(':', ' ')
+ date = date.replace('.', ' ')
+ date = date.split()
+ date[0] = str(int(date[0])-2000)
+ dateString = str(date[0] + date[1] + date[2] + '_' + date[3] + date[4] + date[5])
+
+ return dateString
+
+def get_archive(dateString):
+ os.makedirs('%s_traj_OneD_V' %(dateString))
+ os.makedirs('%s_traj_OneD_V/INPUT' %(dateString))
+ if not os.path.exists('%s_traj_OneD_V/OUTPUT/'%(dateString)):
+ os.makedirs('%s_traj_OneD_V/OUTPUT' %(dateString))
+
+
+def get_backupInput(dateString):
+
+ # Sichern der INPUT-Dateien
+ #==========================================================================
+ dst = '%s_traj_OneD_V/INPUT' %(dateString)
+ src = '2Dto1D_ID.py'
+ shutil.copy(src, dst)
+ print ('--> Gesichert: 2Dto1D_ID.py')
+
+def traj_TwoD2OneD(trajectory):
+ traj=loadtxt(trajectory)
+ traj[:,2]/=100.0
+ traj[:,3]/=100.0
+ lastpos=0
+ lastID=1
+ data=array([[0,0,0,0]])
+ ID=lastID
+ for i in range(shape(traj)[0]):
+ if(traj[i,2]>-2.5 and traj[i,2]<2.5 and traj[i,3]<0):
+ posx=traj[i,2]+2.5
+ posy=0.5*(2.1+2.9)-fabs(traj[i,3])
+ elif(traj[i,2]>-2.5 and traj[i,2]<2.5 and traj[i,3]>0):
+ posx=12.85+2.5-traj[i,2]
+ posy=0.5*(2.1+2.9)-fabs(traj[i,3])
+ elif(traj[i,2]>=2.5):
+ #print(traj[i,2])
+ dis=square(traj[i,3]+2.5)+square(traj[i,2]-2.5)
+ #print(dis)
+ #print((6.25+square(traj[i,3])+square(traj[i,2])-dis)/12.5)
+ theta=arccos((6.25+square(traj[i,3])+square(traj[i,2]-2.5)-dis)/(5*sqrt(square(traj[i,3])+square(traj[i,2]-2.5))))
+ posx=5+theta*2.5
+ disy=sqrt(pow(traj[i,2]-2.5,2)+pow(traj[i,3],2))
+ posy=0.5*(2.1+2.9)-disy
+
+ elif(traj[i,2]<=-2.5):
+ #print(traj[i,2])
+ dis=square(traj[i,3]-2.5)+square(traj[i,2]+2.5)
+ theta=arccos((6.25+square(traj[i,3])+square(traj[i,2]+2.5)-dis)/(5*sqrt(square(traj[i,3])+square(traj[i,2]+2.5))))
+ posx=17.85+theta*2.5
+ disy=sqrt(pow(traj[i,2]+2.5,2)+pow(traj[i,3],2))
+ posy=0.5*(2.1+2.9)-disy
+# if(fabs(lastpos-posx)>25):
+# ID+=1
+# elif (traj[i,0]!=lastID):
+ if (traj[i,0]!=lastID):
+ ID+=1
+ lastpos=posx
+ lastID=traj[i,0]
+ data=append(data,array([[ID,traj[i,1],posx,posy]]),axis=0)
+ data=data[1:]
+ return data
+ #savetxt("traj_ring_adults_06_06.txt",data, fmt='%d %d %.3f %.3f', delimiter='\t')
+
+def add_instantneousV(file, fps=25, timeinterval=10):
+ data=traj_TwoD2OneD(file)
+ deltaT=int(timeinterval/2)
+ V = array([[0,0]])
+ V=delete(V,0,0)
+ for j in unique(data[:,0]):
+ trajectory=data[data[:,0]==j]
+ Frames=trajectory[:,1]
+ T=trajectory[:,1]/fps
+ Sx=trajectory[:,2]
+ Sy=trajectory[:,3]
+ for t in range(len(Frames)):
+ if(len(Frames)>timeinterval):
+ if (t < deltaT):
+ if (Sx[t+deltaT] - Sx[t])<-20:
+ vx = (Sx[t+deltaT] - Sx[t]+25.7)/(T[t+deltaT]-T[t])
+ elif (Sx[t+deltaT] - Sx[t])>20:
+ vx = (Sx[t+deltaT] - Sx[t]-25.7)/(T[t+deltaT]-T[t])
+ else:
+ vx = (Sx[t+deltaT] - Sx[t])/(T[t+deltaT]-T[t])
+ vy = (Sy[t+deltaT] - Sy[t])/(T[t+deltaT]-T[t])
+ elif (t > len(Frames)-deltaT-1):
+ if (Sx[t] - Sx[t-deltaT])<-20:
+ vx = (Sx[t] - Sx[t-deltaT]+25.7)/(T[t]-T[t-deltaT])
+ elif (Sx[t] - Sx[t-deltaT])>20:
+ vx = (Sx[t] - Sx[t-deltaT]-25.7)/(T[t]-T[t-deltaT])
+ else:
+ vx = (Sx[t] - Sx[t-deltaT])/(T[t]-T[t-deltaT])
+ vy = (Sy[t] - Sy[t-deltaT])/(T[t]-T[t-deltaT])
+ else:
+ if (Sx[t+deltaT] - Sx[t-deltaT])<-20:
+ vx = (Sx[t+deltaT] - Sx[t-deltaT]+25.7)/(T[t+deltaT]-T[t-deltaT])
+ elif (Sx[t+deltaT] - Sx[t-deltaT])>20:
+ vx = (Sx[t+deltaT] - Sx[t-deltaT]-25.7)/(T[t+deltaT]-T[t-deltaT])
+ else:
+ vx = (Sx[t+deltaT] - Sx[t-deltaT])/(T[t+deltaT]-T[t-deltaT])
+ vy = (Sy[t+deltaT] - Sy[t-deltaT])/(T[t+deltaT]-T[t-deltaT])
+ else:
+ vx=0
+ vy=0
+ V=append(V,array([[vx,vy]]),axis=0)
+ data=append(data,V,axis=1)
+ return data
+
+
+def traj_TwoD2OneD_V(directory, dateString):
+ files = glob.glob(directory+"/*students*.txt")
+ for file in files:
+ print(file)
+ oneD_traj_file=dateString+"_traj_OneD_V/OUTPUT/OneD_traj_"+file.split("ring_")[1].split("_shifted")[0]+".txt"
+ header="#description: %s\n#framerate: %d\n#ID: the agent ID\n#FR: the current frame\n#X,Y,v: the agents coordinates (in metres) and instantneous velocity (m/s)\n\n#ID\tFR\tX\tY\tVx\tVy"%(file,fps)
+ data=add_instantneousV(file)
+ savetxt(oneD_traj_file, data,fmt='%d %d %.3f %.3f %.3f %.3f', header=header, comments="", delimiter='\t',newline='\r\n')
+
+
+if __name__ == '__main__':
+ args = getParserArgs()
+ Npeds=args.pedNum
+ runName=args.runName
+ fps=25.0
+ dateString = get_dateString()
+ get_archive(dateString)
+ get_backupInput(dateString)
+ traj_TwoD2OneD_V("./_trajectory",dateString)
+# files = glob.glob("./_traj_OneD/stop/*.txt")
+# for file in files:
+# print(file)
+# data=loadtxt(file)
+# figname=dateString+"_TimeSpace/OUTPUT/TS_"+file.split("OneD_traj_")[1].split(".")[0]
+# print(figname)
+# plot_VTimeSpaceDiagram2(data,figname)
+
diff --git a/scripts/Readme.txt b/scripts/Readme.txt
new file mode 100644
index 0000000000000000000000000000000000000000..6b8c6db49153b96bedbbfede64286d66f6f6b8fc
--- /dev/null
+++ b/scripts/Readme.txt
@@ -0,0 +1,29 @@
+_getGait.py: (savitzky_golay)
+ Get gait information directly from 2D trajectory by ditecting each local peaks and valleys from the smoothed curvature information.
+
+txt2xml.py (Also smoothing code inside(Jack & Bernhard)_zero curvature spline)
+ Transfer the txt format trajectory file to xml format, which can be used for JPSreport and JPSvis
+
+_plot_cell_rho.py:
+ Plot geometry, pedestrian as well as the voronoi diagrams. The pedestrian density in each diagram is shown by different colors.
+
+_plot_cell_v.py:
+ Plot geometry, pedestrian as well as the voronoi diagrams. The velocity for each pedestrian in the diagram is shown by different colors.
+
+_plot_FD.py:
+ Plot fundamental diagrams from the output of JPSreport
+
+_plot_N_t.py:
+ Plot the N-t diagram obtained from method A.
+
+_Plot_profiles.py
+ Plot density, velocity and specific profiles based on the output from JPSreport
+
+_SteadyState.py
+ detect the begin and the end of steady state from time series of density and velocity.
+
+txt2txt.py
+ tranfer .txt format trajectory files to meet the requirement of JPSreport
+
+2Dto1D_ID.py
+ Transfer 2D trajectory information to 1D information. This script is only used for single file experiment in a ring.
\ No newline at end of file
diff --git a/scripts/_Plot_N_t.py b/scripts/_Plot_N_t.py
index ed8480a2459ca4d588385c60f7d86a913326d715..52f74bb8716bf9783f51667f53c2c39376388e2d 100644
--- a/scripts/_Plot_N_t.py
+++ b/scripts/_Plot_N_t.py
@@ -1,9 +1,6 @@
-from numpy import *
-import matplotlib
+import numpy as np
import matplotlib.pyplot as plt
from Polygon import *
-import matplotlib.cm as cm
-import pylab
import argparse
import sys
@@ -29,11 +26,11 @@ if __name__ == '__main__':
ax1 = fig.add_subplot(111,aspect='auto')
plt.rc("font", size=30)
plt.rc('pdf',fonttype = 42)
- data_NT = loadtxt("%s/%s"%(pathfile,namefile))
+ data_NT = np.loadtxt("%s/%s"%(pathfile,namefile))
#plt.plot(data_NT[:,0]/fps,data_NT[:,1], 'r-')
plt.plot(data_NT[:,0],data_NT[:,1], 'r-')
- plt.xlabel("t [s]")
- plt.ylabel("N [-]")
+ plt.xlabel("t [s]", size=18)
+ plt.ylabel("N [-]", size=18)
plt.title("%s"%title)
plt.savefig("%s/%s.png"%(pathfile,figname))
plt.close()
diff --git a/scripts/_getGait.py b/scripts/_getGait.py
new file mode 100644
index 0000000000000000000000000000000000000000..d7ac7daa2a00f1b447b8d0955ba198fe4eca3c17
--- /dev/null
+++ b/scripts/_getGait.py
@@ -0,0 +1,203 @@
+from numpy import *
+import matplotlib
+import matplotlib.pyplot as plt
+import argparse
+from numpy import *
+import glob,os
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import matplotlib.cm as cm
+from matplotlib.patches import Arc
+from scipy.signal import argrelextrema
+
+def getParserArgs():
+ parser = argparse.ArgumentParser(description='Combine French data to one file')
+ parser.add_argument("-N", "--pedNum", type=int , default=16, help='The number of pedestrians in the run (default 16)')
+ parser.add_argument("-r", "--runName", default="XP01", help='give the name of each run')
+ args = parser.parse_args()
+ return args
+
+def estimated_autocorrelation(x):
+ n = len(x)
+ variance = x.var()
+ x = x-x.mean()
+ r = correlate(x, x, mode = 'full')[-n:]
+ assert allclose(r, array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)]))
+ result = r/(variance*(arange(n, 0, -1)))
+ return result
+
+def savitzky_golay(y, window_size, order, deriv=0, rate=1):
+ r"""Smooth (and optionally differentiate) data with a Savitzky-Golay filter.
+ The Savitzky-Golay filter removes high frequency noise from data.
+ It has the advantage of preserving the original shape and
+ features of the signal better than other types of filtering
+ approaches, such as moving averages techniques.
+ Parameters
+ ----------
+ y : array_like, shape (N,)
+ the values of the time history of the signal.
+ window_size : int
+ the length of the window. Must be an odd integer number.
+ order : int
+ the order of the polynomial used in the filtering.
+ Must be less then `window_size` - 1.
+ deriv: int
+ the order of the derivative to compute (default = 0 means only smoothing)
+ Returns
+ -------
+ ys : ndarray, shape (N)
+ the smoothed signal (or it's n-th derivative).
+ Notes
+ -----
+ The Savitzky-Golay is a type of low-pass filter, particularly
+ suited for smoothing noisy data. The main idea behind this
+ approach is to make for each point a least-square fit with a
+ polynomial of high order over a odd-sized window centered at
+ the point.
+ Examples
+ --------
+ t = np.linspace(-4, 4, 500)
+ y = np.exp( -t**2 ) + np.random.normal(0, 0.05, t.shape)
+ ysg = savitzky_golay(y, window_size=31, order=4)
+ import matplotlib.pyplot as plt
+ plt.plot(t, y, label='Noisy signal')
+ plt.plot(t, np.exp(-t**2), 'k', lw=1.5, label='Original signal')
+ plt.plot(t, ysg, 'r', label='Filtered signal')
+ plt.legend()
+ plt.show()
+ References
+ ----------
+ .. [1] A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of
+ Data by Simplified Least Squares Procedures. Analytical
+ Chemistry, 1964, 36 (8), pp 1627-1639.
+ .. [2] Numerical Recipes 3rd Edition: The Art of Scientific Computing
+ W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery
+ Cambridge University Press ISBN-13: 9780521880688
+ """
+ import numpy as np
+ from math import factorial
+
+ try:
+ window_size = np.abs(np.int(window_size))
+ order = np.abs(np.int(order))
+ except (ValueError, msg):
+ raise ValueError("window_size and order have to be of type int")
+ if window_size % 2 != 1 or window_size < 1:
+ raise TypeError("window_size size must be a positive odd number")
+ if window_size < order + 2:
+ raise TypeError("window_size is too small for the polynomials order")
+ order_range = range(order+1)
+ half_window = (window_size -1) // 2
+ # precompute coefficients
+ b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
+ m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
+ # pad the signal at the extremes with
+ # values taken from the signal itself
+ firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
+ lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
+ y = np.concatenate((firstvals, y, lastvals))
+ return np.convolve( m[::-1], y, mode='valid')
+
+
+if __name__ == '__main__':
+ args = getParserArgs()
+ Npeds=args.pedNum
+ runName=args.runName
+ print(runName)
+ #ID=0
+ file="./ring_adults_07_05/ring_adults_07_05_hdv_shifted.txt"
+ traj=loadtxt(file)
+ Persons=unique(traj[:,0])
+ freq_size=zeros((1,5))
+ fps = 25
+ for pers in Persons:
+ data=traj[traj[:,0]==pers]
+ v_x=copy(data[5:-5,2])
+ v_y=copy(data[5:-5,3])
+ a_x = copy(v_x[5:-5])
+ a_y = copy(v_y[5:-5])
+ T_Curvature = zeros((len(a_x),2))
+ C = copy(v_y[5:-5])
+ frame_step = 10 ## frame_step should be even
+ window_size = 41
+ order = 3
+ figname=file+"_pers_%d"%(pers)+".png"
+ step_frequency=zeros((1,5))
+ #----------------------claculating velocity, accelaration and curvature--------------------------
+ for i in range(len(v_x)):
+ v_x[i]=fps*(data[i+frame_step,2]-data[i,2])/(frame_step*100.)
+ v_y[i]=fps*(data[i+frame_step,3]-data[i,3])/(frame_step*100.)
+ for j in range(len(a_x)):
+ a_x[j]=fps*(v_x[j+frame_step]-v_x[j])/frame_step
+ a_y[j]=fps*(v_y[j+frame_step]-v_y[j])/frame_step
+ T_Curvature[j,0]=j
+ T_Curvature[j,1]= (v_x[j+5]*a_y[j]-v_y[j+5]*a_x[j])/pow(sqrt(pow(v_x[j+5],2)+pow(v_y[j+5],2)),3)
+
+
+ #-------------------------Smooth the time series of curvatures----------------------------
+ smoothed_curv = savitzky_golay(T_Curvature[:,1], window_size, order)
+ peak_id=argrelextrema(smoothed_curv, greater_equal, mode='wrap')[0]
+ peak=smoothed_curv[peak_id]
+ valley_id=argrelextrema(smoothed_curv, less_equal, mode='wrap')[0]
+ valley=smoothed_curv[valley_id]
+ fig = plt.figure(figsize=(13, 17), dpi=100)
+ fig.subplots_adjust(hspace=0.4)
+ ax1 = fig.add_subplot(211,aspect='equal')
+ plt.rc("font", size=40)
+ ax1.plot(data[:,2]/100.,data[:,3]/100., 'r-', lw=0.4)
+ for i in peak_id:
+ pmax=data[data[:,1]==(i+frame_step/2)]
+ ax1.plot(pmax[:,2]/100.,pmax[:,3]/100., 'bo', lw=0.4)
+ step_frequency=append(step_frequency,pmax,axis=0)
+ for j in valley_id:
+ pmin=data[data[:,1]==(j+frame_step/2)]
+ ax1.plot(pmin[:,2]/100.,pmin[:,3]/100., 'go', lw=0.4)
+ step_frequency=append(step_frequency,pmin,axis=0)
+ ax1.set_ylim(-4,4)
+ ax1.set_xlim(-6,6)
+ fs=40
+ plt.xlabel("$x\, [m]$", fontsize=fs)
+ plt.ylabel("$y\, [m]$", fontsize=fs)
+ ax0 = fig.add_subplot(212,aspect='auto')
+ ax0.plot(T_Curvature[:,0]/fps,T_Curvature[:,1], 'ro', markeredgecolor='r',lw=0.4,markersize=3)
+ ax0.plot(T_Curvature[:,0]/fps,smoothed_curv, 'g-', lw=0.8)
+ ax0.plot(peak_id/fps,peak, 'bo', lw=0.8,markersize=8)
+ ax0.plot(valley_id/fps,valley, 'bo', lw=0.8,markersize=8)
+ plt.xlabel("$t\, [s]$", fontsize=fs)
+ plt.ylabel("$Curvature$", fontsize=fs)
+ plt.savefig(figname)
+ plt.close()
+ step_frequency=delete(step_frequency,0,axis=0)
+ step_frequency=step_frequency[argsort(step_frequency[:,1])]
+ fstep_frequency=file+"_gait_Pers_%d"%(pers)+".dat"
+ freq_size=append(freq_size,diff(step_frequency,axis=0),axis=0)
+ savetxt(fstep_frequency,step_frequency, fmt='%d\t%d\t%.3f\t%.3f\t%.3f',delimiter='\t',newline='\r\n')
+
+ fig = plt.figure(figsize=(15, 12), dpi=100)
+ ax2 = fig.add_subplot(211,aspect='auto')
+ ax2.set_ylim(-2,2)
+ plt.xlabel("$t\, [s]$", fontsize=fs)
+ plt.ylabel("$Autocorrelation$", fontsize=fs)
+ figname=file+"_gait_autocorrelation_Pers_%d"%(pers)+".png"
+ plt.plot(T_Curvature[:,0]/fps,estimated_autocorrelation(T_Curvature[:,1]),'r-')
+ fcurvature_name=file+"_curvature_%d"%(pers)+".dat"
+ savetxt(fcurvature_name,T_Curvature, fmt='%.3f\t%.3f',delimiter='\t',newline='\r\n')
+ plt.savefig(figname)
+ plt.close()
+
+ freq_size=delete(freq_size,0,axis=0)
+ fig = plt.figure(figsize=(15, 12), dpi=100)
+ ax3 = fig.add_subplot(111,aspect='auto')
+ ax3.set_ylim(0,2)
+ ax3.set_xlim(0,5)
+ plt.xlabel("$t\, [s]$", fontsize=fs)
+ plt.ylabel("$step size [m]$", fontsize=fs)
+ figname=file+"_freq_stepsize.png"
+ plt.plot(freq_size[:,1]/fps,sqrt(pow(freq_size[:,2],2)+pow(freq_size[:,3],2))/100.,'ro')
+ freq_size[:,0]=freq_size[:,1]/fps
+ freq_size[:,1]=sqrt(pow(freq_size[:,2],2)+pow(freq_size[:,3],2))/100.
+ freq_size[:,2]=freq_size[:,1]/freq_size[:,0]
+ freq_size=delete(freq_size,[3,4],axis=1)
+ ffreq_size=file+"_freq_size.dat"
+ savetxt(ffreq_size,freq_size, fmt='%.3f\t%.3f\t%.3f',delimiter='\t',newline='\r\n')
+ plt.savefig(figname)
+ plt.close()
diff --git a/scripts/plot_voronoi.py b/scripts/plot_voronoi.py
new file mode 100644
index 0000000000000000000000000000000000000000..6bcdb569707e15fa0fa44c92384114320d7aee11
--- /dev/null
+++ b/scripts/plot_voronoi.py
@@ -0,0 +1,154 @@
+from numpy import *
+import matplotlib
+import matplotlib.pyplot as plt
+from matplotlib.patches import Polygon as pgon
+from Polygon import *
+import matplotlib.cm as cm
+import pylab
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import argparse
+import sys
+from utils import *
+import os
+def getParserArgs():
+ parser = argparse.ArgumentParser(description='Combine French data to one file')
+ parser.add_argument("-f", "--filepath", default="./", help='give the path of source file')
+ parser.add_argument("-n", "--namefile", help='give the name of the source file')
+ parser.add_argument("-x1", "--geominx", type=float, help='give the minmum x of the geometry')
+ parser.add_argument("-x2", "--geomaxx", type=float, help='give the maxmum x of the geometry')
+ parser.add_argument("-y1", "--geominy", type=float, help='give the minmum y of the geometry')
+ parser.add_argument("-y2", "--geomaxy", type=float, help='give the maxmum y of the geometry')
+ args = parser.parse_args()
+ return args
+
+
+try:
+ import xml.etree.cElementTree as ET
+except ImportError:
+ import xml.etree.ElementTree as ET
+
+
+def get_polygon(poly):
+ X = []
+ Y = []
+ # for poly in node.getchildren():
+ # # print "polygon tag: ", poly.tag
+ # # if poly.tag == "obstacle":
+ # # # print poly.getchildren()
+ # # # for pobst in poly.getchildren():
+ # # # # print pobst.tag
+ # # # for q in pobst.getchildren(): # vertex
+ # # # X.append( q.attrib['px'] )
+ # # # Y.append( q.attrib['py'] )
+ # # pass
+ # # else:
+ for p in poly.getchildren(): # vertex
+ X.append(p.attrib['px'])
+ Y.append(p.attrib['py'])
+
+ return X, Y
+
+def plot_geometry(geometry):
+ tree = ET.parse(geometry)
+ root = tree.getroot()
+ for node in root.iter():
+ tag = node.tag
+ # print "subroom tag", tag
+ if tag == "polygon":
+ X, Y = get_polygon(node)
+ plt.plot(X, Y, "k", lw=4)
+ elif tag == "obstacle":
+ # print "obstacle tag",tag
+ for n in node.getchildren():
+ # print "N: ", n
+ X, Y = get_polygon(n)
+ # print "obstacle", X, Y
+ plt.plot(X, Y, "g", lw=4)
+ elif tag == "crossing":
+ X, Y = get_polygon(node)
+ plt.plot(X, Y, "--b", lw=0.9, alpha=0.2)
+ elif tag == "HLine":
+ X, Y = get_polygon(node)
+ plt.plot(X, Y, "--b", lw=0.9, alpha=0.2)
+ elif tag == "transition":
+ X, Y = get_polygon(node)
+ plt.plot(X, Y, "--r", lw=1.6, alpha=0.2)
+
+
+
+if __name__ == '__main__':
+ rho_max = 8.0
+ args = getParserArgs()
+ filepath = args.filepath
+ sys.path.append(filepath)
+ namefile = args.namefile
+ geominX = args.geominx
+ geomaxX = args.geomaxx
+ geominY = args.geominy
+ geomaxY = args.geomaxy
+ fig = plt.figure(figsize=(16*(geomaxX-geominX)/(geomaxY-geominY)+2, 16), dpi=100)
+ ax1 = fig.add_subplot(111,aspect='equal')
+ plt.rc("font", size=30)
+ plt.rc('pdf',fonttype = 42)
+ ax1.set_yticks([int(1*j) for j in range(-2,5)])
+ for label in ax1.get_xticklabels() + ax1.get_yticklabels():
+ label.set_fontsize(30)
+ for tick in ax1.get_xticklines() + ax1.get_yticklines():
+ tick.set_markeredgewidth(2)
+ tick.set_markersize(6)
+ ax1.set_aspect("equal")
+ density=array([])
+ polys = open("%s/polygon%s.dat"%(filepath,namefile)).readlines()
+ for poly in polys:
+ exec("p = %s"%poly)
+ #p=tuple(tuple(i/100.0 for i in inner) for inner in p)
+ xx=1.0/Polygon(p).area()
+ if xx>rho_max:
+ xx=rho_max
+ density=append(density,xx)
+ Maxd=density.max()
+ Mind=density.min()
+ erro=ones(shape(density))*Mind
+ density=rho_max*(density-erro)/(Maxd-Mind)
+ sm = cm.ScalarMappable(cmap=cm.jet)
+ sm.set_array(density)
+ sm.autoscale_None()
+ sm.set_clim(vmin=0,vmax=5)
+ for poly in polys:
+ exec("p = %s"%poly)
+ #p=tuple(tuple(i/100.0 for i in inner) for inner in p)
+ xx=1.0/Polygon(p).area()
+ if xx>rho_max:
+ xx=rho_max
+ ax1.add_patch(pgon(p,facecolor=sm.to_rgba(xx), edgecolor='white',linewidth=2))
+ frameNr = int(namefile.split("_")[-1].split(".")[0])
+
+ print "PWD=", os.getcwd()
+ trajFile = "../Andjiana/hybrid_hall_traj.xml"
+ geoFile = "../Andjiana/hybrid_hall_geo.xml"
+ plot_geometry(geoFile)
+ fps, N, Trajdata = parse_file(trajFile)
+ print Trajdata
+ Trajdata = Trajdata[ Trajdata[:, 1] == frameNr ]
+ #points = loadtxt("%s/points%s.dat"%(filepath,namefile))
+ print "--------------------"
+ print Trajdata
+ ax1.plot(Trajdata[:,2], Trajdata[:,3], "bo", markersize = 20, markeredgewidth=2)
+
+ #ax1.plot(points[:,0],points[:,1],"bo",markersize = 20,markeredgewidth=2)
+ ax1.set_xlim(geominX,geomaxX)
+ ax1.set_ylim(geominY,geomaxY)
+ plt.xlabel("x [m]")
+ plt.ylabel("y [m]")
+ #print density
+ plt.title("%s"%namefile)
+ divider = make_axes_locatable(ax1)
+ cax = divider.append_axes("right", size="2.5%", pad=0.2)
+ cb=fig.colorbar(sm,ax=ax1,cax=cax,format='%.1f')
+ cb.set_label('Density [$m^{-2}$]')
+ plt.savefig("%s/rho_%s.png"%(filepath,namefile))
+ print "figname: %s/rho_%s.png"%(filepath,namefile)
+ plt.show()
+ plt.close()
+
+
diff --git a/scripts/runscript.py b/scripts/runscript.py
new file mode 100644
index 0000000000000000000000000000000000000000..cf181239a142e423dfffc1768da10257cf10d7c2
--- /dev/null
+++ b/scripts/runscript.py
@@ -0,0 +1,12 @@
+import glob
+
+DIR = "../Andjiana/90deg/Output/Fundamental_Diagram/Classical_Voronoi/VoronoiCell"
+
+files = glob.glob("%s/hybrid_hall_traj*.dat"%DIR)
+
+for f in files:
+ ff = f.split(".dat")[0]
+ cmd = python plot_voronoi.py -f DIR -n ff -x1 -2 -x2 14 -y1 -2 -y2 14
+ subprocess.call(["python", ])
+
+
diff --git a/scripts/txt2txt.py b/scripts/txt2txt.py
new file mode 100644
index 0000000000000000000000000000000000000000..16a96d174a7831e172b5f5c6d90872332cbe2d46
--- /dev/null
+++ b/scripts/txt2txt.py
@@ -0,0 +1,34 @@
+#!/usr/bin/python
+# format id frame x y z
+from numpy import *
+import os, sys, glob, math
+import argparse
+
+
+def getParserArgs():
+ parser = argparse.ArgumentParser(description='Tranfer the unit from cm to m and add header for .txt file')
+ parser.add_argument("-p", "--path", default="./", help='give the path of source file')
+ parser.add_argument("-f", "--fps", default="16", type=float, help='give the frame rate of data')
+ args = parser.parse_args()
+ return args
+
+if __name__ == '__main__':
+ args = getParserArgs()
+ path = args.path
+ print(path)
+ sys.path.append(path)
+ fps = args.fps
+ files = glob.glob("%s/*.txt"%(path))
+ for file in files:
+ description="experiment"
+ geometry="geometry.xml"
+ print(os.path.splitext(file)[0])
+ fname = os.path.splitext(file)[0]+".txt"
+ header="#description: %s\n#framerate: %d\n#geometry: %s\n#ID: the agent ID\n#FR: the current frame\n#X,Y,Z: the agents coordinates (in metres)\n\n#ID\tFR\tX\tY\tZ"%(description,fps,geometry)
+ data=loadtxt(file,usecols = (0,1,2,3,4))
+ data[:,2]/=100.
+ data[:,3]/=100.
+ data[:,4]/=100.
+ #np.savetxt("./uo-100-300-300_space.txt",data, fmt= "%d %d %.4f %.4f %.4f", delimiter ='\t', header=header, comments='',newline='\n')
+ savetxt(fname,data, fmt= "%d\t%d\t%.4f\t%.4f\t%.4f", delimiter ='\t', header=header, comments='', newline='\r\n')
+
diff --git a/scripts/xml2txt.py b/scripts/xml2txt.py
new file mode 100644
index 0000000000000000000000000000000000000000..b6f996f50b378cc7ff5dbedb9359405664858f00
--- /dev/null
+++ b/scripts/xml2txt.py
@@ -0,0 +1,70 @@
+"""
+convert jpscore's trajectories from xml to txt format:
+id fr x y z
+input: xml file
+output: xml file
+"""
+
+from sys import argv, exit
+import os
+try:
+ import xml.etree.cElementTree as ET
+except ImportError:
+ import xml.etree.ElementTree as ET
+
+
+def write_header(o, fps, geometry_file):
+ o.write("#xml to txt converted simulation file\n")
+ o.write("#framerate: %0.2f\n" % fps)
+ o.write("#geometry: %s\n" % geometry_file)
+ o.write("#ID: the agent ID\n");
+ o.write("#FR: the current frame\n");
+ o.write("#X,Y,Z: the agents coordinates (in meters)\n");
+ o.write("\n");
+ o.write("#ID\tFR\tX\tY\tZ\n");
+
+
+
+if __name__ == "__main__":
+ if len(argv) < 2:
+ exit('usage: %s filename' % argv[0])
+
+ filename = argv[1]
+ filename1, file_extension = os.path.splitext(filename)
+ if file_extension != ".xml":
+ exit('not an xml file')
+
+ output = filename1 + ".txt"
+ o = open(output, "w")
+ # o.write("# ID frame x y z\n\n")
+ print (">> %s" % output)
+
+ tree = ET.parse(filename)
+ root = tree.getroot()
+
+ for header in root.iter('header'):
+ fps = header.find('frameRate').text
+
+ try:
+ fps = float(fps)
+ except:
+ print ("ERROR: could not read <fps>")
+ exit()
+
+ for g in root.iter('geometry'):
+ file_location = g.find('file').attrib
+ location = file_location['location']
+
+ write_header(o, fps, location)
+ for node in root.iter():
+ tag = node.tag
+ if tag == "frame":
+ frame = node.attrib['ID']
+ for agent in node.getchildren():
+ x = agent.attrib['x']
+ y = agent.attrib['y']
+ z = agent.attrib['z']
+ ID = agent.attrib['ID']
+ o.write("%d\t %d\t %.5f\t %.5f\t %.5f\n" % (int(ID), int(frame), float(x), float(y), float(z)))
+
+ o.close()