Clear Filters
Clear Filters

How can I find out the index of the points here?

20 views (last 30 days)
HI! I have this interesting indexing question.
I have a cloud data set (attached, Cloud.mat). Plotting the X and Y, I get this figure. My question is - I want to find all the indices that falls under these four sections. Can any one please tell me how can I get the indices? Any feedback will be much appreciated.

Accepted Answer

John D'Errico
John D'Errico on 17 Feb 2023
Edited: John D'Errico on 18 Feb 2023
Pretty easy. (Easier yet, IF the lines were drawn at right angles. If they were, then we could do it in literally one line of code.) But I can see they are NOT exactly perpendicular to each other, so this will take a couple more lines of code, but not many.
Anyway, just create a simple triangulation.
triXY = [-1.5 -1.5;1.5 -1.5;-1.5 2;1.5 2;0 0.25]
triXY = 5×2
-1.5000 -1.5000 1.5000 -1.5000 -1.5000 2.0000 1.5000 2.0000 0 0.2500
triangles = [1 2 5;1 3 5;2 4 5;3 4 5];
T = triangulation(triangles,triXY)
T =
triangulation with properties: Points: [5×2 double] ConnectivityList: [4×3 double]
xy = rand(10,2)*3 - 1.5; % some random points.
trimesh(triangles,triXY(:,1),triXY(:,2))
hold on
plot(xy(:,1),xy(:,2),'ro')
axis equal
Next, we consider what tools a triangulation can use.
methods(T)
Methods for class triangulation: barycentricToCartesian edgeAttachments featureEdges isConnected pointLocation vertexAttachments cartesianToBarycentric edges freeBoundary nearestNeighbor size vertexNormal circumcenter faceNormal incenter neighbors triangulation
In there, we see that pointLocation is a method for this class. That may help me.
help pointLocation
--- help for triangulation/pointLocation --- pointLocation Triangle or tetrahedron containing specified point TI = pointLocation(T, QP) returns the index of the triangle/tetrahedron enclosing the query point QP. The matrix QP contains the coordinates of the query points. QP is a mpts-by-ndim matrix where mpts is the number of query points and 2 <= ndim <= 3. TI is a column vector of triangle or tetrahedron IDs corresponding to the row numbers of the triangulation connectivity matrix T.ConnectivityList. The triangle/tetrahedron enclosing the point QP(k,:) is TI(k). Returns NaN for points not located in a triangle or tetrahedron of T. TI = pointLocation(T, QX,QY) and TI = pointLocation (T, QX,QY,QZ) allow the query points to be specified in alternative column vector format when working in 2D and 3D. [TI, BC] = pointLocation(T,...) returns, in addition, the Barycentric coordinates BC. BC is a mpts-by-ndim matrix, each row BC(i,:) represents the Barycentric coordinates of QP(i,:) with respect to the enclosing TI(i). Example 1: % Point Location in 2D T = triangulation([1 2 4; 1 4 3; 2 4 5],[0 0; 2 0; 0 1; 1 1; 2 1]); % Find the triangle that contains the following query point QP = [1, 0.5]; triplot(T,'-b*'), hold on plot(QP(:,1),QP(:,2),'ro') % The query point QP is located in a triangle with index TI = 1 TI = pointLocation(T, QP) Example 2: % Point Location in 3D plus barycentric coordinate evaluation % for a Delaunay triangulation x = rand(10,1); y = rand(10,1); z = rand(10,1); DT = delaunayTriangulation(x,y,z); % Find the tetrahedra that contain the following query points QP = [0.25 0.25 0.25; 0.5 0.5 0.5] [TI, BC] = pointLocation(DT, QP) See also triangulation, triangulation.nearestNeighbor, delaunayTriangulation. Documentation for triangulation/pointLocation doc triangulation/pointLocation Other uses of pointLocation DelaunayTri/pointLocation
pointLocation(T,xy)
ans = 10×1
1 1 4 4 1 1 3 2 1 3
The index returned indicates where the corresponding point fell. You can skip the plots, as you wish.
A second solution is also not difficult. This relies on an affine transformation, if the points in the plane. Think of the data as a coordinate system around an origin at the point (0,0.25). Thus, where the lines cross. A problem is, the lines have the wrong slopes, as they are not perpendicular to each other. We can resolve that simply enough, by scaling y by a factor of 3/3.5. And then we can rotate the coordinate system by 45 degrees.
The final transformation, if we have an array xy of size nx2:
ang = -45; %degree rotation, so a clockwise rotation of 45 degrees
T = @(xy) (xy - [0 0.25]).*[1 3/3.5]*[cosd(ang) sind(ang);-sind(ang) cosd(ang)];
THe function handle T uses capabilities introduced in R2016b. Earlier releases would use bsxfun.
If you look at the transformation, first notice I subtract off 0.25 from y only, so x is unchanged. Effecively, now the lines cross at the origin.
Next, I scaled the problem, so multiplying y by a factor of 3/3.5. That effectively corrects the two lines so they are now perpendicular after the rotation. (This was the error @Walter Roberson made in his suggestion.)
Finally, the 2x2 matrix I multiply by there is a rotation matrix. I've used sind and cosd there so it is clear how it was formed. See what the transformation did to the points listed in triXY.
triXY
triXY = 5×2
-1.5000 -1.5000 1.5000 -1.5000 -1.5000 2.0000 1.5000 2.0000 0 0.2500
T(triXY)
ans = 5×2
-2.1213 0 0 -2.1213 0 2.1213 2.1213 0 0 0
The 5th point now lies at the origin. And each of the other points have been properly scaled, and then rotated. Now, how can we use the above transformation? I'll rotate the points in the array xy, into a new array, then plot the transformation.
newxy = T(xy);
figure
plot(xy(:,1),xy(:,2),'o')
hold on
trimesh(triangles,triXY(:,1),triXY(:,2))
plot(newxy(:,1),newxy(:,2),'x')
plot([xy(:,1),newxy(:,1)]',[xy(:,2),newxy(:,2)]','r-')
axis equal
hold off
grid on
As you can see, the entire system of points has been shifted, stretched just a bit, and then rotated.
How can we now use this? If the final location of a point was originally in the top triangle, then it is now in the first quadrant. In that case, both the new x and the new y are positive numbers. So now we can decide which of the original (quasi-)quadrants a point lies in, merely by testing the signs of the points after rotation.
T(xy) >= 0
ans = 10×2 logical array
0 0 0 0 1 1 1 1 0 0 0 0 1 0 0 1 0 0 1 0
So if both results are 1 there, then the point lies now in the first quadrant. If both results are zero, then the point lies in the 3rd quadrant. And then you can easily distinguish the other quadrants too. A nice trick is to use a binary to decimal encoding.
(T(xy)>0)*[1;2]
ans = 10×1
0 0 3 3 0 0 1 2 0 1
As you can see, that yields an integer, so 0,1,2,3, depending on which of the original regions a point fell in. I'll let you decode which of the regions was which.
  2 Comments
Ashfaq Ahmed
Ashfaq Ahmed on 21 Feb 2023
I think they are the locations (four subsets) within the bounding box. But then the next question arises, how can I connect them with the xy values?
For example, I want to have four variables like, var1, var2, var3, and var4.
var1 will contain the xy values that falls under subsection 1. And same goes for va2, 3 and 4.
Ashfaq Ahmed
Ashfaq Ahmed on 21 Feb 2023
Ok, I think I could solved the problem I mentioned in my comment. I wrote these few lines of codes.
Section_List = pointLocation(T,xy);
sec1 = find(Section_List==1);
sec2 = find(Section_List==2);
sec3 = find(Section_List==3);
sec4 = find(Section_List==4);
TIDE_LOW = Landsat_Time(sec1,:);
TIDE_EBB = Landsat_Time(sec2,:);
TIDE_FLOOD = Landsat_Time(sec3,:);
TIDE_HIGH = Landsat_Time(sec4,:)

Sign in to comment.

More Answers (2)

Torsten
Torsten on 17 Feb 2023
Edited: Torsten on 17 Feb 2023
X = [0,10,0,-10,0,5,0,-5]
X = 1×8
0 10 0 -10 0 5 0 -5
Y = [10,0,-10,0,5,0,-5,0]
Y = 1×8
10 0 -10 0 5 0 -5 0
syms x y
g1 = (y-2)/(x-1.5) == (-1.5-2)/(-1.5-1.5);
g2 = (y-2)/(x-(-1.5)) == (-1.5-2)/(1.5-(-1.5));
g1 = matlabFunction(solve(g1,y))
g1 = function_handle with value:
@(x)x.*(7.0./6.0)+1.0./4.0
g2 = matlabFunction(solve(g2,y))
g2 = function_handle with value:
@(x)x.*(-7.0./6.0)+1.0./4.0
section1 = find(Y >= g1(X) & Y >= g2(X))
section1 = 1×2
1 5
section2 = find(Y < g1(X) & Y > g2(X))
section2 = 1×2
2 6
section3 = find(Y <= g1(X) & Y <= g2(X))
section3 = 1×2
3 7
section4 = find(Y > g1(X) & Y < g2(X))
section4 = 1×2
4 8

Walter Roberson
Walter Roberson on 17 Feb 2023
Subtract the center of mass. Rotate the result by +45 degrees. Now the combinations of sign() of the transformed x and y coordinates tells you which sector you started from.
  1 Comment
John D'Errico
John D'Errico on 18 Feb 2023
Edited: John D'Errico on 18 Feb 2023
Except that the lines drawn are not perpendicular to each other. It would be nice if they were, but they are not.
Check the slopes of the lines. One of them has a slope of
3.5/3
ans = 1.1667
The other slope is
-3.5/3
ans = -1.1667
If they were perpendicular, the slopes would be the negative reciprocal of each other.
You could still perform an appropriate transformation that would make the idea work. But not just a 45 degree rotation.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!