SAGA API v9.10
Loading...
Searching...
No Matches
kdtree.cpp
Go to the documentation of this file.
1
3// //
4// SAGA //
5// //
6// System for Automated Geoscientific Analyses //
7// //
8// Application Programming Interface //
9// //
10// Library: SAGA_API //
11// //
12//-------------------------------------------------------//
13// //
14// kdtree.cpp //
15// //
16// Copyright (C) 2019 by Olaf Conrad //
17// //
18//-------------------------------------------------------//
19// //
20// This file is part of 'SAGA - System for Automated //
21// Geoscientific Analyses'. //
22// //
23// This library is free software; you can redistribute //
24// it and/or modify it under the terms of the GNU Lesser //
25// General Public License as published by the Free //
26// Software Foundation, either version 2.1 of the //
27// License, or (at your option) any later version. //
28// //
29// This library is distributed in the hope that it will //
30// be useful, but WITHOUT ANY WARRANTY; without even the //
31// implied warranty of MERCHANTABILITY or FITNESS FOR A //
32// PARTICULAR PURPOSE. See the GNU Lesser General Public //
33// License for more details. //
34// //
35// You should have received a copy of the GNU Lesser //
36// General Public License along with this program; if //
37// not, see <http://www.gnu.org/licenses/>. //
38// //
39//-------------------------------------------------------//
40// //
41// e-mail: oconrad@saga-gis.org //
42// //
43// contact: Olaf Conrad //
44// Institute of Geography //
45// University of Hamburg //
46// Germany //
47// //
49
50//---------------------------------------------------------
51#include "pointcloud.h"
52
53#include "nanoflann/nanoflann.hpp"
54
55
57// //
58// //
59// //
61
62//---------------------------------------------------------
63class CSG_KDTree_Adaptor
64{
65public:
66 CSG_KDTree_Adaptor(void) { m_pData = NULL; m_zScale = 1.; }
67 virtual ~CSG_KDTree_Adaptor(void) {}
68
69 typedef nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<double, CSG_KDTree_Adaptor>,
70 CSG_KDTree_Adaptor, 2> kd_tree_2d;
71
72 typedef nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<double, CSG_KDTree_Adaptor>,
73 CSG_KDTree_Adaptor, 3> kd_tree_3d;
74
75 //-----------------------------------------------------
76 virtual size_t kdtree_get_point_count (void) const = 0;
77 virtual double kdtree_get_pt (const size_t Index, int Dimension) const = 0;
78
79 template <class BBOX> bool kdtree_get_bbox (BBOX &bb) const
80 {
81 double Extent[3][2];
82
83 if( !Get_Extent(Extent) )
84 {
85 return( false );
86 }
87
88 bb[0].low = Extent[0][0]; bb[0].high = Extent[0][1];
89 bb[1].low = Extent[1][0]; bb[1].high = Extent[1][1];
90
91 if( bb.size() > 2 )
92 {
93 bb[2].low = Extent[2][0]; bb[2].high = Extent[2][1];
94 }
95
96 return( true );
97 }
98
99 //-----------------------------------------------------
100 CSG_Data_Object * Get_Data_Object (void) { return( m_pData ); }
101
102
103protected:
104
105 double m_zScale;
106
107 CSG_Data_Object *m_pData;
108
109
110 virtual bool Get_Extent (double Extent[3][2]) const { return( false ); }
111
112};
113
114//---------------------------------------------------------
115class CSG_KDTree_Adaptor_Points : public CSG_KDTree_Adaptor
116{
117public:
118 CSG_KDTree_Adaptor_Points(CSG_Shapes *pPoints, int zField = -1, double zScale = 1.)
119 {
120 m_pData = m_pPoints = pPoints;
121 m_zField = m_pPoints && zField < m_pPoints->Get_Count() ? zField : -1;
122 m_zScale = zScale;
123 }
124
125 virtual ~CSG_KDTree_Adaptor_Points(void) {}
126
127 virtual size_t kdtree_get_point_count (void) const
128 {
129 return( m_pPoints->Get_Count() );
130 }
131
132 virtual double kdtree_get_pt (const size_t Index, int Dimension) const
133 {
134 CSG_Shape *pPoint = m_pPoints->Get_Shape(Index);
135
136 if( Dimension == 0 ) { return( pPoint->Get_Point().x ); }
137 if( Dimension == 1 ) { return( pPoint->Get_Point().y ); }
138 if( Dimension == 2 ) { return( (m_zField < 0 ? pPoint->Get_Z() : pPoint->asDouble(m_zField)) * m_zScale); }
139
140 return( 0. );
141 }
142
143
144protected:
145
146 int m_zField;
147
148 CSG_Shapes *m_pPoints;
149
150
151 virtual bool Get_Extent (double Extent[3][2]) const
152 {
153 Extent[0][0] = m_pPoints->Get_Extent().Get_XMin();
154 Extent[0][1] = m_pPoints->Get_Extent().Get_XMax();
155 Extent[1][0] = m_pPoints->Get_Extent().Get_YMin();
156 Extent[1][1] = m_pPoints->Get_Extent().Get_YMax();
157 Extent[2][0] = m_zField < 0 ? m_pPoints->Get_ZMin() : m_pPoints->Get_Minimum(m_zField);
158 Extent[2][1] = m_zField < 0 ? m_pPoints->Get_ZMax() : m_pPoints->Get_Maximum(m_zField);
159
160 return( true );
161 }
162
163};
164
165//---------------------------------------------------------
166class CSG_KDTree_Adaptor_PointCloud : public CSG_KDTree_Adaptor
167{
168public:
169 CSG_KDTree_Adaptor_PointCloud(CSG_PointCloud *pPoints, double zScale = 1.)
170 {
171 m_pData = m_pPoints = pPoints;
172 m_zScale = zScale;
173 }
174
175 virtual ~CSG_KDTree_Adaptor_PointCloud(void) {}
176
177 virtual size_t kdtree_get_point_count (void) const
178 {
179 return( m_pPoints->Get_Count() );
180 }
181
182 virtual double kdtree_get_pt (const size_t Index, int Dimension) const
183 {
184 if( Dimension == 0 ) { return( m_pPoints->Get_X((sLong)Index) ); }
185 if( Dimension == 1 ) { return( m_pPoints->Get_Y((sLong)Index) ); }
186 if( Dimension == 2 ) { return( m_pPoints->Get_Z((sLong)Index) * m_zScale ); }
187
188 return( 0. );
189 }
190
191
192protected:
193
194 CSG_PointCloud *m_pPoints;
195
196
197 virtual bool Get_Extent (double Extent[3][2]) const
198 {
199 Extent[0][0] = m_pPoints->Get_Extent().Get_XMin();
200 Extent[0][1] = m_pPoints->Get_Extent().Get_XMax();
201 Extent[1][0] = m_pPoints->Get_Extent().Get_YMin();
202 Extent[1][1] = m_pPoints->Get_Extent().Get_YMax();
203 Extent[2][0] = m_pPoints-> Get_ZMin();
204 Extent[2][1] = m_pPoints-> Get_ZMax();
205
206 return( true );
207 }
208
209};
210
211//---------------------------------------------------------
212class CSG_KDTree_Adaptor_Coordinates : public CSG_KDTree_Adaptor
213{
214public:
215 CSG_KDTree_Adaptor_Coordinates(const double **Points, size_t nPoints)
216 {
217 m_Points = Points;
218 m_nPoints = nPoints;
219 }
220
221 virtual ~CSG_KDTree_Adaptor_Coordinates(void) {}
222
223 virtual size_t kdtree_get_point_count (void) const
224 {
225 return( m_nPoints );
226 }
227
228 virtual double kdtree_get_pt (const size_t Index, int Dimension) const
229 {
230 return( m_Points[Index][Dimension] );
231 }
232
233
234protected:
235
236 size_t m_nPoints;
237
238 const double **m_Points;
239
240};
241
242
244// //
245// //
246// //
248
249//---------------------------------------------------------
254
255//---------------------------------------------------------
258
259//---------------------------------------------------------
261{
262 m_pKDTree = NULL;
263 m_pAdaptor = NULL;
264}
265
266//---------------------------------------------------------
267const char * CSG_KDTree::Get_Version(void)
268{
269 static CSG_String Version(CSG_String::Format("nanoflann %d.%d.%d",
270 (NANOFLANN_VERSION&0xf00)/0x100,
271 (NANOFLANN_VERSION&0x0f0)/0x010,
272 (NANOFLANN_VERSION&0x00f)/0x001)
273 );
274
275 return( Version );
276}
277
278//---------------------------------------------------------
280{
282
283 m_Points .Destroy();
284 m_Indices .Destroy();
285 m_Distances.Destroy();
286
287 return( true );
288}
289
290//---------------------------------------------------------
292{
293 if( i < Get_Match_Count() )
294 {
295 CSG_Shapes *pShapes = m_pAdaptor && m_pAdaptor->Get_Data_Object() ? m_pAdaptor->Get_Data_Object()->asShapes() : NULL;
296
297 if( pShapes )
298 {
299 return( pShapes->Get_Shape(Get_Match_Index(i)) );
300 }
301 }
302
303 return( NULL );
304}
305
306
308// //
309// //
310// //
312
313//---------------------------------------------------------
317//---------------------------------------------------------
322
323//---------------------------------------------------------
328
329//---------------------------------------------------------
335//---------------------------------------------------------
337{
339
340 Create(pPoints, Field);
341}
342
343bool CSG_KDTree_2D::Create(CSG_Shapes *pPoints, int Field)
344{
345 Destroy();
346
347 //-----------------------------------------------------
348 if( Field >= 0 && Field < pPoints->Get_Field_Count() )
349 {
350 m_Points.Create(3, (int)pPoints->Get_Count());
351
352 int n = 0;
353
354 for(int i=0; i<pPoints->Get_Count(); i++)
355 {
356 CSG_Shape *pPoint = pPoints->Get_Shape(i);
357
358 if( !pPoint->is_NoData(Field) )
359 {
360 m_Points[n][0] = pPoint->Get_Point().x;
361 m_Points[n][1] = pPoint->Get_Point().y;
362 m_Points[n][2] = pPoint->asDouble(Field);
363
364 n++;
365 }
366 }
367
368 m_Points.Set_Rows(n); // resize if there are no-data values
369
370 if( n < 1 )
371 {
372 Destroy();
373
374 return( false );
375 }
376
377 m_pAdaptor = new CSG_KDTree_Adaptor_Coordinates(m_Points, m_Points.Get_NRows());
378 m_pKDTree = new CSG_KDTree_Adaptor::kd_tree_2d(2, *m_pAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(10));
379
380 ((CSG_KDTree_Adaptor::kd_tree_2d *)m_pKDTree)->buildIndex();
381
382 return( true );
383 }
384
385 //-----------------------------------------------------
386 if( pPoints->Get_Count() < 1 )
387 {
388 return( false );
389 }
390
391 m_pAdaptor = new CSG_KDTree_Adaptor_Points(pPoints);
392 m_pKDTree = new CSG_KDTree_Adaptor::kd_tree_2d(2, *m_pAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(10));
393
394 ((CSG_KDTree_Adaptor::kd_tree_2d *)m_pKDTree)->buildIndex();
395
396 return( true );
397}
398
399//---------------------------------------------------------
403//---------------------------------------------------------
405{
407
408 Create(pPoints);
409}
410
412{
413 if( pPoints->Get_Count() < 1 )
414 {
415 return( false );
416 }
417
418 Destroy();
419
420 m_pAdaptor = new CSG_KDTree_Adaptor_PointCloud(pPoints);
421 m_pKDTree = new CSG_KDTree_Adaptor::kd_tree_2d(2, *m_pAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(10));
422
423 ((CSG_KDTree_Adaptor::kd_tree_2d *)m_pKDTree)->buildIndex();
424
425 return( true );
426}
427
428//---------------------------------------------------------
433//---------------------------------------------------------
435{
437
438 Create(Points);
439}
440
442{
443 if( Points.Get_NCols() < 2 )
444 {
445 return( false );
446 }
447
448 return( Create((const double **)Points.Get_Data(), Points.Get_NRows()) );
449}
450
451//---------------------------------------------------------
456//---------------------------------------------------------
457CSG_KDTree_2D::CSG_KDTree_2D(const double **Points, size_t nPoints)
458{
460
461 Create(Points, nPoints);
462}
463
464bool CSG_KDTree_2D::Create(const double **Points, size_t nPoints)
465{
466 if( nPoints < 1 )
467 {
468 return( false );
469 }
470
471 Destroy();
472
473 m_pAdaptor = new CSG_KDTree_Adaptor_Coordinates(Points, nPoints);
474 m_pKDTree = new CSG_KDTree_Adaptor::kd_tree_2d(2, *m_pAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(10));
475
476 ((CSG_KDTree_Adaptor::kd_tree_2d *)m_pKDTree)->buildIndex();
477
478 return( true );
479}
480
481//---------------------------------------------------------
483{
484 if( m_pKDTree )
485 {
486 delete((CSG_KDTree_Adaptor::kd_tree_2d *)m_pKDTree);
487
488 m_pKDTree = NULL;
489 }
490
491 return( CSG_KDTree::Destroy() );
492}
493
494
496// //
498
499//---------------------------------------------------------
500size_t CSG_KDTree_2D::Get_Nearest_Points(double Coordinate[2], size_t Count, double Radius)
501{
502 return( Get_Nearest_Points(Coordinate, Count, Radius, m_Indices, m_Distances) );
503}
504
505//---------------------------------------------------------
506size_t CSG_KDTree_2D::Get_Nearest_Points(double Coordinate[2], size_t Count, double Radius, CSG_Array_sLong &Indices, CSG_Vector &Distances)
507{
508 if( Radius > 0. )
509 {
510 nanoflann::SearchParams SearchParams;
511
512 SearchParams.sorted = Count > 0;
513
514 std::vector<std::pair<size_t, double>> Matches;
515
516 ((CSG_KDTree_Adaptor::kd_tree_2d *)m_pKDTree)->radiusSearch(Coordinate, Radius*Radius, Matches, SearchParams);
517
518 if( Count == 0 || Count > Matches.size() )
519 {
520 Count = Matches.size();
521 }
522
523 Indices .Create(Count);
524 Distances.Create(Count);
525
526 for(size_t i=0; i<Count; i++)
527 {
528 Indices [i] = (int)Matches[i]. first ;
529 Distances[i] = sqrt(Matches[i].second);
530 }
531 }
532 else if( Count > 0 )
533 {
534 size_t *_Indices = new size_t[Count];
535
536 Distances.Create(Count);
537
538 Count = Get_Nearest_Points(Coordinate, Count, _Indices, Distances.Get_Data());
539
540 if( Count < (size_t)Distances.Get_N() )
541 {
542 Distances.Set_Rows(Count);
543 }
544
545 Indices.Create(Count);
546
547 for(size_t i=0; i<Count; i++)
548 {
549 Indices[i] = (int)(_Indices[i]);
550 }
551
552 delete[](_Indices);
553 }
554
555 return( Count );
556}
557
558//---------------------------------------------------------
559size_t CSG_KDTree_2D::Get_Nearest_Points(double Coordinate[2], size_t Count, size_t *Indices, double *Distances)
560{
561 Count = ((CSG_KDTree_Adaptor::kd_tree_2d *)m_pKDTree)->knnSearch(Coordinate, Count, Indices, Distances);
562
563 for(size_t i=0; i<Count; i++)
564 {
565 Distances[i] = sqrt(Distances[i]);
566 }
567
568 return( Count );
569}
570
571//---------------------------------------------------------
572bool CSG_KDTree_2D::Get_Nearest_Point(double Coordinate[2], size_t &Index, double &Distance)
573{
574 return( Get_Nearest_Points(Coordinate, 1, &Index, &Distance) == 1 );
575}
576
577//---------------------------------------------------------
578bool CSG_KDTree_2D::Get_Nearest_Point(double Coordinate[2], size_t &Index)
579{
580 double Distance;
581
582 return( Get_Nearest_Points(Coordinate, 1, &Index, &Distance) == 1 );
583}
584
585//---------------------------------------------------------
586bool CSG_KDTree_2D::Get_Nearest_Value(double Coordinate[2], double &Value)
587{
588 size_t Index; double Distance;
589
590 if( Get_Nearest_Points(Coordinate, 1, &Index, &Distance) == 1 )
591 {
592 Value = m_Points.Get_Data() ? Get_Point_Value(Index) : (double)Index;
593
594 return( true );
595 }
596
597 return( false );
598}
599
600//---------------------------------------------------------
602{
603 size_t Index; CSG_Shapes *pShapes = m_pAdaptor && m_pAdaptor->Get_Data_Object() ? m_pAdaptor->Get_Data_Object()->asShapes() : NULL;
604
605 return( pShapes && Get_Nearest_Point(Coordinate, Index) ? pShapes->Get_Shape((sLong)Index) : NULL );
606}
607
608//---------------------------------------------------------
609size_t CSG_KDTree_2D::Get_Duplicates(double Coordinate[2], CSG_Array_sLong &Indices, CSG_Vector &Distances)
610{
611 nanoflann::SearchParams SearchParams;
612
613 SearchParams.sorted = false;
614
615 std::vector<std::pair<size_t, double>> Matches;
616
617 ((CSG_KDTree_Adaptor::kd_tree_2d *)m_pKDTree)->radiusSearch(Coordinate, 0.0000001, Matches, SearchParams);
618
619 Indices .Create(Matches.size());
620 Distances.Create(Matches.size());
621
622 size_t Count = 0;
623
624 for(size_t i=0; i<Matches.size(); i++)
625 {
626 if( Matches[i].second > 0. )
627 {
628 Indices.Dec_Array();
629 }
630 else
631 {
632 Indices[Count++] = (int)Matches[i].first;
633 }
634 }
635
636 return( Count );
637}
638
639//---------------------------------------------------------
640size_t CSG_KDTree_2D::Get_Duplicates(double Coordinate[2])
641{
642 return( Get_Duplicates(Coordinate, m_Indices, m_Distances) );
643}
644
645//---------------------------------------------------------
646size_t CSG_KDTree_2D::Get_Nearest_Points(double x, double y, size_t Count, double Radius)
647{
648 double c[2]; c[0] = x; c[1] = y; return( Get_Nearest_Points(c, Count, Radius) );
649}
650
651size_t CSG_KDTree_2D::Get_Nearest_Points(double x, double y, size_t Count, double Radius, CSG_Array_sLong &Indices, CSG_Vector &Distances)
652{
653 double c[2]; c[0] = x; c[1] = y; return( Get_Nearest_Points(c, Count, Radius, Indices, Distances) );
654}
655
656size_t CSG_KDTree_2D::Get_Nearest_Points(double x, double y, size_t Count, size_t *Indices, double *Distances)
657{
658 double c[2]; c[0] = x; c[1] = y; return( Get_Nearest_Points(c, Count, Indices, Distances) );
659}
660
661bool CSG_KDTree_2D::Get_Nearest_Point(double x, double y, size_t &Index, double &Distance)
662{
663 double c[2]; c[0] = x; c[1] = y; return( Get_Nearest_Point(c, Index, Distance) );
664}
665
666bool CSG_KDTree_2D::Get_Nearest_Point(double x, double y, size_t &Index)
667{
668 double c[2]; c[0] = x; c[1] = y; return( Get_Nearest_Point(c, Index) );
669}
670
671bool CSG_KDTree_2D::Get_Nearest_Value(double x, double y, double &Value)
672{
673 double c[2]; c[0] = x; c[1] = y; return( Get_Nearest_Value(c, Value) );
674}
675
677{
678 double c[2]; c[0] = x; c[1] = y; return( Get_Nearest_Shape(c) );
679}
680
681size_t CSG_KDTree_2D::Get_Duplicates(double x, double y, CSG_Array_sLong &Indices, CSG_Vector &Distances)
682{
683 double c[2]; c[0] = x; c[1] = y; return( Get_Duplicates(c, Indices, Distances) );
684}
685
686size_t CSG_KDTree_2D::Get_Duplicates(double x, double y)
687{
688 double c[2]; c[0] = x; c[1] = y; return( Get_Duplicates(c) );
689}
690
691
693// //
694// //
695// //
697
698//---------------------------------------------------------
703
704//---------------------------------------------------------
709
710//---------------------------------------------------------
719//---------------------------------------------------------
720CSG_KDTree_3D::CSG_KDTree_3D(CSG_Shapes *pPoints, int Field, int zField, double zScale)
721{
723
724 Create(pPoints, Field, zField, zScale);
725}
726
727bool CSG_KDTree_3D::Create(CSG_Shapes *pPoints, int Field, int zField, double zScale)
728{
729 Destroy();
730
731 //-----------------------------------------------------
732 if( Field >= 0 && Field < pPoints->Get_Field_Count() )
733 {
734 m_Points.Create(4, (int)pPoints->Get_Count());
735
736 int n = 0;
737
738 for(int i=0; i<pPoints->Get_Count(); i++)
739 {
740 CSG_Shape *pPoint = pPoints->Get_Shape(i);
741
742 if( !pPoint->is_NoData(Field) )
743 {
744 m_Points[n][0] = pPoint->Get_Point().x;
745 m_Points[n][1] = pPoint->Get_Point().y;
746 m_Points[n][2] = zScale * (zField < 0 ? pPoint->Get_Z() : pPoint->asDouble(zField));
747 m_Points[n][3] = pPoint->asDouble(Field);
748
749 n++;
750 }
751 }
752
753 m_Points.Set_Rows(n); // resize if there are no-data values
754
755 if( n < 1 )
756 {
757 Destroy();
758
759 return( false );
760 }
761
762 m_pAdaptor = new CSG_KDTree_Adaptor_Coordinates(m_Points, m_Points.Get_NRows());
763 m_pKDTree = new CSG_KDTree_Adaptor::kd_tree_3d(3, *m_pAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(10));
764
765 ((CSG_KDTree_Adaptor::kd_tree_3d *)m_pKDTree)->buildIndex();
766
767 return( true );
768 }
769
770 //-----------------------------------------------------
771 if( pPoints->Get_Count() < 1 )
772 {
773 return( false );
774 }
775
776 m_pAdaptor = new CSG_KDTree_Adaptor_Points(pPoints, zField, zScale);
777 m_pKDTree = new CSG_KDTree_Adaptor::kd_tree_3d(3, *m_pAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(10));
778
779 ((CSG_KDTree_Adaptor::kd_tree_3d *)m_pKDTree)->buildIndex();
780
781 return( true );
782}
783
784//---------------------------------------------------------
786{
788
789 Create(pPoints);
790}
791
793{
794 if( pPoints->Get_Count() < 1 )
795 {
796 return( false );
797 }
798
799 Destroy();
800
801 m_pAdaptor = new CSG_KDTree_Adaptor_PointCloud(pPoints);
802 m_pKDTree = new CSG_KDTree_Adaptor::kd_tree_3d(3, *m_pAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(10));
803
804 ((CSG_KDTree_Adaptor::kd_tree_3d *)m_pKDTree)->buildIndex();
805
806 return( true );
807}
808
809//---------------------------------------------------------
814//---------------------------------------------------------
816{
818
819 Create(Points);
820}
821
823{
824 if( Points.Get_NCols() < 3 )
825 {
826 return( false );
827 }
828
829 return( Create((const double **)Points.Get_Data(), Points.Get_NRows()) );
830}
831
832//---------------------------------------------------------
837//---------------------------------------------------------
838CSG_KDTree_3D::CSG_KDTree_3D(const double **Points, size_t nPoints)
839{
841
842 Create(Points, nPoints);
843}
844
845bool CSG_KDTree_3D::Create(const double **Points, size_t nPoints)
846{
847 if( nPoints < 1 )
848 {
849 return( false );
850 }
851
852 Destroy();
853
854 m_pAdaptor = new CSG_KDTree_Adaptor_Coordinates(Points, nPoints);
855 m_pKDTree = new CSG_KDTree_Adaptor::kd_tree_3d(3, *m_pAdaptor, nanoflann::KDTreeSingleIndexAdaptorParams(10));
856
857 ((CSG_KDTree_Adaptor::kd_tree_3d *)m_pKDTree)->buildIndex();
858
859 return( true );
860}
861
862//---------------------------------------------------------
864{
865 if( m_pKDTree )
866 {
867 delete((CSG_KDTree_Adaptor::kd_tree_3d *)m_pKDTree);
868
869 m_pKDTree = NULL;
870 }
871
872 return( CSG_KDTree::Destroy() );
873}
874
875
877// //
879
880//---------------------------------------------------------
881size_t CSG_KDTree_3D::Get_Nearest_Points(double Coordinate[3], size_t Count, double Radius)
882{
883 return( Get_Nearest_Points(Coordinate, Count, Radius, m_Indices, m_Distances) );
884}
885
886//---------------------------------------------------------
887size_t CSG_KDTree_3D::Get_Nearest_Points(double Coordinate[3], size_t Count, double Radius, CSG_Array_sLong &Indices, CSG_Vector &Distances)
888{
889 if( Radius > 0. )
890 {
891 nanoflann::SearchParams SearchParams;
892
893 SearchParams.sorted = Count > 0;
894
895 std::vector<std::pair<size_t, double>> Matches;
896
897 ((CSG_KDTree_Adaptor::kd_tree_3d *)m_pKDTree)->radiusSearch(Coordinate, Radius*Radius, Matches, SearchParams);
898
899 if( Count == 0 || Count > Matches.size() )
900 {
901 Count = Matches.size();
902 }
903
904 Indices .Create(Count);
905 Distances.Create(Count);
906
907 for(size_t i=0; i<Count; i++)
908 {
909 Indices [i] = (int)Matches[i]. first ;
910 Distances[i] = sqrt(Matches[i].second);
911 }
912 }
913 else if( Count > 0 )
914 {
915 size_t *_Indices = new size_t[Count];
916
917 Distances.Create(Count);
918
919 Count = Get_Nearest_Points(Coordinate, Count, _Indices, Distances.Get_Data());
920
921 if( Count < (size_t)Distances.Get_N() )
922 {
923 Distances.Set_Rows(Count);
924 }
925
926 Indices.Create(Count);
927
928 for(size_t i=0; i<Count; i++)
929 {
930 Indices[i] = (int)(_Indices[i]);
931 }
932
933 delete[](_Indices);
934 }
935
936 return( Count );
937}
938
939//---------------------------------------------------------
940size_t CSG_KDTree_3D::Get_Nearest_Points(double Coordinate[3], size_t Count, size_t *Indices, double *Distances)
941{
942 Count = ((CSG_KDTree_Adaptor::kd_tree_3d *)m_pKDTree)->knnSearch(Coordinate, Count, Indices, Distances);
943
944 for(size_t i=0; i<Count; i++)
945 {
946 Distances[i] = sqrt(Distances[i]);
947 }
948
949 return( Count );
950}
951
952//---------------------------------------------------------
953bool CSG_KDTree_3D::Get_Nearest_Point(double Coordinate[3], size_t &Index, double &Distance)
954{
955 return( Get_Nearest_Points(Coordinate, 1, &Index, &Distance) == 1 );
956}
957
958//---------------------------------------------------------
959bool CSG_KDTree_3D::Get_Nearest_Point(double Coordinate[3], size_t &Index)
960{
961 double Distance;
962
963 return( Get_Nearest_Points(Coordinate, 1, &Index, &Distance) == 1 );
964}
965
966//---------------------------------------------------------
967bool CSG_KDTree_3D::Get_Nearest_Value(double Coordinate[3], double &Value)
968{
969 size_t Index; double Distance;
970
971 if( Get_Nearest_Points(Coordinate, 1, &Index, &Distance) == 1 )
972 {
973 Value = m_Points.Get_Data() ? Get_Point_Value(Index) : (double)Index;
974
975 return( true );
976 }
977
978 return( false );
979}
980
981//---------------------------------------------------------
983{
984 size_t Index; CSG_Shapes *pShapes = m_pAdaptor && m_pAdaptor->Get_Data_Object() ? m_pAdaptor->Get_Data_Object()->asShapes() : NULL;
985
986 return( pShapes && Get_Nearest_Point(Coordinate, Index) ? pShapes->Get_Shape((sLong)Index) : NULL );
987}
988
989//---------------------------------------------------------
990size_t CSG_KDTree_3D::Get_Duplicates(double Coordinate[3], CSG_Array_sLong &Indices, CSG_Vector &Distances)
991{
992 nanoflann::SearchParams SearchParams;
993
994 SearchParams.sorted = false;
995
996 std::vector<std::pair<size_t, double>> Matches;
997
998 ((CSG_KDTree_Adaptor::kd_tree_3d *)m_pKDTree)->radiusSearch(Coordinate, 0.0000001, Matches, SearchParams);
999
1000 Indices .Create(Matches.size());
1001 Distances.Create(Matches.size());
1002
1003 size_t Count = 0;
1004
1005 for(size_t i=0; i<Matches.size(); i++)
1006 {
1007 if( Matches[i].second > 0. )
1008 {
1009 Indices.Dec_Array();
1010 }
1011 else
1012 {
1013 Indices[Count++] = (int)Matches[i].first;
1014 }
1015 }
1016
1017 return( Count );
1018}
1019
1020//---------------------------------------------------------
1021size_t CSG_KDTree_3D::Get_Duplicates(double Coordinate[3])
1022{
1023 return( Get_Duplicates(Coordinate, m_Indices, m_Distances) );
1024}
1025
1026//---------------------------------------------------------
1027size_t CSG_KDTree_3D::Get_Nearest_Points(double x, double y, double z, size_t Count, double Radius)
1028{
1029 double c[3]; c[0] = x; c[1] = y; c[2] = z; return( Get_Nearest_Points(c, Count, Radius) );
1030}
1031
1032size_t CSG_KDTree_3D::Get_Nearest_Points(double x, double y, double z, size_t Count, double Radius, CSG_Array_sLong &Indices, CSG_Vector &Distances)
1033{
1034 double c[3]; c[0] = x; c[1] = y; c[2] = z; return( Get_Nearest_Points(c, Count, Radius, Indices, Distances) );
1035}
1036
1037size_t CSG_KDTree_3D::Get_Nearest_Points(double x, double y, double z, size_t Count, size_t *Indices, double *Distances)
1038{
1039 double c[3]; c[0] = x; c[1] = y; c[2] = z; return( Get_Nearest_Points(c, Count, Indices, Distances) );
1040}
1041
1042bool CSG_KDTree_3D::Get_Nearest_Point(double x, double y, double z, size_t &Index, double &Distance)
1043{
1044 double c[3]; c[0] = x; c[1] = y; c[2] = z; return( Get_Nearest_Point(c, Index, Distance) );
1045}
1046
1047bool CSG_KDTree_3D::Get_Nearest_Point(double x, double y, double z, size_t &Index)
1048{
1049 double c[3]; c[0] = x; c[1] = y; c[2] = z; return( Get_Nearest_Point(c, Index) );
1050}
1051
1052bool CSG_KDTree_3D::Get_Nearest_Value(double x, double y, double z, double &Value)
1053{
1054 double c[3]; c[0] = x; c[1] = y; c[2] = z; return( Get_Nearest_Value(c, Value) );
1055}
1056
1057CSG_Shape * CSG_KDTree_3D::Get_Nearest_Shape(double x, double y, double z)
1058{
1059 double c[3]; c[0] = x; c[1] = y; c[2] = z; return( Get_Nearest_Shape(c) );
1060}
1061
1062size_t CSG_KDTree_3D::Get_Duplicates(double x, double y, double z, CSG_Array_sLong &Indices, CSG_Vector &Distances)
1063{
1064 double c[3]; c[0] = x; c[1] = y; c[2] = z; return( Get_Duplicates(c, Indices, Distances) );
1065}
1066
1067size_t CSG_KDTree_3D::Get_Duplicates(double x, double y, double z)
1068{
1069 double c[3]; c[0] = x; c[1] = y; c[2] = z; return( Get_Duplicates(c) );
1070}
1071
1072
1074// //
1075// //
1076// //
1078
1079//---------------------------------------------------------
1080#include "parameters.h"
1081
1082//---------------------------------------------------------
1084{
1085 m_pParameters = NULL;
1086
1087 m_minPoints = 1;
1088 m_maxPoints = 0;
1089 m_Radius = 0.;
1090}
1091
1092//---------------------------------------------------------
1093bool CSG_Parameters_Point_Search::Create(CSG_Parameters *pParameters, const CSG_String &Parent, size_t minPoints)
1094{
1095 if( pParameters == NULL || m_pParameters != NULL )
1096 {
1097 return( false );
1098 }
1099
1100 m_pParameters = pParameters;
1101
1102 if( !Parent.is_Empty() && !(*m_pParameters)(Parent) )
1103 {
1104 m_pParameters->Add_Node("", Parent, _TL("Search Options"), _TL(""));
1105 }
1106
1107 //-----------------------------------------------------
1108 m_pParameters->Add_Choice(Parent,
1109 "SEARCH_RANGE" , _TL("Search Range"),
1110 _TL(""),
1111 CSG_String::Format("%s|%s",
1112 _TL("local"),
1113 _TL("global")
1114 ), 1
1115 );
1116
1117 m_pParameters->Add_Double("SEARCH_RANGE",
1118 "SEARCH_RADIUS" , _TL("Maximum Search Distance"),
1119 _TL("local maximum search distance given in map units"),
1120 1000., 0., true
1121 );
1122
1123 m_pParameters->Add_Choice(Parent,
1124 "SEARCH_POINTS_ALL" , _TL("Number of Points"),
1125 _TL(""),
1126 CSG_String::Format("%s|%s",
1127 _TL("maximum number of nearest points"),
1128 _TL("all points within search distance")
1129 ), 1
1130 );
1131
1132 if( minPoints > 0 )
1133 {
1134 m_pParameters->Add_Int("SEARCH_POINTS_ALL",
1135 "SEARCH_POINTS_MIN" , _TL("Minimum"),
1136 _TL("minimum number of points to use"),
1137 (int)minPoints, 1, true
1138 );
1139 }
1140
1141 m_pParameters->Add_Int("SEARCH_POINTS_ALL",
1142 "SEARCH_POINTS_MAX" , _TL("Maximum"),
1143 _TL("maximum number of nearest points"),
1144 20, 1, true
1145 );
1146
1147 return( true );
1148}
1149
1150
1152// //
1154
1155//---------------------------------------------------------
1157{
1158 if( !m_pParameters || !pParameters || m_pParameters->Get_Identifier().Cmp(pParameters->Get_Identifier()) || !pParameter || !pParameter->asShapes() )
1159 {
1160 return( false );
1161 }
1162
1163 pParameters->Set_Parameter("SEARCH_RADIUS", SG_Get_Rounded_To_SignificantFigures(
1164 5 * sqrt(pParameter->asShapes()->Get_Extent().Get_Area() / pParameter->asShapes()->Get_Count()), 1
1165 ));
1166
1167 return( true );
1168}
1169
1170//---------------------------------------------------------
1172{
1173 if( !m_pParameters || !pParameters || m_pParameters->Get_Identifier().Cmp(pParameters->Get_Identifier()) || !pParameter )
1174 {
1175 return( false );
1176 }
1177
1178 if( pParameter->Cmp_Identifier("SEARCH_RANGE") )
1179 {
1180 pParameters->Set_Enabled("SEARCH_RADIUS" , pParameter->asInt() == 0); // local
1181 pParameters->Set_Enabled("SEARCH_POINTS_MIN", pParameter->asInt() == 0); // when global, no minimum number of points
1182 }
1183
1184 if( pParameter->Cmp_Identifier("SEARCH_POINTS_ALL") )
1185 {
1186 pParameters->Set_Enabled("SEARCH_POINTS_MAX", pParameter->asInt() == 0); // maximum number of points
1187 pParameters->Set_Enabled("SEARCH_DIRECTION" , pParameter->asInt() == 0); // maximum number of points per quadrant
1188 }
1189
1190 return( true );
1191}
1192
1193
1195// //
1197
1198//---------------------------------------------------------
1200{
1201 if( m_pParameters )
1202 {
1203 m_minPoints = (*m_pParameters)("SEARCH_POINTS_MIN")
1204 ? (*m_pParameters)("SEARCH_POINTS_MIN")->asInt () : 0;
1205 m_maxPoints = (*m_pParameters)("SEARCH_POINTS_ALL")->asInt () == 0
1206 ? (*m_pParameters)("SEARCH_POINTS_MAX")->asInt () : 0;
1207 m_Radius = (*m_pParameters)("SEARCH_RANGE" )->asInt () == 0
1208 ? (*m_pParameters)("SEARCH_RADIUS" )->asDouble() : 0.;
1209
1210 return( true );
1211 }
1212
1213 return( false );
1214}
1215
1216//---------------------------------------------------------
1218{
1219 if( bUpdate )
1220 {
1221 Update();
1222 }
1223
1224 return( m_maxPoints == 0 && m_Radius <= 0. );
1225}
1226
1227
1229// //
1230// //
1231// //
1233
1234//---------------------------------------------------------
1239
1240
1242// //
1244
1245//---------------------------------------------------------
1247{
1248 if( !Finalize() || !m_pParameters || !pPoints || pPoints->Get_Count() < 1 || !Update() )
1249 {
1250 return( false );
1251 }
1252
1253 if( Do_Use_All() )
1254 {
1255 m_pPoints = pPoints;
1256 m_zField = zField;
1257
1258 return( true );
1259 }
1260
1261 return( m_Search.Create(pPoints, m_zField = zField) );
1262}
1263
1264//---------------------------------------------------------
1266{
1267 m_pPoints = NULL;
1268 m_zField = -1;
1269
1270 m_Search.Destroy();
1271
1272 return( true );
1273}
1274
1275
1277// //
1279
1280//---------------------------------------------------------
1282{
1283 if( m_pPoints ) // without search engine
1284 {
1285 return( m_pPoints->Get_Count() );
1286 }
1287 else // using search engine
1288 {
1289 return( m_Search.Get_Nearest_Points(x, y, m_maxPoints, m_Radius) );
1290 }
1291}
1292
1293//---------------------------------------------------------
1298
1299
1301// //
1303
1304//---------------------------------------------------------
1305bool CSG_Parameters_Point_Search_KDTree_2D::Get_Point(sLong Index, double &x, double &y, double &z)
1306{
1307 if( m_pPoints ) // without search engine
1308 {
1309 CSG_Shape *pPoint = m_pPoints->Get_Shape(Index);
1310
1311 if( pPoint && !pPoint->is_NoData(m_zField) )
1312 {
1313 x = pPoint->Get_Point().x;
1314 y = pPoint->Get_Point().y;
1315
1316 z = m_zField < 0 ? (double)Index : pPoint->asDouble(m_zField);
1317
1318 return( true );
1319 }
1320 }
1321 else // using search engine
1322 {
1323 if( Index >= 0 && (size_t)Index < m_Search.Get_Match_Count() )
1324 {
1325 if( m_zField >= 0 )
1326 {
1327 double *p = m_Search.Get_Point((size_t)Index);
1328
1329 x = p[0]; y = p[1]; z = p[2];
1330
1331 return( true );
1332 }
1333 else
1334 {
1335 CSG_Shape *pPoint = m_Search.Get_Match_Shape((size_t)Index);
1336
1337 if( pPoint )
1338 {
1339 x = pPoint->Get_Point().x;
1340 y = pPoint->Get_Point().y;
1341
1342 z = (double)Index;
1343
1344 return( true );
1345 }
1346 }
1347 }
1348 }
1349
1350 return( false );
1351}
1352
1353
1355// //
1357
1358//---------------------------------------------------------
1360{
1361 CSG_Array_sLong Index; if( Get_Points(x, y, Index) )
1362 {
1363 Points.Clear();
1364
1365 for(sLong i=0; i<Index.Get_Size(); i++)
1366 {
1367 if( m_zField >= 0 )
1368 {
1369 double *p = m_Search.Get_Point((size_t)Index[i]);
1370
1371 Points.Add(p[0], p[1], p[2]);
1372 }
1373 else
1374 {
1375 CSG_Shape *pPoint = m_Search.Get_Match_Shape((size_t)Index[i]);
1376
1377 if( pPoint )
1378 {
1379 Points.Add(pPoint->Get_Point().x, pPoint->Get_Point().y, (double)Index[i]);
1380 }
1381 }
1382 }
1383
1384 return( true );
1385 }
1386
1387 return( false );
1388}
1389
1390//---------------------------------------------------------
1392{
1393 return( Get_Points(p.x, p.y, Points) );
1394}
1395
1396
1398// //
1400
1401//---------------------------------------------------------
1403{
1404 if( m_pPoints ) // without search engine
1405 {
1406 return( m_pPoints->Get_Count() );
1407 }
1408 else // using search engine
1409 {
1410 CSG_Vector Distances;
1411
1412 return( m_Search.Get_Nearest_Points(x, y, m_maxPoints, m_Radius, Indices, Distances) >= m_minPoints );
1413 }
1414}
1415
1416//---------------------------------------------------------
1418{
1419 return( Get_Points(p.x, p.y, Indices) );
1420}
1421
1422//---------------------------------------------------------
1424{
1425 if( m_pPoints ) // without search engine
1426 {
1427 return( m_pPoints->Get_Count() );
1428 }
1429 else // using search engine
1430 {
1431 return( m_Search.Get_Nearest_Points(x, y, m_maxPoints, m_Radius, Indices, Distances) >= m_minPoints );
1432 }
1433}
1434
1435//---------------------------------------------------------
1437{
1438 return( Get_Points(p.x, p.y, Indices, Distances) );
1439}
1440
1441
1443// //
1444// //
1445// //
1447
1448//---------------------------------------------------------
#define SG_DELETE_SAFE(PTR)
Definition api_core.h:206
signed long long sLong
Definition api_core.h:158
#define _TL(s)
Definition api_core.h:1568
sLong Get_Size(void) const
Definition api_core.h:492
bool Dec_Array(bool bShrink=true)
Definition api_core.h:500
sLong * Create(const CSG_Array_sLong &Array)
virtual bool Destroy(void)
Definition kdtree.cpp:482
virtual size_t Get_Nearest_Points(double Coordinate[2], size_t Count, double Radius)
Definition kdtree.cpp:500
bool Create(CSG_Shapes *pPoints, int Field=-1)
Definition kdtree.cpp:343
virtual size_t Get_Duplicates(double Coordinate[2], CSG_Array_sLong &Indices, CSG_Vector &Distances)
Definition kdtree.cpp:609
virtual double Get_Point_Value(size_t i) const
Definition shapes.h:1280
virtual ~CSG_KDTree_2D(void)
Definition kdtree.cpp:324
virtual bool Get_Nearest_Point(double Coordinate[2], size_t &Index, double &Distance)
Definition kdtree.cpp:572
virtual bool Get_Nearest_Value(double Coordinate[2], double &Value)
Definition kdtree.cpp:586
CSG_KDTree_2D(void)
Definition kdtree.cpp:318
virtual CSG_Shape * Get_Nearest_Shape(double Coordinate[2])
Definition kdtree.cpp:601
virtual bool Destroy(void)
Definition kdtree.cpp:863
virtual ~CSG_KDTree_3D(void)
Definition kdtree.cpp:705
virtual double Get_Point_Value(size_t i) const
Definition shapes.h:1330
virtual CSG_Shape * Get_Nearest_Shape(double Coordinate[3])
Definition kdtree.cpp:982
virtual size_t Get_Duplicates(double Coordinate[3], CSG_Array_sLong &Indices, CSG_Vector &Distances)
Definition kdtree.cpp:990
virtual bool Get_Nearest_Point(double Coordinate[3], size_t &Index, double &Distance)
Definition kdtree.cpp:953
bool Create(CSG_Shapes *pPoints, int Field=-1, int zField=-1, double zScale=1.)
Definition kdtree.cpp:727
virtual bool Get_Nearest_Value(double Coordinate[3], double &Value)
Definition kdtree.cpp:967
CSG_KDTree_3D(void)
Definition kdtree.cpp:699
virtual size_t Get_Nearest_Points(double Coordinate[3], size_t Count, double Radius)
Definition kdtree.cpp:881
CSG_KDTree(void)
Definition kdtree.cpp:250
CSG_Array_sLong m_Indices
Definition shapes.h:1243
size_t Get_Match_Count(void) const
Definition shapes.h:1223
CSG_Matrix m_Points
Definition shapes.h:1247
void _On_Construction(void)
Definition kdtree.cpp:260
void * m_pKDTree
Definition shapes.h:1241
class CSG_KDTree_Adaptor * m_pAdaptor
Definition shapes.h:1239
virtual ~CSG_KDTree(void)
Definition kdtree.cpp:256
CSG_Shape * Get_Match_Shape(size_t i) const
Definition kdtree.cpp:291
CSG_Vector m_Distances
Definition shapes.h:1245
virtual bool Destroy(void)
Definition kdtree.cpp:279
size_t Get_Match_Index(size_t i) const
Definition shapes.h:1224
static const char * Get_Version(void)
Definition kdtree.cpp:267
sLong Get_NRows(void) const
Definition mat_tools.h:525
double ** Get_Data(void) const
Definition mat_tools.h:528
sLong Get_NCols(void) const
Definition mat_tools.h:524
CSG_Shapes * asShapes(void) const
bool Cmp_Identifier(const CSG_String &Identifier) const
int asInt(void) const
Definition parameters.h:287
bool Get_Points(double x, double y, CSG_Points_3D &Points)
Definition kdtree.cpp:1359
sLong Set_Location(double x, double y)
Definition kdtree.cpp:1281
bool Get_Point(sLong Index, double &x, double &y, double &z)
Definition kdtree.cpp:1305
bool Initialize(CSG_Shapes *pPoints, int zField=-1)
Definition kdtree.cpp:1246
bool Do_Use_All(void) const
Definition shapes.h:1372
virtual bool On_Parameter_Changed(class CSG_Parameters *pParameters, class CSG_Parameter *pParameter)
Definition kdtree.cpp:1156
virtual bool Create(class CSG_Parameters *pParameters, const CSG_String &Parent="", size_t minPoints=0)
Definition kdtree.cpp:1093
virtual bool On_Parameters_Enable(class CSG_Parameters *pParameters, class CSG_Parameter *pParameter)
Definition kdtree.cpp:1171
class CSG_Parameters * m_pParameters
Definition shapes.h:1389
const CSG_String & Get_Identifier(void) const
void Set_Enabled(bool bEnabled=true)
bool Set_Parameter(const CSG_String &ID, CSG_Parameter *pValue)
bool Clear(void)
Definition geo_tools.h:328
bool Add(double x, double y, double z)
double Get_Area(void) const
Definition geo_tools.h:513
double Get_XMin(void) const
Definition geo_tools.h:505
virtual TSG_Point Get_Point(int iPoint=0) const =0
virtual double Get_Z(int iPoint=0, int iPart=0, bool bAscending=true) const
Definition shapes.h:192
virtual const CSG_Rect & Get_Extent(void)
Definition shapes.h:810
virtual CSG_Shape * Get_Shape(const CSG_Point &Point, double Epsilon=0.)
Definition shapes.cpp:483
static CSG_String Format(const char *Format,...)
bool is_Empty(void) const
bool is_NoData(int Field) const
double asDouble(int Field) const
sLong Get_Count(void) const
Definition table.h:400
bool Create(const CSG_Vector &Vector)
bool Set_Rows(sLong nRows)
int Get_N(void) const
Definition mat_tools.h:384
double * Get_Data(void) const
Definition mat_tools.h:386
double SG_Get_Rounded_To_SignificantFigures(double Value, int Decimals)
double x
Definition geo_tools.h:129
double y
Definition geo_tools.h:129