Search code examples
matlabimage-processingdiscrete-mathematicspoint-clouds

Finding largest empty rectangle(disorientated) in 2D Point Cloud


As described, I would like to find the largest empty rectangle in a given 2D point cloud. This information will be then used to get parallel lines shown as green lines in the images, which will then be used to get the angle the LiDAR has rotated.

I've tried convex hull, boundary, and alpha shape functions in matlab with this point cloud data. From the looks of it, I think i have to

(1) break the point cloud into two chunks (gift wrap algo?) and then

(2) apply boundary function; With coordinates extracts from boundary function,

(3) I'm thinking to run RANSAC line fitting.

However, RANSAC needs more points to judge the "fitness" of the line. Currently I'm exploring Hough Transform to see if that line detection would work for this case.

Thus the question here is,

(1) am I in the right track to find the green lines? Or is there a better way?

(2) How to get the angular value/gradient from the line when the line is vertical (tan(90) = inf)

The second image depicts the origin may be in the center of the cloud points or at an offset. One can assume there will be parallel lines with varying distance from the center

Parallel lines to be found without Offset

Parallel lines to be found with Offset


Solution

  • You can use kmeans clustering to find two clusters then find two points each from a cluster that distance between them is minimum

    %a function for rotation of points
    function out = rot(theta,data)
        mat = [cos(theta), -sin(theta);
               sin(theta), cos(theta)];
        out = data * mat.';
    end
    a= [rand(505,1)*.4 , rand(505,1) ];
    b= [rand(500,1)*.4 + .6 , rand(500,1)];
    a = [a;b];
    rt = rot(.4,a);
    %kmeans clustering
    km=kmeans(rt,2);
    class1 = rt(km==1,:);
    class2=rt(km==2,:);
    %find min distance
    [srch, d] = dsearchn(class1, class2);
    [mn, I] = min(d);
    
    plot(rt(:,1), rt(:,2),'.')
    hold on
    p1 = class1(srch(I),:);
    p2 = class2(I,:);
    plot(p1(1),p1(2),'*r')
    plot(p2(1),p2(2),'*g')
    xy = [class1(srch(I),:);class2(I,:)]
    plot(xy(:,1),xy(:,2),'-g')
    axis equal