Add files via upload

class_3d_model.py

@@ -0,0 +1,67 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Anton Wang
+"""
+
+# 3D model class
+
+import numpy as np
+
+class a_3d_model:
+    def __init__(self, filepath):
+        self.model_filepath=filepath
+        self.load_obj_file()
+        self.calculate_plane_equations()
+        self.calculate_Q_matrices()
+
+    def load_obj_file(self):
+        with open(self.model_filepath) as file:
+            self.points = []
+            self.faces = []
+            while 1:
+                line = file.readline()
+                if not line:
+                    break
+                strs = line.split(" ")
+                if strs[0] == "v":
+                    self.points.append((float(strs[1]), float(strs[2]), float(strs[3])))
+                if strs[0] == "f":
+                    self.faces.append((int(strs[1]), int(strs[2]), int(strs[3])))
+        self.points=np.array(self.points)
+        self.faces=np.array(self.faces)
+        self.number_of_points=self.points.shape[0]
+        self.number_of_faces=self.faces.shape[0]
+        edge_1=self.faces[:,0:2]
+        edge_2=self.faces[:,1:]
+        edge_3=np.concatenate([self.faces[:,:1], self.faces[:,-1:]], axis=1)
+        self.edges=np.concatenate([edge_1, edge_2, edge_3], axis=0)
+        unique_edges_trans, unique_edges_locs=np.unique(self.edges[:,0]*(10**10)+self.edges[:,1], return_index=True)
+        self.edges=self.edges[unique_edges_locs,:]
+
+    def calculate_plane_equations(self):
+        self.plane_equ_para = []
+        for i in range(0, self.number_of_faces):
+            # solving equation ax+by+cz+d=0, a^2+b^2+c^2=1
+            # set d=-1, give three points (x1, y1 ,z1), (x2, y2, z2), (x3, y3, z3)
+            point_1=self.points[self.faces[i,0]-1, :]
+            point_2=self.points[self.faces[i,1]-1, :]
+            point_3=self.points[self.faces[i,2]-1, :]
+            point_mat=np.array([point_1, point_2, point_3])
+            abc=np.matmul(np.linalg.inv(point_mat), np.array([[1],[1],[1]]))
+            self.plane_equ_para.append(np.concatenate([abc.T, np.array(-1).reshape(1, 1)], axis=1)/(np.sum(abc**2)**0.5))
+        self.plane_equ_para=np.array(self.plane_equ_para)
+        self.plane_equ_para=self.plane_equ_para.reshape(self.plane_equ_para.shape[0], self.plane_equ_para.shape[2])
+
+    def calculate_Q_matrices(self):
+        self.Q_matrices = []
+        for i in range(0, self.number_of_points):
+            point_index=i+1
+            # each point is the solution of the intersection of a set of planes
+            # find the planes for point_index
+            face_set_index=np.where(self.faces==point_index)[0]
+            Q_temp=np.zeros((4,4))
+            for j in face_set_index:
+                p=self.plane_equ_para[j,:]
+                p=p.reshape(1, len(p))
+                Q_temp=Q_temp+np.matmul(p.T, p)
+            self.Q_matrices.append(Q_temp)

class_mesh_simplify.py

@@ -0,0 +1,308 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Anton Wang
+"""
+
+import numpy as np
+import sys
+
+from class_3d_model import a_3d_model
+
+# Mesh simplification calss
+class mesh_simplify(a_3d_model):
+    def __init__(self, input_filepath, threshold, simplify_ratio):
+        if simplify_ratio>1 or simplify_ratio<=0:
+            sys.exit('Error: simplification ratio should be (0<r<=1).')
+        if threshold<0:
+            sys.exit('Error: threshold should be (>=0).')
+        super().__init__(input_filepath)
+        print('Import model: '+str(input_filepath))
+        self.t=threshold
+        self.ratio=simplify_ratio
+
+    # Select all valid pairs.
+    def generate_valid_pairs(self):
+        self.dist_pairs = []
+        for i in range(0, self.number_of_points):
+            current_point_location=i+1
+            current_point=self.points[i,:]
+            current_point_to_others_dist=(np.sum((self.points-current_point)**2,axis=1))**0.5
+            valid_pairs_location=np.where(current_point_to_others_dist<=self.t)[0]+1
+            valid_pairs_location=valid_pairs_location.reshape(len(valid_pairs_location),1)
+            current_valid_pairs=np.concatenate([current_point_location*np.ones((valid_pairs_location.shape[0],1)),valid_pairs_location],axis=1)
+            if i==0:
+                self.dist_pairs=current_valid_pairs
+            else:
+                self.dist_pairs=np.concatenate([self.dist_pairs, current_valid_pairs], axis=0)
+        self.dist_pairs=np.array(self.dist_pairs)
+        find_same=self.dist_pairs[:,1]-self.dist_pairs[:,0]
+        find_same_loc=np.where(find_same==0)[0]
+        self.dist_pairs=np.delete(self.dist_pairs, find_same_loc, axis=0)
+
+        if self.dist_pairs.size > 0:
+            self.valid_pairs=np.concatenate([self.edges,self.dist_pairs],axis=0)
+            self.valid_pairs=np.array(self.valid_pairs, dtype=int)
+        else:
+            self.valid_pairs=self.edges
+
+        find_same=self.valid_pairs[:,1]-self.valid_pairs[:,0]
+        find_same_loc=np.where(find_same==0)[0]
+        self.valid_pairs=np.delete(self.valid_pairs, find_same_loc, axis=0)
+
+        unique_valid_pairs_trans, unique_valid_pairs_loc=np.unique(self.valid_pairs[:,0]*(10**10)+self.valid_pairs[:,1], return_index=True)
+        self.valid_pairs=self.valid_pairs[unique_valid_pairs_loc,:]
+
+    # Compute the optimal contraction target v_opt for each valid pair (v1, v2)
+    # The error v_opt.T*(Q1+Q2)*v_pot of this target vertex becomes the cost of contracting that pair.
+    # Place all the pairs in a heap keyed on cost with the minimum cost pair at the top
+    def calculate_optimal_contraction_pairs_and_cost(self):
+        self.v_optimal = []
+        self.cost = []
+        number_of_valid_pairs=self.valid_pairs.shape[0]
+        for i in range(0, number_of_valid_pairs):
+            current_valid_pair=self.valid_pairs[i,:]
+            v_1_location=current_valid_pair[0]-1
+            v_2_location=current_valid_pair[1]-1
+            # find Q_1
+            Q_1=self.Q_matrices[v_1_location]
+            # find Q_2
+            Q_2=self.Q_matrices[v_2_location]
+            Q=Q_1+Q_2
+            Q_new=np.concatenate([Q[:3,:], np.array([0,0,0,1]).reshape(1,4)], axis=0)
+            if np.linalg.det(Q_new)>0:
+                current_v_opt=np.matmul(np.linalg.inv(Q_new),np.array([0,0,0,1]).reshape(4,1))
+                current_cost=np.matmul(np.matmul(current_v_opt.T, Q), current_v_opt)
+                current_v_opt=current_v_opt.reshape(4)[:3]
+            else:
+                v_1=np.append(self.points[v_1_location,:],1).reshape(4,1)
+                v_2=np.append(self.points[v_2_location,:],1).reshape(4,1)
+                v_mid=(v_1+v_2)/2
+                delta_v_1=np.matmul(np.matmul(v_1.T, Q), v_1)
+                delta_v_2=np.matmul(np.matmul(v_2.T, Q), v_2)
+                delta_v_mid=np.matmul(np.matmul(v_mid.T, Q), v_mid)
+                current_cost=np.min(np.array([delta_v_1, delta_v_2, delta_v_mid]))
+                min_delta_loc=np.argmin(np.array([delta_v_1, delta_v_2, delta_v_mid]))
+                current_v_opt=np.concatenate([v_1,v_2,v_mid],axis=1)[:,min_delta_loc].reshape(4)
+                current_v_opt=current_v_opt[:3]
+            self.v_optimal.append(current_v_opt)
+            self.cost.append(current_cost)
+
+        self.v_optimal=np.array(self.v_optimal)
+        self.cost=np.array(self.cost)
+        self.cost=self.cost.reshape(self.cost.shape[0])
+
+        cost_argsort=np.argsort(self.cost)
+        self.valid_pairs=self.valid_pairs[cost_argsort,:]
+        self.v_optimal=self.v_optimal[cost_argsort,:]
+        self.cost=self.cost[cost_argsort]
+
+        self.new_point=self.v_optimal[0,:]
+        self.new_valid_pair=self.valid_pairs[0,:]
+
+    # Iteratively remove the pair (v1, v2) of least cost from the heap
+    # contract this pair, and update the costs of all valid pairs involving (v1, v2).
+    # until existing points = ratio * original points
+    def iteratively_remove_least_cost_valid_pairs(self):
+        self.new_point_count=0
+        self.status_points=np.zeros(self.number_of_points)
+        self.status_faces=np.zeros(self.number_of_faces)
+        while (self.number_of_points-self.new_point_count)>=self.ratio*(self.number_of_points):
+
+            # current valid pair
+            current_valid_pair=self.new_valid_pair
+            v_1_location=current_valid_pair[0]-1 # point location in self.points
+            v_2_location=current_valid_pair[1]-1
+
+            # update self.points
+            # put the top optimal vertex(point) into the sequence of points
+            self.points[v_1_location,:]=self.new_point.reshape(1,3)
+            self.points[v_2_location,:]=self.new_point.reshape(1,3)
+
+            # set status of points
+            # 0 means no change, -1 means the point is deleted
+            # update v1, v2 to v_opt, then delete v2, keep v1
+            self.status_points[v_2_location]=-1
+
+            # set status of faces
+            # 0 means no change, -1 means the face will be deleted
+            v_1_in_faces_loc=np.where(self.faces==(v_1_location+1))
+            v_2_in_faces_loc=np.where(self.faces==(v_2_location+1))
+            v_1_2_in_one_face_loc = []
+            for item in v_2_in_faces_loc[0]:
+                if np.where(v_1_in_faces_loc[0]==item)[0].size>0:
+                    v_1_2_in_one_face_loc.append(item)
+            v_1_2_in_one_face_loc=np.array(v_1_2_in_one_face_loc)
+            if v_1_2_in_one_face_loc.size>=1:
+                self.status_faces[v_1_2_in_one_face_loc]=-1
+
+            # update self.faces
+            # points of faces involving v1 and v2 are changed accordingly
+            # set v2 to v1
+            self.faces[v_2_in_faces_loc]=v_1_location+1
+
+            # update edges
+            #edge_1=np.delete(self.faces[:,0:2], v_1_2_in_one_face_loc, axis=0)
+            #edge_2=np.delete(self.faces[:,1:], v_1_2_in_one_face_loc, axis=0)
+            #edge_3=np.delete(np.concatenate([self.faces[:,:1], self.faces[:,-1:]], axis=1), v_1_2_in_one_face_loc, axis=0)
+            #self.edges=np.concatenate([edge_1, edge_2, edge_3], axis=0)
+
+            # update self.plane_equ_para
+            v_1_2_in_faces_loc=np.unique(np.append(v_1_in_faces_loc[0], v_2_in_faces_loc[0]))
+            self.update_plane_equation_parameters(v_1_2_in_faces_loc)
+
+            # update self.Q_matrices
+            self.update_Q(current_valid_pair-1, v_1_location)
+
+            # update self.valid_pairs, self.v_optimal, and self.cost
+            self.update_valid_pairs_v_optimal_and_cost(v_1_location)
+            # re-calculate optimal contraction pairs and cost
+            self.update_optimal_contraction_pairs_and_cost(v_1_location)
+
+            if self.new_point_count%100==0:
+                print('Simplification: '+str(100*(self.number_of_points-self.new_point_count)/(self.number_of_points))+'%')
+                print('Remaining: '+str(self.number_of_points-self.new_point_count)+' points')
+                print('\n')
+
+            self.new_point_count=self.new_point_count+1
+
+        print('Simplification: '+str(100*(self.number_of_points-self.new_point_count)/(self.number_of_points+self.new_point_count))+'%')
+        print('Remaining: '+str(self.number_of_points-self.new_point_count)+' points')
+        print('End\n')
+
+    def calculate_plane_equation_for_one_face(self, p1, p2, p3):
+        # input: p1, p2, p3 numpy.array, shape: (3, 1) or (1,3) or (3, )
+        # p1 ,p2, p3 (x, y, z) are three points on a face
+        # plane equ: ax+by+cz+d=0 a^2+b^2+c^2=1
+        # return: numpy.array (a, b, c, d), shape: (1,4)
+        p1=np.array(p1).reshape(3)
+        p2=np.array(p2).reshape(3)
+        p3=np.array(p3).reshape(3)
+        point_mat=np.array([p1, p2, p3])
+        abc=np.matmul(np.linalg.inv(point_mat), np.array([[1],[1],[1]]))
+        output=np.concatenate([abc.T, np.array(-1).reshape(1, 1)], axis=1)/(np.sum(abc**2)**0.5)
+        output=output.reshape(4)
+        return output
+
+    def update_plane_equation_parameters(self, need_updating_loc):
+        # input: need_updating_loc, a numpy.array, shape: (n, ), locations of self.plane_equ_para need updating
+        for i in need_updating_loc:
+            if self.status_faces[i]==-1:
+                self.plane_equ_para[i,:]=np.array([0,0,0,0]).reshape(1,4)
+            else:
+                point_1=self.points[self.faces[i,0]-1, :]
+                point_2=self.points[self.faces[i,1]-1, :]
+                point_3=self.points[self.faces[i,2]-1, :]
+                self.plane_equ_para[i,:]=self.calculate_plane_equation_for_one_face(point_1, point_2, point_3)
+
+    def update_Q(self, replace_locs, target_loc):
+        # input: replace_locs, a numpy.array, shape: (2, ), locations of self.points need updating
+        # input: target_loc, a number, location of self.points need updating
+        face_set_index=np.where(self.faces==target_loc+1)[0]
+        Q_temp=np.zeros((4,4))
+
+        for j in face_set_index:
+            p=self.plane_equ_para[j,:]
+            p=p.reshape(1, len(p))
+            Q_temp=Q_temp+np.matmul(p.T, p)
+
+        for i in replace_locs:
+            self.Q_matrices[i]=Q_temp
+
+    def update_valid_pairs_v_optimal_and_cost(self, target_loc):
+        # input: target_loc, a number, location of self.points need updating
+
+        # processing self.valid_pairs        
+        # replace all the point indexes containing current valid pair with new point index: target_loc+1
+        v_1_loc_in_valid_pairs=np.where(self.valid_pairs==self.new_valid_pair[0])
+        v_2_loc_in_valid_pairs=np.where(self.valid_pairs==self.new_valid_pair[1])
+
+        self.valid_pairs[v_1_loc_in_valid_pairs]=target_loc+1
+        self.valid_pairs[v_2_loc_in_valid_pairs]=target_loc+1
+
+        delete_locs = []
+        for item in v_1_loc_in_valid_pairs[0]:
+            if np.where(v_2_loc_in_valid_pairs[0]==item)[0].size>0:
+                delete_locs.append(item)
+        delete_locs=np.array(delete_locs)
+
+        find_same=self.valid_pairs[:,1]-self.valid_pairs[:,0]
+        find_same_loc=np.where(find_same==0)[0]
+        if find_same_loc.size >=1:
+            delete_locs=np.append(delete_locs, find_same_loc)
+
+        # delete process for self.valid_pairs, self.v_optimal and self.cost
+        self.valid_pairs=np.delete(self.valid_pairs, delete_locs, axis=0)
+        self.v_optimal=np.delete(self.v_optimal, delete_locs, axis=0)
+        self.cost=np.delete(self.cost, delete_locs, axis=0)
+
+        # unique process for self.valid_pairs, self.v_optimal and self.cost
+        unique_valid_pairs_trans, unique_valid_pairs_loc=np.unique(self.valid_pairs[:,0]*(10**10)+self.valid_pairs[:,1], return_index=True)
+        self.valid_pairs=self.valid_pairs[unique_valid_pairs_loc,:]
+        self.v_optimal=self.v_optimal[unique_valid_pairs_loc,:]
+        self.cost=self.cost[unique_valid_pairs_loc]
+
+    def update_optimal_contraction_pairs_and_cost(self, target_loc):
+        # input: target_loc, a number, location of self.points need updating
+        v_target_loc_in_valid_pairs=np.where(self.valid_pairs==target_loc+1)[0]
+        for i in v_target_loc_in_valid_pairs:
+            current_valid_pair=self.valid_pairs[i,:]
+            v_1_location=current_valid_pair[0]-1
+            v_2_location=current_valid_pair[1]-1
+            # find Q_1
+            Q_1=self.Q_matrices[v_1_location]
+            # find Q_2
+            Q_2=self.Q_matrices[v_2_location]
+            Q=Q_1+Q_2
+            Q_new=np.concatenate([Q[:3,:], np.array([0,0,0,1]).reshape(1,4)], axis=0)
+            if np.linalg.det(Q_new)>0:                
+                current_v_opt=np.matmul(np.linalg.inv(Q_new),np.array([0,0,0,1]).reshape(4,1))
+                current_cost=np.matmul(np.matmul(current_v_opt.T, Q), current_v_opt)
+                current_v_opt=current_v_opt.reshape(4)[:3]
+            else:
+                v_1=np.append(self.points[v_1_location,:],1).reshape(4,1)
+                v_2=np.append(self.points[v_2_location,:],1).reshape(4,1)
+                v_mid=(v_1+v_2)/2
+                delta_v_1=np.matmul(np.matmul(v_1.T, Q), v_1)
+                delta_v_2=np.matmul(np.matmul(v_2.T, Q), v_2)
+                delta_v_mid=np.matmul(np.matmul(v_mid.T, Q), v_mid)
+                current_cost=np.min(np.array([delta_v_1, delta_v_2, delta_v_mid]))
+                min_delta_loc=np.argmin(np.array([delta_v_1, delta_v_2, delta_v_mid]))
+                current_v_opt=np.concatenate([v_1,v_2,v_mid],axis=1)[:,min_delta_loc].reshape(4)
+                current_v_opt=current_v_opt[:3]
+            self.v_optimal[i, :]=current_v_opt
+            self.cost[i]=current_cost
+
+        cost_argsort=np.argsort(self.cost)
+        self.valid_pairs=self.valid_pairs[cost_argsort,:]
+        self.v_optimal=self.v_optimal[cost_argsort,:]
+        self.cost=self.cost[cost_argsort]
+
+        self.new_point=self.v_optimal[0,:]
+        self.new_valid_pair=self.valid_pairs[0,:]
+
+    # Generate the simplified 3d model (points/vertices, faces)
+    def generate_new_3d_model(self):
+        point_serial_number=np.arange(self.points.shape[0])+1
+        points_to_delete_locs=np.where(self.status_points==-1)[0]
+        self.points=np.delete(self.points, points_to_delete_locs, axis=0)
+        point_serial_number=np.delete(point_serial_number, points_to_delete_locs)
+        point_serial_number_after_del=np.arange(self.points.shape[0])+1
+
+        faces_to_delete_locs=np.where(self.status_faces==-1)[0]
+        self.faces=np.delete(self.faces, faces_to_delete_locs, axis=0)
+
+        for i in point_serial_number_after_del:
+            point_loc_in_face=np.where(self.faces==point_serial_number[i-1])
+            self.faces[point_loc_in_face]=i
+
+        self.number_of_points=self.points.shape[0]
+        self.number_of_faces=self.faces.shape[0]
+
+    def output(self, output_filepath):
+        with open(output_filepath, 'w') as file_obj:
+            file_obj.write('# '+str(self.number_of_points)+' vertices, '+str(self.number_of_faces)+' faces\n')
+            for i in range(self.number_of_points):
+                file_obj.write('v '+str(self.points[i,0])+' '+str(self.points[i,1])+' '+str(self.points[i,2])+'\n')
+            for i in range(self.number_of_faces):
+                file_obj.write('f '+str(self.faces[i,0])+' '+str(self.faces[i,1])+' '+str(self.faces[i,2])+'\n')
+        print('Output simplified model: '+str(output_filepath))