clc;
close all;
clear;
global i_start;
global j_start;
global base;
global p;
% 1. INPUT INFORMATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i_start = 400; % x coordinate of the upper left image corner
j_start = 400; % y coordinate of the upper left image corner
base = 2; % factor used for up-scaling (marked as lambda_1 in the paper)
p = [0.05 : 0.05 : 0.95, 0.96 : 0.01: 1.05, 1.1 : 0.1: 3]'; % order of statistical moments
L = 60; % Size of the specimen [mm]
t = 10; % base^t = the highest resolution of the image (1024)
n = 16; % base^n = the highest resolution that corresponds to the minimal grain diameter (around 1 micron)
GL0 = 11750; % Threshold value for setting GL values to zero
criterion = 1.55; % criterion for grain-pixels
image_name = 'Horiz plane 4.tif'; % Name of the image file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. UPLOADING AND CREATING DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
GL = imread(image_name);
s = size(GL);
GL = GL(i_start: i_start - 1 + base^t, j_start : j_start - 1 + base^t);
GL = GL - GL0;
GL(find(GL < 0)) = 0;
GL = im2double(GL);
GL_mean = mean2(GL);
ro_ind = GL./GL_mean;
lambda_long = base.^[n:-1:0]; % resolution
l_pix = L./lambda_long; % pixel size
load('GSD_experimental.mat', 'GSD_experimental');
D_experim = GSD_experimental(:,1);
Pr_experim = GSD_experimental(:,2);
% Truncated measurements (min grain diamater is around 50 microns)
D_experim_trunc = D_experim(1:t);
Pr_experim_trunc = 100.*(Pr_experim(1:t) - Pr_experim(t))./(100 - Pr_experim(t));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. ALOCATION OF THE MEMORY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Pr_UMF = zeros(length(lambda_long), 1); % Probability function
c_gamma_theory = zeros(1, n); % Co-dimension function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. COMPUTATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K_p_emp = TM_Analysis_function(ro_ind, t, lambda_long, criterion); % Computation of the empirical moment scaling function
C1 = (K_p_emp(find(p == 1.01)) - K_p_emp(find(p == 0.99)))/0.02; % Numerical calculation of parameter C1
alpha = (K_p_emp(find(p == 1.02)) + K_p_emp(find(p == 0.98)))/(0.02^2*C1); % Numerical calculation of parameter alpha
alpha_p = (1 - 1/alpha)^(-1);
gamma_0 = -C1*alpha_p/alpha;
% Computation of the co-diemnsion function
for k = 1 : n
gamma = log(criterion)/log(lambda_long(k));
if alpha > 1
if gamma > gamma_0
c_gamma_theory(k) = C1*(gamma/(C1*alpha_p) + 1/alpha)^alpha_p;
else
c_gamma_theory(k) = 0;
end
end
if alpha < 1
if gamma < gamma_0
c_gamma_theory(k) = C1*(gamma/(C1*alpha_p) + 1/alpha)^alpha_p;
else
c_gamma_theory(k) = 1e15;
end
end
end
% Computation of the probability function - the maximal resolution is base^16
Pr_UMF(end) = 100;
Pr_UMF(1:n) = 100*(1 - lambda_long(1:n).^(-c_gamma_theory)./lambda_long(1).^(-c_gamma_theory(1)));
% Computation of the truncated probability function - the maximal resolution is base^10
Pr_UMF_trunc = 100*(Pr_UMF(end-t:end) - Pr_UMF(end-t))./(100 - Pr_UMF(end-t));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4. PLOTTING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Central part of the uploaded grey-scale image, that has base^n pixels
figure;
imagesc(ro_ind);
colormap('gray');
colorbar;
set(gca, 'FontWeight', 'bold');
set(gca, 'FontWeight', 'bold');
% Comparison between the experimental GSD data and the proposed model
figure
semilogx(D_experim, Pr_experim, 'ro', 'LineWidth', 3, 'MarkerSize', 10);
hold on;
semilogx(D_experim_trunc, Pr_experim_trunc, 'bs', 'LineWidth', 3, 'MarkerSize', 10);
hold on
semilogx(l_pix, Pr_UMF, 'k--', 'LineWidth', 3);
hold on;
semilogx(l_pix(end-t:end), Pr_UMF_trunc, 'k', 'LineWidth', 3);
xlabel('D [mm]');
ylabel('Percent finer [%]');
ylim([0 100]);
xlim([l_pix(1) L]);
grid on;