Python 计算Mesh顶点法向量

太过爱你忘了你带给我的痛 2022-10-07 15:50 202阅读 0赞

一、将顶点转为点云后利用点云计算法向量的方法

首先将.obj文件读入并将顶点转为点云存储为.pcd文件：

import os
    import numpy as np
    
    
    def read_obj(obj_path):
        with open(obj_path) as file:
            points = []
            faces = []
            while 1:
                line = file.readline()
                if not line:
                    break
                strs = line.split(" ")
                if strs[0] == "v":
                    points.append((float(strs[1]), float(strs[2]), float(strs[3])))
                if strs[0] == "f":
                    faces.append((int(strs[1]), int(strs[2]), int(strs[3])))
        points = np.array(points)
        faces = np.array(faces)
    
        return points, faces
    
    
    def save_pcd(filename, pcd):
        num_points = np.shape(pcd)[0]
        f = open(filename, 'w')
        f.write('# .PCD v0.7 - Point Cloud Data file format \nVERSION 0.7 \nFIELDS x y z \nSIZE 4 4 4 \nTYPE F F F \nCOUNT 1 1 1 \n')
        f.write('WIDTH {} \nHEIGHT 1 \nVIEWPOINT 0 0 0 1 0 0 0 \n'.format(num_points))
        f.write('POINTS {} \nDATA ascii\n'.format(num_points))
        for i in range(num_points):
            new_line = str(pcd[i,0]) + ' ' + str(pcd[i,1]) + ' ' + str(pcd[i,2]) + '\n'
            f.write(new_line)
        f.close()
    
    
    if __name__ == "__main__":
        objfile = "000000.obj"
        points, _ = read_obj(objfile)
        pcdfile = objfile.replace('.obj', '.pcd')
        save_pcd(pcdfile, points)

然后计算每个顶点的法向量：

import open3d as o3d
    import os
    import numpy as np
    
    path = "000000.pcd"
    normalPath = path.replace(".pcd", "_normal.pcd")
    print(path)
    print(normalPath)
    
    print("Load a pcd point cloud, print it, and render it")
    pcd = o3d.io.read_point_cloud(path)
    pcd.paint_uniform_color([0.5, 0.5, 0.5])
    # print(np.asarray(pcd.points))
    # o3d.visualization.draw_geometries([pcd], "Open3D origin points", width=800, height=600, left=50, top=50,
    #                                   point_show_normal=False, mesh_show_wireframe=False,
    #                                   mesh_show_back_face=False)
    
    # print("Downsample the point cloud with a voxel of 0.002")
    # downpcd = pcd.voxel_down_sample(voxel_size=0.002)
    # print(downpcd)
    # o3d.visualization.draw_geometries([downpcd], "Open3D downsample points", width=800, height=600, left=50, top=50,
    #                                   point_show_normal=False, mesh_show_wireframe=False,
    #                                   mesh_show_back_face=False)
    
    print("Compute the normal of the origin point cloud")
    # downpcd.estimate_normals(
    #     search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.01, max_nn=20))
    print(pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamKNN(knn=20)))
    # downpcd.estimate_normals(
    #     search_param=o3d.geometry.KDTreeSearchParamRadius(radius=0.01))
    
    # o3d.visualization.draw_geometries([pcd], "Open3D normal estimation", width=800, height=600, left=50, top=50,
    #                                   point_show_normal=True, mesh_show_wireframe=False,
    #                                   mesh_show_back_face=True)
    # print("Print a normal vector of the 0th point")
    # print(pcd.normals[0])
    
    print("Print the normal vectors of the first 10 points")
    print(np.asarray(pcd.normals)[:10, :])
    # pcd.normals = o3d.utility.Vector3dVector(normals)
    np.save(file="verts_normals.npy", arr=np.asarray(pcd.normals))
    # normals = o3d.np.asarray(pcd.normals)
    # print(normals)
    
    normal_point = o3d.utility.Vector3dVector(pcd.normals)
    normals = o3d.geometry.PointCloud()
    normals.points = normal_point
    normals.paint_uniform_color((0, 1, 0))
    # o3d.visualization.draw_geometries([pcd, normals], "Open3D noramls points", width=800, height=600, left=50, top=50,
    #                                   point_show_normal=False, mesh_show_wireframe=False,
    #                                   mesh_show_back_face=False)
    o3d.io.write_point_cloud(normalPath, normals) 
    print('OK')

二、直接计算三角面片（faces）法向量进而求顶点法向量

# Calculate the normal vector for each face in the mesh,and #
    # then compute the normalized average of the normal vector  #
    # for all faces connected by vertices as its normal vector. #
    import numpy as np
    from numpy.core.numeric import ones
    from OpenGL.GL import *
    
    #This is a function to read obj file
    def read_obj(obj_path):
        with open(obj_path) as file:
            points = []
            faces = []
            while 1:
                line = file.readline()
                if not line:
                    break
                strs = line.split(" ")
                if strs[0] == "v":
                    points.append((float(strs[1]), float(strs[2]), float(strs[3])))
                if strs[0] == "f":
                    faces.append((int(strs[1]), int(strs[2]), int(strs[3])))
        points = np.array(points)
        faces = np.array(faces)
        return points, faces
    
    #This a function to normalize normal vector
    def normalize_v3(arr):
        ''' Normalize a numpy array of 3 component vectors shape=(n,3) '''
        lens = np.sqrt( arr[:,0]**2 + arr[:,1]**2 + arr[:,2]**2 )
        arr[:,0] /= lens
        arr[:,1] /= lens
        arr[:,2] /= lens                
        return arr
    
    # Then to render this in OpenGL
    # def render(va,no):
    #     glEnableClientState( GL_VERTEX_ARRAY )
    #     glEnableClientState( GL_NORMAL_ARRAY )
    #     glVertexPointer( 3, GL_FLOAT, 0, va )
    #     glNormalPointer( GL_FLOAT,    0, no )
    #     glDrawArrays(GL_TRIANGLES,    0, len(va) )
        
    if __name__ == "__main__":
        vertices, faces = read_obj("000001.obj")
        onesT = ones(faces.shape,dtype=np.int) 
        faces = faces-onesT  # Comment this line if faces contains a vertex number starting with 0
        #Create a zeroed array with the same type and shape as our vertices i.e., per vertex normal
        norm = np.zeros( vertices.shape, dtype=vertices.dtype )
        #Create an indexed view into the vertex array using the array of three indices for triangles
        tris = vertices[faces]#(13776,3,3)
        #Calculate the normal for all the triangles, by taking the cross product of the vectors v1-v0, and v2-v0 in each triangle             
        n = np.cross( tris[::,1 ] - tris[::,0]  , tris[::,2 ] - tris[::,0] )
        # n is now an array of normals per triangle. The length of each normal is dependent the vertices, 
        # we need to normalize these, so that our next step weights each normal equally.
        n = normalize_v3(n)
        # now we have a normalized array of normals, one per triangle, i.e., per triangle normals.
        # But instead of one per triangle (i.e., flat shading), we add to each vertex in that triangle, 
        # the triangles' normal. Multiple triangles would then contribute to every vertex, so we need to normalize again afterwards.
        # The cool part, we can actually add the normals through an indexed view of our (zeroed) per vertex normal array
        norm[ faces[:,0] ] += n
        norm[ faces[:,1] ] += n
        norm[ faces[:,2] ] += n
        result = normalize_v3(norm)
        np.save(file="verts_normals1.npy", arr=result)
        print(result)
        # then we create an indexed view into our vertices and normals
        # va = vertices[ faces ]
        # no = norm[ faces ]     
        # render(va,no)

结果对比  
左图方法一结果，右图方法二结果：

Reference

https://sites.google.com/site/dlampetest/python/calculating-normals-of-a-triangle-mesh-using-numpy