基于shp数据实现obj模型的切割（模型单体化思路）

目的

在三维仿真的场景中，无人机影像的建模是一种快速可靠的对大场景地形三维重现的方式，但是无人机影像的建模无法对单个场景的属性进行分析，如对单个楼层、居民区进行属性分析查询。在GIS中，我们通过shp数据描述空间信息，这些数据通常以点线面的形式，即概括了地物的空间关系也包含地物的属性信息。通过shp数据和无人机影像建模的结合，可以将大范围三维场景的模型进行分割，使之成为一个个独立的单元并且赋予属性信息，本文通过一个小实验来实现shp面对obj模型中某个单个建筑的分割。

数据准备

采用网上公开无人机影像数据，通过contextcapture对影像进行建模，生成倾斜摄影数据和obj模型，本文采用obj模型进行切割。

由于模型场景很大，本文采用截取一个单体建筑的方式来对场景进行裁剪，后期需要裁剪所有建筑，则需要迭代计算并优化效率（黑框全选为需要裁剪的建筑）。

shp数据为建筑面数据

算法思路

由于obj模型包含顶点和纹理，而shp数据基本处于二维平面，那么需要考虑顶点落在shp面数据xy平面内的数据，同时，一个面最少由三个顶点构成，如果一个面的部分在shp范围内，部分在范围外，则需要重新分割，生成新的顶点并构建网格。

obj数据的格式

通常简单的obj数据包含顶点信息，纹理信息，面信息，有的obj还包含一个mtl文件
其中常见的obj数据的组成形式为：

• 首行：mtllib *.mtl 表示使用哪个mtl文件,以mtllib开头
• 顶点坐标： v x y z 表示一个顶点的坐标,以v开头
• 纹理坐标: vt u v 表示一个纹理的坐标，以vt开头
• 引用的材质：usemtl * 表示引用mtl文件的哪部分纹理，以usemtl开头
• 面索引：f v_i1/u_i1 v_i2/u_i2 v_i3/u_i3 ，以f开头，分别记录 顶点的序号和纹理的序号，序号从1开始，一个面由三个顶点组成，所以有三个顶点序号和纹理序号
  

mtl文件的格式

mtl记录了纹理的一些配置信息，主要有：

• newmtl *: 创建一个材质，材质名为*，对应obj中的usemtl *
• ka * * *: 环境颜色
• kd * * *：漫反射颜色
• d *: 透明度
• Ns: *: 高光指数
• illum: * : 光照模型
• map_kd: *.jpg:纹理图片的名称

以上可以得出一个结论，一个mtl可以创建多个材质，那么在obj的面索引前面，可以多次引入材质信息，这个材质信息下面对应的面索引的顶点和纹理坐标，使用该材质

解析obj数据

实验采用python语言，通过pywavefront库解析obj模型

model = pywavefront.Wavefront(r"D:\data\osgb\mec\Production_1 (3)\Data\Model\2.obj",collect_faces=True)

解析后的数据存在在model对象上，那么就需要知道model对象对应的成员有哪些：

• materials: 记录了材质信息
  • vertices: 一个坐标数组，五个数据一组，前两位表示纹理坐标，后三位表示顶点坐标
• mtllibs: 模型使用的mtl文件名，以数组形式存储
• vertices: 顶点坐标数组
• meshes: 格网列表，表示面数据索引集合

由于每个材质对应了一部分三角格网数据，格网数据对应了相应的顶点和纹理坐标数据，因此在解析数据的时候以材质进行分类，先记录每个材质的相关信息，在根据材质读取对应的顶点坐标和纹理坐标，进行裁剪和划分。

for mesh_name in model.materials:
    model_mesh = model.materials[mesh_name]
    if model_mesh.texture == None:
        continue
    newmtl = model_mesh.name
    #环境颜色
    Ka = model_mesh.ambient
    #漫反射颜色
    Kd = model_mesh.diffuse
    #透明度
    d = model_mesh.transparency
    #高光指数
    Ns = model_mesh.specular
    #光照模型
    illum = model_mesh.illumination_model
    #材质
    map_Kd = model_mesh.texture.file_name
    file_path = model_mesh.texture.path
    vertices = model_mesh.vertices
    new_faces = []
    new_mtls = []
    split_faces = []
    split_mtls = []
    flag = False
	#获取顶点和纹理坐标
	#检查面的顶点是否在多边形内
	#实现面的分割逻辑
            
    if flag:
        materials = {
            "Ka": Ka,
            "Kd": Kd,
            "d": d,
            "Ns": Ns,
            "illum": illum,
            "map_Kd": map_Kd,
            "file_path": file_path,
            "new_faces": new_faces,
            "new_mtls": new_mtls,
            "split_faces": split_faces,
            "split_mtls": split_mtls,
            "newmtl":newmtl
        } 
        newObj.append(materials)

        

如上所示，首先根据材质读取材质信息，然后检查面的顶点是否在多边形内，如果面的所有顶点都在多边形内，那么只需要保留这个面的顶点和纹理坐标就行；如果面的部分顶点在多边形内，部分不在，则需要进行分割。
由于顶点和纹理信息存储在model_mesh.vertices中，只需要遍历该数据，即可获取对应的顶点坐标和纹理坐标。

  for index in range(int(len(vertices) /5/3) ):
        #计算每一个面
        projected_vertices = []
        new_mtl = []
        for i in range(3):
            projected_vertices.append([vertices[15*index+2 + 5*i],vertices[15*index+3+5*i],vertices[15*index+4+5*i]])
            new_mtl.append([vertices[15*index+0 + 5*i],vertices[15*index+1 + 5*i]])
        #检查面的顶点是否在多边形内
       
        #实现面的分割逻辑
           

进行分割和裁剪面的构建

​ 多边形的构建采用shapely库，由于是对一个shp要素进行裁剪，因此可以直接通过shapely构建一个多边形要素作为该建筑轮廓，实际分割中也可以采取读取shp要素，依次构建的方式。

# 定义多边形
polygon = shapely.geometry.Polygon([(12730175.00000000 ,3571689.50000000),(12730175.00000000 ,3571656.25000000),(12730189.00000000,3571656.25000000),(12730189.00000000 ,3571625.25000000),(12730175.00000000,3571625.00000000),(12730175.00000000,3571628.75000000),(12730113.00000000,3571628.50000000),(12730113.00000000,3571626.00000000),(12730099.00000000,3571626.00000000),(12730099.00000000,3571651.25000000),(12730115.00000000,3571651.25000000),(12730115.00000000,3571647.50000000),(12730159.00000000, 3571647.50000000),(12730159.00000000, 3571698.25000000),(12730175.00000000 ,3571698.25000000 )
])

​ 判断顶点是否在多边形内

#函数: 判断点是否在多边形内
def is_inside_polygon(point,polygon):
    return polygon.contains(shapely.geometry.Point(point))

​ 检查面的顶点是否都在多边形内，如果都在，则将对应顶点和纹理坐标追加到数组中

 #检查面的顶点是否在多边形内
        if all(is_inside_polygon(v,polygon) for v in projected_vertices):
            new_faces.extend([projected_vertices])
            new_mtls.extend([new_mtl])
            flag = True

​ 如果只有部分顶点在多边形内，则对面进行分割，并将分割的面和纹理追加进数组

 elif any(is_inside_polygon(v,polygon) for v in projected_vertices):
            #实现面的分割逻辑
            new_splits,split_mtls = split_face(projected_vertices,new_mtl,polygon)
            split_faces.extend(new_splits)
            split_mtls.extend(split_mtls)
            flag = True

​ 面的分割逻辑需要检查每一条边，判断该边是否和多边形有交点，如果有交点，则分别将交点和在多边形内部点追加进数组

#函数：实现对面的分割
def split_face(face,new_mtls,polygon):
    #存储交点和保留的顶点
    intersections = []
    intersecmtls = []
    inside_vertices = []
    inside_mtls = []
    #检查每条边
    for i in range(len(face)):
        start_vertex = face[i]
        end_vertex = face[(i + 1) % len(face)]
        start_mtl = new_mtls[i]
        end_mtl = new_mtls[(i+1)%len(new_mtls)]
        if is_inside_polygon(start_vertex,polygon):
            inside_vertices.append(start_vertex)
            inside_mtls.append(start_mtl)

        #计算交点
        intersection,intermtl = calculate_intersection((start_vertex,end_vertex),start_mtl,end_mtl,polygon)
        if intersection:
            intersections.append(intersection)
            intersecmtls.append(intermtl)
    
	   #创建新的面
    return new_faces,split_mtls

​ 由于判断交点和内部点是在xy平面上进行的，因此在计算出交点的同时，需要通过线性插值的方式，计算出该点的Z值以及该点的纹理坐标

# 计算线段和多边形交点
def calculate_intersection(edge,start_mtl,end_mtl, polygon):
    line = LineString([edge[0] ,edge[1] ])
    intersection = line.intersection(polygon)

    if intersection.is_empty:
        return None,None
    else:
        if isinstance(intersection,Point):

            # 获取二维交点的坐标
            x, y = intersection.x, intersection.y
        elif isinstance(intersection, LineString):
            # 如果交点是一条线，你可以选择线上的一个点，例如第一个点
            x, y = intersection.coords[0][:2]
        elif isinstance(intersection, MultiPoint):
            # 如果有多个交点，选择一个点，例如第一个点
            x, y = intersection[0].x, intersection[0].y
        else: 
            return None,None

        # 线性插值计算 z 坐标
        z1, z2 = edge[0][2], edge[1][2]
        x1, y1 = edge[0][0], edge[0][1]
        x2, y2 = edge[1][0], edge[1][1]

        # 计算两端点之间的距离
        dist1 = ((x - x1)**2 + (y - y1)**2)**0.5
        dist2 = ((x - x2)**2 + (y - y2)**2)**0.5
        total_dist = dist1 + dist2

        # 根据距离进行线性插值
        (u1,v1) = start_mtl
        (u2,v2) = end_mtl
        z = (z1 * dist2 + z2 * dist1) / total_dist
        u = (u1 * dist2 + u2 * dist1) / total_dist
        v = (v1 * dist2 + v2 * dist1) / total_dist

        return (x, y, z),(u,v)

​ 在得出分割的点和内部点之后，即可根据这些数据构建新的格网面，此时需要分为两种场景，当内部点只有一个的时候，说明有两个分割点，那么可以构建出一个新的多边形面；当内部点有两个的时候，那么分割点也有两个，那么可以构建出两个格网面

def create_new_faces(internal_vertices,intersections,inside_mtls,intersecmtls):

    new_faces = []
    new_mtls = []
    #当一个顶点在内部
    if len(internal_vertices) == 1:
        #新的三角面由一个内部顶点和两个交点组成
        new_faces.append([internal_vertices[0],intersections[0],intersections[1]])
        new_mtls.append([inside_mtls[0],intersecmtls[0],intersecmtls[1]])
    #当两个顶点在内部
    elif len(internal_vertices) == 2:
        #第一个新的三角面由两个内部顶点和一个交点组成
        new_faces.append([internal_vertices[0],internal_vertices[1],intersections[0]])
        new_mtls.append([inside_mtls[0],inside_mtls[1],intersecmtls[0]])

        #第二个新的三角面由一个内部顶点、一个交点和另一个交点组成
        new_faces.append([internal_vertices[1],intersections[0],intersections[1]])
        new_mtls.append([inside_mtls[1],inside_mtls[0],intersecmtls[1]])

    return new_faces,new_mtls

模型的写入

​ 将新的顶点和纹理坐标及面索引重新写入obj同样需要遵守上述的准则，需要写入mtl文件和obj两个文件，同时mtl中所引用的图片得复制到对应的目录下

"""
将顶点和面的数据写入到一个新的OBJ文件中。
:param filename: 新OBJ文件的名称。
:param vertices: 顶点列表，每个顶点是一个(x, y, z)元组。
:param faces: 面列表，每个面是一组顶点索引。
"""
def write_obj_file(folder_path):
    if not os.path.exists(folder_path):
        os.makedirs(folder_path)
    # 构建OBJ和MTL文件的完整路径
    folder_name = os.path.basename(folder_path)
    obj_file_path = os.path.join(folder_path, f"{folder_name}.obj")
    mtl_file_path = os.path.join(folder_path, f"{folder_name}.mtl")
    with open(obj_file_path,'w') as file:
        file.write(f"mtllib {folder_name}.mtl\n")
        # 写入顶点
        for materials in newObj:
            for face in materials['new_faces']:
                file.write(f"v {face[0][0] } {face[0][1] } {face[0][2] }\n")
                file.write(f"v {face[1][0] } {face[1][1] } {face[1][2] }\n")
                file.write(f"v {face[2][0] } {face[2][1] } {face[2][2] }\n")
            for face in materials['split_faces']:
                file.write(f"v {face[0][0] } {face[0][1] } {face[0][2] }\n")
                file.write(f"v {face[1][0] } {face[1][1] } {face[1][2] }\n")
                file.write(f"v {face[2][0] } {face[2][1] } {face[2][2] }\n")
        #写入材质
        for materials in newObj:
            for uv in materials['new_mtls']:
                file.write(f"vt {uv[0][0] } {uv[0][1] } \n")
                file.write(f"vt {uv[1][0] } {uv[1][1] } \n")
                file.write(f"vt {uv[2][0] } {uv[2][1] } \n")
            for uv in materials['split_mtls']:
                file.write(f"vt {uv[0][0] } {uv[0][1] } \n")
                file.write(f"vt {uv[1][0] } {uv[1][1] } \n")
                file.write(f"vt {uv[2][0] } {uv[2][1] } \n")
        #写入面索引
        index = 1
        for materials in newObj:
            file.write(f"usemtl {materials['newmtl']}\n")
            for idx in range(len(materials['new_faces'])):
                face_line = ' '.join([str(3*idx+i + index)+"/"+str(3*idx+i + index) for i in range(3)])
                file.write(f"f {face_line}\n")
            index = index+ 3*len(materials['new_faces'])
            for idx in range(len(materials['split_faces'])):
                face_line = ' '.join([str(3*idx+i + index)+"/"+str(3*idx+i + index) for i in range(3)])
                file.write(f"f {face_line}\n")
            index=len(materials['split_faces'])*3+index
    with open(mtl_file_path,'w') as mtl_file:
        for materials in newObj:
            mtl_file.write(f"newmtl {materials['newmtl']}\n");
            mtl_file.write(f"Ka {materials['Ka'][0]} {materials['Ka'][1]} {materials['Ka'][2]}\n");
            mtl_file.write(f"Kd {materials['Kd'][0]} {materials['Kd'][1]} {materials['Kd'][2]}\n");
            mtl_file.write(f"d {materials['d']}\n");
            mtl_file.write(f"Ns {materials['Ns'][0]}\n");
            mtl_file.write(f"illum {materials['illum']}\n");
            mtl_file.write(f"map_Kd {materials['map_Kd']}\n");
            file_path = materials["file_path"]
            copy_photo(file_path, folder_path)
         

​ 图片的复制是一个简单的io操作，判断要复制的图片是否存在，然后根据地址进行复制即可

#图片复制
def copy_photo(src_file_path, dst_folder_path):
    """
    将照片从源路径复制到目标文件夹。
    :param src_file_path: 源照片的完整路径。
    :param dst_folder_path: 目标文件夹的路径。
    """
    # 确保目标文件夹存在
    if not os.path.exists(dst_folder_path):
        os.makedirs(dst_folder_path)

    # 构建目标文件的路径
    dst_file_path = os.path.join(dst_folder_path, os.path.basename(src_file_path))

    # 复制文件
    shutil.copy(src_file_path, dst_file_path)

    print(f"照片已复制到: {dst_file_path}")

源码下载

源码下载