diff options
Diffstat (limited to 'source/blender/nodes/geometry/nodes/node_geo_point_distribute_poisson_disk.cc')
| -rw-r--r-- | source/blender/nodes/geometry/nodes/node_geo_point_distribute_poisson_disk.cc | 281 |
1 files changed, 281 insertions, 0 deletions
diff --git a/source/blender/nodes/geometry/nodes/node_geo_point_distribute_poisson_disk.cc b/source/blender/nodes/geometry/nodes/node_geo_point_distribute_poisson_disk.cc new file mode 100644 index 00000000000..47764efa15d --- /dev/null +++ b/source/blender/nodes/geometry/nodes/node_geo_point_distribute_poisson_disk.cc @@ -0,0 +1,281 @@ +/* + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* + * Based on Cem Yuksel. 2015. Sample Elimination for Generating Poisson Disk Sample + * ! Sets. Computer Graphics Forum 34, 2 (May 2015), 25-32. + * ! http://www.cemyuksel.com/research/sampleelimination/ + * Copyright (c) 2016, Cem Yuksel <cem@cemyuksel.com> + * All rights reserved. + */ + +#include "BLI_inplace_priority_queue.hh" +#include "BLI_kdtree.h" + +#include "node_geometry_util.hh" + +#include <iostream> +#include <string.h> + +namespace blender::nodes { + +static void tile_point(Vector<float3> *tiled_points, + Vector<size_t> *indices, + const float maximum_distance, + const float3 boundbox, + float3 const &point, + size_t index, + int dimension = 0) +{ + for (int dimension_iter = dimension; dimension_iter < 3; dimension_iter++) { + if (boundbox[dimension_iter] - point[dimension_iter] < maximum_distance) { + float3 point_tiled = point; + point_tiled[dimension_iter] -= boundbox[dimension_iter]; + + tiled_points->append(point_tiled); + indices->append(index); + + tile_point(tiled_points, + indices, + maximum_distance, + boundbox, + point_tiled, + index, + dimension_iter + 1); + } + + if (point[dimension_iter] < maximum_distance) { + float3 point_tiled = point; + point_tiled[dimension_iter] += boundbox[dimension_iter]; + + tiled_points->append(point_tiled); + indices->append(index); + + tile_point(tiled_points, + indices, + maximum_distance, + boundbox, + point_tiled, + index, + dimension_iter + 1); + } + } +} + +/** + * Returns the weight the point gets based on the distance to another point. + */ +static float point_weight_influence_get(const float maximum_distance, + const float minimum_distance, + float distance) +{ + const float alpha = 8.0f; + + if (distance < minimum_distance) { + distance = minimum_distance; + } + + return std::pow(1.0f - distance / maximum_distance, alpha); +} + +/** + * Weight each point based on their proximity to its neighbors + * + * For each index in the weight array add a weight based on the proximity the + * corresponding point has with its neighboors. + **/ +static void points_distance_weight_calculate(Vector<float> *weights, + const size_t point_id, + const float3 *input_points, + const void *kd_tree, + const float minimum_distance, + const float maximum_distance, + InplacePriorityQueue<float> *heap) +{ + KDTreeNearest_3d *nearest_point = nullptr; + int neighbors = BLI_kdtree_3d_range_search( + (KDTree_3d *)kd_tree, input_points[point_id], &nearest_point, maximum_distance); + + for (int i = 0; i < neighbors; i++) { + size_t neighbor_point_id = nearest_point[i].index; + + if (neighbor_point_id >= weights->size()) { + continue; + } + + /* The point should not influence itself. */ + if (neighbor_point_id == point_id) { + continue; + } + + const float weight_influence = point_weight_influence_get( + maximum_distance, minimum_distance, nearest_point[i].dist); + + /* In the first pass we just the weights. */ + if (heap == nullptr) { + (*weights)[point_id] += weight_influence; + } + /* When we run again we need to update the weights and the heap. */ + else { + (*weights)[neighbor_point_id] -= weight_influence; + heap->priority_decreased(neighbor_point_id); + } + } + + if (nearest_point) { + MEM_freeN(nearest_point); + } +} + +/** + * Returns the minimum radius fraction used by the default weight function. + */ +static float weight_limit_fraction_get(const size_t input_size, const size_t output_size) +{ + const float beta = 0.65f; + const float gamma = 1.5f; + float ratio = float(output_size) / float(input_size); + return (1.0f - std::pow(ratio, gamma)) * beta; +} + +/** + * Tile the input points. + */ +static void points_tiling(const float3 *input_points, + const size_t input_size, + void **kd_tree, + const float maximum_distance, + const float3 boundbox) + +{ + Vector<float3> tiled_points(input_points, input_points + input_size); + Vector<size_t> indices(input_size); + + for (size_t i = 0; i < input_size; i++) { + indices[i] = i; + } + + /* Tile the tree based on the boundbox. */ + for (size_t i = 0; i < input_size; i++) { + tile_point(&tiled_points, &indices, maximum_distance, boundbox, input_points[i], i); + } + + /* Build a new tree with the new indices and tiled points. */ + *kd_tree = BLI_kdtree_3d_new(tiled_points.size()); + for (size_t i = 0; i < tiled_points.size(); i++) { + BLI_kdtree_3d_insert(*(KDTree_3d **)kd_tree, indices[i], tiled_points[i]); + } + BLI_kdtree_3d_balance(*(KDTree_3d **)kd_tree); +} + +static void weighted_sample_elimination(const float3 *input_points, + const size_t input_size, + float3 *output_points, + const size_t output_size, + const float maximum_distance, + const float3 boundbox, + const bool do_copy_eliminated) +{ + const float minimum_distance = maximum_distance * + weight_limit_fraction_get(input_size, output_size); + + void *kd_tree = nullptr; + points_tiling(input_points, input_size, &kd_tree, maximum_distance, boundbox); + + /* Assign weights to each sample. */ + Vector<float> weights(input_size, 0.0f); + for (size_t point_id = 0; point_id < weights.size(); point_id++) { + points_distance_weight_calculate( + &weights, point_id, input_points, kd_tree, minimum_distance, maximum_distance, nullptr); + } + + /* Remove the points based on their weight. */ + InplacePriorityQueue<float> heap(weights); + + size_t sample_size = input_size; + while (sample_size > output_size) { + /* For each sample around it, remove its weight contribution and update the heap. */ + size_t point_id = heap.pop_index(); + points_distance_weight_calculate( + &weights, point_id, input_points, kd_tree, minimum_distance, maximum_distance, &heap); + sample_size--; + } + + /* Copy the samples to the output array. */ + size_t target_size = do_copy_eliminated ? input_size : output_size; + for (size_t i = 0; i < target_size; i++) { + size_t index = heap.all_indices()[i]; + output_points[i] = input_points[index]; + } + + /* Cleanup. */ + BLI_kdtree_3d_free((KDTree_3d *)kd_tree); +} + +static void progressive_sampling_reorder(Vector<float3> *output_points, + float maximum_density, + float3 boundbox) +{ + /* Re-order the points for progressive sampling. */ + Vector<float3> temporary_points(output_points->size()); + float3 *source_points = output_points->data(); + float3 *dest_points = temporary_points.data(); + size_t source_size = output_points->size(); + size_t dest_size = 0; + + while (source_size >= 3) { + dest_size = source_size * 0.5f; + + /* Changes the weight function radius using half of the number of samples. + * It is used for progressive sampling. */ + maximum_density *= std::sqrt(2.0f); + weighted_sample_elimination( + source_points, source_size, dest_points, dest_size, maximum_density, boundbox, true); + + if (dest_points != output_points->data()) { + memcpy((*output_points)[dest_size], + dest_points[dest_size], + (source_size - dest_size) * sizeof(float3)); + } + + /* Swap the arrays around. */ + float3 *points_iter = source_points; + source_points = dest_points; + dest_points = points_iter; + source_size = dest_size; + } + if (source_points != output_points->data()) { + memcpy(output_points->data(), source_points, dest_size * sizeof(float3)); + } +} + +void poisson_disk_point_elimination(Vector<float3> const *input_points, + Vector<float3> *output_points, + float maximum_density, + float3 boundbox) +{ + weighted_sample_elimination(input_points->data(), + input_points->size(), + output_points->data(), + output_points->size(), + maximum_density, + boundbox, + false); + + progressive_sampling_reorder(output_points, maximum_density, boundbox); +} + +} // namespace blender::nodes |