（一）由数组到图像：统一matlab、avizo、quant3d、homogenization坐标系
- 1 二维图像与矩阵
- 1.1 matlab中构造图像
- 1.2 avizo与quant3d中显示二维图像
- 2 三维图像与三维数组
- 2.1 matlab中显示三维图像
- 2.2 avizo与quant3d中显示三维图像
- 2.3 总结
- 3 homogenization源码
- 4 以数组为基准得到一致的可视化结果
对于三维数组，
以第一维度索引增加的方向代表三维图像的x轴正向，第二维度索引增加为y轴正向，第三维度索引增加为z轴正向。
在探索三维图像与三维数组的关系之前，要先了解二维图像与矩阵的关系。
1 二维图像与矩阵
1.1 matlab中构造图像
构造一个矩阵，行索引增大为x轴正向，列索引增大为y轴正向（以第一维度索引增加的方向代表三维图像的x轴正向，第二维度索引增加为y轴正向）
%边长为100的正方形图片
a=zeros(100,100);
%从原点出发，沿着x轴方向的矩形，长度为100，宽度为10
a(:,1:10)=1;
%从原点出发，沿着y轴正向的矩形，长度为50，宽度为10
a(1:10,1:50)=1;
%转为8-bit格式，将1变成255，显示为白色，0为背景黑色
a=uint8(255 * a);
%matlab自带的显示图片的函数
imshow(a,[])
%写出tiff
imwrite(a, 'aa.tiff', 'tiff', 'Compression', 'none');结果如下：
由此可见，imshow和imwrite写出的图像是一致的，并且图像和矩阵的坐标系是一致的，即行索引增大为x轴正向，列索引增大为y轴正向，图像的原点位于左上角。
1.2 avizo与quant3d中显示二维图像
把之前用matlab生成的图片aa.tiff分别导入avizo中和quant3d中得到的结果如下：
由此可见，avizo和matlab中原点相同坐标轴不同，quant3d中原点位于图像的左下角。要想在avizo中得到和matlab一致的坐标系，需要将matlab矩阵进行转置；要想在quant3d中得到和matlab一致的坐标系，需要将matlab矩阵进行转置逆时针旋转90度。可如下写出代码：
Aavizo=a';
imwrite(Aavizo, 'Aavizo.tiff', 'tiff', 'Compression', 'none');
Aquant3d=rot90(a);
imwrite(Aquant3d, 'Aquant3d.tiff', 'tiff', 'Compression', 'none');转化后的结果重新导入avizo和quant3d，结果将完全一致
2 三维图像与三维数组
首先构建三维图像，构建a和b两种类型的切片，然后沿着z轴方向堆叠起来，a在z较小的位置，b在z较大的位置。即Za<Zb
%a与上文中图片一致
a=zeros(100,100);
a(:,1:10)=1;
a(1:10,1:50)=1;
%b为图片的中心处有一个正方形
b=zeros(100,100);
b(40:60,40:60)=1;
for i=1:50
A3d(:,:,i)=a;
end
for i=51:100
A3d(:,:,i)=b;
end2.1 matlab中显示三维图像
利用如下代码进行三维可视化
figure();
col=[.7 .7 .8];
hiso = patch(isosurface(A3d,0),'FaceColor',col,'EdgeColor','none');
hiso2 = patch(isocaps(A3d,0),'FaceColor',col,'EdgeColor','none');
axis equal;%axis off;
lighting phong;
isonormals(A3d,hiso);
alpha(0.5);
set(gca,'DataAspectRatio',[1 1 1])
camlight;
hold on;
axis off
%x轴红色
quiver3(0,0,0,110,0,0,'r','LineWidth',2,'MaxHeadSize',0.5);
text(110,0,0,'x','FontSize',12,'Color','r');
%y轴蓝色
quiver3(0,0,0,0,110,0,'b','LineWidth',2,'MaxHeadSize',0.5);
text(0,110,0,'y','FontSize',12,'Color','b');
%z轴绿色
quiver3(0,0,0,0,0,110,'g','LineWidth',2,'MaxHeadSize',0.5);
text(0,0,110,'z','FontSize',12,'Color','g');
view(45,45)结果如下：
由此可见，在matlab的三维展示方法中，会像avizo那样自动的把x轴和y轴进行交换。但是对于z轴而言，就是第三维度索引增加为z轴正向。
2.2 avizo与quant3d中显示三维图像
首先利用matlab导出序列的二维图像，按照顺序，z刻度越小的编号越小
% 获取二值数组的尺寸
[m, n, p] = size(A3d);
% 循环遍历三维数组，将每个切片保存为8-bit TIFF图像
for k = 1:p
% 获取当前切片
currentSlice = A3d(:, :, k);
% 将二值图像转换为8-bit灰度图像（0表示黑色，255表示白色）
grayscaleImage = uint8(255 * currentSlice);
% 生成文件名（可以根据需要定义文件名的格式）
outputFileName = sprintf('output_slice_%03d.tiff', k);
% 保存当前切片为TIFF图像
imwrite(grayscaleImage, ['.\test0929\',outputFileName], 'tiff', 'Compression', 'none');
end结果如下：
将这组序列图片分别导入avizo和quant3d中，结果如下：
2.3 总结
- matlab、avizo和quant3d三个软件中均按照第三维度索引增加的方向为x轴正向，即在输出图像的时候，编号较小的要对应小的z轴刻度值。
- matlab在本文中所用的三维体积成像的方式，会导致和avizo中相似，因而在展示前要将切片进行转置，从而达到交换x轴和y轴的目的。
3 homogenization源码
该开源代码包含了由拓扑结构生成体素，并计算结构的弹性张量，其主要功能如下
为了验证其内置定义的坐标系，利用其本身的GenerateVoxel.m函数即可。
首先在topoly文件夹下新建一个hjp.txt的文件。x、y、z是经过标准化，即最大坐标为1。GRID表示坐标点，例如ID1的坐标为（0，0，0），ID2的坐标为（0.5，0，0）。STRUT是连接GRID的线，例如STRUT1为连接GRID1和GRID2的线。
//Grid ID x y z
GRID 1 0 0 0
GRID 2 0.5 0 0
GRID 3 0.5 0 0
GRID 4 0.5 1 0
GRID 5 0 0 0.5
GRID 6 0 0 1
//Strut ID Start End
STRUT 1 1 2
STRUT 2 3 4
STRUT 3 5 6然后运行命令
[voxel,Density] = GenerateVoxel(10,'topology/hjp.txt',0.1);
% the model has 10 voxels along each axis.
% 0.1 density is the relative density, namely volume fraction观察voxel数组的情况并与预期的拓扑结构做对比
说明开源代码homogenization中，坐标系的方向定义与默认一致，即以第一维度索引增加的方向代表三维图像的x轴正向，第二维度索引增加为y轴正向，第三维度索引增加为z轴正向。
4 以数组为基准得到一致的可视化结果
对于给定数组，以第一维度索引增加的方向代表三维图像的x轴正向，第二维度索引增加为y轴正向，第三维度索引增加为z轴正向。用以下结构作为示例进行完整的展示。
将matlab三维可视化以及将数组输出的图片封装成函数，完整代码如下：
clc
clear
A3d=zeros(100,100,100);
A3d(1:50,1:10,1:10)=1;
A3d(46:55,1:100,1:10)=1;
A3d(1:10,1:10,51:100)=1;
myvolshow(A3d)
MattoAvizo(A3d,'.\testavizo')
MattoQuant3d(A3d,'.\testquant3d')
function myvolshow(A3d)
%以行索引递增为x轴正向
%列索引递增为y正向
%第三维度递增为z轴正向
[m, n, p] = size(A3d);
newA3d=zeros(n,m,p);
for i=1:p
newA3d(:,:,i)=A3d(:,:,i)';
end
figure();
col=[.7 .7 .8];
hiso = patch(isosurface(newA3d,0),'FaceColor',col,'EdgeColor','none');
hiso2 = patch(isocaps(newA3d,0),'FaceColor',col,'EdgeColor','none');
axis equal;
lighting phong;
isonormals(newA3d,hiso);
alpha(0.5);
set(gca,'DataAspectRatio',[1 1 1])
camlight;
hold on;
axis off
%x轴红色
quiver3(0,0,0,1.1*m,0,0,'r','LineWidth',5,'MaxHeadSize',0.5);
text(1.1*m,0,0,'x','FontSize',30,'Color','r');
%y轴蓝色
quiver3(0,0,0,0,1.1*n,0,'b','LineWidth',5,'MaxHeadSize',0.5);
text(0,1.1*n,0,'y','FontSize',30,'Color','b');
%z轴绿色
quiver3(0,0,0,0,0,1.1*p,'g','LineWidth',5,'MaxHeadSize',0.5);
text(0,0,1.1*p,'z','FontSize',30,'Color','g');
view(135,45)
end
function MattoAvizo(A3d,path)
%以数组行索引递增为avizo中的x轴正向
%数组列索引递增为avizo中的y正向
%数组第三维度递增为avizo中的z轴正向
% 获取二值数组的尺寸
[m, n, p] = size(A3d);
% 循环遍历三维数组，将每个切片保存为8-bit TIFF图像
for k = 1:p
% 获取当前切片
currentSlice = A3d(:, :, k)';%%%%%%%转置
% 将二值图像转换为8-bit灰度图像（0表示黑色，255表示白色）
grayscaleImage = uint8(255 * currentSlice);
% 生成文件名（可以根据需要定义文件名的格式）
outputFileName = sprintf('slice_%03d.tiff', k);
% 保存当前切片为TIFF图像
imwrite(grayscaleImage, fullfile(path,outputFileName), 'tiff', 'Compression', 'none');
end
end
function MattoQuant3d(A3d,path)
%以数组行索引递增为Quant3d中的x轴正向
%数组列索引递增为Quant3d中的y正向
%数组第三维度递增为Quant3d中的z轴正向
% 获取二值数组的尺寸
[m, n, p] = size(A3d);
% 循环遍历三维数组，将每个切片保存为8-bit TIFF图像
for k = 1:p
% 获取当前切片
currentSlice = rot90(A3d(:, :, k));%%%%%%%逆时针旋转90
% 将二值图像转换为8-bit灰度图像（0表示黑色，255表示白色）
grayscaleImage = uint8(255 * currentSlice);
% 生成文件名（可以根据需要定义文件名的格式）
outputFileName = sprintf('slice_%03d.tiff', k);
% 保存当前切片为TIFF图像
imwrite(grayscaleImage, fullfile(path,outputFileName), 'tiff', 'Compression', 'none');
end
end可以看到，结果均一致，如下图
