You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
125 lines
3.7 KiB
125 lines
3.7 KiB
function Compare_smooth()
|
|
load('ani_fittingALL_pnorm5.mat')
|
|
term_choose=68;
|
|
aPhys=a_log{term_choose,2};bPhys=b_log{term_choose,2};cPhys=c_log{term_choose,2};
|
|
nz=size(aPhys,1);ny=size(aPhys,2);nx=size(aPhys,3);
|
|
|
|
% ParameterField(aPhys,bPhys,cPhys)
|
|
resolution=40;
|
|
tic
|
|
voxel= voxelField(aPhys,bPhys,cPhys,resolution);
|
|
toc
|
|
save voxel_afterSmooth40.mat voxel
|
|
% tic
|
|
% voxel=voxel_ele(aPhys,bPhys,cPhys,10,1,10);
|
|
% toc
|
|
% sum(voxel,'all')
|
|
end
|
|
|
|
function voxel=voxel_ele(aPhys,bPhys,cPhys,xe,ye,ze)
|
|
resolution=40;
|
|
h=1.0/resolution;
|
|
x2=h/2:h:1-h/2;
|
|
y2=h/2:h:1-h/2;
|
|
z2=h/2:h:1-h/2;
|
|
[x,y,z]=meshgrid(x2+xe-1,y2+ye-1,z2+ze-1);
|
|
[Fa,Fb,Fc]=interpField(aPhys,bPhys,cPhys);
|
|
a_sim = Fa(x,y,z);b_sim = Fb(x,y,z);c_sim = Fc(x,y,z);
|
|
p=cos(2*pi*(x-1/2))+ a_sim.*cos(2*pi*(y-1/2))+ b_sim.*cos(2*pi*(z-1/2))- c_sim;
|
|
voxel = p>=0;
|
|
|
|
a = aPhys(ze,ye,xe);b = bPhys(ze,ye,xe);c = cPhys(ze,ye,xe);
|
|
p2=cos(2*pi*(x-1/2))+ a.*cos(2*pi*(y-1/2))+ b.*cos(2*pi*(z-1/2))- c;
|
|
voxel2 = p2>=0;
|
|
|
|
mean(abs(voxel-voxel2),'all')
|
|
end
|
|
|
|
function voxel= voxelField(aPhys,bPhys,cPhys,resolution)
|
|
[Fa,Fb,Fc]=interpField(aPhys,bPhys,cPhys);
|
|
nz=size(aPhys,1);ny=size(aPhys,2);nx=size(aPhys,3);
|
|
|
|
h=1.0/resolution;
|
|
x2=h/2:h:nx-h/2;
|
|
y2=h/2:h:ny-h/2;
|
|
z2=h/2:h:nz-h/2;
|
|
[x,y,z]=meshgrid(x2,y2,z2);
|
|
a_sim = Fa(x,y,z);b_sim = Fb(x,y,z);c_sim = Fc(x,y,z);
|
|
p=cos(2*pi*(x-1/2))+ a_sim.*cos(2*pi*(y-1/2))+ b_sim.*cos(2*pi*(z-1/2))- c_sim;
|
|
|
|
voxel = p>=0;
|
|
|
|
end
|
|
|
|
|
|
function [Fa,Fb,Fc]=interpField(aPhys,bPhys,cPhys)
|
|
nz=size(aPhys,1);ny=size(aPhys,2);nx=size(aPhys,3);
|
|
l=0.49;
|
|
xx=-1:1:1;
|
|
[dx,dy,dz]=meshgrid(xx,xx,xx);
|
|
dX=[dx(:),dy(:),dz(:)].*l;
|
|
num_p=size(xx,2)^3;
|
|
|
|
count=1;
|
|
a_fun = zeros(nx*ny*nz*num_p,1);
|
|
b_fun = zeros(nx*ny*nz*num_p,1);
|
|
c_fun = zeros(nx*ny*nz*num_p,1);
|
|
in1=zeros(nx*ny*nz*num_p,3);
|
|
for i=1:nx
|
|
for j=1:ny
|
|
for k=1:nz
|
|
for tt=0:num_p-1
|
|
a_fun(count*num_p-tt,1) = aPhys(k,j,i);
|
|
b_fun(count*num_p-tt,1) = bPhys(k,j,i);
|
|
c_fun(count*num_p-tt,1) = cPhys(k,j,i);
|
|
in1(count*num_p-tt,:)=[i-1/2,j-1/2,k-1/2]+dX(tt+1,:);
|
|
end
|
|
count=count+1;
|
|
end
|
|
end
|
|
end
|
|
|
|
Fa=scatteredInterpolant(in1(:,1),in1(:,2),in1(:,3),a_fun,'natural','linear');
|
|
Fb=scatteredInterpolant(in1(:,1),in1(:,2),in1(:,3),b_fun,'natural','linear');
|
|
Fc=scatteredInterpolant(in1(:,1),in1(:,2),in1(:,3),c_fun,'natural','linear');
|
|
% Fa.Method = 'linear';Fb.Method = 'linear';Fc.Method = 'linear';
|
|
end
|
|
|
|
function ParameterField(aPhys,bPhys,cPhys)
|
|
nz=size(aPhys,1);ny=size(aPhys,2);nx=size(aPhys,3);
|
|
figure()
|
|
subplot(311)
|
|
imagesc(squeeze(aPhys(:,1,:))); axis equal; axis tight; axis off;colorbar;caxis([0.5,2])
|
|
subplot(312)
|
|
imagesc(squeeze(bPhys(:,1,:))); axis equal; axis tight; axis off;colorbar;caxis([0.5,2])
|
|
subplot(313)
|
|
imagesc(squeeze(cPhys(:,1,:))); axis equal; axis tight; axis off;colorbar;caxis([-5,1])
|
|
|
|
[Fa,Fb,Fc]=interpField(aPhys,bPhys,cPhys);
|
|
x2=0:0.1:nx;
|
|
y2=0.5;%y2=0:0.1:1;
|
|
z2=0:0.1:nz;
|
|
[x,y,z]=meshgrid(x2,y2,z2);
|
|
a_sim = Fa(x,y,z);%b_sim = Fb(x,y,z);c_sim = Fc(x,y,z);
|
|
|
|
% figure()
|
|
% subplot(311)
|
|
% imagesc(squeeze(a_sim(1,:,:)).'); axis equal; axis tight; axis off;colorbar;caxis([0.5,2])
|
|
% subplot(312)
|
|
% imagesc(squeeze(b_sim(1,:,:)).'); axis equal; axis tight; axis off;colorbar;caxis([0.5,2])
|
|
% subplot(313)
|
|
% imagesc(squeeze(c_sim(1,:,:)).'); axis equal; axis tight; axis off;colorbar;caxis([0.5,2])
|
|
|
|
x2=0.1;
|
|
[x,y,z]=meshgrid(x2,y2,z2);
|
|
a_sim2 = Fa(x,y,z);%b_sim = Fb(x,y,z);c_sim = Fc(x,y,z);
|
|
sum(abs(a_sim2-a_sim),'all')
|
|
% figure()
|
|
% subplot(311)
|
|
% imagesc(squeeze(a_sim(1,:,:)).'); axis equal; axis tight; axis off;colorbar;caxis([0.5,2])
|
|
% subplot(312)
|
|
% imagesc(squeeze(b_sim(1,:,:)).'); axis equal; axis tight; axis off;colorbar;caxis([0.5,2])
|
|
% subplot(313)
|
|
% imagesc(squeeze(c_sim(1,:,:)).'); axis equal; axis tight; axis off;colorbar;caxis([0.5,2])
|
|
|
|
end
|