| 113 | } |
| 114 | |
| 115 | void EigenvaluesFilter::filter(PointView& view) |
| 116 | { |
| 117 | const KD3Index& kdi = view.build3dIndex(); |
| 118 | |
| 119 | for (PointRef p : view) |
| 120 | { |
| 121 | // find neighbors, either by radius or k nearest neighbors |
| 122 | PointIdList ids; |
| 123 | if (m_args->m_radiusArg->set()) |
| 124 | { |
| 125 | ids = kdi.radius(p, m_args->m_radius); |
| 126 | |
| 127 | // if insufficient number of neighbors, eigen solver will fail |
| 128 | // anyway, it may be okay to silently return without setting any of |
| 129 | // the computed features? |
| 130 | if (ids.size() < (size_t)m_args->m_minK) { |
| 131 | log()->get(LogLevel::Info) |
| 132 | << "Skipping point " << p.pointId() << ". Found " << ids.size() |
| 133 | << " neighbors but required " << m_args->m_minK << ".\n"; |
| 134 | continue; |
| 135 | } |
| 136 | } |
| 137 | else |
| 138 | { |
| 139 | ids = kdi.neighbors(p, m_args->m_knn + 1, m_args->m_stride); |
| 140 | } |
| 141 | |
| 142 | // compute covariance of the neighborhood |
| 143 | Matrix3d B = math::computeCovariance(view, ids); |
| 144 | |
| 145 | // Check if the covariance matrix is all zeros |
| 146 | if (B.isZero()) |
| 147 | { |
| 148 | log()->get(LogLevel::Info) |
| 149 | << "Skipping point " << p.pointId() |
| 150 | << ". Covariance matrix is all zeros. This suggests a large " |
| 151 | "number of redundant points. Consider using filters.sample " |
| 152 | "with a small radius to remove redundant points.\n"; |
| 153 | continue; |
| 154 | } |
| 155 | |
| 156 | // perform the eigen decomposition |
| 157 | Eigen::SelfAdjointEigenSolver<Matrix3d> solver(B); |
| 158 | if (solver.info() != Eigen::Success) |
| 159 | throwError("Cannot perform eigen decomposition."); |
| 160 | Vector3d ev = solver.eigenvalues(); |
| 161 | |
| 162 | if (m_args->m_normalize) |
| 163 | { |
| 164 | double sum = ev[0] + ev[1] + ev[2]; |
| 165 | ev /= sum; |
| 166 | } |
| 167 | |
| 168 | p.setField(Id::Eigenvalue0, ev[0]); |
| 169 | p.setField(Id::Eigenvalue1, ev[1]); |
| 170 | p.setField(Id::Eigenvalue2, ev[2]); |
| 171 | } |
| 172 | } |