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

10 次查看（过去 30 天）
Elisa
Elisa 2024-8-1
回答： Madheswaran 2024-8-7
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_size

请先登录，再回答此问题。

回答（1 个）

Madheswaran
Madheswaran 2024-8-7
Hi Elisa,
As you have mentioned there are duplicates in your data given in the files 'clu_sli_1400.mat' and 'I_1400.mat'.
There are several ways to solve this issue. One effective method is to take the simple average of the non-unique values and fit the curve using the unique x-values and the corrosponding y-value mean. Here are the changes that I have made to 'xData' and 'yData' before fitting using the 'fit' function in 'createFit.m':
[xDataNew, ~, indx] = unique(xData);
yDataNew = accumarray(indx, yData, [], @mean);
[fitresult, gof] = fit( xDataNew, yDataNew, ft, opts );
As you can see in the image inserted below, the curve did not fit as expected.
There are other alternatives you can use to handle the duplicates instead of simple averaging, such as
  1. Weighted averaging: Instead of simple averaging, you can try weighted averaging by assigning weights for the duplicates using any criteria (e.g., frequence or any context-dependant criteria).
  2. Penalizing the duplicates: You can introduce some penalty to reduce the impact of duplicate data.
You may need to explore and experiment with different approaches which might be specific to the dataset you have attached to find the most suitable solution.
Hope this helps!

类别

在 Help Center 和 File Exchange 中查找有关 Smoothing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by