| 66 | } |
| 67 | |
| 68 | PointViewSet VoxelCentroidNearestNeighborFilter::run(PointViewPtr view) |
| 69 | { |
| 70 | typedef std::make_signed<std::size_t>::type ssize_t; |
| 71 | double x0 = view->getFieldAs<double>(Dimension::Id::X, 0); |
| 72 | double y0 = view->getFieldAs<double>(Dimension::Id::Y, 0); |
| 73 | double z0 = view->getFieldAs<double>(Dimension::Id::Z, 0); |
| 74 | |
| 75 | // Make an initial pass through the input PointView to index PointIds by |
| 76 | // row, column, and depth. |
| 77 | std::map<std::tuple<ssize_t, ssize_t, ssize_t>, PointIdList> |
| 78 | populated_voxel_ids; |
| 79 | for (PointId id = 0; id < view->size(); ++id) |
| 80 | { |
| 81 | double y = view->getFieldAs<double>(Dimension::Id::Y, id); |
| 82 | double x = view->getFieldAs<double>(Dimension::Id::X, id); |
| 83 | double z = view->getFieldAs<double>(Dimension::Id::Z, id); |
| 84 | ssize_t r = static_cast<ssize_t>(std::floor((y - y0) / m_cell)); |
| 85 | ssize_t c = static_cast<ssize_t>(std::floor((x - x0) / m_cell)); |
| 86 | ssize_t d = static_cast<ssize_t>(std::floor((z - z0) / m_cell)); |
| 87 | populated_voxel_ids[std::make_tuple(r, c, d)].push_back(id); |
| 88 | } |
| 89 | |
| 90 | // Make a second pass through the populated voxels to compute the voxel |
| 91 | // centroid and to find its nearest neighbor. |
| 92 | PointViewPtr output = view->makeNew(); |
| 93 | for (auto const& t : populated_voxel_ids) |
| 94 | { |
| 95 | if (t.second.size() == 1) |
| 96 | { |
| 97 | // If there is only one point in the voxel, simply append it. |
| 98 | output->appendPoint(*view, t.second[0]); |
| 99 | } |
| 100 | else if (t.second.size() == 2) |
| 101 | { |
| 102 | // Else if there are only two, they are equidistant to the |
| 103 | // centroid, so append the one closest to voxel center. |
| 104 | |
| 105 | // Compute voxel center. |
| 106 | double y_center = y0 + (std::get<0>(t.first) + 0.5) * m_cell; |
| 107 | double x_center = x0 + (std::get<1>(t.first) + 0.5) * m_cell; |
| 108 | double z_center = z0 + (std::get<2>(t.first) + 0.5) * m_cell; |
| 109 | |
| 110 | // Compute distance from first point to voxel center. |
| 111 | double x1 = view->getFieldAs<double>(Dimension::Id::X, t.second[0]); |
| 112 | double y1 = view->getFieldAs<double>(Dimension::Id::Y, t.second[0]); |
| 113 | double z1 = view->getFieldAs<double>(Dimension::Id::Z, t.second[0]); |
| 114 | double d1 = pow(x_center - x1, 2) + pow(y_center - y1, 2) + pow(z_center - z1, 2); |
| 115 | |
| 116 | // Compute distance from second point to voxel center. |
| 117 | double x2 = view->getFieldAs<double>(Dimension::Id::X, t.second[1]); |
| 118 | double y2 = view->getFieldAs<double>(Dimension::Id::Y, t.second[1]); |
| 119 | double z2 = view->getFieldAs<double>(Dimension::Id::Z, t.second[1]); |
| 120 | double d2 = pow(x_center - x2, 2) + pow(y_center - y2, 2) + pow(z_center - z2, 2); |
| 121 | |
| 122 | // Append the closer of the two. |
| 123 | if (d1 < d2) |
| 124 | output->appendPoint(*view, t.second[0]); |
| 125 | else |
nothing calls this directly
no test coverage detected