August 2024 | Application Package Repository Telkom University

Month: August 2024

Adjust the fitting for non-unique values data (alternative to the fuction smoothing spline, createFit.f)

/ 2024-08-01 / Matlab

Good morning everyone,
I am trying to extract the fracture aperture from a tomography image. I have written the code for the extraction, and everything works using the fitting smoothing spline (createFit function attached) for clu_sli_12.mat and I_12.mat, but not on clu_sli_1400.mat and I_1400.mat in which we have non-unique data (presented directly in the code below). The variable clu_sli_1400 represents the point cloud used as reference for selecting data in the tomography slice (variable I_1400). We tried using other functions in matlab, without success. Can anybody help suggesting a specific function for fitting this data?
thanks in advance for any help!!!!!!!!
clear;
close all;
clc;
warning off;

load(‘I_1400.mat’)
load(‘clu_sli_1400.mat’)

voxel_size=70.69*10^-3;

%% Iteration on the slices
mean_widths = [];
std_widths = [];
min_widths = [];
max_widths = [];
areas = [];

% for ii = 12 %we have to iterate to ii = 12 and ii = 1400 (the last one not working properly on the fitting)

% fitting data
X = clu_sli(:,1);
Y = clu_sli(:,2);

% Create model fitting
[fitresult, gof] = createFit(X, Y);

% Calcolare i valori della curva di fitting
vv = ppval(fitresult.p, X);

% Calcolare la prima derivata della curva di fitting
d1 = ppval(fnder(fitresult.p, 1), X);

% Calcolare la normale alla curva di fitting
epsilon = 1e-10; % Evita la divisione per zero
nn = 1 ./ (d1 + epsilon);

% Calcolare gli angoli delle normali
theta = atan(-1 ./ d1); % Angolo della normale (in radianti)

% Calcolare le componenti della normale
normal_x = cos(theta); % Componente x della normale
normal_y = sin(theta); % Componente y della normale

% Calcolare l’ampiezza della crepa e visualizzare i segmenti
span = 50;
largh = zeros(size(clu_sli, 1), 1);

% Creare una figura con due subplot
figure;

% Primo subplot: Visualizzare la slice con i pallini rossi
subplot(1, 2, 1);
pcolor(I);
shading interp;
axis equal;
hold on;
scatter(clu_sli(:,1), clu_sli(:,2), ‘r’);
xlabel(‘X’);
ylabel(‘Y’);

% Secondo subplot: Visualizzare la slice con i segmenti verdi
subplot(1, 2, 2);
pcolor(I);
shading interp;
axis equal;
hold on;
scatter(clu_sli(:,1), clu_sli(:,2), ‘r’);
plot(X, vv, ‘b’, ‘LineWidth’, 2);
quiver(X, vv, normal_x, normal_y, 0.5, ‘k’, ‘LineWidth’, 1);

intensity_profiles = cell(size(clu_sli, 1), 1);
vet_in = zeros(size(clu_sli, 1), 2);
vet_fin = zeros(size(clu_sli, 1), 2);

for kk = 1:size(clu_sli, 1)
c = round(X(kk));
r = round(vv(kk));

% Genera i pixel lungo la normale
x_profile = round(c + (-span:span) * normal_x(kk));
y_profile = round(r + (-span:span) * normal_y(kk));

% Assicurarsi che i pixel siano all’interno dell’immagine
valid_idx = x_profile > 0 & x_profile <= size(I, 2) & y_profile > 0 & y_profile <= size(I, 1);
x_profile = x_profile(valid_idx);
y_profile = y_profile(valid_idx);

% Estrarre i valori del profilo dall’immagine
profilo = I(sub2ind(size(I), y_profile, x_profile));
intensity_profiles{kk} = profilo;

% Calcolare l’ampiezza della crepa
largh(kk) = conv_crepa(profilo);

% Calcolare i punti del segmento di larghezza
half_width = largh(kk) / 2;
segment_x = [c – half_width * normal_x(kk), c + half_width * normal_x(kk)];
segment_y = [r – half_width * normal_y(kk), r + half_width * normal_y(kk)];

vet_in(kk,1) = c – half_width * normal_x(kk);
vet_fin(kk,1) = c + half_width * normal_x(kk);

vet_in(kk,2) = r – half_width * normal_y(kk);
vet_fin(kk,2) = r + half_width * normal_y(kk);

% Visualizzare il segmento di larghezza
plot(segment_x, segment_y, ‘g’, ‘LineWidth’, 2);
end

hold off;

% Calcolare la media, la deviazione standard, il minimo e il massimo della larghezza della crepa
mean_width = mean(largh, ‘omitnan’);
std_width = std(largh, ‘omitnan’);
min_width = min(largh, [], ‘omitnan’);
max_width = max(largh, [], ‘omitnan’);

mean_widths = [mean_widths; mean_width];
std_widths = [std_widths; std_width];
min_widths = [min_widths; min_width];
max_widths = [max_widths; max_width];

% Calcolo dell’area all’interno delle curve usando trapz
% Creazione del poligono per calcolare l’area
all_x = [vet_in(:,1); flip(vet_fin(:,1))];
all_y = [vet_in(:,2); flip(vet_fin(:,2))];

% Rimuovi i NaN
valid_idx = ~isnan(all_x) & ~isnan(all_y);
all_x = all_x(valid_idx);
all_y = all_y(valid_idx);

% Assicurarsi che i punti formino un poligono chiuso
if length(all_x) > 1 && (all_x(1) ~= all_x(end) || all_y(1) ~= all_y(end))
all_x(end+1) = all_x(1);
all_y(end+1) = all_y(1);
end

% Calcolo dell’area del poligono usando trapz
area = trapz(all_x, all_y);
areas = [areas; area];

title(‘Area Between Curves per Slice’);

%% calcoli in mm

areas_mm = areas.* voxel_size^2
std_mm=std_widths.*voxel_size
mean_mm=mean_widths.*voxel_size
min_mm=min_widths.*voxel_size
max_mm=max_widths.*voxel_sizeGood morning everyone,
I am trying to extract the fracture aperture from a tomography image. I have written the code for the extraction, and everything works using the fitting smoothing spline (createFit function attached) for clu_sli_12.mat and I_12.mat, but not on clu_sli_1400.mat and I_1400.mat in which we have non-unique data (presented directly in the code below). The variable clu_sli_1400 represents the point cloud used as reference for selecting data in the tomography slice (variable I_1400). We tried using other functions in matlab, without success. Can anybody help suggesting a specific function for fitting this data?
thanks in advance for any help!!!!!!!!
clear;
close all;
clc;
warning off;