Sorting two column vectors into 3D matrix based on position

Question

Using the imfindcircles function in MATLAB to track circles in two images. I start with approximately a grid of circles which deforms. I am trying to sort the two column vector from imfindcircles into matrices so that neighbouring circles are neighbouring elements in the matrices. The first image the circles conform to a grid and the following code works:

[centXsort,IX] = sortrows(centres1,1); %sort by x 
centYsort =zeros(289,2); %preallocate
for i = 1:17:289
    [sortedY,IY] = sortrows(centXsort(i:i+16,:),2); %sort by y within individual column
    centYsort(i:i+16,:) = sortedY; 
end
cent1mat = reshape(centYsort,17,17,2); %reshape into centre matrices

This doesn't work for the second image as some of the circles overlap in the x or y direction, but neighbouring circles never overlap. This means that in the second set of matrices the neighbouring circles aren't neighbouring elements after sorting.

Is there a way to approximate a scatter of points into a matrix?

Answer 1

This answer doesn't work in every single case, but it seems good enough for situations where the points don't vary too wildly.

---

My idea is to start at the grid corners and work our way along the outside diagonals of the matrix, trying to "grab" the nearest points that seem like they fit into the grid-points based any surrounding points we've already captured.

---

You will need to provide:

1. The number of rows (rows) and columns (cols) in the grid.
2. Your data points P arranged in a N x 2 array, rescaled to the unit square on [0,1] x [0,1]. (I assume the you can do this through visual inspection of the corner points of your original data.)
3. A weight parameter edge_weight to tell the algorithm how much the border points should be attracted to the grid border. Some tests show that 3-5 or so are good values.

---

The code, with a test case included:

%// input parameters
rows = 11;     
cols = 11;     
edge_weight = 4;

%// function for getting squared errors between the points list P and a specific point pt
getErr =@(P,pt) sqrt(  sum( bsxfun(@minus,P,pt(:)').^2, 2 )  );    %'

output_grid = zeros(rows,cols,2);   %// output grid matrix
check_grid = zeros(rows,cols);      %// matrix flagging the gridpoints we have covered
[ROW,COL] = meshgrid(...            %// coordinate points of an "ideal grid"
    linspace(0,1,rows),...
    linspace(0,1,cols));


%// create a test case
G = [ROW(:),COL(:)];                       %// the actual grid-points
noise_factor = 0.35;                       %// noise radius allowed
rn = noise_factor/rows;
cn = noise_factor/cols;
row_noise = -rn + 2*rn*rand(numel(ROW),1);
col_noise = -cn + 2*cn*rand(numel(ROW),1);
P = G + [row_noise,col_noise];             %// add noise to get points


%// MAIN LOOP
d = 0;                
while ~isempty(P)                       %// while points remain...
    d = d+1;                            %// increase diagonal depth (d=1 are the outer corners)
    for ii = max(d-rows+1,1):min(d,rows)%// for every row number i...
        i = ii;
        j = d-i+1;                      %// on the dth diagonal, we have d=i+j-1          
        for c = 1:4                     %// repeat for all 4 corners
            if i<rows & j<cols & ~check_grid(i,j) %// check for out-of-bounds/repetitions

                check_grid(i,j) = true;         %// flag gridpoint
                current_gridpoint = [ROW(i,j),COL(i,j)];

                %// get error between all remaining points and the next gridpoint's neighbours
                if i>1
                    errI = getErr(P,output_grid(i-1,j,:));
                else
                    errI = edge_weight*getErr(P,current_gridpoint);
                end
                if check_grid(i+1,j)
                    errI = errI + edge_weight*getErr(P,current_gridpoint);
                end
                if j>1
                    errJ = getErr(P,output_grid(i,j-1,:));
                else
                    errJ = edge_weight*getErr(P,current_gridpoint);
                end
                if check_grid(i,j+1)
                    errJ = errJ + edge_weight*getErr(P,current_gridpoint);
                end

                err = errI.^2 + errJ.^2;


                %// find the point with minimal error, add it to the grid,
                %//     and delete it from the points list
                [~,idx] = min(err);                         
                output_grid(i,j,:) = permute( P(idx,:), [1 3 2] );
                P(idx,:) = [];

            end

            %// rotate the grid 90 degrees and repeat for next corner
            output_grid = cat(3, rot90(output_grid(:,:,1)), rot90(output_grid(:,:,2)));
            check_grid = rot90(check_grid);
            ROW = rot90(ROW);               
            COL = rot90(COL);

        end
    end
end

Code for plotting the resulting points with edges:

%// plotting code
figure(1); clf; hold on;
axis([-0.1 1.1 -0.1 1.1])
for i = 1:size(output_grid,1)
    for j =  1:size(output_grid,2)
        scatter(output_grid(i,j,1),output_grid(i,j,2),'b')
        if i < size(output_grid,1)
            plot(   [output_grid(i,j,1),output_grid(i+1,j,1)],...
                [output_grid(i,j,2),output_grid(i+1,j,2)],...
                'r');
        end
        if j < size(output_grid,2)
            plot(   [output_grid(i,j,1),output_grid(i,j+1,1)],...
                [output_grid(i,j,2),output_grid(i,j+1,2)],...
                'r');
        end
    end
end

Answer 2

I've developed a solution, which works for my case but might not be as robust as required for some. It requires a known number of dots in a 'square' grid and a rough idea of the spacing between the dots. I find the 'AlphaShape' of the dots and all the points that lie along the edge. The edge vector is shifted to start at the minimum and then wrapped around a matrix with the corresponding points are discarded from the list of vertices. Probably not the best idea for large point clouds but good enough for me.

R = 50; % search radius
xy = centres2;
x = centres2(:,1);
y = centres2(:,2);

for i = 1:8
    T = delaunay(xy); % delaunay
    [~,r] = circumcenter(triangulation(T,x,y)); % circumcenters
    T = T(r < R,:); % points within radius
    B = freeBoundary(triangulation(T,x,y)); % find edge vertices
    A = B(:,1);

    EdgeList = [x(A) y(A) x(A)+y(A)]; % find point closest to origin and rotate vector
    [~,I] = min(EdgeList);
    EdgeList = circshift(EdgeList,-I(3)+1);

    n = sqrt(length(xy)); % define zeros matrix
    matX = zeros(n); % wrap x vector around zeros matrix
    matX(1,1:n) = EdgeList(1:n,1);
    matX(2:n-1,n) = EdgeList(n+1:(2*n)-2,1);
    matX(n,n:-1:1) = EdgeList((2*n)-1:(3*n)-2,1);
    matX(n-1:-1:2,1) = EdgeList((3*n)-1:(4*n)-4,1);

    matY = zeros(n); % wrap y vector around zeros matrix
    matY(1,1:n) = EdgeList(1:n,2);
    matY(2:n-1,n) = EdgeList(n+1:(2*n)-2,2);
    matY(n,n:-1:1) = EdgeList((2*n)-1:(3*n)-2,2);
    matY(n-1:-1:2,1) = EdgeList((3*n)-1:(4*n)-4,2);

    centreMatX(i:n+i-1,i:n+i-1) = matX; % paste into main matrix
    centreMatY(i:n+i-1,i:n+i-1) = matY;

    xy(B(:,1),:) = 0; % discard values
    xy = xy(all(xy,2),:);   
    x = xy(:,1);
    y = xy(:,2);   

end

centreMatX(centreMatX == 0) = x;
centreMatY(centreMatY == 0) = y;