图像CT三维整体重建1

Halcom

图像切片序列，层切片，然后整体重建。
视频链接：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;      % 函数表达式
[x,y,z] = 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([1,1,1]); axis equal; camlight; lighting gouraud

ezimplot3的方法实质是跟方法1是一样的。
补充一个：MC(Marching Cube)方法，这是计算机图像学中很有名的方法。
http://www.mathworks.com/matlabc ... 2506-marching-cubes

具体代码如下：

1. %% CT切片
   
2. clc,clear,close all
   
3. file_path =  './试件/';% 图像文件夹路径
   
4. img_path_list = dir(strcat(file_path,'*.jpg')); % 获取该文件夹中所有bmp格式的图像
   
5. img_num = length(img_path_list);  % 获取图像总数量
   
6. if img_num > 0                    % 有满足条件的图像
   
7.     k=1;
   
8.     for j = 1 : 2 : img_num            % 逐一读取图像
   
9.         image_name = img_path_list(j).name;   % 图像名
   
10.         fprintf('正在处理的图像： %d %s\n',j,strcat(file_path,image_name));  % 显示正在处理的图像名
    
11.         img =  imread(strcat(file_path,image_name));
    
12.         if size(img,3)>1
    
13.             img = rgb2gray(img);
    
14.         end
    
15.         img = imresize(img,[512,512]);  % 缩小处理
    
16.         bw=img>10;
    
17.         % 连通域检测，计算特征
    
18.         cc = bwconncomp(bw);         % 连通域操作
    
19.         s  = regionprops(cc, {'centroid','area'});  % 面积特征
    
20.         [A, id] = max([s.Area]);     % 面积最大
    
21.         bw(labelmatrix(cc)~=id)=0;   % 只保留最大二值化块
    
22.         bw = bwfill(bw,'holes');     % 填充孔洞
    
23.         img1 = immultiply(img,bw);   % 交运算
    
24.         D(:,:,k) = img1;
    
25.         k=k+1;
    
26.     end
    
27. end
    
28. %% CT重建
    
29. % x,y,z为空间坐标，x,y为xy轴坐标
    
30. % z为高度：z(:,:,1)=1， z(:,:,2)=2，……，z(:,:,100)=100
    
31. [x,y,z,D1] = reducevolume(D,[4,4,1]);
    
32. % D1 = smooth3(D1);  % 平滑处理
    
33. % 环形曲面
    
34. p1 = patch(isosurface(x,y,z,D1, 5,'verbose'),...
    
35.     'FaceColor','c','EdgeColor','none');
    
36. isonormals(x,y,z,D1,p1);
    
37. % 上下曲面盖
    
38. p2 = patch(isocaps(x,y,z,D1, 5),...
    
39.     'FaceColor','interp','EdgeColor','none');
    
40. view(3);
    
41. axis tight;
    
42. daspect([1,1,0.4]);
    
43. colormap(gray(100));
    
44. camlight;
    
45. lighting gouraud;
    

复制代码

参考：【1】MATLAB绘制3D隐函数曲面的方法总结