| 71 | } |
| 72 | |
| 73 | void LOFFilter::filter(PointView& view) |
| 74 | { |
| 75 | const KD3Index& index = view.build3dIndex(); |
| 76 | |
| 77 | // Increment the minimum number of points, as knnSearch will be returning |
| 78 | // the neighbors along with the query point. |
| 79 | m_minpts++; |
| 80 | |
| 81 | // First pass: Compute the k-distance for each point. |
| 82 | // The k-distance is the Euclidean distance to k-th nearest neighbor. |
| 83 | log()->get(LogLevel::Debug) << "Computing k-distances...\n"; |
| 84 | |
| 85 | typedef std::pair<PointId, size_t> PointIdKPair; |
| 86 | typedef std::pair<PointId, double> PointIdDistPair; |
| 87 | typedef std::map<PointIdKPair, PointIdDistPair> AdjacencyMap; |
| 88 | AdjacencyMap mat; |
| 89 | |
| 90 | for (PointId i = 0; i < view.size(); ++i) |
| 91 | { |
| 92 | PointIdList indices(m_minpts); |
| 93 | std::vector<double> sqr_dists(m_minpts); |
| 94 | index.knnSearch(i, m_minpts, &indices, &sqr_dists); |
| 95 | |
| 96 | for (size_t j = 0; j < m_minpts; ++j) |
| 97 | mat[std::make_pair(i, j)] = std::make_pair(indices[j], std::sqrt(sqr_dists[j])); |
| 98 | |
| 99 | view.setField(Id::NNDistance, i, std::sqrt(sqr_dists[m_minpts - 1])); |
| 100 | } |
| 101 | |
| 102 | // Second pass: Compute the local reachability distance for each point. |
| 103 | // For each neighbor point, the reachability distance is the maximum value |
| 104 | // of that neighbor's k-distance and the distance between the neighbor and |
| 105 | // the current point. The lrd is the inverse of the mean of the reachability |
| 106 | // distances. |
| 107 | log()->get(LogLevel::Debug) << "Computing lrd...\n"; |
| 108 | for (PointId i = 0; i < view.size(); ++i) |
| 109 | { |
| 110 | double M1 = 0.0; |
| 111 | point_count_t n = 0; |
| 112 | for (size_t j = 0; j < m_minpts; ++j) |
| 113 | { |
| 114 | auto val = mat[std::make_pair(i, j)]; |
| 115 | double k = view.getFieldAs<double>(Id::NNDistance, val.first); |
| 116 | double reachdist = (std::max)(k, val.second); |
| 117 | M1 += (reachdist - M1) / ++n; |
| 118 | } |
| 119 | view.setField(Id::LocalReachabilityDistance, i, 1.0 / M1); |
| 120 | } |
| 121 | |
| 122 | // Third pass: Compute the local outlier factor for each point. |
| 123 | // The LOF is the average of the lrd's for a neighborhood of points. |
| 124 | log()->get(LogLevel::Debug) << "Computing LOF...\n"; |
| 125 | for (PointId i = 0; i < view.size(); ++i) |
| 126 | { |
| 127 | double lrdp = view.getFieldAs<double>(Id::LocalReachabilityDistance, i); |
| 128 | double M1 = 0.0; |
| 129 | point_count_t n = 0; |
| 130 | for (size_t j = 0; j < m_minpts; ++j) |