C++ PointViewPtr::calculateBounds方法代码示例

void NitfWriter::writeView(const PointViewPtr view)
{
    //ABELL - Think we can just get this from the LAS file header
    //  when we're done.
    view->calculateBounds(m_bounds);
    LasWriter::writeView(view);
}

void PointView::calculateBounds(const PointViewSet& set, BOX3D& output)
{
    for (auto iter = set.begin(); iter != set.end(); ++iter)
    {
        PointViewPtr buf = *iter;
        buf->calculateBounds(output);
    }
}

void P2gWriter::write(const PointViewPtr view)
{
    for (point_count_t idx = 0; idx < view->size(); idx++)
    {
        double x = view->getFieldAs<double>(Dimension::Id::X, idx);
        double y = view->getFieldAs<double>(Dimension::Id::Y, idx);
        double z = view->getFieldAs<double>(Dimension::Id::Z, idx);
        m_coordinates.push_back(Coordinate{x, y, z});
    }

    view->calculateBounds(m_bounds);
}

void PcdWriter::write(const PointViewPtr view)
{
    pcl::PointCloud<XYZIRGBA>::Ptr cloud(new pcl::PointCloud<XYZIRGBA>);
    BOX3D buffer_bounds;
    view->calculateBounds(buffer_bounds);
    pclsupport::PDALtoPCD(view, *cloud, buffer_bounds);

    pcl::PCDWriter w;

    if (m_compressed)
        w.writeBinaryCompressed<XYZIRGBA>(m_filename, *cloud);
    else
        w.writeASCII<XYZIRGBA>(m_filename, *cloud);
}

TEST(SplitterTest, test_buffer)
{
    Options readerOptions;
    readerOptions.add("filename", Support::datapath("las/1.2-with-color.las"));
    LasReader reader;
    reader.setOptions(readerOptions);

    Options splitterOptions;
    splitterOptions.add("length", 1000);
    splitterOptions.add("buffer", 20);

    SplitterFilter splitter;
    splitter.setOptions(splitterOptions);
    splitter.setInput(reader);

    PointTable table;
    splitter.prepare(table);
    PointViewSet viewSet = splitter.execute(table);

    std::vector<PointViewPtr> views;
    std::map<PointViewPtr, BOX2D> bounds;
    for (auto it = viewSet.begin(); it != viewSet.end(); ++it)
    {
        BOX2D b;
        PointViewPtr v = *it;
        v->calculateBounds(b);
        EXPECT_TRUE(b.maxx - b.minx <= 1040);    
        EXPECT_TRUE(b.maxy - b.miny <= 1040);
        bounds[v] = b;
        views.push_back(v);
    }

    auto sorter = [&bounds](PointViewPtr p1, PointViewPtr p2)
    {
        BOX2D b1 = bounds[p1];
        BOX2D b2 = bounds[p2];

        return b1.minx < b2.minx ?  true :
            b1.minx > b2.minx ? false :
            b1.miny < b2.miny;
    };
    std::sort(views.begin(), views.end(), sorter);

    EXPECT_EQ(views.size(), 24u);
    size_t counts[] = {26, 26, 3, 28, 27, 13, 14, 65, 80, 47, 80, 89, 94,
        77, 5, 79, 65, 34, 63, 67, 74, 69, 36, 5};
    size_t i = 0;
    for (PointViewPtr view : views)
        EXPECT_EQ(view->size(), counts[i++]);
}

void P2gWriter::write(const PointViewPtr view)
{
    std::string z_name = getOptions().getValueOrDefault<std::string>("Z", "Z");


    for (point_count_t idx = 0; idx < view->size(); idx++)
    {
        double x = view->getFieldAs<double>(Dimension::Id::X, idx);
        double y = view->getFieldAs<double>(Dimension::Id::Y, idx);
        double z = view->getFieldAs<double>(Dimension::Id::Z, idx);
        m_coordinates.push_back(boost::tuple<double, double, double>(x, y, z));
    }

    view->calculateBounds(m_bounds);
}

PointViewSet RadiusOutlierFilter::run(PointViewPtr input)
{
    bool logOutput = log()->getLevel() > LogLevel::Debug1;
    if (logOutput)
        log()->floatPrecision(8);
    log()->get(LogLevel::Debug2) << "Process RadiusOutlierFilter...\n";

    // convert PointView to PointXYZ
    typedef pcl::PointCloud<pcl::PointXYZ> Cloud;
    Cloud::Ptr cloud(new Cloud);
    BOX3D bounds;
    input->calculateBounds(bounds);
    pclsupport::PDALtoPCD(input, *cloud, bounds);

    pclsupport::setLogLevel(log()->getLevel());

    // setup the outlier filter
    pcl::RadiusOutlierRemoval<pcl::PointXYZ> ror(true);
    ror.setInputCloud(cloud);
    ror.setMinNeighborsInRadius(m_min_neighbors);
    ror.setRadiusSearch(m_radius);

    pcl::PointCloud<pcl::PointXYZ> output;
    ror.setNegative(true);
    ror.filter(output);

    // filtered to return inliers
    pcl::PointIndicesPtr inliers(new pcl::PointIndices);
    ror.getRemovedIndices(*inliers);

    PointViewSet viewSet;
    if (inliers->indices.empty())
    {
        log()->get(LogLevel::Warning) << "Requested filter would remove all points. Try a larger radius/smaller minimum neighbors.\n";
        viewSet.insert(input);
        return viewSet;
    }

    // inverse are the outliers
    std::vector<int> outliers(input->size()-inliers->indices.size());
    for (PointId i = 0, j = 0, k = 0; i < input->size(); ++i)
    {
        if (i == (PointId)inliers->indices[j])
        {
            j++;
            continue;
        }
        outliers[k++] = i;
    }

    if (!outliers.empty() && (m_classify || m_extract))
    {

        if (m_classify)
        {
            log()->get(LogLevel::Debug2) << "Labeled " << outliers.size() << " outliers as noise!\n";

            // set the classification label of outlier returns as 18
            // (corresponding to ASPRS LAS specification for high noise)
            for (const auto& i : outliers)
            {
                input->setField(Dimension::Id::Classification, i, 18);
            }

            viewSet.insert(input);
        }

        if (m_extract)
        {
            log()->get(LogLevel::Debug2) << "Extracted " << inliers->indices.size() << " inliers!\n";

            // create new PointView containing only outliers
            PointViewPtr output = input->makeNew();
            for (const auto& i : inliers->indices)
            {
                output->appendPoint(*input, i);
            }

            viewSet.erase(input);
            viewSet.insert(output);
        }
    }
    else
    {
        if (outliers.empty())
            log()->get(LogLevel::Warning) << "Filtered cloud has no outliers!\n";

        if (!(m_classify || m_extract))
            log()->get(LogLevel::Warning) << "Must choose --classify or --extract\n";

        // return the input buffer unchanged
        viewSet.insert(input);
    }

    return viewSet;
}

std::vector<PointId> SMRFilter::processGround(PointViewPtr view)
{
    log()->get(LogLevel::Info) << "processGround: Running SMRF...\n";

    // The algorithm consists of four conceptually distinct stages. The first is
    // the creation of the minimum surface (ZImin). The second is the processing
    // of the minimum surface, in which grid cells from the raster are
    // identified as either containing bare earth (BE) or objects (OBJ). This
    // second stage represents the heart of the algorithm. The third step is the
    // creation of a DEM from these gridded points. The fourth step is the
    // identification of the original LIDAR points as either BE or OBJ based on
    // their relationship to the interpolated

    std::vector<PointId> groundIdx;

    BOX2D bounds;
    view->calculateBounds(bounds);
    SpatialReference srs(view->spatialReference());

    // Determine the number of rows and columns at the given cell size.
    m_numCols = ((bounds.maxx - bounds.minx) / m_cellSize) + 1;
    m_numRows = ((bounds.maxy - bounds.miny) / m_cellSize) + 1;

    MatrixXd cx(m_numRows, m_numCols);
    MatrixXd cy(m_numRows, m_numCols);
    for (auto c = 0; c < m_numCols; ++c)
    {
        for (auto r = 0; r < m_numRows; ++r)
        {
            cx(r, c) = bounds.minx + (c + 0.5) * m_cellSize;
            cy(r, c) = bounds.miny + (r + 0.5) * m_cellSize;
        }
    }

    // STEP 1:

    // As with many other ground filtering algorithms, the first step is
    // generation of ZImin from the cell size parameter and the extent of the
    // data. The two vectors corresponding to [min:cellSize:max] for each
    // coordinate – xi and yi – may be supplied by the user or may be easily and
    // automatically calculated from the data. Without supplied ranges, the SMRF
    // algorithm creates a raster from the ceiling of the minimum to the floor
    // of the maximum values for each of the (x,y) dimensions. If the supplied
    // cell size parameter is not an integer, the same general rule applies to
    // values evenly divisible by the cell size. For example, if cell size is
    // equal to 0.5 m, and the x values range from 52345.6 to 52545.4, the range
    // would be [52346 52545].

    // The minimum surface grid ZImin defined by vectors (xi,yi) is filled with
    // the nearest, lowest elevation from the original point cloud (x,y,z)
    // values, provided that the distance to the nearest point does not exceed
    // the supplied cell size parameter. This provision means that some grid
    // points of ZImin will go unfilled. To fill these values, we rely on
    // computationally inexpensive image inpainting techniques. Image inpainting
    // involves the replacement of the empty cells in an image (or matrix) with
    // values calculated from other nearby values. It is a type of interpolation
    // technique derived from artistic replacement of damaged portions of
    // photographs and paintings, where preservation of texture is an important
    // concern (Bertalmio et al., 2000). When empty values are spread through
    // the image, and the ratio of filled to empty pixels is quite high, most
    // methods of inpainting will produce satisfactory results. In an evaluation
    // of inpainting methods on ground identification from the final terrain
    // model, we found that Laplacian techniques produced error rates nearly
    // three times higher than either an average of the eight nearest neighbors
    // or D’Errico’s spring-metaphor inpainting technique (D’Errico, 2004). The
    // spring-metaphor technique imagines springs connecting each cell with its
    // eight adjacent neighbors, where the inpainted value corresponds to the
    // lowest energy state of the set, and where the entire (sparse) set of
    // linear equations is solved using partial differential equations. Both of
    // these latter techniques were nearly the same with regards to total error,
    // with the spring technique performing slightly better than the k-nearest
    // neighbor (KNN) approach.

    MatrixXd ZImin = eigen::createMinMatrix(*view.get(), m_numRows, m_numCols,
                                            m_cellSize, bounds);

    // MatrixXd ZImin_painted = inpaintKnn(cx, cy, ZImin);
    // MatrixXd ZImin_painted = TPS(cx, cy, ZImin);
    MatrixXd ZImin_painted = expandingTPS(cx, cy, ZImin);

    if (!m_outDir.empty())
    {
        std::string filename = FileUtils::toAbsolutePath("zimin.tif", m_outDir);
        eigen::writeMatrix(ZImin, filename, "GTiff", m_cellSize, bounds, srs);

        filename = FileUtils::toAbsolutePath("zimin_painted.tif", m_outDir);
        eigen::writeMatrix(ZImin_painted, filename, "GTiff", m_cellSize, bounds, srs);
    }

    ZImin = ZImin_painted;

    // STEP 2:

    // The second stage of the ground identification algorithm involves the
    // application of a progressive morphological filter to the minimum surface
    // grid (ZImin). At the first iteration, the filter applies an image opening
    // operation to the minimum surface. An opening operation consists of an
    // application of an erosion filter followed by a dilation filter. The
    // erosion acts to snap relative high values to relative lows, where a
    // supplied window radius and shape (or structuring element) defines the
//.........这里部分代码省略.........

TEST(SplitterTest, test_tile_filter)
{
    StageFactory f;

    // create the reader
    Options ops1;
    ops1.add("filename", Support::datapath("las/1.2-with-color.las"));
    LasReader r;
    r.setOptions(ops1);

    Options o;
    Option length("length", 1000);
    o.add(length);

    // create the tile filter and prepare
    SplitterFilter s;
    s.setOptions(o);
    s.setInput(r);

    PointTable table;
    s.prepare(table);
    PointViewSet viewSet = s.execute(table);

    std::vector<PointViewPtr> views;
    std::map<PointViewPtr, BOX2D> bounds;

    for (auto it = viewSet.begin(); it != viewSet.end(); ++it)
    {
        BOX2D b;
        PointViewPtr v = *it;
        v->calculateBounds(b);
        EXPECT_TRUE(b.maxx - b.minx <= 1000);    
        EXPECT_TRUE(b.maxy - b.miny <= 1000);    

        for (auto& p : bounds)
            EXPECT_FALSE(p.second.overlaps(b));

        bounds[v] = b;
        views.push_back(v);
    }

    auto sorter = [&bounds](PointViewPtr p1, PointViewPtr p2)
    {
        BOX2D b1 = bounds[p1];
        BOX2D b2 = bounds[p2];

        return b1.minx < b2.minx ?  true :
            b1.minx > b2.minx ? false :
            b1.miny < b2.miny;
    };
    std::sort(views.begin(), views.end(), sorter);

    EXPECT_EQ(views.size(), 24u);
    size_t counts[] = {24, 25, 2, 26, 27, 10, 82, 68, 43, 57, 7, 71, 73,
        61, 33, 84, 74, 4, 59, 70, 67, 34, 60, 4 };
    for (size_t i = 0; i < views.size(); ++i)
    {
        PointViewPtr view = views[i];
        EXPECT_EQ(view->size(), counts[i]);
    }
}

PointViewSet StatisticalOutlierFilter::run(PointViewPtr input)
{
    bool logOutput = log()->getLevel() > LogLevel::Debug1;
    if (logOutput)
        log()->floatPrecision(8);
    log()->get(LogLevel::Debug2) << "Process StatisticalOutlierFilter...\n";

    // convert PointView to PointXYZ
    typedef pcl::PointCloud<pcl::PointXYZ> Cloud;
    Cloud::Ptr cloud(new Cloud);
    BOX3D bounds;
    input->calculateBounds(bounds);
    pclsupport::PDALtoPCD(input, *cloud, bounds);

    // PCL should provide console output at similar verbosity level as PDAL
    int level = log()->getLevel();
    switch (level)
    {
        case 0:
            pcl::console::setVerbosityLevel(pcl::console::L_ALWAYS);
            break;
        case 1:
            pcl::console::setVerbosityLevel(pcl::console::L_ERROR);
            break;
        case 2:
            pcl::console::setVerbosityLevel(pcl::console::L_WARN);
            break;
        case 3:
            pcl::console::setVerbosityLevel(pcl::console::L_INFO);
            break;
        case 4:
            pcl::console::setVerbosityLevel(pcl::console::L_DEBUG);
            break;
        default:
            pcl::console::setVerbosityLevel(pcl::console::L_VERBOSE);
            break;
    }

    // setup the outlier filter
    pcl::StatisticalOutlierRemoval<pcl::PointXYZ> sor(true);
    sor.setInputCloud(cloud);
    sor.setMeanK(m_meanK);
    sor.setStddevMulThresh(m_multiplier);

    pcl::PointCloud<pcl::PointXYZ> output;
    sor.setNegative(true);
    sor.filter(output);

    // filtered to return inliers
    pcl::PointIndicesPtr inliers(new pcl::PointIndices);
    sor.getRemovedIndices(*inliers);

    log()->get(LogLevel::Debug2) << inliers->indices.size() << std::endl;

    PointViewSet viewSet;
    if (inliers->indices.empty())
    {
        log()->get(LogLevel::Warning) << "Requested filter would remove all points. Try increasing the multiplier.\n";
        viewSet.insert(input);
        return viewSet;
    }

    // inverse are the outliers
    std::vector<int> outliers(input->size()-inliers->indices.size());
    for (PointId i = 0, j = 0, k = 0; i < input->size(); ++i)
    {
        if (i == (PointId)inliers->indices[j])
        {
            j++;
            continue;
        }
        outliers[k++] = i;
    }

    if (!outliers.empty() && (m_classify || m_extract))
    {

        if (m_classify)
        {
            log()->get(LogLevel::Debug2) << "Labeled " << outliers.size() << " outliers as noise!\n";

            // set the classification label of outlier returns as 18
            // (corresponding to ASPRS LAS specification for high noise)
            for (const auto& i : outliers)
            {
                input->setField(Dimension::Id::Classification, i, 18);
            }

            viewSet.insert(input);
        }

        if (m_extract)
        {
            log()->get(LogLevel::Debug2) << "Extracted " << inliers->indices.size() << " inliers!\n";

            // create new PointView containing only outliers
            PointViewPtr output = input->makeNew();
            for (const auto& i : inliers->indices)
            {
                output->appendPoint(*input, i);
//.........这里部分代码省略.........

void NitfWriter::write(const PointViewPtr view)
{
    m_bounds.grow(view->calculateBounds(true));

    LasWriter::write(view);
}

PointViewSet DartSampleFilter::run(PointViewPtr input)
{
    PointViewSet viewSet;
    PointViewPtr output = input->makeNew();

    log()->floatPrecision(2);
    log()->get(LogLevel::Info) << "DartSampleFilter (radius="
                               << m_radius << ")\n";

    BOX3D buffer_bounds;
    input->calculateBounds(buffer_bounds);

    typedef pcl::PointCloud<pcl::PointXYZ> Cloud;
    Cloud::Ptr cloud(new Cloud);
    pclsupport::PDALtoPCD(input, *cloud, buffer_bounds);

    int level = log()->getLevel();
    switch (level)
    {
        case 0:
            pcl::console::setVerbosityLevel(pcl::console::L_ALWAYS);
            break;
        case 1:
            pcl::console::setVerbosityLevel(pcl::console::L_ERROR);
            break;
        case 2:
            pcl::console::setVerbosityLevel(pcl::console::L_WARN);
            break;
        case 3:
            pcl::console::setVerbosityLevel(pcl::console::L_INFO);
            break;
        case 4:
            pcl::console::setVerbosityLevel(pcl::console::L_DEBUG);
            break;
        default:
            pcl::console::setVerbosityLevel(pcl::console::L_VERBOSE);
            break;
    }

    pcl::DartSample<pcl::PointXYZ> ds;
    ds.setInputCloud(cloud);
    ds.setRadius(m_radius);

    std::vector<int> samples;
    ds.filter(samples);

    if (samples.empty())
    {
        log()->get(LogLevel::Warning) << "Filtered cloud has no points!\n";
        return viewSet;
    }

    for (const auto& i : samples)
        output->appendPoint(*input, i);

    double frac = (double)samples.size() / (double)cloud->size();
    log()->get(LogLevel::Info) << "Retaining " << samples.size() << " of "
                               << cloud->size() << " points ("
                               <<  100*frac
                               << "%)\n";

    viewSet.insert(output);
    return viewSet;
}