Hello, I want to limit the number of coordinates generated by the method bwboundaries or boundary while having the same step size.

[b1,l1] = bwboundaries(Image1_seuil,'noholes');

[b2,l2] = bwboundaries(Image2_seuil,'noholes');

hold on

% Coordinates

% I want to have the same array size for x1, x2, y1 and y2

x1=[];

y1=[];

x2=[];

y2=[];

for k1 = 1:length(b1)

boundary1 = b1{k1};

plot(boundary1(:,2), boundary1(:,1), 'r', 'LineWidth', 1)

x1 = [x1; boundary1(:,2)];

y1 = [y1; boundary1(:,1)];

end

for k2 = 1:length(b2)

boundary2 = b2{k2};

plot(boundary2(:,2), boundary2(:,1), 'g', 'LineWidth', 1)

x2 = [x2; boundary2(:,2)];

y2 = [y2; boundary2(:,1)];

end

[Mathieu NOE]

hello

it would help if you could share b1 and b2 as mat file

how do you want to reduce the amount of data if you keep step size ??

[Joumana]

I saved b1 and b2 in example.mat File.

Answer 1

[Mathieu NOE]

hello again

I could reduce the amount of data by a factor 3 roughly ,

more could be done if the large "straigth" lines did not include some local kincks

Code suggestion :

load('example.mat');

x1=[];
y1=[];
x2=[];
y2=[];

c_threshold = 0.5; % curvature threshold parameter

figure(1)
for k1 = 1:length(b1)
   
   boundary1 = b1{k1};
   x = boundary1(:,2);
   y = boundary1(:,1);

    [k] = find_high_curvature(x,y,c_threshold);
    x_out = x(k);
    y_out = y(k);
   plot(x, y, '*r', x_out, y_out, 'ob' ,'LineWidth', 2)     
   hold on
   x1 = [x1; x_out];
   y1 = [y1; y_out];
end


figure(2)
for k2 = 1:length(b2)
   
   boundary2 = b2{k2};
   x = boundary2(:,2);
   y = boundary2(:,1);

    [k] = find_high_curvature(x,y,c_threshold);
    x_out = x(k);
    y_out = y(k);
   plot(x, y, '*r', x_out, y_out, 'ob' ,'LineWidth', 2)     
   hold on
   x2 = [x2; x_out];
   y2 = [y2; y_out];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [highCurvatureIndexes] = find_high_curvature(x,y,threshold)

% Now run along the (x, y)  curve 
% and find the radius of curvature at each location.
numberOfPoints = length(x);
curvature = zeros(1, numberOfPoints);
for t = 1 : numberOfPoints
  if t == 1
    index1 = numberOfPoints;
    index2 = t;
    index3 = t + 1;
  elseif t >= numberOfPoints
    index1 = t-1;
    index2 = t;
    index3 = 1;
  else
    index1 = t-1;
    index2 = t;
    index3 = t + 1;
  end
  % Get the 3 points.
  x1 = x(index1);
  y1 = y(index1);
  x2 = x(index2);
  y2 = y(index2);
  x3 = x(index3);
  y3 = y(index3);
  % Now call Roger's formula:
  % http://www.mathworks.com/matlabcentral/answers/57194#answer_69185
  curvature(t) = 2*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1)) ./ ...
  sqrt(((x2-x1).^2+(y2-y1).^2)*((x3-x1).^2+(y3-y1).^2)*((x3-x2).^2+(y3-y2).^2));
end

% Find high curvature points - 
% indexes where the curvature is greater than threshold
highCurvatureIndexes = find(curvature > threshold);

end

[Joumana]

Thank you for your help, but the length of (x1,y1) is still different from (x2,y2). I need them to have the same size so I can calculate some distances between both images.

[Mathieu NOE]

ok , I didn't understand in first place because you mentionned "step" size and not array size

would it work if we simply pad the shorter array with NaN's ?

how do you compute the distance ?

[Joumana]

Yeah I am sorry, I mixed things up. I use x1,x2,y1 and y2 to calculate distancex, distancey... And I need them to have the same size, so I can have a correct calculation of the distance.

  % Distance x
 for lenx1=1:length(x1)
     for lenx2=1:length(x2)
            distancex=[distancex x2(lenx2)-x1(lenx1)];
     end
 end 

 % Distance y
 for leny1=1:length(y1)
     for leny2=1:length(y2)
            distancey=[distancey y2(leny2)-y1(leny1)];
     end
 end

 for i=1:length(distancex)
     if length(distancex)~= length(distancey)
         break
     end
     distance=[distance ( distancex(i)^2 + distancey(i)^2 )^0.5 ];

 end



%-----------------------Comparaison des valeurs des distances------------

for a=1:length(x1):length(distance)
     Les_distances_min=[Les_distances_min min(distance(a:a+length(x1)-1))];
     
end

[Mathieu NOE]

hello again

the trick I used here is simply to interpolate the longuest data (x2,y2) to the size of the smaller data set (x1,y1)

you can see that the two plots are identical , even though the second on is with (x2,y2) reduced

mission accomplished !

then I looked at the distance computation and I made some suggestion to make it run muuuch faster (less than 0.1 seconds vs many minutes). For loops , especially if they are long and nested (multiple ) , think to preallocate memory instead of let a array grow during iteration (argh !).

tic

load('example.mat');

x1=[];
y1=[];
x2=[];
y2=[];

c_threshold = 0.5; % curvature threshold parameter

for k1 = 1:length(b1)
   
   boundary1 = b1{k1};
   x = boundary1(:,2);
   y = boundary1(:,1);

    [k] = find_high_curvature(x,y,c_threshold);
    x_out = x(k);
    y_out = y(k);
   x1 = [x1; x_out];
   y1 = [y1; y_out];
end


for k2 = 1:length(b2)
   
   boundary2 = b2{k2};
   x = boundary2(:,2);
   y = boundary2(:,1);

    [k] = find_high_curvature(x,y,c_threshold);
    x_out = x(k);
    y_out = y(k);
   x2 = [x2; x_out];
   y2 = [y2; y_out];
end

figure(1)
subplot(1,2,1),plot(x1, y1, '*r', x2, y2, '*b','LineWidth', 1)
title('before resizing x2,y2')
n1 = numel(x1);
n2 = numel(x2);

% reduction of x2,y2 size to match x1,y1 size
x2new = interp1((1:n2),x2,linspace(1,n2,n1))';
y2new = interp1((1:n2),y2,linspace(1,n2,n1))';

subplot(1,2,2),plot(x1, y1, '*r', x2new, y2new, '*c','LineWidth', 1)
title('after resizing x2,y2')

% replace x2,y2 old with new data
x2 = x2new;
y2 = y2new;
n1 = numel(x1);
n2 = numel(x2);

%% Distance computations
% Distance x

% % original code (slow, only distancex alone needs +30 seconds to complete)
% distancex=[];
%  for lenx1=1:length(x1)
%      for lenx2=1:length(x2)
%             distancex=[distancex x2(lenx2)-x1(lenx1)];
%      end
%  end 
 

% prealocate data for faster execution (Elapsed time is 0.020392 seconds!!)
distancex=zeros(1,n1*n2);
k = 0;
 for lenx1=1:n1
     for lenx2=1:n2
            k = k+1;
            distancex(k)= x2(lenx2)-x1(lenx1);
     end
 end 
 
%  % another alternative code (suggestion up to you to try)
%  % Elapsed time is 0.014494 seconds.
%  distX = x2'-x1;
%  distX = reshape(distX',1,[]); % convert 2D array to col vector
%  diff = sum(abs(distancex-distX)) % diff =  0 !! distX and distancex are equal
 
 
%  % Distance y
%  for leny1=1:length(y1)
%      for leny2=1:length(y2)
%             distancey=[distancey y2(leny2)-y1(leny1)];
%      end
%  end


% prealocate data for faster execution (Elapsed time is 0.020392 seconds!!)
distancey=zeros(1,n1*n2);
k = 0;
 for leny1=1:n1
     for leny2=1:n2
            k = k+1;
            distancey(k)= y2(leny2)-y1(leny1);
     end
 end 
 
% original code (slow)
%  for i=1:n1*n2
%      distance=[distance ( distancex(i)^2 + distancey(i)^2 )^0.5 ];
%  end

% much faster and in one line : 
distance = sqrt(distancex.^2 + distancey.^2);

%-----------------------Comparaison des valeurs des distances------------
% original code (slow)
% for a=1:length(x1):length(distance)
%      Les_distances_min=[Les_distances_min min(distance(a:a+length(x1)-1))];
%      
% end

% much faster
ind = 1:n1:n1*n2;
Les_distances_min=zeros(1,numel(ind));% prealocate data for faster execution

for k=1:numel(ind)
    a = ind(k);
     Les_distances_min(k)= min(distance(a:a+n1-1));
     
end

 toc


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [highCurvatureIndexes] = find_high_curvature(x,y,threshold)

% Now run along the (x, y)  curve 
% and find the radius of curvature at each location.
numberOfPoints = length(x);
curvature = zeros(1, numberOfPoints);
for t = 1 : numberOfPoints
  if t == 1
    index1 = numberOfPoints;
    index2 = t;
    index3 = t + 1;
  elseif t >= numberOfPoints
    index1 = t-1;
    index2 = t;
    index3 = 1;
  else
    index1 = t-1;
    index2 = t;
    index3 = t + 1;
  end
  % Get the 3 points.
  x1 = x(index1);
  y1 = y(index1);
  x2 = x(index2);
  y2 = y(index2);
  x3 = x(index3);
  y3 = y(index3);
  % Now call Roger's formula:
  % http://www.mathworks.com/matlabcentral/answers/57194#answer_69185
  curvature(t) = 2*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1)) ./ ...
  sqrt(((x2-x1).^2+(y2-y1).^2)*((x3-x1).^2+(y3-y1).^2)*((x3-x2).^2+(y3-y2).^2));
end

% Find high curvature points - 
% indexes where the curvature is greater than threshold
highCurvatureIndexes = find(curvature > threshold);

end

[Joumana]

Thank you so much, this is perfect! I didn't know how to preallocate memory in matlab, I use it in C but never in matlab. Thank you for your help and your advice. Have a nice day!

[Mathieu NOE]

Thanks - have a great day too

[Joumana]

Hello again, I have a question: if I want to increase the number of coordinates, should I change the value of the threshold? I would like to compare my calculations with different number of coordinates.

[Mathieu NOE]

yeap

you can play with the parameter

c_threshold = 0.5; % curvature threshold parameter

[Joumana]

I tried doing that, but when I go below 0.5 the number of coordinates don't change. And if I go pass 0.5 the number of coordinates is cut down significantly.

[Mathieu NOE]

hello again

I changed one line in the function

change :

highCurvatureIndexes = find(curvature > threshold);

into this

highCurvatureIndexes = find(curvature >= threshold);

now if you put the threshold to zero or a very very low number you get more points , which basically get you back to your original data

I confirm that you can't really go above 0.5 otherwise the number of points is too low

I found myself also that playing with the threshold is not very gradual solution as your data show either a zero curvature (straigth portions) or some (few) positive values and nothing in between

plot of curvature data :

now if you want something more gradual we have to improve the code ; is it what you need ?

otherwise you can simply skip the process of data removal and keep all your original data

[Joumana]

Thank you, I am looking to change the number of coordinates gradually. The reason behind it is to be able to calculate a surface with different sets of coordinates ( for example with 500 corrdinates, 300 coordinates, 200 coordinates ... ) and comparing the results.

[Mathieu NOE]

ok I understand but I am not sure how to force my "curvature" function to give only a specified number of points (100,200,300,...) but maybe there is a solution if a apply also a "decimation" function but I am not sure how to keep the "corners" as the decimation process will act blind

[Joumana]

It dosen't have to be a specified number. But the number of coordinates rises gradually to be able to compare the calculations. I am sorry if I am not giving a clear explanation. Please tell me if you prefer me to explain in another way. And thank you for your help.

[Mathieu NOE]

tried something different , like downsampling your data close enough to the target points

it's simply picking very n points and defining the decimation factor so that the total qty is close as posible to the given N value

this is a different approach that has no more anything to do with my initial curvature computation

tic
clc
clearvars

load('example.mat');
x1=[];
y1=[];
x2=[];
y2=[];

N = 200; % approx total number of points

for k1 = 1:length(b1)
   boundary1 = b1{k1};
   x1 = [x1; boundary1(:,2)];
   y1 = [y1; boundary1(:,1)];
end

for k2 = 1:length(b2)
   boundary2 = b2{k2};
   x2 = [x2; boundary2(:,2)];
   y2 = [y2; boundary2(:,1)];
end

figure(1)
subplot(1,2,1),plot(x1, y1, '*r', x2, y2, '*b','LineWidth', 1)
title('before resizing')

nx = numel(x1);
if nx>N % pick 1 every decim points
    decim = ceil(nx/N);
    ind = 1:decim:nx;
    x1 = x1(ind);
    y1 = y1(ind);
end

nx = numel(x2);
if nx>N % pick 1 every decim points
    decim = ceil(nx/N);
    ind = 1:decim:nx;
    x2 = x2(ind);
    y2 = y2(ind);
end

subplot(1,2,2),plot(x1, y1, '*r', x2, y2, '*b','LineWidth', 1)
title('after resizing')

%% Distance computations
% Distance x
% % original code (slow, only distancex alone needs +30 seconds to complete)
% distancex=[];
%  for lenx1=1:length(x1)
%      for lenx2=1:length(x2)
%             distancex=[distancex x2(lenx2)-x1(lenx1)];
%      end
%  end 
 
% prealocate data for faster execution (Elapsed time is 0.020392 seconds!!)
distancex=zeros(1,n1*n2);
k = 0;
 for lenx1=1:n1
     for lenx2=1:n2
            k = k+1;
            distancex(k)= x2(lenx2)-x1(lenx1);
     end
 end 
 
%  % another alternative code (suggestion up to you to try)
%  % Elapsed time is 0.014494 seconds.
%  distX = x2'-x1;
%  distX = reshape(distX',1,[]); % convert 2D array to col vector
%  diff = sum(abs(distancex-distX)) % diff =  0 !! distX and distancex are equal
 
 
%  % Distance y
%  for leny1=1:length(y1)
%      for leny2=1:length(y2)
%             distancey=[distancey y2(leny2)-y1(leny1)];
%      end
%  end
% prealocate data for faster execution (Elapsed time is 0.020392 seconds!!)
distancey=zeros(1,n1*n2);
k = 0;
 for leny1=1:n1
     for leny2=1:n2
            k = k+1;
            distancey(k)= y2(leny2)-y1(leny1);
     end
 end 
 
% original code (slow)
%  for i=1:n1*n2
%      distance=[distance ( distancex(i)^2 + distancey(i)^2 )^0.5 ];
%  end
% much faster and in one line : 
distance = sqrt(distancex.^2 + distancey.^2);
%-----------------------Comparaison des valeurs des distances------------
% original code (slow)
% for a=1:length(x1):length(distance)
%      Les_distances_min=[Les_distances_min min(distance(a:a+length(x1)-1))];
%      
% end
% much faster
ind = 1:n1:n1*n2;
Les_distances_min=zeros(1,numel(ind));% prealocate data for faster execution
for k=1:numel(ind)
    a = ind(k);
     Les_distances_min(k)= min(distance(a:a+n1-1));
     
end
 toc