[QGIS Commit] r11500 - in trunk/qgis: images/themes/default python/core src/app src/core

svn_qgis at osgeo.org svn_qgis at osgeo.org
Tue Aug 25 15:54:28 EDT 2009


Author: mhugent
Date: 2009-08-25 15:54:28 -0400 (Tue, 25 Aug 2009)
New Revision: 11500

Added:
   trunk/qgis/images/themes/default/mActionReshape.png
   trunk/qgis/src/app/qgsmaptoolreshape.cpp
   trunk/qgis/src/app/qgsmaptoolreshape.h
Modified:
   trunk/qgis/python/core/qgsgeometry.sip
   trunk/qgis/src/app/CMakeLists.txt
   trunk/qgis/src/app/qgisapp.cpp
   trunk/qgis/src/app/qgisapp.h
   trunk/qgis/src/core/qgsgeometry.cpp
   trunk/qgis/src/core/qgsgeometry.h
Log:
[FEATURE] A reshape tool to apply to line/polygon geometries. The part of a geometry between the first and last intersection of the reshape line will be replaced

Added: trunk/qgis/images/themes/default/mActionReshape.png
===================================================================
(Binary files differ)


Property changes on: trunk/qgis/images/themes/default/mActionReshape.png
___________________________________________________________________
Added: svn:mime-type
   + application/octet-stream

Modified: trunk/qgis/python/core/qgsgeometry.sip
===================================================================
--- trunk/qgis/python/core/qgsgeometry.sip	2009-08-25 15:09:16 UTC (rev 11499)
+++ trunk/qgis/python/core/qgsgeometry.sip	2009-08-25 19:54:28 UTC (rev 11500)
@@ -204,6 +204,11 @@
     @return 0 in case of success, 1 if geometry has not been split, error else*/
     int splitGeometry(const QList<QgsPoint>& splitLine, QList<QgsGeometry*>& newGeometries, bool topological, QList<QgsPoint>& topologyTestPoints);
 
+/**Replaces a part of this geometry with another line
+      @return 0 in case of success
+      @note: this function was added in version 1.3*/
+    int reshapeGeometry( const QList<QgsPoint>& reshapeWithLine );
+
  /**Changes this geometry such that it does not intersect the other geometry
        @param other geometry that should not be intersect
        @return 0 in case of success*/

Modified: trunk/qgis/src/app/CMakeLists.txt
===================================================================
--- trunk/qgis/src/app/CMakeLists.txt	2009-08-25 15:09:16 UTC (rev 11499)
+++ trunk/qgis/src/app/CMakeLists.txt	2009-08-25 19:54:28 UTC (rev 11500)
@@ -39,6 +39,7 @@
   qgsmaptoolmovefeature.cpp
   qgsmaptoolmovevertex.cpp
   qgsmaptoolnodetool.cpp
+  qgsmaptoolreshape.cpp
   qgsmaptoolselect.cpp
   qgsmaptoolsimplify.cpp
   qgsmaptoolsplitfeatures.cpp

Modified: trunk/qgis/src/app/qgisapp.cpp
===================================================================
--- trunk/qgis/src/app/qgisapp.cpp	2009-08-25 15:09:16 UTC (rev 11499)
+++ trunk/qgis/src/app/qgisapp.cpp	2009-08-25 19:54:28 UTC (rev 11500)
@@ -171,6 +171,7 @@
 #include "qgsmaptoolnodetool.h"
 #include "qgsmaptoolpan.h"
 #include "qgsmaptoolselect.h"
+#include "qgsmaptoolreshape.h"
 #include "qgsmaptoolsplitfeatures.h"
 #include "qgsmaptoolvertexedit.h"
 #include "qgsmaptoolzoom.h"
@@ -473,6 +474,7 @@
   delete mMapTools.mCaptureLine;
   delete mMapTools.mCapturePolygon;
   delete mMapTools.mMoveFeature;
+  delete mMapTools.mReshapeFeatures;
   delete mMapTools.mSplitFeatures;
   delete mMapTools.mSelect;
   delete mMapTools.mVertexAdd;
@@ -660,6 +662,12 @@
   connect( mActionMoveFeature, SIGNAL( triggered() ), this, SLOT( moveFeature() ) );
   mActionMoveFeature->setEnabled( false );
 
+  mActionReshapeFeatures = new QAction( getThemeIcon( "mActionReshape.png" ), tr( "Reshape Features" ), this );
+  shortcuts->registerAction( mActionReshapeFeatures );
+  mActionReshapeFeatures->setStatusTip( tr( "Reshape Features" ) );
+  connect( mActionReshapeFeatures, SIGNAL( triggered() ), this, SLOT( reshapeFeatures() ) );
+  mActionReshapeFeatures->setEnabled( false );
+
   mActionSplitFeatures = new QAction( getThemeIcon( "mActionSplitFeatures.png" ), tr( "Split Features" ), this );
   shortcuts->registerAction( mActionSplitFeatures );
   mActionSplitFeatures->setStatusTip( tr( "Split Features" ) );
@@ -1062,6 +1070,8 @@
   mMapToolGroup->addAction( mActionCapturePolygon );
   mActionMoveFeature->setCheckable( true );
   mMapToolGroup->addAction( mActionMoveFeature );
+  mActionReshapeFeatures->setCheckable( true );
+  mMapToolGroup->addAction( mActionReshapeFeatures );
   mActionSplitFeatures->setCheckable( true );
   mMapToolGroup->addAction( mActionSplitFeatures );
   mMapToolGroup->addAction( mActionDeleteSelected );
@@ -1175,6 +1185,7 @@
   mEditMenu->addAction( mActionAddIsland );
   mEditMenu->addAction( mActionDeleteRing );
   mEditMenu->addAction( mActionDeletePart );
+  mEditMenu->addAction( mActionReshapeFeatures );
   mEditMenu->addAction( mActionSplitFeatures );
   mEditMenu->addAction( mActionMergeFeatures );
   mEditMenu->addAction( mActionNodeTool );
@@ -1383,6 +1394,7 @@
   mAdvancedDigitizeToolBar->addAction( mActionAddIsland );
   mAdvancedDigitizeToolBar->addAction( mActionDeleteRing );
   mAdvancedDigitizeToolBar->addAction( mActionDeletePart );
+  mAdvancedDigitizeToolBar->addAction( mActionReshapeFeatures );
   mAdvancedDigitizeToolBar->addAction( mActionSplitFeatures );
   mAdvancedDigitizeToolBar->addAction( mActionMergeFeatures );
   mAdvancedDigitizeToolBar->addAction( mActionNodeTool );
@@ -1605,6 +1617,7 @@
   mActionCaptureLine->setIcon( getThemeIcon( "/mActionCaptureLine.png" ) );
   mActionCapturePolygon->setIcon( getThemeIcon( "/mActionCapturePolygon.png" ) );
   mActionMoveFeature->setIcon( getThemeIcon( "/mActionMoveFeature.png" ) );
+  mActionReshapeFeatures->setIcon( getThemeIcon( "/mActionReshape.png" ) );
   mActionSplitFeatures->setIcon( getThemeIcon( "/mActionSplitFeatures.png" ) );
   mActionDeleteSelected->setIcon( getThemeIcon( "/mActionDeleteSelected.png" ) );
   mActionAddVertex->setIcon( getThemeIcon( "/mActionAddVertex.png" ) );
@@ -1735,6 +1748,8 @@
   mActionCapturePolygon->setVisible( false );
   mMapTools.mMoveFeature = new QgsMapToolMoveFeature( mMapCanvas );
   mMapTools.mMoveFeature->setAction( mActionMoveFeature );
+  mMapTools.mReshapeFeatures = new QgsMapToolReshape( mMapCanvas );
+  mMapTools.mReshapeFeatures->setAction( mActionReshapeFeatures );
   mMapTools.mSplitFeatures = new QgsMapToolSplitFeatures( mMapCanvas );
   mMapTools.mSplitFeatures->setAction( mActionSplitFeatures );
   mMapTools.mSelect = new QgsMapToolSelect( mMapCanvas );
@@ -4337,6 +4352,11 @@
   mMapCanvas->setMapTool( mMapTools.mSplitFeatures );
 }
 
+void QgisApp::reshapeFeatures()
+{
+  mMapCanvas->setMapTool( mMapTools.mReshapeFeatures );
+}
+
 void QgisApp::capturePoint()
 {
   if ( mMapCanvas && mMapCanvas->isDrawing() )
@@ -5634,6 +5654,7 @@
         mActionMoveVertex->setEnabled( false );
         mActionAddRing->setEnabled( false );
         mActionAddIsland->setEnabled( false );
+        mActionReshapeFeatures->setEnabled( false );
         mActionSplitFeatures->setEnabled( false );
         mActionSimplifyFeature->setEnabled( false );
         mActionDeleteRing->setEnabled( false );
@@ -5650,6 +5671,7 @@
         {
           mActionCaptureLine->setEnabled( true );
           mActionCaptureLine->setVisible( true );
+          mActionReshapeFeatures->setEnabled( true );
           mActionSplitFeatures->setEnabled( true );
           mActionSimplifyFeature->setEnabled( true );
           mActionDeletePart->setEnabled( true );
@@ -5659,6 +5681,7 @@
         {
           mActionCaptureLine->setEnabled( false );
           mActionCaptureLine->setVisible( false );
+          mActionReshapeFeatures->setEnabled( false );
           mActionSplitFeatures->setEnabled( false );
           mActionSimplifyFeature->setEnabled( false );
           mActionDeletePart->setEnabled( false );
@@ -5679,6 +5702,7 @@
           mActionCapturePolygon->setVisible( true );
           mActionAddRing->setEnabled( true );
           mActionAddIsland->setEnabled( true );
+          mActionReshapeFeatures->setEnabled( true );
           mActionSplitFeatures->setEnabled( true );
           mActionSimplifyFeature->setEnabled( true );
           mActionDeleteRing->setEnabled( true );
@@ -5690,6 +5714,7 @@
           mActionCapturePolygon->setVisible( false );
           mActionAddRing->setEnabled( false );
           mActionAddIsland->setEnabled( false );
+          mActionReshapeFeatures->setEnabled( false );
           mActionSplitFeatures->setEnabled( false );
           mActionSimplifyFeature->setEnabled( false );
           mActionDeleteRing->setEnabled( false );

Modified: trunk/qgis/src/app/qgisapp.h
===================================================================
--- trunk/qgis/src/app/qgisapp.h	2009-08-25 15:09:16 UTC (rev 11499)
+++ trunk/qgis/src/app/qgisapp.h	2009-08-25 19:54:28 UTC (rev 11500)
@@ -497,6 +497,8 @@
     void deleteSelected();
     //! activates the move feature tool
     void moveFeature();
+    //! activates the reshape features tool
+    void reshapeFeatures();
     //! activates the split features tool
     void splitFeatures();
     //! activates the add vertex tool
@@ -725,6 +727,7 @@
     QAction *mActionCapturePolygon;
     QAction *mActionDeleteSelected;
     QAction *mActionMoveFeature;
+    QAction *mActionReshapeFeatures;
     QAction *mActionSplitFeatures;
     QAction *mActionAddVertex;
     QAction *mActionDeleteVertex;
@@ -847,6 +850,7 @@
         QgsMapTool* mCaptureLine;
         QgsMapTool* mCapturePolygon;
         QgsMapTool* mMoveFeature;
+        QgsMapTool* mReshapeFeatures;
         QgsMapTool* mSplitFeatures;
         QgsMapTool* mSelect;
         QgsMapTool* mVertexAdd;

Added: trunk/qgis/src/app/qgsmaptoolreshape.cpp
===================================================================
--- trunk/qgis/src/app/qgsmaptoolreshape.cpp	                        (rev 0)
+++ trunk/qgis/src/app/qgsmaptoolreshape.cpp	2009-08-25 19:54:28 UTC (rev 11500)
@@ -0,0 +1,128 @@
+/***************************************************************************
+    qgsmaptoolreshape.cpp
+    ---------------------------
+    begin                : Juli 2009
+    copyright            : (C) 2009 by Marco Hugentobler
+    email                : marco dot hugentobler at karto dot baug dot ethz dot ch
+ ***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+/* $Id$ */
+
+#include "qgsmaptoolreshape.h"
+#include "qgsgeometry.h"
+#include "qgsmapcanvas.h"
+#include "qgsrubberband.h"
+#include "qgsvectorlayer.h"
+#include <QMessageBox>
+#include <QMouseEvent>
+
+QgsMapToolReshape::QgsMapToolReshape( QgsMapCanvas* canvas ): QgsMapToolCapture( canvas, QgsMapToolCapture::CaptureLine )
+{
+
+}
+
+QgsMapToolReshape::~QgsMapToolReshape()
+{
+
+}
+
+void QgsMapToolReshape::canvasReleaseEvent( QMouseEvent * e )
+{
+  //check if we operate on a vector layer //todo: move this to a function in parent class to avoid duplication
+  QgsVectorLayer *vlayer = dynamic_cast <QgsVectorLayer*>( mCanvas->currentLayer() );
+
+  if ( !vlayer )
+  {
+    QMessageBox::information( 0, tr( "Not a vector layer" ),
+                              tr( "The current layer is not a vector layer" ) );
+    return;
+  }
+
+  if ( !vlayer->isEditable() )
+  {
+    QMessageBox::information( 0, tr( "Layer not editable" ),
+                              tr( "Cannot edit the vector layer. To make it editable, go to the file item "
+                                  "of the layer, right click and check 'Allow Editing'." ) );
+    return;
+  }
+
+  //add point to list and to rubber band
+  int error = addVertex( e->pos() );
+  if ( error == 1 )
+  {
+    //current layer is not a vector layer
+    return;
+  }
+  else if ( error == 2 )
+  {
+    //problem with coordinate transformation
+    QMessageBox::information( 0, tr( "Coordinate transform error" ),
+                              tr( "Cannot transform the point to the layers coordinate system" ) );
+    return;
+  }
+
+  if ( e->button() == Qt::LeftButton )
+  {
+    mCapturing = TRUE;
+  }
+  else if ( e->button() == Qt::RightButton )
+  {
+    mCapturing = FALSE;
+    delete mRubberBand;
+    mRubberBand = 0;
+
+    //find out bounding box of mCaptureList
+    if(mCaptureList.size() < 1)
+    {
+      return;
+    }
+    QgsPoint firstPoint = mCaptureList.at(0);
+    QgsRectangle bbox(firstPoint.x(), firstPoint.y(), firstPoint.x(), firstPoint.y());
+    for(int i = 1; i < mCaptureList.size(); ++i)
+    {
+      bbox.combineExtentWith(mCaptureList.at(i).x(), mCaptureList.at(i).y());
+    }
+
+    //query all the features that intersect bounding box of capture line
+    vlayer->select(QgsAttributeList(), bbox, true, false);
+    QgsFeature f;
+    int reshapeReturn;
+    bool reshapeDone = false;
+
+    vlayer->beginEditCommand( tr( "Reshape" ) );
+    while(vlayer->nextFeature(f))
+    {
+      //query geometry
+      //call geometry->reshape(mCaptureList)
+      //register changed geometry in vector layer
+      QgsGeometry* geom = f.geometry();
+      if(geom)
+      {
+        reshapeReturn = geom->reshapeGeometry(mCaptureList);
+        if(reshapeReturn == 0)
+        {
+          vlayer->changeGeometry(f.id(), geom);
+          reshapeDone = true;
+        }
+      }
+    }
+
+    if(reshapeDone)
+    {
+      vlayer->endEditCommand();
+    }
+    else
+    {
+      vlayer->destroyEditCommand();
+    }
+
+    mCaptureList.clear();
+    mCanvas->refresh();
+  }
+}

Added: trunk/qgis/src/app/qgsmaptoolreshape.h
===================================================================
--- trunk/qgis/src/app/qgsmaptoolreshape.h	                        (rev 0)
+++ trunk/qgis/src/app/qgsmaptoolreshape.h	2009-08-25 19:54:28 UTC (rev 11500)
@@ -0,0 +1,31 @@
+/***************************************************************************
+    qgsmaptoolreshape.h
+    ---------------------
+    begin                : Juli 2009
+    copyright            : (C) 2009 by Marco Hugentobler
+    email                : marco.hugentobler at karto dot baug dot ethz dot ch
+ ***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+/* $Id$ */
+
+#ifndef QGSMAPTOOLRESHAPE_H
+#define QGSMAPTOOLRESHAPE_H
+
+#include "qgsmaptoolcapture.h"
+
+/**A map tool that draws a line and splits the features cut by the line*/
+class QgsMapToolReshape: public QgsMapToolCapture
+{
+  public:
+    QgsMapToolReshape( QgsMapCanvas* canvas );
+    virtual ~QgsMapToolReshape();
+    void canvasReleaseEvent( QMouseEvent * e );
+};
+
+#endif

Modified: trunk/qgis/src/core/qgsgeometry.cpp
===================================================================
--- trunk/qgis/src/core/qgsgeometry.cpp	2009-08-25 15:09:16 UTC (rev 11499)
+++ trunk/qgis/src/core/qgsgeometry.cpp	2009-08-25 19:54:28 UTC (rev 11500)
@@ -17,6 +17,7 @@
 #include <limits>
 #include <cstdarg>
 #include <cstdio>
+#include <cmath>
 
 #include "qgis.h"
 #include "qgsgeometry.h"
@@ -3219,6 +3220,113 @@
   return returnCode;
 }
 
+/**Replaces a part of this geometry with another line*/
+int QgsGeometry::reshapeGeometry( const QList<QgsPoint>& reshapeWithLine )
+{
+  int returnCode = 0;
+  if ( type() == QGis::Point )
+  {
+    return 1; //cannot reshape points
+  }
+
+  GEOSGeometry* reshapeLineGeos = createGeosLineString( reshapeWithLine.toVector() );
+
+   //make sure this geos geometry is up-to-date
+  if ( !mGeos || mDirtyGeos )
+  {
+    exportWkbToGeos();
+  }
+
+  //single or multi?
+  int numGeoms = GEOSGetNumGeometries(mGeos);
+  if(numGeoms == -1)
+  {
+    return 1;
+  }
+
+  bool isMultiGeom = (numGeoms > 1);
+  bool isLine = ( type() == QGis::Line);
+
+  //polygon or multipolygon?
+  if ( !isMultiGeom )
+  {
+    GEOSGeometry* reshapedGeometry;
+    if(isLine)
+    {
+      reshapedGeometry = reshapeLine( mGeos, reshapeLineGeos );
+    }
+    else
+    {
+      reshapedGeometry = reshapePolygon( mGeos, reshapeLineGeos );
+    }
+
+    GEOSGeom_destroy(reshapeLineGeos);
+    if(reshapedGeometry)
+    {
+      GEOSGeom_destroy( mGeos );
+      mGeos = reshapedGeometry;
+      mDirtyWkb = true;
+      return 0;
+    }
+    else
+    {
+      return 1;
+    }
+  }
+  else
+  {
+    //call reshape for each geometry part and replace mGeos with new geometry if reshape took place
+    bool reshapeTookPlace = false;
+
+    GEOSGeometry* currentReshapeGeometry = 0;
+    GEOSGeometry** newGeoms = new GEOSGeometry*[numGeoms];
+
+    for(int i = 0; i < numGeoms; ++i)
+    {
+      if(isLine)
+      {
+        currentReshapeGeometry = reshapeLine( GEOSGetGeometryN(mGeos, i), reshapeLineGeos);
+      }
+      else
+      {
+        currentReshapeGeometry = reshapePolygon( GEOSGetGeometryN(mGeos, i), reshapeLineGeos);
+      }
+
+      if(currentReshapeGeometry)
+      {
+        newGeoms[i] = currentReshapeGeometry;
+        reshapeTookPlace = true;
+      }
+      else
+      {
+        newGeoms[i] = GEOSGeom_clone( GEOSGetGeometryN(mGeos, i) );
+      }
+    }
+     GEOSGeom_destroy(reshapeLineGeos);
+
+    GEOSGeometry* newMultiGeom = GEOSGeom_createCollection(GEOS_MULTIPOLYGON, newGeoms, numGeoms);
+    delete[] newGeoms;
+    if( ! newMultiGeom )
+    {
+      return 3;
+    }
+
+    if(reshapeTookPlace)
+    {
+      GEOSGeom_destroy( mGeos );
+      mGeos = newMultiGeom;
+      mDirtyWkb = true;
+      return 0;
+    }
+    else
+    {
+      GEOSGeom_destroy(newMultiGeom);
+      return 1;
+    }
+  }
+  return 1;
+}
+
 int QgsGeometry::makeDifference( QgsGeometry* other )
 {
   //make sure geos geometry is up to date
@@ -4801,28 +4909,13 @@
   QVector<GEOSGeometry*> testedGeometries;
   GEOSGeometry* intersectGeom = 0;
 
-  //Create a small buffer around the original geometry
-  //and intersect candidate line segments with the buffer.
-  //Then we use the ratio intersection length / segment length to
-  //decide if the line segment belongs to the original geometry or
-  //if it is part of the splitting line
-  double bufferDistance = 0.0000001;
-
   for ( int i = 0; i < GEOSGetNumGeometries( mergedLines ); i++ )
   {
     const GEOSGeometry *testing = GEOSGetGeometryN( mergedLines, i );
-    intersectGeom = GEOSIntersection( mGeos, GEOSBuffer( testing, bufferDistance, DEFAULT_QUADRANT_SEGMENTS ) );
-    double len;
-    GEOSLength( intersectGeom, &len );
-    double testingLen;
-    GEOSLength( testing, &testingLen );
-    double ratio = len / testingLen;
-    //the ratios for geometries that belong to the original line are usually close to 1
-    if ( ratio >= 0.5 && ratio <= 1.5 )
+    if ( lineContainedInLine( testing, mGeos ) == 1 )
     {
       testedGeometries << GEOSGeom_clone( testing );
     }
-    GEOSGeom_destroy( intersectGeom );
   }
 
   mergeGeometriesMultiTypeSplit( testedGeometries );
@@ -4962,6 +5055,331 @@
   return 0;
 }
 
+GEOSGeometry* QgsGeometry::reshapePolygon( const GEOSGeometry* polygon, const GEOSGeometry* reshapeLineGeos )
+{
+  //go through outer shell and all inner rings and check if there is exactly one intersection of a ring and the reshape line
+    int nIntersections = 0;
+    int lastIntersectingRing = -2;
+    const GEOSGeometry* lastIntersectingGeom = 0;
+
+    int nRings = GEOSGetNumInteriorRings(polygon);
+    if(nRings < 0)
+    {
+      return 0;
+    }
+
+    //does outer ring intersect?
+    const GEOSGeometry* outerRing = GEOSGetExteriorRing(polygon);
+    if( GEOSIntersects( outerRing, reshapeLineGeos ) == 1)
+    {
+      ++nIntersections;
+      lastIntersectingRing = -1;
+      lastIntersectingGeom = outerRing;
+    }
+
+    //do inner rings intersect?
+    const GEOSGeometry* innerRings[nRings];
+    for(int i = 0; i < nRings; ++i)
+    {
+      innerRings[i] = GEOSGetInteriorRingN(polygon, i);
+      if(GEOSIntersects( innerRings[i], reshapeLineGeos) == 1)
+      {
+        ++nIntersections;
+        lastIntersectingRing = i;
+        lastIntersectingGeom = innerRings[i];
+      }
+    }
+
+    if(nIntersections != 1) //reshape line is only allowed to intersect one ring
+    {
+      return 0;
+    }
+
+    //we have one intersecting ring, let's try to reshape it
+    GEOSGeometry* reshapeResult = reshapeLine( lastIntersectingGeom, reshapeLineGeos);
+    if(!reshapeResult)
+    {
+      return 0;
+    }
+
+    //if reshaping took place, we need to reassemble the polygon and its rings
+    GEOSGeometry* newRing = 0;
+    const GEOSCoordSequence* reshapeSequence = GEOSGeom_getCoordSeq(reshapeResult);
+    GEOSCoordSequence* newCoordSequence = GEOSCoordSeq_clone(reshapeSequence);
+
+    GEOSGeom_destroy(reshapeResult);
+
+    newRing = GEOSGeom_createLinearRing(newCoordSequence);
+    if(!newRing)
+    {
+      return 0;
+    }
+
+
+    GEOSGeometry* newOuterRing = 0;
+    if(lastIntersectingRing == -1)
+    {
+      newOuterRing = newRing;
+    }
+    else
+    {
+      newOuterRing = GEOSGeom_clone(outerRing);
+    }
+
+    GEOSGeometry** newInnerRings = new GEOSGeometry*[nRings];
+    for(int i = 0; i < nRings; ++i)
+    {
+      if( lastIntersectingRing == i)
+      {
+        newInnerRings[i] = newRing;
+      }
+      else
+      {
+        newInnerRings[i] = GEOSGeom_clone(innerRings[i]);
+      }
+    }
+
+    GEOSGeometry* reshapedPolygon = GEOSGeom_createPolygon( newOuterRing, newInnerRings, nRings);
+    delete[] newInnerRings;
+    if(!reshapedPolygon)
+    {
+      return 0;
+    }
+    return reshapedPolygon;
+}
+
+GEOSGeometry* QgsGeometry::reshapeLine( const GEOSGeometry* line, const GEOSGeometry* reshapeLineGeos )
+{
+  if ( !line || !reshapeLineGeos )
+  {
+    return 0;
+  }
+
+  //make sure there are at least two instersction between line and reshape geometry
+  GEOSGeometry* intersectGeom = GEOSIntersection( line, reshapeLineGeos);
+  bool atLeastTwoIntersections = (GEOSGeomTypeId(intersectGeom) == GEOS_MULTIPOINT && GEOSGetNumGeometries(intersectGeom) > 1);
+  GEOSGeom_destroy(intersectGeom);
+  if(!atLeastTwoIntersections)
+  {
+    return 0;
+  }
+
+
+  bool isRing = false;
+  if(GEOSGeomTypeId(line) == GEOS_LINEARRING)
+  {
+    isRing = true;
+  }
+
+  //begin and end point of original line
+  const GEOSCoordSequence* lineCoordSeq = GEOSGeom_getCoordSeq( line );
+  if ( !lineCoordSeq )
+  {
+    return 0;
+  }
+  unsigned int lineCoordSeqSize;
+  if ( GEOS_DLL GEOSCoordSeq_getSize( lineCoordSeq, &lineCoordSeqSize ) == 0 )
+  {
+    return 0;
+  }
+  if ( lineCoordSeqSize < 2 )
+  {
+    return 0;
+  }
+  //first and last vertex of line
+  double x1, y1, x2, y2;
+  GEOS_DLL GEOSCoordSeq_getX( lineCoordSeq, 0, &x1 );
+  GEOS_DLL GEOSCoordSeq_getY( lineCoordSeq, 0, &y1 );
+  GEOS_DLL GEOSCoordSeq_getX( lineCoordSeq, lineCoordSeqSize - 1, &x2 );
+  GEOS_DLL GEOSCoordSeq_getY( lineCoordSeq, lineCoordSeqSize - 1, &y2 );
+  GEOSGeometry* beginLineVertex = createGeosPoint( QgsPoint( x1, y1 ) );
+  GEOSGeometry* endLineVertex = createGeosPoint( QgsPoint( x2, y2 ) );
+
+//node line and reshape line
+  GEOSGeometry* nodedGeometry = nodeGeometries( reshapeLineGeos, line );
+  if ( !nodedGeometry )
+  {
+    GEOSGeom_destroy( beginLineVertex );
+    GEOSGeom_destroy( endLineVertex );
+    return 0;
+  }
+
+  //and merge them together
+  GEOSGeometry *mergedLines = GEOSLineMerge( nodedGeometry );
+  GEOSGeom_destroy( nodedGeometry );
+  if ( !mergedLines )
+  {
+    GEOSGeom_destroy( beginLineVertex );
+    GEOSGeom_destroy( endLineVertex );
+    return 0;
+  }
+
+  int numMergedLines = GEOSGetNumGeometries( mergedLines );
+  if ( numMergedLines < 2 ) //some special cases. Normally it is >2
+  {
+    GEOSGeom_destroy( beginLineVertex );
+    GEOSGeom_destroy( endLineVertex );
+    if ( numMergedLines == 1 ) //reshape line is from begin to endpoint. So we keep the reshapeline
+    {
+      return GEOSGeom_clone( reshapeLineGeos );
+    }
+    else
+    {
+      return 0;
+    }
+  }
+
+  QList<GEOSGeometry*> resultLineParts; //collection with the line segments that will be contained in result
+  QList<GEOSGeometry*> probableParts; //parts where we can decide on inclusion only after going through all the candidates
+
+  for ( int i = 0; i < numMergedLines; ++i )
+  {
+    const GEOSGeometry* currentGeom;
+
+    currentGeom = GEOSGetGeometryN( mergedLines, i );
+    const GEOSCoordSequence* currentCoordSeq = GEOSGeom_getCoordSeq( currentGeom );
+    unsigned int currentCoordSeqSize;
+    GEOS_DLL GEOSCoordSeq_getSize( currentCoordSeq, &currentCoordSeqSize );
+    if ( currentCoordSeqSize < 2 )
+    {
+      continue;
+    }
+
+    //get the two endpoints of the current line merge result
+    double xBegin, xEnd, yBegin, yEnd;
+    GEOS_DLL GEOSCoordSeq_getX( currentCoordSeq, 0, &xBegin );
+    GEOS_DLL GEOSCoordSeq_getY( currentCoordSeq, 0, &yBegin );
+    GEOS_DLL GEOSCoordSeq_getX( currentCoordSeq, currentCoordSeqSize - 1, &xEnd );
+    GEOS_DLL GEOSCoordSeq_getY( currentCoordSeq, currentCoordSeqSize - 1, &yEnd );
+    GEOSGeometry* beginCurrentGeomVertex = createGeosPoint( QgsPoint( xBegin, yBegin ) );
+    GEOSGeometry* endCurrentGeomVertex = createGeosPoint( QgsPoint( xEnd, yEnd ) );
+
+    //check how many endpoints of the line merge result are on the (original) line
+    int nEndpointsOnOriginalLine = 0;
+    if ( pointContainedInLine( beginCurrentGeomVertex, line ) == 1 )
+    {
+      nEndpointsOnOriginalLine += 1;
+    }
+
+    if ( pointContainedInLine( endCurrentGeomVertex, line ) == 1 )
+    {
+      nEndpointsOnOriginalLine += 1;
+    }
+
+    //check how many endpoints equal the endpoints of the original line
+    int nEndpointsSameAsOriginalLine = 0;
+    if ( GEOSEquals( beginCurrentGeomVertex, beginLineVertex ) == 1 || GEOSEquals( beginCurrentGeomVertex, endLineVertex ) == 1 )
+    {
+      nEndpointsSameAsOriginalLine += 1;
+    }
+    if ( GEOSEquals( endCurrentGeomVertex, beginLineVertex ) == 1 || GEOSEquals( endCurrentGeomVertex, endLineVertex ) == 1 )
+    {
+      nEndpointsSameAsOriginalLine += 1;
+    }
+
+    //check if the current geometry overlaps the original geometry (GEOSOverlap does not seem to work with linestrings)
+    bool currentGeomOverlapsOriginalGeom = false;
+    bool currentGeomOverlapsReshapeLine = false;
+    if ( QgsGeometry::lineContainedInLine( currentGeom, line ) == 1 )
+    {
+      currentGeomOverlapsOriginalGeom = true;
+    }
+    if ( QgsGeometry::lineContainedInLine( currentGeom, reshapeLineGeos ) == 1 )
+    {
+      currentGeomOverlapsReshapeLine = true;
+    }
+
+
+    //logic to decide if this part belongs to the result
+    if ( nEndpointsSameAsOriginalLine == 1 && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
+    {
+      resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
+    }
+    //for closed rings, we take one segment from the candidate list
+    else if(isRing && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom)
+    {
+      probableParts.push_back(GEOSGeom_clone(currentGeom));
+    }
+    else if ( nEndpointsOnOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
+    {
+      resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
+    }
+    else if ( nEndpointsSameAsOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
+    {
+      resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
+    }
+    else if ( currentGeomOverlapsOriginalGeom && currentGeomOverlapsReshapeLine )
+    {
+      resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
+    }
+
+    GEOSGeom_destroy( beginCurrentGeomVertex );
+    GEOSGeom_destroy( endCurrentGeomVertex );
+  }
+
+  //add the longest segment from the probable list for rings (only used for polygon rings)
+  if(isRing)
+  {
+    GEOSGeometry* maxGeom = 0; //the longest geometry in the probabla list
+    GEOSGeometry* currentGeom = 0;
+    double maxLength = -DBL_MAX;
+    double currentLength = 0;
+    for(int i = 0; i < probableParts.length(); ++i)
+    {
+      currentGeom = probableParts.at(i);
+      GEOSLength( currentGeom, &currentLength);
+      if(currentLength > maxLength)
+      {
+        maxLength = currentLength;
+        GEOSGeom_destroy(maxGeom);
+        maxGeom = currentGeom;
+      }
+      else
+      {
+        GEOSGeom_destroy(currentGeom);
+      }
+    }
+    resultLineParts.push_back(maxGeom);
+  }
+
+  GEOSGeom_destroy( beginLineVertex );
+  GEOSGeom_destroy( endLineVertex );
+  GEOSGeom_destroy( mergedLines );
+
+  GEOSGeometry* result = 0;
+  if ( resultLineParts.size() < 1 )
+  {
+    return 0;
+  }
+  if ( resultLineParts.size() == 1 ) //the whole result was reshaped
+  {
+    result = resultLineParts[0];
+  }
+  else //>1
+  {
+    GEOSGeometry* lineArray[resultLineParts.size()];
+    for ( int i = 0; i < resultLineParts.size(); ++i )
+    {
+      lineArray[i] = resultLineParts[i];
+    }
+
+    //create multiline from resultLineParts
+    GEOSGeometry* multiLineGeom = GEOSGeom_createCollection( GEOS_MULTILINESTRING, lineArray, resultLineParts.size() );
+
+    //then do a linemerge with the newly combined partstrings
+    result = GEOSLineMerge( multiLineGeom );
+    GEOSGeom_destroy( multiLineGeom );
+  }
+
+  //now test if the result is a linestring. Otherwise something went wrong
+  if( GEOSGeomTypeId(result) != GEOS_LINESTRING)
+  {
+    GEOSGeom_destroy(result);
+    return 0;
+  }
+  return result;
+}
+
 int QgsGeometry::topologicalTestPointsSplit( const GEOSGeometry* splitLine, QList<QgsPoint>& testPoints ) const
 {
   //Find out the intersection points between splitLineGeos and this geometry.
@@ -5020,7 +5438,7 @@
   return 0;
 }
 
-GEOSGeometry *QgsGeometry::nodeGeometries( const GEOSGeometry *splitLine, GEOSGeometry *geom ) const
+GEOSGeometry *QgsGeometry::nodeGeometries( const GEOSGeometry *splitLine, const GEOSGeometry *geom )
 {
   if ( !splitLine || !geom )
   {
@@ -5028,27 +5446,81 @@
   }
 
   GEOSGeometry *geometryBoundary = 0;
-  bool deleteBoundary = false;
   if ( GEOSGeomTypeId( geom ) == GEOS_POLYGON || GEOSGeomTypeId( geom ) == GEOS_MULTIPOLYGON )
   {
     geometryBoundary = GEOSBoundary( geom );
-    deleteBoundary = true;
   }
   else
   {
-    geometryBoundary = geom;
+    geometryBoundary = GEOSGeom_clone( geom );
   }
 
   GEOSGeometry *splitLineClone = GEOSGeom_clone( splitLine );
   GEOSGeometry *unionGeometry = GEOSUnion( splitLineClone, geometryBoundary );
   GEOSGeom_destroy( splitLineClone );
 
-  if ( deleteBoundary )
-    GEOSGeom_destroy( geometryBoundary );
-
+  GEOSGeom_destroy( geometryBoundary );
   return unionGeometry;
 }
 
+int QgsGeometry::lineContainedInLine( const GEOSGeometry* line1, const GEOSGeometry* line2 )
+{
+  if ( !line1 || !line2 )
+  {
+    return -1;
+  }
+
+  double bufferDistance = 0.00001;
+  GEOSGeometry* bufferGeom = GEOSBuffer( line2, bufferDistance, DEFAULT_QUADRANT_SEGMENTS );
+  if ( !bufferGeom )
+  {
+    return -2;
+  }
+
+  GEOSGeometry* intersectionGeom = GEOSIntersection( bufferGeom, line1 );
+
+  //compare ratio between line1Length and intersectGeomLength (usually close to 1 if line1 is contained in line2)
+  double intersectGeomLength;
+  double line1Length;
+
+  GEOSLength( intersectionGeom, &intersectGeomLength );
+  GEOSLength( line1, &line1Length );
+
+  GEOSGeom_destroy( bufferGeom );
+  GEOSGeom_destroy( intersectionGeom );
+
+  double intersectRatio = line1Length / intersectGeomLength;
+  if ( intersectRatio > 0.9 && intersectRatio < 1.1 )
+  {
+    return 1;
+  }
+  return 0;
+}
+
+int QgsGeometry::pointContainedInLine( const GEOSGeometry* point, const GEOSGeometry* line )
+{
+  if ( !point || !line )
+  {
+    return -1;
+  }
+
+  double bufferDistance = 0.0000001;
+  GEOSGeometry* lineBuffer = GEOSBuffer( line, bufferDistance, 8 );
+  if ( !lineBuffer )
+  {
+    return -2;
+  }
+
+  bool contained = false;
+  if ( GEOSContains( lineBuffer, point ) == 1 )
+  {
+    contained = true;
+  }
+
+  GEOSGeom_destroy( lineBuffer );
+  return contained;
+}
+
 int QgsGeometry::numberOfGeometries( GEOSGeometry* g ) const
 {
   if ( !g )

Modified: trunk/qgis/src/core/qgsgeometry.h
===================================================================
--- trunk/qgis/src/core/qgsgeometry.h	2009-08-25 15:09:16 UTC (rev 11499)
+++ trunk/qgis/src/core/qgsgeometry.h	2009-08-25 19:54:28 UTC (rev 11500)
@@ -252,6 +252,11 @@
                        bool topological,
                        QList<QgsPoint>& topologyTestPoints );
 
+    /**Replaces a part of this geometry with another line
+      @return 0 in case of success
+      @note: this function was added in version 1.3*/
+    int reshapeGeometry( const QList<QgsPoint>& reshapeWithLine );
+
     /**Changes this geometry such that it does not intersect the other geometry
        @param other geometry that should not be intersect
        @return 0 in case of success*/
@@ -439,10 +444,34 @@
      @return 0 in case of success*/
     int topologicalTestPointsSplit( const GEOSGeometry* splitLine, QList<QgsPoint>& testPoints ) const;
 
+    /**Creates a new line from an original line and a reshape line. The part of the input line from the first to the last intersection with the \
+        reshape line will be replaced. The calling function takes ownership of the result.
+    @param origLine the original line
+    @param reshapeLineGeos the reshape line
+    @return the reshaped line or 0 in case of error*/
+    static GEOSGeometry* reshapeLine( const GEOSGeometry* origLine, const GEOSGeometry* reshapeLineGeos );
+
+    /**Creates a new polygon replacing the part from the first to the second intersection with the reshape line. As a polygon ring is always closed,
+        the method keeps the longer part of the existing boundary
+    @param polygon geometry to reshape
+    @param reshapeLineGeos the reshape line
+    @return the reshaped polygon or 0 in case of error*/
+    static GEOSGeometry* reshapePolygon( const GEOSGeometry* polygon, const GEOSGeometry* reshapeLineGeos );
+
     /**Nodes together a split line and a (multi-) polygon geometry in a multilinestring
      @return the noded multiline geometry or 0 in case of error. The calling function takes ownership of the node geometry*/
-    GEOSGeometry* nodeGeometries( const GEOSGeometry *splitLine, GEOSGeometry *poly ) const;
+    static GEOSGeometry* nodeGeometries( const GEOSGeometry *splitLine, const GEOSGeometry *poly );
 
+    /**Tests if line1 is completely contained in line2. This method works with a very small buffer around line2 and therefore is less prone
+        to numerical errors as the corresponding geos method*/
+    static int lineContainedInLine( const GEOSGeometry* line1, const GEOSGeometry* line2 );
+
+    /**Tests if a point is on a line. This method works with a very small buffer and is thus less prone to numerical problems as the direct geos functions.
+      @param point the point to test
+      @param line line to test
+      @return 0 not contained, 1 if contained, <0 in case of error*/
+    static int pointContainedInLine( const GEOSGeometry* point, const GEOSGeometry* line );
+
     /**Returns number of single geometry in a geos geometry. Is save for geos 2 and 3*/
     int numberOfGeometries( GEOSGeometry* g ) const;
 



More information about the QGIS-commit mailing list