Extraction of the spatial coordinates of a curve and calculation of the normal line to each of its points

Hello everyone

I have the following dataset: https://www.dropbox.com/scl/fo/3dthvtj2ex9gmhsu3lkrq/AJt2KVYls3YM89v5SRmSEWk?rlkey=r63d0ej2sv7ku8i4ardla5pvi&dl=0, which gives rise to the data on the panel mz of the figure below:

My idea is to extract the (a,b) set of points that define the whole right side of the whitish curve. Once that it is done, I will need to extract, for each of the points of this curve, the normal line along the track (that it is, the (a,b) coordinates together with the corresponding values of the mz values).

Any ideas on how to do this efficiently?

Additional information (September 3, 2024)

When I talk about the normal to the tangent at each point of the curve mz, I am thinking about something like this:

where the green X symbol corresponds to a certain point of the curve, the black line to the tangent to the curve at that specific (x,y) point, and the dashed yellow line refers to the normal to the tangent line that (x,y) point. I would be interested in extracting the mz value all along that normal line inside the boundaries of the sample. So I would say that some grid data process would be important to make if not too many points actually exist in the given data set.

Answer 1

[Mathieu NOE]

hello

tried a few things here

As far as I have understood correctly your question, the points you are looking for are given by xs,ys,mzs as computed below.

I interpreted this "whole right side of the whitish curve" as the points that have a slightly positive z value so there are slightly on the right side of the white separation line . Adapt the value (and tlerance) to your needs . From the extracted points I obtain a smooth curve using this Fex submission

smoothn - File Exchange - MATLAB Central (mathworks.com)

Zoomed plot showing the raw selected points (in green) and the smoothed curve (in black)

code :

NB : the last portion (commented) is probably not what you need , I have probably misinterpreted the term "normal" in your query

%% load data files 
aa = readmatrix('Space_a.txt');
bb = readmatrix('Space_b.txt');
mz = readmatrix('mz_Data.txt');

% get "white line" border points coordinates 
ztarget = 0.05;
tol = 0.02;
[r,c] =  find(abs(mz - ztarget)<tol);
xl = aa(c);
yl = bb(r);
zl = mz(c,r);

% sort  y data
[yl,ia] = sort(yl);
xl = xl(ia);
zl = zl(ia,:);

%% smoothed curve through the points 
% see Fex : https://www.mathworks.com/matlabcentral/fileexchange/25634-smoothn/
zs = smoothn({xl,yl});
xs =    zs{1};
ys =    zs{2};
% find corresponding z data
mzs = interp2(aa,bb,mz,xs,ys); % notice mzs is equal / very close to target z value (ztarget) which is what we wanted ! 

%% plot
figure(1)
imagesc(aa,bb,mz)
zmin = min(mz(:));
zmax = max(mz(:));
colormap(redwhiteblue(zmin, zmax));
colorbar;
set(gca,'YDir','normal');

hold on 
plot(xl,yl,'*g');
plot(xs,ys,'k','linewidth',2)
hold off

% %% find the normal vectors to the curve
% % see Fex : https://fr.mathworks.com/matlabcentral/fileexchange/32696-2d-line-curvature-and-normals?s_tid=ta_fx_results
% % Vertices = [xs';ys'];
% % Vertices = [xs ys];
% ind = (1:10:numel(xs));
% Vertices = [xs(ind) ys(ind)];
% 
% N=LineNormals2D(Vertices);
% plot([Vertices(:,1) Vertices(:,1)+10*N(:,1)]',[Vertices(:,2) Vertices(:,2)+10*N(:,2)]','c');
% 

[Roderick]

Hi @Mathieu NOE. Thank you very much! I have modified the question to include some additional information about the normal to the tangent at each (x,y) point of the obtained curve. Hope it helps.

[Mathieu NOE]

hello Richard

here some update... so we are now able to trace the normals (using this Fex : 2D Line Curvature and Normals - File Exchange - MATLAB Central (mathworks.com)) as I planned in first place.

I still needed to make sure they lie inside the outer boundary of your shape

now we can see that some normals are discontinuous as a portion of them goes outside the boundary and back inside.

that is reflected in the data (see below) , as for exemple here for the 5th normal, you can clearly sees the jump from the first to the second segment.

with some further code refinement we can decide wheter we keep all segments in the stored data (in cell array : data_out{k} ) of if you consider this or that segment to be removed from the data (up to you to decide)

same discontinuity can also be vizualized in the second plot

I annotated the curves (N1, N2,...) so you have a better understanding of which normal we are talking about

you can define how many normals you want to plot in this line :

N = 6; % specify here number of normals you want to plot

hope it helps !

Code :

%% load data files 
x = readmatrix('Space_a.txt');
y = readmatrix('Space_b.txt');
z = readmatrix('mz_Data.txt');

% % boundary points 
[row,col]=find(abs(z)>eps);% 1/ Find the locations of all non zero entries :
k=boundary(row,col,1);
xb = x(col(k));
yb = y(row(k));

% get "white line" points coordinates 
ztarget = 0.05;
tol = 0.02;
[r,c] =  find(abs(z - ztarget)<tol);
xl = x(c);
yl = y(r);
zl = z(c,r);

% sort  y data
[yl,ia] = sort(yl);
xl = xl(ia);
zl = zl(ia,:);

%% smoothed curve through the points 
% see Fex : https://www.mathworks.com/matlabcentral/fileexchange/25634-smoothn/
zs = smoothn({xl,yl});
xs =    zs{1};
ys =    zs{2};

%% plot
figure(1)
imagesc(x,y,z)
zmin = min(z(:));
zmax = max(z(:));
colormap(redwhiteblue(zmin, zmax));
colorbar;
set(gca,'YDir','normal');
hold on 

plot(xs,ys,'k','linewidth',2) ; % plot the line 
% plot(xb,yb,'m--','linewidth',2); % plot the boudary

%% find the normal vectors to the curve
% see Fex : https://fr.mathworks.com/matlabcentral/fileexchange/32696-2d-line-curvature-and-normals?s_tid=ta_fx_results
N = 6; % specify here number of normals you want to plot 

ns = numel(xs);
r = round(ns/(N-1)); % downsampling factor
ind = (1:r:ns);

% check how far ind(end) is from ns; if too far then simply add ns to the list (probably all the time)
if abs(ind(end) - ns)>round(r/10)
    ind = [ind ns]; % add ns to the vector
end

Vertices = [xs(ind) ys(ind)];
NL=LineNormals2D(Vertices); 
% NL is a unity vector giving the x,y projections of the direction vector (slope) of the normal lines

plot(xs(ind),ys(ind),'gX','markersize',10); % plot the green crosses 

%% constrain the normals to be inside the shape boundary & store xin,yin,zin data in cell array

for k =1:numel(ind)
    yv = Vertices(k,2) + (x-Vertices(k,1))*(NL(k,2)/NL(k,1));
    % constrain the normals to be inside the shape boundary
    [in, on] = inpolygon(x,yv,xb,yb);
    xin = x(in);
    yin = yv(in);
    % find corresponding z data
    zin = interp2(x,y,z,xin,yin); % notice mzs is very close to target z value (ztarget)
    plot(xin,yin,'y*','markersize',1);
    legstr{k} = [ 'N' num2str(k)];
    text(Vertices(k,1)+5,Vertices(k,2),legstr{k});
    
    % store xin,yin,zin data in cell array
    data_out{k} = [xin(:) yin(:) zin(:)]; %#ok<SAGROW>
end

axis equal  % to make sure perpendicular lines are actually showed perpendicular !!
grid off
hold off

%% finally, plot the profile zin 
figure(2);
hold on 

for k =1:numel(ind)
    plot(data_out{k}(:,1),data_out{k}(:,3),'*'); % plot zin vs xin
end
legend(legstr);
xlabel('xin');
hold off

[Mathieu NOE]

after zooming in the main plot I have the feeling that the normals are not 100% perpendicular to the black line and that is probably related to the function LineNormals2D.m I am using from Fex : 2D Line Curvature and Normals - File Exchange - MATLAB Central (mathworks.com))

there's one point I haven't yet completely understood is how to use the second argument (Lines ) and that is probably why I don't get the expected results.

so , assuming I need a bit of time to figure out how to use LineNormals2D.m correctly , I have devised another simpler solution based on picking the points before and after the "green" points and simply computing the normal vector of that linear segment (local tangent)

I believe the results below are more in line whith my / our expectations ....

%% load data files 
x = readmatrix('Space_a.txt');
y = readmatrix('Space_b.txt');
z = readmatrix('mz_Data.txt');

% % boundary points 
[row,col]=find(abs(z)>eps);% 1/ Find the locations of all non zero entries :
k=boundary(row,col,1);
xb = x(col(k));
yb = y(row(k));

% get "white line" points coordinates 
ztarget = 0.05;
tol = 0.02;
[r,c] =  find(abs(z - ztarget)<tol);
xl = x(c);
yl = y(r);
zl = z(c,r);

% sort  y data
[yl,ia] = sort(yl);
xl = xl(ia);
zl = zl(ia,:);

%% smoothed curve through the points 
% see Fex : https://www.mathworks.com/matlabcentral/fileexchange/25634-smoothn/
zs = smoothn({xl,yl});
xs =    zs{1};
ys =    zs{2};

%% plot
figure(1)
imagesc(x,y,z)
zmin = min(z(:));
zmax = max(z(:));
colormap(redwhiteblue(zmin, zmax));
colorbar;
set(gca,'YDir','normal');
hold on 

plot(xs,ys,'k','linewidth',2) ; % plot the line 
% plot(xb,yb,'m--','linewidth',2); % plot the boudary

%% find the normal vectors to the curve
N = 6; % specify here number of normals you want to plot 

ns = numel(xs);
r = round(ns/(N-1)); % downsampling factor
ind = (2:r:ns-1); % you have to start at 2 and end at ns -1 to have one point before and after for the normal vector computation

% check how far ind(end) is from ns; if too far then simply add ns to the list (probably all the time)
if abs(ind(end) - ns)>round(r/10)
    ind = [ind ns-1]; % add ns-1 to the vector
end

% compute normal vector 
for k =1:numel(ind)
    A = [xs(ind(k)-1), ys(ind(k)-1)];
    B = [xs(ind(k)+1), ys(ind(k)+1)];
    % Calculate the vector AB
    AB = B - A;
    % Determine the orientation (choose one of the orthogonal vectors)
    normal_vec = [AB(2), -AB(1)]; % This is one possible normal vector
    % Optionally, normalize the normal vector (adjust its length)
    NV(k,:) = normal_vec / norm(normal_vec); % Make it a unit vector
end

plot(xs(ind),ys(ind),'gX','markersize',10); % plot the green crosses 

%% constrain the normals to be inside the shape boundary & store xin,yin,zin data in cell array

for k =1:numel(ind)
    yv = ys(ind(k)) + (x-xs(ind(k)))*(NV(k,2)/NV(k,1));
    % constrain the normals to be inside the shape boundary
    [in, on] = inpolygon(x,yv,xb,yb);
    xin = x(in);
    yin = yv(in);
    % find corresponding z data
    zin = interp2(x,y,z,xin,yin); % notice mzs is very close to target z value (ztarget)
    plot(xin,yin,'y*','markersize',1);
    legstr{k} = [ 'N' num2str(k)];
    text(xs(ind(k))+5,ys(ind(k)),legstr{k});
    
    % store xin,yin,zin data in cell array
    data_out{k} = [xin(:) yin(:) zin(:)]; %#ok<SAGROW>
end

axis equal  % to make sure perpendicular lines are actually showed perpendicular !!
grid off
hold off

%% finally, plot the profile zin 
figure(2);
hold on 

for k =1:numel(ind)
    plot(data_out{k}(:,1),data_out{k}(:,3),'*'); % plot zin vs xin
end
legend(legstr);
xlabel('xin');
hold off

[Mathieu NOE]

hello again Richard

we can also use the gradient function to compute the slope of the line

once we know the slope of the line we have the slope of the normal vector :

slope_normals = -1/slope_line

%% load data files 
x = readmatrix('Space_a.txt');
y = readmatrix('Space_b.txt');
z = readmatrix('mz_Data.txt');

% % boundary points 
[row,col]=find(abs(z)>eps);% 1/ Find the locations of all non zero entries :
k=boundary(row,col,1);
xb = x(col(k));
yb = y(row(k));

% get "white line" points coordinates 
ztarget = 0.05;
tol = 0.02;
[r,c] =  find(abs(z - ztarget)<tol);
xl = x(c);
yl = y(r);
zl = z(c,r);

% sort  y data
[yl,ia] = sort(yl);
xl = xl(ia);
zl = zl(ia,:);

%% smoothed curve through the points 
% see Fex : https://www.mathworks.com/matlabcentral/fileexchange/25634-smoothn/
zs = smoothn({xl,yl});
xs =    zs{1};
ys =    zs{2};

%% plot
figure(1)
imagesc(x,y,z)
zmin = min(z(:));
zmax = max(z(:));
colormap(redwhiteblue(zmin, zmax));
colorbar;
set(gca,'YDir','normal');
hold on 

plot(xs,ys,'k','linewidth',2) ; % plot the line 
% plot(xb,yb,'m--','linewidth',2); % plot the boudary

%% find the normal vectors to the curve
% see Fex : https://fr.mathworks.com/matlabcentral/fileexchange/32696-2d-line-curvature-and-normals?s_tid=ta_fx_results
N = 6; % specify here number of normals you want to plot 

ns = numel(xs);
r = round(ns/(N-1)); % downsampling factor
ind = (1:r:ns);

% check how far ind(end) is from ns; if too far then simply add ns to the list (probably all the time)
if abs(ind(end) - ns)>round(r/10)
    ind = [ind ns]; % add ns to the vector
end

plot(xs(ind),ys(ind),'gX','markersize',10); % plot the green crosses 

%% constrain the normals to be inside the shape boundary & store xin,yin,zin data in cell array
slope = gradient(ys)./gradient(xs); % "white line" slope computation
for k =1:numel(ind)
    yv = ys(ind(k)) + (x-xs(ind(k)))*(-1/slope(ind(k))); % normal vectors
    % constrain the normals to be inside the shape boundary
    [in, on] = inpolygon(x,yv,xb,yb);
    xin = x(in);
    yin = yv(in);
    % find corresponding z data
    zin = interp2(x,y,z,xin,yin); % notice mzs is very close to target z value (ztarget)
    plot(xin,yin,'y*','markersize',1);
    legstr{k} = [ 'N' num2str(k)];
    text(xs(ind(k))+5,ys(ind(k)),legstr{k});
    
    % store xin,yin,zin data in cell array
    data_out{k} = [xin(:) yin(:) zin(:)]; %#ok<SAGROW>
end

axis equal  % to make sure perpendicular lines are actually showed perpendicular !!
grid off
hold off

%% finally, plot the profile zin 
figure(2);
hold on 

for k =1:numel(ind)
    plot(data_out{k}(:,1),data_out{k}(:,3),'*'); % plot zin vs xin
end
legend(legstr);
xlabel('xin');
hold of

[Roderick]

Hi @Mathieu NOE. Thank you very much for your very detailed answer to my question! It has been very helpful, I really appreciate it.

[Mathieu NOE]

As always, my pleasure

Answer 2

[Aditya Saikumar]

Hi Richard,

From the graph given, the left region has "mz" values less than 0 and the right region has "mz" values greater than 0. We can extract the data points corresponding to the left and the right region along with their "mz" values using "find" function. The following code achieves the same.

% Load data 
a = load('Space_a.txt'); 
b = load('Space_b.txt'); 
mz = load('mz_Data.txt'); 

% Extract the indices corresponding to the right region 
[aIndicesRight, bIndicesRight] = find(mz > 0); 

% Get the noraml values of the right region: 
fprintf("Right:\n"); 
for idx = 1:10 
    i = aIndicesRight(idx); 
    j = bIndicesRight(idx); 
    fprintf("a = %f, b = %f, mz = %f\n", a(i), b(j), mz(i, j)); 
end 
fprintf("\n"); 

Here is the output of the above code:

Please find more information about the “find” function in the following MATLAB Documentation link:

https://www.mathworks.com/help/releases/R2024a/matlab/ref/find.html

I hope this helps!

[Roderick]

Hi @Aditya Saikumar. Thank you very much! I have modified the question to include some additional information about the normal to the tangent at each (x,y) point of the obtained curve. Hope it helps.