Method progressiveFilter

filters/SMRFilter.cpp:669–733

Iteratively open the estimated surface. progressiveFilter can be used to identify both low points and object (i.e., non-ground) points, depending on the inputs.

// identify both low points and object (i.e., non-ground) points, depending on
// the inputs.
std::vector<int> SMRFilter::progressiveFilter(std::vector<double> const& ZImin,
                                              double slope, double max_window)
{
    // "The maximum window radius is supplied as a distance metric (e.g., 21 m),
    // but is internally converted to a pixel equivalent by dividing it by the
    // cell size and rounding the result toward positive infinity (i.e., taking
    // the ceiling value)."
    int max_radius = static_cast<int>(std::ceil(max_window / m_args->m_cell));
    std::vector<double> prevSurface = ZImin;
    std::vector<double> erosion = ZImin;

    // "...the radius of the element at each step [is] increased by one pixel
    // from a starting value of one pixel to the pixel equivalent of the maximum
    // value."
    std::vector<int> Obj(m_rows * m_cols, 0);
    for (int radius = 1; radius <= max_radius; ++radius)
    {
        // "On the first iteration, the minimum surface (ZImin) is opened using
        // a disk-shaped structuring element with a radius of one pixel."
        math::erodeDiamond(erosion, m_rows, m_cols, 1);
        std::vector<double> curOpening = erosion;
        math::dilateDiamond(curOpening, m_rows, m_cols, radius);

        // "An elevation threshold is then calculated, where the value is equal
        // to the supplied slope tolerance parameter multiplied by the product
        // of the window radius and the cell size."
        double threshold = slope * m_args->m_cell * radius;

        // "This elevation threshold is applied to the difference of the minimum
        // and the opened surfaces."

        // Need to provide means of diffing two vectors.
        std::vector<double> diff;
        std::transform(prevSurface.begin(), prevSurface.end(),
                       curOpening.begin(), std::back_inserter(diff),
                       [&](double l, double r) { return std::fabs(l - r); });

        // "Any grid cell with a difference value exceeding the calculated
        // elevation threshold for the iteration is then flagged as an OBJ
        // cell."
        std::vector<int> foo;
        std::transform(diff.begin(), diff.end(), std::back_inserter(foo),
                       [threshold](double x) {
                           return (x > threshold) ? int(1) : int(0);
                       });
        std::transform(Obj.begin(), Obj.end(), foo.begin(), Obj.begin(),
                       [](int a, int b) { return (std::max)(a, b); });

        // "The algorithm then proceeds to the next window radius (up to the
        // maximum), and proceeds as above with the last opened surface acting
        // as the minimum surface for the next difference calculation."
        prevSurface = curOpening;

        size_t ng = std::count(Obj.begin(), Obj.end(), 1);
        size_t g(Obj.size() - ng);
        double p(100.0 * double(ng) / double(Obj.size()));
        log()->floatPrecision(2);
        log()->get(LogLevel::Debug) << "progressiveFilter: radius = " << radius