时空查询之ECQL

前言

  ECQL 是 CQL 的扩展,CQL 是 OGC 标准查询语言,而 ECQL 是 GeoTools 为更好的方便查询,在编程实现时扩展了 CQL,主要扩展在于其移除了 CQL 的一些限制(属性必须在比较运算符的左边,不能创建 Id Filter 进行查询等限制),也和 SQL 更相似。所以可简单认为 CQL 是书面上的标准,而 ECQL 是事实上的标准。

前言

  ECQL 是 CQL 的扩展,CQL 是 OGC 标准查询语言,而 ECQL 是 GeoTools 为更好的方便查询,在编程实现时扩展了 CQL,主要扩展在于其移除了 CQL 的一些限制(属性必须在比较运算符的左边,不能创建 Id Filter 进行查询等限制),也和 SQL 更相似。所以可简单认为 CQL 是书面上的标准,而 ECQL 是事实上的标准。

谓词篇

时间查询主要有以下几个查询谓词:

谓词作用
T TEQUALS Time测试 T 和给定时间相等,相当于 T == Time。
T BEFORE Time测试 T 在给定时间之前,相当于 T < Time。
T BEFORE OR DURING Time Period测试 T 在给定时间段之前或其中,相当于 T <= TimePeriod[1]。
T DURING Time Period测试 T 在给定时间段其中,相当于 TimePeriod[0] <= T <= TimePeriod[1]。
T DURING OR AFTER Time Period测试 T 在给定时间段其中或之后,相当于 TimePeriod[0] <= T。
T AFTER Time测试 T 在给定时间之后,相当于 T > Time。

时间段以 / 分隔符区分前后两个时间,时间格式一般为 yyyy-MM-dd'T'HH:mm:ss.SSS'Z'。

空间查询主要有以下几个查询谓词:

谓词作用
INTERSECTS(A: Geometry, B: Geometry)测试 A 与 B 相交,与 DISJOINT 相反。
DISJOINT(A: Geometry, B: Geometry)测试 A 与 B 不相交,与 INTERSECTS 相反。
CONTAINS(A: Geometry, B: Geometry)测试 A 包含 B,与 WITHIN 相反。
WITHIN(A: Geometry, B: Geometry)测试 B 包含 A,即 A 在 B 中,与 CONTAINS 相反。
TOUCHES(A: Geometry, B: Geometry)测试 A 的边界是否与 B 的边界接触,但内部不相交。
CROSSES(A: Geometry, B: Geometry)测试 A 与 B 是否相交,但不存在包含关系。
OVERLAPS(A: Geometry, B: Geometry)测试 A 与 B 是否重叠,需满足 A 与 B 是同一类型(如都是 POLYGON),并且相交区域同样是 A 和 B 的类型(只能是 POLYGON,不能是 POINT)。
EQUALS(A: Geometry, B: Geometry)测试 A 与 B 完全相等。
RELATE(A: Geometry, B: Geometry, nineIntersectionModel: String)测试 A 与 B 是否满足 DE-9IM 模型,该模型可模拟上述所有情况。
DWITHIN(A: Geometry, B: Geometry, distance: double, units: String)测试 A 与 B 的最短距离是否不超过多少距离,单位有(feet, meters, statute miles, nautical miles, kilometers)。
BEYOND(A: Geometry, B: Geometry, distance: Double, units: String)测试 A 与 B 的最短距离是否超过多少距离。
BBOX(A: Geometry, leftBottomLng: Double, leftBottomLat: Double, rightTopLng: Double, rightTopLat: Double, crs="EPSG:4326")测试 A 是否与给定 box 相交。

Geometry 是指 WKT 格式的数据,主要有以下几种:

类型示例
POINTPOINT(6 10)
LINESTRINGLINESTRING(3 4,10 50,20 25)
POLYGONPOLYGON((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2))
MULTIPOINTMULTIPOINT(3.5 5.6, 4.8 10.5)
MULTILINESTRINGMULTILINESTRING((3 4,10 50,20 25),(-5 -8,-10 -8,-15 -4))
MULTIPOLYGONMULTIPOLYGON(((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2)),((6 3,9 2,9 4,6 3)))
GEOMETRYCOLLECTIONGEOMETRYCOLLECTION(POINT(4 6),LINESTRING(4 6,7 10))

※注: POLYGON 中的边界点必须闭合,即首尾点相同,若存在多个边界,则需要遵循 逆时针,顺时针,逆时针,顺时针... 的点排列顺序,逆时针封闭,顺时针开孔,以形成具有岛和洞的复杂多边形。

  由于 WKT 标准只支持二维的坐标,为支持三维坐标以及齐次线性计算,所以在 PostGIS 中又有 EWKT 标准实现,EWKT 扩展了 WKT,带 Z 结尾用来支持三维坐标,带 M 结尾用来支持齐次线性计算,如 POINTZ(6 10 3)POINTM(6 10 1)POINTZM(6 10 3 1),同时还支持坐标内嵌空间参考系,如 SRID=4326;LINESTRING(-134.921387 58.687767, -135.303391 59.092838)。GeoTools 19.0 之后也默认以 EWKT 进行解析和编码。

查询篇

属性字段查询

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// 查询属性 ATTR1 小于 7 的数据
Filter filter = ECQL.toFilter("ATTR1 < (1 + ((3 / 2) * 4))" );

// 查询属性 ATTR1 小于属性 ATTR2 绝对值的数据
Filter filter = ECQL.toFilter("ATTR1 < abs(ATTR2)" );

// 查询属性 ATTR1 为 test 字符串的数据
Filter filter = ECQL.toFilter("ATTR1 == 'test'" );

// 查询属性 ATTR1 在 10 和 20 之间的数据
Filter filter = ECQL.toFilter( "ATTR1 BETWEEN 10 AND 20" );
Filter filter = ECQL.toFilter( "ATTR1 >= 10 AND ATTR1 <= 20" );

// 多条件查询
Filter filter = ECQL.toFilter("ATTR1 < 10 AND ATTR2 < 2 OR ATTR3 > 10" );

// 查询属性 ATTR1 为 silver 或 oil 或 gold 的数据
Filter filter = ECQL.toFilter("ATTR1 IN ('silver','oil', 'gold' )");

// 以 ID 主键进行查询
Filter filter = ECQL.toFilter("IN ('river.1', 'river.2')");
Filter filter = ECQL.toFilter("IN (300, 301)");

模糊查询

1
2
3
4
5
6
7
8
9
10
11
// 查询属性 ATTR1 包含 abc 字符串的数据
Filter filter = ECQL.toFilter( "ATTR1 LIKE '%abc%'" );

// 查询属性 ATTR1 开头不为 abc 字符串的数据
Filter filter = ECQL.toFilter( "ATTR1 NOT LIKE 'abc%'" );

// 查询属性 cityName 开头为 new 的数据,忽略 new 的大小写
Filter filter = ECQL.toFilter("cityName ILIKE 'new%'");

// 测试字符串是否包含
Filter filter = ECQL.toFilter("'aabbcc' LIKE '%bb%'");

空属性查询

1
2
3
4
5
6
7
8
// 查询有属性 ATTR1 存在的数据
Filter filter = ECQL.toFilter( "ATTR1 EXISTS" );

// 查询属性 ATTR1 不存在的数据
Filter filter = ECQL.toFilter( "ATTR1 DOES-NOT-EXIST" );

// 查询 Name 为 NULL 的数据
Filter filter = ECQL.toFilter("Name IS NULL");

时间查询

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
// 查询时间属性 dtg 等于的数据
Filter filter = ECQL.toFilter( "dtg TEQUALS 2006-11-30T01:30:00Z" );

// 查询时间属性 dtg 在之后的数据
Filter filter = ECQL.toFilter("dtg AFTER 2006-11-30T01:30:00Z");

// 查询时间属性 dtg 在之前的数据
Filter filter = ECQL.toFilter("dtg BEFORE 2006-11-30T01:30:00Z");

// 查询时间属性 dtg 在之间的数据,+3:00 代表 GMT 时间 +3 小时,以 Z 结尾的时间就是 GMT 时间
Filter filter = ECQL.toFilter( "dtg DURING 2006-11-30T00:30:00+03:00/2006-11-30T01:30:00+03:00 ");

// 查询时间属性 dtg 等于的数据
Filter filter = ECQL.toFilter("dtg = 1981-06-20");

// 查询时间属性 dtg 小于等于的数据
Filter filter = ECQL.toFilter("dtg <= 1981-06-20T12:30:01Z");

空间查询

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// 查询空间属性 geom 包含点的数据
Filter filter = ECQL.toFilter( "CONTAINS(geom, POINT(1 2))" );

// 查询空间属性 geom 与 box 相交的数据
Filter filter = ECQL.toFilter( "BBOX(geom, 10,20,30,40)" );

// 查询空间属性 geom 与点最短距离不超过 10 千米的数据
Filter filter = ECQL.toFilter( "DWITHIN(geom, POINT(1 2), 10, kilometers)" );

// 查询空间属性 geom 与线相交的数据(geom 也必须是线)
Filter filter = ECQL.toFilter( "CROSS(geom, LINESTRING(1 2, 10 15))" );

// 查询空间属性 geom 与 GEOMETRYCOLLECTION 相交的数据(geom 也必须是 GEOMETRYCOLLECTION)
Filter filter = ECQL.toFilter( "INTERSECT(geom, GEOMETRYCOLLECTION (POINT (10 10),POINT (30 30),LINESTRING (15 15, 20 20)) )" );

// 查询空间属性 geom 与线相交的数据
Filter filter = ECQL.toFilter( "CROSSES(geom, LINESTRING(1 2, 10 15))" );

// 查询空间属性 geom 与 GEOMETRYCOLLECTION 相交的数据
Filter filter = ECQL.toFilter( "INTERSECTS(geom, GEOMETRYCOLLECTION (POINT (10 10),POINT (30 30),LINESTRING (15 15, 20 20)) )" );

// 查询空间属性 geom 与包含线的数据
Filter filter = ECQL.toFilter("RELATE(geom, LINESTRING (-134.921387 58.687767, -135.303391 59.092838), T*****FF*)");

  在 GeoTools 中,可通过 FilterFactory 来构造 Filter,而不是直接写字符串,具体示例如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2();

// 相当于 Filter filter1 = ECQL.toFilter("ATTR1 = 1 AND ATTR2 < 4" );
List<Filter> filterList = ECQL.toFilterList("ATTR1=1; ATTR2<4");
Filter filter1 = ff.and(filterList);

// 相当于 Filter filter2 = ECQL.toFilter( "BBOX(geom, 10,20,30,40)" );
Filter filter2 = ff.bbox("geom", 10, 20, 30, 40, "EPSG:4326");

// 相当于 Filter filter3 = ECQL.toFilter( "dtg DURING 2006-11-29T00:30:00Z/2006-11-30T00:30:00Z");
Date startTime = ZonedDateTime.of(2006, 11, 29, 0, 30, 0, 0, ZoneOffset.UTC);
Date endTime = Date.from(startTime.plusDays(1).toInstant());
Filter filter3 = ff.between(ff.property("dtg"), ff.literal(startTime), ff.literal(endTime));

后记

  基本可认为 CQL 和 SQL 中查询条件差不多,虽然不支持分组查询等复杂 SQL 特性,但对于一般的时空查询基本够用,CQL 中还有些空间操作函数就不继续写了,如取面积,取缓冲区,取交集,取长度等等,有需要的可自行查询 uDig Common Query Language

参考资料

GeoTools CQL

GeoTools ECQL

GeoServer ECQL Reference / GeoServer 属性查询和空间查询支持 CQL / ECQL过滤器语言

WKT解读

GEOS库学习之三:空间关系、DE-9IM和谓词

GeoMesa踩坑指北

前言

  需要做个 GeoMesa 的微服务,简单熟悉一下 GeoMesa。

前言

  需要做个 GeoMesa 的微服务,简单熟悉一下 GeoMesa。

基础篇

  GeoMesa 可以说是大数据中的 PostGIS,主要用来在存储和处理 GIS 数据时提供相应的索引,从而加快处理速度。GeoMesa 基于 GeoTools,其中最重要的两个概念就是 SimpleFeatureType 和 SimpleFeature,SimpleFeatureType 对应的是关系型数据库中表的描述(表明,表的列字段属性信息等),而 SimpleFeature 对应的是表中每行数据。下面重点谈谈 GeoMesa 中的 SimpleFeatureType 以及其创建索引方式。

  在 GeoMesa 中通常使用 SimpleFeatureTypes.createType 方法进行创建,该方法有两个重载,以没有 namespace 参数的方法为例:

1
2
3
4
def createType(typeName: String, spec: String): SimpleFeatureType = {
val (namespace, name) = parseTypeName(typeName)
createType(namespace, name, spec)
}

先通过 parseTypeName 解析 typeName,以 : 作为分隔符,取最后一个有效(不为空)字符串作为表名(name),其余部分如有效则作为 namespace,否则 namespace 则为 null。spec 参数的通用形式有以下几种:

1
2
3
4
5
6
7
val spec = "name:String,dtg:Date,*geom:Point:srid=4326"

val spec = "name:String,dtg:Date,*geom:Point:srid=4326;geomesa.indices.enabled='z2,id,z3'"

val spec = "name:String:index=true,tags:String:json=true,dtg:Date:default=true,*geom:Point:srid=4326;geomesa.indices.enabled='z2,id,z3'"

val spec = "userId:String,trackId:String,altitude:Double,dtg:Date,*geom:Point:srid=4326;geomesa.index.dtg='dtg',geomesa.table.sharing='true',geomesa.indices='z3:4:3,z2:3:3,id:2:3',geomesa.table.sharing.prefix='\\u0001'"

先使用 ; 分隔符,再使用 , 分隔符,最后使用 : 分隔符。; 分隔符将 spec 分割为两个字符串:前者表示表中的全部列属性信息,列属性经过 , 分隔符分割为多列,列又经过 : 分隔符分割为 列名,列数据类型,列的一些属性(是否是索引,json 数据,默认索引等),而列名首字母 * 代表该字段是用于索引的 geometry 类型,一般采用 WKT 格式进行描述,当然存在数据库时会以字节码进行压缩;后者表示创建表时的 userData,同样经过 , 分隔符分割为多个 userData,userData 的一些默认属性可在 SimpleFeatureTypes.Configs 中看到,其它的可以用户自定义,这里重点说一下 geomesa.indices.enabled 属性,目前 GeoMesa 支持 8 种索引,分别为:

1
2
3
4
5
6
7
8
"attr", // 属性索引
"id", // 主键索引
"s2", // Hilbert 曲线点空间索引
"s3", // Hilbert 曲线点时空索引
"z2", // Z 型曲线点空间索引
"xz2", // Z 型曲线线面空间索引
"z3", // Z 型曲线点时空索引
"xz3" // Z 型曲线线面时空索引

  由于 GeoMesa 中的索引一般存在多个版本,而 geomesa.indices.enabled 默认使用最新的版本,若需要指定版本,需要使用 geomesa.indices该属性是 geomesa 内部属性,不对外开放,通用格式为:

1
s"$name:$version:${mode.flag}:${attributes.mkString(":")}"

name 代表索引类别,version 代表索引版本,mode.flag 代表索引模式(是否支持读写,一般为3,支持读也支持写),attributes 代表是哪些字段需要建立该索引。spec 参数可以只有描述列属性的字段,即不带任何 useData 信息,GeoMesa 会默认添加索引信息,若存在空间和时间字段,则会默认建立 z3(空间字段为点 Point 类型) 或 xz3(空间字段为线面 非Point 类型) 索引,若有多个空间和时间字段,建立索引的字段为第一个空间和第一个时间字段;若只存在空间字段,则会建立 z2 或 xz2 索引;若只有时间字段,则默认建立时间属性索引。当然如没有在 spec 指明索引信息,可以在后续继续添加信息,如下:

1
2
3
4
5
6
7
8
9
10
import org.locationtech.geomesa.utils.interop.SimpleFeatureTypes;

String spec = "name:String,dtg:Date,*geom:Point:srid=4326";
SimpleFeatureType sft = SimpleFeatureTypes.createType("mySft", spec);
// enable a default z3 and a default attribute index
sft.getUserData().put("geomesa.indices.enabled", "z3,attr:name");
// or, enable a default z3 and an attribute index with a Z2 secondary index
sft.getUserData().put("geomesa.indices.enabled", "z3,attr:name:geom");
// or, enable a default z3 and an attribute index with a temporal secondary index
sft.getUserData().put("geomesa.indices.enabled", "z3,attr:name:dtg");

坑篇

导入 OSM 数据问题

  在导入 osm 数据时,若使用 osm-ways 作为 SimpleFeatureType,则 geomesa 会使用数据库存储 node 临时使用,这时其默认使用 H2 Database,若想使用其它数据库,则需要在 lib 导入相应 jdbc 包,若使用 postgresql 数据库,则 geomesa 会触发一个 bug,因为 postgresql 没有 double 类型,只有 double precision 类型,这将导致建表出错。详情见 geomesa/geomesa-convert/geomesa-convert-osm/src/main/scala/org/locationtech/geomesa/convert/osm/OsmWaysConverter.scala 中

1
2
3
4
private def createNodesTable(): Unit = {
val sql = "create table nodes(id BIGINT NOT NULL PRIMARY KEY, lon DOUBLE, lat DOUBLE);"
WithClose(connection.prepareStatement(sql))(_.execute())
}

所以若需要使用 geomesa-convert-osm 导入 osm 数据时,需要进入 geomesa/geomesa-convert/geomesa-convert-osm 文件夹中输入命令

1
mvn dependency:copy-dependencies -DoutputDirectory=./depLib

导出 geomesa-convert-osm 依赖包,将其中的 h2,osm4j,dynsax,trove4j 等一系列库放入 $GEOMESA_HBASE_HOME/lib 中。

s2 索引问题

  s2 索引即 Google S2 Geometry 算法基于 Hilbert 曲线生成一种索引,GeoMesa 的 s2 索引是一个国人提交的,目前 3.2 版本只支持点的时空索引,不支持线面的时空索引,当然官方也在实现自己的 Hilbert 曲线,希望后续 GeoMesa 中会有 h2 索引。Shaun 在导入 osm 数据并启用 s2 索引时,报错,被提示不支持,对比 geomesa-index-api2Index.scala 和 geomesa-index-api2Index.scala 两文件的 defaults 函数可发现 S2Index 直接返回空,而在 geomesa-index-api.scala 中 fromName 函数需要调用 defaults 函数,从而导致 s2 索引不支持,修改 S2Index 的 defaults 函数即可(别忘了在 S2Index 类中首行加上 import org.locationtech.geomesa.utils.geotools.RichSimpleFeatureType.RichSimpleFeatureType)。

后记

  暂时就了解了这么多,等后续熟悉的更多再继续更吧 (ง •_•)ง。

附录

GeoMesa 命令行工具部分参数

Geomesa 命令行参数:

参数描述
-c, --catalog *存放 schema 元数据的catalog 表(相当于数据库)
-f, --feature-nameschema 名(相当于数据库中的表)
-s, --spec要创建 SimpleFeatureType 的说明(即表中列的描述信息,表的 schema,如 "name:String,age:Int,dtg:Date,*geom:Point:srid=4326")
-C, --converter指定转换器,必须为一下之一:1、已经在classpath中的converter 名;2、converter 的配置(一个字符串);3、包括converter的配置的名
–converter-error-mode自定义的转换器的error mode
-t, --threads指定并行度
–input-format指定输入源格式(如csv, tsv, avro, shp, json,)
–no-tracking指定提交的 ingest job何时终止(在脚本中常用)
–run-mode指定运行模式,必须为:local(本地)、distributed (分布式)、distributedcombine(分布式组合)之一
–split-max-size在分布式中,指定切片最大大小(字节)
–src-list输入文件为文本文件,按行输入
–force禁用任何的提示
[files]…指定输入的文件

参考资料:GeoMesa命令行工具---摄取命令

IDEA使用Docker环境开发调试

前言

  IDEA 以前基本没用过,只是简单用过 Android Studio,还基本都忘记了 ( ╯□╰ ),以后应该会用 Scala 做一些大数据方面的东西,而大数据的环境都是 Linux 下的,而 Shaun 日常都是在 Windows 下开发,所以需要用日前做的容器环境来测试调试运行程序,简单记录一下 IDEA 在这方面的使用方法。

前言

  IDEA 以前基本没用过,只是简单用过 Android Studio,还基本都忘记了 ( ╯□╰ ),以后应该会用 Scala 做一些大数据方面的东西,而大数据的环境都是 Linux 下的,而 Shaun 日常都是在 Windows 下开发,所以需要用日前做的容器环境来测试调试运行程序,简单记录一下 IDEA 在这方面的使用方法。

运行篇

  右键项目名(HelloWorld),新建文件(New =》File),指定文件名为 Dockerfile 。写入内容示例如下:

1
2
3
4
FROM stc:2.0
COPY ./target/classes/ /tmp
WORKDIR /tmp
ENTRYPOINT ["scala","HelloWorld"]

点击左上角绿色双箭头,可编辑 Dockerfile(Edit 'Dockerfile') ,指定当前上下文目录(Context folder),Contaier name 等容器启动选项。直接运行 Dockerfile(Run 'Dockerfile'),IDEA 即可自动创建容器,并在容器中运行程序,程序运行完则容器自动停止,若需要运行存在外部依赖的程序,则只能以 jar 包的方式运行。

  设置 IDEA 生成 jar 包如下:在最上面的菜单栏中 File =》Project Structure =》Artifacts =》+ =》JAR =》From modules with dependencies,选择 Main Class,点击右边的文件夹图标即可选择相应类,由于存在外部依赖,所以不能直接用默认的 extract to the target JAR,而是应该选择下面的 copy to the output directory and link via manifest,点击 OK 后,自动或手动选择导出的依赖 jar 包,点击 OK。在最上面的菜单栏中 Build =》Build Artifacts...,可在 out/artifacts/HelloWorld_jar 文件夹中生成所有 jar 包。之后编辑 Dockerfile, 更改 Dockerfile 上下文目录为 out/artifacts/HelloWorld_jar ,指定容器名,在 Command 中输入 java -jar HelloWorld.jar 修改 Dockerfile 中第 2 行命令为 COPY . /tmp,修改第 4 行命令为 CMD ["java", "-jar", "HelloWorld.jar"]。之后运行 Dockerfile 即可在下面 Services 栏对应 Docker 容器 Attached Console 中看到程序运行结果。

调试篇

  除了使用 IDEA 生成 jar 包外,还需要使用 IDEA 的远程调试功能,设置 IDEA 远程调试功能如下:在最上面的菜单栏中 Run =》Edit Configurations... =》+ =》Remote JVM Debug,上方的 Debugger mode 中使用默认的 Attach to remote JVM, 在下面的 Before launch 添加 Launch Docker before debug。在弹窗中选择相应 Dockerfile,在下方的 Custom command 中输入 java -agentlib:jdwp=transport=dt_socket,server=y,suspend=y,address=5005 -jar HelloWorld.jar, 完成后即可使用该配置在 IDEA 调试容器中运行的程序。

后记

  用这种方式使用 IDEA 确实达到了 Shaun 理想的结果,Windows 下开发,Docker 中调试和运行,应付简单的代码调试和运行确实是没问题,但是在复杂的分布式环境下总会碰到一些莫名奇妙的问题,这些问题就是纯粹的经验了。

参考资料

Run a Java application in a Docker container

Debug a Java application using a Dockerfile

大数据环境搭建笔记

前言

  准备开始搞时空数据了,先简单搭一下环境。

前言

  准备开始搞时空数据了,先简单搭一下环境。

准备搭的环境为:jdk-1.8.0,hadoop-3.2.1,hbase-2.2.6,geomesa-hbase_2.11-3.1.0,spark-3.0.1-bin-hadoop3.2,geoserver-2.16.5-bin,geomesa-hbase_2.11-3.2.0-SNAPSHOT,所用的包都已下好并解压到 /home 目录下。

※注hbase-2.2.6 暂不支持最新的 hadoop-3.3.0,Hadoop 也最好使用 jdk-1.8.0,java-11 会有问题。

Hadoop 环境

  首先修改 /etc/hosts 文件中本机 ip 对应的名称为 master,若在容器中安装则需要在 run 开启容器就指定 --hostname master,否则改了也没用,下次启动容器时 hostname 又会回到初始状态,下面开启正式的配置。

修改 /home/hadoop-3.2.1/etc/hadoop/hadoop-env.sh 文件,添加

1
export JAVA_HOME=$JAVA_HOME

修改 /home/hadoop-3.2.1/etc/hadoop/core-site.xml 文件,添加

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
<configuration>

<!-- master 前面配置的主机名称 -->
<!-- <property>
<name>fs.default.name</name>
<value>hdfs://master:9000</value>
</property> -->

<property>
<name>fs.defaultFS</name>
<value>hdfs://master:9000</value>
</property>

<property>
<name>hadoop.tmp.dir</name>
<value>/home/hadoop/data/tmp</value>
</property>

</configuration>

修改 /home/hadoop-3.2.1/etc/hadoop/hdfs-site.xml 文件,添加

1
2
3
4
5
6
7
8
9
10
11
12
13
<configuration>

<property>
<!--指定SecondaryNameNode位置-->
<name>dfs.namenode.secondary.http-address</name>
<value>master:9001</value>
</property>
<property>
<name>dfs.replication</name>
<value>1</value>
</property>

</configuration>

修改 /home/hadoop-3.2.1/etc/hadoop/yarn-site.xml 文件,添加

1
2
3
4
5
6
7
8
9
10
<configuration>

<!-- Site specific YARN configuration properties -->

<property>
<name>yarn.nodemanager.aux-services</name>
<value>mapreduce_shuffle</value>
</property>

</configuration>

修改 /home/hadoop-3.2.1/etc/hadoop/mapred-site.xml 文件,添加

1
2
3
4
5
6
7
8
<configuration>

<property>
<name>mapreduce.framework.name</name>
<value>yarn</value>
</property>

</configuration>

在 /home/hadoop-3.2.1/sbin/start-dfs.sh 和 /home/hadoop-3.2.1/sbin/stop-dfs.sh 文件头添加

1
2
3
4
5
#!/usr/bin/env bash
HDFS_DATANODE_USER=root
HDFS_DATANODE_SECURE_USER=hdfs
HDFS_NAMENODE_USER=root
HDFS_SECONDARYNAMENODE_USER=root

在 /home/hadoop-3.2.1/sbin/start-yarn.sh 和 /home/hadoop-3.2.1/sbin/stop-yarn.sh 文件头添加

1
2
3
4
#!/usr/bin/env bash
YARN_RESOURCEMANAGER_USER=root
HADOOP_SECURE_DN_USER=yarn
YARN_NODEMANAGER_USER=root

设置环境变量,在 /etc/profile 中添加

1
2
3
4
5
#Hadoop Environment Setting
export HADOOP_HOME=/home/hadoop-3.2.1
export JAVA_LIBRARY_PATH=$HADOOP_HOME/lib/native
export LD_LIBRARY_PATH=$JAVA_LIBRARY_PATH
export PATH=$PATH:$HADOOP_HOME/bin:$HADOOP_HOME/sbin

由于容器中默认为 root 用户,所以在 /root/.bashrc 文件末尾添加 source /etc/profile,以开机启用设置的环境变量。

在启动 Hadoop 之前需要执行 hdfs namenode -format 进行格式化,启动命令为 /home/hadoop-3.2.1/sbin/start-all.sh后续若需要清空并重新设置 Hadoop 时,必须先删除 /home/hadoop/ 目录,再重新进行格式化。

HBase 环境

修改 /home/hbase-2.2.6/conf/hbase-env.sh 文件,添加

1
2
3
export JAVA_HOME=$JAVA_HOME
# 使用自带的ZooKeeper管理
export HBASE_MANAGES_ZK=true

修改 /home/hbase-2.2.6/conf/hbase-site.xml 文件,添加

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
<configuration>

<property>
<name>hbase.rootdir</name>
<value>hdfs://master:9000/hbase</value>
</property>
<property>
<name>hbase.dynamic.jars.dir</name>
<value>hdfs://master:9000/hbase/lib</value>
</property>
<property>
<name>hbase.cluster.distributed</name>
<value>true</value>
</property>
<property>
<name>dfs.replication</name>
<value>1</value>
</property>
<property>
<name>hbase.master.maxclockskew</name>
<value>180000</value>
<description>Time difference of regionserver from master</description>
</property>
<property>
<name>hbase.zookeeper.quorum</name>
<value>localhost</value>
</property>
<!-- 修改默认8080 端口-->
<property>
<name>hbase.rest.port</name>
<value>8088</value>
</property>

<!-- 2181 默认端口,尽量不要修改,geomesa-hbase 导入数据时默认连接端口为 2181-->
<property>
<name>hbase.zookeeper.property.clientPort</name>
<value>2181</value>
</property>
<property>
<name>hbase.zookeeper.property.dataDir</name>
<value>/home/hbase/data</value>
</property>
<property>
<name>hbase.unsafe.stream.capability.enforce</name>
<value>false</value>
</property>

<!-- geomesa-hbase -->
<property>
<name>hbase.coprocessor.user.region.classes</name>
<value>org.locationtech.geomesa.hbase.server.coprocessor.GeoMesaCoprocessor</value>
</property>

</configuration>

修改 /home/hbase-2.2.6/conf/regionservers 文件,修改为(原来为 localhost)

1
master

设置环境变量,在 /etc/profile 中添加

1
2
3
#HBase Environment Setting
export HBASE_HOME=/home/hbase-2.2.6
export PATH=$PATH:$HBASE_HOME/bin

配置好之后,执行 start-hbase.sh 启动 HBase。

Spark 环境

修改 /home/spark-3.0.1-bin-hadoop3.2/conf/spark-env.sh 文件,在文件末尾添加

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# 配置JAVA_HOME,一般来说,不配置也可以,但是可能会出现问题,还是配上吧
export JAVA_HOME=$JAVA_HOME
# 一般来说,spark任务有很大可能性需要去HDFS上读取文件,所以配置上
# 如果说你的spark就读取本地文件,也不需要yarn管理,不用配
export HADOOP_CONF_DIR=$HADOOP_HOME/etc/hadoop

# 设置Master的主机名
export SPARK_MASTER_HOST=master
# 提交Application的端口,默认就是这个,万一要改呢,改这里
export SPARK_MASTER_PORT=7077
# 每一个Worker最多可以使用的cpu core的个数,我虚拟机就一个...
# 真实服务器如果有32个,你可以设置为32个
export SPARK_WORKER_CORES=1
# 每一个Worker最多可以使用的内存,我的虚拟机就2g
# 真实服务器如果有128G,你可以设置为100G
export SPARK_WORKER_MEMORY=2g
# master web UI端口默认8080
export SPARK_MASTER_WEBUI_PORT=8090
# worker web UI端口默认8081
export SPARK_WORKER_WEBUI_PORT=8089

复制 /home/spark-3.0.1-bin-hadoop3.2/conf/slaves.template 文件,并重命名为 slaves,将该文件尾修改为

1
2
# 里面的内容原来为localhost,改为master 
master

设置环境变量,在 /etc/profile 中添加

1
2
export SPARK_HOME=/home/spark-3.0.1-bin-hadoop3.2
export PATH=$PATH:$SPARK_HOME/bin:$SPARK_HOME/sbin

将 /home/spark-3.0.1-bin-hadoop3.2/sbin/start-all.sh 重命名为 start-spark-all.sh,将 /home/spark-3.0.1-bin-hadoop3.2/sbin/stop-all.sh 重命名为 stop-spark-all.sh,执行 start-spark-all.sh 启动 Spark。

geomesa-hbase 环境

编译 geomesa

克隆 LocationTech GeoMesa修改 pom.xml,即修改对应依赖的 hadoop 和 hbase 以及 spark 版本(spark 最新的3.0.1版本由 Scala-2.12 编译,而 Geomesa 编译目前采用 Scala-2.11, 所以 Spark 不能使用最新的版本,只能用 2.4.7)。进入 geomesa 根目录,使用命令

1
2
3
4
mvn clean install -DskipTests

# 或仅编译 geomesa-hbase
mvn clean install -pl geomesa-hbase -am -DskipTests

编译 geomesa,中间可能会失败很多次,包下不来,可能需要挂代理或换源,重复使用命令多次即可。

配置 geomesa-hbase

将 /home/geomesa/geomesa-hbase/geomesa-hbase-dist/target/geomesa-hbase_2.11-3.2.0-SNAPSHOT-bin.tar.gz 解压为 /home/geomesa-hbase_2.11-3.2.0-SNAPSHOT,将 /home/geomesa-hbase_2.11-3.2.0-SNAPSHOT/dist/hbase/geomesa-hbase-distributed-runtime-hbase2_2.11-3.2.0-SNAPSHOT.jar 复制到 /home/hbase-2.2.6/lib/ 文件夹中,修改 /home/geomesa-hbase_2.11-3.2.0-SNAPSHOT/conf/dependencies.sh 文件,设置正确的Hadoop 和 hbase 版本,依次执行 /home/geomesa-hbase_2.11-3.2.0-SNAPSHOT/bin/install-dependencies.sh 和 /home/geomesa-hbase_2.11-3.2.0-SNAPSHOT/bin/install-shapefile-support.sh。设置环境变量,在 /etc/profile 中添加

1
2
3
4
5
export GEOMESA_HBASE_HOME=/home/geomesa-hbase_2.11-3.2.0-SNAPSHOT
export GEOMESA_LIB=$GEOMESA_HBASE_HOME/lib
export GEOMESA_CONF_DIR=${GEOMESA_HBASE_HOME}/conf
export CLASSPATH=$CLASSPATH:$GEOMESA_LIB:$GEOMESA_CONF_DIR
export PATH=$PATH:$GEOMESA_HBASE_HOME/bin

测试 geomesa-hbase

启动 Hadoop 和 HBase 之后,可直接使用命令

1
geomesa-hbase ingest --catalog TestGeomesa --feature-name road --input-format shp "/home/shpdata/road.shp"

导入 shp 数据,shp 不能有 id 字段,因为 Geomesa 在创建表时会默认生成一个 id 字段。

也可克隆 geomesa-tutorials ,同样修改其中的 pom.xml 文件,进入 geomesa-tutorials 根目录,使用命令

1
mvn clean install -pl geomesa-tutorials-hbase/geomesa-tutorials-hbase-quickstart -am

编译 geomesa-tutorials,编译完成后,使用命令

1
java -cp geomesa-tutorials-hbase/geomesa-tutorials-hbase-quickstart/target/geomesa-tutorials-hbase-quickstart-3.2.0-SNAPSHOT.jar org.geomesa.example.hbase.HBaseQuickStart --hbase.zookeepers localhost --hbase.catalog geomesaTest

导入数据进 Hbase,导入成功后可通过 hbase shell 进入 hbase,在 hbase shell 中通过 list 查看 hbase 现有的表。

整合 geoserver

导入依赖插件

1
manage-geoserver-plugins.sh -l ${GEOSERVER_HOME}/webapps/geoserver/WEB-INF/lib/ -i

修改 /home/geomesa-hbase_2.11-3.2.0-SNAPSHOT/bin/install-dependencies.sh 中第33行:

1
2
# install_dir="${GEOMESA_HBASE_HOME}/lib"
install_dir="${GEOSERVER_HOME}/webapps/geoserver/WEB-INF/classes"

执行 install-dependencies.sh 安装插件,安装完后将 classes 中的 lib 都移到 ${GEOSERVER_HOME}/webapps/geoserver/WEB-INF/lib中。

后记

  环境搞起来真麻烦,在编译和运行 Geomesa 时总能遇到一些莫名奇妙的问题,Java 系的这一套确实很麻烦,尤其是各种依赖关系,不过最后总算是搞好了,能直接在 geoserver 中看到 geomesa 存在 hbase 里的地图。

参考资料

Centos7系统 Hadoop+HBase+Spark环境搭建

GeoMesa-HBase操作篇——安装

centos7安装geomesa2.0.2_hbase_geoserver2.13.2的方法

hadoop fs 命令使用

GeoMesa HBase Quick Start

Installing GeoMesa HBase

Spark完全分布式集群搭建【Spark2.4.4+Hadoop3.2.1】

附录

最后附上一些常用的端口及说明:

Hbase

配置端口说明
hbase.master.port16000HMaster绑定端口
hbase.master.info.port16010HBase Master的Web UI端口
hbase.regionserver.port16020HBase RegionServer绑定的端口
hbase.regionserver.info.port16030HBase RegionServer的Web UI端口
hbase.zookeeper.property.clientPort2181Zookeeper客户端连接端口
hbase.zookeeper.peerport2888Zookeeper节点内部之间通信的端口
hbase.zookeeper.leaderport3888Zookeeper用来选举主节点的端口
hbase.rest.port8080HBase REST server的端口
hbase.master.port60000HMaster的RPC端口
hbase.master.info.port60010HMaster的http端口
hbase.regionserver.port60020HRegionServer的RPC端口
hbase.regionserver.info.port60030HRegionServer的http端口

Hadoop

配置端口说明
fs.defaultFS9000hdfs访问端口
dfs.namenode.rpc-address9001DataNode会连接这个端口
dfs.datanode.address9866DataNode的数据传输端口
dfs.namenode.http-address9870namenode的web UI 端口
yarn.resourcemanager.webapp.address8088YARN的http端口

Spark

端口说明
8080master的webUI,Tomcat的端口号(已修改为8090)
8081worker的webUI的端口号(已修改为8089)
18080historyServer的webUI的端口号

需开放端口 22,2181,5432,8080,8088,8089,8090,9870,16010,16030。

docker run -dit --privileged=true --name STC2 --hostname master -v E:/Docker/ShareFile:/mnt/sharefile -p 22:22 -p 80:80 -p 2181:2181 -p 5432:5432 -p 8080-8090:8080-8090 -p 9870:9870 -p 16010:16010 -p 16030:16030 stc:2.0 init

设计模式浅谈

前言

  进入职场一年半以来,Shaun 完全独立从 0 到 1 开发了 1.5 个项目(当然也有参与其它项目,但不是 Shaun 独立从 0 到 1 开发的,没多少控制权,就不谈了),一个网页版的高精地图编辑器,半个地图可视化系统,这里面 Shaun 用了不少设计模式,这篇就谈谈 Shaun 用过的和没用过的一些设计模式。

前言

  进入职场一年半以来,Shaun 完全独立从 0 到 1 开发了 1.5 个项目(当然也有参与其它项目,但不是 Shaun 独立从 0 到 1 开发的,没多少控制权,就不谈了),一个网页版的高精地图编辑器,半个地图可视化系统,这里面 Shaun 用了不少设计模式,这篇就谈谈 Shaun 用过的和没用过的一些设计模式。

  以「Head First 设计模式」为参考,Shaun 用 C++ 实现了一遍书中的例子(代理模式及其后面的模式除外),下面进入正文。

模式篇

策略模式

  Shaun 个人认为最能体现面向对象编程思想(抽象封装继承多态)的一种模式,换句话说,只要真正理解和运用面向对象编程,一定会自然而然的用到策略模式。Shaun 在做高精地图编辑器时,需要设计一个渲染模块,渲染模块会包含高亮行为,高亮有两种,一种是直接改变颜色,一种是使用后期处理(OutlinePass 或 UnrealBloomPass 等)进行高亮,这时就需要在渲染类中组合高亮行为。

  策略模式中涉及到的原则有:1、封装变化;2、多用组合,少用继承;3、针对接口编程,不针对实现编程。封装变化这点很考验程序员的经验水平,在写代码之初,往往预料不到变化,所以这一点一般是在编码过程中逐渐完善的,不断进行抽象,从而生成比较合理的基类;第二点一般也是对的,但有时在编码过程中难免会碰到到底是用继承还是组合的问题,这时候可以多想想,组合并不是万能的,有时继承更合适,这时可以请教身边更有经验的程序员,组合的优势在于当子类不想要这个对象时,可以随时丢弃,而继承的优势在于,当子类不想实现这个行为时,可以有默认的行为,而且有些时候只能用继承;针对接口编程没啥好说的,就是抽象。

观察者模式

  这个模式如果在分布式系统中又叫发布订阅模式,该模式常用于消息通知。前端有个 RxJS 的库将这一模式玩出花来了,Shaun 在高精地图编辑器的事件流管理中就使用了该库。在 threejs 中所有渲染对象的都有一个统一的基类 EventDispatcher,该类中就实现了一个简单的观察者模式,用来处理事件监听触发和通知,和真正的观察者相比,区别在于观察者只是一个函数,而不是一个类,其实浏览器的事件监听机制实现方式也和这个类差不多。

  观察者模式中涉及到原则有:松耦合。这里的松耦合是指主题和观察者之间是隔离的,观察者可自行实现自己的更新行为,而主题同样可实现自己的通知机制,两者虽有关联但互不影响。松耦合原则说起来人人会说,但真正能实现松耦合的却不多,实现松耦合的关键在于怎样分离两个系统,两个系统的连接点在哪,这有时很难理清,从而造成逻辑混乱,bug 丛生。

装饰者模式

  利用该模式可以很方便的扩展一些方法和属性,尤其是遇到主体和配件这样的关系时,可以很轻松的添加配件到主体上。Shaun 没用过这个模式,本来在扩展 threejs 一个类时想用,但确实没找到非常明确的主体和配件这样的关系,最后还是简单的使用继承了。

  装饰者模式涉及到的原则有:开放——封闭原则。设计一个类需要考虑对扩展开放,对修改关闭。修改和扩展看似矛盾,但实则可以独立存在,装饰者的配件可以无限加,这是扩展,是开放,而在加配件时无需修改现有代码,这是封闭。当然这一原则并不独属于装饰者模式,应该说是只要用到面向对象的思想开发程序,就一定会用到该原则,否则也就失去了面向对象的意义。但有时这个原则又没必要贯彻彻底,因为对于有些需求可能很难弄清修改和扩展的界限,这时就达到能尽量重用父类的方法就好。

工厂模式

  该模式在稍微大一点的系统中应该都会用到,根据不同的输入,生成不同的对象,这就是工厂模式的本质。至于工厂模式的实现方式一般会根据需求的复杂度来决定:1、只有一个工厂,一类产品,只是为了集中一层 if-else,可用简单工厂模式,甚至一个 builder 函数即可;2、有多个工厂,还是只有一类产品,用工厂模式,多个工厂继承一个工厂父类即可,相当于多个简单工厂组合;3、有多个工厂,多类产品,哪个工厂生产什么产品可能有变化,这时需要用到抽象工厂模式,除正常的继承之外,还需使用组合,组合组成产品的父类,相当于再组合一个工厂。Shaun 在高精地图编辑器中当然是大量使用的工厂模式和简单工厂模式,主要是为了集中 if-else 的处理,比如根据不同的数据类型创建不同的属性栏界面(枚举用下拉框,字符串用文本框,数字用数字栏等),根据不同的路网元素创建对应的渲染器对象以及对应的属性界面等。

  工厂模式涉及到的原则有:依赖倒置原则。尽量依赖抽象,而不是具体类。这其实也是抽象一种作用或好处,即在使用过程中尽量使用最上层的父类,具体类只在创建实例时使用。

单例模式

  写程序的基本都会用到该模式,主要用来创建全局唯一对象,可用来存储和处理系统中常用的各个模块都会用到的一些数据。Shaun 在编辑器中单例模式用了好几个,比如全局唯一的 viewport,用力绘制 3d 图形;全局唯一的路网数据;当然系统中存在太多的单例模式也不好,最好是只有一个,如 Shaun 的编辑器中最好的模式就是创建一个单例的 Editor 类,需要做单例的对象都可以放在该类中,如此保证系统中只有一个单例类,以进行统一管理。

  该模式与面向对象倒是没多大关系了,可以认为是全局变量的优化版,毕竟大的系统中全局变量基本不可避免,这时就可以使用单例模式。

命令模式

  该模式主要用来将函数方法封装成类,这样做的好处就是可以更灵活的执行该方法(将方法放进队列中依次执行,将方法持久化以便系统启动执行),同时也可以保存该方法的一些中间状态,以便撤销操作,回到系统执行该方法前的状态。Shaun 在编辑器中主要用命令模式做撤销重做功能,这基本也是编辑器的必备功能了,可以说没有撤销重做功能的编辑器是不完整的,要实现撤销重做功能除了基本的命令模式之外,还要提供撤销和重做两个栈以保存操作命令。

  该模式与面向对象也没很大关系,只是提供了一个实现一些特殊功能的标准或通用方案。

适配器模式

  该模式正如其名,主要用来转换接口,将一个类的方法用其它类封装一下,以达到兼容其它类接口的目的,同时对外可接口保持不变,该模式通过对象组合实现。Shaun 没使用过该模式,就 Shaun 感觉这个模式应该可以用在维护了好几年的系统上,当新作的东西需要兼容老接口时,可以用适配器模式将新街口封装一下。

  该模式同样只是提供了一种新接口兼容老接口的一种优良方案,当然实际使用过程中可能很难这么完美,但算是一种思路。

外观模式

  该模式算是封装的一种体现。当一个功能需要经过多次函数调用才能完成时,这时可以用另一个方法将这些函数都封装起来,从而简化调用方式。Shaun 用该模式处理整个渲染模块的初始化和资源释放,因为初始化时需要分配好很多东西(光照,viewport,固定图层,地面,天空盒等),而释放时同样需要释放这些东西。该模式同样只能算是提供了一种好的编程实践,实际使用过程可能每个函数都有很多参数,调用顺序可能有变,这时简化调用反而没有必要,让用户自己决定怎样调用更好。

  外观模式涉及到的原则有:最少知识原则。该原则主要用来减少对象依赖,即尽量不将类中组合的对象直接暴露出去,而应该将组合对象的方法再简单封装一下,再将封装后的方法暴露出去,以减少另外的类又依赖类中组合对象的现象。该原则可以适当遵守,因为有时直接使用更方便一点,多次封装之后反而显得逻辑混乱,增加系统的复杂度。

模板方法模式

  该模式是抽象的一种体现。首先抽象出一套固定化的流程,流程中每个步骤的具体行为并一致,有些默认,有些可以重写,父类固定流程,子类负责重写流程中每个步骤,这就时模板方法模式。Shaun 没写过完全符合该模式的代码,只是写了个类似该模式的模块,该模块有三个功能(编辑道路节点,编辑车道节点,编辑车道线),做完前两个功能后,发现这里有一套逻辑是通用的,那就是滑过节点高亮,选择节点,出现 gizmo,拖动 gizmo,完成编辑(当然还有选择节点后不拖动 gizmo 等一套 if-else 中间处理状态),于是 Shaun 把这一套流程抽象出来,固化方法,这三个功能都继承该类,方法该重写的重写,不仅减少了代码量,同时整个流程也更清晰了,很快完成了第三个功能。

  模板方法涉及到的原则有:好莱坞原则。即由父类负责函数调用,而子类负责重写被调用的函数,不用管父类的调用逻辑,也最好不要调用父类的函数。该原则用来理清流程很方便,只需要看父类即可,但实际编程过程中可能也会遇到子类不可避免的会调用父类的一些公共函数的情况,Shaun 觉得只要流程没问题的话,调用父类函数也能接受,并不需要严格遵守模式。

迭代器模式

  迭代器,即对遍历进行封装,一般只能顺序执行,提供 next() 方法获取下一个元素,集合容器的遍历方式一般都会用迭代器进行封装。Shaun 在这一个半项目里没写过迭代器,毕竟这是非常底层的一个模式,语言库本身有的数据结构大多自己实现了迭代器,除非需要设计一个新的集合或容器数据结构,才需要提供相应的迭代器。因为 js 没有 SortedMap 数据结构,为了高效分配路网元素 id,Shaun 利用 object 简单实现了一个,提供了相应的 forEach 方法。

  迭代器模式涉及到的原则有:单一责任原则。即一类只做一件事,这个原则对于涉及最最底层的接口很实用,而大多具体类很难只做一件事。迭代器模式对于顺序访问来说还是非常有用的,毕竟使用迭代器的人不需要管底层到底用的什么数据结构,反正可以顺序遍历即可。

组合模式

  组合模式与其说是一种模式,更不如说就是一颗树,只是树的节点都是可供继承的类。在标准的组合模式中,父类中一定会有全部子类的全部函数,即所有子类的函数要么是继承自父类,要么是重写父类函数的,这其实是违背上面单一责任原则的,因为这必然会造成有些子类不需要这么多函数。而从组合模式会存储孩子节点这点来看,和装饰者模式有点类似,只不过装饰者只会存一个孩子,而组合模式可能会存多个,当然两者做的事是不一样,只是实现手法类似而已。Shaun 没写过标准的组合模式,如果只要符合树形模式都可认为是组合模式,那在高精地图编辑器中,所有路网元素都会继承一个父类,而道路中又包含车道簇,车道簇中包含车道,这也算组合模式。在 threejs 中有个 Object3D 的基类,所有渲染对象都会继承该类,该类中又包含若干孩子,threejs 计算 Model 矩阵时就是一层层孩子的 Model 矩阵乘下去,直到最后的孩子,结果就是最后 Shader 中的 Model 矩阵。

状态模式

  状态机的状态转移可以说是程序设计中最麻烦的部分之一了,这部分如果写不好的话,根本没法修改维护,同时会造成 bug 频发。在高精地图编辑器中鼠标操作有两类模式,一种是选择模式,另一种是编辑模式,选择模式又分为点选和框选,而编辑模式就非常多了,针对路网的不同元素,编辑模式的具体实现都不会一样,Shaun 首先使用 RxJS 封装了一个鼠标操作类(左键右键中键移动等),后续的鼠标操作都继承自该类,可以算是状态模式的父类,后续的鼠标操作就针对相应的需求实现相应的方法即可,当然其中鼠标操作自身也存在状态转移(左键到右键,左键到鼠标移动等),这些一般都是针对特定需求来的,所以这些小的状态转移一般在鼠标操作内部实现,但需要支持随时从编辑模式到选择模式,这意味着编辑模式编辑到一半的东西都需要支持随时释放,恢复编辑前的样子,这算是一个麻烦的地方,有时忘了释放就会出现问题。

  状态模式算是为解决状态转移问题提供一种理想的方案,但其具体实现并不一定要和书上一样,Shaun 在用 C++ 实现时就采用另一套方案,状态类是独立的,控制状态转移的代码都在状态机内,而不是书中这种直接在状态类中控制状态机。好处坏处都有,看具体需求,Shaun 的方式就是状态类和状态机是分离的,状态类不需要管状态机怎么实现的,只需要管当前状态的情况,但需要在状态机中管理状态转移,而书中实现方式状态机的状态转移放到状态类中了,也因此状态类需要依赖状态机。


剩下的模式,Shaun 就没直接写代码实践了,因为大多都需要跨模块实现,有的甚至就是个小项目了,所以就简要谈谈 Shaun 的个人理解

代理模式

  主要可以用来做权限控制,在模块与模块之间的调用过程中,有时不想要一个模块可以访问另一个模块的全部内容,这时可以使用代理模式,创建一个中间模块,避免两个模块直接调用,同时进行访问控制。代理模式在如今的互联网时代不可避免的会用到,或直接或间接,往最小的说,对象组合也可用来实现代理模式。

复合模式

  将多种模式组合在一起使用,比如 MVC 模式,这种模式与其说是模式,更不如说就是一种架构,一种开发包含客户端系统的通用架构,当然每一层都会有很多模式进行组合,从而造成具体实现差异非常大。

反模式

  反模式指的是用“不好的解决方案”去解决一个问题,Shaun 主要想谈谈开发反模式,因为这非常常见。有时候一个解决方案好不好要从多个角度进行衡量,比如现有技术,长期短期,上手难度,开发效率,维护难度等角度,当出现一个新问题时,往往意味着就有解决方案有缺陷,这种缺陷可能很容易弥补,更可能很难,当很难解决时,往往要采用全新的解决方案,这时团队对新解决方案可能都不熟,也没有魄力去采用新解决方案,只能去老解决方案继续强行打补丁,直到最后没法维护,白白浪费了大量的人力和时间,这是非常典型的一种反模式。

桥接模式

  将抽象接口和实现分开,并独立派生,以支持抽象和实现的同时改变,并相互独立,可适用在需要跨平台跨不同设备的系统上。

生成器模式

  有点像是模板方法模式和工厂模式的结合版,使用多个步骤创建一个对象,但步骤没有固定顺序,可适用于流程复杂的规划系统中。

责任链模式

  可以认为是模板方法模式的进阶版,只是模板的步骤方法变成了一个个对象,并且支持步骤的增加和删除,以及更换顺序,一旦某个步骤成功执行,则整个链条终止,可适用于消除显式的 if-else,处理键盘或鼠标事件时,需要针对不同按键触发不同操作,这时可以采用该模式,缺点是链条很长时,要写很多类,导致执行啥很不直观。

蝇量模式

  这个模式算是一种优化内存占用的方案,通过牺牲类的独立性来减少内存,更彻底一点就是不创建类,直接用函数调用来处理就行。

解释器模式

  可用来实现简单语法规则语言的解释器或编译器,每个语法规则都由一个类进行解析,然后组合。

中介者模式

  可认为是状态模式和代理模式的结合版,不过各个状态可以是不同类,由中介者控制系统流转,集中控制逻辑,使被控制对象解耦,但可能会造成中介者本身非常复杂。

备忘录模式

  可用于系统(游戏)存档,存储系统中关键对象的运行状态,通常实现的方案一般是序列化/持久化,为了效率考虑,难的是有时需要增量存档。

原型模式

  js 的原型链应该是原型模式的典型,不仅实现了动态扩展实例,更实现了动态扩展对象,即继承。在高精地图编辑器中,由于需要做自动保存,所以在做序列化和反序列化的同时也简单实现了对象的 clone(),即从当前实例中创建一个完全一样的实例,可认为是 C++ 中的深拷贝。

访问者模式

  相当于加个中间层,从而以最小的代价修改现有系统(一般是增加一个方法),达到外部可以取得系统内部信息的目的。

后记

  曾看过这样一句话:抽象能力决定编程能力,Shaun 个人认为,所谓抽象即提炼事物的共同点,这也是设计模式中反复使用接口的原因,接口即一组具体类的共同点,接口中的函数和变量即为这些具体类共有的,虽然具体行为可以不一样,但行为本身总是存在的。而又有这样一句话:程序等于数据结构加算法,Shaun 的理解是,狭义上的程序确实是这样,一段代码解决一个问题,这是程序,多段代码解决一个问题,这也是程序,多段代码解决多个问题,这亦是程序,一个软件,一个系统,一个平台,都是程序,但显然这些程序不是简单的数据结构和算法就能概括的,其内部必然有一套合理的逻辑进行组织,这套逻辑可以称之为“设计模式”,但这个“设计模式”不仅仅是上面谈的这些模式概念。Shaun 认为好的数据结构和算法确实能使程序达到当前最优,但对于一个大型系统或平台来说,这种最优只能是一种局部最优,这对整个系统的全局最优来说可能是微不足道的,而“设计模式”要解决的是怎样使一个系统达到全局最优,怎么合理组织局部最优。面对现代的超大型系统或平台,传统意义上的设计模式也只能达到局部最优,全局最优基本很少有人能驾驭,都是针对特定的业务需要,不断的试错改进优化,逐渐趋于稳定,但这种稳定可能很难抽象,放进其它的业务中,又得花费大量的人力物力去修改。

  Shaun 个人对现代大型系统架构的理解就是分层分模块,功能太多分模块,模块太多就分层,一层不够分两层,两层不够分三层,三层不够继续分,根据数据流的处理分层,根据功能的不同分模块,模块内部依靠设计模式进行组织,设计模式调度的就是数据结构与算法。Shaun 目前的设计原则就是:每层一个独立的服务控制模块,每个模块一个对外服务功能(或事件或 socket ),同层的各模块之间尽量保持独立,不相互依赖,若各模块存在共同点,则将共同点抽出来,将其作为公共模块或独立为小层,层与层之间通过服务控制模块进行数据流的传输,除服务控制模块之外,模块之间尽量避免相互通信,即每个模块的对外服务功能一般只对本层服务控制模块提供服务,最好是只负责接收数据。如果系统实在太大,就只能保持纵向分层,横向保证各模块间数据流依次传输,并在特定模块节点上进行上下层的数据传输。

  数据结构与算法重不重要?当然重要,数据结构与算法不过关,面试都过不去 ( ╯□╰ ),工作都没有,还何谈什么设计模式,什么架构。设计模式重不重要?当然也重要,不会合理使用设计模式,写出一堆别人没法维护的垃圾代码(当然,这或许是好事 :p ),改个需求要半天,加个功能要半个月,效率太低,这样即使有好的数据结构与算法作用也不大。但是设计模式也不是万能的,针对不同的需求,同一种设计模式有不同的实现方式,所以书中的设计模式也仅供参考而已,与其说设计模式重要,还不如说书中那几个设计原则更重要些。同时一味的追求设计模式也不见得是件好事,设计模式可以参考,但不能生搬硬套,毕竟人是活的,需求也是活的,固定的模式也需要有所改变,总而言之,能以最小的代价解决问题完成需求的模式就是好模式。

积分计算

前言

  最近需要计算一下曲线长度,无法直接得到被积函数的原函数,常规的积分解法牛顿莱布尼茨公式没法使用,所以只能使用数值积分计算积分。

前言

  最近需要计算一下曲线长度,无法直接得到被积函数的原函数,常规的积分解法牛顿莱布尼茨公式没法使用,所以只能使用数值积分计算积分。

  下面主要介绍两种数值积分方法:龙贝格积分(Romberg Quadrature) 和 高斯-克朗罗德积分(Gauss-kronrod Quadrature)。 下面附带的代码只做简单测试,不保证正确性,语言使用 Typescript。

Romberg 篇

  计算积分最简单的当然是使用复化梯形公式,即 \(I=\int_a^b{f(x)dx}=\frac{b-a}{2n}[f(a)+f(b)+2\sum\limits_{i=1}^{n-1}f(x_i)]= T_n, x_i=a+i*h, h=(b-a)/n\) ,若将 n 段每段一分为 2,可得到 \(T_{2n}=T_n/2+\frac{b-a}{2n}\sum\limits_{i=0}^{n-1}f(x_{i+1/2})\) 。考虑数列 \(T=\{T_1,T_2,T_{2^2},...,T_{2^k}\}\),显然该数列必收敛,最后收敛为对应积分,当 \(|T_{2^k}-T_{2^{k-1}}| < ε\)\(ε\) 为精度)时,可取 \(T_{2^k}\) 作为最后的积分结果。但是,直接利用梯形公司求解,收敛很慢,导致计算效率很差,所以需要使用理查德森(Richardson)外推法加速收敛,设 \(T_{2^k}^{(m)}\) 为 m 次加速值,当 m 为 0 时,表示没加速,为原梯形公司,则 \(T_{2^k}^{(m)} = \frac{4^m}{4^m-1}T_{2^{k+1}}^{(m-1)}-\frac{1}{4^m-1}T_{2^k}^{(m-1)}\),当 \(|T_{2^{k+1}}^{(m)}-T_{2^k}^{(m)}| < ε\) 时,则收敛,并取其中一值作为最终的积分值。未经修饰的代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
function rombergIntegrator(f: (x: number) => number, a: number, b: number, tolerance = 1e-8) {
let h = b - a;
let n = 1;
let preTK = ((f(a) + f(b)) * h) / 2;

let tkmArray = [];
let m = 0;
tkmArray[m] = preTK;

while (true) {
m += 1;
console.log(m, tkmArray[m - 1]);

let tk = getTK(f, preTK, n, a, h);
if (Math.abs(tk - preTK) < tolerance) return tk;

let preTKM = tkmArray[0];
let preTK1M = tk;
tkmArray[0] = tk;
for (let i = 1; i <= m; i++) {
let newTKM = getTKM(preTK1M, preTKM, i);
preTKM = tkmArray[i];
preTK1M = newTKM;
tkmArray[i] = newTKM;

if (preTKM !== undefined && Math.abs(preTK1M - preTKM) < tolerance) return preTK1M;
}

preTK = tk;
n *= 2;
h /= 2;
}

function getTK(f: (x: number) => number, preTK: number, n: number, a: number, h: number) {
let sum = 0;
for (let i = 0; i < n; i++) {
let x = a + (i + 0.5) * h;
sum += f(x);
}
return (preTK + h * sum) / 2;
}

function getTKM(preTK1M: number, preTKM: number, m = 0) {
let m4 = 1 << (2 * m); // 4 ** m;
return (m4 * preTK1M - preTKM) / (m4 - 1);
}
}

  由于采用闭型积分规则(积分上下限值参与积分计算),所以以上代码不适合计算两端点值被积函数值无限大的情况(如 1/4 圆弧积分等)。而且该方法不合适求取被积函数在积分区间内导数值变化较大(如被积函数在积分下限附近剧烈波动,在积分上限附近不变化等)的积分,因为该方法是均匀分段,这种情况将导致计算量剧增,这时就需要用到下面的自适应积分。

自适应篇

  自适应积分主要包括两类:全局自适应积分和局部自适应积分,通常情况下全局自适应积分的会比局部自适应积分的表现要好,全局自适应积分一般通过二分递归实现,当满足一定条件时,递归终止,即通过二分分别计算两边的积分,若一边满足一定条件,则不继续划分,从而减少计算量。全局自适应积分中比较经典的有基于辛普森(Simpson)公式的自适应算法,普通辛普森积分公式为:\(I=\int_a^b{f(x)dx}=\frac{b-a}{6}[f(a)+4f(m)+f(b)]= S(a,b), m=(a+b)/2\),复化辛普森公式为 \(I=\int_a^b{f(x)dx}=\frac{h}{3}[f(a)+4\sum\limits_{i=1}^{n}f(x_{2i-1})+2\sum\limits_{i=1}^{n-1}f(x_{2i})+f(b)]= S(a,b)\),其中 \(x_i=a+i*h (i=1,2,3,...,2n),h=\frac{b-a}{2n}\),基于辛普森公式的自适应基本原理如下:令 \(S_2(a,b) = S(a,m)+S(m,b)\),m 为 a,b 中值,若 \(|S(a,b) - S_2(a,b)| < 15ε\),则取 \(S_2(a,b)\)\(S_2(a,b)+(S(a,b)-S_2(a,b))/15\) 作为该区间的积分值,否则,将区间二分递归,同时因为误差会累积,所以每次递归都需要将精度提高两倍,即 \(ε = ε/2\),如此最后的精度才能达到最初的精度。具体 ts 代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
function adaptiveSimpsonIntegrator(f: (x: number) => number, a: number, b: number, epsilon = 1e-8) {
let S = complexSimpson(f, a, b);
return adsimp(f, a, b, epsilon, S);

function adsimp(f: (x: number) => number, a: number, b: number, epsilon = 1e-8, S = 0): number {
let m = a + (b - a) / 2;
let LS = complexSimpson(f, a, m);
let RS = complexSimpson(f, m, b);
const S2 = LS + RS;
const tolerance = 15 * epsilon;
let delta = S - S2;
if (Math.abs(delta) < tolerance) return S2 + delta / 15;
else {
let doubleEPS = epsilon / 2;
return adsimp(f, a, m, doubleEPS, LS) + adsimp(f, m, b, doubleEPS, RS);
}
}

function complexSimpson(f: (x: number) => number, a: number, b: number, n = 1) {
const n2 = n * 2;
const h = (b - a) / n2;
let sum = f(a) + f(b);

for (let i = 1; i < n2; i += 2) {
sum += 4 * f(a + i * h);
}
for (let i = 2; i < n2 - 1; i += 2) {
sum += 2 * f(a + i * h);
}

return (sum * h) / 3;
}
}

  在 D.V. Fedorov 写的「Introduction to Numerical Methods with examples in Javascript」一书中介绍了一种全局自适应方法,即分别使用高阶和低阶的权值分别计算积分,两者之间的差值 \(E\) 作为误差估计,设绝对精度为 \(\delta\) ,相对精度为 \(\epsilon\) ,若 \(|E|<\delta+\epsilon*Q\),Q 为高阶权值计算的积分,则取 Q 作为积分值,否则将积分区间二分,同时使 \(\delta/\sqrt{2}\)\(\epsilon\) 保持不变。具体 ts 代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
function recursiveAdaptiveIntegrator(f: (x: number) => number, a: number, b: number, accuracy = 1e-15) {
return recursiveAdaptiveIntegrate(f, a, b, accuracy);

function recursiveAdaptiveIntegrate(
f: (x: number) => number,
a: number,
b: number,
accuracy = 1e-15,
epsilon = Number.EPSILON,
preFValue?: number[],
): number {
const abscissae = [1 / 6, 2 / 6, 4 / 6, 5 / 6];
const highOrderWeights = [2 / 6, 1 / 6, 1 / 6, 2 / 6];
const lowOrderWeights = [1 / 4, 1 / 4, 1 / 4, 1 / 4];
const isRecompute = [1, 0, 0, 1];
const h = b - a;

const fValue: number[] = [];
if (preFValue === undefined) {
abscissae.forEach((abscissa) => {
const x = a + abscissa * h;
fValue.push(f(x));
});
} else {
for (let k = 0, i = 0; i < abscissae.length; i++) {
if (isRecompute[i]) fValue.push(f(a + abscissae[i] * h));
else fValue.push(preFValue[k++]);
}
}

let highResult = 0;
let lowResult = 0;
for (let i = 0; i < highOrderWeights.length; i++) {
highResult += highOrderWeights[i] * fValue[i];
lowResult += lowOrderWeights[i] * fValue[i];
}
highResult *= h;
lowResult *= h;

const tolerance = accuracy + epsilon * Math.abs(highResult);
let errorEstimate = Math.abs(highResult - lowResult) / 3;
if (errorEstimate < tolerance) {
return highResult;
} else {
accuracy /= Math.sqrt(2);
let m = a + h / 2;
let midIndex = Math.trunc(abscissae.length / 2);
let leftFValue = fValue.slice(0, midIndex);
let rightFValue = fValue.slice(midIndex);
return (
recursiveAdaptiveIntegrate(f, a, m, accuracy, epsilon, leftFValue) +
recursiveAdaptiveIntegrate(f, m, b, accuracy, epsilon, rightFValue)
);
}
}
}

  该方法很巧妙的设计了一组插值点,使得当前计算的函数值正好可以被下次迭代所使用,从而提高性能,同时该方法可以得到 1/4 圆弧长,虽然精度只能到小数点后 8 位,至于 Shaun 写的其它测试函数,都能得到理想精度。

Gauss 篇

  高斯求积法是一种多项式插值积分法,同时由于不计算被积函数在区间两个端点处的值,所以高斯积分法采用的开型积分规则,高斯积分法的衍生方法有很多种,下面主要介绍高斯-勒让德(Gauss-Legendre Quadrature)以及其迭代改良的高斯-克朗罗德法。高斯-勒让德积分法的公式为积分的原始形态,即 \(\int_a^bf(x)dx=\sum\limits_{i=1}^{∞}w_if(x_{i})\approx\sum\limits_{i=1}^{n}w_if(x_{i})\) ,只不过 \(x_i \in [-1,1]\),并且 \(x_i\)\(w_i\) 都通过勒让德多项式求出,所以其原则上只能用在积分区间为 [-1, 1] 上的积分,但是可以将积分从任意区间通过简单的变形变换到 [-1, 1] 上,即 \(\int_a^b{f(x)dx} = \frac{b-a}{2}\int_{-1}^1{f(\frac{b-a}{2}t+\frac{b+a}{2})dt}\) ,从而可以将高斯-勒让德方法扩展到任意积分上。由于每个 n 对应的 \(x_i\)\(w_i\) 都可以查表可得,所以具体代码方面就很简单了,以 n = 4,即插值点个数为 4 为例,ts 代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
function gaussLegendreIntegrate(f: (x: number) => number, a: number, b: number, n: 4 | 8 = 4) {
const weightsAndAbscissae = getWeightsAndAbscissae(n);
const weights = weightsAndAbscissae.weights;
const abscissae = weightsAndAbscissae.abscissae;

const halfH = (b - a) / 2;

let sum = 0;
for (let i = 0; i < abscissae.length; i++) {
let xi = halfH * abscissae[i] + a + halfH;
sum += weights[i] * f(xi);
}

return sum * halfH;

function getWeightsAndAbscissae(n: 4 | 8 = 4) {
switch (n) {
case 8:
return {
weights: [
0.362683783378362,
0.362683783378362,
0.3137066458778873,
0.3137066458778873,
0.2223810344533745,
0.2223810344533745,
0.1012285362903763,
0.1012285362903763,
],
abscissae: [
-0.1834346424956498,
0.1834346424956498,
-0.525532409916329,
0.525532409916329,
-0.7966664774136267,
0.7966664774136267,
-0.9602898564975363,
0.9602898564975363,
],
};
break;

case 4:
default:
return {
weights: [0.6521451548625461, 0.6521451548625461, 0.3478548451374538, 0.3478548451374538],
abscissae: [-0.3399810435848563, 0.3399810435848563, -0.8611363115940526, 0.8611363115940526],
};
break;
}
}
}

  若要提高高斯-勒让德积分法的精度,可通过增加插值点或分多个区间进行积分来实现,但是由于没有误差估计,所以还是没法精确控制精度,对与某些被积函数积分精度高,但对于其它被积函数,积分精度却有限,当然可以简单的引入一些常用的误差估计法,但一般需要重新计算积分,导致效率很低,而高斯-克朗罗德法为其引入了一种基于 Kronrod 点的误差估计法,可充分利用现有计算值,从而达到有效控制精度的同时,性能没有太大的损失。设 \(G(f,n)=\sum\limits_{i=1}^{n}w_if(x_{i})\) 为具有 n 个插值点的高斯-勒让德法计算结果,\(GK(f,n) = \sum\limits_{i=1}^{n}w'_if(x_{i})+\sum\limits_{k=n+1}^{2n+1}w'_kf(x_{k})\) 为高斯-克朗罗德法的计算结果,则 \(|GK(f,n)-G(f,n)|\) 可作为误差估计,有了误差估计,最后再使用全局自适应策略,即可得到精度可控的高斯积分结果。具体 ts 代码如下,以 n = 7 为例:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
function gaussKronrodIntegrator(f: (x: number) => number, a: number, b: number, accuracy = 1e-15) {
return recursiveAdaptiveIntegrate(f, a, b, accuracy);

function recursiveAdaptiveIntegrate(f: (x: number) => number, a: number, b: number, accuracy = 1e-12): number {
const gaussAbscissae = [
0.0,
-4.058451513773971669066064120769615e-1,
4.058451513773971669066064120769615e-1,
-7.415311855993944398638647732807884e-1,
7.415311855993944398638647732807884e-1,
-9.491079123427585245261896840478513e-1,
9.491079123427585245261896840478513e-1,
];
const gaussWeights = [
4.179591836734693877551020408163265e-1,
3.818300505051189449503697754889751e-1,
3.818300505051189449503697754889751e-1,
2.797053914892766679014677714237796e-1,
2.797053914892766679014677714237796e-1,
1.29484966168869693270611432679082e-1,
1.29484966168869693270611432679082e-1,
];

const kronrodAbscissae = gaussAbscissae.concat([
-2.077849550078984676006894037732449e-1,
2.077849550078984676006894037732449e-1,
-5.860872354676911302941448382587296e-1,
5.860872354676911302941448382587296e-1,
-8.648644233597690727897127886409262e-1,
8.648644233597690727897127886409262e-1,
-9.914553711208126392068546975263285e-1,
9.914553711208126392068546975263285e-1,
]);

const kronrodWeights = [
2.094821410847278280129991748917143e-1,
1.903505780647854099132564024210137e-1,
1.903505780647854099132564024210137e-1,
1.406532597155259187451895905102379e-1,
1.406532597155259187451895905102379e-1,
6.309209262997855329070066318920429e-2,
6.309209262997855329070066318920429e-2,
2.044329400752988924141619992346491e-1,
2.044329400752988924141619992346491e-1,
1.690047266392679028265834265985503e-1,
1.690047266392679028265834265985503e-1,
1.04790010322250183839876322541518e-1,
1.04790010322250183839876322541518e-1,
2.293532201052922496373200805896959e-2,
2.293532201052922496373200805896959e-2,
];

const halfH = (b - a) / 2;

let guassResult = 0;
let kronrodResult = 0;
for (let i = 0; i < gaussAbscissae.length; i++) {
let xi = halfH * gaussAbscissae[i] + a + halfH;
let yi = f(xi);
guassResult += gaussWeights[i] * yi;
kronrodResult += kronrodWeights[i] * yi;
}
for (let i = gaussAbscissae.length; i < kronrodAbscissae.length; i++) {
let xi = halfH * kronrodAbscissae[i] + a + halfH;
let yi = f(xi);
kronrodResult += kronrodWeights[i] * yi;
}

if (Math.abs(kronrodResult - guassResult) < accuracy / halfH) return kronrodResult * halfH;
else {
let m = a + (b - a) / 2;
accuracy /= 2;
return recursiveAdaptiveIntegrate(f, a, m, accuracy) + recursiveAdaptiveIntegrate(f, m, b, accuracy);
}
}
}

  简单测试了一下,Shaun 这里写的 gaussKronrodIntegrator 方法最大精度只能到 1e-15,到 16 位就报错递归深度太大了,圆的 1/4 弧长也没法算出来,当然这些问题可通过设置最大递归深度以及处理异常值来解决,Shaun 这里就不继续写了。

后记

  数值积分策略非常多,尤其是针对一些特殊的函数,可能只能使用一些特殊的策略才能计算,Shaun 这里只是介绍了一些比较基础常用的积分方法,能解决大部分积分问题,唯一需要注意一点的就是如何追求性能与精度之间的平衡,因为积分常常涉及到迭代求值,通常而言精度越高,迭代越多,求积时,同时也需要注意被积函数的异常值(如无穷大等),这时可能需要拆分或变换积分区间,并且使用开型积分规则的积分方法进行重新计算。

附录

一些常见的积分变换

\[ \int_a^bf(x) = \int_0^{b-a}f(a+t)dt \\ \int_a^b{f(x)dx} = \frac{b-a}{2}\int_{-1}^1{f(\frac{b-a}{2}t+\frac{b+a}{2})dt} \\ \int_{-∞}^{+∞}{f(x)dx} = \int_{-1}^1{f(\frac{t}{1-t^2})\frac{1+t^2}{(1-t^2)^2}dt} \\ \int_{a}^{+∞}{f(x)dx} = \int_{0}^1{f(a + \frac{t}{1-t})\frac{1}{(1-t)^2}dt} \\ \int_{-∞}^{a}{f(x)dx} = \int_{-1}^0{f(a + \frac{t}{1+t})\frac{-1}{(1+t)^2}dt} \]

参考资料

[1] 积分策略

[2] Gaussian Quadrature Weights and Abscissae

[3] Gauss-Kronrod Quadrature

[4] Gauss-Kronrod Quadrature Nodes and Weights

Geometry增量更新

前言

  优化 DrawCall 是图形学性能优化中老生常谈的问题,而针对 DrawCall 优化有很多方案,大致可分为两种:简化(Simplification)和合并(Consolidation),简化是指减少三角形个数,即将精细的模型变为粗糙的模型以及各种三角形剔除方案(视锥体剔除(Frustum Culling),遮挡剔除(Occlusion Culling)等),而合并自然则是将同一中材质下的多个 geometry 合并成一个 geometry。

前言

  优化 DrawCall 是图形学性能优化中老生常谈的问题,而针对 DrawCall 优化有很多方案,大致可分为两种:简化(Simplification)和合并(Consolidation),简化是指减少三角形个数,即将精细的模型变为粗糙的模型以及各种三角形剔除方案(视锥体剔除(Frustum Culling),遮挡剔除(Occlusion Culling)等),而合并自然则是将同一中材质下的多个 geometry 合并成一个 geometry。

需求篇

  最近一段时间一直在做高精地图道路的编辑,道路的编辑涉及到很多东西,这篇仅简单谈谈编辑性能的问题。为了能更直白的显示高精地图,不能简单的只使用点和线,还需要使用面将路面和路面上的一些路面标识精确的还原出来。在可视化道路时,主要有两种方案:1、将每条道路作为一个单独的 Mesh,即单独控制每条道路的渲染,一条道路至少产生一次 DrawCall,这样可以更方便的对每条道路进行编辑,但在渲染时要求更高的性能,而且道路一多,将不可避免的引起卡顿;2、将所有道路作为一个 Mesh,即直接渲染出一整个路网,这样显著降低了 DrawCall 次数,使渲染更流畅,但问题在于每编辑一次道路时,都需要重新三角化(Tessellation)整个路网,而且在选中一条道路时,为可视化选中效果,同样需要重新三角化该道路,并生成相应的 Mesh,导致编辑卡顿,每编辑完一次都可能需要等待一会儿。所以为了平衡渲染性能和编辑效率,需要有一种折中的方案,即对整个路网的 Geometry 能做到快速的分离与合并,在编辑时将受影响的道路分离出来,而在编辑之后,又将全部道路合并一下,提高显示性能。下面就谈谈 Shaun 对这种方案的一些思考。

编辑篇

  Shaun 为平衡渲染性能和编辑效率,想出的一种方案是增量更新 Geometry,即只删除或增加局部的 Geometry,而其它不受影响的 Geometry 保持原样,如此即可达到快速的分离和合并 Geometry。具体做法如下:

  1. 首先将所有道路的三角化结果合并成一个 Geometry,在合并的同时建立好每条道路的顶点索引以及面的索引(js 中可直接使用 object 进行存储);

  2. 在选择时,根据顶点索引和面索引重建一个 Geometry,再基于该 Geometry 构建一个新的 Mesh 以指示选择效果;

  3. 当删除道路时,需要删除面索引对应的所有面,而顶点索引对应的顶点不需要删除,将顶点索引移到一个用来标识该部分顶点已废弃的容器 F 中;

  4. 当新增道路时,需要先从容器 F 中查找是否有合适的地方放置该道路的顶点,若有,则放置在对应地方,并更新容器 F 中对应元素,若没有,则将该道路的顶点放置在 Geometry 顶点数组的最末尾,放置完顶点之后,同样需要建立该道路的顶点索引和面索引。

※注:至于容器 F 使用怎样的数据结构以及其中的元素该怎样排列,针对不同的顶点索引可以有不同的选择;在新增道路时,同样可以不同的策略来决定放置顶点的位置(可参考操作系统内存分配的模式)。

  由于新增道路时,可能会在容器 F 中产生一些永远无法删除的元素,导致顶点数组空闲碎片。 为抵抗顶点碎片以及减少顶点数目,需要对顶点数据进行压缩(Compaction),即移除没有使用的顶点,将后面的顶点前移,在前移顶点的同时,别忘了需要同时修改面中相应顶点的索引以及更新构建好的每条道路的顶点索引。

后记

  Shaun 这里只是提出了一种想法,最终实现起来发现效果也确实能达到基本需求(针对有很多条道路的大地图,显示性能从原来的十几二十帧到现在的 60 帧,同时选择和编辑也没受到影响),虽然在内存上比原来多增加了近 20%(还有优化的余地),但是为了渲染流畅以及编辑舒服,在如今这个内存越来越不值钱的年代,这种牺牲 Shaun 觉得是能接受的。当然或许有更好的方案,但限于 Shaun 目前的认知,只能暂时想到这一方案了,若有大佬有更好的方案,还望不吝赐教 🙏。

Mapbox显示GeoServer地图

前言

  最近做项目需要用到 Mapbox 这个地图可视化框架,以前也没用过,甲方有自己的地图数据,所以得结合 GeoServer 发布一下,简单记录一下流程。

前言

  最近做项目需要用到 Mapbox 这个地图可视化框架,以前也没用过,甲方有自己的地图数据,所以得结合 GeoServer 发布一下,简单记录一下流程。

发布篇

环境准备

  下载 GeoServer 以及同页面下的 Vector Tiles 插件,将插件中所有 jar 包都复制到 GeoServer 中webapps\geoserver\WEB-INF\lib目录下。在 bin 下执行 startup.bat 启动 geoserver,若需要修改端口,可修改 start.ini 文件中的jetty.port=8080,在浏览器中输入 http://localhost:8080/geoserver/web/,geoserver 中默认账号为admin,密码为geoserver,geoserver 中常用的两个坐标系为 EPSG:4326:wgs84坐标,Mercator 投影,EPSG:900913:wgs84 坐标,Web Mercator 投影,即保证投影为正方形,MapBox 中必须使用EPSG:900913坐标系统,900913 和 3857 是一样的坐标系统,在 PostGIS 中对应的 SRID 是 3857。

设置 Tile Caching

  Tile Layers 中可进行图层和图层组的预览,以及切片,在 seed/truncate 中可以设置切片类型以及自动将切片保存到 \data_dir\gwc 目录中。

  Caching Defaults 需要勾选 Enable TMS Service,以及在 Default Tile Image Formats for:中勾选 application/vnd.mapbox-vector-tile,其它默认即可。

Gridsets 设置新的坐标系统。

地图发布

在 数据 栏下:

  1. 点击工作区,添加新的工作区,命名以及填写 URI,勾选默认工作区。
  2. 点击数据存储,添加新的数据存储,选择数据源,以 Directory of spatial files (shapefiles) 为例,在连接参数下点击 浏览,选择shape文件存放目录,DBF 文件字符集选择 UTF-8 或 GBK。注: shape 文件名中不能有中文。
  3. 点击图层,添加新的资源,添加图层,选择上一步的添加的数据存储名称,点击发布或再次发布,进入发布配置界面,勾选 广告则会在 Layer Preview 中显示,一般不需要勾选,点击Compute from native bounds,GeoServer 会自动计算边框和经纬度信息,然后勾选Linear geometries can contain cicular arcs,使线性几何图形包含环形弧,然后保存。重复当前步骤,直到数据存储中所有图层都发布完毕。
  4. 点击图层组,添加新的图层组,添加图层,然后点击 生成边界,保存,即完成整个地图的发布。

在 Layer Preview 中点击 openLayers 进行地图预览,随意点击地图,若出现乱码,则需要在数据存储中修改 DBF 文件字符集。

Mapbox 访问 GeoServer 地图

  点击 Geoserver 的logo,然后点击 TMS 1.0.0协议,页面跳转后,查找需要访问的外部地址,即对应 TileMap 的 href 属性。MapBox 中访问发布好的切片服务需要在 http://localhost:8080/geoserver/gwc/service/tms/1.0.0/MapBoxTest:Test@EPSG:900913@pbf/{z}/{x}/{y}.pbf,即 href 属性值后面加上 /{z}/{x}/{y}.pbf ,同时注意 在 MapBox style 下 layers 中 source-layer 的值必须为图层名,这里为 "Test" (若使用图层组,则需要找到 Test 图层组下面的图层,使用对应图层名)。简单示例代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
const mapStyle: mapboxgl.Style = {
version: 8,
sources: {
geoserverData: {
type: 'vector',
scheme: 'tms',
tiles: ["http://localhost:8080/geoserver/gwc/service/tms/1.0.0/MapBoxTest:Test@EPSG:900913@pbf/{z}/{x}/{y}.pbf"],
},

// 使用 OSM 数据源作为底图
// OsmTiles: {
// type: "raster",
// tiles: ["http://a.tile.openstreetmap.org/{z}/{x}/{y}.png"],
// tileSize: 256,
// },
},

layers: [
// 背景图层
{
id: 'background',
type: 'background',
paint: {
'background-color': "rgb(0, 0, 0)",
},
// interactive: true,
},

// {
// id: "OsmTiles",
// type: "raster",
// source: "OsmTiles",
// "source-layer": "osmtiles",
// },

// 道路
{
id: 'road',
source: 'geoserverData',
'source-layer': 'Test',

type: 'line',
layout: {
'line-cap': 'round',
'line-join': 'round',
},
paint: {
'line-width': {
base: 1.5,
stops: [
[5, 0.75],
[18, 32],
],
},
'line-color': 'rgb(255, 255, 255)',
},

interactive: true,
},

// 地图标注
{
id: 'label',
source: 'geoserverData',
'source-layer': 'Test',

type: 'symbol',
layout: {
'text-size': {
base: 1,
stops: [
[9, 10],
[20, 16],
],
},
'text-max-angle': 30,
'symbol-spacing': 250,
'text-font': ['Microsoft YaHei'], // 标注使用字体
'symbol-placement': 'line',
'text-padding': 1,
'text-rotation-alignment': 'map',
'text-pitch-alignment': 'viewport',
'text-field': '{name}', // 标注显示属性名
'text-letter-spacing': 0.01,
},
paint: {
'text-color': 'hsl(0, 0%, 0%)',
'text-halo-color': 'hsla(0, 0%, 100%, 0.75)',
'text-halo-width': 1,
'text-halo-blur': 1,
},
interactive: true,
},
]
};

const map = new Map({
container: "map-container", // html container id
// style: "mapbox://styles/mapbox/outdoors-v11", //hosted style id
style: mapStyle,
center: [0, 0], // starting position [经度, 纬度]
zoom: 1, // starting zoom
antialias: true,
// maxZoom: 24,
// minZoom: 1,
// pitch: 0,
// maxPitch: 60,
// // minPitch: 0,
crossSourceCollisions: false,
});

GeoServer 跨域问题

  将 lib 目录中的 jetty-servlets 和 jetty-util 两个 jar 包复制到\webapps\geoserver\WEB-INF\lib目录下,将\webapps\geoserver\WEB-INF\web.xml文件中两个 <!-- Uncomment following filter to enable CORS 注释取消,重启 GeoSever。

后记

  Mapbox 显示的地图可以自定义样式,而且加载速度渲染性能方面也都还可以,最重要的是由于采用前端渲染矢量,所以没有传统瓦片那种缩放模糊的感觉,这点非常好,本来想总结一篇简单的 Mapbox 使用手册,但没时间整理了,还是算了 😅。

参考资料

基于geoserver+mapbox的定制化离线地图技术方案

Docker使用小结

前言

  最近由于项目部署需要,简单学习了 docker 的使用和回顾下 CentOS 中的常用操作。

前言

  最近由于项目部署需要,简单学习了 docker 的使用和回顾下 CentOS 中的常用操作。

Docker 篇

  由于 Shaun 此次需要安装的环境有点偏门,没找到有完全符合要求的镜像,同时也趁着这次机会学一下 docker,所以就还是直接从最开始的装起了。

  首先使用 docker pull centos:7 拉取 CentOS7 的系统镜像,使用 docker images 查看已有的本地镜像信息,使用 docker ps -a 查看当前已有的容器信息,去掉参数 a ,即显示正运行的容器,docker stop [container id | name] 可关闭指定容器,docker start [container id | name] 可打开指定容器,docker rm [container id | name] 可删除指定容器,docker rmi [image id | name] 可删除指定镜像,再删除镜像之前需要先删除依赖该镜像的所有容器。

  使用 docker run -dit -p 80:80 -p 8080:8080 --name CentOS7 centos:7 bash 开启一个新的容器,其中参数的意义为: -i: 交互式操作;-t: 终端;-d: 后台启动;-p: 设置主机的端口映射到容器内的端口;-name: 指定容器名称;最后的 bash 代表使用 bash 终端。在 windows 中直接使用 docker run 运行镜像时会出现 the input device is not a TTY. If you are using mintty, try prefixing the command with 'winpty' 的错误,前面加上 winpty 即可,即 winpty docker run ...。当没有参数 d 时,则直接进入容器,而当存在参数 d 时,由于容器实在后台启动,进入容器时需要执行 docker exec -it [container id | name] bash 才能进入容器,而退出容器可以输入 exit 命令。

  而为了在容器中可以开启后台服务,需要在开启容器时就进行提权,在 Windows 中提权命令为 docker run -dit --privileged=true --name CentOS7 centos:7 init;而在 Linux 中提权命令为 docker run -dit --privileged=true --name CentOS7 centos:7 /usr/sbin/init

  在新开启的容器中添加数据和相应的环境之后,即可使用 docker commit CentOS7 new_image:tag 生成一个新的镜像(生成镜像之前最好关闭容器),该镜像包含已经安装的环境和数据,再使用 docker save -o centos7.tar new_image:tag,可将生成的镜像导出成 tar 包,在其他机器中使用 docker load -i centos7.tar 即可导入该 tar 包,并生成相应的镜像, 从而简单便捷的完成环境和数据迁移部署任务。

  容器有时需要和主机之间传输文件,有两种方案,一种是直接采用共享文件夹的方式,设置某个目录为两系统的共享目录,从而实现文件传输;另一种是使用 docker cp 命令,使用 docker cp src_path container:dst_path 将主机的文件拷贝到容器中;使用 docker cp container:src_path dst_path 将容器的文件拷贝到主机中。

CentOS 篇

安装PostgreSQL12

CentOS安装PostgreSQL

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
# Install the repository RPM:
yum install -y https://download.postgresql.org/pub/repos/yum/reporpms/EL-7-x86_64/pgdg-redhat-repo-latest.noarch.rpm

# Install PostgreSQL:
yum install -y postgresql12-server

# Optionally initialize the database and enable automatic start:
/usr/pgsql-12/bin/postgresql-12-setup initdb
systemctl enable postgresql-12
systemctl start postgresql-12

#开启远程访问
修改/var/lib/pgsql/10/data/postgresql.conf文件,
取消 listen_addresses 的注释,将参数值改为“*”

修改/var/lib/pgsql/10/data/pg_hba.conf文件,增加下
# IPv4 local connections:
host all all 127.0.0.1/32 md5
host all all 0.0.0.0/0 md5

# 给数据库postgres用户分配密码 1
psql -U postgres
alter user postgres with encrypted password '1';

# 重启服务
systemctl restart postgresql-12

由于需要切换账户,所以最好在第一步安装完之后就使用 passwd [username] XXXXXX 设置 root 和 postgres 两个账户的密码。

安装postgis

(https://yum.postgresql.org/12/redhat/rhel-7-x86_64/repoview/)

1
2
3
4
5
6
7
8
9
10
11
12
# 解决依赖
yum install epel-release

yum install -y https://yum.postgresql.org/12/redhat/rhel-7-x86_64/postgis30_12-3.0.2-2.rhel7.x86_64.rpm

# 安装pgrouting
yum install -y https://yum.postgresql.org/12/redhat/rhel-7-x86_64/pgrouting_12-3.0.2-1.rhel7.x86_64.rpm

# 安装ogr_fdw
yum search ogr_fdw # 查询
yum install ogr_fdw # 安装

postgres 移除扩展 drop extension postgis cascade;

安装Java

(https://www.cnblogs.com/bcomll/p/12142747.html)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 搜素java
yum search java | grep -i --color JDK

#安装
yum install java-11-openjdk

# 配置环境变量
vi /etc/profile
# 加入内容
export JAVA_HOME=/usr/lib/jvm/java-11-openjdk-11.0.8.10-0.el7_8.x86_64
export CLASSPATH=.:$JAVA_HOME/jre/lib/rt.jar:$JAVA_HOME/jre/lib/dt.jar:$JAVA_HOME/lib/tool.jar
export PATH=$PATH:$JAVA_HOME/bin

# 重启环境
source /etc/profile

安装GeoServer

(https://blog.csdn.net/weixin_34205076/article/details/88734828)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
mv /tmp/geoserver-2.13.2 /usr/share/geoserver

### 添加环境变量
vi /etc/profile
# 追加
export GEOSERVER_HOME=/usr/share/geoserver
# 重新加载/etc/profile文件
source /etc/profile

# 授权
chown -R root:root /usr/share/geoserver

#### 改造启动脚本
vim /usr/share/geoserver/bin/startup.sh

# 在最上面引入环境变量
source /etc/profile

# 最后执行改为nohup,并将日志输入到 var/log/geoserver.log
nohup "$_RUNJAVA" $JAVA_OPTS $MARLIN_ENABLER -DGEOSERVER_DATA_DIR="$GEOSERVER_DATA_DIR" -Djava.awt.headless=true -DSTOP.PORT=8079 -DSTOP.KEY=geoserver -jar start.jar 1>/dev/null 2>/var/log/geoserver.log &

#### 修改停止脚本
vim /usr/share/geoserver/bin/shutdown.sh
# 在最上面引入环境变量
source /etc/profile

#### 创建服务
vim /lib/systemd/system/geoserver.service

[Unit]
Description=geoserver service
After=network.target

[Service]
Type=forking
LimitNOFILE=65536
ExecStart=/usr/share/geoserver/bin/startup.sh
ExecReload=
ExecStop=/usr/share/geoserver/bin/shutdown.sh
Restart=on-abort

[Install]
WantedBy=multi-user.target


## 开机自启
systemctl enable geoserver
## 开启服务
systemctl start geoserver

安装 nginx

1
2
3
yum install nginx # 下载并安装nginx

systemctl start nginx # 启动nginx服务

  在 /etc/nginx 下可修改 nginx.config 文件,监听端口默认是 80,直接输入本地地址可能并打不开网页,因为直接这样安装的 nginx 的可能“有毒”,在 /usr/share/nginx/html 目录中 index.html 可能并不是一个 html 文件,而只是快捷方式,需要从别的地方拷贝一个真正的 index.html 文件替换该文件才可正常打开网页。

安装 nodejs

(https://github.com/nodesource/distributions)

1
2
3
curl -sL https://rpm.nodesource.com/setup_10.x | bash -

yum install -y nodejs

VSCode 远程篇

  需要安装 VSCode 远程开发插件 https://marketplace.visualstudio.com/items?itemName=ms-vscode-remote.vscode-remote-extensionpack,包括(Remote - WSL,Remote - Containers,Remote - SSH)。

在远程资源资源管理器中 切换到 SSH Targes 标签,点击设置设置,在 C:/Users/用户名/.ssh/config 中输入

1
2
3
Host [随便写]
HostName remote-ip 或 域名
User 远程服务器用户名

  配置SSH,通过命令 ssh-keygen -t rsa -b 4096 生成密钥对(在 C:/Users/用户名/.ssh/ 目录中),将公钥内容拷贝到远程服务器/root/.ssh/authorized_keys 中,修改配置文件 /etc/ssh/sshd_config,取消 #PubkeyAuthentication yes 注释,允许使用基于密钥认证的方式登录。重启 sshd 服务 systemctl restart sshd通过 VSCode ssh 远程连接在结束之后最好将终端全部删除,尤其是最开始的那个 install server 终端。

  设置好SSH之后,通过在 VSCode 中设置 "docker.host":"ssh://your-remote-user@your-remote-machine-fqdn-or-ip-here",可以直接连接在远程服务器上的 docker。该设置最好用于 工作区设置,而不是 用户设置。

  VSCode 远程连接 SSH 有时可能会出现无法连接,一直尝试连接的现象,这时需要粗暴的删除 ~/.vscode-server 目录,重新进行连接,不行的话就只能参考: https://stackoverflow.com/questions/56718453/using-remote-ssh-in-vscode-on-a-target-machine-that-only-allows-inbound-ssh-co,先关闭远程服务器上已存在的所有 vscode-server 进程,通过 https://update.code.visualstudio.com/commit:$COMMIT_ID/server-linux-x64/stable 下载 tar 包,使用 tar -xvzf vscode-server-linux-x64.tar.gz --strip-components 1 后再重新连接。

后记

  据查 docker 依然存在很多缺点,尤其是守护进程,Shaun 这两天在 Windows 中使用 Docker 也时不时的出现 docker 卡死的问题,需要重启 docker 服务, 以后有机会还是使用 Podman 吧。RedHat 系的 yum 也快退休了,以后再需要安装软件可能就直接上毒奶粉(bushi)dnf 了。VSCode 远程开发是真的舒服,文件无缝传输,任意修改,可以尽情享受现代编辑器的方便。

空间中三角形与三角形相交

前言

  一种快速判断空间中三角形与三角形相交的方法,出自论文:Tomas Moller. A Fast Triangle-Triangle Intersection Test. Journal of Graphics Tools, 1997, 2(2):25-30. ,与「快速判断三角形与长方体相交」中那篇论文的作者是同一个人。

前言

  一种快速判断空间中三角形与三角形相交的方法,出自论文:Tomas Moller. A Fast Triangle-Triangle Intersection Test. Journal of Graphics Tools, 1997, 2(2):25-30. ,与「快速判断三角形与长方体相交」中那篇论文的作者是同一个人。

相交篇

  该论文的理论基础部分来自分离轴理论,论文中三角形与三角形的相交测试主要可分为 3 类:1、沿三角形所在平面法线方向的相交测试;2、沿三角形所在两平面交线方向的相交测试;3、共面时的相交测试。下面来逐步分析这些相交测试。

沿法线方向相交

  沿法线相交很简单,直接使用分离轴理论分别判断两三角形在对应两条法线上的投影是否相交,若存在一条法线使投影不相交,则可直接判定两三角形不相交,若投影都相交,则存在两种情况,一种是两三角形共面,一种是两三角形相互跨立,判断这两种情况的依据为计算两段投影之间的距离,具体计算方法为:计算三角形 B 上三个点到三角形 A 所在平面的距离,距离计算的方法可参考「计算几何基础—点到平面的距离」,此距离同时需要保留方向,若三个点的距离都为 0,则两三角形共面;若三个距离都同号,则说明投影不相交,即两三角形不相交;若三个距离存在异号现象,则 B 跨立 A。

沿交线方向相交

  若两三角形相互跨立,则需要判断两三角形是否在两三角形所在平面交线上存在相交。由于两三角形相互跨立,则两三角形必然都与交线相交,则需要分别计算两三角形与交线的交点,根据交点判断两交线段是否相交,若相交,则可直接判定两三角形相交,若不相交,则同样可直接判定两三角形不相交。

  问题的关键现在在于求取两交线段,比较粗暴的方式为:先求出两平面的交线,交线的方向向量为两三角形所在平面法向量的的叉积,交线上的一点通过联立两平面方程进行求取,由于是两个方程求 3 个未知数(x,y,z),所以理论上有无数个解,令交线方向向量绝对值最大的分量对应的未知数为 0,消除一个未知数,还剩两个,可得唯一解,即可求出交线上一点,亦可得到直线参数方程,求两交线段相当于求四个交点,根据直线与线段相交可得到交点,进而得到交线段。

论文中的方法为:求出交线的方向向量 D 后,设 O 为交线上任意一点(不需要求),则交线方程为 \(L= O+tD\),此时求交线段只需要求出 4 个 t。先求交线与三角形 A 的交线段,设三角形 A(V0,V1,V2)三个顶点到三角形 B 所在平面的距离分别为 \(d_0,d_1,d_2\),设 \(d_1\) 与其它两个距离异号,则交线分别与 V0V1 和 V2V1 相交,先求与 V0V1 相交时的 t1,设三角形 B 所在的平面为平面 B,V0 和 V1 在平面 B 上的投影分别为为 K0 和 K1,V0 和 V1 在交线上的投影分别为为 P0 和 P1,交线与 V0V1 的交点为 C1,设 \(P0=O+p_0D,P1=O+p_1D,C1=O+t_1D\),则 \(p_0=(V0-O)·D,p_1=(V1-O)·D\),由论文中图二可知,有两组三角形相似,分别为 V0C1K0和V1C1K1 相似,以及 V0C1P0和V1C1P1 相似,所以 \[ \frac{V0K0}{V1K1} = \frac{V0C1}{V1C1} = \frac{C1P0}{C1P1} \Rightarrow \frac{d_0}{d_1} = \frac{p_0-t_1}{p1-t_1} \Rightarrow t_1 = \frac{d_0p_1-d_1p_0}{d_0-d_1} \] 若点 O 为原点在交线上的投影,则 \((V0-O)·D=V0·D=p_0,  p_1=V1·D\),若将交线投影到交线方向向量绝对值最大的分量对应的坐标轴上,在该投影交线上进行相交测试与在原交线上进行相交测试是等效的,所以此时 \(p_0,p_1\) 即为 V0 和 V1 上对应坐标轴的分量。同理可求出 t2,t3 ,t4,两交线段为 t1t2 和 t3t4。

共面相交

  判断共面相交,相当于判断 2 维中两三角形是否相交,先将三角形投影到 XOY,XOZ,YOZ 平面中投影面积最大的平面(为避免计算面积,可直接令三角形顶点坐标中对应法向向量中绝对值最大的分量为 0,即若法向向量中绝对值最大的分量为 y,则将三角形投影到 XOZ 平面),判断投影后的两三角形是否相交即可,因为此时的两三角形相交,当且仅当其投影三角形相交。2 维中两三角形是否相交,论文中方法是先判断两三角形的边是否相交,即相当于判断 9 组线段是否相交,若存在一组相交,则两三角形相交,若都不相交,则需要判断其中一个三角形是否被另一个三角形包含,具体判断方式取三角形 A 中一顶点,若该顶点在三角形 B 内(判断 2 维中点在多边形内的方法同样可参考「计算几何基础」),则说明三角形 A 在 B 内,同样也需要判断 B 是否在 A 内,若都不在,则两三角形不相交,否则两三角形相交。当然,也可以直接通过分离轴理论来判断两三角形是否相交,毕竟,这就是在 2 维中,而且三角形算是天然的凸多边形。

后记

  三角形是图形学中组成面的基本单元,图形学中面与面的碰撞检测都可以很粗暴的用三角形相交来实现,只不过直接一个个判断效率有点低就是了,所以一般会借助一些基于的树和包围盒空间加速结构,或者针对某种形状的特殊方法,这又是新的方法系列了。