QGIS API Documentation  2.8.2-Wien
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
qgsgeometryanalyzer.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsgeometryanalyzer.cpp - QGIS Tools for vector geometry analysis
3  -------------------
4  begin : 19 March 2009
5  copyright : (C) Carson Farmer
6  email : [email protected]
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #include "qgsgeometryanalyzer.h"
19 
20 #include "qgsapplication.h"
21 #include "qgsfield.h"
22 #include "qgsfeature.h"
23 #include "qgslogger.h"
25 #include "qgsvectorfilewriter.h"
26 #include "qgsvectordataprovider.h"
27 #include "qgsdistancearea.h"
28 #include <QProgressDialog>
29 
31  const QString& shapefileName,
32  double tolerance,
33  bool onlySelectedFeatures,
34  QProgressDialog *p )
35 {
36  if ( !layer )
37  {
38  return false;
39  }
40 
41  QgsVectorDataProvider* dp = layer->dataProvider();
42  if ( !dp )
43  {
44  return false;
45  }
46 
47  QGis::WkbType outputType = dp->geometryType();
48  const QgsCoordinateReferenceSystem crs = layer->crs();
49 
50  QgsVectorFileWriter vWriter( shapefileName, dp->encoding(), layer->pendingFields(), outputType, &crs );
51  QgsFeature currentFeature;
52 
53  //take only selection
54  if ( onlySelectedFeatures )
55  {
56  //use QgsVectorLayer::featureAtId
57  const QgsFeatureIds selection = layer->selectedFeaturesIds();
58  if ( p )
59  {
60  p->setMaximum( selection.size() );
61  }
62 
63  int processedFeatures = 0;
64  QgsFeatureIds::const_iterator it = selection.constBegin();
65  for ( ; it != selection.constEnd(); ++it )
66  {
67  if ( p )
68  {
69  p->setValue( processedFeatures );
70  }
71 
72  if ( p && p->wasCanceled() )
73  {
74  break;
75  }
76  if ( !layer->getFeatures( QgsFeatureRequest().setFilterFid( *it ) ).nextFeature( currentFeature ) )
77  {
78  continue;
79  }
80  simplifyFeature( currentFeature, &vWriter, tolerance );
81  ++processedFeatures;
82  }
83 
84  if ( p )
85  {
86  p->setValue( selection.size() );
87  }
88  }
89  //take all features
90  else
91  {
92  QgsFeatureIterator fit = layer->getFeatures();
93 
94  int featureCount = layer->featureCount();
95  if ( p )
96  {
97  p->setMaximum( featureCount );
98  }
99  int processedFeatures = 0;
100 
101  while ( fit.nextFeature( currentFeature ) )
102  {
103  if ( p )
104  {
105  p->setValue( processedFeatures );
106  }
107  if ( p && p->wasCanceled() )
108  {
109  break;
110  }
111  simplifyFeature( currentFeature, &vWriter, tolerance );
112  ++processedFeatures;
113  }
114  if ( p )
115  {
116  p->setValue( featureCount );
117  }
118  }
119 
120  return true;
121 }
122 
123 void QgsGeometryAnalyzer::simplifyFeature( QgsFeature& f, QgsVectorFileWriter* vfw, double tolerance )
124 {
125  QgsGeometry* featureGeometry = f.geometry();
126  QgsGeometry* tmpGeometry = 0;
127 
128  if ( !featureGeometry )
129  {
130  return;
131  }
132  // simplify feature
133  tmpGeometry = featureGeometry->simplify( tolerance );
134 
135  QgsFeature newFeature;
136  newFeature.setGeometry( tmpGeometry );
137  newFeature.setAttributes( f.attributes() );
138 
139  //add it to vector file writer
140  if ( vfw )
141  {
142  vfw->addFeature( newFeature );
143  }
144 }
145 
146 bool QgsGeometryAnalyzer::centroids( QgsVectorLayer* layer, const QString& shapefileName,
147  bool onlySelectedFeatures, QProgressDialog* p )
148 {
149  if ( !layer )
150  {
151  QgsDebugMsg( "No layer passed to centroids" );
152  return false;
153  }
154 
155  QgsVectorDataProvider* dp = layer->dataProvider();
156  if ( !dp )
157  {
158  QgsDebugMsg( "No data provider for layer passed to centroids" );
159  return false;
160  }
161 
162  QGis::WkbType outputType = QGis::WKBPoint;
163  const QgsCoordinateReferenceSystem crs = layer->crs();
164 
165  QgsVectorFileWriter vWriter( shapefileName, dp->encoding(), layer->pendingFields(), outputType, &crs );
166  QgsFeature currentFeature;
167 
168  //take only selection
169  if ( onlySelectedFeatures )
170  {
171  //use QgsVectorLayer::featureAtId
172  const QgsFeatureIds selection = layer->selectedFeaturesIds();
173  if ( p )
174  {
175  p->setMaximum( selection.size() );
176  }
177 
178  int processedFeatures = 0;
179  QgsFeatureIds::const_iterator it = selection.constBegin();
180  for ( ; it != selection.constEnd(); ++it )
181  {
182  if ( p )
183  {
184  p->setValue( processedFeatures );
185  }
186 
187  if ( p && p->wasCanceled() )
188  {
189  break;
190  }
191  if ( !layer->getFeatures( QgsFeatureRequest().setFilterFid( *it ) ).nextFeature( currentFeature ) )
192  {
193  continue;
194  }
195  centroidFeature( currentFeature, &vWriter );
196  ++processedFeatures;
197  }
198 
199  if ( p )
200  {
201  p->setValue( selection.size() );
202  }
203  }
204  //take all features
205  else
206  {
207  QgsFeatureIterator fit = layer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QgsAttributeList() ) );
208 
209  int featureCount = layer->featureCount();
210  if ( p )
211  {
212  p->setMaximum( featureCount );
213  }
214  int processedFeatures = 0;
215 
216  while ( fit.nextFeature( currentFeature ) )
217  {
218  if ( p )
219  {
220  p->setValue( processedFeatures );
221  }
222  if ( p && p->wasCanceled() )
223  {
224  break;
225  }
226  centroidFeature( currentFeature, &vWriter );
227  ++processedFeatures;
228  }
229  if ( p )
230  {
231  p->setValue( featureCount );
232  }
233  }
234 
235  return true;
236 }
237 
238 
239 void QgsGeometryAnalyzer::centroidFeature( QgsFeature& f, QgsVectorFileWriter* vfw )
240 {
241  QgsGeometry* featureGeometry = f.geometry();
242  QgsGeometry* tmpGeometry = 0;
243 
244  if ( !featureGeometry )
245  {
246  return;
247  }
248 
249  tmpGeometry = featureGeometry->centroid();
250 
251  QgsFeature newFeature;
252  newFeature.setGeometry( tmpGeometry );
253  newFeature.setAttributes( f.attributes() );
254 
255  //add it to vector file writer
256  if ( vfw )
257  {
258  vfw->addFeature( newFeature );
259  }
260 }
261 
263  const QString& shapefileName,
264  bool onlySelectedFeatures,
265  QProgressDialog * )
266 {
267  if ( !layer )
268  {
269  return false;
270  }
271 
272  QgsVectorDataProvider* dp = layer->dataProvider();
273  if ( !dp )
274  {
275  return false;
276  }
277 
278  QGis::WkbType outputType = QGis::WKBPolygon;
279  const QgsCoordinateReferenceSystem crs = layer->crs();
280 
281  QgsFields fields;
282  fields.append( QgsField( QString( "MINX" ), QVariant::Double ) );
283  fields.append( QgsField( QString( "MINY" ), QVariant::Double ) );
284  fields.append( QgsField( QString( "MAXX" ), QVariant::Double ) );
285  fields.append( QgsField( QString( "MAXY" ), QVariant::Double ) );
286  fields.append( QgsField( QString( "CNTX" ), QVariant::Double ) );
287  fields.append( QgsField( QString( "CNTY" ), QVariant::Double ) );
288  fields.append( QgsField( QString( "AREA" ), QVariant::Double ) );
289  fields.append( QgsField( QString( "PERIM" ), QVariant::Double ) );
290  fields.append( QgsField( QString( "HEIGHT" ), QVariant::Double ) );
291  fields.append( QgsField( QString( "WIDTH" ), QVariant::Double ) );
292 
293  QgsVectorFileWriter vWriter( shapefileName, dp->encoding(), fields, outputType, &crs );
294 
295  QgsRectangle rect;
296  if ( onlySelectedFeatures ) // take only selection
297  {
298  rect = layer->boundingBoxOfSelected();
299  }
300  else
301  {
302  rect = layer->extent();
303  }
304 
305  double minx = rect.xMinimum();
306  double miny = rect.yMinimum();
307  double maxx = rect.xMaximum();
308  double maxy = rect.yMaximum();
309  double height = rect.height();
310  double width = rect.width();
311  double cntx = minx + ( width / 2.0 );
312  double cnty = miny + ( height / 2.0 );
313  double area = width * height;
314  double perim = ( 2 * width ) + ( 2 * height );
315 
316  QgsFeature feat;
317  QgsAttributes attrs( 10 );
318  attrs[0] = QVariant( minx );
319  attrs[1] = QVariant( miny );
320  attrs[2] = QVariant( maxx );
321  attrs[3] = QVariant( maxy );
322  attrs[4] = QVariant( cntx );
323  attrs[5] = QVariant( cnty );
324  attrs[6] = QVariant( area );
325  attrs[7] = QVariant( perim );
326  attrs[8] = QVariant( height );
327  attrs[9] = QVariant( width );
328  feat.setAttributes( attrs );
329  feat.setGeometry( QgsGeometry::fromRect( rect ) );
330  vWriter.addFeature( feat );
331  return true;
332 }
333 
334 QList<double> QgsGeometryAnalyzer::simpleMeasure( QgsGeometry* mpGeometry )
335 {
336  QList<double> list;
337  double perim;
338  if ( mpGeometry->wkbType() == QGis::WKBPoint )
339  {
340  QgsPoint pt = mpGeometry->asPoint();
341  list.append( pt.x() );
342  list.append( pt.y() );
343  }
344  else
345  {
346  QgsDistanceArea measure;
347  list.append( measure.measure( mpGeometry ) );
348  if ( mpGeometry->type() == QGis::Polygon )
349  {
350  perim = perimeterMeasure( mpGeometry, measure );
351  list.append( perim );
352  }
353  }
354  return list;
355 }
356 
357 double QgsGeometryAnalyzer::perimeterMeasure( QgsGeometry* geometry, QgsDistanceArea& measure )
358 {
359  double value = 0.00;
360  if ( geometry->isMultipart() )
361  {
362  QgsMultiPolygon poly = geometry->asMultiPolygon();
363  QgsMultiPolygon::iterator it;
364  QgsPolygon::iterator jt;
365  for ( it = poly.begin(); it != poly.end(); ++it )
366  {
367  for ( jt = it->begin(); jt != it->end(); ++jt )
368  {
369  QgsGeometry* geom = QgsGeometry::fromPolyline( *jt );
370  value = value + measure.measure( geom );
371  delete geom;
372  }
373  }
374  }
375  else
376  {
377  QgsPolygon::iterator jt;
378  QgsPolygon poly = geometry->asPolygon();
379  for ( jt = poly.begin(); jt != poly.end(); ++jt )
380  {
381  QgsGeometry* geom = QgsGeometry::fromPolyline( *jt );
382  value = value + measure.measure( geom );
383  delete geom;
384  }
385  }
386  return value;
387 }
388 
389 bool QgsGeometryAnalyzer::convexHull( QgsVectorLayer* layer, const QString& shapefileName,
390  bool onlySelectedFeatures, int uniqueIdField, QProgressDialog* p )
391 {
392  if ( !layer )
393  {
394  return false;
395  }
396  QgsVectorDataProvider* dp = layer->dataProvider();
397  if ( !dp )
398  {
399  return false;
400  }
401  bool useField = false;
402  if ( uniqueIdField == -1 )
403  {
404  uniqueIdField = 0;
405  }
406  else
407  {
408  useField = true;
409  }
410  QgsFields fields;
411  fields.append( QgsField( QString( "UID" ), QVariant::String ) );
412  fields.append( QgsField( QString( "AREA" ), QVariant::Double ) );
413  fields.append( QgsField( QString( "PERIM" ), QVariant::Double ) );
414 
415  QGis::WkbType outputType = QGis::WKBPolygon;
416  const QgsCoordinateReferenceSystem crs = layer->crs();
417 
418  QgsVectorFileWriter vWriter( shapefileName, dp->encoding(), fields, outputType, &crs );
419  QgsFeature currentFeature;
420  QgsGeometry* dissolveGeometry = 0; //dissolve geometry
421  QMultiMap<QString, QgsFeatureId> map;
422 
423  if ( onlySelectedFeatures )
424  {
425  //use QgsVectorLayer::featureAtId
426  const QgsFeatureIds selection = layer->selectedFeaturesIds();
427  QgsFeatureIds::const_iterator it = selection.constBegin();
428  for ( ; it != selection.constEnd(); ++it )
429  {
430 #if 0
431  if ( p )
432  {
433  p->setValue( processedFeatures );
434  }
435  if ( p && p->wasCanceled() )
436  {
437  // break; // it may be better to do something else here?
438  return false;
439  }
440 #endif
441  if ( !layer->getFeatures( QgsFeatureRequest().setFilterFid( *it ) ).nextFeature( currentFeature ) )
442  {
443  continue;
444  }
445  map.insert( currentFeature.attribute( uniqueIdField ).toString(), currentFeature.id() );
446  }
447  }
448  else
449  {
450  QgsFeatureIterator fit = layer->getFeatures();
451  while ( fit.nextFeature( currentFeature ) )
452  {
453 #if 0
454  if ( p )
455  {
456  p->setValue( processedFeatures );
457  }
458  if ( p && p->wasCanceled() )
459  {
460  // break; // it may be better to do something else here?
461  return false;
462  }
463 #endif
464  map.insert( currentFeature.attribute( uniqueIdField ).toString(), currentFeature.id() );
465  }
466  }
467 
468  QMultiMap<QString, QgsFeatureId>::const_iterator jt = map.constBegin();
469  while ( jt != map.constEnd() )
470  {
471  QString currentKey = jt.key();
472  int processedFeatures = 0;
473  //take only selection
474  if ( onlySelectedFeatures )
475  {
476  //use QgsVectorLayer::featureAtId
477  const QgsFeatureIds selection = layer->selectedFeaturesIds();
478  if ( p )
479  {
480  p->setMaximum( selection.size() );
481  }
482  processedFeatures = 0;
483  while ( jt != map.constEnd() && ( jt.key() == currentKey || !useField ) )
484  {
485  if ( p && p->wasCanceled() )
486  {
487  break;
488  }
489  if ( selection.contains( jt.value() ) )
490  {
491  if ( p )
492  {
493  p->setValue( processedFeatures );
494  }
495  if ( !layer->getFeatures( QgsFeatureRequest().setFilterFid( jt.value() ) ).nextFeature( currentFeature ) )
496  {
497  continue;
498  }
499  convexFeature( currentFeature, processedFeatures, &dissolveGeometry );
500  ++processedFeatures;
501  }
502  ++jt;
503  }
504  QList<double> values;
505  if ( !dissolveGeometry )
506  {
507  QgsDebugMsg( "no dissolved geometry - should not happen" );
508  return false;
509  }
510  dissolveGeometry = dissolveGeometry->convexHull();
511  values = simpleMeasure( dissolveGeometry );
512  QgsAttributes attributes( 3 );
513  attributes[0] = QVariant( currentKey );
514  attributes[1] = values[ 0 ];
515  attributes[2] = values[ 1 ];
516  QgsFeature dissolveFeature;
517  dissolveFeature.setAttributes( attributes );
518  dissolveFeature.setGeometry( dissolveGeometry );
519  vWriter.addFeature( dissolveFeature );
520  }
521  //take all features
522  else
523  {
524  int featureCount = layer->featureCount();
525  if ( p )
526  {
527  p->setMaximum( featureCount );
528  }
529  processedFeatures = 0;
530  while ( jt != map.constEnd() && ( jt.key() == currentKey || !useField ) )
531  {
532  if ( p )
533  {
534  p->setValue( processedFeatures );
535  }
536 
537  if ( p && p->wasCanceled() )
538  {
539  break;
540  }
541  if ( !layer->getFeatures( QgsFeatureRequest().setFilterFid( jt.value() ) ).nextFeature( currentFeature ) )
542  {
543  continue;
544  }
545  convexFeature( currentFeature, processedFeatures, &dissolveGeometry );
546  ++processedFeatures;
547  ++jt;
548  }
549  QList<double> values;
550  // QgsGeometry* tmpGeometry = 0;
551  if ( !dissolveGeometry )
552  {
553  QgsDebugMsg( "no dissolved geometry - should not happen" );
554  return false;
555  }
556  dissolveGeometry = dissolveGeometry->convexHull();
557  // values = simpleMeasure( tmpGeometry );
558  values = simpleMeasure( dissolveGeometry );
559  QgsAttributes attributes;
560  attributes[0] = QVariant( currentKey );
561  attributes[1] = QVariant( values[ 0 ] );
562  attributes[2] = QVariant( values[ 1 ] );
563  QgsFeature dissolveFeature;
564  dissolveFeature.setAttributes( attributes );
565  dissolveFeature.setGeometry( dissolveGeometry );
566  vWriter.addFeature( dissolveFeature );
567  }
568  }
569  return true;
570 }
571 
572 
573 void QgsGeometryAnalyzer::convexFeature( QgsFeature& f, int nProcessedFeatures, QgsGeometry** dissolveGeometry )
574 {
575  QgsGeometry* featureGeometry = f.geometry();
576  QgsGeometry* tmpGeometry = 0;
577  QgsGeometry* convexGeometry = 0;
578 
579  if ( !featureGeometry )
580  {
581  return;
582  }
583 
584  convexGeometry = featureGeometry->convexHull();
585 
586  if ( nProcessedFeatures == 0 )
587  {
588  *dissolveGeometry = convexGeometry;
589  }
590  else
591  {
592  tmpGeometry = *dissolveGeometry;
593  *dissolveGeometry = ( *dissolveGeometry )->combine( convexGeometry );
594  delete tmpGeometry;
595  delete convexGeometry;
596  }
597 }
598 
599 bool QgsGeometryAnalyzer::dissolve( QgsVectorLayer* layer, const QString& shapefileName,
600  bool onlySelectedFeatures, int uniqueIdField, QProgressDialog* p )
601 {
602  if ( !layer )
603  {
604  return false;
605  }
606  QgsVectorDataProvider* dp = layer->dataProvider();
607  if ( !dp )
608  {
609  return false;
610  }
611  bool useField = false;
612  if ( uniqueIdField == -1 )
613  {
614  uniqueIdField = 0;
615  }
616  else
617  {
618  useField = true;
619  }
620 
621  QGis::WkbType outputType = dp->geometryType();
622  const QgsCoordinateReferenceSystem crs = layer->crs();
623 
624  QgsVectorFileWriter vWriter( shapefileName, dp->encoding(), layer->pendingFields(), outputType, &crs );
625  QgsFeature currentFeature;
626  QMultiMap<QString, QgsFeatureId> map;
627 
628  if ( onlySelectedFeatures )
629  {
630  //use QgsVectorLayer::featureAtId
631  const QgsFeatureIds selection = layer->selectedFeaturesIds();
632  QgsFeatureIds::const_iterator it = selection.constBegin();
633  for ( ; it != selection.constEnd(); ++it )
634  {
635  if ( !layer->getFeatures( QgsFeatureRequest().setFilterFid( *it ) ).nextFeature( currentFeature ) )
636  {
637  continue;
638  }
639  map.insert( currentFeature.attribute( uniqueIdField ).toString(), currentFeature.id() );
640  }
641  }
642  else
643  {
644  QgsFeatureIterator fit = layer->getFeatures();
645  while ( fit.nextFeature( currentFeature ) )
646  {
647  map.insert( currentFeature.attribute( uniqueIdField ).toString(), currentFeature.id() );
648  }
649  }
650 
651  QgsGeometry *dissolveGeometry = 0; //dissolve geometry
652  QMultiMap<QString, QgsFeatureId>::const_iterator jt = map.constBegin();
653  QgsFeature outputFeature;
654  while ( jt != map.constEnd() )
655  {
656  QString currentKey = jt.key();
657  int processedFeatures = 0;
658  bool first = true;
659  //take only selection
660  if ( onlySelectedFeatures )
661  {
662  //use QgsVectorLayer::featureAtId
663  const QgsFeatureIds selection = layer->selectedFeaturesIds();
664  if ( p )
665  {
666  p->setMaximum( selection.size() );
667  }
668  while ( jt != map.constEnd() && ( jt.key() == currentKey || !useField ) )
669  {
670  if ( p && p->wasCanceled() )
671  {
672  break;
673  }
674  if ( selection.contains( jt.value() ) )
675  {
676  if ( p )
677  {
678  p->setValue( processedFeatures );
679  }
680  if ( !layer->getFeatures( QgsFeatureRequest().setFilterFid( jt.value() ) ).nextFeature( currentFeature ) )
681  {
682  continue;
683  }
684  if ( first )
685  {
686  outputFeature.setAttributes( currentFeature.attributes() );
687  first = false;
688  }
689  dissolveFeature( currentFeature, processedFeatures, &dissolveGeometry );
690  ++processedFeatures;
691  }
692  ++jt;
693  }
694  }
695  //take all features
696  else
697  {
698  int featureCount = layer->featureCount();
699  if ( p )
700  {
701  p->setMaximum( featureCount );
702  }
703  while ( jt != map.constEnd() && ( jt.key() == currentKey || !useField ) )
704  {
705  if ( p )
706  {
707  p->setValue( processedFeatures );
708  }
709 
710  if ( p && p->wasCanceled() )
711  {
712  break;
713  }
714  if ( !layer->getFeatures( QgsFeatureRequest().setFilterFid( jt.value() ) ).nextFeature( currentFeature ) )
715  {
716  continue;
717  }
718  {
719  outputFeature.setAttributes( currentFeature.attributes() );
720  first = false;
721  }
722  dissolveFeature( currentFeature, processedFeatures, &dissolveGeometry );
723  ++processedFeatures;
724  ++jt;
725  }
726  }
727  outputFeature.setGeometry( dissolveGeometry );
728  vWriter.addFeature( outputFeature );
729  }
730  return true;
731 }
732 
733 void QgsGeometryAnalyzer::dissolveFeature( QgsFeature& f, int nProcessedFeatures, QgsGeometry** dissolveGeometry )
734 {
735  QgsGeometry* featureGeometry = f.geometry();
736 
737  if ( !featureGeometry )
738  {
739  return;
740  }
741 
742  if ( nProcessedFeatures == 0 )
743  {
744  size_t geomSize = featureGeometry->wkbSize();
745  *dissolveGeometry = new QgsGeometry();
746  unsigned char* wkb = new unsigned char[geomSize];
747  memcpy( wkb, featureGeometry->asWkb(), geomSize );
748  ( *dissolveGeometry )->fromWkb( wkb, geomSize );
749  }
750  else
751  {
752  *dissolveGeometry = ( *dissolveGeometry )->combine( featureGeometry );
753  }
754 }
755 
756 bool QgsGeometryAnalyzer::buffer( QgsVectorLayer* layer, const QString& shapefileName, double bufferDistance,
757  bool onlySelectedFeatures, bool dissolve, int bufferDistanceField, QProgressDialog* p )
758 {
759  if ( !layer )
760  {
761  return false;
762  }
763 
764  QgsVectorDataProvider* dp = layer->dataProvider();
765  if ( !dp )
766  {
767  return false;
768  }
769 
770  QGis::WkbType outputType = QGis::WKBPolygon;
771  if ( dissolve )
772  {
773  outputType = QGis::WKBMultiPolygon;
774  }
775  const QgsCoordinateReferenceSystem crs = layer->crs();
776 
777  QgsVectorFileWriter vWriter( shapefileName, dp->encoding(), layer->pendingFields(), outputType, &crs );
778  QgsFeature currentFeature;
779  QgsGeometry *dissolveGeometry = 0; //dissolve geometry (if dissolve enabled)
780 
781  //take only selection
782  if ( onlySelectedFeatures )
783  {
784  //use QgsVectorLayer::featureAtId
785  const QgsFeatureIds selection = layer->selectedFeaturesIds();
786  if ( p )
787  {
788  p->setMaximum( selection.size() );
789  }
790 
791  int processedFeatures = 0;
792  QgsFeatureIds::const_iterator it = selection.constBegin();
793  for ( ; it != selection.constEnd(); ++it )
794  {
795  if ( p )
796  {
797  p->setValue( processedFeatures );
798  }
799 
800  if ( p && p->wasCanceled() )
801  {
802  break;
803  }
804  if ( !layer->getFeatures( QgsFeatureRequest().setFilterFid( *it ) ).nextFeature( currentFeature ) )
805  {
806  continue;
807  }
808  bufferFeature( currentFeature, processedFeatures, &vWriter, dissolve, &dissolveGeometry, bufferDistance, bufferDistanceField );
809  ++processedFeatures;
810  }
811 
812  if ( p )
813  {
814  p->setValue( selection.size() );
815  }
816  }
817  //take all features
818  else
819  {
820  QgsFeatureIterator fit = layer->getFeatures();
821 
822  int featureCount = layer->featureCount();
823  if ( p )
824  {
825  p->setMaximum( featureCount );
826  }
827  int processedFeatures = 0;
828 
829  while ( fit.nextFeature( currentFeature ) )
830  {
831  if ( p )
832  {
833  p->setValue( processedFeatures );
834  }
835  if ( p && p->wasCanceled() )
836  {
837  break;
838  }
839  bufferFeature( currentFeature, processedFeatures, &vWriter, dissolve, &dissolveGeometry, bufferDistance, bufferDistanceField );
840  ++processedFeatures;
841  }
842  if ( p )
843  {
844  p->setValue( featureCount );
845  }
846  }
847 
848  if ( dissolve )
849  {
850  QgsFeature dissolveFeature;
851  if ( !dissolveGeometry )
852  {
853  QgsDebugMsg( "no dissolved geometry - should not happen" );
854  return false;
855  }
856  dissolveFeature.setGeometry( dissolveGeometry );
857  vWriter.addFeature( dissolveFeature );
858  }
859  return true;
860 }
861 
862 void QgsGeometryAnalyzer::bufferFeature( QgsFeature& f, int nProcessedFeatures, QgsVectorFileWriter* vfw, bool dissolve,
863  QgsGeometry** dissolveGeometry, double bufferDistance, int bufferDistanceField )
864 {
865  double currentBufferDistance;
866  QgsGeometry* featureGeometry = f.geometry();
867  QgsGeometry* tmpGeometry = 0;
868  QgsGeometry* bufferGeometry = 0;
869 
870  if ( !featureGeometry )
871  {
872  return;
873  }
874 
875  //create buffer
876  if ( bufferDistanceField == -1 )
877  {
878  currentBufferDistance = bufferDistance;
879  }
880  else
881  {
882  currentBufferDistance = f.attribute( bufferDistanceField ).toDouble();
883  }
884  bufferGeometry = featureGeometry->buffer( currentBufferDistance, 5 );
885 
886  if ( dissolve )
887  {
888  if ( nProcessedFeatures == 0 )
889  {
890  *dissolveGeometry = bufferGeometry;
891  }
892  else
893  {
894  tmpGeometry = *dissolveGeometry;
895  *dissolveGeometry = ( *dissolveGeometry )->combine( bufferGeometry );
896  delete tmpGeometry;
897  delete bufferGeometry;
898  }
899  }
900  else //dissolve
901  {
902  QgsFeature newFeature;
903  newFeature.setGeometry( bufferGeometry );
904  newFeature.setAttributes( f.attributes() );
905 
906  //add it to vector file writer
907  if ( vfw )
908  {
909  vfw->addFeature( newFeature );
910  }
911  }
912 }
913 
914 bool QgsGeometryAnalyzer::eventLayer( QgsVectorLayer* lineLayer, QgsVectorLayer* eventLayer, int lineField, int eventField, QgsFeatureIds &unlocatedFeatureIds, const QString& outputLayer,
915  const QString& outputFormat, int locationField1, int locationField2, int offsetField, double offsetScale,
916  bool forceSingleGeometry, QgsVectorDataProvider* memoryProvider, QProgressDialog* p )
917 {
918  if ( !lineLayer || !eventLayer || !lineLayer->isValid() || !eventLayer->isValid() )
919  {
920  return false;
921  }
922 
923  //create line field / id map for line layer
924  QMultiHash< QString, QgsFeature > lineLayerIdMap; //1:n possible (e.g. several linear reference geometries for one feature in the event layer)
925  QgsFeatureIterator fit = lineLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QgsAttributeList() << lineField ) );
926  QgsFeature fet;
927  while ( fit.nextFeature( fet ) )
928  {
929  lineLayerIdMap.insert( fet.attribute( lineField ).toString(), fet );
930  }
931 
932  //create output datasource or attributes in memory provider
933  QgsVectorFileWriter* fileWriter = 0;
934  QgsFeatureList memoryProviderFeatures;
935  if ( !memoryProvider )
936  {
937  QGis::WkbType memoryProviderType = QGis::WKBMultiLineString;
938  if ( locationField2 == -1 )
939  {
940  memoryProviderType = forceSingleGeometry ? QGis::WKBPoint : QGis::WKBMultiPoint;
941  }
942  else
943  {
944  memoryProviderType = forceSingleGeometry ? QGis::WKBLineString : QGis::WKBMultiLineString;
945  }
946  fileWriter = new QgsVectorFileWriter( outputLayer,
947  eventLayer->dataProvider()->encoding(),
948  eventLayer->pendingFields(),
949  memoryProviderType,
950  &( lineLayer->crs() ),
951  outputFormat );
952  }
953  else
954  {
955  memoryProvider->addAttributes( eventLayer->pendingFields().toList() );
956  }
957 
958  //iterate over eventLayer and write new features to output file or layer
959  fit = eventLayer->getFeatures( QgsFeatureRequest().setFlags( QgsFeatureRequest::NoGeometry ) );
960  QgsGeometry* lrsGeom = 0;
961  double measure1, measure2 = 0.0;
962 
963  int nEventFeatures = eventLayer->pendingFeatureCount();
964  int featureCounter = 0;
965  int nOutputFeatures = 0; //number of output features for the current event feature
966  if ( p )
967  {
968  p->setWindowModality( Qt::WindowModal );
969  p->setMinimum( 0 );
970  p->setMaximum( nEventFeatures );
971  p->show();
972  }
973 
974  while ( fit.nextFeature( fet ) )
975  {
976  nOutputFeatures = 0;
977 
978  //update progress dialog
979  if ( p )
980  {
981  if ( p->wasCanceled() )
982  {
983  break;
984  }
985  p->setValue( featureCounter );
986  ++featureCounter;
987  }
988 
989  measure1 = fet.attribute( locationField1 ).toDouble();
990  if ( locationField2 != -1 )
991  {
992  measure2 = fet.attribute( locationField2 ).toDouble();
993  if ( qgsDoubleNear(( measure2 - measure1 ), 0.0 ) )
994  {
995  continue;
996  }
997  }
998 
999  QList<QgsFeature> featureIdList = lineLayerIdMap.values( fet.attribute( eventField ).toString() );
1000  QList<QgsFeature>::const_iterator featureIdIt = featureIdList.constBegin();
1001  for ( ; featureIdIt != featureIdList.constEnd(); ++featureIdIt )
1002  {
1003  if ( locationField2 == -1 )
1004  {
1005  lrsGeom = locateAlongMeasure( measure1, featureIdIt->geometry() );
1006  }
1007  else
1008  {
1009  lrsGeom = locateBetweenMeasures( measure1, measure2, featureIdIt->geometry() );
1010  }
1011 
1012  if ( lrsGeom )
1013  {
1014  ++nOutputFeatures;
1015  addEventLayerFeature( fet, lrsGeom, featureIdIt->geometry(), fileWriter, memoryProviderFeatures, offsetField, offsetScale, forceSingleGeometry );
1016  }
1017  }
1018  if ( nOutputFeatures < 1 )
1019  {
1020  unlocatedFeatureIds.insert( fet.id() );
1021  }
1022  }
1023 
1024  if ( p )
1025  {
1026  p->setValue( nEventFeatures );
1027  }
1028 
1029  if ( memoryProvider )
1030  {
1031  memoryProvider->addFeatures( memoryProviderFeatures );
1032  }
1033  delete fileWriter;
1034  return true;
1035 }
1036 
1037 void QgsGeometryAnalyzer::addEventLayerFeature( QgsFeature& feature, QgsGeometry* geom, QgsGeometry* lineGeom, QgsVectorFileWriter* fileWriter, QgsFeatureList& memoryFeatures,
1038  int offsetField, double offsetScale, bool forceSingleType )
1039 {
1040  if ( !geom )
1041  {
1042  return;
1043  }
1044 
1045  QList<QgsGeometry*> geomList;
1046  if ( forceSingleType )
1047  {
1048  geomList = geom->asGeometryCollection();
1049  }
1050  else
1051  {
1052  geomList.push_back( geom );
1053  }
1054 
1055  QList<QgsGeometry*>::iterator geomIt = geomList.begin();
1056  for ( ; geomIt != geomList.end(); ++geomIt )
1057  {
1058  //consider offset
1059  if ( offsetField >= 0 )
1060  {
1061  double offsetVal = feature.attribute( offsetField ).toDouble();
1062  offsetVal *= offsetScale;
1063  if ( !createOffsetGeometry( *geomIt, lineGeom, offsetVal ) )
1064  {
1065  delete *geomIt;
1066  continue;
1067  }
1068  }
1069 
1070  feature.setGeometry( *geomIt );
1071  if ( fileWriter )
1072  {
1073  fileWriter->addFeature( feature );
1074  }
1075  else
1076  {
1077  memoryFeatures << feature;
1078  }
1079  }
1080 
1081  if ( forceSingleType )
1082  {
1083  delete geom;
1084  }
1085 }
1086 
1087 bool QgsGeometryAnalyzer::createOffsetGeometry( QgsGeometry* geom, QgsGeometry* lineGeom, double offset )
1088 {
1089  if ( !geom || !lineGeom )
1090  {
1091  return false;
1092  }
1093 
1094  QList<QgsGeometry*> inputGeomList;
1095 
1096  if ( geom->isMultipart() )
1097  {
1098  inputGeomList = geom->asGeometryCollection();
1099  }
1100  else
1101  {
1102  inputGeomList.push_back( geom );
1103  }
1104 
1105  QList<GEOSGeometry*> outputGeomList;
1106  QList<QgsGeometry*>::const_iterator inputGeomIt = inputGeomList.constBegin();
1107  GEOSContextHandle_t geosctxt = QgsGeometry::getGEOSHandler();
1108  for ( ; inputGeomIt != inputGeomList.constEnd(); ++inputGeomIt )
1109  {
1110  if ( geom->type() == QGis::Line )
1111  {
1112  //geos 3.3 needed for line offsets
1113 #if defined(GEOS_VERSION_MAJOR) && defined(GEOS_VERSION_MINOR) && \
1114  ((GEOS_VERSION_MAJOR>3) || ((GEOS_VERSION_MAJOR==3) && (GEOS_VERSION_MINOR>=3)))
1115  GEOSGeometry* offsetGeom = GEOSOffsetCurve_r( geosctxt, ( *inputGeomIt )->asGeos(), -offset, 8 /*quadSegments*/, 0 /*joinStyle*/, 5.0 /*mitreLimit*/ );
1116  if ( !offsetGeom || !GEOSisValid_r( geosctxt, offsetGeom ) )
1117  {
1118  return false;
1119  }
1120  if ( !GEOSisValid_r( geosctxt, offsetGeom ) || GEOSGeomTypeId_r( geosctxt, offsetGeom ) != GEOS_LINESTRING || GEOSGeomGetNumPoints_r( geosctxt, offsetGeom ) < 1 )
1121  {
1122  GEOSGeom_destroy_r( geosctxt, offsetGeom );
1123  return false;
1124  }
1125  outputGeomList.push_back( offsetGeom );
1126 #else
1127  outputGeomList.push_back( GEOSGeom_clone_r( geosctxt, ( *inputGeomIt )->asGeos() ) );
1128 #endif
1129  }
1130  else if ( geom->type() == QGis::Point )
1131  {
1132  QgsPoint p = ( *inputGeomIt )->asPoint();
1133  p = createPointOffset( p.x(), p.y(), offset, lineGeom );
1134  GEOSCoordSequence* ptSeq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
1135  GEOSCoordSeq_setX_r( geosctxt, ptSeq, 0, p.x() );
1136  GEOSCoordSeq_setY_r( geosctxt, ptSeq, 0, p.y() );
1137  GEOSGeometry* geosPt = GEOSGeom_createPoint_r( geosctxt, ptSeq );
1138  outputGeomList.push_back( geosPt );
1139  }
1140  }
1141 
1142  if ( !geom->isMultipart() )
1143  {
1144  GEOSGeometry* outputGeom = outputGeomList.at( 0 );
1145  if ( outputGeom )
1146  {
1147  geom->fromGeos( outputGeom );
1148  }
1149  }
1150  else
1151  {
1152  GEOSGeometry** geomArray = new GEOSGeometry*[outputGeomList.size()];
1153  for ( int i = 0; i < outputGeomList.size(); ++i )
1154  {
1155  geomArray[i] = outputGeomList.at( i );
1156  }
1157  GEOSGeometry* collection = 0;
1158  if ( geom->type() == QGis::Point )
1159  {
1160  collection = GEOSGeom_createCollection_r( geosctxt, GEOS_MULTIPOINT, geomArray, outputGeomList.size() );
1161  }
1162  else if ( geom->type() == QGis::Line )
1163  {
1164  collection = GEOSGeom_createCollection_r( geosctxt, GEOS_MULTILINESTRING, geomArray, outputGeomList.size() );
1165  }
1166  geom->fromGeos( collection );
1167  delete[] geomArray;
1168  }
1169  return true;
1170 }
1171 
1172 QgsPoint QgsGeometryAnalyzer::createPointOffset( double x, double y, double dist, QgsGeometry* lineGeom ) const
1173 {
1174  QgsPoint p( x, y );
1175  QgsPoint minDistPoint;
1176  int afterVertexNr;
1177  lineGeom->closestSegmentWithContext( p, minDistPoint, afterVertexNr );
1178 
1179  int beforeVertexNr = afterVertexNr - 1;
1180  QgsPoint beforeVertex = lineGeom->vertexAt( beforeVertexNr );
1181  QgsPoint afterVertex = lineGeom->vertexAt( afterVertexNr );
1182 
1183  //get normal vector
1184  double dx = afterVertex.x() - beforeVertex.x();
1185  double dy = afterVertex.y() - beforeVertex.y();
1186  double normalX = -dy;
1187  double normalY = dx;
1188  double normalLength = sqrt( normalX * normalX + normalY * normalY );
1189  normalX *= ( dist / normalLength );
1190  normalY *= ( dist / normalLength );
1191 
1192  double debugLength = sqrt( normalX * normalX + normalY * normalY ); //control
1193  Q_UNUSED( debugLength );
1194  return QgsPoint( x - normalX, y - normalY ); //negative values -> left side, positive values -> right side
1195 }
1196 
1197 QgsGeometry* QgsGeometryAnalyzer::locateBetweenMeasures( double fromMeasure, double toMeasure, QgsGeometry* lineGeom )
1198 {
1199  if ( !lineGeom )
1200  {
1201  return 0;
1202  }
1203 
1204  QgsMultiPolyline resultGeom;
1205 
1206  //need to go with WKB and z coordinate until QgsGeometry supports M values
1207  const unsigned char* lineWkb = lineGeom->asWkb();
1208 
1209  const unsigned char* ptr = lineWkb + 1;
1210  QGis::WkbType wkbType;
1211  memcpy( &wkbType, ptr, sizeof( wkbType ) );
1212  ptr += sizeof( wkbType );
1213 
1214  if ( wkbType != QGis::WKBLineString25D && wkbType != QGis::WKBMultiLineString25D )
1215  {
1216  return 0;
1217  }
1218 
1219  if ( wkbType == QGis::WKBLineString25D )
1220  {
1221  locateBetweenWkbString( ptr, resultGeom, fromMeasure, toMeasure );
1222  }
1223  else if ( wkbType == QGis::WKBMultiLineString25D )
1224  {
1225  int* nLines = ( int* )ptr;
1226  ptr += sizeof( int );
1227  for ( int i = 0; i < *nLines; ++i )
1228  {
1229  ptr += ( 1 + sizeof( wkbType ) );
1230  ptr = locateBetweenWkbString( ptr, resultGeom, fromMeasure, toMeasure );
1231  }
1232  }
1233 
1234  if ( resultGeom.size() < 1 )
1235  {
1236  return 0;
1237  }
1238  return QgsGeometry::fromMultiPolyline( resultGeom );
1239 }
1240 
1242 {
1243  if ( !lineGeom )
1244  {
1245  return 0;
1246  }
1247 
1248  QgsMultiPoint resultGeom;
1249 
1250  //need to go with WKB and z coordinate until QgsGeometry supports M values
1251  const unsigned char* lineWkb = lineGeom->asWkb();
1252 
1253  const unsigned char* ptr = lineWkb + 1;
1254  QGis::WkbType wkbType;
1255  memcpy( &wkbType, ptr, sizeof( wkbType ) );
1256  ptr += sizeof( wkbType );
1257 
1258  if ( wkbType != QGis::WKBLineString25D && wkbType != QGis::WKBMultiLineString25D )
1259  {
1260  return 0;
1261  }
1262 
1263  if ( wkbType == QGis::WKBLineString25D )
1264  {
1265  locateAlongWkbString( ptr, resultGeom, measure );
1266  }
1267  else if ( wkbType == QGis::WKBMultiLineString25D )
1268  {
1269  int* nLines = ( int* )ptr;
1270  ptr += sizeof( int );
1271  for ( int i = 0; i < *nLines; ++i )
1272  {
1273  ptr += ( 1 + sizeof( wkbType ) );
1274  ptr = locateAlongWkbString( ptr, resultGeom, measure );
1275  }
1276  }
1277 
1278  if ( resultGeom.size() < 1 )
1279  {
1280  return 0;
1281  }
1282  return QgsGeometry::fromMultiPoint( resultGeom );
1283 }
1284 
1285 const unsigned char* QgsGeometryAnalyzer::locateBetweenWkbString( const unsigned char* ptr, QgsMultiPolyline& result, double fromMeasure, double toMeasure )
1286 {
1287  int* nPoints = ( int* ) ptr;
1288  ptr += sizeof( int );
1289  double prevx = 0.0, prevy = 0.0, prevz = 0.0;
1290  double *x, *y, *z;
1291  QgsPolyline currentLine;
1292 
1293  QgsPoint pt1, pt2;
1294  bool measureInSegment; //true if measure is contained in the segment
1295  bool secondPointClipped; //true if second point is != segment endpoint
1296 
1297 
1298  for ( int i = 0; i < *nPoints; ++i )
1299  {
1300  x = ( double* )ptr;
1301  ptr += sizeof( double );
1302  y = ( double* )ptr;
1303  ptr += sizeof( double );
1304  z = ( double* ) ptr;
1305  ptr += sizeof( double );
1306 
1307  if ( i > 0 )
1308  {
1309  measureInSegment = clipSegmentByRange( prevx, prevy, prevz, *x, *y, *z, fromMeasure, toMeasure, pt1, pt2, secondPointClipped );
1310  if ( measureInSegment )
1311  {
1312  if ( currentLine.size() < 1 ) //no points collected yet, so the first point needs to be added to the line
1313  {
1314  currentLine.append( pt1 );
1315  }
1316 
1317  if ( pt1 != pt2 ) //avoid duplicated entry if measure value equals m-value of vertex
1318  {
1319  currentLine.append( pt2 );
1320  }
1321 
1322  if ( secondPointClipped || i == *nPoints - 1 ) //close current segment
1323  {
1324  if ( currentLine.size() > 1 )
1325  {
1326  result.append( currentLine );
1327  }
1328  currentLine.clear();
1329  }
1330  }
1331  }
1332  prevx = *x; prevy = *y; prevz = *z;
1333  }
1334  return ptr;
1335 }
1336 
1337 const unsigned char* QgsGeometryAnalyzer::locateAlongWkbString( const unsigned char* ptr, QgsMultiPoint& result, double measure )
1338 {
1339  int* nPoints = ( int* ) ptr;
1340  ptr += sizeof( int );
1341  double prevx = 0.0, prevy = 0.0, prevz = 0.0;
1342  double *x, *y, *z;
1343 
1344  QgsPoint pt1, pt2;
1345  bool pt1Ok, pt2Ok;
1346 
1347  for ( int i = 0; i < *nPoints; ++i )
1348  {
1349  x = ( double* )ptr;
1350  ptr += sizeof( double );
1351  y = ( double* )ptr;
1352  ptr += sizeof( double );
1353  z = ( double* ) ptr;
1354  ptr += sizeof( double );
1355 
1356  if ( i > 0 )
1357  {
1358  locateAlongSegment( prevx, prevy, prevz, *x, *y, *z, measure, pt1Ok, pt1, pt2Ok, pt2 );
1359  if ( pt1Ok )
1360  {
1361  result.append( pt1 );
1362  }
1363  if ( pt2Ok && ( i == ( *nPoints - 1 ) ) )
1364  {
1365  result.append( pt2 );
1366  }
1367  }
1368  prevx = *x; prevy = *y; prevz = *z;
1369  }
1370  return ptr;
1371 }
1372 
1373 bool QgsGeometryAnalyzer::clipSegmentByRange( double x1, double y1, double m1, double x2, double y2, double m2, double range1, double range2, QgsPoint& pt1,
1374  QgsPoint& pt2, bool& secondPointClipped )
1375 {
1376  bool reversed = m1 > m2;
1377  double tmp;
1378 
1379  //reverse m1, m2 if necessary (and consequently also x1,x2 / y1, y2)
1380  if ( reversed )
1381  {
1382  tmp = m1;
1383  m1 = m2;
1384  m2 = tmp;
1385 
1386  tmp = x1;
1387  x1 = x2;
1388  x2 = tmp;
1389 
1390  tmp = y1;
1391  y1 = y2;
1392  y2 = tmp;
1393  }
1394 
1395  //reverse range1, range2 if necessary
1396  if ( range1 > range2 )
1397  {
1398  tmp = range1;
1399  range1 = range2;
1400  range2 = tmp;
1401  }
1402 
1403  //segment completely outside of range
1404  if ( m2 < range1 || m1 > range2 )
1405  {
1406  return false;
1407  }
1408 
1409  //segment completely inside of range
1410  if ( m2 <= range2 && m1 >= range1 )
1411  {
1412  if ( reversed )
1413  {
1414  pt1.setX( x2 ); pt1.setY( y2 );
1415  pt2.setX( x1 ); pt2.setY( y1 );
1416  }
1417  else
1418  {
1419  pt1.setX( x1 ); pt1.setY( y1 );
1420  pt2.setX( x2 ); pt2.setY( y2 );
1421  }
1422  secondPointClipped = false;
1423  return true;
1424  }
1425 
1426  //m1 inside and m2 not
1427  if ( m1 >= range1 && m1 <= range2 )
1428  {
1429  pt1.setX( x1 ); pt1.setY( y1 );
1430  double dist = ( range2 - m1 ) / ( m2 - m1 );
1431  pt2.setX( x1 + ( x2 - x1 ) * dist );
1432  pt2.setY( y1 + ( y2 - y1 ) * dist );
1433  secondPointClipped = !reversed;
1434  }
1435 
1436  //m2 inside and m1 not
1437  if ( m2 >= range1 && m2 <= range2 )
1438  {
1439  pt2.setX( x2 ); pt2.setY( y2 );
1440  double dist = ( m2 - range1 ) / ( m2 - m1 );
1441  pt1.setX( x2 - ( x2 - x1 ) * dist );
1442  pt1.setY( y2 - ( y2 - y1 ) * dist );
1443  secondPointClipped = reversed;
1444  }
1445 
1446  //range1 and range 2 both inside the segment
1447  if ( range1 >= m1 && range2 <= m2 )
1448  {
1449  double dist1 = ( range1 - m1 ) / ( m2 - m1 );
1450  double dist2 = ( range2 - m1 ) / ( m2 - m1 );
1451  pt1.setX( x1 + ( x2 - x1 ) * dist1 );
1452  pt1.setY( y1 + ( y2 - y1 ) * dist1 );
1453  pt2.setX( x1 + ( x2 - x1 ) * dist2 );
1454  pt2.setY( y1 + ( y2 - y1 ) * dist2 );
1455  secondPointClipped = true;
1456  }
1457 
1458  if ( reversed ) //switch p1 and p2
1459  {
1460  QgsPoint tmpPt = pt1;
1461  pt1 = pt2;
1462  pt2 = tmpPt;
1463  }
1464 
1465  return true;
1466 }
1467 
1468 void QgsGeometryAnalyzer::locateAlongSegment( double x1, double y1, double m1, double x2, double y2, double m2, double measure, bool& pt1Ok, QgsPoint& pt1, bool& pt2Ok, QgsPoint& pt2 )
1469 {
1470  bool reversed = false;
1471  pt1Ok = false;
1472  pt2Ok = false;
1473  double tolerance = 0.000001; //work with a small tolerance to catch e.g. locations at endpoints
1474 
1475  if ( m1 > m2 )
1476  {
1477  double tmp = m1;
1478  m1 = m2;
1479  m2 = tmp;
1480  reversed = true;
1481  }
1482 
1483  //segment does not match
1484  if (( m1 - measure ) > tolerance || ( measure - m2 ) > tolerance )
1485  {
1486  pt1Ok = false;
1487  pt2Ok = false;
1488  return;
1489  }
1490 
1491  //match with vertex1
1492  if ( qgsDoubleNear( m1, measure, tolerance ) )
1493  {
1494  if ( reversed )
1495  {
1496  pt2Ok = true;
1497  pt2.setX( x2 ); pt2.setY( y2 );
1498  }
1499  else
1500  {
1501  pt1Ok = true;
1502  pt1.setX( x1 ); pt1.setY( y1 );
1503  }
1504  }
1505 
1506  //match with vertex2
1507  if ( qgsDoubleNear( m2, measure, tolerance ) )
1508  {
1509  if ( reversed )
1510  {
1511  pt1Ok = true;
1512  pt1.setX( x1 ); pt1.setY( y1 );
1513  }
1514  else
1515  {
1516  pt2Ok = true;
1517  pt2.setX( x2 ); pt2.setY( y2 );
1518  }
1519  }
1520 
1521 
1522  if ( pt1Ok || pt2Ok )
1523  {
1524  return;
1525  }
1526 
1527  //match between the vertices
1528  if ( qgsDoubleNear( m1, m2 ) )
1529  {
1530  pt1.setX( x1 );
1531  pt1.setY( y1 );
1532  pt1Ok = true;
1533  return;
1534  }
1535  double dist = ( measure - m1 ) / ( m2 - m1 );
1536  if ( reversed )
1537  {
1538  dist = 1 - dist;
1539  }
1540 
1541  pt1.setX( x1 + dist * ( x2 - x1 ) );
1542  pt1.setY( y1 + dist * ( y2 - y1 ) );
1543  pt1Ok = true;
1544 }