How to omit/edit some elemental coordinates of a matrix

Hello

I have a rectanglular grid (28.66*10)m with the element dimension (1.433*1)m. On the other hand, there is a slope with vertices (0,10), (10,10), (18.66,5) and (28.66,5) that should be mapped onto the rectangular grid. The confusing part is how to omit the coordinates of elements which fall out of the slope area (outside of solid lines) and also preserve the small triangular (with three corners) and also trapzoidal elements (with 5 corners) in the sloped border (e.g. green and yellow one )? I have written the code which gives me the x and y coordinates of 4 corners of the initial rectangular grid but have difficulty with the rest. Any idea/suggestion is highly appreciated.

The order of element numbering is shown in the picture.

nye=10;
nxe=20;
T=is the matrix including the centre points of the grid elements in the order of numbering
dy=1;
nm=0;
%The slope angle is 30(deg);
n=0;
dx=1.433;
for y=1:nye
    n=n+1;
    for x=1:nxe

    nm=nm+1;
    coord(1, 1)=T(x,1)-(0.5*dx); %four corners of an element in the rectangular random field in the x direction
    if x==1 
    coord(1,1)=0;
    end
    coord(2, 1)=T(x,1)+(0.5*dx);
    coord(3, 1)=coord(2, 1);
    coord(4, 1)=coord(1, 1);

coord(1, 2)=T(x+(nxe*(n-1)),2)-(0.5*dy); %four corners of an element in the rectangular random field in the y direction
    coord(2, 2)=coord(1, 2);  
    coord(3, 2)=T(x+(nxe*(n-1)),2)+(0.5*dy);
    coord(4, 2)=coord(3, 2);
      
        for i=1:4
            g_coordx(nm, i) = coord(i, 1); %in the x direction(horizontal)
            g_coordy(nm, i) = coord(i, 2); %in the y direction (vertical)
        end
    end
    end

[DGM]

What is the condition which dictates whether an element is inside the ROI? If it's the inclusion of the element's center, then you may be able to use something like poly2mask().

[Pooneh Shah Malekpoor]

There is no condition. I should consider the remains of the elements intersected by the slope (within solid region or to the left of the sloped line) and omit the right-hand-side elements

Answer 1

[Voss]

nye=10;
nxe=20;
dx=1.433;
dy=1;
% T=is the matrix including the centre points of the grid elements in the order of numbering

Here I construct T:

xc = dx*((1:nxe)-0.5);
yc = dy*((1:nye)-0.5);
[xc,yc] = ndgrid(xc,yc);
T = [xc(:) yc(:)];

Here I construct g_coordx and g_coordy from T, replacing your code above:

g_coordx = T(:,1) + [-1 1 1 -1]*dx/2;
g_coordy = T(:,2) + [-1 -1 1 1]*dy/2;

Store the points defining the boundary between the two regions:

pts = [0 10; 10 10; 18.66 5; 28.66 5];

Now, plot the boundary and the grid corners:

figure
plot(pts(:,1),pts(:,2),'k','LineWidth',2)
hold on
plot(g_coordx,g_coordy,'g.')

Now some logic to determine which grid corners are "beyond" (to the upper-right of) the boundary.

First, dealing with the sloping line:

m = (pts(3,2)-pts(2,2))/(pts(3,1)-pts(2,1)); % slope of the line
f = @(x)m*(x-pts(2,1))+pts(2,2);             % function representing the equation of the line

% for any x-coordinate, if the corresponding y-coordinate along the line
% (i.e., f(x)) is less than the actual y-coordinate of the point, then the
% point is above the line.
% additionally, constrain this to those points whose x-coordinate is less
% than pts(3,1).
is_above_slope = f(g_coordx) < g_coordy & g_coordx < pts(3,1);

A plot to visualize this:

figure
plot(pts(:,1),pts(:,2),'k','LineWidth',2)
hold on
plot(g_coordx,g_coordy,'g.')
plot(g_coordx(is_above_slope),g_coordy(is_above_slope),'r.')

Now determine which points are above the flat part of the border:

is_above_flat = g_coordx >= pts(3,1) & g_coordy > pts(3,2);

Another plot for this part:

figure
plot(pts(:,1),pts(:,2),'k','LineWidth',2)
hold on
plot(g_coordx,g_coordy,'g.')
plot(g_coordx(is_above_flat),g_coordy(is_above_flat),'r.')

Another plot to put those two together:

figure
plot(pts(:,1),pts(:,2),'k','LineWidth',2)
hold on
plot(g_coordx,g_coordy,'g.')
plot(g_coordx(is_above_slope | is_above_flat),g_coordy(is_above_slope | is_above_flat),'r.')

Finally, remove those points from g_coordx and g_coordy:

g_coordx(is_above_slope | is_above_flat) = [];
g_coordy(is_above_slope | is_above_flat) = [];

And a final plot to check that we have the correct points remaining:

figure
plot(pts(:,1),pts(:,2),'k','LineWidth',2)
hold on
plot(g_coordx,g_coordy,'g.')