Method classifyGround

filters/SMRFilter.cpp:296–398

}

void SMRFilter::classifyGround(PointViewPtr view, std::vector<double>& ZIpro)
{
    // "While many authors use a single value for the elevation threshold, we
    // suggest that a second parameter be used to increase the threshold on
    // steep slopes, transforming the threshold to a slope-dependent value. The
    // total permissible distance is then equal to a fixed elevation threshold
    // plus the scaling value multiplied by the slope of the DEM at each LIDAR
    // point. The rationale behind this approach is that small horizontal and
    // vertical displacements yield larger errors on steep slopes, and as a
    // result the BE/OBJ threshold distance should be more permissive at these
    // points."
    MatrixXd gsurfs(m_rows, m_cols);
    MatrixXd thresh(m_rows, m_cols);
    {
        MatrixXd ZIproM = Map<MatrixXd>(ZIpro.data(), m_rows, m_cols);
        MatrixXd scaled = ZIproM / m_args->m_cell;

        MatrixXd gx = math::gradX(scaled);
        MatrixXd gy = math::gradY(scaled);
        gsurfs = (gx.cwiseProduct(gx) + gy.cwiseProduct(gy)).cwiseSqrt();
        std::vector<double> gsurfsV(gsurfs.data(),
                                    gsurfs.data() + gsurfs.size());

        //ABELL - We can eliminate this copy if we're OK with not writing
        //  both the filled and non-filled array to output.
        std::vector<double> gsurfs_fillV = gsurfsV;
        knnfill(view, gsurfs_fillV);
        gsurfs = Map<MatrixXd>(gsurfs_fillV.data(), m_rows, m_cols);
        thresh =
            (m_args->m_threshold + m_args->m_scalar * gsurfs.array()).matrix();

        if (!m_args->m_dir.empty())
        {
            std::string fname =
                FileUtils::toAbsolutePath("gx.tif", m_args->m_dir);
            math::writeMatrix(gx, fname, "GTiff", m_args->m_cell, m_bounds, m_srs);

            fname = FileUtils::toAbsolutePath("gy.tif", m_args->m_dir);
            math::writeMatrix(gy, fname, "GTiff", m_args->m_cell, m_bounds, m_srs);

            fname = FileUtils::toAbsolutePath("gsurfs.tif", m_args->m_dir);
            math::writeMatrix(gsurfs, fname, "GTiff", m_args->m_cell, m_bounds, m_srs);

            fname = FileUtils::toAbsolutePath("gsurfs_fill.tif", m_args->m_dir);
            MatrixXd gsurfs_fill =
                Map<MatrixXd>(gsurfs_fillV.data(), m_rows, m_cols);
            math::writeMatrix(gsurfs_fill, fname, "GTiff", m_args->m_cell, m_bounds, m_srs);

            fname = FileUtils::toAbsolutePath("thresh.tif", m_args->m_dir);
            math::writeMatrix(thresh, fname, "GTiff", m_args->m_cell, m_bounds, m_srs);
        }
    }

    point_count_t ng(0);
    point_count_t g(0);
    for (PointRef p : *view)
    {
        double x = p.getFieldAs<double>(Id::X);