|
1 | 1 | #include <data/VoxelGrid.h>
|
| 2 | + |
| 3 | +#include <iostream> |
| 4 | + |
| 5 | +void VoxelGrid::initialize(float resolution, const Eigen::Vector4f& min, const Eigen::Vector4f& max) { |
| 6 | + clear(); |
| 7 | + |
| 8 | + resolution_ = resolution; |
| 9 | + sizex_ = std::ceil((max.x() - min.x()) / resolution_); |
| 10 | + sizey_ = std::ceil((max.y() - min.y()) / resolution_); |
| 11 | + sizez_ = std::ceil((max.z() - min.z()) / resolution_); |
| 12 | + |
| 13 | + voxels_.resize(sizex_ * sizey_ * sizez_); |
| 14 | + // ensure that min, max are always inside the voxel grid. |
| 15 | + float ox = min.x() - 0.5 * (sizex_ * resolution - (max.x() - min.x())); |
| 16 | + float oy = min.y() - 0.5 * (sizey_ * resolution - (max.y() - min.y())); |
| 17 | + float oz = min.z() - 0.5 * (sizez_ * resolution - (max.z() - min.z())); |
| 18 | + offset_ = Eigen::Vector4f(ox, oy, oz, 1); |
| 19 | + |
| 20 | + // center_.head(3) = 0.5 * (max - min) + min; |
| 21 | + // center_[3] = 0; |
| 22 | + |
| 23 | + occlusions_.resize(sizex_ * sizey_ * sizez_); |
| 24 | +} |
| 25 | + |
| 26 | +void VoxelGrid::clear() { |
| 27 | + for (auto idx : occupied_) { |
| 28 | + voxels_[idx].count = 0; |
| 29 | + voxels_[idx].labels.clear(); |
| 30 | + } |
| 31 | + |
| 32 | + for (auto idx : occluded_) { |
| 33 | + voxels_[idx].count = 0; |
| 34 | + voxels_[idx].labels.clear(); |
| 35 | + } |
| 36 | + |
| 37 | + occupied_.clear(); |
| 38 | + occluded_.clear(); |
| 39 | +} |
| 40 | + |
| 41 | +void VoxelGrid::insert(const Eigen::Vector4f& p, uint32_t label) { |
| 42 | + Eigen::Vector4f tp = p - offset_; |
| 43 | + int32_t i = std::floor(tp.x() / resolution_); |
| 44 | + int32_t j = std::floor(tp.y() / resolution_); |
| 45 | + int32_t k = std::floor(tp.z() / resolution_); |
| 46 | + |
| 47 | + if ((i >= int32_t(sizex_)) || (j >= int32_t(sizey_)) || (k >= int32_t(sizez_))) return; |
| 48 | + if ((i < 0) || (j < 0) || (k < 0)) return; |
| 49 | + |
| 50 | + int32_t gidx = index(i, j, k); |
| 51 | + if (gidx < 0 || gidx >= int32_t(voxels_.size())) return; |
| 52 | + |
| 53 | + occupied_.push_back(gidx); |
| 54 | + |
| 55 | + // float n = voxels_[gidx].count; |
| 56 | + voxels_[gidx].labels[label] += 1; //(1. / (n + 1)) * (n * voxels_[gidx].point + p); |
| 57 | + voxels_[gidx].count += 1; |
| 58 | + occlusionsValid_ = false; |
| 59 | +} |
| 60 | + |
| 61 | +bool VoxelGrid::isOccluded(int32_t i, int32_t j, int32_t k) const { return occlusions_[index(i, j, k)] > -1; } |
| 62 | + |
| 63 | +bool VoxelGrid::isFree(int32_t i, int32_t j, int32_t k) const { return occlusions_[index(i, j, k)] == -1; } |
| 64 | + |
| 65 | +void VoxelGrid::insertOcclusionLabels() { |
| 66 | + if (!occlusionsValid_) updateOcclusions(); |
| 67 | + |
| 68 | + for (uint32_t i = 0; i < sizex_; ++i) { |
| 69 | + for (uint32_t j = 0; j < sizey_; ++j) { |
| 70 | + for (uint32_t k = 0; k < sizez_; ++k) { |
| 71 | + // heuristic: find label from above. |
| 72 | + if (occlusions_[index(i, j, k)] != index(i, j, k)) { |
| 73 | + int32_t n = 1; |
| 74 | + while ((k + n < sizez_) && isOccluded(i, j, k + n) && voxels_[index(i, j, k + n)].count == 0) n += 1; |
| 75 | + if (k + n < sizez_ && voxels_[index(i, j, k + n)].count > 0) { |
| 76 | + voxels_[index(i, j, k)].count = voxels_[index(i, j, k + n)].count; |
| 77 | + voxels_[index(i, j, k)].labels = voxels_[index(i, j, k + n)].labels; |
| 78 | + } |
| 79 | + } |
| 80 | + } |
| 81 | + } |
| 82 | + } |
| 83 | +} |
| 84 | + |
| 85 | +void VoxelGrid::updateOcclusions() { |
| 86 | + for (uint32_t i = 0; i < sizex_; ++i) { |
| 87 | + for (uint32_t j = 0; j < sizey_; ++j) { |
| 88 | + for (uint32_t k = 0; k < sizez_; ++k) { |
| 89 | + occlusions_[index(i, j, k)] = occludedBy(i, j, k); |
| 90 | + if (occlusions_[index(i, j, k)] != index(i, j, k)) occluded_.push_back(index(i, j, k)); |
| 91 | + } |
| 92 | + } |
| 93 | + } |
| 94 | + |
| 95 | + occlusionsValid_ = true; |
| 96 | +} |
| 97 | + |
| 98 | +int32_t VoxelGrid::occludedBy(int32_t i, int32_t j, int32_t k, std::vector<Eigen::Vector3i>* visited) const { |
| 99 | + float NextCrossingT[3], DeltaT[3]; /** t for next intersection with voxel boundary of axis, t increment for axis |
| 100 | + **/ |
| 101 | + int32_t Step[3], Out[3], Pos[3]; /** voxel increment for axis, index of of outside voxels, current position **/ |
| 102 | + float dir[3]; /** ray direction **/ |
| 103 | + |
| 104 | + if (visited != nullptr) visited->clear(); |
| 105 | + |
| 106 | + Pos[0] = i; |
| 107 | + Pos[1] = j; |
| 108 | + Pos[2] = k; |
| 109 | + |
| 110 | + /** calculate direction vector assuming sensor at (0,0,_heightOffset) **/ |
| 111 | + Eigen::Vector3f startpoint = voxel2position(Pos[0], Pos[1], Pos[2]); |
| 112 | + |
| 113 | + double halfResolution = 0.5 * resolution_; |
| 114 | + |
| 115 | + // startpoint[0] += halfResolution; |
| 116 | + // startpoint[1] += halfResolution; |
| 117 | + // startpoint[2] += halfResolution; |
| 118 | + |
| 119 | + dir[0] = -startpoint[0]; |
| 120 | + dir[1] = -startpoint[1]; |
| 121 | + dir[2] = -startpoint[2]; |
| 122 | + |
| 123 | + /** initialize variables for traversal **/ |
| 124 | + for (uint32_t axis = 0; axis < 3; ++axis) { |
| 125 | + if (dir[axis] < 0) { |
| 126 | + NextCrossingT[axis] = -halfResolution / dir[axis]; |
| 127 | + DeltaT[axis] = -resolution_ / dir[axis]; |
| 128 | + Step[axis] = -1; |
| 129 | + Out[axis] = 0; |
| 130 | + } else { |
| 131 | + NextCrossingT[axis] = halfResolution / dir[axis]; |
| 132 | + DeltaT[axis] = resolution_ / dir[axis]; |
| 133 | + Step[axis] = 1; |
| 134 | + Out[axis] = size(axis); |
| 135 | + } |
| 136 | + } |
| 137 | + |
| 138 | + Eigen::Vector3f endpoint(0, 0, 0); |
| 139 | + Eigen::Vector3i endindexes = position2voxel(endpoint); |
| 140 | + int32_t i_end = endindexes[0]; |
| 141 | + int32_t j_end = endindexes[1]; |
| 142 | + int32_t k_end = endindexes[2]; |
| 143 | + |
| 144 | + const int32_t cmpToAxis[8] = {2, 1, 2, 1, 2, 2, 0, 0}; |
| 145 | + int32_t iteration = 0; |
| 146 | + for (;;) // loop infinitely... |
| 147 | + { |
| 148 | + if (Pos[0] < 0 || Pos[1] < 0 || Pos[2] < 0) break; |
| 149 | + if (Pos[0] >= int32_t(sizex_) || Pos[1] >= int32_t(sizey_) || Pos[2] >= int32_t(sizez_)) break; |
| 150 | + |
| 151 | + int32_t idx = index(Pos[0], Pos[1], Pos[2]); |
| 152 | + bool occupied = voxels_[idx].count > 0; |
| 153 | + if (visited != nullptr) visited->push_back(Eigen::Vector3i(Pos[0], Pos[1], Pos[2])); |
| 154 | + |
| 155 | + if (occupied) return idx; |
| 156 | + |
| 157 | + int32_t bits = ((NextCrossingT[0] < NextCrossingT[1]) << 2) + ((NextCrossingT[0] < NextCrossingT[2]) << 1) + |
| 158 | + ((NextCrossingT[1] < NextCrossingT[2])); |
| 159 | + int32_t stepAxis = cmpToAxis[bits]; /* branch-free looping */ |
| 160 | + |
| 161 | + Pos[stepAxis] += Step[stepAxis]; |
| 162 | + NextCrossingT[stepAxis] += DeltaT[stepAxis]; |
| 163 | + |
| 164 | + /** note the first condition should never happen, since we want to reach a point inside the grid **/ |
| 165 | + if (Pos[stepAxis] == Out[stepAxis]) break; //... until outside, and leaving the loop here! |
| 166 | + if (Pos[0] == i_end && Pos[1] == j_end && Pos[2] == k_end) break; // .. or the sensor origin is reached. |
| 167 | + |
| 168 | + ++iteration; |
| 169 | + } |
| 170 | + |
| 171 | + return -1; |
| 172 | +} |
0 commit comments