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

160 阅读10分钟

目的

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

裁剪顶点.png

数据准备

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

无人机影像的空三建模.png
由于模型场景很大,本文采用截取一个单体建筑的方式来对场景进行裁剪,后期需要裁剪所有建筑,则需要迭代计算并优化效率(黑框全选为需要裁剪的建筑)。

image.png shp数据为建筑面数据

image.png

算法思路

由于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}")

image.png

image.png