文章阅读目录大纲
https://github.com/rsharp-lang/NRRD
NRRD(Nearly Raw Raster Data)是一种用于存储类似于热图成像数据的文件格式。其实我们可以将NRRD看作为类似于bitmap之类的未压缩的原始光栅图像文件。只要我们有对应的解码方式，我们就可以像查看普通图片文件一样查看NRRD文件。
就像普通的图片文件一样，NRRD文件之中一般存储的是一个二维矩阵数据或者三维矩阵数据。但是NRRD文件之中的光栅数据和普通图像文件之中的光栅数据不太一样：
- 普通图像文件之中的光栅数据一般是RGB颜色值，既每一个像素点上的颜色是确定的一个值
- NRRD图像文件之中的光栅数据一般是一个数值，我们需要经过一定的线性变换，将这些光栅数值映射为对应的颜色，才可以将NRRD中存储的图像显示出来
在上面所提到的线性变化转换过程，其实就是一个热图绘制的过程。我们一般按照不同的颜色谱做线性变换映射，就可以得到对应的不同颜色系列下的NRRD热图成像渲染结果。对于NRRD图像文件的热图成像渲染原理，其实是和质谱成像的渲染原理一摸一样的（对于质谱成像渲染而言，其主要的原理也就是将对应的扫描点上的目标离子的intensity值取出，构建出一个和NRRD文件中的光栅矩阵数据一摸一样的矩阵数据，基于这个矩阵数据进行线性变换映射到对应的颜色值完成热图成像可视化操作）。
文件格式介绍
和一般的原始数据文件一样，在NRRD文件之中，分为文件头和数据区域两个部分。在NRRD文件之中，文件头是由ASCII编码的文本组成的，我们可以很容易的通过诸如vscode之类的代码编辑器打开NRRD文件，以文本的形式查看文件的文件头信息，例如：
从上面的截图之中，我们可以看得出一些信息，如下所示：
- type: unsigned int 光栅原始数据是以无符号的32位整形数类型进行存储的，既接下来做原始数据矩阵的读取我们需要按照32位整形数的格式读取出来
- dimension: 3 光栅矩阵数据是一个三维的点云扫描模型
- sizes: 528 320 456 矩阵的三个维度大小信息
- endian: little 光栅原始数据是以小端编码的32位无符号整形数存储的
- encoding: gzip 我们需要将原始数据进行gzip流解压缩之后在进行矩阵的读取操作！
解析上面截图所示的文件头信息非常简单，我们只需要使用.NET系统框架所提供的StreamReader进行字符串读取操作即可。
Private Function loadNrrdHeader() As Header
Dim read As New StreamReader(file, System.Text.Encoding.ASCII)
Dim magic As String = read.ReadLine
Dim line As Value(Of String) = ""
Dim metadata As NamedValue(Of String)
Dim header As New Header With {
.magicNumber = [Enum].Parse(GetType(MagicNumber), magic),
.metadata = New Dictionary(Of String, String)
}
Dim readSize As New List(Of String)
Do While Not (line = read.ReadLine).StringEmpty
Call readSize.Add(line)
If line.First = "#"c Then
comments.Add(line)
Else
metadata = line.GetTagValue(":", trim:=True)
header.add(metadata.Name, metadata.Value)
End If
Loop
' 20230216
' evaluate the actual text offset
' the streamReader has a bug about buffer size
' it could cause the incorrect stream position
scan0 = Aggregate str As String
In readSize
Let chars = str.Length + 2
Into Sum(chars)
Return header
End Function假设你使用的是微软的.NET环境中的StreamReader进行NRRD文件数据读取操作，在这里面有一个坑需要注意：因为StreamReader对象是按照一定的Buffer大小进行流数据的读取操作的，所以在读完NRRD的文件头信息之后，因为Buffer读取的原因，Stream流的位置很有可能已经偏移到光栅矩阵数据区内部了。所以我们还需要统计一下我们所读取的文件头的字符数量，来判断出真正的光栅原始数据在流之中的位置scan0，再进行后续的光栅矩阵原始数据读取操作。
读取光栅矩阵原始数据
在解析完NRRD的文件头数据之后，我们一般会需要有四个比较重要的信息用于后续的光栅矩阵原始数据的读取操作，从NRRD的文件头元数据之中我们必须要了解到如下参数信息：
- type: 我们使用何种数据类型的读取方式进行原始的Byte转换为对应的光栅矩阵数值。在这里的数据类型是和c语言之中的数据类型保持一致的。对于type字段而言，一般可能会出现的数据类型有：int8, uint8, int16, uint16, int32, uint32, int64, uint64, int, float, double, uchar, block
- encoding: 光栅原始数据是以什么压缩格式存储的？在NRRD文件之中，支持的数据压缩格式有：raw, text, ascii, hex这几种类型其实都表示一个意思：无压缩，此时我们可以直接进行原始数据的读取操作；gzip, bzip2分别表示光栅矩阵原始数据数据流的压缩格式分别位gzip流压缩以及bzip2流压缩，此时我们需要将从文件中读取出来的原始数据进行一次对应的解压缩操作，基于解压缩出来的流程据再进行光栅矩阵原始数据的读取操作
- endian：数字的大小端对齐方式，如果在进行数字矩阵读取操作的时候，没有设置正确的大小端对齐方式，那么我们解码出来的数字矩阵一般会是无序的乱码
- size: 从size属性里面我们一般可以了解到我们的矩阵是如何从原始数据中读取出来的一维数据上构建出来的
因为进行流数据的加压缩操作比较简单，假设我们已经可以正确的将光栅数据给解压缩出来了，我们在这里就主要重点介绍一下对原始数据的读取操作过程。对于NRRD的光栅矩阵数据而言，不管这个矩阵数据是多少维度的，其最终在原始数据文件之中的存储形式依然是必须要以一维线性向量的格式进行存储，所以基于NRRD的文件头信息，我们就可以很轻松的计算出我们所需要读取的一维向量的长度：
Dim sizes As Integer() = options.sizes
Dim totalLen As Integer = 1
For Each int As Integer In sizes
totalLen *= int
Next只需要将所有的维度值全部相乘起来即可知道矩阵展开后的一维向量的长度有多长了，在上面我们计算得到了最终的一维向量的总长度值totalLen。接下来我们只需要按照设定的endian大小端对齐方式，读取出所计算出来的长度的一维向量即可：
If type = Types.block Then
' Don't do anything special, just return the slice containing all blocks.
Return buffer.ReadBytes(totalLen * options.blockSize)
Else
buffer.ByteOrder = options.endian
End If
Select Case type
Case Types.int8 : Return buffer.ReadSBytes(totalLen)
Case Types.uint8 : Return buffer.ReadBytes(totalLen)
Case Types.int16 : Return buffer.ReadInt16s(totalLen)
Case Types.uint16 : Return buffer.ReadUInt16s(totalLen)
Case Types.int32 : Return buffer.ReadInt32s(totalLen)
Case Types.uint32 : Return buffer.ReadUInt32s(totalLen)
Case Types.float : Return buffer.ReadSingles(totalLen)
Case Types.double : Return buffer.ReadDoubles(totalLen)
Case Else
Throw New NotImplementedException(type.Description)
End Select接着我们就可以在这个所读取的一维向量的结果基础上，解析为二维矩阵或者三维矩阵即可。针对于二维矩阵，我们只需要按照第一个维度长度进行一维向量的切割即可，例如：
Dim dimensionSize As Integer() = metadata.sizes
Dim rawDbls As Double() = CLRVector.asNumeric(rawdata)
Dim plant2D As Double()() = rawDbls.Split(dimensionSize(0))
If plant2D.Length <> dimensionSize(1) Then
Throw New InvalidProgramException
End If对于三维的矩阵数据，我们要做的无非就是在二维折叠的结果基础上，再做一次第二个维度的折叠即可，例如：
Dim dimensionSize As Integer() = metadata.sizes
Dim rawDbls As Double() = CLRVector.asNumeric(rawdata)
Dim fold1 = rawDbls.Split(dimensionSize(0))
Dim fold2 = fold1.Split(dimensionSize(1))
If fold2.Length <> dimensionSize(2) Then
Throw New InvalidOperationException
End If是不是非常的简单？!!!
使用NRRD程序包进行医学影像数据读取
基于上面的所描述的原始数据文件的读取操作过程，我将代码打包成了一个R#程序包，开源在github上，叫做NRRD。基于这个程序包，我们可以非常简单的进行NRRD光栅矩阵对象的读取操作，例如只需要下面的几行代码：
require(NRRD);
# open raw NRRD data file
nrrd = NRRD::nrrdRead("..\\data\\stent.nrrd");
# view of the file header data
header = as.list(NRRD::metadata(nrrd));
# and get raster data for image rendering
raster = NRRD::getRaster(nrrd);那现在我们很轻松的基于R#脚本就拥有了从NRRD文件之中所读取出来的光栅矩阵数据之后，如果我们需要将得到光栅矩阵数据进行可视化，该怎样做呢？其实，如果我们了解过热图成像或者质谱成像的原理的话，实际上对于这个光栅矩阵的原始数据进行成像的原理应该就会很清楚了。在我们拿到这个矩阵之后，可以将矩阵的行和列看作为二维图像空间之中的x和y坐标信息，然后对应的矩阵中的单元格值可以映射为一个对应的颜色，即可将从NRRD文件之中拿到的光栅矩阵数据给可视化出来。将光栅矩阵中的数值映射为对应的颜色值的方法原理，大家可以参考一下《【热图数据可视化】颜色插值计算原理》的内容介绍，一摸一样。
实现上面的热图成像的代码，大家可以参考rasterHeatmap
require(graphics2D);
# for a 2d image
bitmap(file = "./raster_image_heatmap.png", size = [1024,1024], fill ="black");
# draw heatmap
graphics2D::rasterHeatmap(raster, colorName = "viridis");
dev.off();对于一个三维扫描的数据，其最本质上仍然是基于二维的矩阵堆叠出来的。如果我们需要将其某一断层显示出来，那么我们可以通过getRasterLayer
# for a specific layer from the 3d scan data
bitmap(file = `./stent/raster__${i}.png`, size = [256,256], fill ="black");
# draw heatmap
graphics2D::rasterHeatmap(
x = raster → NRRD::getRasterLayer(i),
colorName = "viridis:turbo"
);
dev.off();
- 【MZKit】简单自动化组织分区 - 2023年11月5日
- 【MZKit教程】质谱成像原始数据文件查看 - 2023年6月29日
- 生物序列图嵌入算法 - 2023年6月29日
No responses yet