Drawing planes in geology

Dear all

I need help to add an algoritm or modified the exsit one to get the figures annexed in file:the following is my code:

clc; clear; close all;
%% ---------------- 1. Load Data ------------------------------------
filename = 'ThreeFaultModel_Modified.xlsx';
data = xlsread(filename);

x = data(:,1);             
g_Bouguer = data(:,2);     
x = x(:); g_Bouguer = g_Bouguer(:);

%% ---------------- 2. User Inputs ---------------------------------
N_layers = input('Enter number of layers: ');

rho_profiles = zeros(1, N_layers);
for i = 1:N_layers
    rho_profiles(i) = input(['Density of layer ' num2str(i) ': ']);
end

rho_basement = input('Basement density: ');
rho_contrasts = rho_profiles - rho_basement;

%% ---------------- 3. Inversion ------------------------
G = 6.67e-3;

z_inv = bsxfun(@rdivide, g_Bouguer, (2*pi*G*rho_contrasts));

%% ---------------- 4. Detect ONE Main Fault ------------------------
VG = gradient(z_inv(:,1), x);

% اختيار أقوى تغير (fault)
[~, idx_fault] = max(abs(VG));

MainFault.x = x(idx_fault);
MainFault.dip = atan2d(VG(idx_fault),1);         % تقدير تقريبي
MainFault.throw = max(z_inv(:,1)) * 0.3;          % تقدير بسيط
MainFault.type = 'Normal';

if MainFault.dip > 0
    MainFault.type = 'Reverse';
end
if abs(MainFault.dip) > 85
    MainFault.type = 'Vertical';
end

%% ---------------- 5. Apply ONE Fault to All Layers ----------------
z_layers = z_inv;

for L = 1:N_layers
    if strcmp(MainFault.type,'Reverse')
        z_layers(idx_fault:end, L:end) = z_layers(idx_fault:end, L:end) + MainFault.throw;
        
    elseif strcmp(MainFault.type,'Normal')
        z_layers(idx_fault:end, L:end) = z_layers(idx_fault:end, L:end) - MainFault.throw;
        
    else
        % Vertical → no vertical shift
    end
end

%% ---------------- 6. Plot Depth ------------------------
figure; hold on;
for i = 1:N_layers
    plot(x, z_layers(:,i),'LineWidth',2);
end
set(gca,'YDir','reverse');
grid on;
title('Depth with Single Fault');
xlabel('Distance (km)');
ylabel('Depth (km)');
legend(arrayfun(@(k) sprintf('Layer %d',k),1:N_layers,'UniformOutput',false));

%% ---------------- 7. Plot Single Fault Plane ------------------------
z_max = max(z_layers(:));

if strcmp(MainFault.type,'Vertical')
    plot([MainFault.x MainFault.x],[0 z_max],'r--','LineWidth',2);
else
    dx = z_max / tand(MainFault.dip);
    
    if MainFault.dip > 0
        x2 = MainFault.x - dx;
    else
        x2 = MainFault.x + dx;
    end
    
    plot([MainFault.x x2],[0 z_max],'r--','LineWidth',2);
end

%% ---------------- 8. Dip Line from Surface ------------------------
figure; hold on;

for i = 1:N_layers
    plot(x, z_layers(:,i),'LineWidth',2);
end

set(gca,'YDir','reverse');
grid on;
title('Single Fault Dip Line');
xlabel('Distance (km)');
ylabel('Depth (km)');

z1 = max(z_layers(:,1));
dip = MainFault.dip;
x0 = MainFault.x;

if abs(dip) > 89
    x1 = x0;
else
    x1 = x0 + z1 / tand(dip);
end

plot([x0 x1],[0 z1],'r--','LineWidth',2);
text((x0+x1)/2,0.3,'Main Fault','Color','r');

%% ---------------- 9. Geological Section ------------------------
figure; hold on;

colors = lines(N_layers);

for i = N_layers:-1:1
    if i == 1
        top = zeros(size(x));
    else
        top = z_layers(:,i-1);
    end
    
    bot = z_layers(:,i);
    
    fill([x;flipud(x)], [top;flipud(bot)], colors(i,:), ...
        'EdgeColor','k','FaceAlpha',0.6);
end

% رسم الصدع
if strcmp(MainFault.type,'Vertical')
    plot([MainFault.x MainFault.x],[0 z_max],'k--','LineWidth',2);
else
    dx = z_max / tand(MainFault.dip);
    
    if MainFault.dip > 0
        x2 = MainFault.x - dx;
    else
        x2 = MainFault.x + dx;
    end
    
    plot([MainFault.x x2],[0 z_max],'k--','LineWidth',2);
end

set(gca,'YDir','reverse');
xlabel('Distance (km)');
ylabel('Depth (km)');
title('Geological Cross Section with Single Fault');

legend(arrayfun(@(i) sprintf('Layer %d',i),1:N_layers,'UniformOutput',false));

we use the following paramters:

Enter number of layers: 3

Density of layer 1: 2.20

Density of layer 2: 2.25

Density of layer 3: 2.30

Basement density: 2.67

[dpb]

%% ---------------- 1. Load Data ------------------------------------
filename = 'ThreeFaultModel_Modified.xlsx';
data = xlsread(filename);
x = data(:,1);             
g_Bouguer = data(:,2);     
x = x(:); g_Bouguer = g_Bouguer(:);

%% ---------------- 2. User Inputs ---------------------------------
N_layers = 3;                           % input('Enter number of layers: ');
rho_profiles=[2.20 2.25 2.30];          % zeros(1, N_layers);
                                        %    for i = 1:N_layers
                                        %        rho_profiles(i) = input(['Density of layer ' num2str(i) ': ']);
                                        %    end

rho_basement =  2.67;                   % input('Basement density: ');
rho_contrasts = rho_profiles - rho_basement;

%% ---------------- 3. Inversion ------------------------
G = 6.67e-3;
z_inv = bsxfun(@rdivide, g_Bouguer, (2*pi*G*rho_contrasts));

%% ---------------- 4. Detect ONE Main Fault ------------------------
VG = gradient(z_inv(:,1), x);

% اختيار أقوى تغير (fault)
[~, idx_fault] = max(abs(VG));
MainFault.x = x(idx_fault);
MainFault.dip = atan2d(VG(idx_fault),1);         % تقدير تقريبي
MainFault.throw = max(z_inv(:,1)) * 0.3;          % تقدير بسيط
MainFault.type = 'Normal';
if MainFault.dip > 0
    MainFault.type = 'Reverse';
end
if abs(MainFault.dip) > 85
    MainFault.type = 'Vertical';
end

%% ---------------- 5. Apply ONE Fault to All Layers ----------------
z_layers = z_inv;
for L = 1:N_layers
    if strcmp(MainFault.type,'Reverse')
        z_layers(idx_fault:end, L:end) = z_layers(idx_fault:end, L:end) + MainFault.throw;
    elseif strcmp(MainFault.type,'Normal')
        z_layers(idx_fault:end, L:end) = z_layers(idx_fault:end, L:end) - MainFault.throw;
    else
        % Vertical → no vertical shift
    end
end

%% ---------------- 6. Plot Depth ------------------------
figure; hold on;
for i = 1:N_layers
    plot(x, z_layers(:,i),'LineWidth',2);
end
set(gca,'YDir','reverse');
grid on;
title('Depth with Single Fault');
xlabel('Distance (km)');
ylabel('Depth (km)');
legend(arrayfun(@(k) sprintf('Layer %d',k),1:N_layers,'UniformOutput',false));

%% ---------------- 7. Plot Single Fault Plane ------------------------
z_max = max(z_layers(:));
if strcmp(MainFault.type,'Vertical')
    plot([MainFault.x MainFault.x],[0 z_max],'r--','LineWidth',2);
else
    dx = z_max / tand(MainFault.dip);
    if MainFault.dip > 0
        x2 = MainFault.x - dx;
    else
        x2 = MainFault.x + dx;
    end
    plot([MainFault.x x2],[0 z_max],'r--','LineWidth',2);
end

%% ---------------- 8. Dip Line from Surface ------------------------
figure; hold on;
for i = 1:N_layers
    plot(x, z_layers(:,i),'LineWidth',2);
end
set(gca,'YDir','reverse');
grid on;
title('Single Fault Dip Line');
xlabel('Distance (km)');
ylabel('Depth (km)');
z1 = max(z_layers(:,1));
dip = MainFault.dip;
x0 = MainFault.x;
if abs(dip) > 89
    x1 = x0;
else
    x1 = x0 + z1 / tand(dip);
end
plot([x0 x1],[0 z1],'r--','LineWidth',2);
text((x0+x1)/2,0.3,'Main Fault','Color','r');

%% ---------------- 9. Geological Section ------------------------
figure; hold on;
colors = lines(N_layers);
for i = N_layers:-1:1
    if i == 1
        top = zeros(size(x));
    else
        top = z_layers(:,i-1);
    end
    bot = z_layers(:,i);
    fill([x;flipud(x)], [top;flipud(bot)], colors(i,:), ...
        'EdgeColor','k','FaceAlpha',0.6);
end
% رسم الصدع
if strcmp(MainFault.type,'Vertical')
    plot([MainFault.x MainFault.x],[0 z_max],'k--','LineWidth',2);
else
    dx = z_max / tand(MainFault.dip);
    if MainFault.dip > 0
        x2 = MainFault.x - dx;
    else
        x2 = MainFault.x + dx;
    end
    plot([MainFault.x x2],[0 z_max],'k--','LineWidth',2);
end
set(gca,'YDir','reverse');
xlabel('Distance (km)');
ylabel('Depth (km)');
title('Geological Cross Section with Single Fault');
legend(arrayfun(@(i) sprintf('Layer %d',i),1:N_layers,'UniformOutput',false));

So, what is the specific Q?

[Star Strider]

Please save the figure in .png format instead, and describe what you want to do in some detail. (The .jpg format image, even when enlarged -- and it has to be enlarged -- is so indistinct that it is impossible to read it, even if it were in a European language using the appropriate alphabet.)

[Moustafa]

Please see the attached file I need code makes what should be in (After) under each figure after you after processing figures I made it manually to understand the idea

[Torsten]

I can't understand what you try to do. Please post the translation of the google translator from your native language - your English skills are a problem.

[dpb]

It appears the issue is to draw the best-fit lines at the fault zone boundaries.

@Moustafa, if you have the Signal Processing TB, look at <findchangepts> to locate the breakpoints in each trace. Then use the means of those for each set to solve for the best fit line with polyfit and can then evaluate those and plot over whatever range want.

[Steven Lord]

You could also use the ischange function in MATLAB, possibly in conjunction with the detrend or trenddecomp functions or the Find and Remove Trends task.

[dpb]

Actually, with these data, it may be that simply

% simulate and plot a representative shift 
z=[2*ones(1,3) 0.5*ones(1,5) 3*ones(1,8)];
x=0:numel(z)-1;
hAx=subplot(2,1,1);
plot(x,z,'linewidth',2)
ylim([0 5])
hAx.YDir='reverse';

% find the breakpoints
dzz=diff(sign(diff(z)));    % take first, second difference, scale to +/-1 range w/ sign()
iChangePts=find(dzz);       % the breakpoints are in pairs 

ix = 1×4

2 3 7 8

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

hold on
plot(ix,2,'xr','markersize',20) % show where are on axes

My previous comment about using +/-2 is appropriate logic for finding crossing points -- here it finds the extent of the section but not the actual inflection point pairs at each end...

[Moustafa]

فى الملف المرفق رسم لنتائج الكود قبل عمل التعديلات (الرسومات المعدلة يدويا باسفل كل ؤسمة من النتائج) المطلوب عملها بكود ماتلاب بطريقة تلقائية

The attached file contains a graph of the code results before modifications (manually modified graphs are below each result tag) that need to be performed automatically using MATLAB code.

[Torsten]

Well, I didn't mean that you post your question in Arabic, but that you enter your question in the google translator in Arabic and post the English translation. But I did this for you ...

But I think it will be hard to code automatically the modifications to be made after processing the data.

[dpb]

" it will be hard to code automatically the modifications"

Especially since there's no definition of what the modifications wanted are--I surmised the first set of lines is a best fit of the breakpoints which wouldn't appear to be terribly difficult.

However, the lines shown towards the right hand side are certainly puzzling as to how those were generated in the positions in which they are shown, certainly it won't be possible to write an algorithm without a specification.

Answer 1

[Mathieu NOE]

well , hello again @Moustafa

seems I already answered this question last year : How to shift lines to their correction positions (I need to correct figure) - MATLAB Answers - MATLAB Central

and the code I provided is still working IMHO :

clc; close all; clear;
% Load the Excel data
filename = 'ThreeFaultModel_Modified.xlsx';
data = xlsread(filename, 'Sheet1');
% Extract relevant columns
x = data(:, 1); % Distance (x in km)
Field = data(:, 2); % Earth filed (Field in unit)
%==============================================================%
% % Input number of layers and densities
% num_blocks = input('Enter the number of blocks: ');
% block_densities = zeros(1, num_blocks);
% for i = 1:num_blocks
%     block_densities(i) = input(['Enter the density of block ', num2str(i), ' (kg/m^3): ']);
% end
% EDIT to run code here
num_blocks = 5;
block_densities = [2.40, 2.45, 2.50, 2.55, 2.60];
%==============================================================%
% Constants
G = 0.00676; 
Lower_density = 2.67; % in kg/m^3
%==============================================================%
% Calculate inverted depth profile for each layer
z_inv = zeros(length(x), num_blocks);
for i = 1:num_blocks
    density_contrast = block_densities(i) - Lower_density;
    if density_contrast ~= 0
        z_inv(:, i) = Field ./ (2 * pi * G * density_contrast);
    else
        z_inv(:, i) = NaN; % Avoid division by zero
    end
end
%==============================================================%
% Compute vertical gradient (VG) of inverted depth (clean)
VG = diff(z_inv(:, 1)) ./ diff(x);
%==============================================================%
% Set fault threshold and find f indices based on d changes
f_threshold = 0.5; % Threshold for identifying significant d changes
f_indices = find(abs(diff(z_inv(:, 1))) > f_threshold);
%==============================================================%
% Initialize f locations and dip arrays
%==============================================================%
f_locations = x(f_indices);  % Automatically determined f locations
f_dip_angles = nan(size(f_indices)); % Placeholder for calculated dip 
% Calculate dip  for each identified f
for i = 1:length(f_indices)
    idx = f_indices(i);
    if idx < length(x)
        f_dip_angles(i) = atand(abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / (x(idx + 1) - x(idx)));
    else
        f_dip_angles(i) = atand(abs(z_inv(idx, 1) - z_inv(idx - 1, 1)) / (x(idx) - x(idx - 1)));
    end
end
%==============================================================%
%                    Displacement of faults 
%==============================================================%
D_faults = zeros(size(f_dip_angles));
for i = 1:length(f_indices)
    idx = f_indices(i);
    dip_angle_rad = deg2rad(f_dip_angles(i)); % Convert dip to radians
    D_faults(i) = abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / sin(dip_angle_rad);
end
% Assign displacement values
D1 = D_faults(1); % NF displacemen
D2 = D_faults(2); % VF displacement
D3 = D_faults(3); % RF displacement
%==============================================================%
%                 Processing Data for Interpretation
%==============================================================%
A = [x Field z_inv];  % New Data Obtained
col_names = {'x', 'Field'};
for i = 1:num_blocks
    col_names{end+1} = ['z', num2str(i)];
end
dataM = array2table(A, 'VariableNames', col_names);
t1 = dataM;
[nr, nc] = size(t1);
t1_bottoms = t1;
for jj = 3:nc
    for ii = 1:nr-1
        if t1_bottoms{ii, jj} ~= t1_bottoms{ii+1, jj}
            t1_bottoms{ii, jj} = NaN;
        end
    end
end
%==============================================================%
%                     Identifying NaN rows
%==============================================================%
nans = isnan(t1_bottoms{:, 3:end});
nan_rows = find(any(nans, 2));
xc = t1_bottoms{nan_rows, 1}; % Corrected x-coordinates
yc = zeros(numel(nan_rows), 1); % y-coordinates for NaN rows
for ii = 1:numel(nan_rows)
    idx = find(~nans(nan_rows(ii), :), 1, 'last');
    if isempty(idx)
        yc(ii) = 0;
    else
        yc(ii) = t1_bottoms{nan_rows(ii), idx+2};
    end
end
%==============================================================%
% Plot the interpreted d profiles
figure(3)
% Plot the fault lines in black (matching the image)
plot(x, zeros(size(x)), 'k', 'LineWidth', 3) % top horizontal lines
set(gca, 'YDir', 'reverse')  % Reverse Y-axis for depth representation
box on
grid on
xlabel('Distance (km)'); 
ylabel('D (km)'); 
title('New Interpretation of Profile Data')
hold on
yl = [0 ceil(max(z_inv(:)))];  % Get Y-axis limits
% fault lines are estimated with polyfit
ind_nan = [];
for kk = 1:numel(f_indices)
    xtmp = x;
    xa = x(f_indices(kk));
    xb = x(f_indices(kk)+1);
    ya = z_inv(f_indices(kk),:);
    yb = z_inv(f_indices(kk)+1,:);
    p{kk} = polyfit([ya yb],[xa*ones(1,5) xb*ones(1,5)],1);
    xf = polyval(p{kk},yl);
    % Plot the fault lines in red (matching the image)
    plot(xf, yl, 'r', 'LineWidth', 3)
    ind1 = find(xtmp>min(xf)-1, 1 );
    ind2 = find(xtmp<max(xf)+1, 1, 'last' );
    ind_nan = [ind_nan [ind1: ind2]];
end
x(ind_nan) = [];
z_inv(ind_nan,:) = [];
% complete horizontal lines datas
xnew = x;
z_inv_new = z_inv;
for kk = 1:numel(f_indices)
    xf = polyval(p{kk},yl);
    ind1(kk) = find(x>min(xf)-1, 1 );
    ind2(kk) = find(x<max(xf)+1, 1, 'last' );
    ya = z_inv(ind1(kk),:);
    yb = z_inv(ind2(kk),:);
   
    % replace ends of z_inv data to match red lines 
    % define new end values 
    xa = polyval(p{kk},ya);
    ytmp = diag(ya);
    ytmp(abs(ytmp)<eps) = NaN;
    % add new end values
    xnew = [xnew; xa(:)];
    z_inv_new = [z_inv_new; ytmp];
    % define new start values 
    xb = polyval(p{kk},yb);
    ytmp = diag(yb);
    ytmp(abs(ytmp)<eps) = NaN;
    % add new start values
    xnew = [xnew; xb(:)];
    z_inv_new = [z_inv_new; ytmp]; 
   
end
%% finally , the plot !!
% have to plot each individual segment otherwise the plot is messy
[z_inv_unic,ia,ic] = uniquetol(z_inv_new(:),1e-6);
ii = isnan(z_inv_unic);
z_inv_unic(ii) = [];
ia(ii) = [];
for kk = 1:numel(z_inv_unic)
    [r,c] = find(z_inv_new == z_inv_unic(kk));
    for m = 1:numel(c)
        plot(xnew(r),z_inv_new(r,c),'-', 'LineWidth', 1)
    end
end

[Moustafa]

When I run the modified code got the folloqing message:

Undefined function 'uniquetol' for input arguments of type 'double'.

Error in Tako6M3 (line 161)

[z_inv_unic,ia,ic] = uniquetol(z_inv_new(:),1e-6);

[Torsten]

The MATLAB function "uniquetol" was introduced in 2015. If you have a MATLAB version < 2015, try

https://eidors3d.sourceforge.net/doc/index.html?eidors/overloads/uniquetol.html

clc; close all; clear;
% Load the Excel data
filename = 'ThreeFaultModel_Modified.xlsx';
data = xlsread(filename, 'Sheet1');
% Extract relevant columns
x = data(:, 1); % Distance (x in km)
Field = data(:, 2); % Earth filed (Field in unit)
%==============================================================%
% % Input number of layers and densities
% num_blocks = input('Enter the number of blocks: ');
% block_densities = zeros(1, num_blocks);
% for i = 1:num_blocks
%     block_densities(i) = input(['Enter the density of block ', num2str(i), ' (kg/m^3): ']);
% end
% EDIT to run code here
num_blocks = 5;
block_densities = [2.40, 2.45, 2.50, 2.55, 2.60];
%==============================================================%
% Constants
G = 0.00676; 
Lower_density = 2.67; % in kg/m^3
%==============================================================%
% Calculate inverted depth profile for each layer
z_inv = zeros(length(x), num_blocks);
for i = 1:num_blocks
    density_contrast = block_densities(i) - Lower_density;
    if density_contrast ~= 0
        z_inv(:, i) = Field ./ (2 * pi * G * density_contrast);
    else
        z_inv(:, i) = NaN; % Avoid division by zero
    end
end
%==============================================================%
% Compute vertical gradient (VG) of inverted depth (clean)
VG = diff(z_inv(:, 1)) ./ diff(x);
%==============================================================%
% Set fault threshold and find f indices based on d changes
f_threshold = 0.5; % Threshold for identifying significant d changes
f_indices = find(abs(diff(z_inv(:, 1))) > f_threshold);
%==============================================================%
% Initialize f locations and dip arrays
%==============================================================%
f_locations = x(f_indices);  % Automatically determined f locations
f_dip_angles = nan(size(f_indices)); % Placeholder for calculated dip 
% Calculate dip  for each identified f
for i = 1:length(f_indices)
    idx = f_indices(i);
    if idx < length(x)
        f_dip_angles(i) = atand(abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / (x(idx + 1) - x(idx)));
    else
        f_dip_angles(i) = atand(abs(z_inv(idx, 1) - z_inv(idx - 1, 1)) / (x(idx) - x(idx - 1)));
    end
end
%==============================================================%
%                    Displacement of faults 
%==============================================================%
D_faults = zeros(size(f_dip_angles));
for i = 1:length(f_indices)
    idx = f_indices(i);
    dip_angle_rad = deg2rad(f_dip_angles(i)); % Convert dip to radians
    D_faults(i) = abs(z_inv(idx + 1, 1) - z_inv(idx, 1)) / sin(dip_angle_rad);
end
% Assign displacement values
D1 = D_faults(1); % NF displacemen
D2 = D_faults(2); % VF displacement
D3 = D_faults(3); % RF displacement
%==============================================================%
%                 Processing Data for Interpretation
%==============================================================%
A = [x Field z_inv];  % New Data Obtained
col_names = {'x', 'Field'};
for i = 1:num_blocks
    col_names{end+1} = ['z', num2str(i)];
end
dataM = array2table(A, 'VariableNames', col_names);
t1 = dataM;
[nr, nc] = size(t1);
t1_bottoms = t1;
for jj = 3:nc
    for ii = 1:nr-1
        if t1_bottoms{ii, jj} ~= t1_bottoms{ii+1, jj}
            t1_bottoms{ii, jj} = NaN;
        end
    end
end
%==============================================================%
%                     Identifying NaN rows
%==============================================================%
nans = isnan(t1_bottoms{:, 3:end});
nan_rows = find(any(nans, 2));
xc = t1_bottoms{nan_rows, 1}; % Corrected x-coordinates
yc = zeros(numel(nan_rows), 1); % y-coordinates for NaN rows
for ii = 1:numel(nan_rows)
    idx = find(~nans(nan_rows(ii), :), 1, 'last');
    if isempty(idx)
        yc(ii) = 0;
    else
        yc(ii) = t1_bottoms{nan_rows(ii), idx+2};
    end
end
%==============================================================%
% Plot the interpreted d profiles
figure(3)
% Plot the fault lines in black (matching the image)
plot(x, zeros(size(x)), 'k', 'LineWidth', 3) % top horizontal lines
set(gca, 'YDir', 'reverse')  % Reverse Y-axis for depth representation
box on
grid on
xlabel('Distance (km)'); 
ylabel('D (km)'); 
title('New Interpretation of Profile Data')
hold on
yl = [0 ceil(max(z_inv(:)))];  % Get Y-axis limits
% fault lines are estimated with polyfit
ind_nan = [];
for kk = 1:numel(f_indices)
    xtmp = x;
    xa = x(f_indices(kk));
    xb = x(f_indices(kk)+1);
    ya = z_inv(f_indices(kk),:);
    yb = z_inv(f_indices(kk)+1,:);
    p{kk} = polyfit([ya yb],[xa*ones(1,5) xb*ones(1,5)],1);
    xf = polyval(p{kk},yl);
    % Plot the fault lines in red (matching the image)
    plot(xf, yl, 'r', 'LineWidth', 3)
    ind1 = find(xtmp>min(xf)-1, 1 );
    ind2 = find(xtmp<max(xf)+1, 1, 'last' );
    ind_nan = [ind_nan [ind1: ind2]];
end
x(ind_nan) = [];
z_inv(ind_nan,:) = [];
% complete horizontal lines datas
xnew = x;
z_inv_new = z_inv;
for kk = 1:numel(f_indices)
    xf = polyval(p{kk},yl);
    ind1(kk) = find(x>min(xf)-1, 1 );
    ind2(kk) = find(x<max(xf)+1, 1, 'last' );
    ya = z_inv(ind1(kk),:);
    yb = z_inv(ind2(kk),:);
   
    % replace ends of z_inv data to match red lines 
    % define new end values 
    xa = polyval(p{kk},ya);
    ytmp = diag(ya);
    ytmp(abs(ytmp)<eps) = NaN;
    % add new end values
    xnew = [xnew; xa(:)];
    z_inv_new = [z_inv_new; ytmp];
    % define new start values 
    xb = polyval(p{kk},yb);
    ytmp = diag(yb);
    ytmp(abs(ytmp)<eps) = NaN;
    % add new start values
    xnew = [xnew; xb(:)];
    z_inv_new = [z_inv_new; ytmp]; 
   
end
%% finally , the plot !!
% have to plot each individual segment otherwise the plot is messy
[z_inv_unic,ia,ic] = uniquetol_repl(z_inv_new(:),1e-6);
ii = isnan(z_inv_unic);
z_inv_unic(ii) = [];
ia(ii) = [];
for kk = 1:numel(z_inv_unic)
    [r,c] = find(z_inv_new == z_inv_unic(kk));
    for m = 1:numel(c)
        plot(xnew(r),z_inv_new(r,c),'-', 'LineWidth', 1)
    end
end

function [z,ii,jj] = uniquetol_repl(x,tol,varargin)
%UNIQUETOL Unique element within a tolerance.
%   [Y,I,J] = UNIQUETOL(X,TOL) is very similar to UNIQUE, but allows an
%   additional tolerance input, TOL. TOL can be taken as the total absolute
%   difference between similar elements. TOL must be a none negative
%   scalar. If not provided, TOL is assumed to be 0, which makes UNIQUETOL
%   identical to UNIQUE.
%
%   UNIQUETOL(...,'ROWS')
%   UNIQUETOL(...,'FIRST')
%   UNIQUETOL(...,'LAST')
%   These expressions are identical to the UNIQUE counterparts.
%
%   See also UNIQUE.
 
 
 if size(x,1) == 1, x = x(:); end
 if nargin < 2 || isempty(tol) || tol == 0
     [z,ii,jj] = unique(x,varargin{:});
     return;
 end
 [y,ii,jj] = unique(x,varargin{:});
 if size(x,2) > 1
     [~,ord] = sort(sum(x.^2,1),2,'descend');
     [y,io] = sortrows(y,ord);
     [~,jo] = sort(io);
     ii = ii(io);
     jj = jo(jj);
 end
 d = sum(abs(diff(y,1,1)),2);
 
 isTol = [true;d > tol];
 z = y(isTol,:);
 bin = cumsum(isTol); % [n,bin] = histc(y,z);
 jj = bin(jj);
 ii = ii(isTol);
 
end

[Moustafa]

You solve partially the problem with R2025a instead R2014b

[Torsten]

I don't understand your comment. "uniquetol" has been replaced by "uniquetol_repl" in the above code - thus it should work with R2014b. Or are there other problems asides "uniquetol" that prevent the above code from running under R2014b ?