三維跨孔電磁波CT數(shù)據(jù)可視化框架搭建

1、三維CT可視化結(jié)果

對(duì)三維CT反演結(jié)果進(jìn)行180°旋轉(zhuǎn)，動(dòng)態(tài)展示三維CT反演結(jié)果。

三維CT平面圖

三維測(cè)線布置

CT數(shù)據(jù)解譯結(jié)果。

文章來(lái)源地址http://www.zghlxwxcb.cn/news/detail-694911.html

2、matlab代碼

2.1、CT數(shù)據(jù)格式整理并保存

close all
clear
clc

[inFileName,PathName] = uigetfile('*.txt',...
    '選擇CT數(shù)據(jù)文件','MultiSelect','on');

filename = strcat(PathName,inFileName);

% 將CT數(shù)據(jù)讀取出來(lái)；
data = importdata(filename);

% 將每一層的數(shù)據(jù)提取出來(lái)，以第二列數(shù)據(jù)為標(biāo)準(zhǔn)（同一深度）；
x = data(:,2);
% 找出每一層的起點(diǎn)；
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;

% 將深度記錄下來(lái)
np_depth = x(list(1:end-1));

% 找出每一行的最大列數(shù)；
row = max(data(:,1));
row = ceil(row);
% 一共有col行；

Nper = zeros(col,row);
% 將吸收系數(shù)值填充進(jìn)Nper矩陣；
% 注意吸收系數(shù)填充的位置，起點(diǎn)靠左還是靠右；
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];

% 保存吸收系數(shù)及深度
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);


% 將矩陣左右翻轉(zhuǎn)；
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','視吸收系數(shù)[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;  %設(shè)置仰角為30度。
% for az=0:1:1080  %讓方位角從0變到360，繞z軸一周
%     view(az,el);
%     drawnow;
% end

% az= 345;   %設(shè)置方位角為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

到了這里，關(guān)于三維跨孔電磁波CT數(shù)據(jù)可視化框架搭建的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！