BlackPianoCat
commented
3 years ago @AnonymFish If you want only Gaussian Surfaces simulation, this code is very good for that. If you want non-Gaussian surfaces simulation you need Pearson distribution, which I am not sure if it is the same as the distribution of MATLAB. If you use JohnsonSU distribution, the algorithm doesn't always converge. I tried to run it in python with both JohnsonSU and Pearson (the Pearson of python is not the same function as the Pearson of MATLAB), but it doesn't converge.
%Created: 2020.03.08
%Generation of non-Gaussian surfaces with predetermined
%1. Mean value (mu)
% 2. Standard deviation (rms)
% 3. Skewness (skew)
% 4. Kurtosis (kurt)
% 5. Correlation lengths (ksix, ksiy)
% 6. Roughness exponent (aplha)
% The input values of kurtosis and skewness should satisfy the inequality
% kurtosis > 1 + skewness^2
% The method is described in Yang et al. CMES, vol.103, no.4, pp.251-279,2014
% We use pearson translation to transform random noise to non-gaussian series and
% inverse Fourier transform to generate Gaussian surfaces.
% The method is implemented in three steps
% 1. Generation of a Gaussian self-affine surface z_gs with rms, ksix,ksiy,alpha
% 2. Generation of non-Gaussian noise z_ngn with mu,rms,skew,kurt
% 3. Arranging z_ngn according to the spatial arrangement of z_g to get a non-Gaussian self-affine surface z_ngs
function [z_ngs]=SAimage_fft_2(N_total,rms,skewness,kurtosis,corlength_x,corlength_y,alpha)
%INPUT
tic;
N=N_total; %surface side N
numpoints=N;
for n=1:N/2,
hhcf_y(n)=0.0;
end
Pyy=0.0;
% open the output files
foutsur=fopen('surface.txt','w');
% input moments and spatial parameters
mu=0;
rms=rms;
skew=skewness;
kurt=kurtosis;
ksix=corlength_x;
ksiy=corlength_y;
alpha=alpha;
%1st step: Generation of a Gaussian surface
%Determine the autocorrelation function R(tx,ty)
R(1:(numpoints+1),1:(numpoints+1))=0;
txmin=-numpoints/2;
txmax=numpoints/2;
dtx=(txmax-txmin)/(numpoints);
tymin=-numpoints/2;
tymax=numpoints/2;
dty=(tymax-tymin)/(numpoints);
tx=[txmin:dtx:txmax];
ty=[tymin:dty:tymax];
for tx=txmin:dtx:txmax
for ty=tymin:dty:tymax
% R(tx+txmax+1,ty+tymax+1)=((w^2)*exp(-(abs(tx)/ksi).^(2*a)))*((w^2)*exp(-(abs(ty)/ksi).^(2*a)));
R(tx+txmax+1,ty+tymax+1)=((rms^2)*exp(-(abs(sqrt((tx/ksix)^2+(ty/ksiy)^2))).^(2*alpha)));
% R(tx+txmax+1,ty+tymax+1)=((rms^2)*exp(-(abs(sqrt((tx/ksix)^2+(ty/ksiy)^2))).^(2*alpha)))*besselj(0,(2*pi*sqr(tx^2+ty^2))/lamda);
end
end
%According to the Wiener-Khinchine theorem FR is the power spectrum of the
%desired profile
FR = fft2(R,numpoints,numpoints);
AMPR=sqrt(dtx^2+dty^2)*abs(FR);
% 2nd step: Generate a white noise, normalize it and take its Fourier transform
X=rand(numpoints,numpoints);
aveX=mean(mean((X)));
dif2X=(X-aveX).^2;
stdX=sqrt(mean(mean(dif2X)));
X=X/stdX;
XF = fft2(X,numpoints,numpoints);
% 3nd step: Multiply the two Fourier transforms
YF=XF.*sqrt(AMPR);
% 4th step: Perform the inverse Fourier transform of YF and get the desired
% surface
zaf = ifft2(YF,numpoints,numpoints);
z=real(zaf);
% figure,surf(z);
avez=mean(mean(z));
dif2z=(z-avez).^2;
stdz=sqrt(mean(mean(dif2z)));
z=((z-avez)*rms)/stdz;
% Define the fraction of the surface to be analysed
xmin=1;
xmax=N;
ymin=1;
ymax=N;
z_gs=z(xmin:xmax,ymin:ymax);
Nh=xmax-xmin+1;
% figure, surf(z_gs);
% himage1=mat2gray(h1,[min(min(h1)),max(max(h1))]);
image_gs=mat2gray(z_gs);
% figure, imshow(image_gs)
clear R AMPR z
% return;
%2nd step: Generation of a non-Gaussian noise NxN
z_ngn=pearsrnd(mu,rms,skew,kurt,N,N);
%3rd step: Combination of z_gs with z_ngn to output a z_ngs
v_gs=z_gs(:);
v_ngn=z_ngn(:);
[vs_gs,Igs]=sort(v_gs);
[vs_ngn,Ingn]=sort(v_ngn);
for iv=1:N*N
ivs=Igs(iv);
v_ngs(ivs)=vs_ngn(iv);
end
z_ngs=reshape(v_ngs,N,N)';
figure, surf(z_ngs,'linestyle','none')
image_ngs=mat2gray(z_ngs);
bw=imbinarize(image_ngs);
figure, imshow(bw);
return;
% 1. CORRELATION FUNCTIONS
%a. 1-D height-height correlation function
inpsur=z_ngs;
for ndif=1:N/2,
surf1=inpsur(1:N,1:(N-ndif));
surf2=inpsur(1:N,(ndif+1):N);
difsur2=(surf1-surf2).^2;
hhcf1d(ndif)=sqrt(mean(mean(difsur2)));
rdif(ndif)=ndif;
end
%subplot(2,2,2);
figure,
loglog(rdif,hhcf1d);
grid on;
xlabel('r(nm)');
ylabel('G(r) (nm)');
title('1-D height-height correlation function');
toc;
fclose(foutsur);Let me know if you succeeded to solve the problem with non-Gaussian surfaces in python.
By the way, I am not a professor. I am just a master of Mathematical Modeling.
AnonymFish
Dear prof, I am a freshman in this area, so may i acess to your matlab code about simulating? Thanks so much for your help.