diff --git a/methods/Method_D.cpp b/methods/Method_D.cpp
index 784e1dc690e9a7c93f2a5a8ee522f70214897e5a..55127f0de139113deba2e03d68f4a028286b489f 100644
--- a/methods/Method_D.cpp
+++ b/methods/Method_D.cpp
@@ -28,7 +28,9 @@
 
 #include "Method_D.h"
 #include <cmath>
-
+#include<map>
+#include<iostream>
+#include<vector>
 //using std::string;
 //using std::vector;
 //using std::ofstream;
@@ -135,7 +137,7 @@ bool Method_D::Process (const PedData& peddata,const std::string& scriptsLocatio
           {
                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> 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);
@@ -164,24 +166,28 @@ bool Method_D::Process (const PedData& peddata,const std::string& scriptsLocatio
                               YInFrame[1]+= offset;
                          }
                     }
-                    vector<polygon_2d> polygons = GetPolygons(XInFrame, YInFrame, VInFrame, IdInFrame);
+                    std::vector<std::pair<polygon_2d, int> > polygons_id = GetPolygons(XInFrame, YInFrame, VInFrame, IdInFrame);
+                    vector<polygon_2d> polygons;
+                    for (auto p: polygons_id)
+                         polygons.push_back(p.first);
+
                     if(!polygons.empty())
                     {
-                         OutputVoronoiResults(polygons, str_frid, VInFrame);
+                         OutputVoronoiResults(polygons, str_frid, VInFrame); // TODO polygons_id
                          if(_calcIndividualFD)
                          {
                               if(false==_isOneDimensional)
                               {
-                                   GetIndividualFD(polygons,VInFrame, IdInFrame, _areaForMethod_D->_poly, str_frid);
+                                   GetIndividualFD(polygons,VInFrame, IdInFrame, _areaForMethod_D->_poly, str_frid); // TODO polygons_id
                               }
                          }
                          if(_getProfile)
                          { //	field analysis
-                              GetProfiles(str_frid, polygons, VInFrame);
+                              GetProfiles(str_frid, polygons, VInFrame); // TODO polygons_id
                          }
                          if(_outputVoronoiCellData)
                          { // output the Voronoi polygons of a frame
-                              OutputVoroGraph(str_frid, polygons, NumPeds, XInFrame, YInFrame,VInFrame);
+                              OutputVoroGraph(str_frid, polygons_id, NumPeds, XInFrame, YInFrame,VInFrame);
                          }
                     }
                     else
@@ -193,7 +199,7 @@ bool Method_D::Process (const PedData& peddata,const std::string& scriptsLocatio
           }
           else
           {
-               //Log->Write("INFO: \tThe number of the pedestrians is less than 2 !!");
+               Log->Write("INFO: \tThe number of the pedestrians is less than 2 !!");
           }
      }
      fclose(_fVoronoiRhoV);
@@ -248,23 +254,27 @@ bool Method_D::OpenFileIndividualFD()
      }
 }
 
-vector<polygon_2d> Method_D::GetPolygons(vector<double>& XInFrame, vector<double>& YInFrame, vector<double>& VInFrame, vector<int>& IdInFrame)
+std::vector<std::pair<polygon_2d, int> > Method_D::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)));
-     vector<polygon_2d>  polygons = vd.getVoronoiPolygons(XInFrame, YInFrame, VInFrame,IdInFrame, boundpoint);
+     std::vector<std::pair<polygon_2d, int> > polygons_id;
+     polygons_id = vd.getVoronoiPolygons(XInFrame, YInFrame, VInFrame,IdInFrame, boundpoint);
+     polygon_2d poly ;
      if(_cutByCircle)
      {
-          polygons = vd.cutPolygonsWithCircle(polygons, XInFrame, YInFrame, _cutRadius,_circleEdges);
+          polygons_id = vd.cutPolygonsWithCircle(polygons_id, XInFrame, YInFrame, _cutRadius,_circleEdges);
      }
-     polygons = vd.cutPolygonsWithGeometry(polygons, _geoPoly, XInFrame, YInFrame);
+     polygons_id = vd.cutPolygonsWithGeometry(polygons_id, _geoPoly, XInFrame, YInFrame);
 
-     for(auto && p:polygons)
+     for(auto && p:polygons_id)
      {
-          ReducePrecision(p);
+          poly = p.first;
+          ReducePrecision(poly);
+          // TODO update polygon_id?
      }
-     return polygons;
+     return polygons_id;
 }
 /**
  * Output the Voronoi density and velocity in the corresponding file
@@ -395,10 +405,11 @@ void Method_D::GetProfiles(const string& frameId, const vector<polygon_2d>& poly
      fclose(Prf_density);
 }
 
-void Method_D::OutputVoroGraph(const string & frameId, vector<polygon_2d>& polygons, int numPedsInFrame, vector<double>& XInFrame, vector<double>& YInFrame,const vector<double>& VInFrame)
+void Method_D::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;
      string voronoiLocation=_projectRootDir+VORO_LOCATION+"VoronoiCell/";
+     polygon_2d poly;
 
 
 #if defined(_WIN32)
@@ -407,14 +418,14 @@ void Method_D::OutputVoroGraph(const string & frameId, vector<polygon_2d>& polyg
      mkdir(voronoiLocation.c_str(), 0777);
 #endif
 
-
      string polygon=voronoiLocation+"/polygon"+_trajName+"_id_"+_measureAreaId+"_"+frameId+".dat";
      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 && poly:polygons)
+          for(auto && p:polygons_id)
           {
+               poly = p.first;
                for(auto&& point:poly.outer())
                {
                     point.x(point.x()*CMtoM);
@@ -428,7 +439,8 @@ void Method_D::OutputVoroGraph(const string & frameId, vector<polygon_2d>& polyg
                          point.y(point.y()*CMtoM);
                     }
                }
-               polys << dsv(poly) << endl;
+               polys << "["<< p.second << ", " <<dsv(poly) << "]" << endl;
+               //polys  <<dsv(poly)<< endl;
           }
      }
      else
diff --git a/methods/Method_D.h b/methods/Method_D.h
index a115b780b80c2ee01d780c3fae5ede9493ac0a2f..21be1cabc62c9aa252399257e6c14da595e8068a 100644
--- a/methods/Method_D.h
+++ b/methods/Method_D.h
@@ -103,14 +103,14 @@ private:
      int _stopFrame;
 
 
-     std::vector<polygon_2d> GetPolygons(std::vector<double>& XInFrame, std::vector<double>& YInFrame,
+     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);
      double GetVoronoiDensity(const std::vector<polygon_2d>& polygon, const polygon_2d & measureArea);
      double GetVoronoiDensity2(const std::vector<polygon_2d>& polygon, double* XInFrame, double* YInFrame, const polygon_2d& measureArea);
      double GetVoronoiVelocity(const std::vector<polygon_2d>& polygon, const std::vector<double>& Velocity, const polygon_2d & measureArea);
      void GetProfiles(const std::string& frameId, const std::vector<polygon_2d>& polygons, const std::vector<double>& velocity);
-     void OutputVoroGraph(const std::string & frameId, std::vector<polygon_2d>& polygons, int numPedsInFrame,std::vector<double>& XInFrame,
+     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);
      /**
diff --git a/methods/VoronoiDiagram.cpp b/methods/VoronoiDiagram.cpp
index 89ca23fdd94ef1b0b752cfea46f84328a0e1d457..4662efbb0b5d9bf637457eb6a489de1f51e34ad3 100644
--- a/methods/VoronoiDiagram.cpp
+++ b/methods/VoronoiDiagram.cpp
@@ -44,9 +44,9 @@ VoronoiDiagram::VoronoiDiagram()
 VoronoiDiagram::~VoronoiDiagram()
 {
 }
-
+//typedef std::vector<std::pair<polygon_2d>, int> >  poly_id;
 // Traversing Voronoi edges using cell iterator.
-vector<polygon_2d> VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame, vector<double>& YInFrame,
+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();
@@ -59,16 +59,16 @@ vector<polygon_2d> VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame,
      {
           points.push_back(point_type2(round(XInFrame[i]), round(YInFrame[i])));
           XInFrame_temp.push_back(round(XInFrame[i]));
-		  YInFrame_temp.push_back(round(YInFrame[i]));
-		  VInFrame_temp.push_back(VInFrame[i]);
-		  IdInFrame_temp.push_back(IdInFrame[i]);
+          YInFrame_temp.push_back(round(YInFrame[i]));
+          VInFrame_temp.push_back(VInFrame[i]);
+          IdInFrame_temp.push_back(IdInFrame[i]);
      }
 
      VD voronoidiagram;
      construct_voronoi(points.begin(), points.end(), &voronoidiagram);
      int Ncell = 0;
-     std::vector<polygon_2d> polygons;
-
+     //std::vector<polygon_2d> polygons;
+     std::vector<std::pair<polygon_2d, int> >    polygons_id;
      double Bd_Box_minX = -Bound_Max;
      double Bd_Box_minY = -Bound_Max;
      double Bd_Box_maxX = Bound_Max;
@@ -83,7 +83,7 @@ vector<polygon_2d> VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame,
      for (voronoi_diagram<double>::const_cell_iterator it = voronoidiagram.cells().begin();
                it != voronoidiagram.cells().end(); ++it)
      {
-    	  polygon_2d poly;
+          polygon_2d poly;
           vector<point_type2> polypts;
           point_type2 pt_s;
           point_type2 pt_e;
@@ -170,7 +170,6 @@ vector<polygon_2d> VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame,
                }
                edge = edge->next();
           } while (edge != cell.incident_edge());
-          Ncell++;
 
           if (infinite_s && infinite_e)
           {
@@ -190,9 +189,14 @@ vector<polygon_2d> VoronoiDiagram::getVoronoiPolygons(vector<double>& XInFrame,
           }
 
           boost::geometry::correct(poly);
-          polygons.push_back(poly);
+          //polygons.push_back(poly);
+          //cout << "poly is: " << typeid(poly).name() << '\n'
+          int id_ped = IdInFrame[Ncell];
+          std::pair<polygon_2d, int>  poly_id = std::make_pair(poly, id_ped);
+          polygons_id.push_back(poly_id);
+          Ncell++;
      }
-     return polygons;
+     return polygons_id;
 }
 
 point_type2 VoronoiDiagram::retrieve_point(const cell_type& cell)
@@ -359,7 +363,7 @@ vector<point_type2> VoronoiDiagram::add_bounding_points(const point_type2& pt1,
                else
                {
                     double tempx = pt1.x()
-                            		                      + (pt2.x() - pt1.x()) * (pt.y() - pt1.y()) / (pt2.y() - pt1.y());
+                                                              + (pt2.x() - pt1.x()) * (pt.y() - pt1.y()) / (pt2.y() - pt1.y());
                     if (tempx < pt.x())
                     {
                          bounding_vertex.push_back(point_type2(minX, maxY));
@@ -391,7 +395,7 @@ vector<point_type2> VoronoiDiagram::add_bounding_points(const point_type2& pt1,
                else
                {
                     double tempy = pt1.y()
-                            		                      + (pt2.y() - pt1.y()) * (pt.x() - pt1.x()) / (pt2.x() - pt1.x());
+                                                              + (pt2.y() - pt1.y()) * (pt.x() - pt1.x()) / (pt2.x() - pt1.x());
                     if (tempy < pt.y())
                     {
                          bounding_vertex.push_back(point_type2(minX, minY));
@@ -433,15 +437,18 @@ point_type2 VoronoiDiagram::getIntersectionPoint(const point_2d& pt0, const poin
      return interpts;
 }
 
-std::vector<polygon_2d> VoronoiDiagram::cutPolygonsWithGeometry(const std::vector<polygon_2d>& polygon,
+std::vector<std::pair<polygon_2d, int> > VoronoiDiagram::cutPolygonsWithGeometry(const std::vector<std::pair<polygon_2d, int> > & polygon,
           const polygon_2d& Geometry, const vector<double>& xs, const vector<double>& ys)
 {
-     vector<polygon_2d> intersetionpolygons;
+     // vector<polygon_2d> intersetionpolygons;
+     std::vector<std::pair<polygon_2d, int> > intersetionpolygons;
      int temp = 0;
      for (const auto & polygon_iterator:polygon)
      {
           polygon_list v;
-          intersection(Geometry, polygon_iterator, v);
+          polygon_2d p = polygon_iterator.first;
+
+          intersection(Geometry, p, v);
           BOOST_FOREACH(auto & it, v)
           {
                if (within(point_2d(xs[temp], ys[temp]), it))
@@ -452,7 +459,7 @@ std::vector<polygon_2d> VoronoiDiagram::cutPolygonsWithGeometry(const std::vecto
                     simplify(it,simplified,J_EPS);
 
                     correct(simplified);
-                    intersetionpolygons.push_back(simplified);
+                    intersetionpolygons.push_back(std::make_pair(simplified, polygon_iterator.second));
                }
           }
           temp++;
@@ -461,14 +468,15 @@ std::vector<polygon_2d> VoronoiDiagram::cutPolygonsWithGeometry(const std::vecto
 }
 
 
-std::vector<polygon_2d> VoronoiDiagram::cutPolygonsWithCircle(const std::vector<polygon_2d>& polygon,
+std::vector<std::pair<polygon_2d, int> >  VoronoiDiagram::cutPolygonsWithCircle(const std::vector<std::pair<polygon_2d, int> > & polygon,
           const vector<double>& xs, const vector<double>& ys, double radius, int edges)
 {
-     std::vector<polygon_2d> intersetionpolygons;
+     std::vector<std::pair<polygon_2d, int> >  intersetionpolygons;
      int temp = 0;
      for (const auto & polygon_iterator:polygon)
      {
           polygon_2d circle;
+          polygon_2d p = polygon_iterator.first;
           {
                for (int angle = 0; angle <=edges; angle++)
                {
@@ -479,7 +487,7 @@ std::vector<polygon_2d> VoronoiDiagram::cutPolygonsWithCircle(const std::vector<
           }
           correct(circle);
           polygon_list v;
-          intersection(circle, polygon_iterator, v);
+          intersection(circle, p, v);
           BOOST_FOREACH(auto & it, v)
           {
                if (within(point_2d(xs[temp], ys[temp]), it))
@@ -489,7 +497,7 @@ std::vector<polygon_2d> VoronoiDiagram::cutPolygonsWithCircle(const std::vector<
                     polygon_2d simplified;
                     simplify(it,simplified,J_EPS);
                     correct(simplified);
-                    intersetionpolygons.push_back(simplified);
+                    intersetionpolygons.push_back(std::make_pair(simplified, polygon_iterator.second));
                }
           }
           temp++;
diff --git a/methods/VoronoiDiagram.h b/methods/VoronoiDiagram.h
index d8f0c58d1f8a172b6d3391c42ab3e98313b4b6f4..88ec2091e796b1ac6d9f5eb07fe6a90ce2521869 100644
--- a/methods/VoronoiDiagram.h
+++ b/methods/VoronoiDiagram.h
@@ -79,10 +79,10 @@ public:
      VoronoiDiagram();
      virtual ~VoronoiDiagram();
 
-     std::vector<polygon_2d> getVoronoiPolygons(std::vector<double>& XInFrame, std::vector<double>& YInFrame, std::vector<double>& VInFrame,
+    std::vector<std::pair<polygon_2d, int> > getVoronoiPolygons(std::vector<double>& XInFrame, std::vector<double>& YInFrame, std::vector<double>& VInFrame,
                std::vector<int>& IdInFrame, const double Bound_Max);
-     std::vector<polygon_2d> cutPolygonsWithGeometry(const std::vector<polygon_2d>& polygon, const polygon_2d& Geometry, const std::vector<double>& xs, const std::vector<double>& ys);
-     std::vector<polygon_2d> cutPolygonsWithCircle(const std::vector<polygon_2d>& polygon, const std::vector<double>& xs, const std::vector<double>& ys, double radius, int edges);
+     std::vector<std::pair<polygon_2d, int> >  cutPolygonsWithGeometry(const std::vector<std::pair<polygon_2d, int> > & polygon, const polygon_2d& Geometry, const std::vector<double>& xs, const std::vector<double>& ys);
+     std::vector<std::pair<polygon_2d, int> >  cutPolygonsWithCircle(const std::vector<std::pair<polygon_2d, int> > & polygon, const std::vector<double>& xs, const std::vector<double>& ys, double radius, int edges);
 };
 
-#endif /* VORONOIDIAGRAM_H_ */
\ No newline at end of file
+#endif /* VORONOIDIAGRAM_H_ */
diff --git a/scripts/_Plot_cell_rho.py b/scripts/_Plot_cell_rho.py
index e469951241c91291d09a0fdb1bd64e32cc1703e6..d2edf9f7e1bb2e931df46dc2fcba57ccb8fe0673 100644
--- a/scripts/_Plot_cell_rho.py
+++ b/scripts/_Plot_cell_rho.py
@@ -2,7 +2,7 @@ import os
 import sys
 import logging
 import argparse
-#from xml.dom import minidom
+from xml.dom import minidom
 try:
     import xml.etree.cElementTree as ET
 except ImportError:
@@ -117,7 +117,7 @@ def parse_xml_traj_file(filename):
         xmldoc = minidom.parse(filename)
     except:
         logging.critical('could not parse file. exit')
-        exit(FAILURE)
+        exit(0)
     N = int(xmldoc.getElementsByTagName('agents')[0].childNodes[0].data)
     fps = xmldoc.getElementsByTagName('frameRate')[0].childNodes[0].data #type unicode
     fps = round(float(fps))