Method to get as large and few as number boxes

oriLat=[48.3860120000000;49.3686790000000;48.0145280000000;46.6825390000000;47.4787370000000;46.5073030000000;46.5176190000000;44.9243550000000;45.2906840000000;42.3015220000000;41.9182910000000;44.6990460000000;45.7893010000000;45.2741950000000;43.7351650000000;43.8312240000000;41.7326470000000;41.3863160000000;44.9774490000000;45.2345130000000;47.3537330000000;44.8174190000000;43.3218860000000;41.5511100000000;39.4791860000000;37.8949860000000;37.1004480000000;39.4531150000000;38.0365030000000;38.6352170000000;35.2258250000000;31.1253830000000;26.7715300000000;24.8661310000000;29.9185130000000;30.6849560000000;30.0184100000000;28.9338120000000;29.4307800000000;27.9083070000000;25.8403780000000;29.7380790000000;28.9821380000000;31.7520090000000;31.3322390000000;34.5540170000000;40.2609740000000;46.2591090000000;48.3860120000000];
oriLon=[-124.725839000000;-94.9521110000000;-89.4892260000000;-91.9618890000000;-87.9292690000000;-87.3667670000000;-84.1179250000000;-87.8434330000000;-86.9777800000000;-87.8347690000000;-86.5978990000000;-86.2484740000000;-84.7727650000000;-83.3851040000000;-83.9477400000000;-82.6336410000000;-83.4538320000000;-82.4605990000000;-74.9927560000000;-70.8444300000000;-68.2697100000000;-66.9498950000000;-70.5538540000000;-69.9649820000000;-75.5930680000000;-75.3386230000000;-75.9796080000000;-76.0123120000000;-76.3220930000000;-77.2467040000000;-75.5336270000000;-81.4020960000000;-80.0321200000000;-80.6511890000000;-83.6792190000000;-88.0083960000000;-89.8450650000000;-89.4009660000000;-94.6703890000000;-97.0033250000000;-97.4226360000000;-101.400636000000;-103.281190000000;-106.417940000000;-111.074825000000;-120.622575000000;-124.363414000000;-123.547659000000;-124.725839000000];
N = 100 ;   % can be varied
lon = linspace(min(oriLon),max(oriLon),N) ;
lat = linspace(min(oriLat),max(oriLat),N) ;
%
[Lon,Lat] = meshgrid(lon,lat) ;
% pick points lying inside the given region 
idx = inpolygon(Lon(:), Lat(:),oriLon, oriLat) ;
Lon(~idx) = NaN ; Lat(~idx) = NaN ;
plot( oriLon, oriLat, 'b')
hold on
plot(Lon,Lat,'.r')

The above gives me multiple uniform grids fitting within a polygon. From above I would like to reduce number of recangulars by replacing them with larger grids/ boxes. Each box could have different size. Usually near the edge, the size would be smaller.

What would be the approach to go about this?

[John D'Errico]

Sorry. Breaking a general domain with such convoluted boundaries into rectangles of maximum size is quite difficult to do well. It is even difficult to do poorly. And sorry, but you will find no simple tool that solves it for you.

[Walter Roberson]

I programmed something along those lines 15-ish years ago. I do not recall the algorithm I used, though.

[Pete sherer]

any suggestions on lgic?

[Walter Roberson]

You asked to minimize the number of rectangles used. That makes it an optimization problem -- and literature implies that it is NP-hard .

You cannot just optimize locally for such a system. Imagine for example if there were a single "unoccupied" square (such as a lake) on your map near (-100,40) then that could break up what might have been a wide rectangle near that area, and suddenly the optimum might flip to instead using tall rectangles.

So if you need the best (fewest rectangles) possible for the situation, then you are going to have undertake some difficult optimization.

If you do not need the best then you might be able to use heuristics to get a "not unreasonable" fit. For example you might proceed and at each step look for the largest rectangle you could inscribe on the current map, mark it as occupied and iterate. If you find multiple rectangles with the same largest area and they do not overlap then you might be able to save steps by using them all; if the multiple largest rectangles do overlap then making decisions between them might take more work. This all will not result in an "optimal" solution -- "greedy" solutions often turn out to be non-optimal.

I have to wonder why you want to do this? The task you have set is not, for example, equivalent to determining the least number of "cuts" to cut out the shape from enclosing material.

[Steven Lord]

What's your ultimate goal in creating these boxes? What are you hoping to use these boxes to do?

Do you require these boxes to stay strictly within the boundaries of the shape or is some coverage of the area outside the shape acceptable?

Are the boxes required to be rectangles? Are general quadrilaterals allowed? Are polygons with other number of sides allowed (triangulate the area, tile it with hexagons, cover it with N-sided polygons for large N that approximate circles, etc.)? Are arbitrarily shaped regions allowed (jigsaw puzzle pieces or electoral district drawing / gerrymandering)?

What's your success criteria? What's most important from the list below?

• maximum amount of coverage / minimize uncovered area
• minimize amount of area outside the shape that is covered (if covering outside the area is allowed)
• minimize number of boxes
• maximize box size

[Walter Roberson]

For example, imagine that you have a 2 x 10 area in which the top right is not available. Would you want that covered by a (1 x 9 + 1 x 10) or as a (2 x 9 + 1 x 1) ? Or as a 2 x 10 with it being "forgiven" that it goes outside the boundary?

Answer 1

[Fabio Freschi]

The following attempt is far from being a perfect answer to your question. All perplexities of @Walter Roberson are also mine, so I have drastically simplified the problem assuming

1. the use of square boxes instead of rectangular ones
2. the possibility of randomly starting the creation of new boxes
3. the number of boxes is not minimum

If you are optimistic and benevolent, the solution found can be said to be "sub-optimal", but I think it's a starting point you can use to find a better solution.

clear variables, close all
% perimeter
oriLat = [48.3860120000000;49.3686790000000;48.0145280000000;46.6825390000000;47.4787370000000;46.5073030000000;46.5176190000000;44.9243550000000;45.2906840000000;42.3015220000000;41.9182910000000;44.6990460000000;45.7893010000000;45.2741950000000;43.7351650000000;43.8312240000000;41.7326470000000;41.3863160000000;44.9774490000000;45.2345130000000;47.3537330000000;44.8174190000000;43.3218860000000;41.5511100000000;39.4791860000000;37.8949860000000;37.1004480000000;39.4531150000000;38.0365030000000;38.6352170000000;35.2258250000000;31.1253830000000;26.7715300000000;24.8661310000000;29.9185130000000;30.6849560000000;30.0184100000000;28.9338120000000;29.4307800000000;27.9083070000000;25.8403780000000;29.7380790000000;28.9821380000000;31.7520090000000;31.3322390000000;34.5540170000000;40.2609740000000;46.2591090000000;48.3860120000000];
oriLon = [-124.725839000000;-94.9521110000000;-89.4892260000000;-91.9618890000000;-87.9292690000000;-87.3667670000000;-84.1179250000000;-87.8434330000000;-86.9777800000000;-87.8347690000000;-86.5978990000000;-86.2484740000000;-84.7727650000000;-83.3851040000000;-83.9477400000000;-82.6336410000000;-83.4538320000000;-82.4605990000000;-74.9927560000000;-70.8444300000000;-68.2697100000000;-66.9498950000000;-70.5538540000000;-69.9649820000000;-75.5930680000000;-75.3386230000000;-75.9796080000000;-76.0123120000000;-76.3220930000000;-77.2467040000000;-75.5336270000000;-81.4020960000000;-80.0321200000000;-80.6511890000000;-83.6792190000000;-88.0083960000000;-89.8450650000000;-89.4009660000000;-94.6703890000000;-97.0033250000000;-97.4226360000000;-101.400636000000;-103.281190000000;-106.417940000000;-111.074825000000;-120.622575000000;-124.363414000000;-123.547659000000;-124.725839000000];

% discretization size
N1 = 80;
N2 = 60;
x = linspace(min(oriLon),max(oriLon),N1) ;
y = linspace(min(oriLat),max(oriLat),N2) ;

% create meshgrid
[x,y] = meshgrid(x,y) ;

% pick points lying inside the given region
idx = inpolygon(x(:),y(:),oriLon, oriLat);
% reshape idx
idx = reshape(idx,size(x));

% coloring matrix (-1: out, 0: in & never touched, >0 in & touched)
col = zeros(size(idx));
col(~idx) = -NaN;

% plot
h = figure; hold on
plot(oriLon, oriLat, 'b')

% rng(2);

for icol = 1:N1*N2
    % select a random point to start
    idx0 = find(col(:) == 0);
    if isempty(idx0)
        break
    end
    [i0,j0] = ind2sub(size(col),idx0(randi(length(idx0))));

    % move north-est
    chk = true;
    k_ne = 0;
    while chk
        k_ne = k_ne+1;

        if any(any(col(i0:i0+k_ne,j0:j0+k_ne) ~= 0))
            k_ne = k_ne-1;
            % col(i0:i0+k_ne,j0:j0+k_ne) = icol;
            chk = false;
        end
    end

    % move south-west
    chk = true;
    k_sw = 0;
    while chk
        k_sw = k_sw+1;

        if any(any(col(i0-k_sw:i0+k_ne,j0-k_sw:j0+k_ne) ~= 0))
            k_sw = k_sw-1;
            col(i0-k_sw:i0+k_ne,j0-k_sw:j0+k_ne) = icol;
            chk = false;
        end
    end

    % plot colored points
    scatter(reshape(x(i0-k_sw:i0+k_ne,j0-k_sw:j0+k_ne),[],1),...
        reshape(y(i0-k_sw:i0+k_ne,j0-k_sw:j0+k_ne),[],1),...
        15,reshape(col(i0-k_sw:i0+k_ne,j0-k_sw:j0+k_ne),[],1),'filled');
    if k_ne+k_sw > 1
        xl = [x(i0-k_sw,j0-k_sw)
            x(i0+k_ne,j0-k_sw)
            x(i0+k_ne,j0+k_ne)
            x(i0-k_sw,j0+k_ne)
            x(i0-k_sw,j0-k_sw)];
        yl = [y(i0-k_sw,j0-k_sw)
            y(i0+k_ne,j0-k_sw)
            y(i0+k_ne,j0+k_ne)
            y(i0-k_sw,j0+k_ne)
            y(i0-k_sw,j0-k_sw)];
        plot(xl,yl,'-','LineWidth',2);
    end

end

fprintf('final number of boxes: %d\n',icol)

final number of boxes: 335

Hope it helps

[Fabio Freschi]

The first improvement you could try is to loop over the four edges of the square to check if they can be moved north-est-south-west, creating a rectangle with a larger area. A mor sistematic and optimal approach could be based on the search of the maximum area inscribed reectangle (link):

1. Find the largest rectangle inscribed in your picture
2. identify the sub-areas identified by the intersection of the original picture and the rectangle
3. Iterate

[Fabio Freschi]

I added the possibility of creation rectangles from squares, as suggested in my provious comment, with some cleanup and small improvements: the number of rectangles is now reduced to 25%-30% with respect to the previous method.

clear variables, close all
% perimeter
oriLat = [48.3860120000000;49.3686790000000;48.0145280000000;46.6825390000000;47.4787370000000;46.5073030000000;46.5176190000000;44.9243550000000;45.2906840000000;42.3015220000000;41.9182910000000;44.6990460000000;45.7893010000000;45.2741950000000;43.7351650000000;43.8312240000000;41.7326470000000;41.3863160000000;44.9774490000000;45.2345130000000;47.3537330000000;44.8174190000000;43.3218860000000;41.5511100000000;39.4791860000000;37.8949860000000;37.1004480000000;39.4531150000000;38.0365030000000;38.6352170000000;35.2258250000000;31.1253830000000;26.7715300000000;24.8661310000000;29.9185130000000;30.6849560000000;30.0184100000000;28.9338120000000;29.4307800000000;27.9083070000000;25.8403780000000;29.7380790000000;28.9821380000000;31.7520090000000;31.3322390000000;34.5540170000000;40.2609740000000;46.2591090000000;48.3860120000000];
oriLon = [-124.725839000000;-94.9521110000000;-89.4892260000000;-91.9618890000000;-87.9292690000000;-87.3667670000000;-84.1179250000000;-87.8434330000000;-86.9777800000000;-87.8347690000000;-86.5978990000000;-86.2484740000000;-84.7727650000000;-83.3851040000000;-83.9477400000000;-82.6336410000000;-83.4538320000000;-82.4605990000000;-74.9927560000000;-70.8444300000000;-68.2697100000000;-66.9498950000000;-70.5538540000000;-69.9649820000000;-75.5930680000000;-75.3386230000000;-75.9796080000000;-76.0123120000000;-76.3220930000000;-77.2467040000000;-75.5336270000000;-81.4020960000000;-80.0321200000000;-80.6511890000000;-83.6792190000000;-88.0083960000000;-89.8450650000000;-89.4009660000000;-94.6703890000000;-97.0033250000000;-97.4226360000000;-101.400636000000;-103.281190000000;-106.417940000000;-111.074825000000;-120.622575000000;-124.363414000000;-123.547659000000;-124.725839000000];

% discretization size
N1 = 80;
N2 = 60;
x = linspace(min(oriLon),max(oriLon),N1) ;
y = linspace(min(oriLat),max(oriLat),N2) ;

% create meshgrid
[x,y] = meshgrid(x,y) ;

% pick points lying inside the given region
idx = inpolygon(x(:),y(:),oriLon, oriLat);
% reshape idx
idx = reshape(idx,size(x));

% coloring matrix (-1: out, 0: in & never touched, >0 in & touched)
col = zeros(size(idx));
col(~idx) = -NaN;

% plot
h = figure; hold on
plot(oriLon, oriLat, 'b')

rng(2);

for icol = 1:N1*N2
    % select a random point to start
    idx0 = find(col(:) == 0);
    if isempty(idx0)
        break
    end
    [i0,j0] = ind2sub(size(col),idx0(randi(length(idx0))));

    % initial values
    k_n = 0;
    k_e = 0;
    k_s = 0;
    k_w = 0;

    % move north-est
    chk = true;
    while chk
        k_n = k_n+1;
        k_e = k_e+1;

        if any(col(i0:i0+k_e,j0:j0+k_n) ~= 0,'all')
            k_n = k_n-1;
            k_e = k_e-1;
            chk = false;
        end
    end

    % try only north
    chk = true;
    while chk
        k_n = k_n+1;

        if any(col(i0:i0+k_e,j0:j0+k_n) ~= 0,'all')
            k_n = k_n-1;
            chk = false;
        end
    end

    % try only east
    chk = true;
    while chk
        k_e = k_e+1;

        if any(col(i0:i0+k_e,j0:j0+k_n) ~= 0,'all')
            k_e = k_e-1;
            chk = false;
        end
    end

    % move south-west
    chk = true;
    while chk
        k_s = k_s+1;
        k_w = k_w+1;

        if any(col(i0-k_w:i0+k_e,j0-k_s:j0+k_n) ~= 0,'all')
            k_s = k_s-1;
            k_w = k_w-1;
            chk = false;
        end
    end

    % try only south
    chk = true;
    while chk
        k_s = k_s+1;

        if any(col(i0-k_w:i0+k_e,j0-k_s:j0+k_n) ~= 0,'all')
            k_s = k_s-1;
            chk = false;
        end
    end

    % try only west
    chk = true;
    while chk
        k_w = k_w+1;

        if any(col(i0-k_w:i0+k_e,j0-k_s:j0+k_n) ~= 0,'all')
            k_w = k_w-1;
            chk = false;
        end
    end

    % set new box id
    col(i0-k_w:i0+k_e,j0-k_s:j0+k_n) = icol;

    % plot colored points
    scatter(reshape(x(i0-k_w:i0+k_e,j0-k_s:j0+k_n),[],1),...
        reshape(y(i0-k_w:i0+k_e,j0-k_s:j0+k_n),[],1),...
        15,reshape(col(i0-k_w:i0+k_e,j0-k_s:j0+k_n),[],1),'filled');
    if k_n+k_s > 0 && k_e+k_w > 0
        xl = [x(i0-k_w,j0-k_s)
            x(i0+k_e,j0-k_s)
            x(i0+k_e,j0+k_n)
            x(i0-k_w,j0+k_n)
            x(i0-k_w,j0-k_s)];
        yl = [y(i0-k_w,j0-k_s)
            y(i0+k_e,j0-k_s)
            y(i0+k_e,j0+k_n)
            y(i0-k_w,j0+k_n)
            y(i0-k_w,j0-k_s)];

        plot(xl,yl,'-','LineWidth',2);
    end

end

fprintf('final number of boxes: %d\n',icol)

final number of boxes: 91

[Pete sherer]

Hi All,

Thans very much for the suggestions and thought processes. Yes the box could go over the coundary area. I will look at the code Fabio's provided and fine-tuning it from there.

Thanks again for all of the helpful tips.