fractal dimensions calcuation #2153
Unanswered
Zhanlf0226
asked this question in
Ask Anything
Replies: 0 comments
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
-
hi everyone,
I wanna calculate the mineral fractal dimensions in a mylonitic sample, strangely, the result is <1 by using three linear fit methods, it means somewhere wrong. Enclosed is my data and my code. Could you help me check it what's the problem? Thank you very much!
code:
%% Import Script for EBSD Data
% This script was automatically created by the import wizard. You should
% run the whole script or parts of it in order to import your data. There
% is no problem in making any changes to this script.
%% Specify Crystal and Specimen Symmetries
% crystal symmetry
CS = {...
'notIndexed',...
crystalSymmetry('321', [4.9 4.9 5.5], 'X||a*', 'Y||b', 'Z||c*', 'mineral', 'Quartz-new', 'color', [0.53 0.81 0.98]),...
crystalSymmetry('12/m1', [8.6 13 7.2], [90,116,90]degree, 'X||a', 'Y||b*', 'Z||c', 'mineral', 'Orthoclase', 'color', [0.56 0.74 0.56]),...
crystalSymmetry('12/m1', [8.3 13 7.1], [90,116.05,90]degree, 'X||a', 'Y||b*', 'Z||c', 'mineral', 'Anorthoclase', 'color', [0.85 0.65 0.13])};
% plotting convention
setMTEXpref('xAxisDirection','east');
setMTEXpref('zAxisDirection','outOfPlane');
%% Specify File Names
% path to files
pname = 'D:\matlab2023b\mtex-5.11.2\mtex-5.11.2\dong_sample_2024\GL2208-1';
% which files to be imported
fname = [pname '\trailleft.cpr'];
%% Import the Data
% create an EBSD variable containing the data
ebsd = EBSD.load(fname, CS, 'interface', 'crc', 'convertEuler2SpatialReferenceFrame');
%% Compute grains and smooth
[grains, ebsd.grainId] = calcGrains(ebsd, 'angle', 5*degree);
grains = grains(~grains.isBoundary);
grains = smooth(grains, 5);
figure
plot(grains)
title('Grains before removing not indexed data')
%% Plot shape factor
figure
plot(grains, grains.shapeFactor)
mtexColorbar('title','shape factor')
title('Shape Factor Mapping after removing not indexed data and boundary grains')
%%
% Extract grains of 'Quartz-new'
grainsQ = grains('Quartz-new');
% Check data range
disp(['Equivalent Radius range: ', num2str(min(grainsQ.equivalentRadius)), ' to ', num2str(max(grainsQ.equivalentRadius))]);
disp(['Perimeter range: ', num2str(min(grainsQ.perimeter)), ' to ', num2str(max(grainsQ.perimeter))]);
% Remove grains with zero or negative equivalent radius or perimeter
validIdx = (grainsQ.equivalentRadius > 0) & (grainsQ.perimeter > 0);
grainsQ = grainsQ(validIdx);
% Log-transform the data
logEquivalentRadius = log(grainsQ.equivalentRadius);
logPerimeter = log(grainsQ.perimeter);
% Plot scatter plot
figure
scatter(grainsQ.equivalentRadius, grainsQ.perimeter)
xlabel('Equivalent Radius')
ylabel('Perimeter')
% Set axes to logarithmic scale
set(gca, 'XScale', 'log', 'YScale', 'log')
xlim([10, 10^4])
ylim([10^2, 10^5])
% 1. Linear fit using polyfit
p = polyfit(logEquivalentRadius, logPerimeter, 1);
hold on
plot([10,10^4],exp(p(2) + p(1) * log([10,10^4])),'LineWidth',3)
hold off
% Set legend
legend('Data', 'Polyfit', 'Fitlm', 'Regress')
% Output fractal dimensions
fractalDimension_polyfit = p(1);
trailleft.zip
Beta Was this translation helpful? Give feedback.
All reactions