（一）dcm格式医用CT图像处理：avizo与matlab

---

1 对dcm格式三维数据的读取与展示

dcm格式和普通的图片格式对于avizo和matlab的坐标系而言无明显区别，因而可以参考文章1和文章2 。读取方法如下：

clc
clear
dcmpath = '.\xu-invivo-transformed';
[image3d,pixnsize3d]=DCMtoMat(dcmpath);
myvolshow(image3d)

function [image3d,pixnsize3d]=DCMtoMat(dcmpath)
%输入：dcmpath是储存一些列二维dcm图片的文件夹的地址
%输出image3d为matlab坐标系下的三维图像
%以行索引递增为x轴正向
%列索引递增为y正向
%第三维度递增为z轴正向
%输出pixnsize3d每个像素x、y、z三个方向的尺寸，单位为mm
    files_list = dir(fullfile(dcmpath,'*.dcm'));
    bone=[];
    for i=1:length(files_list)
        tmp=dicomread(fullfile(dcmpath, files_list(i).name));
        bone(:,:,i) = tmp';%需要转置从而调整x、y的位置来符合matlab的坐标系
    end
    info=dicominfo(fullfile(dcmpath, files_list(i).name));
    pixnsize3d=zeros(1,3);
    tmp=info.PixelSpacing;
    pixnsize3d(1)=tmp(1)/1e6;
    pixnsize3d(2)=tmp(2)/1e6;
    pixnsize3d(3)=info.SpacingBetweenSlices/1e6;%原单位默认是nm，转为mm
    image3d=bone;
end
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

如下图，可以看到坐标系是一致的。

2 对感兴趣的区域(mask区域)进行分析

有时3d的医用CT图像会搭配有手工标注的raw格式文件作为mask，下面的代码展示了如何对感兴趣的mask区域进行提取和分析：

clc
clear
dcmpath = '.\xu-invivo-transformed';
rawfile='.\xu-invivo mask\xu-invivo.mask.raw';
[image3d,pixnsize3d]=DCMtoMat(dcmpath);
mask=fread(fopen(rawfile));
mask=reshape(mask,size(image3d));
newimage3d=image3d.*mask;
myvolshow(newimage3d)
%提取mask区域内的HU值的分布
Hu=image3d(:);
mask=mask(:);
Hu=Hu(mask>0);
figure
histogram(Hu,'FaceColor', [0.2, 0.4, 0.6], 'FaceAlpha', 0.5)
legend('Data1')
xlabel('HU value')
ylabel('Num of voxels')
set(gca,'fontname','Times New Roman','FontSize',16);

为了进一步完善，新增加了两个函数

1. minRectangularRegion：用于提取三维数组中包含所有非零元素的最小长方体区域
2. myvolshowV2：与myvolshow的区别在于增加了可以显示用红色显示mask位置

clc
clear
dcmpath = '.\xu-invivo-transformed';
rawfile='.\xu-invivo mask\xu-invivo.mask.raw';
[image3d,pixnsize3d]=DCMtoMat(dcmpath);
%读取raw格式的mask文件
mask=fread(fopen(rawfile));
mask=reshape(mask,size(image3d));
myvolshowV2(image3d,mask)%展示mask在原文件中的位置
newimage3d=image3d.*mask;
%为了方便展示，将包含mask的最小长方体提取出来
newimage3d=minRectangularRegion(newimage3d);
myvolshow(newimage3d)%仅展示mask的区域
%提取mask区域内的HU值的分布
Hu1=image3d(:);
mask=mask(:);
Hu1=Hu1(mask>0);
%%

dcmpath = '.\xu-exvivo-original';
rawfile='.\xu-exvivo-mask\xu-exvivo.mask.raw';
[image3d,pixnsize3d]=DCMtoMat(dcmpath);
%读取raw格式的mask文件
mask=fread(fopen(rawfile));
mask=reshape(mask,size(image3d));
myvolshowV2(image3d,mask)%展示mask在原文件中的位置
newimage3d=image3d.*mask;
%为了方便展示，将包含mask的最小长方体提取出来
newimage3d=minRectangularRegion(newimage3d);
myvolshow(newimage3d)%仅展示mask的区域
%提取mask区域内的HU值的分布
Hu2=image3d(:);
mask=mask(:);
Hu2=Hu2(mask>0);


figure
histogram(Hu1,'FaceColor', [0.2, 0.4, 0.6], 'FaceAlpha', 1)
hold on
histogram(Hu2, 'FaceColor', [0.8, 0.1, 0.3], 'FaceAlpha', 0.7);
legend('In-vivo','Ex-vivo')
xlabel('HU value')
ylabel('Num of voxels')
set(gca,'fontname','Times New Roman','FontSize',16);

function [image3d,pixnsize3d]=DCMtoMat(dcmpath)
%读取dcm格式文件
%输入：dcmpath是储存一些列二维dcm图片的文件夹的地址
%输出image3d为matlab坐标系下的三维图像
%以行索引递增为x轴正向
%列索引递增为y正向
%第三维度递增为z轴正向
%输出pixnsize3d每个像素x、y、z三个方向的尺寸，单位为mm
    files_list = dir(fullfile(dcmpath,'*.dcm'));
    bone=[];
    for i=1:length(files_list)
        tmp=dicomread(fullfile(dcmpath, files_list(i).name));
        bone(:,:,i) = tmp';%需要转置从而调整x、y的位置来符合matlab的坐标系
    end
    info=dicominfo(fullfile(dcmpath, files_list(i).name));
    pixnsize3d=zeros(1,3);
    tmp=info.PixelSpacing;
    pixnsize3d(1)=tmp(1)/1e6;
    pixnsize3d(2)=tmp(2)/1e6;
    pixnsize3d(3)=info.SpacingBetweenSlices/1e6;%原单位默认是nm，转为mm
    image3d=bone;
end

function newA3d=minRectangularRegion(A3d)
% 提取三维数组中包含所有非零元素的最小长方体区域
% 查找非零元素的索引
[rows, cols, depths] = ind2sub(size(A3d), find(A3d));
% 计算最小矩形区域的边界
minRow = min(rows);
maxRow = max(rows);
minCol = min(cols);
maxCol = max(cols);
minDepth = min(depths);
maxDepth = max(depths);
% 提取最小矩形区域
newA3d = A3d(minRow:maxRow, minCol:maxCol, minDepth:maxDepth);
end

function myvolshow(A3d)
%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 myvolshowV2(A3d,mask)
%与myvolshow的区别在于增加了可以显示用红色显示mask位置
%A3d和显示的结果和mask均以行索引递增为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)
    hold on;
    [x,y,z]=ind2sub(size(mask),find(mask(:)));
    plot3(x,y,z,'square','Markersize',4,'MarkerFaceColor','r','Color','r');
end

结果如下：