| 71 | } |
| 72 | |
| 73 | PointViewSet MADFilter::run(PointViewPtr view) |
| 74 | { |
| 75 | using namespace Dimension; |
| 76 | |
| 77 | PointViewPtr output = view->makeNew(); |
| 78 | |
| 79 | auto estimate_median = [](std::vector<double> vals) |
| 80 | { |
| 81 | std::nth_element(vals.begin(), vals.begin()+vals.size()/2, vals.end()); |
| 82 | return *(vals.begin()+vals.size()/2); |
| 83 | }; |
| 84 | |
| 85 | std::vector<double> z(view->size()); |
| 86 | for (PointId j = 0; j < view->size(); ++j) |
| 87 | z[j] = view->getFieldAs<double>(m_dimId, j); |
| 88 | |
| 89 | double median = estimate_median(z); |
| 90 | log()->get(LogLevel::Debug) << getName() << |
| 91 | " estimated median value: " << median << std::endl; |
| 92 | |
| 93 | std::transform(z.begin(), z.end(), z.begin(), |
| 94 | [median](double v) { return std::fabs(v - median); }); |
| 95 | double mad = estimate_median(z)*m_madMultiplier; |
| 96 | log()->get(LogLevel::Debug) << getName() << " mad " << mad << std::endl; |
| 97 | |
| 98 | for (PointId j = 0; j < view->size(); ++j) |
| 99 | { |
| 100 | if (z[j]/mad < m_multiplier) |
| 101 | output->appendPoint(*view, j); |
| 102 | } |
| 103 | |
| 104 | double low_fence = median - m_multiplier * mad; |
| 105 | double hi_fence = median + m_multiplier * mad; |
| 106 | |
| 107 | log()->get(LogLevel::Debug) << getName() << " cropping " << m_dimName |
| 108 | << " in the range (" << low_fence |
| 109 | << "," << hi_fence << ")" << std::endl; |
| 110 | |
| 111 | PointViewSet viewSet; |
| 112 | viewSet.insert(output); |
| 113 | return viewSet; |
| 114 | } |
| 115 | |
| 116 | } // namespace pdal |