Vaa3D笔记（part1）

安装过程

预备步骤：安装Qt6

先下载官网的online安装程序（Qt6只能online安装）
安装时注意要安装Qt 5 Compatibility Module
如果后续需要添加插件，进入Qt目录后，运行 MaintenanceTool.exe 即可。
如果下载速度过慢，可以使用镜像源安装，在cmd中，

cd online安装程序的路径
.\qt-unified-windows-x64-4.6.0-online.exe --mirror https://mirrors.ustc.edu.cn/qtproject

mirror后面加常用的国内镜像源都可以。

清华：https://pypi.tuna.tsinghua.edu.cn/simple

阿里云：http://mirrors.aliyun.com/pypi/simple/

中国科技大学 https://pypi.mirrors.ustc.edu.cn/simple/

华中理工大学：http://pypi.hustunique.com/

山东理工大学：http://pypi.sdutlinux.org/ 

豆瓣：http://pypi.douban.com/simple/

参考：
qt 安装太慢的问题如何解决？
pip install 使用国内镜像

方法一：直接安装发行版

以Vaa3D v1.1.2为例，进入以下网址，
https://github.com/Vaa3D/release/releases/tag/v1.1.2
下载对应操作系统的安装包解压即可。

方法二：下载源码后自行编译构建

具体步骤可完全参考，
https://github.com/Vaa3D/Vaa3D_Wiki/wiki/Build-Vaa3D-Qt6-version-on-Windows

注意Qt的shadow build会在源代码目录 .\v3d_external\v3d_main\v3d 的上级目录中创建 build-v3d_qt6-Desktop_Qt_6_5_1_MinGW_64_bit-Release 文件夹，将可执行程序exe和动态链接库dll拖入.\v3d_external\bin中即可。
注意最后构建出的exe文件叫Vaa3D-X.exe，不是老版的Vaa3D.exe了。
后续要导入库的话，需要将github上库的源码下载下来，编译后放入plugins文件夹下。

Vaa3D 插件

Vaa3D 插件有两种使用方式，一是在GUI菜单栏中的Plug-In下直接对当前打开图像使用，然后Ctrl+S保存或另存为。二是使用命令行。

常用指令	功能
Vaa3D-x.exe /h	查看插件列表
Vaa3D-x.exe /h /x ada_threshold	查看某一插件信息（可加.dll或目录）
Vaa3D-x.exe /x ada_threshold	调用插件（后加输入参数）

当然，也可以用其他编程语言（如python）批量执行命令行，但注意linux shell和windows cmd的语法是不同的。

通用概念

通过学习源码，先在此总结一些通用的类，以方便查阅。
一些预先定义的类或数据结构：

类	含义
V3DLONG	可理解为整数
V3DPluginCallback2	可理解为插件整体的一个回调类指针，其下定义了很多方法、属性
V3DPluginArgList	参数列表的指针，输入输出参数都采用这个类定义的对象
Image4DSimple	图像指针
V3D_UINT8	8位无符号整型
V3D_UINT16	16位无符号整型
V3D_FLOAT32	32位浮点数
v3dhandle	图像窗口指针
LocationSimple	可理解为一个坐标对象，包含x、y、z三个属性
LandmarkList	可理解为定义了Landmark类的列表，其中Landmark类和LocationSimple类类似。

一些预先定义的函数：

函数	功能
v3d_msg	Vaa3D窗口的弹出消息框接口，输入文本参数即可

此外还有一些常用于继承的类，这些类都定义在 v3d_basicdatatype.h 或 v3d_interface.h 中，但因为本人还未看这两个源码，所以很多类的介绍只能模糊说明。

因为V3D处理的图片中数据类型有三种，算法必须要用模板 template 去定义，可以利用这个方法快速找到算法核心代码。

ada_threshold

.\plugins\image_thresholding\Simple_Adaptive_Thresholding\ada_threshold.dll
官方介绍如下：

Simple adaptive thresholding: for each voxel, compute a threshold which is the average intensity of its neighboring voxels and then subtract the threshold from the current voxel’s intensity value. If the result is < 0, then set it as 0. The neighborhood is defined along 6 axial directions in 3D, with N samples of each direction (N – the ‘number of sampling points’ in the parameter setting dialog), and M voxels between every nearest pair of samples (M – the ‘sampling iterval’ in the parameter setting dialog).

翻译如下：

简单自适应阈值化：对于每个体素，计算一个阈值，该阈值是其相邻体素的平均强度，然后从当前体素的强度值中减去该阈值。如果结果小于0，则将其设置为0。邻域沿着3D空间的6个轴向定义，在参数设置对话框中的每个方向上有N个采样点（N是“采样点数量”），并且在每对最近的采样点之间有M个体素（M是“采样间隔”）。
voxel（volume pixel）：体积元素、立体像素、体素

下面通过插件的源码进行分析学习。

Image4DSimple *subject = callback.loadImage(inimg_file);

这一行将图片导入，命名为subject。

Image4DSimple outimg;
outimg.setData((unsigned char *)pData, sz0, sz1, sz2, sz3, subject->getDatatype());

这两行将图片重新保存到outimg中输出。

V3DLONG sz0 = subject->getXDim();
V3DLONG sz1 = subject->getYDim();
V3DLONG sz2 = subject->getZDim();
V3DLONG sz3 = subject->getCDim();

V3DLONG sz_data[4]; sz_data[0]=sz0; sz_data[1]=sz1; sz_data[2]=sz2; sz_data[3]=1;

这一部分可以看出，sz0、sz1、sz2、sz3分别表示图像的三个维度和色彩维度的大小。

了解了这些，看算法函数BinaryProcess。

template <class T>
void BinaryProcess(T *apsInput, T * aspOutput, V3DLONG iImageWidth, V3DLONG iImageHeight, V3DLONG iImageLayer, V3DLONG h, V3DLONG d)
{
	V3DLONG i, j,k,n,count;
	double t, temp;

	V3DLONG mCount = iImageHeight * iImageWidth;
	for (i=0; i<iImageLayer; i++)
	{
		for (j=0; j<iImageHeight; j++)
		{
			for (k=0; k<iImageWidth; k++)
			{
				V3DLONG curpos = i * mCount + j*iImageWidth + k;
				V3DLONG curpos1 = i* mCount + j*iImageWidth;
				V3DLONG curpos2 = j* iImageWidth + k;
				temp = 0;
				count = 0;
				for(n =1 ; n <= d  ;n++)
				{
					if (k>h*n) {temp += apsInput[curpos1 + k-(h*n)]; count++;}
					if (k+(h*n)< iImageWidth) { temp += apsInput[curpos1 + k+(h*n)]; count++;}
                    if (j>h*n) {temp += apsInput[i* mCount + (j-(h*n))*iImageWidth + k]; count++;}//
					if (j+(h*n)<iImageHeight) {temp += apsInput[i* mCount + (j+(h*n))*iImageWidth + k]; count++;}//
					if (i>(h*n)) {temp += apsInput[(i-(h*n))* mCount + curpos2]; count++;}//
					if (i+(h*n)< iImageLayer) {temp += apsInput[(i+(h*n))* mCount + j* iImageWidth + k ]; count++;}
				}
				t =  apsInput[curpos]-temp/(count);
				aspOutput[curpos]= (t > 0)? t : 0;
			}
		}
	}
}

这里的h和d在另一个函数中定义：

V3DLONG h;
V3DLONG d;
if (method_code == 1)
{
	h = 5;
	d = 3;
}
else
{
	if( method_code == 2)
	{
		AdaTDialog dialog(callback, parent);
		if (dialog.exec()!=QDialog::Accepted)
		return;
		else
		{
			h = dialog.Ddistance->text().toLong()-1;
			d = dialog.Dnumber->text().toLong()-1;
			printf("d% h,d% d \n ",h,d);
		}
	}
}

综合两段代码可以看出，h代表步长，d代表采样点个数或者说采样点范围，简单自适应阈值化的算法原理是，在每个点，沿着6个方向以h为步长，各取最多d个点求这些点的平均值作为当前点的背景值，用当前点的原始值减去背景值作为阈值化处理后的结果。method1为参数取默认值（h=5,d=3）的方法，method2为输入指定参数的方法。
在Vaa3D中实验后，发现此算法的method1因为步长、采样点个数较少会使得中间的细胞体消失（因为步长和背景点个数较少），不过也可以以这个为思路设计定位细胞体。

APP1（from neuron_tracing）

论文相关

Hanchuan P,Fuhui L,Gene M. Automatic 3D neuron tracing using all-path pruning.[J]. Bioinformatics (Oxford, England),2011,27(13).
如果仅需要了解算法，论文重点在下载后的P3-P7中。
论文中常用缩写含义如下：

缩写	含义
APP	all-path pruning，全路径修剪
SC	structural component，神经元重建对象中的所有结构
ICR	initial over-complete reconstruction，初始完全重建图，是
MCMR subgraph algorithm	maximal-covering minimal redundant subgraph algorithm，最大覆盖下最小冗余子图算法，是用于修剪重建图的一种策略

APP1 算法属于神经元的重建算法，重建算法的目标是最终生成准确的神经元结构数字化的模型。以APP1算法为例，在论文《Automatic 3D neuron tracing using all-path pruning - PMC》的2.1节中提到：

Typically, a neuron reconstruction is described as a tree graph, which has a root/seed reconstruction node, many leaf nodes, branching nodes and other inter-nodes.

可知APP1算法，将神经元结构表示为一个重建根节点（种子节点），若干个叶子节点、分支节点和中间节点的组合。
论文中将图像的平均强度值$t_a$作为全局阈值来定义图像的前景（大于这个阈值），假设所有值大于$t_a$的体素都是神经元的一部分，否则则为图像的背景。这是一个十分保守的设定以确保不会丢失信息，在定义完全局阈值后，定义了333的中值滤波器来消除噪声。
之后，在前景上定义无向图$G(V,E)$，将顶点定义为前景的所有像素点并记录其位置，将边设定为只存在于两个相邻点之间（这里相邻的定义为，三个坐标差值最多为1且不全为0），且边的权重定义为如下公式：
$$e(v_0,v_1)=||v_0-v_1||(\frac{g_I(v_0)+g_I(v_1)}{2})$$$$g_I(p)=exp({\lambda }_I(1-\frac{I(p)}{I_{max}}{)}^2)$$
其中，第一个公式中的范数采用欧几里得距离，$g_I(\cdot )$定义为和强度成反比的函数（点越亮，权重越小），${\lambda }_I$一般取10。如果要求种子节点到其他节点的最短路径则采用 Dijkstra algorithm，考虑边都是只存于邻近点之间，十分稀疏，该算法的复杂度为$O(N\log N)$（$N$为节点数）。

In the resultant shortest path map, we search vertices that have no child, and call such a set as the leaf set. Obviously, we can back trace a path from every leaf vertex to the seed location, which is also called the root vertex. All these paths share many common pieces. We organize the entire solution as a tree graph. Since this step detects all ordered paths from the root to every image foreground voxels, and thus include no false negative, we call this tree graph the all-path reconstruction, which is an ICR.

这一段文字需要单独拉出来理解，首先the resultant shortest path map指的是由根节点出发的Dijkstra algorithm生成的图，由这个算法的原理我们知道，得到的这个图是一棵树，接下来说这棵树上没有子节点的点就是叶子节点，注意前面边的权重被定义为点之间的欧几里得距离乘以一个和两个点亮度成反比的函数，也就是说神经元边上那些更暗的点之间的距离会更远，因此最短路径走的时候就不会沿着神经元图像的边界去走，而是会走在神经元中间，再延申到边界，这就保证了我们这棵树的子节点能够实际上等价于神经元的边界。
由以上步骤得到的树图，被称为初始完全重建图（ICR），至此，算法第一步完成。

算法的第二步是对已有的 ICR 进行修剪，采用的是叫 MCMR (maximal-covering minimal redundant) subgraph search 的一种策略，论文中提供的具体的算法有Dark-leaf pruning (DLP)、Covered-leaf pruning (CLP)、 Inter-node pruning (INP) and refinement这三种。

DLP的原理就是把亮度低于某个阈值的叶子节点剪掉，重复此操作直到所有叶子节点的亮度都大于这个阈值。

CLP的原理核心思想是一种覆盖法，对于每个节点，以该节点为中心定义一个可调节半径的球，将这个球的半径从0逐渐增大直到这个球中有至少0.1%的体素的值低于$t_a$（即，使这个球在前景范围内尽可能大），将所有节点及其最大覆盖半径（covering radius）作为一个SC保存下，这个SC也可以理解为是每个节点的最大球的覆盖范围。
根据覆盖半径的定义，定义$\Omega (\cdot )$为覆盖体积或覆盖质量，代表对应结点最大覆盖球的体积（包含体素的数量），或者质量（包含体素的值之和），文章中也提到用质量的定义会更加合理。之后定义节点 $a$ 被节点 $b$ 显著覆盖的条件为：
$$\frac{\Omega (a\cap b)}{\Omega (a)}\ge 0.9$$
我们可以认为，那些能被其他点覆盖的点，是可以移除的点，但我们仍需要定义一个移除点的顺序。注意到，当一个叶子节点不被任何其他节点覆盖的时候，这个节点能一直留存下来，否则这个叶子节点可以安全地被剔除，因此可以像DLP一样，每次检查所有叶子节点，将可以剪掉的叶子节点都剪掉，重复此操作直到不再有叶子节点被剪掉。

INP认为，虽然CLP已经把ICR剪得很简洁了，但还可以再减少重建图的复杂度，其核心思想是，在不损害图的连通性和完整性的基础上，移除叶子节点到分支节点到根节点之间冗余的中间节点。其算法为，先用一个队列 L 存储所有子节点，从一个子节点 $a$ 出发，如果其父节点 $b$ 能被 $a$ 显著覆盖，则删除节点 $b$ 并用其父节点代替 $a$ 的父节点。如果 $b$ 未被删除，则从 $b$ 出发重复此操作，直到到达分支节点或者父节点，如果到达的是分支节点，把该分支节点先压入队列中，对队列中的所有节点重复以上步骤，直到队列为空。最终，我们成功将ICR修剪为最精简的树。

代码概览

这里参考 app1_connector.cpp 的代码，该文件导入了以下头文件。

#include "vn_app1.h"
#include "vn_imgpreprocess.h"
#include "basic_surf_objs.h"
#include "marker_radius.h"
#include "v3dneuron_gd_tracing.h"
#include "fastmarching_dt.h"
#include "volimg_proc.h"

其中 vn_app1.h 导入了 vn.h 导入了 vn_imgpreprocess.h。

首先 vn_imgpreprocess.h 和 vn_imgpreprocess.cpp 中定义了多种采样、值映射函数，具体在使用到之后再说。

basic_surf_objs.h 这个头文件，其实我在文件夹内只找到 my_surf_objs.h，里面定义了marker文件、swc文件的读取函数。

vn类

vn.h 中定义了 PARA_VN 类，VN是V3DNeuron的缩写。

struct PARA_VN //VN - V3DNeuron
{
    Image4DSimple * p4dImage;
    int xc0, xc1, yc0, yc1, zc0, zc1; //the six bounding box boundaries
    LandmarkList landmarks;
    PARA_VN() 
    {
        p4dImage = NULL;
        xc0 = xc1 = yc0 = yc1 = zc0 = zc1 = 0;
        landmarks.clear();
    }
    bool initialize(V3DPluginCallback2 &callback);
}

上述代码省略了初始化函数中的内容，可以看出PARA_VN中主要包含了一个存储图像的对象，六个代表图像边界的数值，以及一个列表对象 landmarks，初始化函数中有用到对象的size方法获取列表长度，removeAt方法删除指定位置元素，该列表中每个元素有三个属性x、y、z，故猜测该列表用于存储标记点坐标。
初始化函数主要功能是根据当前窗口图像初始化该类的所有属性（包括landmarks），并删除超出边界的landmarks中的元素。
目前为止可以看出 vn.h 头文件创建了神经元类。

vn_app1.h 中定义了 PARA_APP1 继承自 PARA_VN。

int  bkg_thresh; //for initial reconstruction generation
int  visible_thresh; //for dark pruning
int  channel;
int  downsample_factor; //when set to be 0, then set b_256cube to be 1 which means downsample based on fit to a 256px cube automatically
int  b_256cube;

bool b_menu;

QString inimg_file, inmarker_file, outswc_file;

PARA_APP1()
{
    bkg_thresh = -1; //-1 for auto-thresholding, i.e. using the ave of image channel; //10; //30; change to 10 by PHC 2012-10-11
    visible_thresh = 30;
    channel = 0;
    downsample_factor = 2;
    b_256cube = 1; //whether or not preprocessing to downsample to a 256xYxZ cube UINT8 for tracing

    b_menu = true;

    inimg_file = "";
    inmarker_file = "";
    outswc_file = "";
}

bool app1_dialog();

bool proc_app1(V3DPluginCallback2 &callback, PARA_APP1 &p, const QString & versionStr);

其中 app1_dialog 函数在此处已完整定义，proc_app1 函数则只是声明，完整定义在 app1_connector.cpp 中。
该类的初始化属性中，包括两个阈值参数、下采样因子、输入的图片文件名、标记文件名和输出的swc文件（用于显示重建神经元的图像）。app1_dialog 函数用于生成算法的GUI对话框。

markerRadius 函数

前面论文中提到的CLP法，需要计算每个节点的最大覆盖半径，该方法的实现在 marker_radius.h 中。
marker_radius.h 中定义了 markerRadius 函数（分为四种方法），该函数输入图像、图像大小、标记点、背景阈值四个参数，输出标记点的最大覆盖半径。

template<class T> double markerRadius_accurate(T* &inimg1d, V3DLONG * sz, MyMarker & marker, double thresh);

函数内，令 r 从1 开始，每次+1，最大到max_r（图像长宽高中最小值的一半）遍历，每次遍历统计覆盖球中总体素个数和背景体素个数，当背景体素个数除以总体素个数大于0.0001时，或者覆盖球超出图像范围时结束循环并输出最大半径。

v3dneuron_gd_tracing 函数

v3dneuron_gd_tracing.h 以及 v3dneuron_gd_tracing.cpp 中定义了一个主要接口函数，

NeuronTree v3dneuron_GD_tracing(unsigned char ****p4d, V3DLONG sz[4], 
								LocationSimple & p0, vector<LocationSimple> & pp, 
								CurveTracePara & trace_para, double trace_z_thickness);

和非常多的enabling function。

附录

Vaa3D常用文件格式

格式	介绍
V3DRAW	Vaa3D源生图像数据文件
V3DPBD	Vaa3D源生图像数据文件
TIF	医学影像常用图像文件，可存储2维或3维图像
SWC	是一种打包文件，常用于存储重建后的神经元
MARKER	图像处理时生成的临时文件，记录了标记点，不能单独打开

二维tif图转为三维tif图

这里使用python3.7实现。
需要安装 pylibtiff 库和 tifffile 库。
注意 pylibtiff 不要直接pip安装，可以在https://www.lfd.uci.edu/~gohlke/pythonlibs/​中找到与python版本对应的 whl 安装包，下载后，再在pip安装 wheel 库后，再pip安装下载好的安装包。
另外 pylibtiff 库在处理三维数据时，可能存在问题，保存三维 tiff 数据可以用 tifffile 库（注意3个f）解决。

# tif2to3.py
import os
import numpy as np
import glob
from libtiff import TIFF
import tifffile

# 读取所有要合并的图片的文件名
tiflist = []
filepath = r'E:\assignment\raw_data\20230704\ZSeries-07042023-0844-10867'
for file in glob.glob(os.path.join(filepath,'*.tif')):
    tiflist.append(file)
tiflist.sort()

# 用numpy库连接生成三维数组
img3d = []
for i in range(len(tiflist)):
    tif = TIFF.open(tiflist[i],mode='r')
    image = tif.read_image()
    img3d.append(image.copy())

# 用tifffile生成3Dtif图片
img3d = np.array(img3d)
tifffile.imsave(filepath+'.tif',img3d)

参考：
OpenCV—Python PyLibTiff_psd 图像基本操作以及图像格式转换
Python将多张2D TIFF图片转为一个3D TIFF文件 by 笃谷