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

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.

Answer 1

[John D'Errico]

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.

[Ashfaq Ahmed]

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]

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,:)

Answer 2

[Torsten]

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

Answer 3

[Walter Roberson]

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.

[John D'Errico]

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.