diff --git a/.github/OcMesher.png b/.github/OcMesher.png new file mode 100644 index 0000000..75cf8a5 Binary files /dev/null and b/.github/OcMesher.png differ diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a80cafa --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +__pycache__ +lib +results \ No newline at end of file diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..6882739 --- /dev/null +++ b/LICENSE @@ -0,0 +1,28 @@ +BSD 3-Clause License + +Copyright (c) 2023, Princeton Vision & Learning Lab + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 index 0000000..8044322 --- /dev/null +++ b/README.md @@ -0,0 +1,53 @@ +## OcMesher: View-Dependent Octree-based Mesh Extraction in Unbounded Scenes for Procedural Synthetic Data + +Implementation source-code for OcMesher, which extracts a mesh for an unbounded scene represented by signed distance functions (SDFs). Even though the scene is unbounded, the mesh is memory-efficient, and highly detailed from a given set of camera views. OcMesher is used by default in [Infinigen](https://github.com/princeton-vl/infinigen) to speed up video generation, improve rendering quality, and export terrain meshes to external simulators. + + + +If you use OcMesher in your work, please cite our academic paper: + +

+ + View-Dependent Octree-based Mesh Extraction in Unbounded Scenes for Procedural Synthetic Data + +

+

+ Zeyu Ma, + Alexander Raistrick, + Lahav Lipson, + Jia Deng
+

+ +``` +@article{ocmesher2023view, + title={View-Dependent Octree-based Mesh Extraction in Unbounded Scenes for Procedural Synthetic Data}, + author={Ma, Zeyu and Raistrick, Alexander and Lipson, Lahav and Deng, Jia}, + year={2023} +} +``` + +Please view the video [here](https://youtu.be/YA1c5L0Ncuw) for more qualitative results + +## Getting Started + +:bulb: Note: OcMesher is installed by default in Infinigen as of v1.2.0 - if you wish to use OcMesher with Infinigen please follow the Installation instructions on the infinigen repo. Use the instructions below only if you want a standalone installation & demo. + +### Standalone Installation + +``` +git clone https://github.com/princeton-vl/OcMesher.git +cd OcMesher +bash install.sh +conda create --name ocmesher python=3.10 +conda activate ocmesher +pip install -r requirements.txt +``` + + +### Demo + +``` +python demo.py +``` + +This example uses one camera and the Perlin Noise from the Python library `vnoise` and outputs the resulting mesh in `results/demo.obj`. diff --git a/__init__.py b/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/demo.py b/demo.py new file mode 100644 index 0000000..13da5b1 --- /dev/null +++ b/demo.py @@ -0,0 +1,38 @@ +# Copyright (c) Princeton University. +# This source code is licensed under the BSD 3-Clause license found in the LICENSE file in the root directory of this source tree. + +# Authors: Zeyu Ma + +import os +import sys +import numpy as np +from ocmesher import OcMesher + +import vnoise +noise = vnoise.Noise() + +def f(XYZ): + scale = 2 + h = noise.noise2(XYZ[:, 0] / scale, XYZ[:, 1] / scale, grid_mode=False, octaves=4) + + return XYZ[:, 2] - h + + +cam_poses = [np.array([ + [1, 0, 0, 0], + [0, 0, 1, 0], + [0, -1, 0, 3], + [0, 0, 0, 1], +])] +Ks = [np.array([ + [2000, 0, 640], + [0, 2000, 360], + [0, 0, 1] +])] +Hs = [720] +Ws = [1280] + +mesher = OcMesher((cam_poses, Ks, Hs, Ws), pixels_per_cube=16) +meshes, in_view_tags = mesher([f]) +os.makedirs("results", exist_ok=True) +meshes[0].export("results/demo.obj") \ No newline at end of file diff --git a/install.sh b/install.sh new file mode 100644 index 0000000..cb1886c --- /dev/null +++ b/install.sh @@ -0,0 +1,27 @@ +#!/bin/bash + +shopt -s expand_aliases +set -e + +OS=$(uname -s) +ARCH=$(uname -m) + +if [ "${OS}" = "Linux" ]; then + alias gx1="g++ -O3 -c -fpic -fopenmp " + alias gx2="g++ -O3 -shared -fopenmp " +elif [ "${OS}" = "Darwin" ]; then + if [ "${ARCH}" = "arm64" ]; then + compiler="/opt/homebrew/opt/llvm/bin/clang++" + else + compiler="/usr/local/opt/llvm/bin/clang++" + fi + alias gx1="${compiler} -O3 -c -fpic -fopenmp " + alias gx2="${compiler} -O3 -shared -fopenmp " +else + echo "Unsupported OS" + exit -1 +fi + +mkdir -p ocmesher/lib +gx1 -o ocmesher/lib/core.o ocmesher/source/core.cpp +gx2 -o ocmesher/lib/core.so ocmesher/lib/core.o diff --git a/ocmesher/__init__.py b/ocmesher/__init__.py new file mode 100644 index 0000000..2fda981 --- /dev/null +++ b/ocmesher/__init__.py @@ -0,0 +1 @@ +from .core import OcMesher \ No newline at end of file diff --git a/ocmesher/core.py b/ocmesher/core.py new file mode 100644 index 0000000..d0c3769 --- /dev/null +++ b/ocmesher/core.py @@ -0,0 +1,230 @@ +# Copyright (c) Princeton University. +# This source code is licensed under the BSD 3-Clause license found in the LICENSE file in the root directory of this source tree. + +# Authors: Zeyu Ma + +from pathlib import Path +import sys +import gin +import numpy as np +import trimesh +from tqdm import tqdm + +from .utils.interface import AC, POINTER, AsDouble, AsFloat, AsInt, AsBool, c_bool, c_double, c_float, c_int32, load_cdll, register_func +from .utils.timer import Timer + +@gin.configurable +class OcMesher: + def __init__(self, + cameras, + pixels_per_cube=8, + inv_scale=10, + min_dist=1, + center=(0, 0, 0), size=1e3, + memory_limit_mb=1000, + bisection_iters=15, + enclosed=True, + simplify_occluded=True, + visible_relax_iter=2, + coarse_count=500000, + ): + dll = load_cdll(str(Path(__file__).parent.resolve()/"lib"/"core.so")) + self.float_type = c_double + self.np_float_type = np.float64 + self.AF = AsDouble + self.sdf_float_type = c_float + self.sdf_np_float_type = np.float32 + self.sdf_AF = AsFloat + + self.memory_limit_mb = memory_limit_mb + + cam_poses, Ks, Hs, Ws = cameras + self.n_cameras = len(cam_poses) + self.cameras = np.zeros(23 * self.n_cameras, dtype=self.np_float_type) + for i in range(self.n_cameras): + self.cameras[23 * i: 23 * (i+1)] = np.concatenate([ + np.linalg.inv(cam_poses[i])[:3, :4].reshape(-1), + Ks[i].reshape(-1), [Hs[i]], [Ws[i]] + ]).astype(self.np_float_type) + + self.inview_pixels_per_cube = self.np_float_type(pixels_per_cube) + self.inv_scale = self.np_float_type(inv_scale) + self.min_dist = self.np_float_type(min_dist) + + self.center = np.array(center, self.np_float_type) + self.size = self.np_float_type(size) + + self.bisection_iters = bisection_iters + self.enclosed = enclosed + self.simplify_occluded = simplify_occluded + self.visible_relax_iter = visible_relax_iter + self.coarse_count = coarse_count + + register_func(self, dll, "run_coarse", [ + POINTER(self.float_type), self.float_type, + c_int32, POINTER(self.float_type), + self.float_type, self.float_type, self.float_type, + c_int32, c_int32, c_int32, + ], c_int32) + register_func(self, dll, "fine_group", [], c_int32) + register_func(self, dll, "fine_iteration", [POINTER(self.sdf_float_type)], c_int32) + register_func(self, dll, "fine_iteration_output", [POINTER(self.float_type)]) + register_func(self, dll, "vis_filter", [c_bool, c_int32], c_int32) + register_func(self, dll, "final_iteration", [], c_int32) + register_func(self, dll, "final_iteration_occluded", [], c_int32) + register_func(self, dll, "final_iteration2", [POINTER(self.float_type)]) + register_func(self, dll, "final_iteration3", [POINTER(self.sdf_float_type)], c_int32) + register_func(self, dll, "final_iteration3_occluded", [POINTER(self.sdf_float_type)]) + register_func(self, dll, "final_remaining", [POINTER(c_int32)]) + register_func(self, dll, "get_verts_center", [c_int32, POINTER(self.float_type)]) + register_func(self, dll, "update_verts", [c_int32, POINTER(self.sdf_float_type), POINTER(self.sdf_float_type), POINTER(self.float_type)]) + register_func(self, dll, "get_lr_verts", [c_int32, POINTER(self.float_type), POINTER(self.float_type)]) + register_func(self, dll, "finalize_verts", [c_int32, POINTER(self.sdf_float_type), POINTER(self.sdf_float_type), POINTER(self.float_type)]) + register_func(self, dll, "construct_faces", [c_int32, POINTER(self.float_type), POINTER(c_int32)]) + register_func(self, dll, "get_extra_verts_center", [POINTER(self.float_type), POINTER(self.float_type)]) + register_func(self, dll, "update_extra_verts", [ + POINTER(self.sdf_float_type), POINTER(self.sdf_float_type), + POINTER(self.sdf_float_type), POINTER(self.sdf_float_type), + POINTER(self.float_type), POINTER(self.float_type), + ]) + register_func(self, dll, "get_lr_extra_verts", [POINTER(self.float_type), POINTER(self.float_type), POINTER(self.float_type), POINTER(self.float_type)]) + register_func(self, dll, "finalize_extra_verts", [ + POINTER(self.sdf_float_type), POINTER(self.sdf_float_type), POINTER(self.float_type), + POINTER(self.sdf_float_type), POINTER(self.sdf_float_type), POINTER(self.float_type), + ]) + register_func(self, dll, "get_faces", [POINTER(c_int32)]) + register_func(self, dll, "get_in_view_tag", [c_int32, POINTER(c_bool)]) + + + def kernel_caller(self, kernels, XYZ_all): + n_XYZ = len(XYZ_all) + if n_XYZ == 0: return np.zeros((0, len(kernels)), dtype=self.sdf_np_float_type) + step = 10000000 + sdfs = [] + for i in range(0, n_XYZ, step): + XYZ = XYZ_all[i: i+step] + sdfs_i = [] + if self.enclosed: + out_bound = np.zeros(len(XYZ), dtype=bool) + for i in range(3): + out_bound |= XYZ[:, i] <= self.center[i] - self.size / 2 + out_bound |= XYZ[:, i] >= self.center[i] + self.size / 2 + for kernel in kernels: + sdf = kernel(XYZ) + if self.enclosed: sdf[out_bound] = 1 + sdfs_i.append(sdf) + sdfs.append(np.stack(sdfs_i, -1).astype(self.sdf_np_float_type)) + return np.concatenate(sdfs, 0) + + def __call__(self, kernels): + n_elements = len(kernels) + # octree only considering cameras, not sdf + with Timer("coarse step part1"): + n_blocks = self.run_coarse( + self.AF(self.center), self.size, + self.n_cameras, self.AF(self.cameras), + self.inview_pixels_per_cube, + self.inv_scale, self.min_dist, + self.coarse_count, self.memory_limit_mb, n_elements + ) + # start considering sdf + with Timer("coarse step part2"), tqdm(total=n_blocks) as pbar: + while True: + inc = self.fine_group() + if inc == 0: break + pbar.update(inc) + n = self.fine_iteration(POINTER(self.sdf_float_type)()) + while n > 0: + positions = AC(np.zeros((n, 3), dtype=self.np_float_type)) + self.fine_iteration_output(self.AF(positions)) + sdf = AC(self.kernel_caller(kernels, positions).min(axis=-1)) + n = self.fine_iteration(self.sdf_AF(sdf)) + with Timer("filter visible blocks"): + n_vis_block = self.vis_filter(self.simplify_occluded, self.visible_relax_iter) + + with Timer("fine step"), tqdm(total=n_vis_block) as pbar: + nv = np.zeros(1, dtype=np.int32) + while True: + n = self.final_iteration(AsInt(nv)) + if n == 0: break + positions = AC(np.zeros((n, 3), dtype=self.np_float_type)) + self.final_iteration2(self.AF(positions)) + sdf = AC(self.kernel_caller(kernels, positions)) + inc = self.final_iteration3(self.sdf_AF(sdf)) + pbar.update(inc) + n = self.final_iteration_occluded(AsInt(nv)) + if n != 0: + positions = AC(np.zeros((n, 3), dtype=self.np_float_type)) + self.final_iteration2(self.AF(positions)) + sdf = AC(self.kernel_caller(kernels, positions)) + self.final_iteration3_occluded(self.sdf_AF(sdf)) + nv = np.zeros(n_elements, dtype=np.int32) + self.final_remaining(AsInt(nv)) + del positions, sdf + + with Timer("construct mesh"): + meshes = [] + in_view_tags = [] + for e in range(n_elements): + k_e = kernels[e:e+1] + centers = np.zeros((nv[e], 3), dtype=self.np_float_type) + self.get_verts_center(e, self.AF(centers)) + center_sdf = self.kernel_caller(k_e, centers) + cubes = AC(np.zeros((nv[e] * 8, 3), dtype=self.np_float_type)) + self.update_verts(e, POINTER(self.sdf_float_type)(), POINTER(self.sdf_float_type)(), self.AF(cubes)) + for _ in tqdm(range(self.bisection_iters)): + sdf = self.kernel_caller(k_e, cubes) + self.update_verts(e, self.sdf_AF(AC(sdf)), self.sdf_AF(AC(center_sdf)), self.AF(cubes)) + cubes_r = AC(np.zeros((nv[e] * 8, 3), dtype=self.np_float_type)) + self.get_lr_verts(e, self.AF(cubes), self.AF(cubes_r)) + sdf_l = self.kernel_caller(k_e, cubes) + sdf_r = self.kernel_caller(k_e, cubes_r) + del cubes, cubes_r, centers, center_sdf + vertices = np.zeros((nv[e], 3), dtype=self.np_float_type) + self.finalize_verts(e, self.sdf_AF(sdf_l), self.sdf_AF(sdf_r), self.AF(vertices)) + del sdf_l, sdf_r + cnts = np.zeros(3, dtype=np.int32) + self.construct_faces(e, self.AF(vertices), AsInt(cnts)) + nve, nvf, nf = cnts + edge_vertices_c = AC(np.zeros((nve, 3), dtype=self.np_float_type)) + face_vertices_c = AC(np.zeros((nvf, 3), dtype=self.np_float_type)) + self.get_extra_verts_center(self.AF(edge_vertices_c), self.AF(face_vertices_c)) + ecenter_sdf = self.kernel_caller(k_e, edge_vertices_c) + fcenter_sdf = self.kernel_caller(k_e, face_vertices_c) + edge_vertices_lr = AC(np.zeros((nve * 2, 3), dtype=self.np_float_type)) + face_vertices_lr = AC(np.zeros((nvf * 4, 3), dtype=self.np_float_type)) + self.update_extra_verts( + POINTER(self.sdf_float_type)(), POINTER(self.sdf_float_type)(), + POINTER(self.sdf_float_type)(), POINTER(self.sdf_float_type)(), + self.AF(edge_vertices_lr), self.AF(face_vertices_lr), + ) + for _ in range(self.bisection_iters): + e_sdf = self.kernel_caller(k_e, edge_vertices_lr) + f_sdf = self.kernel_caller(k_e, face_vertices_lr) + self.update_extra_verts( + self.sdf_AF(e_sdf), self.sdf_AF(f_sdf), + self.sdf_AF(ecenter_sdf), self.sdf_AF(fcenter_sdf), + self.AF(edge_vertices_lr), self.AF(face_vertices_lr), + ) + del edge_vertices_c, face_vertices_c, ecenter_sdf, fcenter_sdf + edge_vertices_r = AC(np.zeros((nve * 2, 3), dtype=self.np_float_type)) + face_vertices_r = AC(np.zeros((nvf * 4, 3), dtype=self.np_float_type)) + self.get_lr_extra_verts(self.AF(edge_vertices_lr), self.AF(edge_vertices_r), self.AF(face_vertices_lr), self.AF(face_vertices_r)) + esdf_l = self.kernel_caller(k_e, edge_vertices_lr) + esdf_r = self.kernel_caller(k_e, edge_vertices_r) + fsdf_l = self.kernel_caller(k_e, face_vertices_lr) + fsdf_r = self.kernel_caller(k_e, face_vertices_r) + del edge_vertices_lr, edge_vertices_r, face_vertices_lr, face_vertices_r + edge_vertices = np.zeros((nve, 3), dtype=self.np_float_type) + face_vertices = np.zeros((nvf, 3), dtype=self.np_float_type) + self.finalize_extra_verts(self.sdf_AF(esdf_l), self.sdf_AF(esdf_r), self.AF(edge_vertices), self.sdf_AF(fsdf_l), self.sdf_AF(fsdf_r), self.AF(face_vertices)) + del esdf_l, esdf_r, fsdf_l, fsdf_r + faces = AC(np.zeros((nf, 3), dtype=np.int32)) + self.get_faces(AsInt(faces)) + vertices = np.concatenate((vertices, edge_vertices, face_vertices)) + in_view_tag = np.zeros(vertices.shape[0], dtype=bool) + self.get_in_view_tag(e, AsBool(in_view_tag)) + in_view_tags.append(in_view_tag) + meshes.append(trimesh.Trimesh(vertices=vertices, faces=faces, process=False)) + print(f"element {e} has vertices #{meshes[-1].vertices.shape[0]} faces #{meshes[-1].faces.shape[0]}") + return meshes, in_view_tags \ No newline at end of file diff --git a/ocmesher/source/core.cpp b/ocmesher/source/core.cpp new file mode 100644 index 0000000..271c7f6 --- /dev/null +++ b/ocmesher/source/core.cpp @@ -0,0 +1,1095 @@ +// Copyright (c) Princeton University. +// This source code is licensed under the BSD 3-Clause license found in the LICENSE file in the root directory of this source tree. + +// Authors: Zeyu Ma + +#include "core.h" + +namespace coarse { + vector nodes; + std::priority_queue > nodes_heap; + vector nodes_vector; +} + +namespace fine { + int start_node, end_node; + vector > cubes_queue; + vector cubes; + vector cubes_index; + vector vertices; + vector vertices_index; + vector output_vertices; + vector output_vertices_index; +} + +namespace solid { + vector cubes; + set cubes_set; + set visible_set, occluded_set; +} + +namespace final { + queue new_nodes; + vector v; + vector visible_nodes_cube; + vector occluded_nodes_id; + int gl, start_node, end_node, size0; + map vertices; + vector bipolar_edges_s; + vector > bipolar_edges; + vector > bipolar_edges_vindices; + map, int> bipolar_edges_vertices; + vector vertices_cnt; + vector > bipolar_edges_vertices_vector; + vector > bipolar_edges_computed_vertices; + vector nodes; + vector > in_view_tag; + vector searched; +} + +namespace computing { + vector faces; + vector > edge_vertices; + vector edge_vertices_in_view_tag; + vector > face_vertices; + vector face_vertices_in_view_tag; + map, int> face_vertices_map; +} + +extern "C" { + int run_coarse( + T *center, T size, + int n_cams, T *cams, + T pixels_per_cube, + T occ_scale, + T min_dist, + int coarse_count, + int memory_limit_mb, + int n_elements + ) { + using namespace coarse; + params::center = center; + params::size = size; + params::n_cams = n_cams; + params::cams = cams; + params::pixels_per_cube = pixels_per_cube; + params::occ_scale = occ_scale; + params::min_dist = min_dist; + params::coarse_count = coarse_count; + params::memory_limit_mb = memory_limit_mb; + params::n_elements = n_elements; + node root; + mark_leaf_node(root); + memset(root.c.coords, 0, 3 * sizeof(int)); + root.c.L = 0; + nodes.clear(); + nodes.push_back(root); + nodes_heap = std::priority_queue >(); + nodes_heap.push(mp(projected_size(root.c), 0)); + + int t = 0; + while (!nodes_heap.empty() && nodes.size() < coarse_count) { + pair top = nodes_heap.top(); + if (top.first < occ_scale) break; + nodes_heap.pop(); + int i0 = nodes.size(); + expand_octree(nodes, top.second); + for (int i = 0; i < 8; i++) nodes_heap.push(mp(projected_size(nodes[i0 + i].c), i0 + i)); + } + fine::end_node = 0; + solid::cubes.clear(); + solid::cubes_set.clear(); + nodes_vector.clear(); + while (!nodes_heap.empty()) { + pair top = nodes_heap.top(); + mark_grid_node(nodes[top.second], max(0, int_log(top.first / params::occ_scale))); + nodes_vector.push_back(top.second); + nodes_heap.pop(); + } + return nodes_vector.size(); + } + + int fine_group() { + using namespace coarse; + using namespace fine; + cubes_queue.clear(); + cubes_index.clear(); + vertices_index.clear(); + if (end_node == nodes_vector.size()) { + cubes.clear(); + vertices.clear(); + nodes_vector.clear(); + return 0; + } + start_node = end_node; + int cubes_size = 0, vertices_size = 0; + for (;;) { + int top = nodes_vector[end_node++]; + int s = grid_node_level(nodes[top]); + cubes_index.push_back(cubes_size); + vertices_index.push_back(vertices_size); + cubes_size += cubex(1<>20) * (sizeof(bool)+6*sizeof(int)+sizeof(vertex)+(3+params::n_elements)*sizeof(T)) > params::memory_limit_mb) break; + if (end_node == nodes_vector.size()) break; + } + cubes = vector(cubes_size, 0); + vertices = vector(vertices_size, -1); + + for (int i = 0; i < end_node - start_node; i++) { + int ss = 1 << grid_node_level(nodes[nodes_vector[start_node + i]]); + for (int j = 0; j < ss; j++) + for (int k = 0; k < ss; k++) + for (int f = 0; f < 6; f++) { + int coords[3]; + coords[f / 2] = (f & 1) * (ss - 1); + coords[(f/2+1) % 3] = j; + coords[(f/2+2) % 3] = k; + int ci = cube_index(coords[0], coords[1], coords[2], ss), ici = cubes_index[i] + ci; + if (cubes[ici]) continue; + cubes[ici] = true; + cubes_queue.push_back(mp(i, make_int3(coords[0], coords[1], coords[2]))); + } + } + return end_node - start_node; + } + + int fine_iteration(sdfT *sdf) { + using namespace coarse; + using namespace fine; + if (sdf != NULL) { + for (int i = 0; i < output_vertices.size(); i++) { + assert(!std::isnan(sdf[i])); + vertices[output_vertices_index[i]] = sdf[i]>=0? 1: 2; + } + int cqs = cubes_queue.size(); + for (int j = 0; j < cqs; j++) { + int i = cubes_queue[j].first; + int s = grid_node_level(nodes[nodes_vector[start_node + i]]), ss = 1<>1)&1); + int sign1 = vertices[vertices_index[i] + cube_index(vcoords[0], vcoords[1], vcoords[2], ss+1)]; + bool border1 = 0, border2 = 0; + for (int p = 0; p < 3; p++) border1 |= vcoords[p]==0 || vcoords[p]==ss; + vcoords[e/4] = coords[e/4]; + int sign2 = vertices[vertices_index[i] + cube_index(vcoords[0], vcoords[1], vcoords[2], ss+1)]; + for (int p = 0; p < 3; p++) border2 |= vcoords[p]==0 || vcoords[p]==ss; + if (sign1 != sign2) { + flag = true; + if (border1 && border2) { + for (int k = 0; k < 4; k++) { + int inter_coords[3]; + for (int p = 0; p < 3; p++) assign(inter_coords[p], cubei.coords[p], 2<>1)&1)); + cube new_cube = search(&nodes[0], inter_coords, cubei.L + s + 1); + if (!is_boundary(new_cube)) solid::cubes_set.insert(cube_to_key(new_cube)); + } + } + else { + for (int c = 0; c < 4; c++) { + int ccoords[3]; + ccoords[e/4] = vcoords[e/4]; + ccoords[(e/4+1) % 3] = vcoords[(e/4+1) % 3] + (c&1) - 1; + ccoords[(e/4+2) % 3] = vcoords[(e/4+2) % 3] + ((c>>1)&1) - 1; + int ci = cube_index(ccoords[0], ccoords[1], ccoords[2], ss), ici = cubes_index[i] + ci; + if (cubes[ici]) continue; + cubes[ici] = true; + cubes_queue.push_back(mp(i, make_int3(ccoords[0], ccoords[1], ccoords[2]))); + cube c0; + for (int p = 0; p < 3; p++) assign(c0.coords[p], cubei.coords[p], 1<::iterator iter = cubes_set.begin(); iter != cubes_set.end(); iter++) { + cube c; + key_to_cube(c, *iter); + cubes.push_back(c); + } + cubes_set.clear(); + + vector visible(cubes.size(), false); + T factor = 10; + vector canvas; + for (int k = 0; k < n_cams; k++) { + T *current_cam = cams + k * (12 + 9 + 2); + int H = int(current_cam[21] / factor), W = int(current_cam[22] / factor); + if (simplify_occluded) { + canvas = vector(H * W, std::numeric_limits::infinity()); + #pragma omp parallel for + for (int i = 0; i < cubes.size(); i++) { + T image_coords[3]; + projected_coords(cubes[i], k, image_coords, NULL); + if (image_coords[2] >= 0) { + int x = floor(image_coords[0] / factor); + int y = floor(image_coords[1] / factor); + if (x >= 0 && y >= 0 && x < W && y < H) { + #pragma omp critical + { + if (image_coords[2] <= canvas[x * H + y]) { + canvas[x * H + y] = image_coords[2]; + } + } + } + } + } + } + #pragma omp parallel for + for (int i = 0; i < cubes.size(); i++) { + T image_coords[3]; + projected_coords(cubes[i], k, image_coords, NULL); + if (image_coords[2] >= 0) { + int x = floor(image_coords[0] / factor); + int y = floor(image_coords[1] / factor); + if (x >= -relax_iters && y >= -relax_iters && x < W + relax_iters && y < H + relax_iters) { + if (simplify_occluded) { + for (int dx = -relax_iters; dx <= relax_iters; dx++) + for (int dy = -relax_iters; dy <= relax_iters; dy++) { + int nx = x + dx, ny = y + dy; + if (nx >= 0 && ny >= 0 && nx < W && ny < H) { + if (image_coords[2] <= canvas[nx * H + ny]) { + visible[i] = true; + } + } + } + } + else visible[i] = true; + } + } + } + } + visible_set.clear(); + occluded_set.clear(); + set new_visible_set, old_visible_set; + for (int i = 0; i < visible.size(); i++) { + if (visible[i]) visible_set.insert(cube_to_key(cubes[i])); + else occluded_set.insert(cube_to_key(cubes[i])); + } + visible.clear(); + cubes.clear(); + for (int i = 0; i < relax_iters; i++) { + for (set::iterator iter = visible_set.begin(); iter != visible_set.end(); iter++) { + cube c; + key_to_cube(c, *iter); + int coords[3]; + for (int f = 0; f < 6; f++) { + assign(coords[f/3], c.coords[f/3], 2, -1 + 4*(f&1)); + assign(coords[(f/3+1)%3], c.coords[(f/3+1)%3], 2, 1); + assign(coords[(f/3+2)%3], c.coords[(f/3+2)%3], 2, 1); + cube new_cube = search(&coarse::nodes[0], coords, c.L + 1); + if (!is_boundary(new_cube)) { + if (occluded_set.count(cube_to_key(new_cube))) { + occluded_set.erase(cube_to_key(new_cube)); + new_visible_set.insert(cube_to_key(new_cube)); + } + } + } + } + old_visible_set.insert(visible_set.begin(), visible_set.end()); + visible_set = new_visible_set; + new_visible_set.clear(); + } + visible_set.insert(old_visible_set.begin(), old_visible_set.end()); + old_visible_set.clear(); + new_visible_set.clear(); + + node root; + memset(root.nxts, -1, 8 * sizeof(int)); + memset(root.c.coords, 0, 3 * sizeof(int)); + root.c.L = 0; + final::nodes.clear(); + final::nodes.push_back(root); + final::visible_nodes_cube.clear(); + final::occluded_nodes_id.clear(); + for (set::iterator iter = visible_set.begin(); iter != visible_set.end(); iter++) { + cube c; + key_to_cube(c, *iter); + final::visible_nodes_cube.push_back(c); + } + for (set::iterator iter = occluded_set.begin(); iter != occluded_set.end(); iter++) { + cube c; + key_to_cube(c, *iter); + final::occluded_nodes_id.push_back(divide_to_cube(final::nodes, c)); + } + final::bipolar_edges.clear(); + final::bipolar_edges_s.clear(); + final::bipolar_edges_vertices.clear(); + final::in_view_tag.clear(); + final::bipolar_edges_vindices.clear(); + for (int i = 0; i < params::n_elements; i++) { + final::bipolar_edges.push_back(vector()); + final::bipolar_edges_s.push_back(0); + final::bipolar_edges_vindices.push_back(vector()); + final::in_view_tag.push_back(vector()); + } + final::vertices_cnt = vector(5, 0); + final::bipolar_edges.push_back(vector()); + final::start_node = 0; + final::gl = 0; + for (;;) { + int total_nodes = 0; + for (int i = 0; i < final::visible_nodes_cube.size(); i++) { + int is = max(0, int_log(projected_size(final::visible_nodes_cube[i])) - final::gl); + total_nodes += max(0, cubex(1<>20) * sizeof(node) < params::memory_limit_mb * 0.6) break; + final::gl++; + } + return final::visible_nodes_cube.size(); + } + + int final_iteration() { + using namespace final; + assert((nodes.size() >> 20) * sizeof(node) < params::memory_limit_mb * 0.8); + if (start_node == visible_nodes_cube.size()) { + coarse::nodes.clear(); + solid::visible_set.clear(); + visible_nodes_cube.clear(); + return 0; + } + int total_nodes = 0, total_verts = 0; + for (end_node = start_node; end_node < visible_nodes_cube.size(); end_node++) { + int is = int_log(projected_size(visible_nodes_cube[end_node])); + total_nodes += (1 << (3*(is-gl))) * 8 / 7; + total_verts += cubex((1<>20) * sizeof(node) + (total_verts>>20) * (sizeof(vertex)+sizeof(T*)+sizeof(key_cube)+sizeof(int)+2*sizeof(key_cube*)+(3+params::n_elements)*sizeof(T)) > params::memory_limit_mb) break; + } + assert(end_node != start_node); + + vertices.clear(); + size0 = nodes.size(); + assert(new_nodes.empty()); + for (int i = start_node; i < end_node; i++) { + new_nodes.push(divide_to_cube(nodes, visible_nodes_cube[i])); + while (!new_nodes.empty()) { + int ind = new_nodes.front(); + new_nodes.pop(); + T s = projected_size(nodes[ind].c); + if (s > 1< 0); + return vertices.size(); + } + + int final_iteration_occluded() { + using namespace final; + for (int i = 0; i < occluded_nodes_id.size(); i++) { + v.resize(8); + enumerate_vertices(&v[0], nodes[occluded_nodes_id[i]]); + for (int j = 0; j < 8; j++) if (!vertices.count(cube_to_key(v[j]))) vertices[cube_to_key(v[j])] = vertices.size(); + } + return vertices.size(); + } + + void final_iteration2(T *xyz) { + using namespace params; + using namespace final; + for (map::iterator iter = vertices.begin(); iter != vertices.end(); iter++) { + vertex v; + key_to_cube(v, iter->first); + compute_coords(xyz + iter->second*3, v.coords, v.L); + } + } + + int final_iteration3(sdfT *sdf) { + using namespace final; + using namespace solid; + int start = visible_nodes_cube.size() - start_node; + for (int i = 0; i < vertices.size() * params::n_elements; i++) assert(!std::isnan(sdf[i])); + for (int i = size0; i < nodes.size(); i++) + if (leaf_node(nodes[i])) find_edges(nodes[i], vertices, sdf, bipolar_edges); + vertices.clear(); + for (int i = 0; i < params::n_elements + 1; i++) { + vector &bei = bipolar_edges[i]; + int s = i==params::n_elements?0:bipolar_edges_s[i]; + std::sort(bei.begin() + s, bei.end()); + bei.erase(std::unique(bei.begin() + s, bei.end()), bei.end()); + int e = bei.size(); + searched.resize((e-s)*4); + #pragma omp parallel for + for (int j = s; j < e; j++) { + cube e0; + key_to_cube(e0, bei[j].second); + int dir = bei[j].first; + int coords[3], L = e0.L + 1; + int ks[4] = {0, 1, 3, 2}; + bool flag = true, flag2 = true; + for (int ik = 0; ik < 4; ik++) { + int k = dir > 0? ks[ik]: ks[3 - ik]; + int c0 = abs(dir) - 1, c1 = (c0 + 1) % 3, c2 = (c1 + 1) % 3; + assign(coords[c0], e0.coords[c0], 2, 1); + assign(coords[c1], e0.coords[c1], 2, -1 + 2*(k&1)); + assign(coords[c2], e0.coords[c2], 2, -1 + 2*((k>>1)&1)); + cube &c = searched[(j-s)*4 + ik]; + c = search(&nodes[0], coords, L); + if (is_boundary(c)) { + flag = false; + break; + } + else if (is_exterior(c)) { + flag2 = false; + } + } + if (!flag) { + for (int k = 0; k < 4; k++) + mark_boundary(searched[(j-s)*4 + k]); + } + else if (!flag2) { + for (int k = 0; k < 4; k++) + mark_exterior(searched[(j-s)*4 + k]); + } + } + if (i == params::n_elements) continue; + int j1 = s, j2 = e - 1; + while (j1 < j2) { + while (j1 < e && is_regular(searched[4*(j1-s)])) j1++; + while (j2 >= s && !is_regular(searched[4*(j2-s)])) j2--; + if (j1 < j2) { + std::swap(bei[j1], bei[j2]); + for (int k = 0; k < 4; k++) + std::swap(searched[4*(j1-s) + k], searched[4*(j2-s) + k]); + j1++; + j2--; + } + } + vector &ivt = in_view_tag[i]; + for (int j = s; j < j1; j++) { + for (int k = 0; k < 4; k++) { + cube c = searched[(j-s)*4 + k]; + int vid; + key_cube key = cube_to_key(c); + if (bipolar_edges_vertices.count(mp(i, key))) vid = bipolar_edges_vertices[mp(i, key)]; + else { + vid = bipolar_edges_vertices[mp(i, key)] = vertices_cnt[i]++; + ivt.push_back(occluded_set.count(key) == 0); + } + bipolar_edges_vindices[i].push_back(vid); + } + } + bipolar_edges_s[i] = j1; + j2 = e - 1; + while (j1 < j2) { + while (j1 < e && is_exterior(searched[4*(j1-s)])) j1++; + while (j2 >= bipolar_edges_s[i] && is_boundary(searched[4*(j2-s)])) j2--; + if (j1 < j2) { + std::swap(bei[j1], bei[j2]); + for (int k = 0; k < 4; k++) + std::swap(searched[4*(j1-s) + k], searched[4*(j2-s) + k]); + j1++; + j2--; + } + } + bei.erase(bei.begin() + j1, bei.end()); + } + + for (int i = 0; i < size0; i++) + for (int j = 0; j < 8; j++) + if (nodes[i].nxts[j] >= size0) nodes[i].nxts[j] = -1; + nodes.erase(nodes.begin() + size0, nodes.end()); + + for (int i = start_node; i < end_node; i++) { + assert(new_nodes.empty()); + int nodes_id = divide_to_cube(nodes, visible_nodes_cube[i]); + memset(nodes[nodes_id].nxts, -1, 3 * sizeof(int)); + new_nodes.push(nodes_id); + while (!new_nodes.empty()) { + int ind = new_nodes.front(); + new_nodes.pop(); + T s = projected_size(nodes[ind].c); + if (s > 1< &ben = bipolar_edges[params::n_elements]; + for (int i = 0; i < ben.size(); i++) { + int dir = ben[i].first; + cube e0; + key_to_cube(e0, ben[i].second); + int coords[3], L = e0.L + 1; + for (int k = 0; k < 4; k++) { + int c0 = abs(dir) - 1, c1 = (c0 + 1) % 3, c2 = (c1 + 1) % 3; + assign(coords[c0], e0.coords[c0], 2, 1); + assign(coords[c1], e0.coords[c1], 2, -1 + 2*(k&1)); + assign(coords[c2], e0.coords[c2], 2, -1 + 2*((k>>1)&1)); + cube c = search(&coarse::nodes[0], coords, L); + if (!is_boundary(c)) { + key_cube key = cube_to_key(c); + if (occluded_set.count(key) || visible_set.count(key)) continue; + visible_set.insert(key); + visible_nodes_cube.push_back(c); + } + } + } + ben.clear(); + + start_node = end_node; + return start - (visible_nodes_cube.size() - start_node); + } + + void final_iteration3_occluded(sdfT *sdf) { + using namespace final; + using namespace solid; + for (int i = 0; i < vertices.size() * params::n_elements; i++) assert(!std::isnan(sdf[i])); + for (int i = 0; i < occluded_nodes_id.size(); i++) + find_edges(nodes[occluded_nodes_id[i]], vertices, sdf, bipolar_edges); + vertices.clear(); + occluded_nodes_id.clear(); + } + + void final_remaining(int *nv) { + using namespace final; + bipolar_edges_computed_vertices.clear(); + bipolar_edges_vertices_vector.clear(); + for (int i = 0; i < params::n_elements; i++) { + vector &bei = bipolar_edges[i]; + int s = bipolar_edges_s[i]; + std::sort(bei.begin() + s, bei.end()); + bei.erase(std::unique(bei.begin() + s, bei.end()), bei.end()); + int e = bei.size(); + searched.resize((e-s)*4); + vector tags((e-s)*4); + #pragma omp parallel for + for (int j = s; j < e; j++) { + cube e0; + key_to_cube(e0, bei[j].second); + int dir = bei[j].first; + int coords[3], L = e0.L + 1; + int ks[4] = {0, 1, 3, 2}; + bool flag = false; + for (int ik = 0; ik < 4; ik++) { + int k = dir > 0? ks[ik]: ks[3 - ik]; + int c0 = abs(dir) - 1, c1 = (c0 + 1) % 3, c2 = (c1 + 1) % 3; + assign(coords[c0], e0.coords[c0], 2, 1); + assign(coords[c1], e0.coords[c1], 2, -1 + 2*(k&1)); + assign(coords[c2], e0.coords[c2], 2, -1 + 2*((k>>1)&1)); + cube &c = searched[(j-s)*4 + ik]; + c = search(&nodes[0], coords, L, 1); + if (is_boundary(c)) { + flag = true; + break; + } + else { + tags[(j-s)*4 + ik] = solid::occluded_set.count(cube_to_key(c)) == 0 && !is_exterior(search(&nodes[0], coords, L)); + } + } + if (flag) { + for (int k = 0; k < 4; k++) + mark_boundary(searched[(j-s)*4 + k]); + } + } + int j1 = s, j2 = e - 1; + while (j1 < j2) { + while (j1 < e && is_regular(searched[4*(j1-s)])) j1++; + while (j2 >= s && !is_regular(searched[4*(j2-s)])) j2--; + if (j1 < j2) { + std::swap(bei[j1], bei[j2]); + for (int k = 0; k < 4; k++) { + std::swap(searched[4*(j1-s) + k], searched[4*(j2-s) + k]); + std::swap(tags[4*(j1-s) + k], tags[4*(j2-s) + k]); + } + j1++; + j2--; + } + } + vector &ivt = in_view_tag[i]; + for (int j = s; j < j1; j++) { + for (int k = 0; k < 4; k++) { + cube c = searched[(j-s)*4 + k]; + int vid; + key_cube key = cube_to_key(c); + if (bipolar_edges_vertices.count(mp(i, key))) vid = bipolar_edges_vertices[mp(i, key)]; + else { + vid = bipolar_edges_vertices[mp(i, key)] = vertices_cnt[i]++; + ivt.push_back(tags[(j-s)*4 + k]); + } + bipolar_edges_vindices[i].push_back(vid); + } + } + bei.erase(bei.begin() + j1, bei.end()); + bipolar_edges_computed_vertices.push_back(vector(vertices_cnt[i])); + bipolar_edges_vertices_vector.push_back(vector(vertices_cnt[i])); + nv[i] = vertices_cnt[i]; + } + for (map, int>::iterator iter = bipolar_edges_vertices.begin(); iter != bipolar_edges_vertices.end(); iter++) { + int i = iter->first.first; + vector &becvi = bipolar_edges_computed_vertices[i]; + vector &bevvi = bipolar_edges_vertices_vector[i]; + cube c; + key_to_cube(c, iter->first.second); + compute_center(becvi[iter->second].c, c); + becvi[iter->second].l = 0; + becvi[iter->second].r = 0.5 * params::size / (1<second] = iter->first.second; + } + bipolar_edges_vertices.clear(); + nodes.clear(); + solid::occluded_set.clear(); + } + + void get_verts_center(int e, T *positions) { + using namespace final; + vector &becv = bipolar_edges_computed_vertices[e]; + for (int i = 0; i < becv.size(); i++) { + memcpy(positions + 3 * i, becv[i].c, 3 * sizeof(T)); + } + } + + void get_extra_verts_center(T *epositions, T *fpositions) { + using namespace computing; + for (int i = 0; i < edge_vertices.size(); i++) { + memcpy(epositions + 3 * i, edge_vertices[i].second.c, 3 * sizeof(T)); + } + for (int i = 0; i < face_vertices.size(); i++) { + memcpy(fpositions + 3 * i, face_vertices[i].second.c, 3 * sizeof(T)); + } + } + + void update_verts(int e, sdfT *sdf, sdfT *center_sdf, T *positions) { + using namespace final; + vector &becv = bipolar_edges_computed_vertices[e]; + #pragma omp parallel for + for (int i = 0; i < becv.size(); i++) { + if (sdf != NULL) { + T mid = (becv[i].l + becv[i].r) / 2; + bool bipolar = 0; + for (int j = 0; j < 8; j++) + if ((sdf[i * 8 + j] >= 0) != (center_sdf[i] >= 0)) { + bipolar = 1; + break; + } + if (bipolar) becv[i].r = mid; + else becv[i].l = mid; + } + T mid = (becv[i].l + becv[i].r) / 2; + for (int j = 0; j < 8; j++) { + int cid = i * 8 + j; + for (int k = 0; k < 3; k++) + positions[cid*3 + k] = becv[i].c[k] + (((j>>k)&1)*2-1) * mid; + } + } + } + + void update_extra_verts(sdfT *esdf, sdfT *fsdf, sdfT *ecenter_sdf, sdfT *fcenter_sdf, T *epositions, T*fpositions) { + using namespace computing; + int edge_n = 2; + #pragma omp parallel for + for (int i = 0; i < edge_vertices.size(); i++) { + if (esdf != NULL) { + T mid = (edge_vertices[i].second.l + edge_vertices[i].second.r) / 2; + bool bipolar = 0; + for (int j = 0; j < edge_n; j++) + if ((esdf[i * edge_n + j] >= 0) != (ecenter_sdf[i] >= 0)) { + bipolar = 1; + break; + } + if (bipolar) edge_vertices[i].second.r = mid; + else edge_vertices[i].second.l = mid; + } + T mid = (edge_vertices[i].second.l + edge_vertices[i].second.r) / 2; + for (int j = 0; j < edge_n; j++) { + int cid = i * edge_n + j; + for (int k = 0; k < 3; k++) { + epositions[cid*3 + k] = edge_vertices[i].second.c[k]; + if (k == edge_vertices[i].first) epositions[cid*3 + k] += (j*2-1) * mid; + } + } + } + int face_n = 4; + #pragma omp parallel for + for (int i = 0; i < face_vertices.size(); i++) { + if (fsdf != NULL) { + T mid = (face_vertices[i].second.l + face_vertices[i].second.r) / 2; + bool bipolar = 0; + for (int j = 0; j < face_n; j++) + if ((fsdf[i * face_n + j] >= 0) != (fcenter_sdf[i] >= 0)) { + bipolar = 1; + break; + } + if (bipolar) face_vertices[i].second.r = mid; + else face_vertices[i].second.l = mid; + } + T mid = (face_vertices[i].second.l + face_vertices[i].second.r) / 2; + for (int j = 0; j < face_n; j++) { + int cid = i * face_n + j; + for (int k = 0; k < 3; k++) { + fpositions[cid*3 + k] = face_vertices[i].second.c[k]; + if (k == (face_vertices[i].first+1)%3) fpositions[cid*3 + k] += (first_digit(j)*2-1) * mid; + else if (k == (face_vertices[i].first+2)%3) fpositions[cid*3 + k] += (second_digit(j)*2-1) * mid; + } + } + } + } + + void get_lr_verts(int e, T *cube_l, T *cube_r) { + using namespace final; + vector &becv = bipolar_edges_computed_vertices[e]; + for (int i = 0; i < becv.size(); i++) { + for (int j = 0; j < 8; j++) { + int cid = i * 8 + j; + for (int k = 0; k < 3; k++) { + cube_l[cid*3 + k] = becv[i].c[k] + (((j>>k)&1)*2-1) * becv[i].l; + cube_r[cid*3 + k] = becv[i].c[k] + (((j>>k)&1)*2-1) * becv[i].r; + } + } + } + } + + void get_lr_extra_verts(T *epos_l, T *epos_r, T *fpos_l, T *fpos_r) { + using namespace computing; + for (int i = 0; i < edge_vertices.size(); i++) { + for (int j = 0; j < 2; j++) { + int cid = i * 2 + j; + for (int k = 0; k < 3; k++) { + epos_l[cid*3 + k] = edge_vertices[i].second.c[k]; + if (k == edge_vertices[i].first) epos_l[cid*3 + k] += (j*2-1) * edge_vertices[i].second.l; + epos_r[cid*3 + k] = edge_vertices[i].second.c[k]; + if (k == edge_vertices[i].first) epos_r[cid*3 + k] += (j*2-1) * edge_vertices[i].second.r; + } + } + } + for (int i = 0; i < face_vertices.size(); i++) { + for (int j = 0; j < 4; j++) { + int cid = i * 4 + j; + for (int k = 0; k < 3; k++) { + fpos_l[cid*3 + k] = face_vertices[i].second.c[k]; + if (k == (face_vertices[i].first+1)%3) fpos_l[cid*3 + k] += (first_digit(j)*2-1) * face_vertices[i].second.l; + else if (k == (face_vertices[i].first+2)%3) fpos_l[cid*3 + k] += (second_digit(j)*2-1) * face_vertices[i].second.l; + fpos_r[cid*3 + k] = face_vertices[i].second.c[k]; + if (k == (face_vertices[i].first+1)%3) fpos_r[cid*3 + k] += (first_digit(j)*2-1) * face_vertices[i].second.r; + else if (k == (face_vertices[i].first+2)%3) fpos_r[cid*3 + k] += (second_digit(j)*2-1) * face_vertices[i].second.r; + } + } + } + } + + // todo consider more than corners when a vetex cube has complex side face + void finalize_verts(int e, sdfT *sdf_l, sdfT *sdf_r, T *verts) { + using namespace final; + vector &becv = bipolar_edges_computed_vertices[e]; + for (int i = 0; i < becv.size(); i++) { + T v[3]={0}; + int w=0; + for (int j = 0; j < 8; j++) + if ((sdf_l[i * 8 + j] >= 0) != (sdf_r[i * 8 + j] >= 0)) { + w++; + for (int k = 0; k < 3; k++) { + v[k] += becv[i].c[k] + (((j>>k)&1)*2-1) * (becv[i].l + becv[i].r) / 2; + } + } + if (w == 0) { + for (int k = 0; k < 3; k++) verts[i * 3 + k] = becv[i].c[k]; + } + else { + for (int k = 0; k < 3; k++) verts[i * 3 + k] = v[k] / w; + } + } + bipolar_edges_computed_vertices[e].clear(); + } + + void finalize_extra_verts(sdfT *esdf_l, sdfT *esdf_r, T *everts, sdfT *fsdf_l, sdfT *fsdf_r, T *fverts) { + using namespace computing; + for (int i = 0; i < edge_vertices.size(); i++) { + T v[3]={0}; + int w=0; + for (int j = 0; j < 2; j++) + if ((esdf_l[i * 2 + j] >= 0) != (esdf_r[i * 2 + j] >= 0)) { + w++; + T mid = (edge_vertices[i].second.l + edge_vertices[i].second.r) / 2; + for (int k = 0; k < 3; k++) { + v[k] += edge_vertices[i].second.c[k]; + if (k == edge_vertices[i].first) v[k] += (j*2-1) * mid; + } + } + assert(w != 0); + for (int k = 0; k < 3; k++) everts[i * 3 + k] = v[k] / w; + } + edge_vertices.clear(); + for (int i = 0; i < face_vertices.size(); i++) { + T v[3]={0}; + int w=0; + for (int j = 0; j < 4; j++) + if ((fsdf_l[i * 4 + j] >= 0) != (fsdf_r[i * 4 + j] >= 0)) { + w++; + T mid = (face_vertices[i].second.l + face_vertices[i].second.r) / 2; + for (int k = 0; k < 3; k++) { + v[k] += face_vertices[i].second.c[k]; + if (k == (face_vertices[i].first+1)%3) v[k] += (first_digit(j)*2-1) * mid; + else if (k == (face_vertices[i].first+2)%3) v[k] += (second_digit(j)*2-1) * mid; + } + } + if (w == 0) { + for (int k = 0; k < 3; k++) fverts[i * 3 + k] = face_vertices[i].second.c[k]; + } + else { + for (int k = 0; k < 3; k++) fverts[i * 3 + k] = v[k] / w; + } + } + face_vertices.clear(); + } + + void get_in_view_tag(int e, bool *output) { + using namespace final; + using namespace computing; + int cnt = 0; + vector &ivt = in_view_tag[e]; + for (int i = 0; i < ivt.size(); i++) { + output[cnt++] = ivt[i]; + } + ivt.clear(); + for (int i = 0; i < edge_vertices_in_view_tag.size(); i++) { + output[cnt++] = edge_vertices_in_view_tag[i]; + } + edge_vertices_in_view_tag.clear(); + for (int i = 0; i < face_vertices_in_view_tag.size(); i++) { + output[cnt++] = face_vertices_in_view_tag[i]; + } + face_vertices_in_view_tag.clear(); + } + + void construct_faces(int e, T *final_vertices, int *cnt) { + using namespace final; + using namespace computing; + vector &edge_e = bipolar_edges[e]; + vector &vertices_ids = bipolar_edges_vindices[e]; + vector &unique_vertices = bipolar_edges_vertices_vector[e]; + faces.clear(); + edge_vertices.clear(); + edge_vertices_in_view_tag.clear(); + face_vertices.clear(); + face_vertices_in_view_tag.clear(); + face_vertices_map.clear(); + int nv = unique_vertices.size(); + for (int i = 0; i < edge_e.size(); i++) { + T computed_edge[6]; + cube c; + key_to_cube(c, edge_e[i].second); + compute_coords(computed_edge, c.coords, c.L); + int dir = abs(edge_e[i].first) - 1; + c.coords[dir]++; + compute_coords(computed_edge + 3, c.coords, c.L); + T computed_faces[12 * 4]; + int computed_faces_L[4], computed_faces_dir[4]; + bool intersect[4]={0}; + bool condition1 = 1; + for (int j = 0; j < 4; j++) { + int v1 = vertices_ids[i * 4 + j], v2 = vertices_ids[i * 4 + (j+1) % 4]; + if (v1 == v2) continue; + cube c1, c2; + key_to_cube(c1, unique_vertices[v1]); + key_to_cube(c2, unique_vertices[v2]); + if (c1.L < c2.L) { + cube t = c1; + c1 = c2; + c2 = t; + } + bool found = 0; + for (int ax = 0; ax < 3; ax++) { + for (int dir = 0; dir <= 1; dir++) { + if (c1.coords[ax] + dir == (c2.coords[ax] + 1 - dir) << (c1.L - c2.L)) { + for (int d1 = 0; d1 <= 1; d1++) + for (int d2 = 0; d2 <= 1; d2++) { + cube c1_ = c1; + c1_.coords[ax] += dir; + c1_.coords[(ax+1)%3] += d1; + c1_.coords[(ax+2)%3] += d2; + compute_coords(computed_faces + 12 * j + 3 * (d1 * 2 + d2), c1_.coords, c1_.L); + } + computed_faces_L[j] = c1.L; + computed_faces_dir[j] = ax; + found = 1; + break; + } + } + if (found) break; + } + assert(found); + if (c1.L == c2.L) { + intersect[j] = 1; + continue; + } + // tri_seg_intersect should be strict + bool i1 = tri_seg_intersect(computed_faces + 12 * j + 6, computed_faces + 12 * j, computed_faces + 12 * j + 3, final_vertices + v1 * 3, final_vertices + v2 * 3); + bool i2 = tri_seg_intersect(computed_faces + 12 * j + 3, computed_faces + 12 * j + 9, computed_faces + 12 * j + 6, final_vertices + v1 * 3, final_vertices + v2 * 3); + intersect[j] = i1 || i2; + condition1 &= intersect[j]; + } + if (!condition1) { + computed_vertex v; + for (int j = 0; j < 3; j++) + v.c[j] = (computed_edge[j] + computed_edge[j + 3]) / 2; + v.l = 0; + v.r = 0.5 * params::size / (1< +#include +#include +#include +#include +#include +#include +#include +#include +#define mp std::make_pair +#define pair std::pair +#define max std::max +#define min std::min +#define INT_MAX 2147483647 +#define cubex(x) (x) * (x) * (x) +#define cube_index(x, y, z, s) (x)*(s)*(s)+(y)*(s)+(z) +#define xpp(p) p.first +#define ypp(p) p.second.first +#define zpp(p) p.second.second +#define vector std::vector +#define set std::set +#define map std::map +#define queue std::queue +#define cube_to_key(cube) mp(cube.coords[0], mp(cube.coords[1], mp(cube.coords[2], cube.L))) +#define key_to_cube(c, key) {c.coords[0] = (key).first; c.coords[1] = (key).second.first; c.coords[2] = (key).second.second.first; c.L = (key).second.second.second; } +#define assign(x, y, a, b) {assert((y) < (INT_MAX-(max(0,b)))/(a)); x=(y)*(a)+(b);} +#define leaf_node(node) (node.nxts[0] <= -2) +#define mark_leaf_node(node) node.nxts[0] = -2 +#define mark_grid_node(node, gl) node.nxts[0] = -2-gl +#define grid_node_level(node) (-node.nxts[0] - 2) +#define is_regular(cube) (cube.L >= 0) +#define is_boundary(cube) (cube.L == -1) +#define is_exterior(cube) (cube.L == -2) +#define mark_boundary(cube) cube.L = -1 +#define mark_exterior(cube) cube.L = -2 +#define add_faces(faces, v1, v2, v3) {assert(v1 != v2 && v2 != v3 && v3 != v1); faces.push_back(v1); faces.push_back(v2); faces.push_back(v3);} +#define first_digit(j) (j&1) +#define second_digit(j) ((j>>1)&1) +#define make_int3(x,y,z) mp(x, mp(y, z)) +typedef double T; +typedef float sdfT; + +struct computed_vertex { + T c[3], l, r; +}; + +struct cube { + int coords[3], L; +}; + +typedef cube vertex; +typedef pair > int3; +typedef pair > > key_cube; +typedef pair key_edge; + + +struct node { + cube c; + int nxts[8]; +}; + +int int_log(T x) { + return int(max(T(0), (T)ceil(log2(x)))); +} + +namespace params { + int n_cams, memory_limit_mb, coarse_count, n_elements; + T *center, *cams; + T size, pixels_per_cube, occ_scale, min_dist; +} + +void enumerate_vertices(vertex *v, node n) { + assert(leaf_node(n)); + int s=grid_node_level(n), ss = 1<>= 1; + v[vid].L--; + } + } +} + +void compute_coords(T *coords, int *icoords, int L) { + for (int j = 0; j < 3; j++) + coords[j] = params::center[j] - params::size / 2 + params::size * icoords[j] / (1 << L); +} + +void compute_center(T *coords, cube v) { + for (int j = 0; j < 3; j++) + coords[j] = params::center[j] - params::size / 2 + params::size * (v.coords[j] + 0.5) / (1 << v.L); +} + +void projected_coords(cube c, int k, T *icoords, T *r) { + using namespace params; + T Pw[3], Pc[3]; + T *current_cam = cams + k * (12 + 9 + 2); + compute_center(Pw, c); + for (int i = 0; i < 3; i++) { + Pc[i] = current_cam[i * 4 + 3]; + for (int j = 0; j < 3; j++) { + Pc[i] += Pw[j] * current_cam[i * 4 + j]; + } + } + if (r != NULL) { + *r = sqrt(Pc[0] * Pc[0] + Pc[1] * Pc[1] + Pc[2] * Pc[2]); + *r = max(*r, min_dist); + } + if (icoords != NULL) { + for (int i = 0; i < 3; i++) { + icoords[i] = 0; + for (int j = 0; j < 3; j++) { + icoords[i] += Pc[j] * current_cam[12 + i * 3 + j]; + } + } + icoords[0] /= icoords[2]; + icoords[1] /= icoords[2]; + } +} + +T projected_size(cube c, int k) { + using namespace params; + T *current_cam = cams + k * (12 + 9 + 2); + T r; + projected_coords(c, k, NULL, &r); + T W = current_cam[22]; + T pix_ang = atan(W / 2 / current_cam[12]) * 2 / W; + T ang = pix_ang * pixels_per_cube; + return size / (1 << c.L) / r / ang; +} + +T projected_size(cube c) { + T max_size = 0; + for (int k = 0; k < params::n_cams; k++) { + T size_k = projected_size(c, k); + max_size = max(max_size, size_k); + } + return max_size; +} + +void expand_octree(vector &nodes, int index) { + assert(leaf_node(nodes[index])); + int base_size = nodes.size(); + for (int i = 0; i < 8; i++) { + node *current = &nodes[index]; + current->nxts[i] = nodes.size(); + node node0; + mark_leaf_node(node0); + node0.c.L = current->c.L + 1; + for (int k = 0; k < 3; k++) { + int offset = (i>>k) & 1; + assign(node0.c.coords[k], current->c.coords[k], 2, offset); + } + nodes.push_back(node0); + } +} + +void partial_expand_octree(vector &nodes, int index, int instance) { + int base_size = nodes.size(); + assert(!leaf_node(nodes[index])); + for (int i = 0; i < 8; i++) { + node *current = &nodes[index]; + if ((instance>>i) & 1) { + if (current->nxts[i] == -1) { + current->nxts[i] = nodes.size(); + node node0; + memset(node0.nxts, -1, 8 * sizeof(int)); + node0.c.L = current->c.L + 1; + for (int k = 0; k < 3; k++) { + int offset = (i>>k) & 1; + assign(node0.c.coords[k], current->c.coords[k], 2, offset); + } + nodes.push_back(node0); + } + } + } +} + +cube search(node *nodes, int *coords, int L, bool include_exterior=0) { + assert(L > 0); + cube res; + for (int i = 0; i < 3; i++) + if (!(coords[i] > 0 && coords[i] < 1<> (L - current.c.L - 1))) nxt += (1 << p); + current_id = current.nxts[nxt]; + if (current_id == -1) { + if (!include_exterior) { + mark_exterior(res); + } + else { + for (int k = 0; k < 3; k++) { + int offset = (nxt>>k) & 1; + assign(res.coords[k], current.c.coords[k], 2, offset); + } + res.L = current.c.L + 1; + } + return res; + } + current = nodes[current_id]; + } + if (!leaf_node(current)) { + mark_boundary(res); + return res; + } + int gl = grid_node_level(current); + if (L <= current.c.L + gl) { + mark_boundary(res); + return res; + } + for (int i = 0; i < 3; i++) + if (coords[i] & ((1<<(L-current.c.L-gl))-1) == 0) { + mark_boundary(res); + return res; + } + res.L = current.c.L + gl; + for (int p = 0; p < 3; p++) res.coords[p] = coords[p] >> (L - res.L); + return res; +} + +int divide_to_cube(vector &nodes, cube c) { + node current = nodes[0]; + int current_id = 0; + for (;;) { + assert(!leaf_node(current)); + bool flag = 1; + for (int p = 0; p < 3; p++) flag &= current.c.coords[p] == c.coords[p]; + if (flag && current.c.L == c.L) { + mark_leaf_node(nodes[current_id]); + return current_id; + } + int nxt = 0; + for (int p = 0; p < 3; p++) + if (2 * current.c.coords[p] < c.coords[p] >> (c.L - current.c.L - 1)) nxt += (1 << p); + partial_expand_octree(nodes, current_id, 1< &vertices, sdfT *sdf, vector > &bipolar_edges) { + int s = grid_node_level(n), ss=1<::infinity(), sdf2_min=std::numeric_limits::infinity(); + for (int e = 0; e < params::n_elements; e++) { + sdf1_min = min(sdf1[e], sdf1_min); + sdf2_min = min(sdf2[e], sdf2_min); + if ((sdf1[e] >= 0) != (sdf2[e] >= 0)) { + cube c0; + for (int p = 0; p < 3; p++) + assign(c0.coords[p], n.c.coords[p], ss, coords[p]); + c0.L = n.c.L + s; + int dir = edir + 1; + if (sdf1[e] < 0) dir *= -1; + bipolar_edges[e].push_back(mp(dir, cube_to_key(c0))); + } + } + if ((sdf1_min >= 0) != (sdf2_min >= 0)) { + cube c0; + for (int p = 0; p < 3; p++) + assign(c0.coords[p], n.c.coords[p], ss, coords[p]); + c0.L = n.c.L + s; + int dir = edir + 1; + if (sdf1_min < 0) dir *= -1; + bipolar_edges[params::n_elements].push_back(mp(dir, cube_to_key(c0))); + } + } +} + +int compute_boundary(cube c, cube bound) { + int b[3][2]; + for (int i = 0; i < 3; i++) + for (int p = 0; p < 2; p++) { + b[i][p] = c.coords[i] + p == (bound.coords[i]+p) << (c.L-bound.L); + } + int instance=0; + for (int i = 0; i < 8; i++) { + for (int j = 0; j < 3; j++) + if (b[j][(i>>j)&1]) { + instance |= 1 << i; + break; + } + } + return instance; +} + +inline T det(T matrix[3][3]) { + T a = matrix[0][0]; + T b = matrix[0][1]; + T c = matrix[0][2]; + T d = matrix[1][0]; + T e = matrix[1][1]; + T f = matrix[1][2]; + T g = matrix[2][0]; + T h = matrix[2][1]; + T i = matrix[2][2]; + return a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g); +} + +bool tri_seg_intersect(T *t1, T *t2, T *t3, T *s1, T *s2) { + // note (t1,t3) of the tri is allowed to intersect + T m[3][3]; + for (int i = 0; i < 3; i++) { + m[0][i] = t1[i] - s1[i]; + m[1][i] = t2[i] - s1[i]; + m[2][i] = t3[i] - s1[i]; + } + T det1 = det(m); + for (int i = 0; i < 3; i++) { + m[0][i] = t1[i] - s2[i]; + m[1][i] = t2[i] - s2[i]; + m[2][i] = t3[i] - s2[i]; + } + T det2 = det(m); + if (!(det1 > 0 && det2 < 0 || det1 < 0 && det2 > 0)) return 0; + + for (int i = 0; i < 3; i++) { + m[0][i] = t1[i] - s1[i]; + m[1][i] = t2[i] - s1[i]; + m[2][i] = s2[i] - s1[i]; + } + det1 = det(m); + for (int i = 0; i < 3; i++) { + m[0][i] = t2[i] - s1[i]; + m[1][i] = t3[i] - s1[i]; + } + det2 = det(m); + if (!(det1 > 0 && det2 > 0 || det1 < 0 && det2 < 0)) return 0; + for (int i = 0; i < 3; i++) { + m[0][i] = t3[i] - s1[i]; + m[1][i] = t1[i] - s1[i]; + } + det1 = det(m); + if (!(det1 >= 0 && det2 > 0 || det1 <= 0 && det2 < 0)) return 0; + return 1; +} diff --git a/ocmesher/utils/__init__.py b/ocmesher/utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/ocmesher/utils/interface.py b/ocmesher/utils/interface.py new file mode 100644 index 0000000..a6c0e03 --- /dev/null +++ b/ocmesher/utils/interface.py @@ -0,0 +1,31 @@ +# Copyright (c) Princeton University. +# This source code is licensed under the BSD 3-Clause license found in the LICENSE file in the root directory of this source tree. + +# Authors: Zeyu Ma + +import sys +from ctypes import CDLL, POINTER, RTLD_LOCAL, c_double, c_float, c_int32, c_bool +from pathlib import Path + +from numpy import ascontiguousarray as AC + + +# note: size of x should not exceed maximum +def AsInt(x): + return x.ctypes.data_as(POINTER(c_int32)) +def AsDouble(x): + return x.ctypes.data_as(POINTER(c_double)) +def AsFloat(x): + return x.ctypes.data_as(POINTER(c_float)) +def AsBool(x): + return x.ctypes.data_as(POINTER(c_bool)) + +def register_func(me, dll, name, argtypes=[], restype=None, caller_name=None): + if caller_name is None: caller_name = name + setattr(me, caller_name, getattr(dll, name)) + func = getattr(me, caller_name) + func.argtypes = argtypes + func.restype = restype + +def load_cdll(path): + return CDLL(Path(sys.path[-1]) / path, mode=RTLD_LOCAL) \ No newline at end of file diff --git a/ocmesher/utils/timer.py b/ocmesher/utils/timer.py new file mode 100644 index 0000000..045aba6 --- /dev/null +++ b/ocmesher/utils/timer.py @@ -0,0 +1,37 @@ +# Copyright (c) Princeton University. +# This source code is licensed under the BSD 3-Clause license found in the LICENSE file in the root directory of this source tree. + +# Authors: Zeyu Ma + +from datetime import datetime +import os + +import psutil + + +class Timer: + + def __init__(self, desc, disable_timer=False): + self.disable_timer = disable_timer + if self.disable_timer: + return + self.name = f'[{desc}]' + + def __enter__(self): + if self.disable_timer: + return + self.start = datetime.now() + + + def __exit__(self, exc_type, exc_val, traceback): + if self.disable_timer: + return + self.end = datetime.now() + self.duration = self.end - self.start # timedelta + if exc_type is None: + process = psutil.Process(os.getpid()) + print(f'{self.name} finished in {str(self.duration)} with memory usage {process.memory_info().rss / 1024**3} GB') + else: + print(f'{self.name} failed with {exc_type}') + + diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..638f725 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,6 @@ +numpy +vnoise +trimesh +gin-config +tqdm +psutil \ No newline at end of file