WO2022104681A1 - 一种柱坐标系等值面提取方法 - Google Patents

一种柱坐标系等值面提取方法 Download PDF

Info

Publication number
WO2022104681A1
WO2022104681A1 PCT/CN2020/130374 CN2020130374W WO2022104681A1 WO 2022104681 A1 WO2022104681 A1 WO 2022104681A1 CN 2020130374 W CN2020130374 W CN 2020130374W WO 2022104681 A1 WO2022104681 A1 WO 2022104681A1
Authority
WO
WIPO (PCT)
Prior art keywords
processing unit
isosurface
data
vertex
coordinate system
Prior art date
Application number
PCT/CN2020/130374
Other languages
English (en)
French (fr)
Inventor
卢光辉
蔡洪斌
王忠杰
张珀熙
王小凤
陈简
Original Assignee
电子科技大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 电子科技大学 filed Critical 电子科技大学
Priority to PCT/CN2020/130374 priority Critical patent/WO2022104681A1/zh
Publication of WO2022104681A1 publication Critical patent/WO2022104681A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T19/00—Manipulating 3D models or images for computer graphics

Definitions

  • the invention belongs to the technical field of three-dimensional data field visualization, and in particular relates to a method for extracting isosurfaces in a cylindrical coordinate system.
  • the Marching Cubes algorithm is a classic isosurface extraction algorithm, which is simple, robust, and fast. However, it is ambiguous and may generate isosurfaces with incorrect topology.
  • the Marching Tetrahedra algorithm cuts the minimum processing unit of Marching Cubes, and cuts the hexahedron into several tetrahedra, trying to solve the ambiguity problem, but it cannot completely correct the ambiguity problem and reduce the extraction speed.
  • the Marching Cubes 33 algorithm solves the ambiguity problem by revising the division template of Marching Cubes, and retains the advantages of robustness and speed, but the computational efficiency still needs to be improved to meet the visualization requirements in the case of large data volumes.
  • the present invention provides a method for extracting an isosurface of a cylindrical coordinate system, comprising the following main steps:
  • Step 1 read in the data, read in the three-dimensional cylindrical regular grid data.
  • Step 2 data preprocessing, determine whether the data is completely processed, if it is completely processed, go to step 5 to perform surface subdivision, otherwise read a processing unit to determine whether it is a valid unit, and if it is a valid unit, perform isosurface extraction , otherwise skip the processing unit and read the next processing unit.
  • This step mainly includes:
  • Step 2.1 obtaining the minimum scalar value S min and the maximum scalar value S max among the eight scalar values of the processing unit.
  • Step 2.2 if the target value S target of the isosurface is not in the range [S min , S max ], then the processing unit is an invalid unit, skip the processing unit, read the next processing unit, and repeat step 2; otherwise, is a valid unit, go to step 3.
  • Step 3 judging the division of a single processing unit, marking the processing unit, and judging the division of the processing unit according to the mark.
  • This step mainly includes:
  • Step 3.1 according to the S target marking each vertex on the processing unit, determine the positional relationship between the vertex and the isosurface. If the vertex scalar value S ⁇ S target , then the vertex is outside the isosurface, marked with "+”, otherwise it is inside the isosurface, marked with "-". Use an 8-bit binary number to store the tag variable of the processing unit, "+” is recorded as 1, "-" is recorded as 0.
  • Step 3.2 establishing a situation look-up table and a subdivision look-up table.
  • the index in the case lookup table contains the case number to which the processing unit belongs and information about the triangle patch vertices.
  • the index in the subdivision lookup table contains the number of triangular patches in the processing unit and the edge number of the vertex of the triangular patch.
  • Step 3.3 according to the marked variable to the situation lookup table, look up the basic subdivision situation category of the processing unit, and obtain the situation index. According to the situation, the index is searched in the subdivision lookup table to find the subdivision index, and the intersection mode of the processing unit and the isosurface is obtained.
  • Step 4 Calculate the coordinates of the vertices of the isosurface.
  • step 3 the intersection method between the processing unit and the isosurface is obtained. Therefore, the coordinates of the intersection between the processing unit and the isosurface are calculated, that is, the vertex coordinates of the isosurface, so as to obtain the equivalent value. Triangular mesh of faces. Return to step 2 after completing this step.
  • This step mainly includes:
  • Step 4.1 use linear interpolation to calculate the cylindrical coordinates of the vertices of the isosurface.
  • Step 4.2 using the conversion formula from cylindrical coordinates to rectangular coordinates to convert the cylindrical coordinates of the vertices into rectangular coordinates.
  • Step 4.3 connect the vertices in the order defined in the subdivision index to get a triangular mesh.
  • step 5 the triangular mesh is refined, and the extracted isosurface is subdivided, so that the edges and corners caused by the sparse data are smoother.
  • the beneficial effects of the invention are as follows: an effective isosurface extraction method is provided for the data uniformly sampled in the cylindrical coordinate system, the drawing efficiency is improved, and the correctness of the isosurface topology structure is guaranteed; for sparse data, the combination of The tessellation technology enhances the visualization effect while reducing the huge resource consumption caused by computing.
  • Fig. 1 shows the basic flow of a cylindrical coordinate system isosurface extraction method of the present invention
  • Fig. 2 shows the concrete flow of a kind of cylindrical coordinate system isosurface extraction method of the present invention
  • Fig. 3 shows the data format of the regular grid data under the cylindrical coordinate system of a cylindrical coordinate system isosurface extraction method of the present invention
  • Fig. 4 shows the subdivision of all processing units in the isosurface extraction algorithm under the cylindrical coordinate system of a cylindrical coordinate system isosurface extraction method of the present invention
  • Fig. 5 shows the sub-flow chart in the isosurface extraction algorithm under the cylindrical coordinate system of a cylindrical coordinate system isosurface extraction method of the present invention - determine the situation;
  • Fig. 6 shows the sub-flow chart in the isosurface extraction algorithm under the cylindrical coordinate system of a cylindrical coordinate system isosurface extraction method of the present invention - surface subdivision;
  • FIG. 7 shows the vertex addition and update strategy of the tessellation algorithm in the isosurface extraction algorithm under the cylindrical coordinate system of a cylindrical coordinate system isosurface extraction method of the present invention.
  • Step 1 Read in the data, and read in the three-dimensional cylindrical regular grid data. As shown in Figure 3, every 8 adjacent points constitute a processing unit, and the data contained in the processing unit includes the scalar values of the 8 vertices and the Cylindrical coordinates.
  • Step 2 data preprocessing, determine whether the data is completely processed, if it is completely processed, go to step 5 to perform surface subdivision, otherwise read a processing unit to determine whether there is an intersection with the isosurface, and if there is an intersection, perform the equivalent Face extraction, otherwise skip this processing unit and read the next processing unit.
  • This step mainly includes:
  • Step 2.1 obtaining the minimum scalar value S min and the maximum scalar value S max among the eight scalar values of the processing unit.
  • Step 2.2 if the target value S target of the isosurface is not in the range [S min , S max ], then the processing unit and the isosurface have an intersection, skip the processing unit, read the next processing unit, and repeat step 2 ; otherwise, it is a valid unit, go to step 3.
  • Step 3 judging the division of a single processing unit, marking the processing unit, and judging the division of the processing unit according to the mark.
  • This step mainly includes:
  • Step 3.1 according to the S target marking each vertex on the processing unit, determine the positional relationship between the vertex and the isosurface. If the vertex scalar value S ⁇ S target , then the vertex is outside the isosurface, marked with "+”, otherwise it is inside the isosurface, marked with "-". Use an 8-bit binary number to store the tag variable of the processing unit, "+” is recorded as 1, "-" is recorded as 0.
  • Step 3.2 establishing a situation look-up table and a subdivision look-up table.
  • Each vertex has 2 possible states, so a processing unit containing 8 vertices has 256 cases.
  • the 256 cases can be reduced to 128, and it is guaranteed that none of the subdivision cases are lost.
  • the case lookup table is an int array with a size of 128.
  • the subscript of the array element corresponds to the tag variable of the processing unit.
  • Each element is composed of 4-digit hexadecimal numbers, including the case number to which the corresponding processing unit belongs and related Information about the vertices of the triangular patch.
  • the 128 cases are reduced to 15 base cases.
  • Marching Cubes 33 expands the 15 basic cases to 33 to ensure the correctness of the extracted isosurface.
  • the subdivision lookup table contains all the possibilities of these 33 subdivision cases. It is an int array of size 2182, each element is composed of 4-bit hexadecimal numbers, representing a triangular patch, including the triangular patch in processing. Location information in the cell.
  • the index in the subdivision lookup table contains the number of triangular facets of the processing unit and the edge number where the vertex is located.
  • Step 3.3 according to the index of the situation where the processing unit is found in the situation lookup table, if the processing unit belongs to situations 1, 2, 5, 8, 9, 11 and 14, then there is no ambiguity, and it can be directly searched in the subdivision Find the corresponding subdivision index in the table, that is, the composition of the triangular patch. As shown in Figure 4 and Figure 5. For other cases, make the following analysis:
  • the processing unit belongs to case 3, then there is a facet ambiguity.
  • the face number that needs further testing can be obtained from the index found in the situation lookup table. Assuming that the four vertices of the face are A, B, C, and D respectively, and their scalar values are S A , S B , S C , and S D , respectively, then, use formula (1) to perform face ambiguity test on this face .
  • the processing unit belongs to case 6, then there is a facet ambiguity.
  • the serial number of the test surface is obtained according to the situation index number. If there is no ambiguity surface, the segmentation index will be returned directly; otherwise, an internal ambiguity will be generated. Similarly, the diagonal line of the test will be obtained according to the situation index number, and the internal ambiguity test will be performed and the segmentation index will be returned. sub-index.
  • the processing unit belongs to case 7, then there are three facet ambiguities.
  • the face ambiguity test is performed on the three faces respectively. If all faces are ambiguous faces, then there will be an internal ambiguity. After the internal ambiguity test, the subdivision index will be returned. Otherwise, the corresponding subdivision index will be returned according to the face ambiguity test.
  • the processing unit belongs to case 10, then there are two face ambiguities in opposite positions, and one of the two face ambiguities must be non-ambiguous. Perform a face ambiguity test on the two faces respectively, if both faces are not ambiguous faces, then an internal ambiguity will be generated, the internal ambiguity test will be performed, and the corresponding subdivision index will be returned according to the result.
  • processing unit belongs to case 12, then there are two facet ambiguities, and there must be one facet that is not ambiguous. Perform a face ambiguity test on the two faces respectively. If there is an ambiguity face, the subdivision index will be returned directly; otherwise, an internal ambiguity will be generated, and an internal ambiguity test will be performed.
  • step 4 the coordinates of the vertices of the isosurface are calculated, and in step 3, the intersection mode of the processing unit and the isosurface is obtained, thereby calculating the coordinates of the intersection point between the processing unit and the isosurface, that is, the vertex coordinates of the isosurface. Connect vertices in the order defined in the mesh index, resulting in a triangular mesh. Return to step 2 after completing this step.
  • This step mainly includes:
  • Step 4.1 use linear interpolation to calculate the cylindrical coordinates of the vertices of the isosurface. Assuming that the scalar values of the two endpoints of the side where the intersection is located are S 0 , S 1 , and the cylindrical coordinates are v 0 , v 1 , then the cylindrical coordinates of the points on the isosurface are
  • Step 4.2 using formula (5) to convert the vertex cylindrical coordinates into rectangular coordinates.
  • Step 4.3 connect the vertices in the order defined in the subdivision index to get a triangular mesh.
  • step 5 the triangular mesh is refined, and the extracted isosurface is subdivided, so that the edges and corners caused by the sparse data are smoother.
  • the flow chart shown in FIG. 6 shows the operation flow of tessellation in the present invention.
  • This step mainly includes:
  • Step 5.1 input grid data, use half-edge structure to store data.
  • Step 5.2 Obtain edge data and add new vertices to the calculation strategy. Determine whether the vertex to be added is located on the grid boundary, if so, assuming that the vertex coordinates of the grid boundary where the new vertex is located are v 0 , v 1 , according to the adding strategy shown in Figure 7(b), calculate the coordinate v of the newly added vertex as
  • Step 5.3 read the vertex data and update the original vertex. Determine whether the original vertex is located at the grid boundary. If so, suppose the original vertex coordinate is v 0 , and the mesh boundary vertex coordinates of the original vertex are v 1 and v 2 . According to the update strategy shown in Figure 7(d), the original vertex v 0 After the update, the coordinate v is
  • step 5.3 Read the next vertex data, if back to the starting point, go to step 5.4; otherwise, repeat step 5.3
  • Step 5.4 regenerate the triangular mesh.

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Graphics (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Image Generation (AREA)

Abstract

一种柱坐标系等值面提取方法。该方法包括：读入数据，数据预处理，判断单个处理单元的剖分情况，计算等值面顶点坐标，曲面细化三角网格。该方法利用具有33种三角剖分情况的Marching Cubes 33算法，解决了传统等值面提取算法Marching Cubes的二义性问题，使得生成的三角网格具有正确的拓扑结构。通过改变Marching Cubes 33的最小处理单元的几何结构，使其能够应用在柱坐标系下采样的柱形规则网格数据，生成质量较高的等值面。结合曲面细分技术，解决了稀疏数据生成粗糙等值面的问题，生成高质量的柱形等值面。同时，对每个最小处理单元进行预判断，跳过了三维数据场中存在的大量无效体素，缩短了绘制时间。

Description

一种柱坐标系等值面提取方法 技术领域
本发明属于三维数据场可视化技术领域，尤其涉及一种柱坐标系等值面提取方法。
背景技术
随着计算机硬件的不断进步，计算机图形学和图像处理等技术的持续发展，使得对科学数据的实时分析成为了可能。利用可视化技术对大规模数据进行分析及显示，便于快速提取数据中具有典型意义的特征与结果，加强用户对数据的分析和理解，为各应用领域提供了强大的分析工具。等值面提取方法作为可视化技术的一个分支，提供了一种机制来理解三维数据场的层次结构，加快了数据的处理速度，改变了长期以来数据处理的低效性，避免了人工处理过程中信息的大量丢失。
可视化技术目前已经得到了广泛的发展，但是主要的数据研究对象仍局限在规则网格数据。然而，在计算流体力学、计算电磁学等科学仿真实验中，用于模拟的网格不同于可视化方法中的网格，通常为一种结构化网格，生成的数据类型为非规则网格数据。从目前的研究现状来看，对于这类非规则网格数据，常用的方式是对原始数据进行格式转换后使用传统的可视化方法。例如，针对将非规则网格数据转换为规则网格数据这一过程，付则鑫等人在《基于立方体网格插值算法的气象雷达数据的三维重建》一文中针对锥型的气象雷达数据，提出立方体网格插值算法；薛源等人在《基于网格吸附插值法的不规则尾流场构建方法》中针对不规则的尾流场数据，提出网格吸附插值算法。除此之外，杨超在《虚拟战场中电磁环境三维建模与绘制方法研究》中，针对柱形规则网格电磁数据，提出柱坐标系下的体绘制算法。但是，网格转换法在数据预处理上需要消耗大量的时间成本，得到的结果也无法完全保留原始数据的特征，体绘制算法相对于等值面提取算法来说，计算复杂度更高，对计算机的硬件环境有一定的要求。
Marching Cubes算法是一种经典的等值面提取算法，具有简单、健壮、快速等特点，但存在二义性，可能生成拓扑结构不正确的等值面。Marching Tetrahedra 算法对Marching Cubes的最小处理单元进行切割，由六面体切割为若干个四面体，试图解决二义性问题，但并不能完全修正二义性问题，并且降低了提取速度。Marching Cubes 33算法通过订正Marching Cubes的剖分模板，解决了二义性问题，并保留了健壮、快速等优势，但计算效率仍有待提高，以满足大数据量情况下的可视化要求。
发明内容
为了克服现有技术的不足，本发明提供了一种柱坐标系等值面提取方法，包括以下主要步骤：
步骤1，读入数据，读入三维柱形规则网格数据。
步骤2，数据预处理，判断数据是否处理完全，若处理完全，则进入步骤5进行曲面细分，否则读取一个处理单元，判断是否为有效单元，若为有效单元，则进行等值面提取，否则跳过该处理单元，读取下一个处理单元。
本步骤主要包括：
步骤2.1，获取处理单元8个标量值中的最小标量值S min和最大标量值S max。
步骤2.2，若等值面的目标值S target不在范围[S min,S max]内，那么该处理单元为无效单元，跳过该处理单元，读取下一处理单元，重复步骤2；否则，为有效单元，进入步骤3。
步骤3，判断单个处理单元的剖分情况，标记处理单元，根据标记判断该处理单元的剖分情况。
本步骤主要包括：
步骤3.1，根据S target标记处理单元上的每个顶点，判断顶点与等值面的位置关系。若顶点标量值S≥S target，那么该顶点位于等值面外部，标记为“+”，否则位于等值面内部，标记为“-”。使用一个8位二进制数保存处理单元的标记变量，“+”记为1，“-”记为0。
步骤3.2，建立情况查找表和剖分查找表。情况查找表中的索引包含了该处理单元所属的情况号以及有关三角面片顶点的信息。剖分查找表中的索引包含了处理单元中的三角面片个数、三角面片顶点所在的边序号。
步骤3.3，根据标记变量到情况查找表中查找处理单元的基本剖分情况类别， 获取情况索引。根据情况索引到剖分查找表中查找剖分索引，得到处理单元与等值面的相交方式。
步骤4，计算等值面顶点坐标，步骤3中得到处理单元与等值面的相交方式，由此，计算处理单元与等值面的交点坐标，即等值面的顶点坐标，从而得到等值面的三角网格。完成该步骤后返回至步骤2。
本步骤主要包括：
步骤4.1，利用线性插值计算等值面的顶点柱坐标。
步骤4.2，利用柱坐标到直角坐标的转换公式将顶点的柱坐标转换为直角坐标。
步骤4.3，按照剖分索引中定义的顺序连接顶点，得到三角网格。
步骤5，曲面细化三角网格，对提取的等值面进行曲面细分，使得稀疏数据导致的棱角部分更加平滑。
本发明的有益效果是：为柱坐标系下均匀采样的数据提供了一种有效的等值面提取方式，并提高了绘制效率以及保证了等值面拓扑结构的正确性；对于稀疏数据，结合曲面细分技术，在减少计算带来的巨大资源消耗的同时，增强了可视化效果。
附图说明
图1示出了本发明一种柱坐标系等值面提取方法的基本流程；
图2示出了本发明一种柱坐标系等值面提取方法的具体流程；
图3示出了本发明一种柱坐标系等值面提取方法的柱坐标系下的规则网格数据的数据格式；
图4示出了本发明一种柱坐标系等值面提取方法的柱坐标系下等值面提取算法中的所有处理单元剖分情况；
图5示出了本发明一种柱坐标系等值面提取方法的柱坐标系下等值面提取算法中的子流程图——确定所属情况；
图6示出了本发明一种柱坐标系等值面提取方法的柱坐标系下等值面提取算法中的子流程图——曲面细分；
图7示出了本发明一种柱坐标系等值面提取方法的柱坐标系下等值面提取算法中的曲面细分算法的顶点添加和更新策略。
具体实施方式
下面结合附图和实施例对本发明优先实施方式进一步说明。
图1所示的流程图给出了本发明实施的基本流程，图2所示的流程图给出了本发明整个实施的具体过程：
步骤1，读入数据，读入三维柱形规则网格数据，如图3所示，相邻的每8个点构成一个处理单元，处理单元中包含的数据包括8个顶点的标量值和柱坐标。
步骤2，数据预处理，判断数据是否处理完全，若处理完全，则进入步骤5进行曲面细分，否则读取一个处理单元，判断是否与等值面有交点，若有交点，则进行等值面提取，否则跳过该处理单元，读取下一个处理单元。
本步骤主要包括：
步骤2.1，获取处理单元8个标量值中的最小标量值S min和最大标量值S max。
步骤2.2，若等值面的目标值S target不在范围[S min,S max]内，那么该处理单元与等值面存在交点，跳过该处理单元，读取下一处理单元，重复步骤2；否则，为有效单元，进入步骤3。
步骤3，判断单个处理单元的剖分情况，标记处理单元，根据标记判断该处理单元的剖分情况。
本步骤主要包括：
步骤3.1，根据S target标记处理单元上的每个顶点，判断顶点与等值面的位置关系。若顶点标量值S≥S target，那么该顶点位于等值面外部，标记为“+”，否则位于等值面内部，标记为“-”。使用一个8位二进制数保存处理单元的标记变量，“+”记为1，“-”记为0。
步骤3.2，建立情况查找表和剖分查找表。每个顶点有2种可能的状态，那么包含8个顶点的处理单元有256种情况，根据互补对称性，可以将256种情况缩减至128种，并保证了不丢失任何一种剖分情况。情况查找表是一个大小为128的int数组，数组元素的下标与处理单元的标记变量对应，每个元素由4位16进制数组成，包含了与其对应的处理单元所属的情况号以及有关三角面片顶点的信息。根据旋转对称性，将128种情况缩减为15种基本情况。为了解决二 义性问题Marching Cubes 33将15种基本情况扩充为33种，保证提取的等值面的正确性。剖分查找表包含了这33种剖分情况的所有可能，是一个大小为2182的int数组，每个元素由4位16进制数组成，代表一个三角面片，包含了三角面片在处理单元中的位置信息。剖分查找表中的索引包含了处理单元的三角面片个数、顶点所在的边序号。
步骤3.3，根据标记变量到情况查找表中的找到处理单元的情况索引，若处理单元属于情况1、2、5、8、9、11和14，那么它不存在歧义，可以直接在剖分查找表中找到对应的剖分索引，即三角面片构成方式。如图4、图5所示。针对其他情况，作以下分析：
a.若处理单元属于情况3，那么存在一个面歧义。根据在情况查找表中找到的索引可以得到需要进一步测试的面序号。假设该面的4个顶点分别为A、B、C、D，它们的标量值分别为S A、S B、S C、S D，然后，使用公式(1)对该面进行面歧义测试。
Figure PCTCN2020130374-appb-000001
若S≥0，那么该面存在歧义，否则不存在歧义，根据测试结果在剖分查找表中找到对应的剖分索引。
b.若处理单元属于情况4，那么存在一个内部歧义。根据在情况查找表中找到的索引可以确定需要进行内部歧义测试的对角线。假设处理单元8个顶点的标量值为
Figure PCTCN2020130374-appb-000002
其中
Figure PCTCN2020130374-appb-000003
为测试对角线上两顶点的标量值，若F t为A 0到A 1方向上对处理单元的任意一水平切平面，t∈(0,1)，那么，关于t有如下关系式
at 2+bt+c>0        (2)
Figure PCTCN2020130374-appb-000004
当a≥0时，该处理单元不存在内部歧义；当a<0时，若-b/2a<0，不存在歧义，若-b/2a≥0，且A 0到A 1存在歧义面，那么该处理单元存在内部歧义，否则不存在。由此在剖分情况查找表中确定剖分索引。
c.若处理单元属于情况6，那么存在一个面歧义。根据情况索引号得到测试面的序号，若不存在歧义面，则直接返回剖分索引；否则将又产生一个内部歧义，同样根据情况索引号得到测试的对角线，进行内部歧义测试并返回剖分索引。
d.若处理单元属于情况7，那么存在三个面歧义。对三个面分别进行面歧义测试，如果所有面均是歧义面，那么又会存在一个内部歧义，再进行内部歧义测试后返回剖分索引，否则根据面歧义测试情况返回对应的剖分索引。
e.若处理单元属于情况10，那么存在两个位置对立的面歧义，并且两个面歧义中一定存在一个不是歧义面。对两个面分别进行面歧义测试，若两个面均不是歧义面，那么将产生一个内部歧义，进行内部歧义测试，根据结果返回对应的剖分索引。
f.若处理单元属于情况12，那么存在两个面歧义，并且一定存在一个不是歧义面。对两个面分别进行面歧义测试，若存在一个歧义面，直接返回剖分索引；否则将会产生一个内部歧义，进行内部歧义测试。
g.若处理单元属于情况13，那么所有的面均存在面歧义。对所有面进行面歧义测试，若存在四个歧义面，并且非歧义的两个面相邻，此时将产生一个内部歧义，进行内部歧义测试，否则，根据面歧义结果直接返回剖分索引。
步骤4，计算等值面顶点坐标，步骤3中得到处理单元与等值面的相交方式，由此，计算处理单元与等值面的交点坐标，即等值面的顶点坐标。按照剖分索引中定义的顺序连接顶点，得到三角网格。完成该步骤后返回至步骤2。
本步骤主要包括：
步骤4.1，利用线性插值计算等值面的顶点柱坐标。假设交点所在边的两个端点的标量值为S 0、S 1，柱坐标为v 0、v 1，那么等值面上点的柱坐标为
Figure PCTCN2020130374-appb-000005
步骤4.2，利用公式(5)将顶点柱坐标转换为直角坐标。
Figure PCTCN2020130374-appb-000006
其中，
Figure PCTCN2020130374-appb-000007
为等值面顶点的柱坐标，(x,y,z)为等值面顶点的直角坐标，也是最终的绘制坐标。
步骤4.3，按照剖分索引中定义的顺序连接顶点，得到三角网格。
步骤5，曲面细化三角网格，对提取的等值面进行曲面细分，使得稀疏数据导致的棱角部分更加平滑。图6所示的流程图给出了本发明中曲面细分操作流程。
本步骤主要包括：
步骤5.1，输入网格数据，使用半边结构存储数据。
步骤5.2，获取边数据，计算策略新增顶点。判断将要添加顶点是否位于网格边界，若是，假设新增顶点所在网格边界顶点坐标为v 0、v 1，根据图7(b)所示的添加策略，计算新增的顶点坐标v为
Figure PCTCN2020130374-appb-000008
若将要添加的顶点位于网格内部，假设新增顶点的所在边顶点坐标为v 0、v 1，其他周围网格顶点坐标为v 2、v 3，根据图7(a)所示的添加策略，计算新增的顶点坐标v为
Figure PCTCN2020130374-appb-000009
读取下一个边数据，如果回到起始边，则进入步骤5.3；否则，重复步骤5.2。
步骤5.3，读取顶点数据，更新原始顶点。判断原始顶点是否位于网格边界，若是，假设原始顶点坐标为v 0，原始顶点所在网格边界顶点坐标为v 1、v 2，根据图7(d)所示的更新策略，原始顶点v 0经过更新后坐标v为
Figure PCTCN2020130374-appb-000010
若原始顶点位于网格内部，假设原始顶点坐标为v 0，周围网格顶点坐标为v i(i＝1,2,...,n)，周围网格顶点数n＝6，根据图7(c)所示的更新策略，原始顶点v 0经过更新后坐标v为
Figure PCTCN2020130374-appb-000011
读取下一个顶点数据，如果回到起点，则进入步骤5.4；否则，重复步骤5.3
步骤5.4，重新生成三角网格。

Claims (5)

  1. 一种柱坐标系等值面提取方法，其特征在于，包括以下步骤：,
    步骤1，读入数据，读入三维柱形规则网格数据。
    步骤2，数据预处理，判断数据是否全部处理完成，以及处理单元是否与等值面有交点，跳过三维数据场中的与等值面不存在交点的单元。
    步骤3，判断单个处理单元的剖分情况，标记处理单元，利用Marching Cubes 33算法中提出的渐进判定法判断该处理单元的剖分情况，使得生成的等值面具有正确的拓扑结构。
    步骤4，计算等值面的顶点坐标，由此生成直角坐标系下三维数据的等值面。
    步骤5，曲面细化三角网格，对提取的等值面进行曲面细分。
  2. 根据权利要求1中所述的一种柱坐标系等值面提取方法，其特征在于所述的步骤1中，柱坐标系下的数据以柱形规则网格的形式存在，相邻的8个点作为一个处理单元，每个处理单元分开处理，最终的等值面由所有处理单元的三角面片拼接得到。
  3. 根据权利要求1中所述的一种柱坐标系等值面提取方法，其特征在于所述的步骤2中，数据预处理。所述的步骤2进一步包括：
    步骤2.1，获取处理单元8个标量值中的最小标量值S min和最大标量值S max。
    步骤2.2，若等值面的目标值S target不在范围[S min,S max]内，那么该处理单元与等值面没有交点，跳过该处理单元，读取下一处理单元，重复步骤2；否则，为有效单元，进入步骤3。
  4. 根据权利要求1中所属的一种柱坐标系等值面提取方法，其特征在于所述的步骤3中，判断处理单元的剖分情况。所述的步骤3进一步包括：
    步骤3.1，根据S target标记处理单元上的每个顶点，判断顶点与等值面的位置关系。若顶点标量值S≥S target，那么该顶点位于等值面外部，标记为1，否则位于等值面内部，标记为0。使用一个8位二进制数保存处理单元的标记变量。
    步骤3.2，建立情况查找表和剖分查找表。情况查找表中的索引包含了该处理单元所属的情况索引以及有关三角面片顶点的信息，是一个大小为128的int数组，数组下标与标记变量对应，每个元素由4位16进制数组成。剖分查找表中的索引包含了处理单元与等值面的相交信息，包括处理单元中的三角面片个数、 三角面片顶点所在的边序号，是一个大小为2182的int数组，每个元素由4位16进制数组成，代表一个三角面片，包含了三角面片顶点在处理单元中的位置信息。
    步骤3.3，根据标记变量到情况查找表中查找处理单元的基本剖分情况，获取情况索引。根据情况索引，利用渐进判定法判断处理单元是否存在面歧义或内部歧义，并结合剖分查找表，解决传统Marching Cubes方法的二义性问题，得到剖分索引，得到处理单元与等值面的相交方式。
  5. 根据权利要求1中所属的一种柱坐标系等值面提取方法，其特征在于所述的步骤4中，计算等值面顶点坐标。所述的步骤4进一步包括：
    步骤4.1，根据处理单元与等值面的交点所在边的端点的标量值、柱坐标，利用线性插值计算等值面的顶点柱坐标。
    步骤4.2，利用柱坐标到直角坐标的转换公式，将顶点柱坐标转换为直角坐标。
    步骤4.3，按照剖分索引中定义的顺序连接顶点，得到三角网格。完成该步骤后返回至步骤2.
PCT/CN2020/130374 2020-11-20 2020-11-20 一种柱坐标系等值面提取方法 WO2022104681A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
PCT/CN2020/130374 WO2022104681A1 (zh) 2020-11-20 2020-11-20 一种柱坐标系等值面提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2020/130374 WO2022104681A1 (zh) 2020-11-20 2020-11-20 一种柱坐标系等值面提取方法

Publications (1)

Publication Number Publication Date
WO2022104681A1 true WO2022104681A1 (zh) 2022-05-27

Family

ID=81708224

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2020/130374 WO2022104681A1 (zh) 2020-11-20 2020-11-20 一种柱坐标系等值面提取方法

Country Status (1)

Country Link
WO (1) WO2022104681A1 (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115270614A (zh) * 2022-07-18 2022-11-01 郑州轻工业大学 一种泥水循环系统多物理场数字孪生体可视化生成方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110007933A1 (en) * 2009-07-10 2011-01-13 Microsoft Corporation 3D Image Processing
CN102509358A (zh) * 2011-11-24 2012-06-20 武汉大学 一种基于三棱柱剖分法的动态三维等值面提取方法
CN103295269A (zh) * 2013-06-26 2013-09-11 电子科技大学 一种电磁环境体数据等值面提取方法
CN109940894A (zh) * 2019-04-01 2019-06-28 上海大学 一种基于有限支撑半径控制的卷积曲面混合建模方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110007933A1 (en) * 2009-07-10 2011-01-13 Microsoft Corporation 3D Image Processing
CN102509358A (zh) * 2011-11-24 2012-06-20 武汉大学 一种基于三棱柱剖分法的动态三维等值面提取方法
CN103295269A (zh) * 2013-06-26 2013-09-11 电子科技大学 一种电磁环境体数据等值面提取方法
CN109940894A (zh) * 2019-04-01 2019-06-28 上海大学 一种基于有限支撑半径控制的卷积曲面混合建模方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115270614A (zh) * 2022-07-18 2022-11-01 郑州轻工业大学 一种泥水循环系统多物理场数字孪生体可视化生成方法
CN115270614B (zh) * 2022-07-18 2024-05-28 郑州轻工业大学 一种泥水循环系统多物理场数字孪生体可视化生成方法

Similar Documents

Publication Publication Date Title
Sheffer et al. Robust spherical parameterization of triangular meshes
Zhao et al. Mathematical morphology-based generalization of complex 3D building models incorporating semantic relationships
Li et al. Geometric structure simplification of 3D building models
Anderson et al. Smooth, volume-accurate material interface reconstruction
WO2021203711A1 (zh) 一种基于几何重建模型的等几何分析方法
Hu et al. Geometric optimization of building information models in MEP projects: algorithms and techniques for improving storage, transmission and display
WO2022104681A1 (zh) 一种柱坐标系等值面提取方法
Moenning et al. Fast marching farthest point sampling for implicit surfaces and point clouds
Jaillet et al. Fast Quadtree/Octree adaptive meshing and re-meshing with linear mixed elements
Guo et al. Line-based 3d building abstraction and polygonal surface reconstruction from images
Wu et al. Extracting POP: Pairwise orthogonal planes from point cloud using RANSAC
Zhang et al. A geometry and texture coupled flexible generalization of urban building models
CN111402422A (zh) 三维表面重建方法、装置和电子设备
CN113536416A (zh) 一种基于室内空间布局约束的场景模型补全方法
CN112487124A (zh) 一种使用VBA将CorelDraw地图中点状要素转换到SuperMap的方法
CN112070895A (zh) 一种高质量的实时等值面网络生成方法
Shen et al. An adaptive triangulation optimization algorithm based on empty circumcircle
Nielson MC/sup*: star functions for marching cubes
WO2022041169A1 (zh) 一种高质量的实时等值面网络生成方法
Shakaev et al. View-Dependent Level of Detail for Real-Time Rendering of Large Isosurfaces
Ma et al. Application of point cloud segmentation algorithm in high-precision virtual assembly
CN111027244B (zh) 一种百亿级颗粒模型的构建方法
CN107633543A (zh) 考虑局部拓扑结构的线条形状对应方法
Lopes et al. Interactive approaches to contouring and isosurfacing for geovisualization
Wang et al. Adaptive isosurface reconstruction using a volumetric-divergence-based metric

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 20961972

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20961972

Country of ref document: EP

Kind code of ref document: A1