https://github.com/xieguigang/sciBASIC

等高线指的是地形图上高程相等的相邻各点所连成的闭合曲线。把地面上海拔高度相同的点连成的闭合曲线，并垂直投影到一个水平面上，并按比例缩绘在图纸上，就得到等高线。

等高线计算方法

等高线可以通过一种名称为Marching Squares的算法来实现。Marching squares是计算机图形学中一个用于矩阵网格点数据生成等值面的一个算法。在Marching squares算法之中，其主要思想是将每个2*2网格，其中4个顶点网格点根据等值线阈值处理为0，1标量。

• 1 代表数值大于阈值
• 0 代表数值小于阈值

算法主要流程为：

1. 遍历每个小网格，从16中固定情形中选择划线类别
2. 利用线性插值，结合网格点数值找寻交点
3. 划线

假设在一个矩形范围内，有一些离散的高度数据，形如(x, y, height)这样的数据集，通过等高线的绘制方式，我们可以得到一张这样的等高线图。

在Marching squares算法中，square表示的是方块（格子、网格）。因此，我们会首先要把平面分成很多小的网格。

Friend Function padData(data As Double()(), levels As Double()) As Double()()
    Dim rows = data.Length
    Dim cols = data(0).Length

    ' superMin is a value guaranteed to never be included in any isoline.
    ' The idea is to surround the data with low values, forcing polygon
    ' sides to be created.
    Dim superMin = levels(0)

    For i = 1 To levels.Length - 1
        superMin = stdNum.Min(superMin, levels(i))
    Next

    superMin -= 1

    Dim padded = MAT(Of Double)(rows + 2, cols + 2)

    For i = 0 To cols + 2 - 1
        padded(0)(i) = superMin
        padded(rows + 1)(i) = superMin
    Next

    For i = 0 To rows + 2 - 1
        padded(i)(0) = superMin
        padded(i)(cols + 1) = superMin
    Next

    For i = 0 To rows - 1
        For j = 0 To cols - 1
            padded(i + 1)(j + 1) = data(i)(j)
        Next
    Next

    Return padded
End Function

在我们所划分的网格中，每一个顶点都有一个插值过后的数值。这些数据值就是高度值，假设我们现在要画一条高度值为某一个threshold的等高线。我们可以把threshold值和网格上面的数值进行比较，大于等于该threshold值的网格顶点记为1，用实心黑点标记；反之，网格顶点为0，不标记。

Friend Function makeContour(data As Double()(), level As Double) As IsoCell()()
    ' Pad data to guarantee iso GeneralPaths will be closed shapes.
    Dim numRows = data.Length
    Dim numCols = data(0).Length

    ' Create array indicating iso cell neighbor info.
    Dim contours As IsoCell()() = MAT(Of IsoCell)(numRows - 1, numCols - 1)

    For r = 0 To numRows - 1 - 1
        For c = 0 To numCols - 1 - 1
            ' Determine if neighbors are above or below threshold.
            Dim ll = If(data(r)(c) > level, 0, 1)
            Dim lr = If(data(r)(c + 1) > level, 0, 2)
            Dim ur = If(data(r + 1)(c + 1) > level, 0, 4)
            Dim ul = If(data(r + 1)(c) > level, 0, 8)
            Dim nInfo = ll Or lr Or ur Or ul
            Dim isFlipped = False

            ' Deal with ambiguous cases.
            If nInfo = 5 OrElse nInfo = 10 Then
                ' Calculate average of 4 surrounding values.
                Dim center = (data(r)(c) + data(r)(c + 1) + data(r + 1)(c + 1) + data(r + 1)(c)) / 4

                If nInfo = 5 AndAlso center < level Then
                    isFlipped = True
                ElseIf nInfo = 10 AndAlso center < level Then
                    isFlipped = True
                End If
            End If

            ' Store neighbor info as one number.
            contours(r)(c) = New IsoCell()
            contours(r)(c).flipped = isFlipped
            contours(r)(c).neighborInfo = nInfo
        Next
    Next

    Return contours
End Function

进行完顶点的标记之后，我们最后只需要进行划线，产生一个闭合的多边形，就可以得到某一个高度值对应的等高线多边形了。在我们的网格中相邻的四个顶点就是一个cell，对于每一个cell，根据顶点的标记结果，一共有2^4 = 16种情况。

Private Sub interpolateCrossing(isoData As IsoCell()(), data As Double()(), r As Integer, c As Integer, threshold As Double)
    Dim a, b As Double
    Dim cell = isoData(r)(c)
    Dim ll = data(r)(c)
    Dim lr = data(r)(c + 1)
    Dim ul = data(r + 1)(c)
    Dim ur = data(r + 1)(c + 1)

    ' Left side of iso cell.
    Select Case cell.neighborInfo
        Case 1, 3, 5, 7, 8, 10, 12, 14
            a = ll
            b = ul
            cell.left = (threshold - a) / (b - a) ' frac from LL
        Case Else
    End Select

    ' Bottom side of iso cell.
    Select Case cell.neighborInfo
        Case 1, 2, 5, 6, 9, 10, 13, 14
            a = ll
            b = lr
            cell.bottom = (threshold - a) / (b - a) ' frac from LL
        Case Else
    End Select

    ' Top side of iso cell.
    Select Case cell.neighborInfo
        Case 4, 5, 6, 7, 8, 9, 10, 11
            a = ul
            b = ur
            cell.top = (threshold - a) / (b - a) ' frac from UL
        Case Else
    End Select

    ' Right side of iso cell.
    Select Case cell.neighborInfo
        Case 2, 3, 4, 5, 10, 11, 12, 13
            a = lr
            b = ur
            cell.right = (threshold - a) / (b - a) ' frac from LR
        Case Else
    End Select
End Sub

这16种情况对应于轮廓线的绘制也有16种情况。这里在每个cell的每条边都插入了Midpoint。通过连接cell的中点，就可以绘制出轮廓线多边形，这些绘制出来的多边形就是我们的等高线。

最后呢，我们只需要从低高度值按照一定的步长一层层的往高处叠加绘制我们的等高线多边形，就可以得到一幅比较完整的等高线图的绘制结果。

Public Overridable Function mkIsos(data As Double()(), levels As Double()) As GeneralPath()
    ' Pad data to guarantee iso GeneralPaths will be closed shapes.
    Dim dataP = padData(data, levels)
    Dim isos = New GeneralPath(levels.Length - 1) {}

    For i = 0 To levels.Length - 1
        ' Create contour for this level using Marching Squares algorithm.
        Dim contour = makeContour(dataP, levels(i))
        ' Convert contour to GeneralPath.
        isos(i) = mkIso(contour, dataP, levels(i))
    Next

    Return isos
End Function

在这里我们所实现的等高线图的绘制算法，实际上就是在不同的高度水平上将基于Marching Squares算法得到的特征区域的多边形轮廓堆叠在一个二维平面上。如果我们将高度信息给忽略掉，只实现一个高度水平的轮廓计算，理论上我们是可以实现类似于Contour Tracing算法相似的多边形轮廓获取的效果。

如果我们将Marching Squares算法扩展到三围空间，cell变成了cube。这个时候，对应的算法就是Marching cubes这个名字。Marching cubes算法常用于生成一些医学图像。

mzkit中使用等高线可视化代谢物分布

在一级质谱数据之中，一个数据点是由三种维度的值所构成的：m/z，rt，intensity。假若我们将m/z值当作为Y轴，rt值当作为X轴，这个时候我们就可以确定一个坐标点的。再加上对应的intensity响应度值，那么我们就可以通过这三个维度的数据来进行等高线的绘制了。在mzkit程序之中，大家可以通过原始数据文件的概览功能进行一级扫描原始数据的等高线图的查看。

整个过程非常的简单，就只需要在原始数据文件上，通过鼠标右键，选择等高线图绘制菜单就可以完成了。一个原始数据文件的一级质谱扫描的等高线图的可视化效果是和下面类似的图像。

质谱数据一级扫描结果的等高线图

大家对默认的等高线图的样式不太满意的话，可以在图片上点击右键。通过出现的参数调整窗口进行一些可视化样式的调整。

R#语言中进行等高线图的绘制

我们在R#语言之中进行等高线的绘制，只需要将矩阵对象准备好，然后调用对应的绘图函数就好了。例如

imports "charts" from "graphics";

bitmap(file = `${dirname(@script)}/1_Contour.png`) {
    `${dirname(@script)}/1_Contour.csv`
    |> read.csv(row.names = NULL)
    |> contourPlot(colorSet = "Jet")
    ;
}

contourPlot函数是在R#语言之中进行等高线绘图的函数。在这个函数之中，其接受一个数据框对象。对于数据框的数据结构，大家可以通过下面的代码进行查看：

> head(read.csv("1_Contour.csv"))
#         x          y          data
# <mode>  <integer>  <integer>  <double>
# [1, ]   231        21          17.471372
# [2, ]   252        21          16.2102122
# [3, ]   196        28          16.0810415
# [4, ]   217        35          18.067797
# [5, ]   336        35          18.6838193
# [6, ]   490        49          18.1078316

> str(read.csv("1_Contour.csv"))
# 'data.frame':   1807 obs. of  3 variables:
#  $ x    : int [1:1807] 231 252 196 217 336 490 ...
#  $ y    : int [1:1807] 21 21 28 35 35 49 ...
#  $ data : num [1:1807]  17.471372  16.2102122  16.0810415  18.067797  18.6838193  18.1078316 ...

可以看到，我们只需要在数据框之中定义好坐标点，以及坐标点对应的数值就好了。

偏红色的表示对应的数值较高，为等高线的山峰；偏蓝色的区域表示对应的数值较低，为等高线的山谷

• Author
• Recent Posts

谢桂纲

高级数据科学家 at 苏州帕诺米克

Working on Engineered bacteria CAD design on its genome from scratch. Writing scientific computing software for Tianhe & Sunway TaihuLight supercomputer. Do scientific computing programming in R/R# language, he is also the programming language designer of the R# language on the .NET runtime.

Latest posts by 谢桂纲 (see all)

• 机器学习驱动的生物标志物发现与疾病预测集成工具包 - 2025年10月7日
• CCL对象检测算法 - 2025年8月11日
• Boids鸟群模拟算法讲解 - 2025年8月10日