Method setPlaneFit

filters/PlaneFitFilter.cpp:120–168

}

void PlaneFitFilter::setPlaneFit(PointView& view, const PointId& i)
{
    // Find k-nearest neighbors of i.
    const KD3Index& kdi = view.build3dIndex();
    PointIdList ni = kdi.neighbors(i, m_knn + 1);

    // Normal based only on neighbors, so exclude first point.
    PointIdList neighbors(ni.begin() + 1, ni.end());

    // Covariance and normal are based off demeaned coordinates, so we record
    // the centroid to properly offset the coordinates when computing point to
    // plance distance.
    Vector3d centroid = math::computeCentroid(view, neighbors);

    // Compute covariance of the neighbors.
    Matrix3d B = math::computeCovariance(view, neighbors);

    // Check if the covariance matrix is all zeros
    if (B.isZero())
    {
        log()->get(LogLevel::Info)
            << "Skipping point " << i
            << ". Covariance matrix is all zeros. This suggests a large "
               "number of redundant points. Consider using filters.sample "
               "with a small radius to remove redundant points.\n";
        return;
    }

    // Perform the eigen decomposition, using the eigenvector of the smallest
    // eigenvalue as the normal.
    Eigen::SelfAdjointEigenSolver<Matrix3d> solver(B);
    if (solver.info() != Eigen::Success)
        throwError("Cannot perform eigen decomposition.");
    Vector3d normal = solver.eigenvectors().col(0);

    // Compute point to plane distance of the query point.
    double d = absDistance(view, i, centroid, normal);

    // Compute mean point to plane distance of neighbors.
    double d_sum(0.0);
    for (PointId const& j : neighbors)
    {
        d_sum += absDistance(view, j, centroid, normal);
    }
    double d_bar(d_sum / m_knn);

    // Compute and set the plane fit criterion.
    view.setField(Id::PlaneFit, i, d / (d + d_bar));
}

} // namespace pdal