|
| 1 | +/* |
| 2 | + * Copyright (C) 2024, Graz University of Technology |
| 3 | + * This code is licensed under the MIT license (see LICENSE.txt in this folder for details) |
| 4 | + */ |
| 5 | + |
| 6 | + #pragma once |
| 7 | + |
| 8 | + #include "../auxiliary.h" |
| 9 | + |
| 10 | + #include <cooperative_groups.h> |
| 11 | + namespace cg = cooperative_groups; |
| 12 | + |
| 13 | +__device__ __inline__ uint64_t constructSortKey(uint32_t tile_id, float depth) |
| 14 | +{ |
| 15 | + uint64_t key = tile_id; |
| 16 | + key <<= 32; |
| 17 | + key |= *((uint32_t*)&depth); |
| 18 | + return key; |
| 19 | +} |
| 20 | + |
| 21 | +// Given a ray and a Gaussian primitive, compute the intersection depth. |
| 22 | +__device__ __inline__ bool getIntersectPoint( |
| 23 | + const int W, const int H, |
| 24 | + const float fx, const float fy, |
| 25 | + const float2 scale, |
| 26 | + const glm::vec2 pixel_center, |
| 27 | + const float* view2gaussian, |
| 28 | + float& depth |
| 29 | +){ |
| 30 | + |
| 31 | + // Fisrt compute two homogeneous planes, See Eq. (8) |
| 32 | + float3 Tu = {view2gaussian[0], view2gaussian[1], view2gaussian[2]}; |
| 33 | + float3 Tv = {view2gaussian[3], view2gaussian[4], view2gaussian[5]}; |
| 34 | + float3 Tw = {view2gaussian[6], view2gaussian[7], view2gaussian[8]}; |
| 35 | + float3 k = {-Tu.x + pixel_center.x * Tw.x, -Tu.y + pixel_center.x * Tw.y, -Tu.z + pixel_center.x * Tw.z}; |
| 36 | + float3 l = {-Tv.x + pixel_center.y * Tw.x, -Tv.y + pixel_center.y * Tw.y, -Tv.z + pixel_center.y * Tw.z}; |
| 37 | + // cross product of two planes is a line (i.e., homogeneous point), See Eq. (10) |
| 38 | + float3 p = crossProduct(k, l); |
| 39 | + |
| 40 | + if (p.z == 0.0) return false; // there is not intersection |
| 41 | + // TODO: no intersection if distance < scale |
| 42 | + |
| 43 | + // 3d homogeneous point to 2d point on the splat |
| 44 | + float2 s = {p.x / p.z, p.y / p.z}; |
| 45 | + // 3d distance. Compute Mahalanobis distance in the canonical splat' space |
| 46 | + float rho3d = (s.x * s.x + s.y * s.y); |
| 47 | + |
| 48 | + depth = (s.x * Tw.x + s.y * Tw.y) + Tw.z; // splat depth |
| 49 | + return true; |
| 50 | +} |
| 51 | + |
| 52 | + |
| 53 | +template<bool TILE_BASED_CULLING = false, bool LOAD_BALANCING = true> |
| 54 | +__global__ void duplicateWithKeys_extended( |
| 55 | + int P, |
| 56 | + int W, int H, |
| 57 | + const float focal_x, const float focal_y, |
| 58 | + const float2* __restrict__ points_xy, |
| 59 | + const float* __restrict__ depths, |
| 60 | + const float2* __restrict__ scales, |
| 61 | + const float* __restrict__ view2gaussians, |
| 62 | + const uint32_t* __restrict__ offsets, |
| 63 | + const int* __restrict__ radii, |
| 64 | + const float2* __restrict__ rects, |
| 65 | + uint64_t* __restrict__ gaussian_keys_unsorted, |
| 66 | + uint32_t* __restrict__ gaussian_values_unsorted, |
| 67 | + dim3 grid) |
| 68 | +{ |
| 69 | + auto block = cg::this_thread_block(); |
| 70 | + auto warp = cg::tiled_partition<WARP_SIZE>(block); |
| 71 | + |
| 72 | + // Since the projection of the quadratic surface on the image is non-convex, |
| 73 | + // there is no explicit solution for computing the pixel with the maximum weight on the image, |
| 74 | + // and tile-based culling is not performed. |
| 75 | + constexpr bool EVAL_MAX_CONTRIB_POS = false; |
| 76 | + constexpr bool PER_TILE_DEPTH = true; |
| 77 | + |
| 78 | +#define RETURN_OR_INACTIVE() if constexpr(LOAD_BALANCING) { active = false; } else { return; } |
| 79 | + uint32_t idx = cg::this_grid().thread_rank(); |
| 80 | + bool active = true; |
| 81 | + if (idx >= P) { |
| 82 | + RETURN_OR_INACTIVE(); |
| 83 | + idx = P - 1; |
| 84 | + } |
| 85 | + |
| 86 | + const int radius = radii[idx]; |
| 87 | + if (radius <= 0) { |
| 88 | + RETURN_OR_INACTIVE(); |
| 89 | + } |
| 90 | + |
| 91 | + // If the thread exceeds the Gaussian index, the Gaussian projection is zero, |
| 92 | + // and there are no Gaussians to process in the current warp, return. |
| 93 | + if constexpr(LOAD_BALANCING) |
| 94 | + if (__ballot_sync(WARP_MASK, active) == 0) |
| 95 | + return; |
| 96 | + |
| 97 | + // Find this Gaussian's offset in buffer for writing keys/values. |
| 98 | + uint32_t off_init = (idx == 0) ? 0 : offsets[idx - 1]; |
| 99 | + |
| 100 | + const int offset_to_init = offsets[idx]; |
| 101 | + const float global_depth_init = depths[idx]; |
| 102 | + |
| 103 | + const float2 xy_init = points_xy[idx]; |
| 104 | + const float2 rect_dims_init = rects[idx]; |
| 105 | + |
| 106 | + __shared__ float2 s_xy[BLOCK_SIZE]; |
| 107 | + __shared__ float2 s_rect_dims[BLOCK_SIZE]; |
| 108 | + __shared__ float s_radius[BLOCK_SIZE]; |
| 109 | + s_xy[block.thread_rank()] = xy_init; |
| 110 | + s_rect_dims[block.thread_rank()] = rect_dims_init; |
| 111 | + s_radius[block.thread_rank()] = radius; |
| 112 | + |
| 113 | + uint2 rect_min_init, rect_max_init; |
| 114 | +#if FAST_INFERENCE |
| 115 | + if (radius > MAX_BILLBOARD_SIZE) |
| 116 | + getRectOld(xy_init, radius, rect_min_init, rect_max_init, grid); |
| 117 | + else |
| 118 | + getRect(xy_init, rect_dims_init, rect_min_init, rect_max_init, grid); |
| 119 | +# else |
| 120 | + getRectOld(xy_init, radius, rect_min_init, rect_max_init, grid); |
| 121 | +#endif |
| 122 | + |
| 123 | + __shared__ float s_view2gaussians[BLOCK_SIZE * 9]; |
| 124 | + __shared__ float2 s_scales[BLOCK_SIZE]; |
| 125 | + |
| 126 | + if (PER_TILE_DEPTH) |
| 127 | + { |
| 128 | + s_scales[block.thread_rank()] = scales[idx]; |
| 129 | + for (int ii = 0; ii < 9; ii++) |
| 130 | + s_view2gaussians[9 * block.thread_rank() + ii] = view2gaussians[idx * 9 + ii]; |
| 131 | + } |
| 132 | + |
| 133 | + constexpr uint32_t SEQUENTIAL_TILE_THRESH = 32U; // all tiles above this threshold will be computed cooperatively |
| 134 | + const uint32_t rect_width_init = (rect_max_init.x - rect_min_init.x); |
| 135 | + const uint32_t tile_count_init = (rect_max_init.y - rect_min_init.y) * rect_width_init; |
| 136 | + |
| 137 | + // Generate no key/value pair for invisible Gaussians |
| 138 | + if (tile_count_init == 0) { |
| 139 | + RETURN_OR_INACTIVE(); |
| 140 | + } |
| 141 | + auto tile_function = [&](const int W, const int H, |
| 142 | + const float fx, const float fy, |
| 143 | + float2 xy, |
| 144 | + int x, int y,// tile ID |
| 145 | + const float2 scale, |
| 146 | + const float* view2gaussian, |
| 147 | + const float global_depth, |
| 148 | + float& depth) |
| 149 | + { |
| 150 | + const glm::vec2 tile_min(x * BLOCK_X, y * BLOCK_Y); |
| 151 | + const glm::vec2 tile_max((x + 1) * BLOCK_X - 1, (y + 1) * BLOCK_Y - 1); // 像素坐标 |
| 152 | + |
| 153 | + glm::vec2 max_pos; |
| 154 | + if constexpr (PER_TILE_DEPTH) |
| 155 | + { |
| 156 | + glm::vec2 target_pos = {max(min(xy.x, tile_max.x), tile_min.x), max(min(xy.y, tile_max.y), tile_min.y)}; |
| 157 | + |
| 158 | + // Or select the tile's center pixel as the target_pos. |
| 159 | + // const glm::vec2 tile_center = (tile_min + tile_max) * 0.5f; |
| 160 | + // glm::vec2 target_pos = tile_center; |
| 161 | + |
| 162 | + bool intersect = getIntersectPoint( |
| 163 | + W, H, fx, fy, scale, target_pos, view2gaussian, depth); // Compute the intersection point of the quadratic surface. |
| 164 | + if (intersect) |
| 165 | + depth = max(0.0f, depth); |
| 166 | + else // If there is no intersection, sort by the Gaussian centroid. |
| 167 | + depth = global_depth; |
| 168 | + } |
| 169 | + else |
| 170 | + { |
| 171 | + depth = global_depth; |
| 172 | + } |
| 173 | + |
| 174 | + // Since the quadratic surface is non-convex, tile-based culling is not performed. |
| 175 | + // return (!TILE_BASED_CULLING) || max_opac_factor <= opacity_factor_threshold; |
| 176 | + return true; |
| 177 | + }; |
| 178 | + |
| 179 | + if (active) |
| 180 | + { |
| 181 | + const float2 scale_init = { |
| 182 | + s_scales[block.thread_rank()].x, |
| 183 | + s_scales[block.thread_rank()].y}; |
| 184 | + |
| 185 | + float view2gaussian_init[9]; |
| 186 | + for (int ii = 0; ii < 9; ii++) |
| 187 | + view2gaussian_init[ii] = s_view2gaussians[9 * block.thread_rank() + ii]; |
| 188 | + |
| 189 | + for (uint32_t tile_idx = 0; tile_idx < tile_count_init && (!LOAD_BALANCING || tile_idx < SEQUENTIAL_TILE_THRESH); tile_idx++) |
| 190 | + { |
| 191 | + const int y = (tile_idx / rect_width_init) + rect_min_init.y; |
| 192 | + const int x = (tile_idx % rect_width_init) + rect_min_init.x; |
| 193 | + |
| 194 | + float depth; |
| 195 | + bool write_tile = tile_function( |
| 196 | + W, H, focal_x, focal_y, |
| 197 | + xy_init, x, y, scale_init, view2gaussian_init, global_depth_init, depth); |
| 198 | + if (write_tile) |
| 199 | + { |
| 200 | + if (off_init < offset_to_init) |
| 201 | + { |
| 202 | + const uint32_t tile_id = y * grid.x + x; |
| 203 | + gaussian_values_unsorted[off_init] = idx; |
| 204 | + gaussian_keys_unsorted[off_init] = constructSortKey(tile_id, depth); |
| 205 | + } |
| 206 | + else |
| 207 | + { |
| 208 | +#ifdef DUPLICATE_OPT_DEBUG |
| 209 | + printf("Error (sequential): Too little memory reserved in preprocess: off=%d off_to=%d idx=%d\n", off_init, offset_to_init, idx); |
| 210 | +#endif |
| 211 | + } |
| 212 | + off_init++; |
| 213 | + } |
| 214 | + } |
| 215 | + } |
| 216 | + |
| 217 | +#undef RETURN_OR_INACTIVE |
| 218 | + |
| 219 | + if (!LOAD_BALANCING) // Coordinate to handle the unprocessed tasks of other threads within the same warp. |
| 220 | + return; |
| 221 | + |
| 222 | + const uint32_t idx_init = idx; // Current thread idx. |
| 223 | + const uint32_t lane_idx = cg::this_thread_block().thread_rank() % WARP_SIZE; |
| 224 | + const uint32_t warp_idx = cg::this_thread_block().thread_rank() / WARP_SIZE; |
| 225 | + unsigned int lane_mask_allprev_excl = 0xFFFFFFFFU >> (WARP_SIZE - lane_idx); |
| 226 | + |
| 227 | + const int32_t compute_cooperatively = active && tile_count_init > SEQUENTIAL_TILE_THRESH; // Determine whether additional idle threads are needed for computation. |
| 228 | + const uint32_t remaining_threads = __ballot_sync(WARP_MASK, compute_cooperatively); |
| 229 | + if (remaining_threads == 0) |
| 230 | + return; |
| 231 | + |
| 232 | + uint32_t n_remaining_threads = __popc(remaining_threads); // The number of threads required for collaborative computation. |
| 233 | + for (int n = 0; n < n_remaining_threads && n < WARP_SIZE; n++) |
| 234 | + { |
| 235 | + int i = __fns(remaining_threads, 0, n+1); // find lane index of next remaining thread |
| 236 | + |
| 237 | + uint32_t idx_coop = __shfl_sync(WARP_MASK, idx_init, i); |
| 238 | + uint32_t off_coop = __shfl_sync(WARP_MASK, off_init, i); |
| 239 | + |
| 240 | + const uint32_t offset_to = __shfl_sync(WARP_MASK, offset_to_init, i); |
| 241 | + const float global_depth = __shfl_sync(WARP_MASK, global_depth_init, i); |
| 242 | + |
| 243 | + const float2 xy = s_xy[warp.meta_group_rank() * WARP_SIZE + i]; |
| 244 | + const float2 rect_dims = s_rect_dims[warp.meta_group_rank() * WARP_SIZE + i]; |
| 245 | + const float rad = s_radius[warp.meta_group_rank() * WARP_SIZE + i]; |
| 246 | + const float2 scale = { |
| 247 | + s_scales[warp.meta_group_rank() * WARP_SIZE + i].x, |
| 248 | + s_scales[warp.meta_group_rank() * WARP_SIZE + i].y}; |
| 249 | + float view2gaussian[9]; |
| 250 | + for (int ii = 0; ii < 9; ii++) |
| 251 | + view2gaussian[ii] = s_view2gaussians[9 * (warp.meta_group_rank() * WARP_SIZE + i) + ii]; |
| 252 | + |
| 253 | + uint2 rect_min, rect_max; |
| 254 | +#if FAST_INFERENCE |
| 255 | + if (radius > MAX_BILLBOARD_SIZE) |
| 256 | + getRectOld(xy, rad, rect_min, rect_max, grid); |
| 257 | + else |
| 258 | + getRect(xy, rect_dims, rect_min, rect_max, grid); |
| 259 | +#else |
| 260 | + getRectOld(xy, rad, rect_min, rect_max, grid); |
| 261 | +#endif |
| 262 | + |
| 263 | + const uint32_t rect_width = (rect_max.x - rect_min.x); |
| 264 | + const uint32_t tile_count = (rect_max.y - rect_min.y) * rect_width; |
| 265 | + const uint32_t remaining_tile_count = tile_count - SEQUENTIAL_TILE_THRESH; |
| 266 | + const int32_t n_iterations = (remaining_tile_count + WARP_SIZE - 1) / WARP_SIZE; |
| 267 | + for (int it = 0; it < n_iterations; it++) |
| 268 | + { |
| 269 | + int tile_idx = it * WARP_SIZE + lane_idx + SEQUENTIAL_TILE_THRESH; // it*32 + local_warp_idx + 32 |
| 270 | + int active_curr_it = tile_idx < tile_count; |
| 271 | + |
| 272 | + int y = (tile_idx / rect_width) + rect_min.y; |
| 273 | + int x = (tile_idx % rect_width) + rect_min.x; |
| 274 | + |
| 275 | + float depth; |
| 276 | + bool write_tile = tile_function( |
| 277 | + W, H, focal_x, focal_y, |
| 278 | + xy, x, y, scale, view2gaussian, global_depth, depth |
| 279 | + ); |
| 280 | + |
| 281 | + const uint32_t write = active_curr_it && write_tile; |
| 282 | + |
| 283 | + uint32_t n_writes, write_offset; |
| 284 | + if constexpr (!TILE_BASED_CULLING) |
| 285 | + { |
| 286 | + n_writes = WARP_SIZE; |
| 287 | + write_offset = off_coop + lane_idx; |
| 288 | + } |
| 289 | + else |
| 290 | + { |
| 291 | + const uint32_t write_ballot = __ballot_sync(WARP_MASK, write); |
| 292 | + n_writes = __popc(write_ballot); |
| 293 | + |
| 294 | + const uint32_t write_offset_it = __popc(write_ballot & lane_mask_allprev_excl); |
| 295 | + write_offset = off_coop + write_offset_it; |
| 296 | + } |
| 297 | + |
| 298 | + if (write) |
| 299 | + { |
| 300 | + if (write_offset < offset_to) |
| 301 | + { |
| 302 | + const uint32_t tile_id = y * grid.x + x; |
| 303 | + gaussian_values_unsorted[write_offset] = idx_coop; |
| 304 | + gaussian_keys_unsorted[write_offset] = constructSortKey(tile_id, depth); |
| 305 | + } |
| 306 | + #ifdef DUPLICATE_OPT_DEBUG |
| 307 | + else |
| 308 | + { |
| 309 | + printf("Error (parallel): Too little memory reserved in preprocess: off=%d off_to=%d idx=%d tile_count=%d it=%d | x=%d y=%d rect=(%d %d - %d %d)\n", |
| 310 | + write_offset, offset_to, idx_coop, tile_count, it, x, y, rect_min.x, rect_min.y, rect_max.x, rect_max.y); |
| 311 | + } |
| 312 | + #endif |
| 313 | + } |
| 314 | + off_coop += n_writes; |
| 315 | + } |
| 316 | + |
| 317 | + __syncwarp(); |
| 318 | + } |
| 319 | + } |
0 commit comments