图像CT三维局部重建2

Halcom

图像CT序列，提取每一张图像的某一部分，将该部分重建。
视频链接：http://pan.baidu.com/s/1qYG8ajM
录制的视频是算法底层原理讲解，底层代码实现，方便大家真正掌握算法实质，开发出更加出色的算法。录制视频的初衷是：避免读者朋友利用大把时间学习已有常见算法，本系列视频旨在让读者朋友快速深入了解这些常见算法原理，有更多的时间去研究更加高大上算法（价值）。

具体链接在halcom.cn论坛，联系人QQ：3283892722
该论坛是一个学习交流平台，我会逐一的和大家分享学习。
欢迎大家录制视频，并提交给我，我来设置视频，你可在论坛进行打赏分享。
视频专用播放器：http://halcom.cn/forum.php?mod=viewthread&tid=258&extra=page%3D1

运行环境：win7+32bit+matlab2014a

1. % 分割石头
   
2. clc,clear,close all
   
3. D0 = 5; % 阻止的频率点与频域中心的距离
   
4. n = 1;  % 阶次
   
5. file_path =  './试件/';% 图像文件夹路径
   
6. img_path_list = dir(strcat(file_path,'*.jpg')); % 获取该文件夹中所有bmp格式的图像
   
7. img_num = length(img_path_list);  % 获取图像总数量
   
8. if img_num > 0                    % 有满足条件的图像
   
9.     k=1;
   
10.     for j = 1 : 2 : img_num                   % 逐一读取图像
    
11.         image_name = img_path_list(j).name;   % 图像名
    
12.         fprintf('正在处理的图像： %d %s\n',j,strcat(file_path,image_name));  % 显示正在处理的图像名
    
13.         img =  imread(strcat(file_path,image_name));
    
14.         if size(img,3)>1
    
15.             img = rgb2gray(img);
    
16.         end
    
17.         img = imresize(img,[512,512]);  % 缩小处理
    
18.         bw=img>10;
    
19.         % 连通域检测，计算特征
    
20.         cc = bwconncomp(bw);         % 连通域操作
    
21.         s  = regionprops(cc, {'centroid','area'});  % 面积特征
    
22.         [A, id] = max([s.Area]);     % 面积最大
    
23.         bw(labelmatrix(cc)~=id)=0;   % 只保留最大二值化块
    
24.         bw = bwfill(bw,'holes');     % 填充孔洞
    
25.         img1 = immultiply(img,bw);   % 交运算
    
26. %         figure(1),imshow(img1); title('目标图像');
    
27.         
    
28.         img_med = medfilt2(img1);
    
29.         img_btw = freqfilter_btw_Hp(img_med, D0, n);     % 巴特沃斯高通滤波器
    
30. %         figure(2), imshow(img_btw,[]); title('巴特沃斯高通滤波图像');
    
31.         
    
32.         % 分割石头
    
33.         stone_bw = img_btw>=12;
    
34.         stone_bw = bwareaopen(stone_bw,100);  % 去掉面积小于100的小块
    
35. %         figure(3),        imshow(stone_bw); title('石头成分图像');
    
36.         
    
37.         % 分割空气
    
38.         air_bw = img_med<60 & img_med>0;   % 分割空气分割
    
39.         % 边缘不存在空气
    
40.         se = strel('disk',5);  % 形态学操作
    
41.         bw_air = imerode(bw, se);
    
42.         air_bw = immultiply(air_bw,bw_air);   % 交运算
    
43.         air_bw = bwareaopen(air_bw,50);    % 去掉面积小于50的小块
    
44.         air_bw = bwfill(air_bw,'holes');   % 填充空洞
    
45.         air_part = immultiply(img1,air_bw);   % 交运算
    
46. %         figure(4),        imshow(air_bw); title('空气成分图像');
    
47.         
    
48.         % 分割凝胶部分
    
49.         gel_bw = img_btw<12;                  % 分割凝胶部分
    
50.         gel_bw = immultiply(gel_bw,bw_air);   % 交运算
    
51.         gel_bw = gel_bw-air_bw;           % 去除空气部分
    
52.         gel_bw = bwareaopen(gel_bw,50);   % 去掉面积小于50的小块
    
53.         gel_part = immultiply(img1,gel_bw);   % 交运算
    
54. %         figure(5),        imshow(gel_bw); title('凝胶成分图像');
    
55.         
    
56.         % 修正石头
    
57.         stone_bwz = bw-gel_bw-air_bw;
    
58.         stone_bwz = bwareaopen(stone_bwz,100);   % 去掉面积小于100的小块
    
59.         stone_part = immultiply(img1,stone_bwz);   % 交运算
    
60. %         figure(6),        imshow(stone_bwz); title('石头成分图像');
    
61.         
    
62.         D(:,:,k) = img1;
    
63.         stone(:,:,k) = stone_part;
    
64.         gel(:,:,k) = gel_part;
    
65.         air(:,:,k) = air_part;
    
66.         k=k+1;
    
67.     end
    
68. end
    
69. save D.mat D stone gel air
    
70. 
    
71. %% 空气部分
    
72. clear all;
    
73. load('D.mat');
    
74. figure(10),
    
75. [x,y,z,air1] = reducevolume(air,[8,8,1]);
    
76. p1 = patch(isosurface(x,y,z,air1, 5,'verbose'),...
    
77.     'FaceColor','c','EdgeColor','none');
    
78. isonormals(x,y,z,air1,p1);
    
79. p2 = patch(isocaps(x,y,z,air1, 5),...
    
80.     'FaceColor','c','EdgeColor','none');
    
81. view(3); axis tight;
    
82. daspect([1,1,0.4]);
    
83. colormap('jet');
    
84. camlight;
    
85. lighting gouraud;
    
86. hold off
    
87. 
    
88. % 获取点坐标
    
89. [t, p, pxZoom, pyZoom, pzZoom] = get_3D_points(x,y,z,air1);
    
90. % triangulation
    
91. figure(103),plot3(p(:,1)*pxZoom,p(:,2)*pyZoom,p(:,3)*pzZoom,'b.','MarkerSize',3);
    
92. xlabel('x'); ylabel('y'); zlabel('z');
    
93. title('空气部分');
    
94. figure(104),hold on;
    
95. trisurf(t, p(:,1)*pxZoom,p(:,2)*pyZoom,p(:,3)*pzZoom, 'edgecolor','none')
    
96. view(-37.5,30);
    
97. % axis tight;
    
98. xlabel('x'); ylabel('y'); zlabel('z');
    
99. camlight;
    
100. lighting gouraud;
     
101. colormap('lines');
     
102. box on;
     
103. hold off;

复制代码

参考：

【1】图像CT三维整体重建1