| 102 | } |
| 103 | |
| 104 | void ReciprocityFilter::setReciprocity(PointView& view, const PointId& i) |
| 105 | { |
| 106 | // Find k-nearest neighbors of i. |
| 107 | const KD3Index& kdi = view.build3dIndex(); |
| 108 | PointIdList ni = kdi.neighbors(i, m_knn + 1); |
| 109 | |
| 110 | // Initialize number of unidirectional neighbors to 0. |
| 111 | point_count_t uni(0); |
| 112 | |
| 113 | // Visit each neighbor of i, finding its k-nearest neighbors. If i is |
| 114 | // not a nearest neighbor of one of its neighbors, increment uni. |
| 115 | for (PointId const& j : ni) |
| 116 | { |
| 117 | // The query point itself will always show up as a neighbor and can |
| 118 | // be skipped. |
| 119 | if (j == i) |
| 120 | continue; |
| 121 | |
| 122 | // Find k-nearest neighbors of j. |
| 123 | PointIdList nj = kdi.neighbors(j, m_knn + 1); |
| 124 | |
| 125 | // If i is not a neighbor of j, increment uni. |
| 126 | if (std::find(nj.begin(), nj.end(), i) == nj.end()) |
| 127 | ++uni; |
| 128 | } |
| 129 | |
| 130 | // Compute reciprocity as percentage of neighbors that do NOT contain |
| 131 | // id as a neighbor. |
| 132 | double reciprocity = 100.0 * uni / m_knn; |
| 133 | view.setField(Id::Reciprocity, i, reciprocity); |
| 134 | } |
| 135 | |
| 136 | } // namespace pdal |