前言

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

准备搭的环境为：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

Hadoop

Spark

需开放端口 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 用过的和没用过的一些设计模式。

　　以「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

前言

　　优化 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 发布一下，简单记录一下流程。

发布篇

环境准备

　　下载 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 的使用和回顾下 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. ，与「快速判断三角形与长方体相交」中那篇论文的作者是同一个人。

相交篇

　　该论文的理论基础部分来自分离轴理论，论文中三角形与三角形的相交测试主要可分为 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 维中，而且三角形算是天然的凸多边形。

后记

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

前言

　　由于本篇主要是谈谈基础，所以一些快速运算方法一般不在本篇探讨范围之内，一些特殊的快速手法等后续专门独立开篇再谈。

基础篇

解线性方程组：令 \(A\) 为 \(n×n\) 的矩阵，\(α\) 和 \(β\) 为 n 维的列向量，设 \(A = \begin{bmatrix} \alpha_1 & \alpha_2 & \alpha_3 & ... & \alpha_n \end{bmatrix}\)，对于线性方程组 \(Ax = β\)，初等数学中最常规的就是消元法了，但在线性代数中，有两种解法，一种是等式两边同乘 \(A\) 的逆阵，得到 \(x = A^{-1}β\)；另一种是克莱姆法则，可得 \(x_i = |A_i| / |A|\)，其中 \(|A|\) 为 矩阵 \(A\) 对应的行列式， \(|A_i|\) 为将矩阵 \(A\) 中第 \(i\) 列换成 \(β\) 后对应的行列式，如 \(x_3 = \left| \begin{array}{cccc} \alpha_1 & \alpha_2 & \beta & ... & \alpha_n \end{array} \right| / |A|\)。

向量内积：又称向量点积，向量点乘。设 a，b 为空间中两个 n 维向量 \(a=(x_1,x_2,x_3,...,x_n)\)，\(b=(y_1,y_2,y_3,...,y_n)\)，则 \(a·b=|a|*|b|*cos(α)=\sum\limits_{i=1}^{n}x_iy_i\)，其中\(|a|= \sqrt{(x_1)^2+(x_2)^2+(x_3)^2+...+(x_n)^2},\alpha\) 为两向量之间夹角。向量内积一般用来计算投影（令 \(|b|=1\)，则 \(a \cdot b=|a|cos(α)\)，即为向量 a 在 b 上的投影），和两向量之间角度。

向量外积：又称向量叉积，向量叉乘，\(|a×b|=||a|*|b|*sin(\alpha)|\)。向量外积一般只针对二维向量和三维向量，对于二维向量，\(a×b=a.x*b.y-b.x*a.y\)，可以认为是一个值（其实也是一个向量），对于三维向量 \(a×b=(a.y*b.z-b.y*a.z, a.x*b.z-b.x*a.z, a.x*b.y-b.x*a.y)\)，是一个向量，且该向量一定垂直于 a，b 两向量构成的平面，所以三维向量的外积一般可以用来计算平面的法向量。向量外积的几何意义为两向量构成平行四边形的面积，二维向量外积直接取绝对值即为面积，不取绝对值则可以用来判断两向量构成三角形的点是以顺时针排列（小于 0）的还是逆时针排列（大于 0）的，三维向量外积取所得向量的模即为面积。

向量混合积： 设 \(a\)，\(b\)，\(c\) 为空间中三个向量，则其混合积 \((a,b,c)= (a×b)·c = -(c×b)·a = -(a×c)·b\)，\((a,b,c) = (b,c,a) = (c,a,b)\)。若 \(a\)，\(b\)，\(c\) 都为 3 维的向量，矩阵 \(A = [a, b, c]\)，则 \(|A| = (a, b, c)\) ，其中 \(a,b,c\) 是列向量或行向量都可以，因为 \(|A| = |A^T|\)，该混合积的几何意义为这三个向量组成的平行六面体的体积。

直线参数方程：设 P 为直线上任意一点，若 P1 和 P2 为直线上已知两点，则 \(P=P1 + (P2-P1)*t, t\in[-∞,+∞]\)。当然直线的参数方程有很多种形式，之所以用两点式，是因为两点式除了能表示直线，同样能表示射线（\(t\ge 0\)），也能表示线段（\(0\le t\le 1\)）。

平面方程：设 P 为平面上任意一点，O 为平面上已知一点，n 为平面法向量，则平面方程为：\((P-O)·n=0\) 。

三角形内的点：设 P0，P1，P2 为三角形 3 个顶点，若 P 为三角形内一点，同时可认为 (P1-P0) 和 (P2-P0) 为一组基向量，则 P 满足 \(P=P0+u*(P1-P0)+v*(P2-P0),u \in [0,1],v \in [0,1],u+v\le 1\)。

距离篇

　　一般的距离计算都是向量计算，要么点乘，要么叉乘。

点到直线的距离

　　计算点到直线的距离通用的解法是使用向量叉乘，设 P 为直线外一点，Q1 和 Q2 为直线上两点，则可得到两向量 \(a=P-Q1, b=Q2-Q1\)，则 P 到直线的距离为 \(d=|a×b|/|b|\)，叉乘是面积，而面积又等于底乘高，\(|b|\) 为底，\(d\) 为高。当然点到直线的距离还有其他的一些方法，如 利用公式，利用点积再使用勾股定理，直接利用点积计算直线法向量上的投影等，这些方法都有一定的局限性。

点到线段的距离

　　设 d 为点 P 到线段的距离，表示为点 P 到线段上最近点的距离，设 P1，P2 分别为线段两端点，计算 \(a=(P-P1)·(P2-P1) / |P2-P1|\)，若 \(0≤a≤1\)，则 d 为点 P 到线段所在直线的距离，若 \(a<0\)，则 d 为点 P 到点 P1 的距离，若 \(a>1\)，则 d 为点 P 到点 P2 的距离。

三维中点到三角形的距离

　　先计算点 P 到三角形所在平面的投影点 \(P'\)，若 \(P'\) 在三角形内，则只需要求点 P 到三角形所在平面的距离，否则需要分别求点 P 到三角形三条边的距离（点到线段的距离），取三个距离中最小值即为点到三角形的距离。

点到平面的距离

　　直接利用点乘计算，平面外一点与平面内一点构成的向量到平面法向量上的投影，即为点到平面的距离。设 P 为平面外一点，O 为平面内一点，n 为平面法向量，则点 P 到平面的距离为 \(d=(P-O) \cdot n\)。

直线到直线的距离

　　设 P1 和 P2 为直线 1 上两点，Q1 和 Q2 为直线 2 上两点，则直线 1 与直线 2 之间的距离计算流程为：

1. 先判断两直线是否平行，即 \((P2-P1) = k*(Q2-Q1), k \neq 0\) 有解。若平行，则相当于计算点到直线的距离；若不平行，则进行下一步。
2. 再判断两直线是否共面，即 P1，P2，Q1，Q2 四点共面，先计算 \(n=(P2-P1)×(Q2-Q1)\)，再分别计算 \(n·(Q1-P1), n·(Q2-P1), n·(Q1-P2), n·(Q2-P2)\)，若这四个值都为 0 ，则两直线共面（之所以需要判断 4 个值，是为防出现 3 点共线情况，当然也可以先判断 3 点共线，再取不共线的一点构成向量与 n 做点积进行判断），若两直线共面，则两直线必然相交；若两直线不共面，则两直线间距离为 \(d=(Q1-P1) \cdot n\)。

直线到平面的距离

　　设 P1 和 P2 为直线上两点，n 为平面法向量，则直线与平面之间的距离计算流程为：

1. 先判断直线与平面是否平行，即若 \((P2-P1) \cdot n = 0\) ，则直线与平面平行，则直线到平面的距离为 点 P1 到平面的距离；
2. 若直线与平面不平行，则直线与平面必相交。

相交篇

　　一般的相交求交点都是联立方程组，但有些只需要做相交测试的，可以利用一些特殊方法快速求出来。有些相交直接求很麻烦或很难，可以反向求不相交的情况，而在实际编程中，一般都有很多条件语句，以便快速返回提高效率。

直线与直线相交

　　设直线 1 方程为 \(P=P1 + (P2-P1)*t\)，直线 2 方程为 \(P=Q1 + (Q2-Q1)*u\)，联立两方程得 \(P1 + (P2-P1)*t = Q1 + (Q2-Q1)*u\)，即 \(P1-Q1 = (Q2-Q1)*u - (P2-P1)*t\)，设 \(v_0=Q2-Q1,v_1=P1-P2,v_2=P1-Q1\)，即 \(v_0*u+v_1*t=v_2\)，现在的问题就是解这个方程了，一种是直接把向量分解成单一维度，列方程组；一种是等式两边分别同时点乘 \(v_0\) 和 \(v_1\)，可得 \[ \begin{cases} (v_0 \cdot v_0)*u+(v1 \cdot v_0)*t=v_2 \cdot v_0 , \\ (v_0 \cdot v_1)*u+(v1 \cdot v_1)*t=v_2 \cdot v_1 \end{cases} \] 解得： \[ \begin{cases} u=((v1·v1)(v2·v0)-(v1·v0)(v2·v1)) / ((v0·v0)(v1·v1) - (v0·v1)(v1·v0)) , \\ t = ((v0·v0)(v2·v1)-(v0·v1)(v2·v0)) / ((v0·v0)(v1·v1) - (v0·v1)(v1·v0)) \end{cases} \] 若方程有解，则两直线相交，由于两点式方程可以很简单的转换为线段和射线，所以该方法同样可以判断两线段相交，两射线相交，射线与直线与线段相交等。针对二维直线，还有一种解法是将原方程写成矩阵形式，利用克莱姆法则进行求解，不过总的来说，最终的解都是一种形式。

直线与平面相交

　　求直线与平面相交，直接联立直线方程和平面方程即可，得 \((P1-O+(P2-P1)*t)·n=0\) ，即 \(t=(P1-O)·n/((P1-P2)·n)\)。若 t 有解，则直线与平面相交。同样也可以用该方法判断线段或射线与平面相交。

直线与三角形相交

　　二维中判断直线和三角形相交相当于判断直线和线段相交，而在三维中则同样需要联立直线和三角形方程，设直线方程为 \(P=P1+(P2-P1)*t\)，三角形方程为 \(P=Q1+(Q2-Q1)*u+(Q3-Q1)*v\)， 则联立后方程为 \(P1-Q1=(P1-P2)*t+(Q2-Q1)*u+(Q3-Q1)*v\)，令 \(V_0=P1-Q1,V_1=P1-P2,V_2=Q2-Q1,V_3=Q3-Q1\)，则 \(V_0=V_1*t+V_2*u+V_3*v\)，可以分解向量求解，也可以使用克莱姆法则得： \[ t = \frac{\begin{vmatrix} V_0 & V_2 & V_3 \end{vmatrix}}{\begin{vmatrix} V_1 & V_2 & V_3 \end{vmatrix}} \qquad u = \frac{\begin{vmatrix} V_1 & V_0 & V_3 \end{vmatrix}}{\begin{vmatrix} V_1 & V_2 & V_3 \end{vmatrix}} \qquad v = \frac{\begin{vmatrix} V_1 & V_2 & V_0 \end{vmatrix}}{\begin{vmatrix} V_1 & V_2 & V_3 \end{vmatrix}} \qquad \] 由于三阶行列式也可以用向量混合积来求值，所以 \[ t = \frac{-V_0 × V_3 · V_2 }{V_1 × V_2 · V_3} \qquad u = \frac{V_0 × V_3 · V_1 }{V_1 × V_2 · V_3} \qquad v = \frac{V_1 × V_2 · V_0 }{V_1 × V_2 · V_3} \qquad \] 若方程有解，且满足三角形条件，则相交。

二维中凸多边形与凸多边形相交

　　曾在「快速判断三角形与长方体相交」中写过判断两凸多边形是否相交直接使用分离轴理论（separating axis theorem， AST）即可，简而言之就是，取两多边形任意一条边，计算两多边形在该边法向量上的投影是否相交，若存在一条边，使投影不相交，则两凸多边形不相交。

二维中点在多边形内

　　判断点在多边形内有很多种方法，利用叉乘，面积等方法虽然思想简单粗暴但一般计算量较大，且有一定的局限性，仅限凸多边形，但有一种相对快速且能应对各种简单多边形的方法——射线法，射线法的本质是判断射线与线段相交，即从已知点处引一条沿 X 轴正向的射线，若射线与多边形边的相交条数为奇数，则该点在多边形内，该法的缺陷在于若点在边上则需要单独判断。判断该射线与多边形的边是否相交也比较简单，设射线起点为 P0，多边形边的两个端点分别为 P1 和 P2，则射线与边相交需满足的条件为：1、\(min(P1.y, P2.y)<=P0.y<=max(P1.y, P2.y)\)；2、\(P0.x <= (P2.x-P1.x)*(P0.y-P1.y)/(P2.y-P1.y)+P1.x\)。

圆与三角形相交

　　判断圆与三角形相交，即判断圆心到三角形各边的距离是否小于圆的半径，若存在一条边，使圆心到其的距离小于半径，则圆与三角形相交，该距离计算不是点到直线的距离，而是点到线段的距离。

直线与长方体相交

　　判断直线与长方体相交，首先需要将坐标原点平移至长方体中心，再计算过原点且垂直于直线的法向量 n ，随后计算原点到直线的距离 d，若满足 \(|d| <= h_x|n_x| + h_y|n_y| + h_z|n_z|\) ，则 直线与长方体相交。法向量 n 的求法为：设 P 为直线上一点，l 为直线方向向量，则 \(n=-(P·l)*l+P\)。

球与 AABB 相交

判断球与 AABB 相交，首先需要将坐标原点平移至球心，设平移后的 AABB 最小顶点为 V1，最大顶点为 V2，球半径为 r。最简单粗暴的当然是若原点不在 AABB 内，则直接求原点到 8 个面的距离，取最小值，若最小值比半径大，则不相交。另一种方法是将这个问题转化为求解不等式，若存在一点 \(P(x, y, z)\)，使 P 在球内，同时 P 在 AABB 中，即： \[ \begin{cases} x^2+y^2+z^2 ≤ r^2 , \\ V1.x ≤ x ≤ V2.x , \\ V1.y ≤ y ≤ V2.y , \\ V1.z ≤ z ≤ V2.z \end{cases} \] 若该不等式有解，则球与 AABB 相交，否则不相交。解该不等式也简单，由于是存在而不是任意，所以只需要求 \(min(x^2+y^2+z^2)≤r^2 \Rightarrow min(x^2) + min(y^2)+ min(z^2)≤r^2\)，则若 xyz 分量可以取 0，则对应分量取 0，否则取 \(x = min(|V1.x|, |V2.x|), y = min(|V1.y|, |V2.y|),z = min(|V1.z|, |V2.z|)\)，若满足不等式，则相交。

该方法同样可用来求 OBB 与球相交，只需要先利用坐标系变换将 OBB 转成 AABB。

反射篇

　　令入射向量为 \(I\)，法向量为 \(N\)，反射向量为 \(R\)，入射向量与反射向量构成的平面与镜面交线的方向向量为 \(T\)，这四个向量都为单位向量。首先以 \(I\) 和 \(R\) 组成一个菱形，则 \(N\) 和 \(T\) 则为该菱形对角线的方向向量，若已知 \(I\) 和 \(R\)，则 \(N = (R - I).normalize()\)，\(T = (I + R).normalize()\) ，\(normalize()\) 指向量归一化；若已知 \(I\) 和 \(N\)，则 \(R = I + 2*|I \cdot N|*N\)，\(2*|I \cdot N|*N\) 为菱形 \(N\) 方向的对角线向量。

　　镜面反射公式为 \((r',g',b') = (r,g,b) + (1-r,1-g,1-b)*t\)，当 \(t\) 越大则越接近白色，表现为越亮；漫反射公式为 \((r',g',b') = (r,g,b)*t\)，当 \(t\) 越小，则越接近黑色，表现为越暗。 其中 \(t\in[0,1]\)，\(t\) 为入射向量到平面法向量上的投影，即点积。

后记

　　该篇不出意外的话也会是一个长期支持篇，等以后有碰到其他的一些计算几何知识再持续更新吧。

前言

　　接触并使用高精地图和 OpenDRIVE 已有一年的时间，简要写写 Shaun 对 OpenDRIVE 的一些认知。

基础篇

OpenDRIVE 目前最新的版本的 1.6，下面主要结合 1.4，1.5 和1.6 版本一起看。

　　在 OpenDRIVE 中主要有两种坐标系统，一种是常见的 X/Y/Z 空间坐标系统，另一种是 S/T/H 坐标系统。其中 X/Y/Z 坐标系统常与地理信息的各种坐标系统一起使用，S/T/H 坐标系统是针对 OpenDRIVE 中道路参考线设定的一种局部坐标系统，在 OpenDRIVE 中称前者为 “inertial co-ordinates”，后者为 “track co-ordinates”（1.6 中直接为 Reference Line System）。除此之外，还有个局部坐标系 U/V/Z，该坐标系统系统是相对于 S/T/H 坐标系统平移旋转而来的。

对于旋转，以逆时针为正，heading 是指绕 z/h 轴旋转，pitch 是指绕 y/t 轴旋转，roll 是指绕 x/s 轴旋转。

对于曲率，逆时针延申的曲线曲率为正，顺时针延申的曲线曲率为负。

在 OpenDRIVE 中对于道路和车道的描述，有以下几个重要的概念：

• Reference Line。用来指示一条道路的骨架，是 S/T/H 坐标系统的依据，道路中各车道线需要参考这条线。
• Lanes。用来描述各车道以及各车道所属 Lane Section。
• Lane Offset。其意义在 OpenDRIVE 标准看起来很清除，但实际用起来非常模糊，offset 到底是只能偏移一个车道，或半个车道或多个车道，所属哪个 Lane Section，其意义不明，故下文解析篇将直接忽略该属性。
• Lane Section。可以简单理解为子道路，一个 lane section 中包含多条车道，一条道路包含多个 lane section。一个 lane section 中车道数是一个常数，所以对于存在 m 变 n 车道的一条道路，至少要划分为两个 lane section。
• Superelevation、Crossfall 和 Shape。用来表示路面倾斜程度。Superelevation 是整个路面侧向太高，即路面倾斜程度，Crossfall 在 1.6 中已被废弃，原因为 1.6 完善了 Shape，可以完全取代 CrossFall，甚至能做的更好，Shape 主要用来描述路面两侧车道的倾斜程度，通过三次曲线和线性插值可以精确到车道各点倾斜程度。
• Road Linkage。道路之间的连接关系，通过前继（predecessor）和后继（successor）建立道路之间连接关系，若一条路存在多个前继或后继，则其对应前继或后继应该为 Junction ，即目前的标准中暂不支持道路多前后继，1.5 中只支持车道多前后继，道路多前后继依然不支持。
• Junction。交汇处，通俗意义上的路口，主要包含 incomingRoad 和 connectingRoad。
• Junction Group。可以简单理解为交通环岛。
• Neighbors。相邻关系，和 Road Linkage 类似，一个是前后连接关系，一个侧面相邻关系。
• Surface。车道或道路表面材质，有 OpenCRG 则优先使用，没有则由应用程序自定义字符串。
• ....... 等等。还有一些冷门的元素暂时没用到，就不介绍了。

下面就是真正的解析内容了。

解析篇

　　一些简单的就不介绍了，就介绍解析时需要注意的一些重点元素。

Road Geometry

　　geometry 信息可以说是 OpenDRIVE 中最重要的信息，没有之一。OpenDRIVE 中最重要的元素为 Road，而 Road 中最重要的是 Reference Line，而 geometry 正是用来描述 Reference Line 几何线条形状的信息。

　　首先 geometry 标签中共有 5 个属性，分别为 s （该段 geometry 沿参考线起始位置），x（该段 geometry 在 inertial system 下起始横坐标），y（该段 geometry 在 inertial system 下起始纵坐标），hdg（该段 geometry 在 inertial system 下起始弧度），length（该段 geometry 长度）。其次，geometry 共有 5 种线型，分别为 straight lines（直线），spirals（螺线），arcs（圆弧线），cubic polynomials（三次曲线），parametric cubic polynomials（参数化三次曲线），这 5 种线型分别由 geometry 下 5 种标签控制。解析 geometry 的要点在于：先不用管 geometry 标签中的属性，直接解析对应线型标签，需要满足两点：1、起始点坐标一定为 (0, 0)；2、起始点斜率，即导数也一定为 0。依据这两点正确解析完线型之后，再根据 hdg 旋转线型，根据 (x, y) 将线型平移到正确位置。 下面就具体看看这 5 种线型：

1. line，直线。没有任何属性，直接根据 geometry 中 (x, y) 和 hdg 就可得到直线方程。
2. spiral，螺线。有两个属性 curvStart（起始曲率），curvEnd（结束曲率）。螺线解析的代码已经由 OpenDRIVE 官方直接给出了，里面涉及到的数值计算方法就不详解了，直接看官方提供 API 的输入输出，输入有两个：s（从原点开始，螺线延展的长度），cDot（螺线的曲率关于 s 的一阶导数）；输出有 3 个：x（横坐标），y（纵坐标），t（该点的切线弧度）。解析螺线最大的问题应该就是如何得到这两个输入参数，由螺线的性质可以得到，螺线的曲率一定随着螺线的长度均匀变化的，换句话说，对于一条已知螺线，cDot 一定是常数。则 \(cDot = (curvEnd - curvStart) / length\)，其中 length 为 geometry 中的长度属性，下一步需要求出 s，由于已知螺线在原点处的曲率一定为 0，则 (curv - 0) / (s - 0) = cDot，即 s = curv / cDot。由此可得到螺线的各点坐标和切向方向，但由于解析线型需要满足上面说的两点，所以需要将螺线坐标以起始点进行平移和旋转，以满足起始点为原点，起始点切线弧度为 0，最后再将完成平移和旋转后的点根据 geometry 的属性进行旋转和平移以得到真正的坐标点。
3. arc，圆弧线。只有一个属性 curvature（曲率），可根据曲率直接得到圆的半径，在根据曲率的正负可得到该曲线是以顺时针延申还是逆时针延申，再根据解析线型需要满足的两点和 geometry 的属性可得到真正的 inertial system 中的点。
4. poly3，三次曲线，在 1.6 中已被废弃。精度要求低一点可直接插值计算，要求高一点则需要利用 length ，数值计算和二分法直接求出最大的 u，然后插值。
5. paramPoly3，参数化三次曲线。最有名的两个参数化三次曲线就是 Bezier 曲线和 Hermite 曲线，为满足解析线型需要满足的两点，一定有 \(a_u=0, a_v=0,b_v=0\)，需要注意的是一般参数化曲线的参数取值范围为 \([0,1]\)，即该标签最后一个属性 pRange="normalized"，若碰到特殊情况 pRange="arcLength"，则需要将 \(a_u,b_u,\dots\) 等属性转化为参数取值范围为 \([0,1]\) 时对应的属性。

最重要的元素 Geometry 的解析就是这样了，下面谈谈 Shaun 对基于三次曲线的一些元素的理解。

　　比较重要的基于三次曲线的元素主要有 elevation（控制路面的高度），superelevation（道路侧面抬高弧度），crossfall（路面两侧弧度，已被废弃），shape（路面两侧形状，特殊，下文详细介绍），laneOffset（车道偏移量），border 和 width（特殊，下文详细介绍），这些元素所使用的三次曲线一般都基于道路参考线（特殊除外），即该三次曲线的横坐标一般为 s，纵坐标即三次曲线的值则为各元素的信息，这些三次曲线一般是分段描述的，即道路参考线上两个关键点的 s 之间必会生成对应的一段三次曲线，每一段三次曲线的横坐标取值范围都为 \([0,poly3Length]\)，其中 ploy3Length 为该段三次曲线的长度。这些三次曲线段全部合起来则构成沿道路参考线的一条三次曲线，根据 s 计算三次曲线的取值得到对应的信息。

LaneOffset

　　laneOffset 用来指示的是所有车道沿参考线法线方向的偏移量，所有车道包括中心车道（centerLane），由于中心车道没有宽度，所以也可以叫道路中心线，这和道路参考线是两个概念，若 laneOffset 都为 0，则道路中心线和道路参考线重合。

Shape

　　路面两侧车道的高度，该元素虽然也是使用三次曲线进行描述的，但是该三次曲线的横坐标为 S/T/H 坐标中的 t，不是 s，不同的 t 得到不同的高度，该元素下可能存在多段三次曲线，这些三次曲线段是独立的，只是描述其属性 s 所在位置横截面的 shape，而没有完全指定 s 的横截面的 shape，则由两临近 s 的 三次曲线计算相同的 t 对应的高度然后线型插值得到（1.6 标准中插值的公式 Shaun 觉得有问题）。

border 和 width

　　之所以把这个两个元素放在一起写，是因为这两个元素描述的本质上是一个东西，都是车道边界所在的位置。同时，Shaun 还是想吐槽一下，作为一个既定的标准，完全不应该把这两个元素同时放出来，只需要放出一个即可，既然为标准就应该做到完全统一，至于具体用哪个是应用程序的事，但是出来的东西必须得唯一。border 是指车道边界到道路中心线在道路参考线法线上的投影，有正负之分，一般左边车道为正，右边车道为负，而 width 是指车道的宽度，这两者完全可以相互转换。虽然描述这两者的三次曲线是基于道路参考线 s 的，但是这个 s 是相对于 laneSection 标签中 s 属性的，即其真正沿道路参考线的 s' 应该为 s' = s + laneSection.s。

　　至于Junction 的解析好像没什么需要注意的，就不写了，至于 1.5 中的 Virtual Junction，引入新的道路前后继描述方式，也新加了一个 virtual connection，解析时到也没有需要注意的，唯一需要注意是在应用程序中该如何利用这些信息。

　　对于 Signal 和 Object，Shaun 觉得标准中强制规定 s > 0 同样是一件非常不合理的事情，既然能超出道路长度，只准正向超出，不能反向超出，这有点不讲道理。 Signal 和 Object 同样需要注意在应用程序中的用法，至于解析方面，需要注意的应该就是 Object 中 outline 下面的 cornerRoad 和 cornerLocal。

cornerRoad 和 cornerLocal

　　这两个元素描述的本质上也是一个东西，都是面域对象多边形边界上的轮廓点，虽然 1.5 中引入了复杂多边形的概念，即有内轮廓和外轮廓，但 cornerRoad 和 cornerLocal 还是一样的解析。cornerRoad 是直接相对于道路参考线的坐标，即将 X/Y/Z 坐标直接转换为 S/T/H 坐标得到的，所以直接解析转换就可得到该点的真正坐标。而 cornerLocal 的解析则相对要麻烦些，首先根据 object 标签中 s 和 t 属性计算得到 Object 的位置，以道路参考线上 s 所在的位置建立 S/T/H 坐标系统，将该坐标系统平移到 Object 所在位置，即该 S/T/H 系统以Object 的位置为原点，将 cornerLocal 上 u,v,z 属性带入该坐标系统计算出真正的坐标。由于 Object 中的 roll,pitch,hdg 属性，所以还需要对坐标以 Object 的位置为中心进行相应旋转才能得到真正 X/Y/Z 坐标系统中坐标。

好了，以上就是 Shaun 觉得在解析 OpenDRIVE 时需要注意的一些地方了。

后记

　　整个地图行业本来就存在很强的壁垒，国内更是如此，而高精地图作为专为自动驾驶服务的地图，国内估计更是少有人接触过。

　　OpenDRIVE 虽作为一种高精地图标准，但离那种被广泛认可关注的标准还有很长一段距离，虽然它发展了十几年，但奈何整个行业才算是起步阶段，所以之前其一直发展的十分缓慢，之前用 OpenDRIVE 较多的应该是交通仿真领域，这个领域同样存在很强的壁垒，要在这种领域得到整个行业的认可是一件非常困难的事，因为其本来就存在一套自己的格式，而想要应用别人的标准，就势必需要逐渐抛弃自己的格式，这对企业来说需要下很大的决心。

　　就 Shaun 目前所知，行业内虽然有许多企业使用 OpenDRIVE 标准，但基本都是自己魔改后的标准，原因在于 OpenDRIVE 标准还没做到真正成为标准的地步，虽然其一些基本元素具备，但很有很多元素是完全缺失或不完善的，更重要的是，其可视化程序竟然不开源，在整个计算机领域内，还从没见过一种标准没有其开源实现的东西，没有相应的源码，只靠文字来理解难免会产生歧义，各大厂商自然就会选择实现自己的 OpenDRIVE，很难统一。历史上一些经典的论文和算法，都有其相应的开源实现，没有开源实现的论文和算法基本都淹没在了历史长河之中。所以没有开源实现，在计算机领域内很难真正推广开，就很难成为真正的标准。综上，OpenDRIVE 的发展任重而道远啊！希望现在在 ASAM 手中能发展的更快更完善些吧。

前言

　　前段时间做了一件很有意思的事，使用一张普通纹理图片做成了一个天空盒（SkyBox），其实是半个，因为最终形成的天空盒是上下对称的，看起来有天空之境的感觉。

天空盒篇

　　常见的天空盒一般采用立方体包围盒做 geometry，再使用 CubeMap 将 6 张纹理图贴到立方体包围盒上，但是如果只用一张纹理图片，要映射到立方体包围盒上，无论怎么展 UV，都会造成部分棱和顶点 UV 聚集或突变，从而在这样的地方造成贴图扭曲或错位的不协调现象。

　　既然立方体包围盒不可行，就只有采用另一种 geometry —— 球 了，但是直接用球作为 geometry 会造成两种问题：1、球的极点处会有扭曲现象，因为球极点处的 UV 是聚集的，同一个位置，V 都为 1 或 0，而 U 则为 [0, 1] 都有，即在该位置处一个 V 对应多个 U，自然会造成扭曲；2、由于球是一个圆周，所以 U 为 0 的位置和 U 为 1 的位置是同一个，对于一张普通纹理图，这就会造成错位，能看到非常明显的接缝线。所以需要重新展开球面的 UV，而常见的展 UV 方法主要有以下几种：

1. 对于平面，uv 坐标自然通过线性归一化解决，取左下角坐标和右上角坐标，得到平面的宽度和长度，再用当前坐标减去左下角坐标，最后除以平面宽度和长度得到 uv 坐标。
2. 对于球面，自然是取水平方位角 \(θ, \theta \in [-\pi, \pi]\) 和 垂直仰角 \(φ, φ\in [0, \pi]\) 作为 uv 坐标。对于 threejs 而言，方位角在 XZ 平面上，可计算 \(\theta = Math.atan2(z, x)\)，仰角 \(φ = Math.acos(-y)\)，其中 \(x,y,z\) 为当前坐标归一化后的值，则 \(u = (θ +\pi)/(2*\pi)\)，\(v=φ/\pi\)。
3. 对于一般的曲面，可能会依据当前顶点的法线计算 uv 坐标。对法线向量 x,y,z 三个分量的值进行排序，取最小的两个作为 uv，或者将最小的两个除以最大的一个作为 uv（CubeMap 中通常用这种进行 6 个面的纹理映射，最大的那个决定映射 6 个面中的哪一面），或者取固定的两个分量作为 uv ，或者得到最小的两个分量索引取对应顶点坐标的分量值作为 uv 等等，若 uv 取值为 [-1,1]，则还需要做 \(uv = uv *0.5+0.5\) 。根据法线计算 uv 一般需要根据不同的使用情况选择不同的 uv 计算策略。

　　Shaun 很快的解决了第二个问题，想要将一张普通纹理图片贴在整个球面上而不留下接缝线是不现实的，而铺在半个球面上，另一半球面做对称，这个是可行的，而且看起来的效果也不错，计算 uv 就很简单了，原始的 v 不需要改变，只需要将 u 从 [0, 0.5] 映射为 [0, 1]，从 [0.5, 1] 映射为 [1, 0]，即 \[ f(u)=\begin{cases} 2*u,　0 \le u \le 0.5 \\ 2*(1-u),　0.5 \le u \le 1 \end{cases} \] ，接缝线问题使用对称性巧妙的解决了，但第一个问题，极点扭曲现象，还是没法解决，接下来才是真正的难点 ╯︿╰。

　　为了解决极点扭曲问题，首先需要知道极点扭曲的根本原因是同一个位置对应了多个 uv，所以要么将这些 uv 给散开，要么抛弃一部分 uv，散开 uv 的方式 Shaun 没想出来，抛弃一部分 uv 倒是有一种简单的方式，不过抛弃 uv 也意味着会损失一部分纹理，就这个天空盒而言，抛弃一部分 uv 是能接受的，不然平面到球面必然有形变。具体抛弃方法为，将球面铺平，变为一个大圆面，即忽略球面顶点坐标的 z 值，求出球的 AABB，将 AABB 铺平，即为大圆的外接正方形，用计算平面 uv 的方式重新球上各顶点的 uv 坐标，这种方式的确能解决极点扭曲的问题，但又带来了一个更为严重的问题。

　　这个问题就是，球面上半部分显示很正常，但是越靠近底部大圆的部分，纹理拉伸的越厉害，造成天空盒四周都出现很严重的纹理拉伸现象。出现纹理拉伸现象的原因也很好理解，那就是越靠近底部大圆的顶点，uv 坐标之间的间隔越小，即同样的 uv 间隔，顶点之间的距离变大了，在进行纹理插值时，自然会导致纹理拉伸现象。知道问题出现的原因了，那怎么解决了？这又是一个新问题 😔。

　　导致纹理拉伸现象的原因是线性映射，那能不能对计算好的 uv 进行非线性映射，从而抵消 OpenGL 线性纹理插值的影响，非线性映射最重要的是找到合适的非线性函数，常见非线性函数一般有幂函数（伽马变换就是一种幂函数变换，幂 < 1 拉伸小值，幂 > 1 拉伸大值），对数函数，指数等，对于 Shaun 这里的情况，常规的非线性函数肯定是不行的，只能自己想一个函数。

　　注意到，这是在球面上，要取一个非线性映射函数，自然需要从圆的弧长入手，先计算圆的弧长。计算弧长的本质是勾股定理 \(ds=\sqrt{(dx)^2 + (dy)^2} = \sqrt{1+(dy/dx)^2} * dx\)，对于圆 \(x^2 + y^2 = r^2,　y \ge 0\)，即 \(y=\sqrt{r^2-x^2}\)，有 \(dy/dx = -x/ \sqrt{r^2-x^2}\)，则对于圆的弧长 \(L=\int \sqrt{1+x^2/(r^2-x^2)}dx = r \int 1/ \sqrt{1-(x/r)^2}d(x/r) = r*arcsin(x/r)|\)，即对于圆在第一象限的弧长可以为 \(L=r*arcsin(x/r),　0\le x \le r\) 。

　　由于 v 是均匀的，所以只需要对 u 进行拉伸即可，离圆心越近的点，则越接近 0.5， 离圆心越远的点则越偏离 0.5，非线性拉伸公式为 $ u = α * (u-0.5) + 0.5$，其中 \(α = L / (\pi*r/2),　L中的x为点到圆心的距离\)，使用这种拉伸后，天空盒四周的纹理拉伸现象确实不见了，但是又引入了新的问题，天空盒顶部出现了局部拉伸现象，出现微弱的纹理模糊，虽然区域不大，依靠一些手段可以让用户看不到这块区域，但是不能自欺欺人，这一块存在总让 Shaun 觉得很不舒服，于是，就有了下面的终极解决方案。

　　Shaun 最终想出解决方案是：既然单纯的拉伸不能完美解决问题，那还是只能从问题根源入手，完全重新计算 uv 坐标，这次计算 uv 坐标，还是需要借助上文的 \(α\)。具体计算 uv 坐标的方式如下：

1. 先计算球面上顶点相对大圆圆心的角度，即 \(angle = Math.atan2(y - center.y, x - center.x)\)，其中 center 即为大圆圆心（0, 0）。
2. 根据顶点到圆心的距离得到 \(\alpha\)，计算 \(α_x = α * cos(angle)， α_y=α * sin(angle)\) 。
3. 计算 uv：\(u = α_x * 0.5 + 0.5， v=α_y*0.5+0.5\) 。

　　使用这个方式计算 uv，天空盒的全部问题都解决了，天空盒没有任何拉伸扭曲等令人看起来不协调的地方，至于对称也说的过去，Shaun 个人感觉挺漂亮的 (￣▽￣)"。自此天空盒的事就算是告一段落了。

后记

　　做这个天空盒，确实花费了 Shaun 了不少力气，在做的那两天，满脑子都是为什么会拉伸扭曲，以及如何解决拉伸扭曲，最终想出了这套方案，简单优雅，最后的效果也是完美达到了 Shaun 的预期。不过一般人应该也用不到 Shaun 这套方案，这只是 Shaun 自己想做做而已，搞不出来没关系，搞出来当然是好的。