From: ageay <ageay>
Date: Tue, 5 Jun 2012 15:42:37 +0000 (+0000)
Subject: 2D Convex hull computation.
X-Git-Tag: V6_main_FINAL~669
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=d8b0fb1760354306e67d091ca617fbe4c0870618;p=tools%2Fmedcoupling.git

2D Convex hull computation.
---

diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx
index ae6a410d3..ec57ec0b9 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -3052,6 +3052,56 @@ void MEDCouplingUMesh::checkButterflyCells(std::vector<int>& cells, double eps)
     }
 }
 
+/*!
+ * This method is typically requested to unbutterfly 2D linear cells in \b this.
+ *
+ * This method expects that space dimension is equal to 2 and mesh dimension is equal to 2 too. If it is not the case an INTERP_KERNEL::Exception will be thrown.
+ * This method works only for linear 2D cells. If there is any of non linear cells (INTERP_KERNEL::NORM_QUAD8 for example) an INTERP_KERNEL::Exception will be thrown too.
+ * 
+ * For each 2D linear cell in \b this, this method builds the convex envelop (or the convex hull) of the current cell.
+ * This convex envelop is computed using Jarvis march algorithm.
+ * The coordinates and the number of cells of \b this remain unchanged on invocation of this method.
+ * Only connectivity of some cells could be modified if those cells were not representing a convex envelop. If a cell already equals its convex envelop (regardless orientation)
+ * its connectivity will remain unchanged. If the computation leads to a modification of nodal connectivity of a cell its geometric type will be modified to INTERP_KERNEL::NORM_POLYGON.
+ *
+ * @return a newly allocated array containing cellIds that have been modified if any. If no cells have been impacted by this method NULL is returned.
+ */
+DataArrayInt *MEDCouplingUMesh::convexEnvelop2D() throw(INTERP_KERNEL::Exception)
+{
+  if(getMeshDimension()!=2 || getSpaceDimension()!=2)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convexEnvelop2D  works only for meshDim=2 and spaceDim=2 !");
+  checkFullyDefined();
+  const double *coords=getCoords()->getConstPointer();
+  int nbOfCells=getNumberOfCells();
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nodalConnecIndexOut=DataArrayInt::New();
+  nodalConnecIndexOut->alloc(nbOfCells+1,1);
+  std::vector<int> nodalConnecOut;
+  int *workIndexOut=nodalConnecIndexOut->getPointer();
+  *workIndexOut=0;
+  const int *nodalConnecIn=_nodal_connec->getConstPointer();
+  const int *nodalConnecIndexIn=_nodal_connec_index->getConstPointer();
+  std::set<INTERP_KERNEL::NormalizedCellType> types;
+  std::vector<int> isChanged;
+  for(int i=0;i<nbOfCells;i++,workIndexOut++)
+    {
+      std::size_t pos=nodalConnecOut.size()+1;
+      if(BuildConvecEnvelopOf2DCellJarvis(coords,nodalConnecIn+nodalConnecIndexIn[i],nodalConnecIn+nodalConnecIndexIn[i+1],nodalConnecOut))
+        isChanged.push_back(i);
+      types.insert((INTERP_KERNEL::NormalizedCellType)nodalConnecOut[pos]);
+      workIndexOut[1]=(int)nodalConnecOut.size();
+    }
+  if(isChanged.empty())
+    return 0;
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nodalConnecOut2=DataArrayInt::New();
+  nodalConnecOut2->alloc((int)nodalConnecOut.size(),1);
+  std::copy(nodalConnecOut.begin(),nodalConnecOut.end(),nodalConnecOut2->getPointer());
+  setConnectivity(nodalConnecOut2,nodalConnecIndexOut,false);
+  _types=types;
+  DataArrayInt *ret=DataArrayInt::New(); ret->alloc((int)isChanged.size(),1);
+  std::copy(isChanged.begin(),isChanged.end(),ret->getPointer());
+  return ret;
+}
+
 /*!
  * This method is expected to be applied on a mesh with spaceDim==3 and meshDim==3. If not an exception will be thrown.
  * This method analyzes only linear extruded 3D cells (NORM_HEXA8,NORM_PENTA6,NORM_HEXGP12...)
@@ -5878,6 +5928,101 @@ void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair<i
     }
 }
 
+/*!
+ * This method compute the convex hull of a single 2D cell. This method tries to conserve at maximum the given input connectivity. In particular, if the orientation of cell is not clockwise
+ * as in MED format norm. If definitely the result of Jarvis algorithm is not matchable with the input connectivity, the result will be copied into \b nodalConnecOut parameter and
+ * the geometric cell type set to INTERP_KERNEL::NORM_POLYGON.
+ * This method excepts that \b coords parameter is expected to be in dimension 2. [\b nodalConnBg, \b nodalConnEnd) is the nodal connectivity of the input
+ * cell (geometric cell type included at the position 0). If the meshdimension of the input cell is not equal to 2 an INTERP_KERNEL::Exception will be thrown.
+ * 
+ * @return false if the input connectivity represents already the convex hull, true if the input cell needs to be reordered.
+ */
+bool MEDCouplingUMesh::BuildConvecEnvelopOf2DCellJarvis(const double *coords, const int *nodalConnBg, const int *nodalConnEnd, std::vector<int>& nodalConnecOut) throw(INTERP_KERNEL::Exception)
+{
+  std::size_t sz=std::distance(nodalConnBg,nodalConnEnd);
+  if(sz>=4)
+    {
+      const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)*nodalConnBg);
+      if(cm.getDimension()==2)
+        {
+          const int *node=nodalConnBg+1;
+          int startNode=*node++;
+          double refX=coords[2*startNode];
+          for(;node!=nodalConnEnd;node++)
+            {
+              if(coords[2*(*node)]<refX)
+                {
+                  startNode=*node;
+                  refX=coords[2*startNode];
+                }
+            }
+          std::vector<int> tmpOut; tmpOut.reserve(sz); tmpOut.push_back(startNode);
+          refX=1e300;
+          double tmp1;
+          double tmp2[2];
+          double angle0=-M_PI/2;
+          //
+          int nextNode=-1;
+          int prevNode=-1;
+          double resRef;
+          double angleNext;
+          while(nextNode!=startNode)
+            {
+              nextNode=-1;
+              resRef=1e300;
+              for(node=nodalConnBg+1;node!=nodalConnEnd;node++)
+                {
+                  if(*node!=tmpOut.back() && *node!=prevNode)
+                    {
+                      tmp2[0]=coords[2*(*node)]-coords[2*tmpOut.back()]; tmp2[1]=coords[2*(*node)+1]-coords[2*tmpOut.back()+1];
+                      double angleM=INTERP_KERNEL::EdgeArcCircle::GetAbsoluteAngle(tmp2,tmp1);
+                      double res;
+                      if(angleM<=angle0)
+                        res=angle0-angleM;
+                      else
+                        res=angle0-angleM+2.*M_PI;
+                      if(res<resRef)
+                        {
+                          nextNode=*node;
+                          resRef=res;
+                          angleNext=angleM;
+                        }
+                    }
+                }
+              if(nextNode!=startNode)
+                {
+                  angle0=angleNext-M_PI;
+                  if(angle0<-M_PI)
+                    angle0+=2*M_PI;
+                  prevNode=tmpOut.back();
+                  tmpOut.push_back(nextNode);
+                }
+            }
+          std::vector<int> tmp3(2*(sz-1));
+          std::vector<int>::iterator it=std::copy(nodalConnBg+1,nodalConnEnd,tmp3.begin());
+          std::copy(nodalConnBg+1,nodalConnEnd,it);
+          if(std::search(tmp3.begin(),tmp3.end(),tmpOut.begin(),tmpOut.end())!=tmp3.end())
+            {
+              nodalConnecOut.insert(nodalConnecOut.end(),nodalConnBg,nodalConnEnd);
+              return false;
+            }
+          if(std::search(tmp3.rbegin(),tmp3.rend(),tmpOut.begin(),tmpOut.end())!=tmp3.rend())
+            {
+              nodalConnecOut.insert(nodalConnecOut.end(),nodalConnBg,nodalConnEnd);
+              return false;
+            }
+          else
+            {
+              nodalConnecOut.push_back((int)INTERP_KERNEL::NORM_POLYGON);
+              nodalConnecOut.insert(nodalConnecOut.end(),tmpOut.begin(),tmpOut.end());
+              return true;
+            }
+        }
+    }
+  else
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::BuildConvecEnvelopOf2DCellJarvis : invalid 2D cell connectivity !");
+}
+
 /*!
  * Given a 2D mesh conn by (conn2D,connI2D) it returns a single polygon
  */
diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx
index 163a8ccf9..95982a9f9 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.hxx
@@ -138,6 +138,7 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT void getCellsContainingPoint(const double *pos, double eps, std::vector<int>& elts) const;
     MEDCOUPLING_EXPORT void getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, std::vector<int>& elts, std::vector<int>& eltsIndex) const;
     MEDCOUPLING_EXPORT void checkButterflyCells(std::vector<int>& cells, double eps=1e-12) const;
+    MEDCOUPLING_EXPORT DataArrayInt *convexEnvelop2D() throw(INTERP_KERNEL::Exception);
     MEDCOUPLING_EXPORT void findAndCorrectBadOriented3DExtrudedCells(std::vector<int>& cells) throw(INTERP_KERNEL::Exception);
     MEDCOUPLING_EXPORT void getBoundingBoxForBBTree(std::vector<double>& bbox) const;
     MEDCOUPLING_EXPORT MEDCouplingUMesh *buildExtrudedMesh(const MEDCouplingUMesh *mesh1D, int policy);
@@ -190,9 +191,8 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT static bool IsPolygonWellOriented(bool isQuadratic, const double *vec, const int *begin, const int *end, const double *coords);
     MEDCOUPLING_EXPORT static bool IsPolyhedronWellOriented(const int *begin, const int *end, const double *coords);
     MEDCOUPLING_EXPORT static void TryToCorrectPolyhedronOrientation(int *begin, int *end, const double *coords) throw(INTERP_KERNEL::Exception);
-/// @cond INTERNAL
     MEDCOUPLING_EXPORT static MEDCouplingUMesh *Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2) throw(INTERP_KERNEL::Exception);
-/// @endcond
+    MEDCOUPLING_EXPORT static bool BuildConvecEnvelopOf2DCellJarvis(const double *coords, const int *nodalConnBg, const int *nodalConnEnd, std::vector<int>& nodalConnecOut) throw(INTERP_KERNEL::Exception);
   private:
     MEDCouplingUMesh();
     MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCopy);
diff --git a/src/MEDCoupling_Swig/MEDCoupling.i b/src/MEDCoupling_Swig/MEDCoupling.i
index ed9a8de46..96dde2318 100644
--- a/src/MEDCoupling_Swig/MEDCoupling.i
+++ b/src/MEDCoupling_Swig/MEDCoupling.i
@@ -272,6 +272,7 @@ using namespace INTERP_KERNEL;
 %newobject ParaMEDMEM::MEDCouplingUMesh::getCellIdsLyingOnNodes;
 %newobject ParaMEDMEM::MEDCouplingUMesh::buildSetInstanceFromThis;
 %newobject ParaMEDMEM::MEDCouplingUMesh::getCellIdsCrossingPlane;
+%newobject ParaMEDMEM::MEDCouplingUMesh::convexEnvelop2D;
 %newobject ParaMEDMEM::MEDCouplingUMeshCellByTypeEntry::__iter__;
 %newobject ParaMEDMEM::MEDCouplingUMeshCellEntry::__iter__;
 %newobject ParaMEDMEM::MEDCouplingExtrudedMesh::New;
@@ -1081,6 +1082,7 @@ namespace ParaMEDMEM
     MEDCouplingFieldDouble *getAspectRatioField() const throw(INTERP_KERNEL::Exception);
     MEDCouplingFieldDouble *getWarpField() const throw(INTERP_KERNEL::Exception);
     MEDCouplingFieldDouble *getSkewField() const throw(INTERP_KERNEL::Exception);
+    DataArrayInt *convexEnvelop2D() throw(INTERP_KERNEL::Exception);
     static MEDCouplingUMesh *Build0DMeshFromCoords(DataArrayDouble *da) throw(INTERP_KERNEL::Exception);
     static MEDCouplingUMesh *MergeUMeshes(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception);
     static MEDCouplingUMesh *MergeUMeshesOnSameCoords(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception);