update

AttitudeToVector.m

@@ -0,0 +1,13 @@
+% name：AttitudToVector.m
+% description：get the vector of attitude
+% input：attitude[x,y,z,attribute] 
+% return：vector[x,y,z,dx,dy,dz]
+% author：Linze Du, Yongqiang Tong, Baoyi Zhang, Umair Khan, Lifang Wang and Hao Deng.  
+% version：ver1.0.0[2022.9.23]
+
+function [vector]=AttitudeToVector(attitude)
+    dipAngl=deg2rad(attitude(:,5));
+    dipOrin=deg2rad(attitude(:,4));
+    dipVct=[cos(dipAngl).*sin(dipOrin),cos(dipAngl).*cos(dipOrin),-sin(dipAngl)];
+    vector=[attitude(:,1:3),dipVct];
+end

BaseFunc.m

@@ -0,0 +1,19 @@
+
+% name：BsaeFunc.m
+% description：calculating bsaefuncton value
+% input：k（float）：approximate minima in calculation(useless in current version) 
+%         meshgrid [X,Y,Z]
+%         X_：constants points
+%         Y_：constants points
+%         Z_：constants points
+%         type（int）：set basefunction type(useless in current version)
+% return：basefunction value
+% author：Linze Du, Yongqiang Tong, Baoyi Zhang, Umair Khan, Lifang Wang and Hao Deng.
+% version：ver1.0.0[2022.9.23]
+function basefunc=BaseFunc(X,Y,Z,X_,Y_,Z_,type)
+
+basefunc=((sqrt((X-X_).^2 + (Y-Y_).^2 +(Z-Z_).^2 )).^3);
+% Gaussian
+% MQ
+% IMQ
+end

DipvectorToGradientvector.m

@@ -0,0 +1,18 @@
+% name：DipvectorToGradientvct.m
+% description：get the gradient vector from attitude vector
+% input：k（float）：attitude_vector[x,y,z,dx,dy,dz]
+% return：gradient_vector[x,y,z,gx,gy,gz]
+% author：Linze Du, Yongqiang Tong, Baoyi Zhang, Umair Khan, Lifang Wang and Hao Deng.  
+% version：ver1.0.0[2022.9.23]
+function [gradient_vector]=DipvectorToGradientvector(attitude_vector)
+
+    stkOrin=DipvectorToStrikvector(attitude_vector);%走向矢量
+    gradient_vector=attitude_vector;
+    gradient_vector(:,4:6)=cross(attitude_vector(:,4:6),stkOrin(:,4:6));%叉积求得梯度矢量
+
+end
+
+
+
+
+

DipvectorToStrikvector.m

@@ -0,0 +1,18 @@
+% name：DipvectorToStrikvector .m
+% description：get the strike vector from attitude vector
+% input：attitude_vector[x,y,z,dx,dy,dz]
+% return：strike_vector[x,y,z,gx,gy,gz]
+% author：Linze Du, Yongqiang Tong, Baoyi Zhang, Umair Khan, Lifang Wang and Hao Deng.  
+% version：ver1.0.0[2022.9.23]
+
+function strike_vector= DipvectorToStrikvector(attitude_vector)
+
+    dipOrin=attitude_vector(:,4:6);
+    dipOrin(:,end)=0;
+    stkOrin=cross(dipOrin,attitude_vector(:,4:6));
+    Norm=sqrt(stkOrin(:,1).^2+stkOrin(:,2).^2);
+    stkNorm=[Norm,Norm,Norm];
+    stkOrin=stkOrin./stkNorm;
+    strike_vector=[attitude_vector(:,1:3),stkOrin];
+
+end

Experiment_AdaHRBF.m

@@ -0,0 +1,208 @@
+% loading fault attribute points
+attributes_fault = xlsread('Data\FaultChange2.xls');
+% loading digital elevation model
+DEM = xlsread('Data\DEM100m.xlsx');
+% loading stratigraphic attribute points in two sub-domains
+attributes_north = xlsread('Data\NStrataPoint.xls');
+attributes_south = xlsread('Data\SStrataPoint.xls');
+% loading stratigraphic attribute points in two sub-domains
+attitude_north = xlsread('Data\NFattitude.xlsx');
+attitude_south = xlsread('Data\SFattitude.xlsx');
+% approximate minima in partial derivative calculation 
+k=1e-4;
+% set basefunction type(cubic function) 
+basefunctype=1;
+
+% geting initial gradient magnitude from attitude points
+attitude_north=AttitudeToVector(attitude_north);
+attitude_south=AttitudeToVector(attitude_south);
+gradient_north=DipvectorToGradientvector(attitude_north);
+gradient_south=DipvectorToGradientvector(attitude_south);
+% geting the scope of study area and initial the mesh grid
+dataMax_n=max([attributes_north(:,[1,2,3]);attitude_north(:,[1,2,3]);attributes_south(:,[1,2,3]);attitude_south(:,[1,2,3])]);
+dataMin_n=min([attributes_north(:,[1,2,3]);attitude_north(:,[1,2,3]);attributes_south(:,[1,2,3]);attitude_south(:,[1,2,3])]);
+dataMax_s=max([attributes_north(:,[1,2,3]);attitude_north(:,[1,2,3]);attributes_south(:,[1,2,3]);attitude_south(:,[1,2,3])]);
+dataMin_s=min([attributes_north(:,[1,2,3]);attitude_north(:,[1,2,3]);attributes_south(:,[1,2,3]);attitude_south(:,[1,2,3])]);
+[meshgrid_NX,meshgrid_NY,meshgrid_NZ] = meshgrid(dataMin_n(1,1):100:dataMax_n(1,1),dataMin_n(1,2):100:dataMax_n(1,2),0:25:dataMax_n(1,3));
+[meshgrid_SX,meshgrid_SY,meshgrid_SZ] = meshgrid(dataMin_s(1,1):100:dataMax_s(1,1),dataMin_s(1,2):100:dataMax_s(1,2),0:25:dataMax_s(1,3));
+%%
+%HRBF: interpolating an initial 3D scalar field function in northern sub-domain
+[alph_n,bravo_n,charlie_n]=GetParameters(1e-4,1,attributes_north,gradient_north);
+valueGrid_north=GetValueGrid(k,basefunctype,meshgrid_NX,meshgrid_NY,meshgrid_NZ,attributes_north,gradient_north,alph_n,bravo_n,charlie_n);
+
+%HRBF: interpolating an initial 3D scalar field function in southern sub-domain 
+[alphs,betas,charlies]=GetParameters(1e-4,1,attributes_south,gradient_south);
+valueGrid_south=GetValueGrid(k,basefunctype,meshgrid_SX,meshgrid_SY,meshgrid_SZ,attributes_south,gradient_south,alphs,betas,charlies);
+%%
+%adaHRBF: interpolating an optimal 3D scalar field in northern sub-domain
+[opt_gradient_north,opt_alph_n,opt_bravo_n,opt_charlie_n]=OptGradientMagnitudes(attributes_north,gradient_north,100,1e-4,300);
+opt_valueGrid_north=GetValueGrid(k,basefunctype,meshgrid_NX,meshgrid_NY,meshgrid_NZ,attributes_north,opt_gradient_north,opt_alph_n,opt_bravo_n,opt_charlie_n);
+
+%adaHRBF: interpolating an optimal 3D scalar field in southern sub-domain
+[opt_gradient_south,opt_alph_s,opt_bravo_s,opt_charlie_s]=OptGradientMagnitudes(attributes_south,gradient_south,100,1e-4,300);
+opt_valueGrid_south=GetValueGrid(k,basefunctype,meshgrid_SX,meshgrid_SY,meshgrid_SZ,attributes_south,opt_gradient_south,opt_alph_s,opt_bravo_s,opt_charlie_s);
+%%
+%RBF: interpolating a 3D scalar field for fault
+[alph,~,charlie]=GetParameters(1e-4,1,attributes_fault);
+valueGrid_fault=GetValueGrid(k,basefunctype,meshgrid_NX,meshgrid_NY,meshgrid_NZ,attributes_fault,nan,alph,nan,charlie);
+%%
+% removing unnecessary grids
+index_n=find(valueGrid_fault>=0);
+index_s=find(valueGrid_fault<0);
+
+valueGrid_north(index_n)=nan;
+valueGrid_south(index_s)=nan;
+
+opt_valueGrid_north(index_n)=nan;
+opt_valueGrid_south(index_s)=nan;
+
+index=find(meshgrid_NY<=-1.21*meshgrid_NX+3379598.04 | meshgrid_NY>=-1.21*meshgrid_NX+3389045.69 | meshgrid_NY>=0.83*meshgrid_NX+2028031.65 | meshgrid_NY<=0.83*meshgrid_NX+2011701.6);
+valueGrid_north(index)=nan;
+valueGrid_fault(index)=nan;
+opt_valueGrid_north(index)=nan;
+index=find(meshgrid_SY<=-1.21*meshgrid_SX+3379598.04 | meshgrid_SY>=-1.21*meshgrid_SX+3389045.69 | meshgrid_SY>=0.83*meshgrid_SX+2028031.65 | meshgrid_SY<=0.83*meshgrid_SX+2011701.6);
+valueGrid_south(index)=nan;
+valueGrid_fault(index)=nan;
+opt_valueGrid_south(index)=nan;
+
+% removing grids outside the study area in the northern subdomain
+bar = waitbar(0,'Removing grids outside the northern study area...'); 
+for i=1:size(meshgrid_NZ,1) 
+    for j=1:size(meshgrid_NZ,2)
+        for k=1:size(meshgrid_NZ,3)
+            high = find(DEM(:,1)==meshgrid_NX(i,j,k) & DEM(:,2) == meshgrid_NY(i,j,k));
+            if (meshgrid_NZ(i,j,k)>=DEM(high,3) | meshgrid_NZ(i,j,k)<0)
+                valueGrid_north(i,j,k) = nan;
+                opt_valueGrid_north(i,j,k) = nan;
+            end
+        end
+    end
+    waitbar(i/(size(meshgrid_NZ,1)),bar,'Removing grids outside the northern study area...')    
+end
+close(bar)
+%removing grids outside the study area in the southern subdomain
+bar = waitbar(0,'Removing grids outside the southern study area...'); 
+for i=1:size(meshgrid_SZ,1)
+    for j=1:size(meshgrid_SZ,2)
+        for k=1:size(meshgrid_SZ,3)
+            high = find(DEM(:,1)==meshgrid_SX(i,j,k) & DEM(:,2) == meshgrid_SY(i,j,k));
+            if (meshgrid_SZ(i,j,k)>=DEM(high,3) | meshgrid_SZ(i,j,k)<0)
+                valueGrid_south(i,j,k) = nan;
+                opt_valueGrid_south(i,j,k) = nan;
+            end
+        end
+    end
+    waitbar(i/(size(meshgrid_SZ,1)),bar,'Removing grids outside the southern study area...')    
+end
+close(bar)
+%%
+% dsiplaying
+value_list=[4872,7233,8429,9145,9674,10602,11717,8429,9145,9674,10602,11142];
+color_list={'#ffe7ff','#fedfff','#ecffb0','#e2e2e2','#d0cfe1','#c5c4ca','#ffdfc0','#ecffb0','#e2e2e2','#d0cfe1','#c5c4ca','#ffe7cf'};
+area=[662000 677000 2568000 2582000 0 1150];
+color_axis=[4000,13000];
+
+% dsiplaying scalar field in the northern subdomain
+fig1=figure('Name','HRBFscatters');
+figure(fig1);
+% scalar field with initial gradient magnitude
+subplot(1,2,1);
+caxis(color_axis)
+colorbar;
+axis(area);
+axis equal;
+set(gca,'XLim',[area(1) area(2)]);
+set(gca,'YLim',[area(3) area(4)]);
+set(gca,'ZLim',[area(5) area(6)]);
+scatter3(meshgrid_NX(:),meshgrid_NY(:),meshgrid_NZ(:),25,valueGrid_north(:),'.');
+hold on;
+% scalar field with optimized gradient magnitude
+subplot(1,2,2);
+caxis(color_axis)
+colorbar;
+axis(area);
+axis equal;
+set(gca,'XLim',[area(1) area(2)]);
+set(gca,'YLim',[area(3) area(4)]);
+set(gca,'ZLim',[area(5) area(6)]);
+scatter3(meshgrid_NX(:),meshgrid_NY(:),meshgrid_NZ(:),25,opt_valueGrid_north(:),'.');
+hold on;
+%%
+% displaying stratigraphic interfaces in the northern subdomain
+fig2=figure('Name','HRBFisosurfaces');         
+figure(fig2);
+% stratigraphic interfaces with initial gradient magnitude
+subplot(1,2,1);
+light('position',[-0.5,-0.5,0.5],'style','local','color','[0.5,0.5,0.5]')
+grid on; box on;
+axis(area);
+axis equal;
+set(gca,'XLim',[area(1) area(2)]);
+set(gca,'YLim',[area(3) area(4)]);
+set(gca,'ZLim',[area(5) area(6)]);
+for i=1:7
+    pic = patch(isosurface(meshgrid_NX,meshgrid_NY,meshgrid_NZ,valueGrid_north,value_list(i)));
+    isonormals(meshgrid_NX,meshgrid_NY,meshgrid_NZ,valueGrid_north,pic)
+    set(pic,'FaceColor',color_list{i},'EdgeColor','none');
+    hold on;
+end
+% dsiplaying faults interfaces
+pic_f = patch(isosurface(meshgrid_NX,meshgrid_NY,meshgrid_NZ,valueGrid_fault,0));
+isonormals(meshgrid_NX,meshgrid_NY,meshgrid_NZ,valueGrid_fault,pic_f)
+set(pic_f,'FaceColor','red','EdgeColor','none');
+
+% stratigraphic interfaces with optimized gradient magnitude
+subplot(1,2,2);
+light('position',[-0.5,-0.5,0.5],'style','local','color','[0.5,0.5,0.5]')
+grid on; box on;
+axis(area);
+axis equal;
+set(gca,'XLim',[area(1) area(2)]);
+set(gca,'YLim',[area(3) area(4)]);
+set(gca,'ZLim',[area(5) area(6)]);
+for i=1:7
+    pic = patch(isosurface(meshgrid_NX,meshgrid_NY,meshgrid_NZ,opt_valueGrid_north,value_list(i)));
+    isonormals(meshgrid_NX,meshgrid_NY,meshgrid_NZ,opt_valueGrid_north,pic)
+    set(pic,'FaceColor',color_list{i},'EdgeColor','none');
+    hold on;
+end
+% dsiplaying faults interfaces
+pic_f = patch(isosurface(meshgrid_NX,meshgrid_NY,meshgrid_NZ,valueGrid_fault,0));
+isonormals(meshgrid_NX,meshgrid_NY,meshgrid_NZ,valueGrid_fault,pic_f)
+set(pic_f,'FaceColor','red','EdgeColor','none');
+
+%%
+% dsiplaying scalar field in the southern subdomain
+figure(fig1);
+% scalar field with initial gradient magnitude
+subplot(1,2,1);
+scatter3(meshgrid_SX(:),meshgrid_SY(:),meshgrid_SZ(:),25,valueGrid_south(:),'.');
+% scalar field with optimized gradient magnitude
+subplot(1,2,2);
+scatter3(meshgrid_SX(:),meshgrid_SY(:),meshgrid_SZ(:),25,opt_valueGrid_south(:),'.');
+hold off;
+%%
+% dsiplaying stratigraphic interfaces in the southern subdomain
+figure(fig2);
+% stratigraphic interfaces with initial gradient magnitude
+subplot(1,2,1);
+for i=8:12
+    pic = patch(isosurface(meshgrid_SX,meshgrid_SY,meshgrid_SZ,valueGrid_south,value_list(i)));
+    isonormals(meshgrid_SX,meshgrid_SY,meshgrid_SZ,valueGrid_south,pic)
+    set(pic,'FaceColor',color_list{i},'EdgeColor','none');
+    hold on;
+end
+%stratigraphic interfaces with optimized gradient magnitude
+subplot(1,2,2);
+for i=8:12
+    pic = patch(isosurface(meshgrid_SX,meshgrid_SY,meshgrid_SZ,opt_valueGrid_south,value_list(i)));
+    isonormals(meshgrid_SX,meshgrid_SY,meshgrid_SZ,opt_valueGrid_south,pic)
+    set(pic,'FaceColor',color_list{i},'EdgeColor','none');
+    hold on;
+end
+hold off;
+
+
+
+
+

GetMatrixA.m

@@ -0,0 +1,104 @@
+% name：GetMatrixA.m
+% description：get the dividend of coefficients solving linear system
+% input：k（float）：approximate minima in calculation(useless in current version) 
+%         baseFuncType（int）：set basefunction type(useless in current version) 
+%         varargin{1}：[x,y,z,attribute]relatively butial depth for RBF,(ada)HRBF
+%         varargin{2}：[x,y,z,gx,gy,gz]gradient magnitude for (ada)HRBF
+% return：matrix A 
+% author：
+% version：ver1.0.0[2022.9.23]
+
+function A = GetMatrixA(k,baseFuncType,varargin)
+if nargin == 3
+% RBF
+    attribute_points = varargin{1};
+    clear varargin;
+elseif nargin == 4
+% HRBF
+    attribute_points = varargin{1};
+    gradient_points = varargin{2};
+    r_m=size(gradient_points,1);
+    clear varargin;
+else
+   error('GetMatrixA:The number of parameters is incorrect.')
+end
+
+r_num= size(attribute_points,1);
+
+%A（1，1）
+Phi=zeros(r_num,r_num);
+for i=1:r_num
+    xi=ones(r_num,1)*attribute_points(i,1);
+    yi=ones(r_num,1)*attribute_points(i,2);
+    zi=ones(r_num,1)*attribute_points(i,3);
+    Phi(i,1:r_num)=(BaseFunc(attribute_points(:,1),attribute_points(:,2),attribute_points(:,3),xi,yi,zi,baseFuncType))';
+end
+
+%A（1，3）
+Charlie=ones(r_num,4);
+Charlie(:,2:4)=attribute_points(1:r_num,1:3);
+
+if nargin == 3
+%construct matrix A in RBF
+    A = zeros(r_num+4,r_num+4);
+    A(1:r_num,1:r_num)=Phi;
+    A(1:r_num,r_num+1:r_num+4)=Charlie;
+    A(r_num+1:r_num+4,1:r_num)=Charlie';
+    return
+end
+
+%A（1，2）
+Nabla_Phi=zeros(r_num,3*r_m);
+for i=1:r_m    
+    xi=ones(r_num,1)*gradient_points(i,1);
+    yi=ones(r_num,1)*gradient_points(i,2);
+    zi=ones(r_num,1)*gradient_points(i,3);
+    [Nabla_Phi(:,3*i-2),Nabla_Phi(:,3*i-1),Nabla_Phi(:,3*i)]=GetPartialDerivative(attribute_points(:,1),attribute_points(:,2),attribute_points(:,3), xi,yi,zi,k,baseFuncType);
+end
+
+%A（2，1）
+Nabla_Phi_T=zeros(3*r_m,r_num);
+for i=1:r_m    
+    xi=ones(1,r_num)*gradient_points(i,1);
+    yi=ones(1,r_num)*gradient_points(i,2);
+    zi=ones(1,r_num)*gradient_points(i,3);
+    [Nabla_Phi_T(3*i-2,:),Nabla_Phi_T(3*i-1,:),Nabla_Phi_T(3*i,:)]=GetPartialDerivative( xi,yi,zi,attribute_points(:,1)',attribute_points(:,2)',attribute_points(:,3)',k,baseFuncType);
+end
+
+%A（2，2）
+Nabla2_Phi_T=zeros(3*r_m,3*r_m);
+for i=1:r_m
+    for j=1:r_m
+        if i~=j
+            e3=GetSecPartialDerivative(gradient_points(i,1),gradient_points(i,2),gradient_points(i,3),gradient_points(j,1),gradient_points(j,2),gradient_points(j,3),0.01,1);
+        else
+            e3=[0,0,0;0,0,0;0,0,0];
+        end 
+        Nabla2_Phi_T(3*(i-1)+(1:3),3*(j-1)+(1:3))=double(e3);
+    end
+end
+clear e3;
+
+%A（2，3）
+Nabla_Charlie=zeros(3*r_m,4);
+for i=1:3*r_m
+    one_pos=mod(i-1,3);
+    Nabla_Charlie(i,end+one_pos-2)=1;
+end
+
+if nargin == 4
+% construct matrix A in (ada）HRBF
+    A = zeros(r_num+3*r_m+4,r_num+3*r_m+4);
+    A(1:r_num,1:r_num)=Phi;
+    A(1:r_num,r_num+1:r_num+3*r_m)=Nabla_Phi;
+    A(1:r_num,r_num+3*r_m+1:r_num+3*r_m+4)=Charlie;
+    A((r_num+1):(r_num+3*r_m),1:r_num)=Nabla_Phi_T;
+    A((r_num+1):(r_num+3*r_m),(r_num+1):(r_num+3*r_m))= Nabla2_Phi_T;
+    A((r_num+1):(r_num+3*r_m),(r_num+3*r_m+1):(r_num+3*r_m+4))=Nabla_Charlie;
+    A((r_num+3*r_m+1):(r_num+3*r_m+4),1:r_num)=Charlie';
+    A((r_num+3*r_m+1):(r_num+3*r_m+4),(r_num+1):(r_num+3*r_m))=Nabla_Charlie';
+    return
+end
+
+end
+