Different answers for the same algorithm

I am getting different values of the width upon calling the findwidth function or writing its algorithm (exactly the same) in the main script. Is there a reason to justify that and know what is true?

The main script looks like this:

for i=1:10
    :
for j=1:13
for p=1:length(c1)
     : 
  
  P1new=[P11(concavev(i):n,:);P11(1:c1(p),:)];
    width(i,j,p)=findwidth(P1new)
end
end
end

%the width function :
function width=findwidth(P11)

[k,av] = convhull(P11);
CP=[P11(k,1),P11(k,2)];
for x=1:(length(CP)-1)
   z=1:(length(CP)-1);
   for  z=1:(length(CP)-1); 
  D(x,z)=Seg_point_dis(CP(z,:),CP(z+1,:),CP(x,:))
   end
  end
  D1=max(D,[],2)    
width=min(D1);

end

[darova]

What is going on in this script? It's hard to understand

[farah fadel]

I am iterating over some vertices(i), edges(j), and number of intersection(p) to define a polygon, the issue is when I want to find its width (varies across each i,j,p)using the function findwidth, I will have a value of the width different than that if i just copied the findwidth function algorithm and aded it to the script

[darova]

Can't help. Sorry

[Aditya Patil]

Can you provide the entire script? Most likely, the issue is due to either use of random variables somewhere, or repeated use of variables.

[farah fadel]

%The purpose of this code is to decompose a polygon into two polygones
P11=[0.45 0.75;2.37 1.49;4.15 0.32;3.63 1.48;5.58 1.78;7.45 3.21;6.12 3.58;4.75 6.15;4.32 5.94;3.75 4.55;2.45 6.44;1.55 5.45;2.51 3.67;0.45 0.75]
n=length(P11)-1;%closed polygon
[ncc,concavev]=CheckConcavity(P11);%check the number of convex polygons(ncc) and their location (concavev)
%in the above given polygon ncc=5,n=13,concavev=[2;4;7;10;13]
for i=1:ncc %iterate over the number of concave vertices
    for j=1:n%iterate over all the edges
     
    e(j,:)=P11(j,:)-P11(j+1,:);%calcualting the edges of the polygon
    a(j)=e(j,2)/e(j,1);%slope of the edges
    b(j)=bin(P11(concavev(i),1),P11(concavev(i),2),a(j)) ;%intercept of the edges
    pp = [a(j),b(j)]; %1st line represent an edge
        % 2nd line represent the drawn from the concave vertex, and parallel to the edge
        for k=setdiff(1:n,[concavev(i),(concavev(i)-1)]) 
        x2 = [P11(k) P11(k+1)];
        y2 = [P11(k,2) P11(k+1,2)];
        %fit linear polynomial
        p2= polyfit(x2,y2,1);
        %calculate intersection of the twolines
        x_intersect(k) = fzero(@(x) polyval(pp-p2,x),1);
        y_intersect(k) = polyval(p2,x_intersect(k));
        %these intersections are the possible options to decompose the polygon at
        xi(j,k)=x_intersect(k);  yi(j,k)=y_intersect(k);
        end
       
        [n1,c1] =find(xi(j,:))%finding the nonzero terms, place where intersection took place  
        for p=1:length(c1)%number of intersection per edge
        P1(p,:)=[xi(j,c1(p)),yi(j,c1(p))];
          if c1(p)<concavev(i)
            P1new=[P11(concavev(i):n,:);P11(1:c1(p),:);P1(p,:);P11(concavev(i),:)];%decompsed polygon
          end
  %now we need to calculate the width of all possible polygons we could get to choose the one with min width 
[k1,av] = convhull(P1new);%calculate the convexhull of the pollygon
CP=[P1new(k1,1),P1new(k1,2)];%define the new convex polygon
for x=1:(size(CP,2)-1)%vertex
   for  z=1:(size(CP,2)-1) %edge
  D(x,z)=Seg_point_dis(CP(z,:),CP(z+1,:),CP(x,:));
  %finding all the possible distances between the vertices and the edges combination
   end
end
  D1=max(D,[],2) ;%The maximum dis btwn vertex and edge spans the polygon
width1=min(D1);% the width will be the maximum Span
          
            widthmethod2(i,j,p)=width1%from script
            widthmethod1(i,j,p)=findwidth(P1new)%from function
        end
    end
end

[farah fadel]

I tried to check for the varibles issues, but still nothing changed.(I changed the variable x in the function but still same result)

[Guillaume]

You haven't provided the code for Seg_point_dis. Possibly the problem comes from there.

In any case, there are many issues with the code of findwidth:

• No comments whatsoever. Coupled with meaningless variable names, it's impossible to know what it's meant to be doing. Even you in 6 months time won't be able to understand it.
• You're using length on a 2D matrix. length returns either the number of rows or the number of columns (whichever is greater) Most of the times it will work fine but if your convex hull is only one point your code will error. I know it's unlikely but why leave the chance when using size(CP, 2) is guaranteed to give you the number of rows all the time.
• You create a vector z=1:(length(CP)-1); and on the very next line throw it away and use z as a for loop iterator.
• You're building a matrix D in a double loop. At each step of the inner loop, you calculate the maximum of D while you're still building it.

[farah fadel]

Thank you for the feedback, I fixed the issues you told me about,but still the same problem.

This is the function Seg_point_dis if it would give any help

function dist=Seg_point_dis(a,b,x)
d_ab = norm(a-b);
d_ax = norm(a-x);
d_bx = norm(b-x);

if dot(a-b,x-b)*dot(b-a,x-a)>=0
    A = [a,1;b,1;x,1];
    dist = abs(det(A))/d_ab;        
else
    dist = min(d_ax, d_bx);
end

[farah fadel]

Based on the results, the answers obtained when writing the algorithm instead of calling the function seems more reasonable

[Aditya Patil]

Can you provide both the variants of the code so that we can compare the results and find out what the issue is. Also, it would be helpful if you comment, align and rename variables to make the code more readable.

[farah fadel]

I edited the scipt in the comment above, if there is anything still unclear please tell me.

Results:

The location of the nonempty var are the same but the values differ.

widthmethod2(1,:,1)=[3.4923 3.4923 3.4923 4.2720 4.0440 5.1376 5.1376 4.9211 3.4305 3.4305 4.5582 4.5582 4.8062]

widthmethod1(1;:,1)=[4.1941 ;5.1838 ; 4.4855 ; 1.4271 ;4.1062 ; 0.4544 ;4.5847 ;3.9177; 3.5832 ; 4.7346 ;4.3275 4.5915 ; 4.3693]

[Aditya Patil]

The arguments to bin function seem incorrect. Also, can you provide the modified code with changes that you made as per above comments by me and Guillaume?

[farah fadel]

the bin functions take the x,y and the slope respectively to output the intercept.

Yes sure, do I send them by email?

[Aditya Patil]

Can you attach the files here?

[farah fadel]

Done, thank you

[Guillaume]

I've only checked the first few lines of your code and stopped at CheckConc.

You're still using length on a 2-column matrix (both in the script and in CheckConc. I recommend that you never use length. For a vector, use numel. For a matrix, use size and explicitly specify which dimension you want to know the size of, rather than assuming that length will return the size of the dimension you are hoping for. So, use size(P1, 1) and size(v, 1)

The reason I stopped at CheckConc is that the function is wrong. It only performs N-1 checks for a N sided polygon so is missing one vertex in the concavity check. I've not tried to understand the math you're doing. There are more efficient algorithm to check that concavity. See stackoverflow and math.stackexhange. If you're guaranteed that your polygon is not self-intersecting that latter one would be best. It can be implemented without any loop in matlab.

Answer 1

[Aditya Patil]

Inside findwidth function, you have modified the variable D, which is a global variable. This creates the issue.

Instead, pass D as a parameter to the function, as follows ( parameter named E in the code below)

function width=findwidth(P11,ncc, E) 
    if ncc~=0 
        [k1,av] = convhull(P11); 
        CP=[P11(k1,1),P11(k1,2)]; 
    else 
        CP=P11; 
    end 
    for q=1:(length(CP)-1)%:(length(CP)-1) 
        for  z=1:(length(CP)-1)  
            E(q,z)=Seg_point_dis(CP(z,:),CP(z+1,:),CP(q,:)); 
        end 
    end 
    D1=max(E,[],2) ; 
    width=min(D1); 
end 

[farah fadel]

I am still getting different answer