o3d_shapes (转)

这一章介绍怎么用顶点数组创建一个3D 模型,如何创建一个shape 对象,缓冲(buffers) ,域(fields) 等等。

 


    由于要定义这个3D 模型的每个顶点,然后存入顶点数组,所以这章不会画出一个比较复杂的3D 模型,我们只是画一个立方体来说明如何创建一个3D 模型,如果对于一个复杂的3D 模型还是一个个顶点画的话,只能说太牛逼了,那时候就要用到3dmax ,maya 等软件来建模了,经过传换成一定格式的文件后再在程序中导入。

 


    上一章中直接用了一个createSphere ()函数创建了一个球体,这个函数是o3d 自带的,其实也有自带的创建立方体的函数,当然这章是不会用的,否则就没东西可以讲了。我们要做的就是自己写一个。

 

 

    function createCube(material, size){

   

    这里的两个参数,前一个是这个立方体要用到的材质,材质这块下一章会讲到。Size 顾名思义就是这个立方体的大小了。

 


    在o3d 中首先要做的便是创建一个shape 对象,一个图元(primitive )对象和一个streamBank 对象,

 


var cubeShape = g_pack . createObject ( 'Shape' );
    var cubePrimitive = g_pack . createObject ( 'Primitive' );
    var streamBank = g_pack . createObject ( 'StreamBank' );

 

cubePrimitive . material = material ;      // 定义这个形状的材质
    cubePrimitive . owner = cubeShape ;            // 定义shape 和图元的关系
    cubePrimitive . streamBank = streamBank ; ·· // 定义这个shape 的streamBank

 

    其中shape 是许多图元的集合,而一个图元则是一个点或者多个点通过一定规律连接在一起组成的图形 ( 不像纸上画画最重要的是把每条线画出来,在O3d 或者说是传统的openGL 和D3D 中画一个模型是把这个模型的每个顶点的坐标定义好,然后系统便会自己把这些点连接起来) 其中点与点之间的连接可以有多种规律,比如说比较常用的每三个点通过连线组成三角形,多个点连成多个连续的三角形,每四个点组成的四边形,多个点连成的多个连续的四边形等等。

   

             下面就要定义一下图元中各个点之间连接的规律了

cubePrimitive . primitiveType = g_o3d . Primitive . TRIANGLELIST ; // 三角形绘制

cubePrimitive . numberPrimitives = 12 ;         // 这个图元中有12 个三角形
    cubePrimitive . numberVertices = 8 ;        // 有8 个顶点

 

    创建一个绘制元素

cubePrimitive . createDrawElement ( g_pack , null );

 

  接下来填充这8 个顶点的数据,在js 中这些数据都是存放在一个一维数组中的:

var positionArray = [
    - size/2 , - size/2 ,   size/2 ,   // 顶点0
      size/2 , - size/2 ,   size/2 ,   // 顶点1
    - size/2 ,   size/2 ,   size/2 ,   // 顶点2
      size/2 ,   size/2 ,   size/2 ,   // 顶点3
    - size/2 ,   size/2 , - size/2 ,   // 顶点4
      size/2 ,   size/2 , - size/2 ,   // 顶点5
    - size/2 , - size/2 , - size/2 ,   // 顶点6
      size/2 , - size/2 , - size/2   // 顶点7
  ];

 

 这么大串size/2 (还记得size 是这个函数传进来的参数吧)都是每个顶点的坐标。


 

现在就有个问题,就是size 要多大才能够在屏幕当中显示一个身材合适的立方体呢,一开始可能会假象这个size 取成10 吧,10 这个数字不大不小蛮好的,可是你会发现屏幕中除了一片灰色就什么都没有了,因为10 对于现在这个情况来说太大了,连摄像机都在那个立方体里面了(当然你不会看到立方体里面那个面,因为一般情况下o3d 是只绘制一个面的)。

 

在o3d 中到底用的是什么单位导致10 这个在现实世界中这么小的一个数字能够变得这么大。咳咳,其实貌似在这样的3D 世界中是没有单位之说的( 至少我还不知道) ,因为这个物体的大小是相对于镜头的位置的,如果镜头离物体近了,size 就算很小那物体还是庞然大物,因此想要调整物体大小,可以调整size 的大小,也可以把镜头拉远(还记得第一章中设置照相机的位置吧),也可以把物体放远点,这些都是跟真实世界中差不多的。

 

O3d 中这些顶点的坐标数据都会先放入一个缓冲(buffers) 中,缓冲是一个对象(object) ,用来暂时存储顶点的各种数据的,像位置坐标,色彩,贴图坐标。下面就先创建缓冲对象:

    // 创建顶点缓冲来存放顶点数据

    var positionsBuffer = g_pack.createObject('VertexBuffer');

   


    然后创建域(field) ,域是一小块buffer 。

var positionsField = positionsBuffer . createField ( 'FloatField' , 3 );

positionsBuffer . set ( positionArray );     // 设置该缓冲所用的数组

   


    在一开始有说到要创建一个streamBank 对象,但是没解释为什么,这个真是说来话长啊,因为o3d 使用的是可编程渲染管线(即shader ),而非固定函数渲染管线(fixed function pipeline) ,shader (又可以叫做着色器)可以分为vertex shader (顶点着色器)和pixel shader (像素着色器)。有了pixel shader 之后,你能决定里面的每个像素该怎么填充。因为GPU 浮点运算能力比CPU 比显卡强得多,所以把这些计算交给显卡做再合适不过了。

   


    而这次绘制正方形,每个顶点也是通过vertex shader 来渲染的,是不是有点大材小用了,其实我也这么觉得,不过我们还好暂时还不用去写这个shader ,可以让Js 的接口来做这件事,这就是这个streamBank 做的了,streamBank 把域(field) 中的数据作为vertex shader 的输入。可以是js 中的顶点数据和底层shader 交互的通道。下面这段代码就是把域(field )和vertex shader 联系起来。

streamBank . setVertexStream (
    g_o3d . Stream . POSITION , // semantic: This stream stores vertex positions
    0 ,                     // semantic index: First (and only) position stream
    positionsField ,         // field: the field this stream uses.
    0 );                     // start_index: How many elements to skip in the field.

( 注释的英文我就不翻译了- -)

   

          最后你可以创建一个index buffer ,这个是用来重排上面vertex buffer 里各个点的得排列顺序的。你可以不加这个,不过显示出来的是不是一个立方体就不一定了- - 。

var indicesArray = [
      0 , 1 , 2 ,   // face 1
      2 , 1 , 3 ,
      2 , 3 , 4 ,   // face 2
      4 , 3 , 5 ,
      4 , 5 , 6 ,   // face 3
      6 , 5 , 7 ,
      6 , 7 , 0 ,   // face 4
      0 , 7 , 1 ,
      1 , 7 , 3 ,   // face 5
      3 , 7 , 5 ,
      6 , 0 , 4 ,   // face 6
      4 , 0 , 2
    ];

var indexBuffer = g_pack . createObject ( 'IndexBuffer' );

indexBuffer . set ( indicesArray );

cubePrimitive . indexBuffer = indexBuffer ;

Index buffer 和vertex buffer 的创建方式和存储方式是差不多的. 是两个差不多的Object ,都是用来存储数据的。

   

          函数最后就是把这个已经创建好的shape 返回了。

return cubeShape;

}

        

          函数写好了,得试下先,把前一章用来创建一个球体的语句替换掉,替换成

     var shape = createSphere(material,0.5);

   

         然后就会发现原来那个很丑的圆就变成一个很丑的正方形了。为了体现出这是一个立方体,而不是一个正方形,我们还要做透视投影变换,这个得用shading language 了,

<textarea id="effect">

  // World View Projection matrix that will transform the input vertices

  // to screen space.

  float4x4 worldViewProjection : WorldViewProjection;

 

  // input parameters for our vertex shader

  struct VertexShaderInput {

    float4 position : POSITION;

  };

 

  // input parameters for our pixel shader

  struct PixelShaderInput {

    float4 position : POSITION;

  };

 

  /**

   * The vertex shader simply transforms the input vertices to screen space.

   */

  PixelShaderInput vertexShaderFunction(VertexShaderInput input) {

    PixelShaderInput output;

 

    // Multiply the vertex positions by the worldViewProjection matrix to

    // transform them to screen space.

    output.position = mul(input.position, worldViewProjection);

    return output;

  }

 

  /**

   * This pixel shader just returns the color red.

   */

  float4 pixelShaderFunction(PixelShaderInput input): COLOR {

    return float4(1, 0, 0, 1);  // Red.

  }

 

  // Here we tell our effect file *which* functions are

  // our vertex and pixel shaders.

 

  // #o3d VertexShaderEntryPoint vertexShaderFunction

  // #o3d PixelShaderEntryPoint pixelShaderFunction

  // #o3d MatrixLoadOrder RowMajor

</textarea>

 

    这段就是传说中的shading language 了,它放在一个textarea 中方便js 程序获取,可以看到里面有两个函数,vertexShaderFunction 和pixelShaderFunction ,分别为vertex shader 和pixelShader 的入口,这个怎么工作,通俗点讲vertexShaderFunction 就是把每个顶点的数据传进去后经过一定的计算后再返回进行绘制。pixelShaderFunction 也是这样,在这里vertexShaderFunction 就一句代码:output.position = mul(input.position, worldViewProjection); 就是进行了一次投影变换. 而pixelShaderFunction 就只是把这个立方体的颜色变成了红色。最后三句不要看成是注释啊,前两句申明了vertex shader 的入口函数是vertexShaderFunction , pixelShader 的入口函数是pixelShaderFunction 。

   

          接下来就要用这段代码作为这个立方体的shader string, 还记不记得在前面创建过一个effect ,在后面加上下面两句代码

    var shaderString = document.getElementById('effect').value;

  effect.loadFromFXString(shaderString);               // 载入shading language

   

            然后就会发现这个图形已经有立方体的样子了(我们已经能够脑淫出这是个立方体了- - )。

   

            换个视角看会更像,还记得上一章设置摄像机参数吧。

viewInfo.drawContext.view = g_math.matrix4.lookAt([2, 2, 2],  // 将摄像机放置在(2 ,2 ,2 )位置上

                                                      [0, 0, 0],  // 目标仍是立方体

                                                      [0, 1, 0]); // up

   

 

以下代码中line = plt.contour(Xnew, Ynew, znew, levels=levels, cmap='jet', linewidths=1.5)报错: TypeError: Shapes of x (300, 200) and z (200, 300) do not match 原代码: def ekma_plot(ekma_dir, ekma_cases, run_id, current_date): yyyy, mm, dd = get_ymd(current_date) df_combine = pd.DataFrame([]) max_o3 = [] matched_o3 = [] c_nox = [] ekma_cases_new = ekma_cases.copy() for f in ekma_cases: f_path = ekma_dir + '/' + f + '/output/speciesConcentrations.output' if not os.path.exists(f_path): ekma_cases_new.remove(f) continue if os.stat(f_path).st_size < 1.0: ekma_cases_new.remove(f) continue df_i = pd.read_csv(f_path, sep='\\s+') df_i = df_i.iloc[:24] df_i['t'] = df_i['t'].apply(lambda x: datetime(yyyy, mm, dd) + timedelta(hours=round(x / 3600 - 16))) df_i['hour'] = df_i['t'].apply(lambda x: x.hour) df_i = df_i.groupby('hour').mean() max_o3 = np.append(max_o3, df_i['O3'].max()) matched_o3 = np.append(matched_o3, df_i['O3'].iloc[12]) c_nox = np.append(c_nox, np.nanmean(df_i['NO2'].iloc[9:17])) # +np.nanmean(df_i[5].iloc[9:17]))) $ NO2+NO r_nox = pd.DataFrame([i.split('_')[0] for i in ekma_cases_new]) r_nox = r_nox[0].str.extract(r'(\d+)') # 从NOX20%中提取数字,下面的VOCs同理 r_avocs = pd.DataFrame([i.split('_')[1] for i in ekma_cases_new]) r_avocs = r_avocs[0].str.extract(r'(\d+)') max_o3 = list([i / 2.5e10 for i in max_o3]) matched_o3 = list([i / 2.5e10 for i in matched_o3]) r_nox_ = np.unique(r_nox) r_avocs_ = np.unique(r_avocs) # # # rr_avocs = np.sort([int(i) for i in r_avocs_]) # rr_nox = np.sort([int(i) for i in r_nox_]) # max_o3_ = np.zeros((len(r_avocs_), len(r_nox_))) # for jj in range(0, len(r_nox_)): # for j in range(0, len(r_avocs_)): # a = [max_o3[i] for i in range(0, len(max_o3)) if # (r_avocs.loc[i, 0] == r_avocs_[j]) & (r_nox.loc[i, 0] == r_nox_[jj])] # if len(a) == 1: # max_o3_[j, jj] = a[0] # else: # max_o3_[j, jj] = float('nan') # # fekma = interpolate.interp2d([int(i) for i in r_avocs[0]], [int(i) for i in r_nox[0]], np.reshape(max_o3, (-1, 1)), # kind='cubic') # # xnew = np.arange(0, 200, 1) # ynew = np.arange(0, 300, 1) # znew = fekma(xnew, ynew) ############################# # 构建网格数据(确保唯一值已排序) rr_avocs = np.sort([int(i) for i in r_avocs_]) # x坐标 (AVOCS百分比) rr_nox = np.sort([int(i) for i in r_nox_]) # y坐标 (NOX百分比) # 初始化网格值矩阵(使用NaN占位) max_o3_grid = np.full((len(rr_avocs), len(rr_nox)), np.nan) # 填充网格值(按排序后的坐标索引) for idx, case in enumerate(ekma_cases_new): avoc_val = int(r_avocs.loc[idx, 0]) nox_val = int(r_nox.loc[idx, 0]) i = np.where(rr_avocs == avoc_val)[0][0] # x坐标索引 j = np.where(rr_nox == nox_val)[0][0] # y坐标索引 max_o3_grid[i, j] = max_o3[idx] # 填充值 # 处理NaN值(用0填充,可根据需求改为插值) max_o3_grid = np.nan_to_num(max_o3_grid, nan=0.0) # 创建插值器(使用三次样条) interp = RegularGridInterpolator( (rr_avocs, rr_nox), # 网格坐标 (x, y) max_o3_grid, # 网格值 (二维数组) method='cubic', # 三次插值 bounds_error=False, # 允许外插 fill_value=0.0 # 外插值设为0 ) # 生成新网格 xnew = np.arange(0, 200, 1) ynew = np.arange(0, 300, 1) X, Y = np.meshgrid(xnew, ynew, indexing='ij') # 创建坐标网格 points = np.stack([X.ravel(), Y.ravel()], axis=1) # 换为(N,2)格式 # 执行插值 znew = interp(points).reshape(len(xnew), len(ynew)) ############################# # 实测数据 datafile_path = 'D:\\work\\环科院气象\\unzip\\工程化\\project2\\atchem2\\data\\6.csv' df_real = pd.read_csv(datafile_path) # 将从第 8 列开始的所有列换为数值类型 df_real[df_real.columns[7:]] = df_real[df_real.columns[7:]].apply(pd.to_numeric, errors='coerce') lim_nan = df_real[df_real.columns[7:]].apply(lambda x: x < 0) df_real[lim_nan] = np.nan df_real['Date'] = pd.to_datetime(df_real['开始时间'], format="%Y-%m-%d %H:%M:%S") # df_real['Date'] = pd.to_datetime(df_real['开始时间'], format="%Y/%m/%d %H:%M") df_real['YMD'] = df_real['Date'].apply(lambda x: x.strftime("%Y-%m")) df_real = df_real.fillna(method='bfill') # df_origin = df_origin.interpolate(method='linear') #df_real[df_real.columns[7:]] = df_real[df_real.columns[7:]].interpolate(method='linear') df_real[df_real.columns[7:]] = df_real[df_real.columns[7:]].interpolate(method='bfill') df_real[df_real.columns[7:]] = df_real[df_real.columns[7:]].interpolate(method='ffill') all_map_file = 'D:\\work\\环科院气象\\unzip\\工程化\\project2\\atchem2\\static\\variablesMapping.txt' df_map = pd.read_csv(all_map_file, sep="\\s+") df_cst = pd.read_csv(CSTFILE, header=None) lim_cst1 = df_map['MCM'].apply(lambda x: x in df_cst[0].values) lim_cst2 = df_map['MCM'].apply(lambda x: x not in ['NO', 'NO2', 'NOXINJ', 'CH4', 'PAN', 'CO']) df_var_voc = df_map['VARIABLE'][lim_cst1 & lim_cst2] df_var_nox = df_map['VARIABLE'][df_map['MCM'].apply(lambda x: x in ['NO', 'NO2'])] datelist = df_real['YMD'].unique() all_day_voc = 0 all_day_nox = 0 for i in datelist: i = datelist[0] df_day = df_real[df_real['YMD'] == i] df_day_voc = df_day[df_var_voc.values] df_day_nox = df_day[df_var_nox.values] voc = df_day_voc.iloc[9:17].mean().sum() nox = df_day_nox.iloc[9:17].mean().sum() all_day_voc = all_day_voc + voc all_day_nox = all_day_nox + nox c_vocs = np.asarray([int(i) * all_day_voc / len(datelist) / 100 for i in r_avocs[0]]) c_nox_real = np.asarray([int(i) * all_day_nox / len(datelist) / 100 for i in r_nox[0]]) xnew = np.arange(0, 200, 1) * max(c_vocs) / 200 ynew = np.arange(0, 300, 1) * max(c_nox_real) / 300 Xnew, Ynew = np.meshgrid(xnew, ynew) fig = plt.figure(figsize=(12, 10)) ax1 = fig.add_subplot(111) levels = np.arange(0, znew.max() + 10, 10) line = plt.contour(Xnew, Ynew, znew, levels=levels, cmap='jet', linewidths=1.5)
09-24
import numpy as np import open3d as o3d import matplotlib.pyplot as plt import os from matplotlib.gridspec import GridSpec from tornado.gen import Return # 设置中文字体为宋体 plt.rcParams["font.family"] = ["SimSun"] plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题 class DirectionalCompressionPointCloudProcessor: def __init__(self, random_seed=42, output_dir="output_models", compression_factor=0.5, compression_direction=np.array([0, 0, 1])): """初始化点云处理器,支持沿任意方向向量压缩""" self.random_seed = random_seed self.output_dir = output_dir # 输出目录 self.compression_factor = compression_factor # 压缩因子,1.0表示不压缩,值越小压缩程度越高 self.compression_direction = self._normalize_vector(compression_direction) # 压缩方向向量(已归一化) np.random.seed(random_seed) # 创建输出目录(如果不存在) os.makedirs(self.output_dir, exist_ok=True) # 初始化数据结构 self.data = None self.points = None self.labels = None self.unique_labels = {} # 可视化几何体存储 self.geometries = [] self.geometries_pcd = [] # 处理参数 self.cluster_index = 0 self.min_points_threshold = 20 self.small_cluster_threshold = 10 self.outlier_std_ratio = 1.2 # 平滑处理参数 self.smoothing_iterations = 5 # 平滑迭代次数 self.smoothing_lambda = 0.5 # 平滑强度参数 # 法线计算参数 - 使用固定KNN参数 self.knn_param = 20 # 固定KNN参数 self.consistent_knn = 10 # 法线一致性检查的邻点数 # 可视化参数 self.visualization_pause_time = 0 # 可视化暂停时间,0表示无限等待 self.coord_offset_ratio = 0.3 # 坐标系偏移比例 self.coord_position_choice = 0 # 0:左下, 1:右下, 2:左上, 3:右上 # 三视图参数 self.viewpoint_size = 100 # 三视图点大小 self.save_views = True # 是否保存三视图 def _normalize_vector(self, vector): """归一化向量""" norm = np.linalg.norm(vector) if norm == 0: return np.array([0, 0, 1]) # 默认Z轴方向 return vector / norm def set_compression_direction(self, direction): """设置压缩方向向量""" self.compression_direction = self._normalize_vector(direction) print(f"已设置压缩方向向量: {self.compression_direction}") def load_data(self, file_path): """加载点云数据""" try: self.data = np.genfromtxt( file_path, delimiter=',', skip_header=1 ) self.points = self.data[:, :3] self.labels = self.data[:, 3].astype(int) for label in np.unique(self.labels): self.unique_labels[label] = self.points[self.labels == label] # 计算点云整体边界框,用于放置坐标系 self._calculate_overall_bounding_box() print(f"成功加载数据,共 {len(self.points)} 个点,{len(self.unique_labels)} 个类簇") return True except Exception as e: print(f"数据加载失败: {str(e)}") return False def _calculate_overall_bounding_box(self): """计算整个点云的边界框,用于确定坐标系位置""" if self.points is not None and len(self.points) > 0: self.overall_min = np.min(self.points, axis=0) self.overall_max = np.max(self.points, axis=0) self.overall_center = np.mean(self.points, axis=0) self.overall_size = self.overall_max - self.overall_min else: self.overall_min = np.array([0, 0, 0]) self.overall_max = np.array([100, 100, 100]) self.overall_center = np.array([50, 50, 50]) self.overall_size = np.array([100, 100, 100]) def _filter_outliers(self, pcd): """过滤离群点""" cl, ind = pcd.remove_statistical_outlier( nb_neighbors=15, std_ratio=self.outlier_std_ratio ) return cl, ind def _orient_normals_externally(self, pcd): """将法线方向统一向外""" points = np.asarray(pcd.points) normals = np.asarray(pcd.normals) # 使用包围盒中心方法 bbox = pcd.get_axis_aligned_bounding_box() center = bbox.get_center() to_center = points - center dot_products = np.sum(normals * to_center, axis=1) flip_mask = dot_products > 0 # 指向中心的需要翻 normals[flip_mask] *= -1 pcd.normals = o3d.utility.Vector3dVector(normals) return pcd def _create_colored_coordinate_frame(self, size=1.0): """创建自定义颜色的坐标系:x轴蓝色,y轴绿色,z轴红色""" frame = o3d.geometry.TriangleMesh.create_coordinate_frame(size=size) # 设置坐标轴颜色 frame.vertex_colors[0] = [0, 0, 1] # x轴顶点1 - 蓝色 frame.vertex_colors[1] = [0, 0, 1] # x轴顶点2 - 蓝色 frame.vertex_colors[2] = [0, 1, 0] # y轴顶点1 - 绿色 frame.vertex_colors[3] = [0, 1, 0] # y轴顶点2 - 绿色 frame.vertex_colors[4] = [1, 0, 0] # z轴顶点1 - 红色 frame.vertex_colors[5] = [1, 0, 0] # z轴顶点2 - 红色 # 根据选择计算不同的坐标系位置 offset = self.overall_size * self.coord_offset_ratio if self.coord_position_choice == 0: # 左下 coord_position = self.overall_min - offset elif self.coord_position_choice == 1: # 右下 coord_position = np.array([self.overall_max[0] + offset[0], self.overall_min[1] - offset[1], self.overall_min[2] - offset[2]]) elif self.coord_position_choice == 2: # 左上 coord_position = np.array([self.overall_min[0] - offset[0], self.overall_max[1] + offset[1], self.overall_min[2] - offset[2]]) else: # 右上 coord_position = self.overall_max + offset # 平移坐标系到计算位置 frame.translate(coord_position) return frame def _laplacian_smoothing(self, mesh, iterations=5, lambda_factor=0.5): """对网格进行拉普拉斯平滑处理,减少尖锐边缘和凹陷""" if not mesh: return mesh # 兼容旧版本Open3D,不使用clone()方法 # 创建网格的深拷贝以避免修改原始数据 smoothed_mesh = o3d.geometry.TriangleMesh() smoothed_mesh.vertices = o3d.utility.Vector3dVector(np.asarray(mesh.vertices)) smoothed_mesh.triangles = o3d.utility.Vector3iVector(np.asarray(mesh.triangles)) smoothed_mesh.vertex_normals = o3d.utility.Vector3dVector(np.asarray(mesh.vertex_normals)) smoothed_mesh.triangle_normals = o3d.utility.Vector3dVector(np.asarray(mesh.triangle_normals)) # 获取顶点和邻接信息 vertices = np.asarray(smoothed_mesh.vertices) adjacency_list = smoothed_mesh.adjacency_list # 执行多轮平滑 for _ in range(iterations): new_vertices = vertices.copy() for i in range(len(vertices)): # 确保 i 在 adjacency_list 的有效范围内 if i < len(adjacency_list) and adjacency_list[i]: neighbors = adjacency_list[i] # 计算相邻顶点的平均值 neighbor_avg = np.mean(vertices[neighbors], axis=0) # 应用拉普拉斯平滑公式 new_vertices[i] = vertices[i] + lambda_factor * (neighbor_avg - vertices[i]) # 更新顶点位置 smoothed_mesh.vertices = o3d.utility.Vector3dVector(new_vertices) # 重新计算法线 smoothed_mesh.compute_vertex_normals() return smoothed_mesh def _directional_compression(self, mesh, compression_factor=None, use_max_edge=True): """ 基于单一边缘平面的方向压缩算法: 将方向向量作为法向量,以最边缘的一个顶点确定投影平面, 所有顶点都将相对于这个单一平面进行比例移动 参数: compression_factor: 压缩因子,0-1之间的值表示压缩比例 use_max_edge: 如果为True,使用正方向最边缘的平面;如果为False,使用负方向最边缘的平面 """ if not mesh: return mesh # 如果未指定压缩因子,使用类中定义的默认值 if compression_factor is None: compression_factor = self.compression_factor # 确保压缩因子在有效范围内 if compression_factor <= 0 or compression_factor > 1: raise ValueError("压缩因子必须是0到1之间的值(不包括0,包括1)") # 获取网格顶点 vertices = np.asarray(mesh.vertices) # 确保压缩方向向量已归一化 plane_normal = self.compression_direction / np.linalg.norm(self.compression_direction) # 计算所有顶点在压缩方向上的投影 projections = np.dot(vertices, plane_normal) # 找到投影的最大值和最小值,确定沿压缩方向的边缘位置 min_proj = np.min(projections) max_proj = np.max(projections) original_thickness = max_proj - min_proj if original_thickness <= 0: # 避免厚度为零的情况 return mesh # 选择单一边缘平面作为投影基准 if use_max_edge: # 使用正方向最边缘的平面 edge_proj = max_proj # 平面方程:plane_normal · x + d = 0,其中d = -edge_proj d = -edge_proj # 计算每个顶点到边缘平面的带符号距离(负值表示在平面内侧) distances = np.dot(vertices, plane_normal) + d else: # 使用负方向最边缘的平面 edge_proj = min_proj d = -edge_proj distances = np.dot(vertices, plane_normal) + d # 处理顶点位置 new_vertices = [] for i, vertex in enumerate(vertices): dist = distances[i] # 应用压缩因子:新距离 = 原始距离 × 压缩因子 new_dist = dist * compression_factor # 计算需要移动的距离 move_distance = new_dist - dist # 计算移动向量 move_vector = move_distance * plane_normal # 计算新顶点位置 new_vertex = vertex + move_vector new_vertices.append(new_vertex) # 更新顶点 mesh.vertices = o3d.utility.Vector3dVector(new_vertices) # 平滑处理 mesh = self._laplacian_smoothing( mesh, iterations=self.smoothing_iterations, lambda_factor=self.smoothing_lambda ) # 重新计算法线 mesh.compute_vertex_normals() return mesh def _save_3d_model(self, mesh, cluster_index, label): """保存三维模型为多种格式""" if not mesh: return # 构建文件名(包含类簇信息) base_filename = f"Ore_cluster_{label}_sub_{cluster_index}" # 保存为PLY格式(包含颜色信息) ply_path = os.path.join(self.output_dir, f"{base_filename}_DirCompress.ply") o3d.io.write_triangle_mesh(ply_path, mesh) print(f"已保存PLY模型: {ply_path}") # 保存为STL格式(广泛用于3D打印) stl_path = os.path.join(self.output_dir, f"{base_filename}_DirCompress.stl") o3d.io.write_triangle_mesh(stl_path, mesh) print(f"已保存STL模型: {stl_path}") # 保存为OBJ格式(在许多3D软件中兼容) obj_path = os.path.join(self.output_dir, f"{base_filename}_DirCompress.obj") o3d.io.write_triangle_mesh(obj_path, mesh) print(f"已保存OBJ模型: {obj_path}") def _process_sub_cluster(self, sub_cluster, cluster_index, label): """处理单个子簇并确保生成的网格为闭合体模型""" # 添加微小噪声打破共面性 noise = np.random.normal(0, 1e-6, sub_cluster.shape) sub_cluster += noise # 创建点云对象 pcd = o3d.geometry.PointCloud() pcd.points = o3d.utility.Vector3dVector(sub_cluster) # 过滤点数太少的子簇 if len(pcd.points) < self.small_cluster_threshold: print(f"子簇点数不足 {self.small_cluster_threshold},跳过处理") return None, None # 不过滤离群点 filtered_pcd = pcd print(f"使用KNN参数计算法向量: {self.knn_param}个最近邻") # 计算法向量 - 使用固定KNN参数 filtered_pcd.estimate_normals( o3d.geometry.KDTreeSearchParamKNN(knn=self.knn_param), fast_normal_computation=False # 禁用快速计算,提高精度 ) # 法线方向一致化 filtered_pcd.orient_normals_consistent_tangent_plane(k=self.consistent_knn) # 设置颜色 color = plt.get_cmap("Accent")(cluster_index % 8)[:3] filtered_pcd.paint_uniform_color(color) # 泊松曲面重建 # try: with o3d.utility.VerbosityContextManager( o3d.utility.VerbosityLevel.Debug ) as cm: mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson( filtered_pcd, depth=5, width= 0, scale=1.1, linear_fit=True ) # 后处理:确保网格为闭合流形(兼容性更强) if mesh is not None: # 清理网格 mesh.remove_degenerate_triangles() mesh.remove_duplicated_triangles() mesh.remove_duplicated_vertices() mesh.remove_unreferenced_vertices() # 沿指定方向压缩模型 mesh = self._directional_compression(mesh) # mesh = self.compress_mesh(mesh) # 计算法线(确保可视化正常) mesh.compute_vertex_normals() mesh.paint_uniform_color(color) # 保存三维模型 self._save_3d_model(mesh, cluster_index, label) return filtered_pcd, mesh # except Exception as e: # print(f"泊松重建失败: {str(e)}") # return filtered_pcd, None def process_clusters(self): """处理所有类簇""" if not self.unique_labels: print("没有可处理的类簇数据") return # 遍历每个类簇 for label, cluster_points in self.unique_labels.items(): if label == -1: # 跳过噪声点 continue if len(cluster_points) < self.min_points_threshold: print(f"类簇 {label} 点数不足 {self.min_points_threshold},跳过处理") continue print(f"处理类簇 {label},包含 {len(cluster_points)} 个点") current_clusters = self._get_sub_clusters(label, cluster_points) for sub_cluster in current_clusters: # 传递label参数用于保存文件 pcd, mesh = self._process_sub_cluster(sub_cluster, self.cluster_index, label) if pcd: self.geometries_pcd.append(pcd) if mesh: self.geometries.append(mesh) self.cluster_index += 1 print(f"共处理 {self.cluster_index} 个有效类簇") def _get_sub_clusters(self, label, cluster_points): """获取子簇,直接返回原始类簇""" return [cluster_points] def visualize_orthographic_views(self): """生成并显示三视图(正视图、侧视图、俯视图)""" if self.points is None or len(self.points) == 0: print("没有点云数据可生成三视图") return # 创建图形和子图 fig = plt.figure(figsize=(15, 10)) gs = GridSpec(2, 2, figure=fig) # 正视图 (X-Y平面) ax1 = fig.add_subplot(gs[0, 0]) # 侧视图 (Y-Z平面) ax2 = fig.add_subplot(gs[0, 1]) # 俯视图 (X-Z平面) ax3 = fig.add_subplot(gs[1, 0]) # 设置标题 ax1.set_title('正视图 (X-Y平面)') ax2.set_title('侧视图 (Y-Z平面)') ax3.set_title('俯视图 (X-Z平面)') # 设置坐标轴标签 ax1.set_xlabel('X') ax1.set_ylabel('Y') ax2.set_xlabel('Y') ax2.set_ylabel('Z') ax3.set_xlabel('X') ax3.set_ylabel('Z') # 确保各视图比例一致 ax1.set_aspect('equal') ax2.set_aspect('equal') ax3.set_aspect('equal') # 获取颜色映射 cmap = plt.get_cmap("Accent") # 按类别绘制点 for i, (label, points) in enumerate(self.unique_labels.items()): if label == -1: # 跳过噪声点 continue color = cmap(i % 8)[:3] # 正视图 (X-Y) ax1.scatter(points[:, 0], points[:, 1], c=[color], s=self.viewpoint_size, alpha=0.6, label=f'类别 {label}') # 侧视图 (Y-Z) ax2.scatter(points[:, 1], points[:, 2], c=[color], s=self.viewpoint_size, alpha=0.6) # 俯视图 (X-Z) ax3.scatter(points[:, 0], points[:, 2], c=[color], s=self.viewpoint_size, alpha=0.6) # 添加图例 ax1.legend(bbox_to_anchor=(1.05, 1), loc='upper left') # 调整布局 plt.tight_layout() # 保存三视图(如果启用) if self.save_views: views_path = os.path.join(self.output_dir, "point_cloud_orthographic_views.png") plt.savefig(views_path, dpi=300, bbox_inches='tight') print(f"已保存三视图: {views_path}") # 显示三视图 plt.show() def visualize_results(self): """可视化处理结果,包括3D视图和三视图""" # 先显示三视图 self.visualize_orthographic_views() if not self.geometries_pcd and not self.geometries: print("没有可可视化的3D结果") return # 创建自定义颜色的坐标系 coord_frame_size = np.max(self.overall_size) * 0.1 colored_coord_frame = self._create_colored_coordinate_frame(size=coord_frame_size) # 可视化点云(显示法线) if self.geometries_pcd: print("显示点云与法线方向 (按Q键关闭窗口)") vis = o3d.visualization.Visualizer() vis.create_window(window_name="点云与法线方向", width=1024, height=768) # 添加自定义坐标系 vis.add_geometry(colored_coord_frame) # 添加点云 for pcd in self.geometries_pcd: vis.add_geometry(pcd) # 设置渲染参数 opt = vis.get_render_option() opt.point_size = 3 opt.line_width = 1 opt.show_coordinate_frame = False # 运行可视化 vis.run() vis.destroy_window() # 可视化点云和重建网格 if self.geometries: print("显示点云和泊松重建网格 (按Q键关闭窗口)") vis = o3d.visualization.Visualizer() vis.create_window(window_name="点云和泊松重建网格", width=1024, height=768) # 添加自定义坐标系 vis.add_geometry(colored_coord_frame) # 添加几何体 for geom in self.geometries + self.geometries_pcd: vis.add_geometry(geom) opt = vis.get_render_option() opt.mesh_show_back_face = True opt.point_size = 2 opt.show_coordinate_frame = False vis.run() vis.destroy_window() def compress_mesh(self, mesh): # 检查网格是否有效 if not mesh.has_vertices(): raise ValueError("STL 文件中没有顶点数据。") # 2. 定义压缩方向 compression_direction = np.array([-0.9848, -0.1736, 0.0]) compression_direction = compression_direction / np.linalg.norm(compression_direction) # 3. 定义压缩比例 scale_factor = 0.04 # 压缩比例 # 4. 获取顶点数据 vertices = np.asarray(mesh.vertices) # 新增:计算顶点中心点,用于后续平移校正 center = np.mean(vertices, axis=0) # 新增:将顶点平移到以原点为中心 vertices_centered = vertices - center # 5. 对顶点进行投影与压缩 dot_products = np.dot(vertices_centered, compression_direction) projections = np.outer(dot_products, compression_direction) orthogonal_components = vertices_centered - projections compressed_projections = projections * scale_factor new_vertices_centered = orthogonal_components + compressed_projections #将顶点平移回原来的位置 new_vertices = new_vertices_centered + center # 6. 更新网格顶点 mesh.vertices = o3d.utility.Vector3dVector(new_vertices) # 恢复法线计算(很重要) mesh.compute_vertex_normals() return mesh if __name__ == "__main__": # 定义压缩方向向量 compression_direction = np.array([-0.9848, -0.1736, 0]) # 短轴方向 # 创建处理器实例,指定输出目录和压缩参数 processor = DirectionalCompressionPointCloudProcessor( output_dir="3d_models_New_output", compression_factor= 0.04, # 压缩因子,0.3表示压缩到原来的30% compression_direction=compression_direction ) # 加载数据(请替换为实际文件路径) file_path = 'midpoints_with_cluster_index_Ore2Part_filtered_10_602030.csv' if processor.load_data(file_path): # 可以在处理前动态更改压缩方向 # processor.set_compression_direction(np.array([0, 1, 0])) # 例如改为Y轴方向 # 处理类簇(处理过程中会自动保存模型) processor.process_clusters() # 可视化结果(现在包括三视图) processor.visualize_results() else: print("无法继续处理,数据加载失败") 如何保证压缩后的网格模型尽量减少三角网交叉现象
09-16
### 使用 Open3D 实现 Alpha Shapes Alpha Shapes 是一种用于表示点集形状的方法,可以看作是对 Delaunay 三角剖分的一种推广。通过调整参数 α 的大小,可以在不同层次上捕捉几何结构的细节。 为了使用 Open3D 库实现 Alpha Shapes,在 Python 中可以通过以下方式完成: ```python import open3d as o3d import numpy as np # 创建一个示例点云数据 points = np.array([ [0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 1] ]) # 将 NumPy 数组换为 Open3D PointCloud 对象 pcd = o3d.geometry.PointCloud() pcd.points = o3d.utility.Vector3dVector(points) # 计算 Alpha Shape alpha = 1.5 # 调整此值来改变 Alpha Shape 的紧密程度 mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(pcd, alpha=alpha) # 可视化结果 o3d.visualization.draw_geometries([mesh]) ``` 上述代码展示了如何创建简单的三维点集合,并将其作为输入传递给 `create_from_point_cloud_alpha_shape` 函数以生成相应的 Alpha Shape 结构[^1]。注意这里使用的 `alpha` 参数决定了最终模型的紧凑度;较小的α值会产生更精细的结果,而较大的α则会使表面更加平滑简化。 对于颜色绑定方面,虽然提供的参考资料提到了关于立方体的颜色处理方法[^2],但这部分并不适用于当前讨论的主题——即基于 Open3D 的 Alpha Shapes 构建过程。如果确实需要对生成的 mesh 进行着色,则可以根据具体需求设置顶点或面片的颜色属性。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值