文章阅读目录大纲
https://github.com/xieguigang/sciBASIC
等高线指的是地形图上高程相等的相邻各点所连成的闭合曲线。把地面上海拔高度相同的点连成的闭合曲线，并垂直投影到一个水平面上，并按比例缩绘在图纸上，就得到等高线。
等高线计算方法
等高线可以通过一种名称为Marching Squares的算法来实现。Marching squares是计算机图形学中一个用于矩阵网格点数据生成等值面的一个算法。在Marching squares算法之中，其主要思想是将每个2*2网格，其中4个顶点网格点根据等值线阈值处理为0，1标量。
10
算法主要流程为：
- 遍历每个小网格，从16中固定情形中选择划线类别
- 利用线性插值，结合网格点数值找寻交点
- 划线
假设在一个矩形范围内，有一些离散的高度数据，形如
这样的数据集，通过等高线的绘制方式，我们可以得到一张这样的等高线图。(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
> 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 ...可以看到，我们只需要在数据框之中定义好坐标点，以及坐标点对应的数值就好了。
偏红色的表示对应的数值较高，为等高线的山峰；偏蓝色的区域表示对应的数值较低，为等高线的山谷
- 【MZKit】简单自动化组织分区 - 2023年11月5日
- 【MZKit教程】质谱成像原始数据文件查看 - 2023年6月29日
- 生物序列图嵌入算法 - 2023年6月29日
One response
[…] 实际上，我们通过Marching Squares算法进行绘制等高线图也可以实现轮廓扫描获取的功能 […]