三维跨孔电磁波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:n
    if 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)-1
    begin = list(i);
    bend = list(i+1);
    len = bend - begin;
    np_o = np(begin:bend-1);
    np_y = y(begin);
    % 判断矩阵是在左边还是在右边；
    if np_y == 0
        Nper(i,1:len) = np_o;
    else
        Nper(i,end-len+1:end) = np_o;
    end
end

Nper = [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 
clc

load('*.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:1000
    data_new(:,:,i) = np;
end

xslice = [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 == 1
           imwrite(imind,cm,fname,'gif', 'Loopcount',inf);
       else
           imwrite(imind,cm,fname,'gif','WriteMode','append','DelayTime', 0.25);
       end
    end
end