Use inpolygon command for multiple polygon areas

Dear all,

I have some difficulty with assigning points into multiple square polygons. I have data files containing (x,y,z) coordinates of every particle. With the command inpolygon I locate the particles that are inside a specified polygon:

%% Loading the data
rhoPart = 2540;

files = dir(fullfile(uigetdir,'*.data*'));
[~,Index] = natsort({files.name});
files = files(Index);
expData = cell(length(files),1);
k = 1;
for i = 1:length(files)
    fid = fopen(fullfile(files(i).folder,files(i).name),'r');
    %% Reading the data
    % Read all the data from the file
    dataRead = textscan(fid,'%f %f %f %f %f %f %f %f %f %f %f %f %f %f','HeaderLines',1);
    frewind(fid);
    % Write headerline N, time, xmin, ymin, zmin, xmax, ymax, zmax
    runData{k} = strsplit(fgetl(fid), ' ');
    % Write only the x, y, and z components of the particles, particle radius,
    % z component+ particle radius and volume of the particle
    expData{k} = [dataRead{1}(:,1) dataRead{2}(:,1) dataRead{3}(:,1) dataRead{7}(:,1) dataRead{3}(:,1)+dataRead{7}(:,1) rhoPart*(4/3)*pi*(dataRead{7}(:,1).^3)];
    % Write only the vx,vy,vz of the particles and magnitude
    velData{k} = [dataRead{4}(:,1) dataRead{5}(:,1) dataRead{6}(:,1) sqrt(dataRead{4}(:,1).^2 + dataRead{5}(:,1).^2 + dataRead{6}(:,1).^2)];
    fclose(fid);
    k = k + 1;
end

%% Classify (into a structure)
% Number of simulation runs
N_run = 1;
% Number of .data files to be extracted of every simulation run
N_files = 2745;

k = 1;
for i = 1:N_run
    for j = 1:N_files
        runData_r{j,i} = str2double(runData{k});
        expData_r{j,i} = expData{k};
        velData_r{j,i} = velData{k};
        k = k + 1;
    end
end

% Polygon
xc = expData_r{1800,1}(:,1);
zc = expData_r{1800,1}(:,1);
xv = [-0.0184 -0.0061 -0.0061 -0.0184 -0.0184];
zv = [0.2857 0.2857 0.2755 0.2755 0.2857];
in = inpolygon(xc,zc,xv,zv);
figure
plot(xv,zv,'LineWidth',2)
hold on
plot(xc(in),zc(in),'r+')

I get the following results:

But I want to do this procedure for a whole mesgrid. Thus, locate the particles of every gridbin (i.e. multiple polygons). How can I do that?

Answer 1

[KSSV]

Refer this: https://in.mathworks.com/matlabcentral/answers/592294-data-picking-from-grid?s_tid=srchtitle

[Tessa Kol]

Thank you so much. With the code below I got the following graph:

%% Velocity
close all

% Number of bins over the x-axis
Nx = 50;
% Number of bins over the z-axis
Nz = 50;
% Bins defined over the width of the silo
bins_x = linspace(-0.3,0.3,Nx);
% Bins defined over the height of the silo
bins_z = linspace(0,0.5,Nz);

[X,Z] = meshgrid(bins_x,bins_z);

xc = expData_r{1800,1}(:,1);
zc = expData{1800,1}(:,3);
plot(X,Z,'r',X',Z','r')
hold on
for i = 1:Nz-1
    for j = 1:Nx-1
        P = [X(i,j) Z(i,j) ; 
             X(i,j+1) Z(i,j+1) ; 
             X(i+1,j+1) Z(i+1,j+1) ;
             X(i+1,j) Z(i+1,j)] ;
         
        idx = inpolygon(xc,zc,P(:,1),P(:,2)) ; 
        iwant = [xc(idx) zc(idx)] ; 
        plot(xc(idx),zc(idx),'.')
        drawnow
    end
end

Now every (x,z) coordinate of the particle is devided into a bin. But every particle has also a certain velocity which are listed in velData.

How can I find (index) the corresponding velocities (vx, vz) of these coordinates?

[KSSV]

As you have the indices of the points, the respective indices can be used for velcoities too.

[KSSV]

If you are looking for the points which are not in the grid..you need to interpolate...read about interp2.

[Tessa Kol]

I think I got it wrong though (I added iwant2):

%% Velocity
close all

% Number of bins over the x-axis
Nx = 50;
% Number of bins over the z-axis
Nz = 50;
% Bins defined over the width of the silo
bins_x = linspace(-0.3,0.3,Nx);
% Bins defined over the height of the silo
bins_z = linspace(0,0.6,Nz);

[X,Z] = meshgrid(bins_x,bins_z);

xc = expData_r{1800,1}(:,1);
zc = expData_r{1800,1}(:,3);
vx = velData_r{1800,1}(:,1);
vz = velData_r{1800,1}(:,3);
plot(X,Z,'r',X',Z','r')
hold on
for i = 1:Nz-1
    for j = 1:Nx-1
        P = [X(i,j) Z(i,j) ; 
             X(i,j+1) Z(i,j+1) ; 
             X(i+1,j+1) Z(i+1,j+1) ;
             X(i+1,j) Z(i+1,j)] ;
         
        idx = inpolygon(xc,zc,P(:,1),P(:,2)) ; 
        iwant = [xc(idx) zc(idx)] ;
        iwant2 = [vx(idx) vz(idx)];
        plot(xc(idx),zc(idx),'.')
        plot(vx(idx), vz(idx),'o')
        drawnow
    end
end

Ultimately I want to make a quiver plot of the average velocity in a grid. So my goal is to take the mean of vx in every grid bin and to tak the mean of vz in every grid bin. The same thing for the coordinates.

And then I can do this:

quiver(xc_avr, zc_avr, vx_avr, vz_avr)

[KSSV]

You can use mean to get the average velcoities right?

[Tessa Kol]

Yes, I check if the indexing went correctly and it seems fine to me. So now I have to do indeed with mean. I worked with the code below:

vz_mean = cellfun(@(x)mean(x(:,2)),iwant2);
vx_mean = cellfun(@(x)mean(x(:,1)),iwant2);

However, I can still not use quiver since X and Z do not match with vz_mean and vx mean. Because in the for loop iwant2 is calculated for 49 x 49 bins and X and Z are 50 x 50 bins. How can I resolve this?

[KSSV]

Why do you think the indices are wrong?

[Tessa Kol]

Everything is solved by now. I succesfully made a quiver plot. Thank you for you help!

Answer 2

[Matt J]

Use discretize(),

https://www.mathworks.com/help/matlab/ref/double.discretize.html

Apply it separately to all of your xv's and then on all your zv's.