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

首先将.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