From 3c6763c909403af374c821a2de0b97e7fed78fd7 Mon Sep 17 00:00:00 2001
From: Otmane Lahlou <otmane.lahlou@c-s.fr>
Date: Mon, 19 Mar 2012 17:28:26 +0100
Subject: [PATCH] ENH: use OGR to compute distance between geometries

---
 Code/Common/otbDataNode.h   |  11 +++-
 Code/Common/otbDataNode.txx | 104 ++++++++++++++++++++++++++++++++++++
 2 files changed, 114 insertions(+), 1 deletion(-)

diff --git a/Code/Common/otbDataNode.h b/Code/Common/otbDataNode.h
index 11985fb767..4efa6242f5 100644
--- a/Code/Common/otbDataNode.h
+++ b/Code/Common/otbDataNode.h
@@ -26,6 +26,8 @@
 
 #include <iostream>
 
+#include "ogrsf_frmts.h"
+
 namespace otb
 {
 /**
@@ -91,6 +93,8 @@ public:
   typedef typename LineType::Pointer                                           LinePointerType;
   typedef typename LineType::ConstPointer                                      LineConstPointerType;
   typedef Polygon<ValuePrecisionType>                                          PolygonType;
+  typedef typename PolygonType::VertexListType                                 VertexListType;
+  typedef typename VertexListType::ConstPointer                                VertexListConstPointerType;
   typedef typename PolygonType::Pointer                                        PolygonPointerType;
   typedef typename PolygonType::ConstPointer                                   PolygonConstPointerType;
   typedef ObjectList<PolygonType>                                              PolygonListType;
@@ -258,6 +262,10 @@ public:
     */
   std::vector<std::string> GetFieldList() const;
 
+  /** \return the distance to a point */
+  double EuclideanDistance(const DataNode* node);
+
+  double EuclideanDistance(const PointType point);
 /**
  * Clear all fields.
  */
@@ -271,6 +279,8 @@ protected:
   /** PrintSelf method */
   void PrintSelf(std::ostream& os, itk::Indent indent) const;
 
+  OGRGeometry* ConvertDataNodeToOGRGeometry(const DataNode* dataNode);
+
 private:
   DataNode(const Self&); //purposely not implemented
   void operator =(const Self&); //purposely not implemented
@@ -296,7 +306,6 @@ private:
 
   /** The fields map */
 //   FieldMapType         m_FieldMap;
-
 };
 } // end namespace
 
diff --git a/Code/Common/otbDataNode.txx b/Code/Common/otbDataNode.txx
index d020ec7ce2..d7878caf66 100644
--- a/Code/Common/otbDataNode.txx
+++ b/Code/Common/otbDataNode.txx
@@ -532,6 +532,110 @@ DataNode<TPrecision, VDimension, TValuePrecision>
 {
   return m_NodeType == FEATURE_COLLECTION;
 }
+
+
+template <class TPrecision, unsigned int VDimension, class TValuePrecision>
+OGRGeometry *
+DataNode<TPrecision, VDimension, TValuePrecision>
+::ConvertDataNodeToOGRGeometry(const DataNode* dataNode)
+{
+  switch(dataNode->GetNodeType())
+    {
+    case FEATURE_POINT:
+    {
+    OGRPoint *ogrPoint = (OGRPoint *) OGRGeometryFactory::createGeometry(wkbPoint); ;
+    ogrPoint->setX(dataNode->GetPoint()[0]);
+    ogrPoint->setY(dataNode->GetPoint()[1]);
+    return ogrPoint;
+    }
+    break;
+    case FEATURE_LINE:
+    {
+    //Build the ogrObject
+    OGRLineString *             ogrLine = (OGRLineString *) OGRGeometryFactory::createGeometry(wkbLineString);
+    VertexListConstPointerType vertexList = dataNode->GetLine()->GetVertexList();
+    typename VertexListType::ConstIterator vIt = vertexList->Begin();
+
+    while (vIt != vertexList->End())
+      {
+      OGRPoint ogrPoint;
+      ogrPoint.setX(vIt.Value()[0]);
+      ogrPoint.setY(vIt.Value()[1]);
+      if (Dimension > 2)
+        {
+        ogrPoint.setZ(vIt.Value()[2]);
+        }
+      ogrLine->addPoint(&ogrPoint);
+      ++vIt;
+      }
+    return ogrLine;
+    }
+    break;
+    case FEATURE_POLYGON:
+    {
+    OGRPolygon * polygon = (OGRPolygon *) OGRGeometryFactory::createGeometry(wkbPolygon);
+    OGRLinearRing *  ogrExternalRing = (OGRLinearRing *) OGRGeometryFactory::createGeometry(wkbLinearRing);
+
+    VertexListConstPointerType vertexList = dataNode->GetPolygonExteriorRing()->GetVertexList();
+
+    typename VertexListType::ConstIterator vIt = vertexList->Begin();
+
+    while (vIt != vertexList->End())
+      {
+      OGRPoint ogrPoint;
+      ogrPoint.setX(vIt.Value()[0]);
+      ogrPoint.setY(vIt.Value()[1]);
+      if ( Dimension > 2)
+        {
+        ogrPoint.setZ(vIt.Value()[2]);
+        }
+
+      ogrExternalRing->addPoint(&ogrPoint);
+      ++vIt;
+      }
+    polygon->addRing(ogrExternalRing);
+    // Close the polygon
+    polygon->closeRings();
+    OGRGeometryFactory::destroyGeometry(ogrExternalRing);
+    return polygon;
+    }
+    break;
+    }
+
+  return NULL;  
+}
+
+template <class TPrecision, unsigned int VDimension, class TValuePrecision>
+double
+DataNode<TPrecision, VDimension, TValuePrecision>
+::EuclideanDistance(const DataNode* node)
+{
+  // Convert the nodes to OGRGeometries
+  OGRGeometry * dstGeomtery     = this->ConvertDataNodeToOGRGeometry(node);
+  OGRGeometry * currentGeometry = this->ConvertDataNodeToOGRGeometry(this);
+  
+  // Compute the distance
+  return currentGeometry->Distance(dstGeomtery);
+}
+
+template <class TPrecision, unsigned int VDimension, class TValuePrecision>
+double
+DataNode<TPrecision, VDimension, TValuePrecision>
+::EuclideanDistance(const PointType point)
+{
+  // Convert Point to point to ogrPoint
+  OGRPoint ogrPointSrc;
+  ogrPointSrc.setX(point[0]);
+  ogrPointSrc.setY(point[1]);
+
+  // Convert the current datanode to an OGRGeometry
+  OGRGeometry * currentGeometry = this->ConvertDataNodeToOGRGeometry(this);
+
+  // return the distance
+  return currentGeometry->Distance(&ogrPointSrc);
+}
+
+
 } // end namespace otb
 
 #endif
-- 
GitLab