1 /******************************************************************************
2 * Copyright (c) 2019, Bradley J Chambers (brad.chambers@gmail.com)
3 *
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following
8 * conditions are met:
9 *
10 * * Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * * Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in
14 * the documentation and/or other materials provided
15 * with the distribution.
16 * * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the
17 * names of its contributors may be used to endorse or promote
18 * products derived from this software without specific prior
19 * written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
24 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
25 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
26 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
27 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
28 * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
29 * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
30 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
31 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
32 * OF SUCH DAMAGE.
33 ****************************************************************************/
34
35 // Adapted from MIT-licensed implemenation provided by
36 // https://github.com/intel-isl/Open3D/pull/1038.
37
38 #include "DBSCANFilter.hpp"
39
40 #include <pdal/KDIndex.hpp>
41
42 #include <string>
43 #include <unordered_set>
44
45 namespace pdal
46 {
47
48 using namespace Dimension;
49
50 static StaticPluginInfo const s_info{
51 "filters.dbscan", "DBSCAN Clustering.",
52 "http://pdal.io/stages/filters.dbscan.html"};
53
CREATE_STATIC_STAGE(DBSCANFilter,s_info)54 CREATE_STATIC_STAGE(DBSCANFilter, s_info)
55
56 std::string DBSCANFilter::getName() const
57 {
58 return s_info.name;
59 }
60
DBSCANFilter()61 DBSCANFilter::DBSCANFilter() : Filter()
62 {
63 }
64
addArgs(ProgramArgs & args)65 void DBSCANFilter::addArgs(ProgramArgs& args)
66 {
67 args.add("min_points", "Min points per cluster", m_minPoints,
68 static_cast<uint64_t>(6));
69 args.add("eps", "Epsilon", m_eps, 1.0);
70 args.add("dimensions", "Dimensions to cluster", m_dimStringList,
71 {"X", "Y", "Z"});
72 }
73
addDimensions(PointLayoutPtr layout)74 void DBSCANFilter::addDimensions(PointLayoutPtr layout)
75 {
76 layout->registerDim(Id::ClusterID);
77 }
78
prepared(PointTableRef table)79 void DBSCANFilter::prepared(PointTableRef table)
80 {
81 const PointLayoutPtr layout(table.layout());
82
83 if (m_dimStringList.size())
84 {
85 for (std::string& s : m_dimStringList)
86 {
87 Id id = layout->findDim(s);
88 if (id == Id::Unknown)
89 throwError("Invalid dimension '" + s +
90 "' specified for "
91 "'dimensions' option.");
92 m_dimIdList.push_back(id);
93 }
94 }
95 }
96
filter(PointView & view)97 void DBSCANFilter::filter(PointView& view)
98 {
99 // Construct KDFlexIndex for radius search.
100 KDFlexIndex kdfi(view, m_dimIdList);
101 kdfi.build();
102
103 // First pass through point cloud precomputes neighbor indices and
104 // initializes ClusterID to -2.
105 std::vector<PointIdList> neighbors(view.size());
106 for (PointId idx = 0; idx < view.size(); ++idx)
107 {
108 neighbors[idx] = kdfi.radius(idx, m_eps);
109 view.setField(Id::ClusterID, idx, -2);
110 }
111
112 // Second pass through point cloud performs DBSCAN clustering.
113 int64_t cluster_label = 0;
114 for (PointId idx = 0; idx < view.size(); ++idx)
115 {
116 // Point has already been labeled, so move on to next.
117 if (view.getFieldAs<int64_t>(Id::ClusterID, idx) != -2)
118 continue;
119
120 // Density of the neighborhood does not meet minimum number of points
121 // constraint, label as noise.
122 if (neighbors[idx].size() < m_minPoints)
123 {
124 view.setField(Id::ClusterID, idx, -1);
125 continue;
126 }
127
128 // Initialize some bookkeeping to track which neighbors have been
129 // considered for addition to the current cluster.
130 std::unordered_set<PointId> neighbors_next(neighbors[idx].begin(),
131 neighbors[idx].end());
132 std::unordered_set<PointId> neighbors_visited;
133 neighbors_visited.insert(idx);
134
135 // Unlabeled point encountered; assign cluster label.
136 view.setField(Id::ClusterID, idx, cluster_label);
137
138 // Consider all neighbors.
139 while (!neighbors_next.empty())
140 {
141 // Select first neighbor and move it to the visited set.
142 PointId p = *neighbors_next.begin();
143 neighbors_next.erase(neighbors_next.begin());
144 neighbors_visited.insert(p);
145
146 // Reassign cluster label to neighbor previously marked as noise.
147 if (view.getFieldAs<int64_t>(Id::ClusterID, p) == -1)
148 view.setField(Id::ClusterID, p, cluster_label);
149
150 // Neighbor has already been labeled, so move on to next.
151 if (view.getFieldAs<int64_t>(Id::ClusterID, p) != -2)
152 continue;
153
154 // Assign cluster label to neighbor.
155 view.setField(Id::ClusterID, p, cluster_label);
156
157 // If density of neighbor's neighborhood is sufficient, add it's
158 // neighbors to the set of neighbors to consider if they are not
159 // already.
160 if (neighbors[p].size() >= m_minPoints)
161 {
162 for (PointId q : neighbors[p])
163 {
164 if (neighbors_visited.count(q) == 0)
165 neighbors_next.insert(q);
166 }
167 }
168 }
169
170 cluster_label++;
171 }
172 }
173
174 } // namespace pdal
175