diff --git a/Analysis.cpp b/Analysis.cpp
index 4b2c821629dbc3d818cb604f9c05ffe1bb00cb1f..ff25423d7fd385d22e5e6e713323f1934d61363d 100644
--- a/Analysis.cpp
+++ b/Analysis.cpp
@@ -119,6 +119,11 @@ std::string Analysis::GetFilename (const std::string& str)
 void Analysis::InitArgs(ArgumentParser* args)
 {
      string s = "Parameter:\n";
+     _building = new Building();
+     _building->LoadGeometry(args->GetGeometryFilename());
+     // create the polygons
+     _building->InitGeometry();
+     // _building->AddSurroundingRoom();
 
      if(args->GetIsMethodA()) {
           _DoesUseMethodA = true;
@@ -177,6 +182,11 @@ void Analysis::InitArgs(ArgumentParser* args)
           _geoPoly = ReadGeometry(args->GetGeometryFilename(), _areaForMethod_I);
      }
 
+     if( _DoesUseMethodD &&  _DoesUseMethodI)
+     {
+          Log->Write("Warning:\t Using both method D and I is not safe!");
+          // because ReadGeomtry() may be called twice
+     }
      _deltaF = args->GetDelatT_Vins();
      _cutByCircle = args->GetIsCutByCircle();
      _getProfile = args->GetIsGetProfile();
@@ -200,12 +210,7 @@ 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();
-
+     Log->Write("INFO:\tReadGeometry with %s", geometryFile.string().c_str());
      double geo_minX  = FLT_MAX;
      double geo_minY  = FLT_MAX;
      double geo_maxX  = -FLT_MAX;
diff --git a/methods/Method_I.cpp b/methods/Method_I.cpp
index 4f498a1822fb9e7bc9f2520a111175841f38f075..ff96188f99e4fa1c3d300154d23c0ebdc6bc5cfd 100644
--- a/methods/Method_I.cpp
+++ b/methods/Method_I.cpp
@@ -136,6 +136,7 @@ bool Method_I::Process(const PedData& peddata,const fs::path& scriptsLocation, c
           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> ZInFrame = peddata.GetZInFrame(frameNr, ids, zPos_measureArea);
           vector<double> VInFrame = peddata.GetVInFrame(frameNr, ids, zPos_measureArea);
           //vector int to_remove
           //------------------------------Remove peds outside geometry------------------------------------------
@@ -143,10 +144,11 @@ bool Method_I::Process(const PedData& peddata,const fs::path& scriptsLocation, c
           {
                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 );
+                    Log->Write("Warning:\tPedestrian at <x=%.4f, y=%.4f, , z=%.4f> is not in geometry and not considered in analysis!", XInFrame[i]*CMtoM, YInFrame[i]*CMtoM, ZInFrame[i]*CMtoM );
                     IdInFrame.erase(IdInFrame.begin() + i);
                     XInFrame.erase(XInFrame.begin() + i);
                     YInFrame.erase(YInFrame.begin() + i);
+                    ZInFrame.erase(ZInFrame.begin() + i);
                     VInFrame.erase(VInFrame.begin() + i);
                     Log->Write("Warning:\t Pedestrian removed");
                     i--;
@@ -189,7 +191,7 @@ bool Method_I::Process(const PedData& peddata,const fs::path& scriptsLocation, c
                          if(!_isOneDimensional)
                          {
                               // GetIndividualFD(polygons,VInFrame, IdInFrame,  str_frid); // TODO polygons_id
-                              GetIndividualFD(polygons,VInFrame, IdInFrame,  str_frid, XInFrame, YInFrame); //
+                              GetIndividualFD(polygons,VInFrame, IdInFrame,  str_frid, XInFrame, YInFrame, ZInFrame); //
                          }
                     }
                     if(_getProfile)
@@ -266,7 +268,7 @@ bool Method_I::Process(const PedData& peddata,const fs::path& scriptsLocation, c
                }
                else
                {
-                    fprintf(_fIndividualFD,"#framerate (fps):\t%.2f\n\n#Frame	\t	PersID	\t	x/m \t y/m \t Individual density(m^(-2)) \t   Individual velocity(m/s)  \t Voronoi Polygon\n",_fps);
+                    fprintf(_fIndividualFD,"#framerate (fps):\t%.2f\n\n#Frame	\t	PersID	\t	x/m \t y/m \t z/m \t Individual density(m^(-2)) \t   Individual velocity(m/s)  \t Voronoi Polygon\n",_fps);
                }
                return true;
           }
@@ -547,11 +549,11 @@ bool Method_I::Process(const PedData& peddata,const fs::path& scriptsLocation, c
           }
      }
 
-void Method_I::GetIndividualFD(const vector<polygon_2d>& polygon, const vector<double>& Velocity, const vector<int>& Id, const string& frid, vector<double>& XInFrame, vector<double>& YInFrame)
+void Method_I::GetIndividualFD(const vector<polygon_2d>& polygon, const vector<double>& Velocity, const vector<int>& Id, const string& frid, vector<double>& XInFrame, vector<double>& YInFrame, vector<double>& ZInFrame)
 {
      double uniquedensity=0;
      double uniquevelocity=0;
-     double x, y;
+     double x, y, z;
      int uniqueId=0;
      int temp=0;
      for (const auto & polygon_iterator:polygon)
@@ -562,11 +564,13 @@ void Method_I::GetIndividualFD(const vector<polygon_2d>& polygon, const vector<d
           uniqueId=Id[temp];
           x = XInFrame[temp]*CMtoM;
           y = YInFrame[temp]*CMtoM;
-          fprintf(_fIndividualFD,"%s\t %d\t %.4f\t %.4f\t %.4f\t %.4f\t %s\n",
+          z = ZInFrame[temp]*CMtoM;
+          fprintf(_fIndividualFD,"%s\t %d\t %.4f\t %.4f\t %.4f\t %.4f\t %.4f\t %s\n",
                   frid.c_str(),
                   uniqueId,
                   x,
                   y,
+                  z,
                   uniquedensity,
                   uniquevelocity,
                   polygon_str.c_str()
diff --git a/methods/Method_I.h b/methods/Method_I.h
index 380f65460c1468c83b22ef0a1a895e160490f50f..01ec100843c7d3e951eede75814b31a83eaf87a0 100644
--- a/methods/Method_I.h
+++ b/methods/Method_I.h
@@ -104,7 +104,7 @@ private:
      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,const std::vector<double>& VInFrame);
      void GetIndividualFD(const std::vector<polygon_2d>& polygon, const std::vector<double>& Velocity, const std::vector<int>& Id, const std::string& frid);
-     void GetIndividualFD(const std::vector<polygon_2d>& polygon, const std::vector<double>& Velocity, const std::vector<int>& Id, const std::string& frid, std::vector<double>& XInFrame, std::vector<double>& YInFrame);
+     void GetIndividualFD(const std::vector<polygon_2d>& polygon, const std::vector<double>& Velocity, const std::vector<int>& Id, const std::string& frid, std::vector<double>& XInFrame, std::vector<double>& YInFrame, std::vector<double>& ZInFrame);
      /**
       * Reduce the precision of the points to two digits
       * @param polygon
diff --git a/methods/PedData.cpp b/methods/PedData.cpp
index 2e4123dfd742ceb0fdb4e0e522771424f44f45af..66b88d2ec228d2ed16837a561ca4b9976cf313b1 100644
--- a/methods/PedData.cpp
+++ b/methods/PedData.cpp
@@ -430,7 +430,7 @@ bool PedData::InitializeVariables(TiXmlElement* xRootNode)
           for(TiXmlElement* xAgent = xFrame->FirstChildElement("agent"); xAgent;
               xAgent = xAgent->NextSiblingElement("agent"))
           {
-               //get agent id, x, y
+               //get agent id, x, y, z
                double x= atof(xAgent->Attribute("x"));
                double y= atof(xAgent->Attribute("y"));
                double z= atof(xAgent->Attribute("z"));
@@ -587,6 +587,27 @@ vector<double> PedData::GetZInFrame(int frame, const vector<int>& ids) const
      }
      return ZInFrame;
 }
+vector<double> PedData::GetZInFrame(int frame, const vector<int>& ids, double zPos) const
+{
+     vector<double> ZInFrame;
+     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)
+               {
+                    ZInFrame.push_back(_zCor(id,frame));
+               }
+          }
+          else
+          {
+               ZInFrame.push_back(_zCor(id,frame));
+          }
+
+     }
+     return ZInFrame;
+}
 
 vector<int> PedData::GetIdInFrame(int frame, const vector<int>& ids) const
 {
diff --git a/methods/PedData.h b/methods/PedData.h
index d71576de0ec03c6ac76150b2074709dd9f1f642f..3fa72be24da4d96cff2bbb81fc3f4c7e00026711 100644
--- a/methods/PedData.h
+++ b/methods/PedData.h
@@ -77,6 +77,7 @@ public:
      std::vector<double> GetXInFrame(int frame, const std::vector<int>& ids) const;
      std::vector<double> GetYInFrame(int frame, const std::vector<int>& ids) const;
      std::vector<double> GetZInFrame(int frame, const std::vector<int>& ids) const;
+     std::vector<double> GetZInFrame(int frame, const std::vector<int>& ids, double zPos) const;
      std::vector<double> GetVInFrame(int frame, const std::vector<int>& ids, double zPos) const;
      bool ReadData(const fs::path& projectRootDir,const fs::path& outputDir, const fs::path& path, const fs::path& filename, const FileFormat& _trajformat, int deltaF, std::string vComponent, const bool IgnoreBackwardMovement);
      fs::path GetOutputLocation() const;