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

PointViewSet LocateFilter::run(PointViewPtr inView)
{
    PointViewSet viewSet;
    if (!inView->size())
        return viewSet;

    PointId minidx, maxidx;
    double minval = (std::numeric_limits<double>::max)();
    double maxval = std::numeric_limits<double>::lowest();

    for (PointId idx = 0; idx < inView->size(); idx++)
    {
        double val = inView->getFieldAs<double>(m_dimId, idx);
        if (val > maxval)
        {
            maxval = val;
            maxidx = idx;
        }
        if (val < minval)
        {
            minval = val;
            minidx = idx;
        }
    }

    PointViewPtr outView = inView->makeNew();

    if (Utils::iequals("min", m_minmax))
        outView->appendPoint(*inView.get(), minidx);
    if (Utils::iequals("max", m_minmax))
        outView->appendPoint(*inView.get(), maxidx);

    viewSet.insert(outView);
    return viewSet;
}

TEST(Ilvis2ReaderTest, testReadHigh)
{
    Option filename("filename",
        Support::datapath("ilvis2/ILVIS2_TEST_FILE.TXT"));
    Options options(filename);
    options.add("mapping","high");
    std::shared_ptr<Ilvis2Reader> reader(new Ilvis2Reader);
    reader->setOptions(options);

    PointTable table;

    reader->prepare(table);
    PointViewSet viewSet = reader->execute(table);
    EXPECT_EQ(viewSet.size(), 1u);
    PointViewPtr view = *viewSet.begin();

    EXPECT_EQ(view->size(), 3u);

    checkPoint(*view.get(), 0, 42504.48313,
             78.307672,-58.785213,1956.777
            );

    checkPoint(*view.get(), 1, 42504.48512,
             78.307592, 101.215097, 1956.588
            );

    checkPoint(*view.get(), 2, 42504.48712,
             78.307512, -58.78459, 2956.667
            );
}

PointViewSet DecimationFilter::run(PointViewPtr inView)
{
    PointViewSet viewSet;
    PointViewPtr outView = inView->makeNew();
    decimate(*inView.get(), *outView.get());
    viewSet.insert(outView);
    return viewSet;
}

void BpfWriter::writeView(const PointViewPtr dataShared)
{
    setAutoXForm(dataShared);

    // Avoid reference count overhead internally.
    const PointView* data(dataShared.get());

    // We know that X, Y and Z are dimensions 0, 1 and 2.
    m_dims[0].m_offset = m_xXform.m_offset;
    m_dims[1].m_offset = m_yXform.m_offset;
    m_dims[2].m_offset = m_zXform.m_offset;

    switch (m_header.m_pointFormat)
    {
    case BpfFormat::PointMajor:
        writePointMajor(data);
        break;
    case BpfFormat::DimMajor:
        writeDimMajor(data);
        break;
    case BpfFormat::ByteMajor:
        writeByteMajor(data);
        break;
    }
    m_header.m_numPts += data->size();
}

MetadataNode InfoKernel::dumpQuery(PointViewPtr inView) const
{
    int count;
    std::string location;

    // See if there's a provided point count.
    StringList parts = Utils::split2(m_queryPoint, '/');
    if (parts.size() == 2)
    {
        location = parts[0];
        count = atoi(parts[1].c_str());
    }
    else if (parts.size() == 1)
    {
        location = parts[0];
        count = inView->size();
    }
    else
        count = 0;
    if (count == 0)
        throw pdal_error("Invalid location specificiation. "
            "--query=\"X,Y[/count]\"");

    auto seps = [](char c){ return (c == ',' || c == '|' || c == ' '); };

    std::vector<std::string> tokens = Utils::split2(location, seps);
    std::vector<double> values;
    for (auto ti = tokens.begin(); ti != tokens.end(); ++ti)
    {
        double d;
        if (Utils::fromString(*ti, d))
            values.push_back(d);
    }

    if (values.size() != 2 && values.size() != 3)
        throw pdal_error("--points must be two or three values");

    PointViewPtr outView = inView->makeNew();

    std::vector<PointId> ids;
    if (values.size() >= 3)
    {
        KD3Index kdi(*inView);
        kdi.build();
        ids = kdi.neighbors(values[0], values[1], values[2], count);
    }
    else
    {
        KD2Index kdi(*inView);
        kdi.build();
        ids = kdi.neighbors(values[0], values[1], count);
    }

    for (auto i = ids.begin(); i != ids.end(); ++i)
        outView->appendPoint(*inView.get(), *i);

    return outView->toMetadata();
}

TEST(PLangTest, PLangTest_array)
{
    PointTable table;
    PointViewPtr view = makeTestView(table, 40);

    plang::Array array;
    array.update(view);
    verifyTestView(*view.get(), 4);

}

PointViewSet MergeFilter::run(PointViewPtr in)
{
    PointViewSet viewSet;

    // If the SRS of all the point views aren't the same, print a warning
    // unless we're explicitly overriding the SRS.
    if (getSpatialReference().empty() &&
      (in->spatialReference() != m_view->spatialReference()))
        log()->get(LogLevel::Warning) << getName() << ": merging points "
            "with inconsistent spatial references." << std::endl;
    m_view->append(*in.get());
    viewSet.insert(m_view);
    return viewSet;
}

TEST(SbetReaderTest, testRead)
{
    Option filename("filename", Support::datapath("sbet/2-points.sbet"), "");
    Options options(filename);
    std::shared_ptr<SbetReader> reader(new SbetReader);
    reader->setOptions(options);

    PointTable table;

    reader->prepare(table);
    PointViewSet viewSet = reader->execute(table);
    EXPECT_EQ(viewSet.size(), 1u);
    PointViewPtr view = *viewSet.begin();

    EXPECT_EQ(view->size(), 2u);

    checkPoint(*view.get(), 0,
               1.516310028360710e+05, 5.680211852972264e-01,
               -2.041654392303940e+00, 1.077152953296560e+02,
               -2.332420866600025e+00, -3.335067504871401e-01,
               -3.093961631767838e-02, -2.813407149321339e-02,
               -2.429905393889139e-02, 3.046773230278662e+00,
               -2.198414736922658e-02, 7.859639737752390e-01,
               7.849084719295495e-01, -2.978807916450262e-01,
               6.226807982589819e-05, 9.312162756440178e-03,
               7.217812320996525e-02);
    checkPoint(*view.get(), 1,
               1.516310078318641e+05, 5.680211834722869e-01,
               -2.041654392034053e+00, 1.077151424357507e+02,
               -2.336228229691271e+00, -3.324663118952635e-01,
               -3.022948961008987e-02, -2.813856631423094e-02,
               -2.425215669392169e-02, 3.047131105236811e+00,
               -2.198416007932108e-02, 8.397590491636475e-01,
               3.252165276637165e-01, -1.558883225990844e-01,
               8.379685112283802e-04, 7.372886784718076e-03,
               7.179027672314571e-02);
}

void LasWriter::write(const PointViewPtr view)
{
    setAutoOffset(view);

    size_t pointLen = m_lasHeader.pointLen();

    // Make a buffer of at most a meg.
    std::vector<char> buf(std::min((size_t)1000000, pointLen * view->size()));

    const PointView& viewRef(*view.get());

    //ABELL - Removed callback handling for now.
    point_count_t remaining = view->size();
    PointId idx = 0;
    while (remaining)
    {
        point_count_t filled = fillWriteBuf(viewRef, idx, buf);
        idx += filled;
        remaining -= filled;

#ifdef PDAL_HAVE_LASZIP
        if (m_lasHeader.compressed())
        {
            char *pos = buf.data();
            for (point_count_t i = 0; i < filled; i++)
            {
                memcpy(m_zipPoint->m_lz_point_data.data(), pos, pointLen);
                if (!m_zipper->write(m_zipPoint->m_lz_point))
                {
                    std::ostringstream oss;
                    const char* err = m_zipper->get_error();
                    if (err == NULL)
                        err = "(unknown error)";
                    oss << "Error writing point: " << std::string(err);
                    throw pdal_error(oss.str());
                }
                pos += pointLen;
            }
        }
        else
            m_ostream->write(buf.data(), filled * pointLen);
#else
        m_ostream->write(buf.data(), filled * pointLen);
#endif
    }

    m_numPointsWritten = view->size() - remaining;
}

void LasWriter::writeView(const PointViewPtr view)
{
    Utils::writeProgress(m_progressFd, "READYVIEW",
        std::to_string(view->size()));
    m_scaling.setAutoXForm(view);

    point_count_t pointLen = m_lasHeader.pointLen();

    // Since we use the LASzip API, we can't benefit from building
    // a buffer of multiple points, so loop.
    if (m_compression == LasCompression::LasZip)
    {
        PointRef point(*view, 0);
        for (PointId idx = 0; idx < view->size(); ++idx)
        {
            point.setPointId(idx);
            processPoint(point);
        }
    }
    else
    {
        // Make a buffer of at most a meg.
        m_pointBuf.resize((std::min)((point_count_t)1000000,
                    pointLen * view->size()));

        const PointView& viewRef(*view.get());

        point_count_t remaining = view->size();
        PointId idx = 0;
        while (remaining)
        {
            point_count_t filled = fillWriteBuf(viewRef, idx, m_pointBuf);
            idx += filled;
            remaining -= filled;

            if (m_compression == LasCompression::LazPerf)
                writeLazPerfBuf(m_pointBuf.data(), pointLen, filled);
            else
                m_ostream->write(m_pointBuf.data(), filled * pointLen);
        }
    }
    Utils::writeProgress(m_progressFd, "DONEVIEW",
        std::to_string(view->size()));
}

point_count_t PgReader::readPgPatch(PointViewPtr view, point_count_t numPts)
{
    point_count_t numRemaining = m_patch.remaining;
    PointId nextId = view->size();
    point_count_t numRead = 0;

    size_t offset = (m_patch.count - m_patch.remaining) * packedPointSize();
    char *pos = (char *)(m_patch.binary.data() + offset);

    while (numRead < numPts && numRemaining > 0)
    {
        writePoint(*view.get(), nextId, pos);
        pos += packedPointSize();
        numRemaining--;
        nextId++;
        numRead++;
    }
    m_patch.remaining = numRemaining;
    return numRead;
}

PointViewSet GroupByFilter::run(PointViewPtr inView)
{
    PointViewSet viewSet;
    if (!inView->size())
        return viewSet;

    for (PointId idx = 0; idx < inView->size(); idx++)
    {
        uint64_t val = inView->getFieldAs<uint64_t>(m_dimId, idx);
        PointViewPtr& outView = m_viewMap[val];
        if (!outView)
            outView = inView->makeNew();
        outView->appendPoint(*inView.get(), idx);
    }

    // Pull the buffers out of the map and stick them in the standard
    // output set.
    for (auto bi = m_viewMap.begin(); bi != m_viewMap.end(); ++bi)
        viewSet.insert(bi->second);
    return viewSet;
}

// Test reprojecting UTM 15 to DD with a filter
TEST(ReprojectionFilterTest, ReprojectionFilterTest_test_1)
{
    const char* epsg4326_wkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]";

    PointTable table;

    const double postX = -93.351563;
    const double postY = 41.577148;
    const double postZ = 16.000000;

    {
        const SpatialReference out_ref(epsg4326_wkt);

        Options ops1;
        ops1.add("filename", Support::datapath("las/utm15.las"));
        LasReader reader;
        reader.setOptions(ops1);

        Options options;
        options.add("out_srs", out_ref.getWKT());

        ReprojectionFilter reprojectionFilter;
        reprojectionFilter.setOptions(options);
        reprojectionFilter.setInput(reader);

        reprojectionFilter.prepare(table);
        PointViewSet viewSet = reprojectionFilter.execute(table);
        EXPECT_EQ(viewSet.size(), 1u);
        PointViewPtr view = *viewSet.begin();

        double x, y, z;
        getPoint(*view.get(), x, y, z);

        EXPECT_FLOAT_EQ(x, postX);
        EXPECT_FLOAT_EQ(y, postY);
        EXPECT_FLOAT_EQ(z, postZ);
    }
}

MetadataNode InfoKernel::dumpPoints(PointViewPtr inView) const
{
    MetadataNode root;
    PointViewPtr outView = inView->makeNew();

    // Stick points in a inViewfer.
    std::vector<PointId> points = getListOfPoints(m_pointIndexes);
    for (size_t i = 0; i < points.size(); ++i)
    {
        PointId id = (PointId)points[i];
        if (id < inView->size())
            outView->appendPoint(*inView.get(), id);
    }

    MetadataNode tree = outView->toMetadata();
    std::string prefix("point ");
    for (size_t i = 0; i < outView->size(); ++i)
    {
        MetadataNode n = tree.findChild(std::to_string(i));
        n.add("PointId", points[i]);
        root.add(n.clone("point"));
    }
    return root;
}

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
//.........这里部分代码省略.........