图像CT三维整体重建1
图像切片序列，层切片，然后整体重建。视频链接：http://pan.baidu.com/s/1pLNvBtD
录制的视频是算法底层原理讲解，底层代码实现，方便大家真正掌握算法实质，开发出更加出色的算法。录制视频的初衷是：避免读者朋友利用大把时间学习已有常见算法，本系列视频旨在让读者朋友快速深入了解这些常见算法原理，有更多的时间去研究更加高大上算法（价值）。
具体链接在halcom.cn论坛，联系人QQ：3283892722
该论坛是一个学习交流平台，我会逐一的和大家分享学习。
欢迎大家录制视频，并提交给我，我来设置视频，你可在论坛进行打赏分享。
视频专用播放器：http://halcom.cn/forum.php?mod=viewthread&tid=258&extra=page%3D1
运行环境：win7+32bit+matlab2014a
3D重建函数说明：
l.isosurface 等值面函数
调用格式：fv = isosurface(X,Y,Z,V,isovalue)
作用：返回某个等值面（由isovalue指定）的表面（faces）和顶点（vertices）数据，存放在结构体fv中（fv由vertices、faces两个域构成）。如果是画隐函数 v = f(x,y,z) = 0 的三维图形，那么等值面的数值为isovalue = 0。
2.patch函数
调用格式：patch(X,Y,C) 以平面坐标(X, Y)为顶点，构造平面多边形，C是RGB颜色向量
patch(X,Y,Z,C)以空间3-D坐标(X, Y,Z)为顶点，构造空间3D曲面，C是RGB颜色向量
patch(fv) 通过包含vertices、faces两个域的结构体fv来构造3D曲面，fv可以直接由等值面函数isosurface得到
例如：patch(isosurface(X,Y,Z,V,0))
3.isonormals等值面法线函数
调用格式：isonormals(X,Y,Z,V,p)
实现功能：计算等值面V的顶点法线，将patch曲面p的法线设置为计算得到的法线（p是patch返回得到的句柄）。如果不设置法线的话，得到曲面在过渡地带看起来可能不是很光滑
解决办法一：isosurface patch isonormals
实现原理：先定义3元显函数v =f(x, y, z), 则 v = 0 定义的等值面就是z = g(x,y)的3D曲面。利用isosurface函数获取v= 0 的等值面，将得到的等值面直接输入给patch函数，得出patch句柄p，并画出patch曲面的平面视角图形。对p用isonormals函数设置曲面顶点数据的法线，最后设置颜色、亮度、3D视角，得到3D曲面。
代码如下：
f = @(x,y,z) x.*y.*z.*log(1 x.^2 y.^2 z.^2)-10; % 函数表达式
= meshgrid(-10:.2:10,-10:.2:10,-10:.2:10); % 画图范围
v = f(x,y,z);
h = patch(isosurface(x,y,z,v,0));
isonormals(x,y,z,v,h)
set(h,'FaceColor','r','EdgeColor','none');
xlabel('x');ylabel('y');zlabel('z');
alpha(1)
grid on; view(); axis equal; camlight; lighting gouraud
ezimplot3的方法实质是跟方法1是一样的。
补充一个：MC(Marching Cube)方法，这是计算机图像学中很有名的方法。
http://www.mathworks.com/matlabc ... 2506-marching-cubes
具体代码如下：
%% CT切片
clc,clear,close all
file_path ='./试件/';% 图像文件夹路径
img_path_list = dir(strcat(file_path,'*.jpg')); % 获取该文件夹中所有bmp格式的图像
img_num = length(img_path_list);% 获取图像总数量
if img_num > 0 % 有满足条件的图像
k=1;
for j = 1 : 2 : img_num % 逐一读取图像
image_name = img_path_list(j).name; % 图像名
fprintf('正在处理的图像： %d %s\n',j,strcat(file_path,image_name));% 显示正在处理的图像名
img =imread(strcat(file_path,image_name));
if size(img,3)>1
img = rgb2gray(img);
end
img = imresize(img,);% 缩小处理
bw=img>10;
% 连通域检测，计算特征
cc = bwconncomp(bw); % 连通域操作
s= regionprops(cc, {'centroid','area'});% 面积特征
= max(); % 面积最大
bw(labelmatrix(cc)~=id)=0; % 只保留最大二值化块
bw = bwfill(bw,'holes'); % 填充孔洞
img1 = immultiply(img,bw); % 交运算
D(:,:,k) = img1;
k=k+1;
end
end
%% CT重建
% x,y,z为空间坐标，x,y为xy轴坐标
% z为高度：z(:,:,1)=1， z(:,:,2)=2，……，z(:,:,100)=100
= reducevolume(D,);
% D1 = smooth3(D1);% 平滑处理
% 环形曲面
p1 = patch(isosurface(x,y,z,D1, 5,'verbose'),...
'FaceColor','c','EdgeColor','none');
isonormals(x,y,z,D1,p1);
% 上下曲面盖
p2 = patch(isocaps(x,y,z,D1, 5),...
'FaceColor','interp','EdgeColor','none');
view(3);
axis tight;
daspect();
colormap(gray(100));
camlight;
lighting gouraud;
参考：【1】MATLAB绘制3D隐函数曲面的方法总结
页:
[1]