From 7c9dc5682a9dac4b2f3857c96b803766b76d6d91 Mon Sep 17 00:00:00 2001 From: tatsy Date: Wed, 11 Jan 2017 19:11:26 +0900 Subject: [PATCH] Revising simplification algorithm to support boundary vertices. --- .gitignore | 2 +- demos/demo_simplify.py | 7 +- pygeop/__init__.py | 10 +- pygeop/geom2d/__init__.py | 12 +-- pygeop/geom3d/__init__.py | 16 ++-- pygeop/geom3d/objmesh.py | 191 ++++++++++++++++++++++++++++++++++++++ pygeop/geom3d/simplify.py | 58 +++++++++--- pygeop/geom3d/trimesh.py | 95 ++++++++++++------- pygeop/geom3d/vertex.py | 3 +- 9 files changed, 320 insertions(+), 74 deletions(-) create mode 100644 pygeop/geom3d/objmesh.py diff --git a/.gitignore b/.gitignore index 87a1b78..8e1224c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,6 @@ # User defined data/ignore/* - +.idea/* # Byte-compiled / optimized / DLL files __pycache__/ diff --git a/demos/demo_simplify.py b/demos/demo_simplify.py index eeab8ae..881e3f1 100644 --- a/demos/demo_simplify.py +++ b/demos/demo_simplify.py @@ -1,6 +1,11 @@ +import os +import sys import argparse -from pygeop.geom3d import * +sys.path.append(os.path.join(os.path.abspath(os.path.dirname(__file__)), os.pardir)) + +from pygeop.geom3d import TriMesh, simplify + def main(): parser = argparse.ArgumentParser(description='Demo for mesh simplification.') diff --git a/pygeop/__init__.py b/pygeop/__init__.py index f2e6aec..f0135f3 100644 --- a/pygeop/__init__.py +++ b/pygeop/__init__.py @@ -1,9 +1 @@ -import os - -current_dir = os.path.dirname(os.path.realpath(__file__)) -modules = [ m for m in os.listdir(current_dir) if os.path.isfile(os.path.join(current_dir, m)) ] -modules = filter(lambda m : m != '__init__.py' and not m.startswith('.'), modules) -for m in modules: - base = m.split('.')[0] - cmd = 'from .%s import *' % base - exec(cmd) +from .exception import * \ No newline at end of file diff --git a/pygeop/geom2d/__init__.py b/pygeop/geom2d/__init__.py index f2e6aec..098415f 100644 --- a/pygeop/geom2d/__init__.py +++ b/pygeop/geom2d/__init__.py @@ -1,9 +1,3 @@ -import os - -current_dir = os.path.dirname(os.path.realpath(__file__)) -modules = [ m for m in os.listdir(current_dir) if os.path.isfile(os.path.join(current_dir, m)) ] -modules = filter(lambda m : m != '__init__.py' and not m.startswith('.'), modules) -for m in modules: - base = m.split('.')[0] - cmd = 'from .%s import *' % base - exec(cmd) +from .line import * +from .point import * +from .utils import * \ No newline at end of file diff --git a/pygeop/geom3d/__init__.py b/pygeop/geom3d/__init__.py index f2e6aec..0fa3895 100644 --- a/pygeop/geom3d/__init__.py +++ b/pygeop/geom3d/__init__.py @@ -1,9 +1,7 @@ -import os - -current_dir = os.path.dirname(os.path.realpath(__file__)) -modules = [ m for m in os.listdir(current_dir) if os.path.isfile(os.path.join(current_dir, m)) ] -modules = filter(lambda m : m != '__init__.py' and not m.startswith('.'), modules) -for m in modules: - base = m.split('.')[0] - cmd = 'from .%s import *' % base - exec(cmd) +from .face import * +from .halfedge import * +from .simplify import simplify +from .trimesh import * +from .vector import * +from .vertex import * +from .objmesh import ObjMesh diff --git a/pygeop/geom3d/objmesh.py b/pygeop/geom3d/objmesh.py new file mode 100644 index 0000000..7cfb9b5 --- /dev/null +++ b/pygeop/geom3d/objmesh.py @@ -0,0 +1,191 @@ +import re +import numpy as np + +from .vector import Vector + +def parse_face(items): + pat = re.compile('([0-9]+)\/([0-9]+)\/([0-9]+)') + mat = pat.search(items[0]) + if mat is None: + pat = re.compile('([0-9]+)\/\/([0-9]+)') + mat = pat.search(items[0]) + + if mat is None: + pat = re.compile('([0-9]+)\/([0-9]+)') + mat = pat.search(items[0]) + + if mat is None: + pat = re.compile('([0-9]+)') + mat = pat.search(items[0]) + + indices = [ int(pat.search(it).group(1)) - 1 for it in items ] + return indices + +class ObjMesh(object): + def __init__(self, filename=None): + self._vertices = np.array([], dtype=np.float32) + self._normals = np.array([], dtype=np.float32) + self._texcoords = np.array([], dtype=np.float32) + self._indices = np.array([], dtype=np.uint32) + + if filename is not None: + self.load(filename) + + # Property for vertices + def _get_vertices(self): + return self._vertices + + def _set_vertices(self, vertices): + self._vertices = np.array(vertices, dtype=np.float32).flatten() + self.compute_normals() + + vertices = property(_get_vertices, _set_vertices) + + # Property for normals + def _get_normals(self): + return self._normals + + def _set_normals(self, normals): + self._normals = np.array(normals, dtype=np.float32).flatten() + + normals = property(_get_normals, _set_normals) + + # Property for texcoords + def _get_texcoords(self): + return self._texcoords + + def _set_texcoords(self, texcoords): + self._texcoords = np.array(texcoords, dtype=np.float32).flatten() + + texcoords = property(_get_texcoords, _set_texcoords) + + # Property for indices + def _get_indices(self): + return self._indices + + def _set_indices(self, indices): + self._indices = np.array(indices, dtype=np.uint32).flatten() + + indices = property(_get_indices, _set_indices) + + # Load mesh + def load(self, filename): + with open(filename, 'r') as f: + lines = [ l.strip() for l in f ] + lines = filter(lambda l : l != '' and not l.startswith('#'), lines) + + vertices = [] + normals = [] + texcoords = [] + indices = [] + for l in lines: + it = [ x for x in re.split('\s+', l.strip()) ] + if it[0] == 'v': + it = [ float(i) for i in it[1:] ] + vertices.append(it) + + elif it[0] == 'vt': + texcoords.append((float(it[1]), float(it[2]))) + + elif it[0] == 'vn': + it = [ float(i) for i in it[1:] ] + normals.append(it) + + elif it[0] == 'f': + it = it[1:] + indices.append(parse_face(it)) + else: + print('Unknown identifier: {}'.format(it[0])) + + if len(indices) > 0: + self.indices = np.array(indices).flatten() + + if len(vertices) > 0: + self.vertices = np.array(vertices).flatten() + + if len(normals) > 0: + self.normals = np.array(normals).flatten() + + if len(texcoords) > 0: + self.texcoords = np.array(texcoords).flatten() + + def save(self, filename): + assert len(self.vertices) > 0 + assert self.vertices.ndim == 1 + assert self.vertices.dtype == np.float32 + assert self.indices.ndim == 1 + assert self.indices.dtype == np.uint32 + + with open(filename, 'w') as fp: + # Write positions + for i in range(0, self.vertices.size, 3): + x = self.vertices[i + 0] + y = self.vertices[i + 1] + z = self.vertices[i + 2] + fp.write('v {0:.8f} {1:.8f} {2:.8f}\n'.format(x, y, z)) + + # Write normals + has_normal = False + if self.normals.size > 0: + has_normal = True + for i in range(0, self.normals.size, 3): + x = self.normals[i + 0] + y = self.normals[i + 1] + z = self.normals[i + 2] + fp.write('vn {0:.8f} {1:.8f} {2:.8f}\n'.format(x, y, z)) + + # Write texcoords + has_texcoord = False + if self.texcoords.size >0: + has_texcoord = True + for i in range(0, self.texcoords.size, 2): + x = self.texcoords[i + 0] + y = self.texcoords[i + 1] + fp.write('vt {0:.8f} {1:.8f}\n'.format(x, y)) + + # Write indices + for i in range(0, len(self.indices), 3): + i0 = self.indices[i + 0] + 1 + i1 = self.indices[i + 1] + 1 + i2 = self.indices[i + 2] + 1 + + if has_normal and has_texcoord: + fp.write('f {0}/{0}/{0} {1}/{1}/{1} {2}/{2}/{2}\n'.format(i0, i1, i2)) + + elif has_texcoord: + fp.write('f {0}/{0} {1}/{1} {2}/{2}\n'.format(i0, i1, i2)) + + elif has_normal: + fp.write('f {0}//{0} {1}//{1} {2}//{2}\n'.format(i0, i1, i2)) + + else: + fp.write('f {0} {1} {2}\n'.format(i0, i1, i2)) + + def compute_normals(self): + vectors = [ Vector(self.vertices[i + 0], self.vertices[i + 1], self.vertices[i + 2]) + for i in range(0, self.vertices.size, 3) ] + normals = [ Vector(0.0, 0.0, 0.0) for i in range(self.n_vertices()) ] + + for i in range(0, self.indices.size, 3): + i0 = self.indices[i + 0] + i1 = self.indices[i + 1] + i2 = self.indices[i + 2] + v0 = vectors[i0] + v1 = vectors[i1] + v2 = vectors[i2] + normal = (v1 - v0).cross(v2 - v0) + normals[i0] += normal + normals[i1] += normal + normals[i2] += normal + + for n in normals: + n.normalize() + + self.normals = [ (n.x, n.y, n.z) for n in normals ] + + + def n_vertices(self): + return self.vertices.size // 3 + + def n_normals(self): + return self.normals.size// 3 diff --git a/pygeop/geom3d/simplify.py b/pygeop/geom3d/simplify.py index 2373948..95fdf33 100644 --- a/pygeop/geom3d/simplify.py +++ b/pygeop/geom3d/simplify.py @@ -6,9 +6,7 @@ import numpy as np from ..exception import PygpException -from .trimesh import TriMesh from .vector import Vector -from .halfedge import Halfedge class PriorityQueue(object): def __init__(self): @@ -71,7 +69,9 @@ def progress_bar(x, total, width=50): bar[tick] = '>' bar = ''.join(bar) - print('[ {0:6.2f} % ] [ {1} ]'.format(100.0 * ratio, bar), end='\r', flush=True) + percent = min(100.0, 100.0 * ratio) + print('[ {0:6.2f} % ] [ {1} ]'.format(percent, bar), + end='\n' if tick >= width else '\r', flush=True) def compute_QEM(Q1, Q2, v1, v2): Q = np.identity(4) @@ -87,6 +87,7 @@ def compute_QEM(Q1, Q2, v1, v2): return qem, Vector(v_bar[0], v_bar[1], v_bar[2]) def simplify(mesh, ratio=0.5, remains=-1, show_progress=True): + EPS = 1.0e-6 start_time = time.clock() nv = mesh.n_vertices() @@ -98,8 +99,11 @@ def simplify(mesh, ratio=0.5, remains=-1, show_progress=True): n_remove = nv - remains # Compute matrix Q + if show_progress: + print('Computing matrix Q') + Qs = [ np.zeros((4, 4)) for i in range(nv) ] - for t in mesh.faces: + for i, t in enumerate(mesh.faces): vs = list(t.vertices()) assert len(vs) == 3 @@ -116,9 +120,21 @@ def simplify(mesh, ratio=0.5, remains=-1, show_progress=True): Qs[vs[1].index] += w * Q Qs[vs[2].index] += w * Q + if show_progress: + progress_bar(i, len(mesh.faces)) + + if show_progress: + progress_bar(len(mesh.faces), len(mesh.faces)) + # Push QEMs + if show_progress: + print('Computing QEMs') + pque = PriorityQueue() - for he in mesh.halfedges: + for i, he in enumerate(mesh.halfedges): + # if he.vertex_from.is_boundary or he.vertex_to.is_boundary: + # continue + i1 = he.vertex_from.index i2 = he.vertex_to.index v1 = he.vertex_from.position @@ -127,7 +143,13 @@ def simplify(mesh, ratio=0.5, remains=-1, show_progress=True): Q2 = Qs[i2] qem, v_bar = compute_QEM(Q1, Q2, v1, v2) pque.push(QEMNode(qem, i1, i2, v_bar)) - + + if show_progress: + progress_bar(i, len(mesh.halfedges)) + + if show_progress: + progress_bar(len(mesh.halfedges), len(mesh.halfedges)) + removed = 0 uftree = UnionFindTree(nv) while removed < n_remove: @@ -135,6 +157,7 @@ def simplify(mesh, ratio=0.5, remains=-1, show_progress=True): try: qn = pque.pop() except IndexError: + print('Target number is not reached!') break ii, jj, v_bar = qn.ii, qn.jj, qn.vec @@ -190,9 +213,13 @@ def simplify(mesh, ratio=0.5, remains=-1, show_progress=True): continue # Collapse halfedge - mesh.collapse_halfedge(v_j, v_i, v_bar) - uftree.merge(v_i.index, v_j.index) - assert v_i.index == uftree.root(v_j.index) + try: + mesh.collapse_halfedge(v_j, v_i, v_bar) + uftree.merge(v_i.index, v_j.index) + assert v_i.index == uftree.root(v_j.index) + except: + print('Error!') + pass # Check triangle shapes is_update = True @@ -206,11 +233,13 @@ def simplify(mesh, ratio=0.5, remains=-1, show_progress=True): e0 = v1 - v0 e1 = v2 - v0 - a1 = math.acos(e0.dot(e1) / (e0.norm() * e1.norm())) + c1 = e0.dot(e1) / (e0.norm() * e1.norm()) + a1 = math.acos(max(-1.0, min(c1, 1.0))) e2 = v1 - v3 e3 = v2 - v3 - a2 = math.acos(e2.dot(e3) / (e2.norm() * e3.norm())) + c2 = e2.dot(e3) / (e2.norm() * e3.norm()) + a2 = math.acos(max(-1.0, min(c2, 1.0))) if a1 + a2 > math.pi: mesh.flip_halfedge(he) @@ -228,6 +257,9 @@ def simplify(mesh, ratio=0.5, remains=-1, show_progress=True): ps = [ v.position for v in vs ] norm = (ps[1] - ps[0]).cross(ps[2] - ps[0]) w = norm.norm() + if w < EPS: + continue + norm = norm / w d = -norm.dot(ps[0]) @@ -242,13 +274,15 @@ def simplify(mesh, ratio=0.5, remains=-1, show_progress=True): assert mesh.vertices[v1.index] is not None assert mesh.vertices[v2.index] is not None + # if v1.is_boundary: continue + # if v2.is_boundary: continue if v1.degree() <= 3: continue if v2.degree() <= 3: continue Q1 = Qs[v1.index] Q2 = Qs[v2.index] qem, v_bar = compute_QEM(Q1, Q2, v1.position, v2.position) - pque.push(QEMNode(qem, i, j, v_bar)) + pque.push(QEMNode(qem, v1.index, v2.index, v_bar)) # Progress removed += 1 diff --git a/pygeop/geom3d/trimesh.py b/pygeop/geom3d/trimesh.py index 2012bd1..0290a49 100644 --- a/pygeop/geom3d/trimesh.py +++ b/pygeop/geom3d/trimesh.py @@ -2,10 +2,10 @@ import math from ..exception import PygpException -from .vector import Vector from .vertex import Vertex from .halfedge import Halfedge from .face import Face +from .objmesh import ObjMesh class TriMesh(object): def __init__(self, filename=''): @@ -18,30 +18,21 @@ def __init__(self, filename=''): self.load(filename) def load(self, filename): - with open(filename, 'r') as fp: - self.clear() - temp_vertices = [] - temp_indices = [] - for l in fp: - l = l.strip() - if l.startswith('v '): - v = [ float(it) for it in re.split('\s+', l)[1:] ] - temp_vertices.append((v[0], v[1], v[2])) - - if l.startswith('f '): - f = [ int(it) - 1 for it in re.split('\s+', l)[1:] ] - temp_indices.extend([ f[0], f[1], f[2] ]) - - unique_vertices = {} - for i in temp_indices: - v = temp_vertices[i] - - if v not in unique_vertices: - unique_vertices[v] = len(self.vertices) - self.vertices.append(Vertex(v[0], v[1], v[2])) - self.vertices[-1].index = unique_vertices[v] - - self.indices.append(unique_vertices[v]) + obj = ObjMesh(filename) + + unique_vertices = {} + for i in obj.indices: + vx = obj.vertices[i * 3 + 0] + vy = obj.vertices[i * 3 + 1] + vz = obj.vertices[i * 3 + 2] + v = (vx, vy, vz) + + if v not in unique_vertices: + unique_vertices[v] = len(self.vertices) + self.vertices.append(Vertex(v[0], v[1], v[2])) + self.vertices[-1].index = unique_vertices[v] + + self.indices.append(unique_vertices[v]) self._make_halfedge() @@ -78,21 +69,28 @@ def collapse_halfedge(self, v_from, v_to, update_position=None): reverse_halfedge = target_halfedge.opposite + # Boundary halfedge + is_boundary = v_from.is_boundary and v_to.is_boundary + if target_halfedge.face is None: + target_halfedge, reverse_halfedge = reverse_halfedge, target_halfedge + # Update v_to's halfedge - v_to.halfedge = target_halfedge.next.opposite.next + target_halfedge.vertex_to.halfedge = target_halfedge.next.opposite.next # Update halfedges of surrounding vertices target_halfedge.next.vertex_to.halfedge = target_halfedge.next.opposite - reverse_halfedge.next.vertex_to.halfedge = reverse_halfedge.next.opposite + if not is_boundary: + reverse_halfedge.next.vertex_to.halfedge = reverse_halfedge.next.opposite # Update topology he0 = target_halfedge.next.opposite he1 = target_halfedge.next.next.opposite he0.opposite, he1.opposite = he1, he0 - he2 = reverse_halfedge.next.opposite - he3 = reverse_halfedge.next.next.opposite - he2.opposite, he3.opposite = he3, he2 + if not is_boundary: + he2 = reverse_halfedge.next.opposite + he3 = reverse_halfedge.next.next.opposite + he2.opposite, he3.opposite = he3, he2 for he in v_to.halfedges(): he.vertex_from = v_to @@ -100,7 +98,8 @@ def collapse_halfedge(self, v_from, v_to, update_position=None): # Remove faces self.faces[target_halfedge.face.index] = None - self.faces[reverse_halfedge.face.index] = None + if not is_boundary: + self.faces[reverse_halfedge.face.index] = None # Delete halfedge self.halfedges[target_halfedge.index] = None @@ -111,6 +110,13 @@ def collapse_halfedge(self, v_from, v_to, update_position=None): if update_position is not None: self.vertices[v_to.index].position = update_position + def collapse_halfedge_boundary(self, target_halfedge): + reverse_halfedge = target_halfedge.opposite + if target_halfedge.face is None: + target_halfedge, reverse_halfedge = reverse_halfedge, target_halfedge + + + def flip_halfedge(self, he): rev = he.opposite @@ -239,6 +245,7 @@ def _make_halfedge(self): table[self.vertices[self.indices[i + 1]].index].append(he1) table[self.vertices[self.indices[i + 2]].index].append(he2) + # Set opposite halfedges for he0 in self.halfedges: for he1 in table[he0.vertex_to.index]: if he0.vertex_from == he1.vertex_to and \ @@ -248,8 +255,32 @@ def _make_halfedge(self): he1.opposite = he0 break - assert he0.opposite is not None + # Opposite halfedge not found + # Mark vertices as border vertices + if he0.opposite is None: + he0.vertex_from.is_boundary = True + he0.vertex_to.is_boundary = True + + he1 = Halfedge() + he1.vertex_from = he0.vertex_to + he1.vertex_to = he0.vertex_from + he1.opposite = he0 + he0.opposite = he1 + he1.vertex_from.halfedge = he1 + + self.halfedges.append(he1) + + # Process border vertices + for v in self.vertices: + if v.is_boundary: + he = v.halfedge + while True: + if he.opposite.next is None: + he.opposite.next = v.halfedge + break + + he = he.opposite.next for i, he in enumerate(self.halfedges): he.index = i diff --git a/pygeop/geom3d/vertex.py b/pygeop/geom3d/vertex.py index 347d62c..e437ffa 100644 --- a/pygeop/geom3d/vertex.py +++ b/pygeop/geom3d/vertex.py @@ -7,6 +7,7 @@ def __init__(self, x=0.0, y=0.0, z=0.0): self.z = z self.index= -1 self.halfedge = None + self.is_boundary = False def copy(self): v = Vertex(self.x, self.y, self.z) @@ -46,7 +47,7 @@ def vertices(self): return ( he.vertex_to for he in self.halfedges() ) def faces(self): - return ( he.face for he in self.halfedges() ) + return ( he.face for he in self.halfedges() if he.face is not None ) def __eq__(self, v): return self.x == v.x and self.y == v.y and self.z == v.z