三维跨孔电磁波CT数据可视化框架搭建

利用matlab实现对跨孔电磁波CT实测数据反演，并搭建了三维CT数据可视化框架，可装填实测CT反演数据。

1、三维CT可视化结果

对三维CT反演结果进行180°旋转，动态展示三维CT反演结果。

三维CT平面图

三维测线布置

CT数据解译结果。

2、matlab代码

2.1、CT数据格式整理并保存

close all
clear
clc[inFileName,PathName] = uigetfile('*.txt',...'选择CT数据文件','MultiSelect','on');filename = strcat(PathName,inFileName);% 将CT数据读取出来；
data = importdata(filename);% 将每一层的数据提取出来，以第二列数据为标准（同一深度）；
x = data(:,2);
% 找出每一层的起点；
k = 1;
n = length(x);
for i = 2:nif x(i) ~= x(i-1)list(k) = i;k = k + 1;end
end
list(k) = n+1;
list = [1,list];
col = k;% 将深度记录下来
np_depth = x(list(1:end-1));% 找出每一行的最大列数；
row = max(data(:,1));
row = ceil(row);
% 一共有col行；Nper = zeros(col,row);
% 将吸收系数值填充进Nper矩阵；
% 注意吸收系数填充的位置，起点靠左还是靠右；
np = data(:,3);
y = data(:,1);
for i = 1:length(list)-1begin = list(i);bend = list(i+1);len = bend - begin;np_o = np(begin:bend-1);np_y = y(begin);% 判断矩阵是在左边还是在右边；if np_y == 0Nper(i,1:len) = np_o;elseNper(i,end-len+1:end) = np_o;end
endNper = [np_depth,Nper];% 保存吸收系数及深度
filename = strcat(inFileName(1:end-4),'.mat');
save(filename,'Nper');contourf(Nper,30);colormap(jet);
set(gca,'ydir','reverse');
axis equal;

2.2、三维可视化

close all
clear 
clcload('*.mat');
np1 = Nper(:,2:end);
np1_depth = Nper(:,1);
clear Nper
load('*.mat');
np2 = Nper(:,2:end);
np2_depth = Nper(:,1);
clear Nper
load('*.mat');
np3 = Nper(:,2:end);
np3_depth = Nper(:,1);% 将矩阵左右翻转；
np1 = fliplr(np1);% np1(np1 == 0) = nan;
% contourf(np1,100,'LineStyle','none');
% colormap(jet);colorbar;
% shading interp
% caxis([0.1,0.7]);
% set(gca,'ydir','reverse');
% axis equal;clear np
np = zeros(50,51);
np(1:50,1:18) = np1;
np(7:50,19:35) = np2;
np(5:48,36:51) = np3;np(np == 0) = nan;
contourf(np,80,'LineStyle','none');
colormap(hsv);colorbar;
set(get(colorbar,'title'),'string','视吸收系数[Nper/m]','fontsize',14);
shading flat
caxis([0,0.65]);
set(gca,'ydir','reverse');
xlabel('水平距离/m');
ylabel('深度/m');
axis equal;data_new = zeros(50,51,100);
for i = 1:1000data_new(:,:,i) = np;
endxslice = [10,40];
yslice = [];
zslice = 1:499:1000;
slice(data_new,xslice,yslice,zslice);
colormap('hsv');
h = colorbar;
set(get(h,'title'),'string','beta(Nper/m)');
caxis([0,0.65]);
colorbar('eastoutside');
shading interp
set(gca,'xticklabel',[]);
set(gca,'yticklabel',[]);
set(gca,'zticklabel',[]);
% axis off
alpha(0.8);
view(345,-15);
% 
spinningGIF('zk45-48.gif');% el=-45;  %设置仰角为30度。
% for az=0:1:1080  %让方位角从0变到360，绕z轴一周
%     view(az,el);
%     drawnow;
% end% az= 345;   %设置方位角为0
% for el=0:1:360*1000   %仰角从0变到360
%     view(az,el);
%     drawnow;
% end% spinningGIF(fname): makes a spinning GIF of the current plot and saves it
% Usage: make your 3D plot (using plot3(...) or scatter3(...) etc.) and
% then call SpinningGIF with the file name that you want
function spinningGIF(fname)
%     axis off
%     view(0,10)center = get(gca, 'CameraTarget');pos = get(gca, 'CameraPosition');radius = norm(center(1:2) - pos(1:2));angles = pi:0.02*pi:2*pi;for ii=1:length(angles)angle = angles(ii);set(gca, 'CameraPosition', [center(1) + radius * cos(angle),...center(2) + radius * sin(angle),...pos(3)]);drawnow;frame = getframe(1);im = frame2im(frame);[imind,cm] = rgb2ind(im,256);if ii == 1imwrite(imind,cm,fname,'gif', 'Loopcount',inf);elseimwrite(imind,cm,fname,'gif','WriteMode','append','DelayTime', 0.25);endend
end