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

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

Answer 1

[Madheswaran]

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!