-
Notifications
You must be signed in to change notification settings - Fork 31
/
pypoisson.pyx
200 lines (133 loc) · 6.63 KB
/
pypoisson.pyx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
#cython: language_level=3str
cimport numpy as np
import numpy as np
from libcpp.vector cimport vector
from libcpp.string cimport string
from libc.stdlib cimport malloc, free
from libcpp cimport bool
np.import_array()
cdef extern from "PoissonRecon_v6_13/src/PoissonReconLib.h":
cdef int PoissonReconLibMain(int argc, char* argv[])
cdef vector[double] double_data
cdef vector[int] int_data
cdef vector[double] mem_data
def poisson_reconstruction(points, normals, depth=8, full_depth=5, scale=1.1,
samples_per_node=1.0, cg_depth=0.0,
enable_polygon_mesh=False, enable_density=False,
nthreads=0, verbose=False):
"""
Python Binding of Screened Poisson Surface Reconstruction (Version 6.13)
more information, see http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version6.13/
Usage:
faces, vertices = poisson_reconstruction(points, normals, depth=10)
Parameters
----------
points: array-like
list of oriented vertices with the x-, y-, and z-coordinates of the positions
normals: array-like
list of the x-, y-, and z-coordinates of the normals
depth: Integer
This integer is the maximum depth of the tree that will be used for surface reconstruction.
Running at depth d corresponds to solving on a voxel grid whose resolution is no larger than 2^d x 2^d x 2^d.
Note that since the reconstructor adapts the octree to the sampling density, the specified reconstruction depth is only an upper bound.
The default value for this parameter is 8.
full_depth: Integer
This integer specifies the depth beyond depth the octree will be adapted.
At coarser depths, the octree will be complete, containing all 2^d x 2^d x 2^d nodes.
The default value for this parameter is 5.
scale: float
This floating point value specifies the ratio between the diameter of the cube used for reconstruction and the diameter of the samples' bounding cube.
The default value is 1.1.
samples_per_node: float
This floating point value specifies the minimum number of sample points that should fall within an octree node as the octree construction is adapted to sampling density.
For noise-free samples, small values in the range [1.0 - 5.0] can be used.
For more noisy samples, larger values in the range [15.0 - 20.0] may be needed to provide a smoother, noise-reduced, reconstruction.
The default value is 1.0.
cg_depth: Integer
This integer is the depth up to which a conjugate-gradients solver will be used to solve the linear system.
Beyond this depth Gauss-Seidel relaxation will be used.
The default value for this parameter is 0.
enable_polygon_mesh: Bool
Enabling this flag tells the reconstructor to output a polygon mesh (rather than triangulating the results of Marching Cubes).
The default value for this parameter is False.
enable_density: Bool
Enabling this flag tells the reconstructor to output the estimated depth values of the iso-surface vertices
The default value for this parameter is False.
nthreads: int
Number of omp threads to use. Default = 0 = all available threads.
verbose: Bool
Enable verbose mode.
Returns
-------
faces: array-like
faces of the reconstructed mesh
vertices: array-like
vertices of the reconstructed mesh
"""
return _poisson_reconstruction(np.ascontiguousarray(np.float64(points)),
np.ascontiguousarray(np.float64(normals)),
depth, full_depth, scale, samples_per_node,
cg_depth, enable_polygon_mesh,
enable_density, nthreads, verbose)
cdef _poisson_reconstruction(np.float64_t[:, ::1] points,
np.float64_t[:, ::1] normals,
int depth=8,
int full_depth=5,
double scale=1.10,
double samples_per_node=1.0,
double cg_depth=0.0,
bool enable_polygon_mesh=False,
bool enable_density=False,
int nthreads=0,
bool verbose=False):
cdef:
char **c_argv
string arg_depth = str(depth).encode()
string arg_full_depth = str(full_depth).encode()
string arg_scale = str(scale).encode()
string arg_samples_per_node = str(samples_per_node).encode()
string arg_cg_depth = str(cg_depth).encode()
string arg_nthreads = str(nthreads).encode()
int_data.clear()
double_data.clear()
mem_data.clear()
point_nrows, point_ncols = np.shape(points)
normal_nrows, normal_ncols = np.shape(normals)
mem_data.resize(point_ncols * point_nrows + normal_ncols * normal_nrows)
for i in range(point_nrows):
for j in range(point_ncols):
mem_data[j + i*(point_ncols + normal_ncols)] = points[i,j]
mem_data[j + point_ncols + i *(point_ncols + normal_ncols)] = normals[i,j]
args = [b"PoissonRecon", b"--in", b"none", b"--out", b"none", b"--depth", arg_depth.c_str(),
b"--fullDepth", arg_full_depth.c_str(), b"--scale", arg_scale.c_str(),
b"--samplesPerNode", arg_samples_per_node.c_str(),
b"--cgDepth", arg_cg_depth.c_str()]
if verbose == True:
args += [b"--verbose"]
if nthreads > 0:
args += [b"--threads", arg_nthreads.c_str()]
if enable_polygon_mesh:
args += [b"--polygonMesh"]
if enable_density:
args += [b"--density"]
c_argv = <char**> malloc(sizeof(char*) * len(args))
for idx, s in enumerate(args):
c_argv[idx] = s
try:
PoissonReconLibMain(len(args), c_argv)
finally:
free(c_argv)
face_cols, vertex_cols = 3, 3
face_rows = int_data.size() // face_cols
vertex_rows = double_data.size() // vertex_cols
cdef int *ptr_faces = &int_data[0]
cdef double *ptr_vertices = &double_data[0]
faces = np.zeros((face_rows*face_cols,), dtype=np.int32 )
vertices = np.zeros((vertex_rows*vertex_cols,), dtype=np.float64)
for i in range(face_rows*face_cols):
faces[i] = ptr_faces[i]
for i in range(vertex_rows*vertex_cols):
vertices[i] = ptr_vertices[i]
int_data.clear()
double_data.clear()
return faces.reshape(face_rows,face_cols), vertices.reshape(vertex_rows,vertex_cols)