diff --git a/CMakeLists.txt b/CMakeLists.txt
index 637671e72c0ed8769fd211bf6e1b0f209f43892e..1f6ce515dfcd1e28476c6aed5d43e957b4b5bb32 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -180,6 +180,7 @@ set(methods
   methods/Method_B.cpp
   methods/Method_C.cpp
   methods/Method_D.cpp
+  methods/Method_I.cpp
   )
 set(source_files
   Analysis.cpp
@@ -211,6 +212,7 @@ set (  header_files
   methods/Method_B.h
   methods/Method_C.h
   methods/Method_D.h
+  methods/Method_I.h
   IO/OutputHandler.h
   general/ArgumentParser.h
   general/Macros.h
diff --git a/methods/Method_I.cpp b/methods/Method_I.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3f65c2c3011cd177cde7dc75a005636bbc996229
--- /dev/null
+++ b/methods/Method_I.cpp
@@ -0,0 +1,816 @@
+/**
+ * \file        Method_I.cpp
+ * \date        Feb 07, 2019
+ * \version     v0.8
+ * \copyright   <2009-2019> Forschungszentrum Jülich GmbH. All rights reserved.
+ *
+ * \section License
+ * This file is part of JuPedSim.
+ *
+ * JuPedSim is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * any later version.
+ *
+ * JuPedSim is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with JuPedSim. If not, see <http://www.gnu.org/licenses/>.
+ *
+ * \section Description
+ * In this file functions related to method I are defined.
+ *
+ *
+ **/
+
+#include "Method_I.h"
+#include <cmath>
+#include<map>
+#include<iostream>
+#include<vector>
+#include <tuple>
+//using std::string;
+//using std::vector;
+//using std::ofstream;
+using namespace std;
+
+
+
+Method_I::Method_I()
+{
+     _grid_size_X = 0.10;
+     _grid_size_Y = 0.10;
+     _fps=16;
+     _outputVoronoiCellData = false;
+     _getProfile = false;
+     _geoMinX = 0;
+     _geoMinY = 0;
+     _geoMaxX = 0;
+     _geoMaxY = 0;
+     _cutByCircle = false;
+     _cutRadius = -1;
+     _circleEdges = -1;
+     _fIndividualFD = nullptr;
+     _calcIndividualFD = false;
+     _fVoronoiRhoV = nullptr;
+     _areaForMethod_I = nullptr;
+     _plotVoronoiCellData=false;
+     _isOneDimensional=false;
+     _startFrame =-1;
+     _stopFrame = -1;
+}
+
+Method_I::~Method_I()
+{
+
+}
+// auto outOfRange = [&](int number, int start, int end) -> bool
+// {
+//      return (number < start || number > end);
+// };
+bool Method_I::Process (const PedData& peddata,const fs::path& scriptsLocation, const double& zPos_measureArea)
+{
+     bool return_value = true;
+     _scriptsLocation = scriptsLocation;
+     _outputLocation = peddata.GetOutputLocation();
+     _peds_t = peddata.GetPedsFrame();
+     _trajName = peddata.GetTrajName();
+     _projectRootDir = peddata.GetProjectRootDir();
+     _measureAreaId = boost::lexical_cast<string>(_areaForMethod_I->_id);
+     _fps =peddata.GetFps();
+     int mycounter = 0;
+     int minFrame = peddata.GetMinFrame();
+     Log->Write("INFO:\tMethod D: frame rate fps: <%.2f>, start: <%d>, stop: <%d> (minFrame = %d)", _fps, _startFrame, _stopFrame, minFrame);
+     if(_startFrame != _stopFrame)
+     {
+          if(_startFrame==-1)
+          {
+               _startFrame = minFrame;
+          }
+          if(_stopFrame==-1)
+          {
+               _stopFrame = peddata.GetNumFrames()+minFrame;
+          }
+          for(std::map<int , std::vector<int> >::iterator ite=_peds_t.begin();ite!=_peds_t.end();)
+          {
+               if((ite->first + minFrame)<_startFrame || (ite->first + minFrame) >_stopFrame)
+               {
+                    mycounter++;
+                    ite = _peds_t.erase(ite);
+               }
+               else
+               {
+                    ++ite;
+               }
+
+
+          }
+     }
+
+     if(!OpenFileMethodD())
+     {
+          return_value = false;
+     }
+     if(_calcIndividualFD)
+     {
+          if (!OpenFileIndividualFD())
+          {
+               return_value = false;
+          }
+     }
+     Log->Write("------------------------Analyzing with Method D-----------------------------");
+     //for(int frameNr = 0; frameNr < peddata.GetNumFrames(); frameNr++ )
+     //for(std::map<int , std::vector<int> >::iterator ite=_peds_t.begin();ite!=_peds_t.end();ite++)
+     for(auto ite: _peds_t)
+     {
+          int frameNr = ite.first;
+          int frid =  frameNr + minFrame;
+          //padd the frameid with 0
+          std::ostringstream ss;
+          ss << std::setw(5) << std::setfill('0') << std::internal << frid;
+          const std::string str_frid = ss.str();
+          if(!(frid%50))
+          {
+               Log->Write("INFO:\tframe ID = %d",frid);
+          }
+          vector<int> ids=_peds_t[frameNr];
+          vector<int> IdInFrame = peddata.GetIdInFrame(frameNr, ids, zPos_measureArea);
+          vector<double> XInFrame = peddata.GetXInFrame(frameNr, ids, zPos_measureArea);
+          vector<double> YInFrame = peddata.GetYInFrame(frameNr, ids, zPos_measureArea);
+          vector<double> VInFrame = peddata.GetVInFrame(frameNr, ids, zPos_measureArea);
+          //vector int to_remove
+          //------------------------------Remove peds outside geometry------------------------------------------
+          for( int i=0;i<(int)IdInFrame.size();i++)
+          {
+               if(false==within(point_2d(round(XInFrame[i]), round(YInFrame[i])), _geoPoly))
+               {
+                    Log->Write("Warning:\tPedestrian at <x=%.4f, y=%.4f> is not in geometry and not considered in analysis!", XInFrame[i]*CMtoM, YInFrame[i]*CMtoM );
+                    IdInFrame.erase(IdInFrame.begin() + i);
+                    XInFrame.erase(XInFrame.begin() + i);
+                    YInFrame.erase(YInFrame.begin() + i);
+                    VInFrame.erase(VInFrame.begin() + i);
+                    i--;
+               }
+          }
+          int NumPeds = IdInFrame.size();
+          //---------------------------------------------------------------------------------------------------------------
+          if(NumPeds>3)
+          {
+               if(_isOneDimensional)
+               {
+                    CalcVoronoiResults1D(XInFrame, VInFrame, IdInFrame, _areaForMethod_I->_poly,str_frid);
+               }
+               else
+               {
+                    if(IsPointsOnOneLine(XInFrame, YInFrame))
+                    {
+                         if(fabs(XInFrame[1]-XInFrame[0])<dmin)
+                         {
+                              XInFrame[1]+= offset;
+                         }
+                         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())
+                    {
+                         OutputVoronoiResults(polygons, str_frid, VInFrame); // TODO polygons_id
+                         if(_calcIndividualFD)
+                         {
+                              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);
+                         }
+                    }
+                    else
+                    {
+                         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);
+                    }
+               }
+          } // if N >3
+          else
+          {
+               Log->Write("INFO: \tThe number of the pedestrians is small (%d). Frame = %d (minFrame = %d)\n", NumPeds, frid, minFrame);
+          }
+     }
+     fclose(_fVoronoiRhoV);
+     if(_calcIndividualFD)
+     {
+          fclose(_fIndividualFD);
+     }
+     return return_value;
+}
+
+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());
+// _outputLocation.string() +  voroLocation+"rho_v_Voronoi_"+_trajName+"_id_"+_measureAreaId+".dat";
+     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)
+          {
+               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;
+     }
+}
+
+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
+     {
+          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;
+     }
+}
+
+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)
+     {
+          poly = p.first;
+          ReducePrecision(poly);
+          // TODO update polygon_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);
+}
+
+/*
+ * 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::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++;
+     }
+     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);
+          }
+          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;
+
+     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";
+     }
+
+     fs::path polygonPath=voroLocPath / "polygon";
+     if(!fs::exists(polygonPath))
+     {
+          if(!fs::create_directory(polygonPath))
+          {
+               Log->Write("ERROR:\tcan not create directory <%s>", polygonPath.string().c_str());
+               exit(EXIT_FAILURE);
+          }
+     }
+     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)
+          {
+               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;
+          }
+     }
+     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))
+          {
+               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())
+     {
+          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();
+}
+
+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_string2(polygon_iterator);
+              // string measureArea_str =
+              // polygon_to_string(measureArea); // maybe used for debugging
+              string v_str = polygon_to_string2(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::SetCalculateIndividualFD(bool individualFD)
+{
+     _calcIndividualFD = individualFD;
+}
+
+void Method_I::SetStartFrame(int startFrame)
+{
+     _startFrame=startFrame;
+}
+
+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::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::SetGeometryFileName(const fs::path& geometryFile)
+{
+     _geometryFileName=geometryFile;
+}
+
+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::SetCalculateProfiles(bool calcProfile)
+{
+     _getProfile =calcProfile;
+}
+
+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::SetMeasurementArea (MeasurementArea_B* area)
+{
+     _areaForMethod_I = area;
+}
+
+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);
+     }
+}
+
+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++)
+          {
+               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;
+}
+
+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)
+               {
+                    voronoi_distance.push_back(dist[j]);
+                    XLeftNeighbor.push_back(XLeftTemp[j]);
+                    XRightNeighbor.push_back(XRightTemp[j]);
+                    break;
+               }
+          }
+     }
+
+     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)
+          {
+               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 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)
+     {
+          for(unsigned int i=2; i<XInFrame.size();i++)
+          {
+               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++)
+          {
+               double dist=fabs(slope*XInFrame[i] - YInFrame[i] + intercept)/sqrt(slope*slope +1);
+               if(dist > dmin)
+               {
+                    isOnOneLine=false;
+                    break;
+               }
+          }
+     }
+     return isOnOneLine;
+}
diff --git a/methods/Method_I.h b/methods/Method_I.h
new file mode 100644
index 0000000000000000000000000000000000000000..26373caaecbc696a28c066b95099423d496cd944
--- /dev/null
+++ b/methods/Method_I.h
@@ -0,0 +1,119 @@
+/**
+ * \file        Method_I.h
+ * \date        Feb 07, 201
+ * \version     v0.8
+ * \copyright   <2009-2019> Forschungszentrum Juelich GmbH. All rights reserved.
+ *
+ * \section License
+ * This file is part of JuPedSim.
+ *
+ * JuPedSim is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * any later version.
+ *
+ * JuPedSim is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with JuPedSim. If not, see <http://www.gnu.org/licenses/>.
+ *
+ * \section Description
+ *
+ *
+ **/
+
+#ifndef Method_I_H_
+#define Method_I_H_
+
+#include "PedData.h"
+#include "../Analysis.h"
+#include "VoronoiDiagram.h"
+
+
+
+//handle more than two person are in one line
+#define dmin 200
+#define offset 200
+
+
+class Method_I
+{
+public:
+     Method_I();
+     virtual ~Method_I();
+     bool Process (const PedData& peddata,const fs::path& scriptsLocation, const double& zPos_measureArea);
+     void SetCalculateIndividualFD(bool individualFD);
+     void Setcutbycircle(double radius,int edges);
+     void SetGeometryPolygon(polygon_2d geometryPolygon);
+     void SetGeometryFileName(const fs::path& geometryFile);
+     void SetGeometryBoundaries(double minX, double minY, double maxX, double maxY);
+     void SetGridSize(double x, double y);
+     void SetCalculateProfiles(bool calcProfile);
+     void SetOutputVoronoiCellData(bool outputCellData);
+     void SetPlotVoronoiGraph(bool plotVoronoiGraph);
+     void SetPlotVoronoiIndex(bool plotVoronoiIndex);
+     void SetMeasurementArea (MeasurementArea_B* area);
+     void SetDimensional (bool dimension);
+     void SetTrajectoriesLocation(const fs::path& trajectoryPath);
+     void SetStartFrame(int startFrame);
+     void SetStopFrame(int stopFrame);
+
+private:
+     std::map<int , std::vector<int> > _peds_t;
+     std::string _measureAreaId;
+     MeasurementArea_B* _areaForMethod_I;
+     fs::path _trajName;
+     fs::path _projectRootDir;
+     fs::path _outputLocation;
+     fs::path _scriptsLocation;
+     bool _calcIndividualFD;
+     polygon_2d _areaIndividualFD;
+     bool _getProfile;
+     bool _outputVoronoiCellData;
+     bool _plotVoronoiCellData;
+     bool _plotVoronoiIndex;
+     bool _isOneDimensional;
+     bool _cutByCircle;       //Adjust whether cut each original voronoi cell by a circle
+     double _cutRadius;
+     int _circleEdges;
+     polygon_2d _geoPoly;
+     double _geoMinX;  // LOWest vertex of the geometry (x coordinate)
+     double _geoMinY;  //  LOWest vertex of the geometry (y coordinate)
+     double _geoMaxX; // Highest vertex of the geometry
+     double _geoMaxY;
+     FILE* _fVoronoiRhoV;
+     FILE* _fIndividualFD;
+     double _grid_size_X;      // the size of the grid
+     double _grid_size_Y;
+     float _fps;
+     bool OpenFileMethodD();
+     bool OpenFileIndividualFD();
+     fs::path _geometryFileName;
+     fs::path _trajectoryPath;
+     int _startFrame;
+     int _stopFrame;
+
+
+     std::vector<std::pair<polygon_2d, int> >  GetPolygons(std::vector<double>& XInFrame, std::vector<double>& YInFrame,
+                                                           std::vector<double>& VInFrame, std::vector<int>& IdInFrame);
+     void OutputVoronoiResults(const std::vector<polygon_2d>&  polygons, const std::string& frid, const std::vector<double>& VInFrame);
+     std::tuple<double,double> GetVoronoiDensityVelocity(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<std::pair<polygon_2d, int> >& polygons, int numPedsInFrame,std::vector<double>& XInFrame,
+                          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, const std::string& frid);
+     /**
+      * Reduce the precision of the points to two digits
+      * @param polygon
+      */
+     void CalcVoronoiResults1D(std::vector<double>& XInFrame, std::vector<double>& VInFrame, std::vector<int>& IdInFrame, const polygon_2d & measureArea, const std::string& frid);
+     void ReducePrecision(polygon_2d& polygon);
+     bool IsPedInGeometry(int frames, int peds, double **Xcor, double **Ycor, int  *firstFrame, int *lastFrame); //check whether all the pedestrians are in the geometry
+     double getOverlapRatio(const double& left, const double& right, const double& measurearea_left, const double& measurearea_right);
+     bool IsPointsOnOneLine(std::vector<double>& XInFrame, std::vector<double>& YInFrame);
+};
+
+#endif /* Method_I_H_ */