Using a for loop to change a matrix

Hello everyone

Okay, so I have a 40x3 matrix with nodal point IDs. This can then be broken down into x and y coordinates per nodal point. Each line in the 40x3 matrix represents a triangle. I'm trying to make node identification go in a counter-clockwise direction. I've used the equation in the code below to calculate the polygon (triangle) area. For those values that are positive, the direction is counter-clockwise. For negative values it is clockwise. Its roughly 50/50.

What I'm trying to do is use a for loop, so some other method, to make the nodal identification direction change. So instead of having x1 y1, x2 y2, x3 y3, I want x1 y1, x3 y3, x2 y2. I've been mulling over this for a few days, looking online to see if anyone has had similar ideas but I can't find anything. For what I have down below, I keep getting an error. And obviously that wouldn't be the full way of doing it, but it was a way I was going to try.

Any help would be greatly appreciated.

My code:

%%Determing boundary nodes
Tri.triangles=R.ci;

Tri.triX=x(Tri.triangles); Tri.triY=y(Tri.triangles);
x1=Tri.triX(:,1); y1=Tri.triY(:,1);
x2=Tri.triX(:,2); y2=Tri.triY(:,2);
x3=Tri.triX(:,3); y3=Tri.triY(:,3);

% direction of triangle circulation, +ive is counterclockwise
for i=1:length(Tri.triangles);
Tri.area(i)=0.5*((x1(i)*y2(i))-(x2(i)*y3(i))+(x2(i)*y3(i))-(x3(i)*y2(i)));
T_area=Tri.area';
t_neg(i)=T_area(i)<0; T_neg=t_neg';
end
T_neg=double(T_neg);  

bb=x1(T_neg(x1)>0);

Error:

Subscript indices must either be real positive integers or logicals.

T_area=

-0.338828523177654
-0.565042113652453
-0.903958166425582
0.00953685300191864
-0.271408101485577
-0.670282627630513
0.370416248624679
0.260899886256084
-0.225398410111666
-0.00953673134790734
-0.826227724959608
0.826595614606049
-0.418818252452184
0.117110523569863
0.0329695575637743
-0.0329722759779543
-1.57312464073766
-0.658750543370843
0.658923690207303
-0.155193869140930
-0.0884571871720254
0.275831339182332
0.0546458313474432
-0.325210246955976
-0.117114751192275
0.286314578843303
0.494672991335392
0.569307863712311
0.871622048784047
0.121534811332822
1.57387580676004
0.169117909390479
-0.324359477497637
-0.185630145482719
-0.402470884844661
0.402549884282053
-0.286296253907494
0.0546600596280769
0.403549019247294
0.0615939649869688

Answer 1

[dpb]

Given a function that returns the vector of areas given the node and coordinate arrays, then

ix=triArea(nodes,coordinates)<0;          % the rows that fail
nodes(ix,[2:3])=fliplr(nodes(ix,[2:3]));  % swap the offending ids

should do the trick, methinks...

ADDENDUM

OK, am back from meeting and had chance to look into the guts a little more. The area formula as modified in your recent comment still isn't correct; the vector cross formula for area is A=|AB x AC|/2. Algebraically if name the points ABC-->123 this reduces to

AB=[x2−x1,y2−y1]; AC=[x3−x1,y3−y1];  % vectors 
A=[(x2−x1)(y3−y1)−(x3−x1)(y2−y1]/2;  % Area of ΔABC

You can write this in Matlab as a function handle succinctly as:

>> fn_triArea=@(x,y) ((x(:,2)-x(:,1)).*(y(:,3)-y(:,1))-(x(:,3)-x(:,1)).*(y(:,2)-y(:,1)))/2;

In use it's simply

>> A=fn_triArea(x,y)
A =
 1.0e-03 *
  0.3445
  0.2559
  0.2817
  0.2110
  0.3301
  0.2789
>>

For comparison,

>> [A polyarea(x',y').']
ans =
 1.0e-03 *
  0.3445    0.3445
  0.2559    0.2559
  0.2817    0.2817
  0.2110    0.2110
  0.3301    0.3301
  0.2789    0.2789
>>

we note the same values as the builtin polyarea function which gives confirmation we got the formula right. I did a second and third comparison as well of the area of the total rectangle bounding the triangle less the three triangles' outside the inscribed one and the direct cross calculation as well--all were in agreement. I didn't check your algebra but you've still got a transcription or sign error somewhere.

Note the x, y coordinate arrays you've given all result in positive areas, if we pick a couple at random and flip coordinates what happens?

>> ix=randi(length(x),2,1)  % two random points to flip coordinates
ix =
   2
   5
>> x(ix,[2:3])=fliplr(x(ix,[2:3]));  % and flip 'em...
>> y(ix,[2:3])=fliplr(y(ix,[2:3]));
>> [fn_triArea(x,y) A polyarea(x.',y.').']
ans =
 1.0e-03 *
  0.3445    0.3445    0.3445
 -0.2559    0.2559    0.2559
  0.2817    0.2817    0.2817
  0.2110    0.2110    0.2110
 -0.3301    0.3301    0.3301
  0.2789    0.2789    0.2789
>>

Note the two selected triangles now do have negative areas and again compare numerically with the previous results.

Of course, you also need to swap the node numbers to match but if you're storing the full x,y arrays and computing based on them and their row indices instead of actual node ID then the area computation doesn't even need the node; it's just an external bookkeeping addon.

But, there's no indexing error and the logic described will fix your problem implemented correctly.

HTH...

PS: The above results are from your x,y arrays you pasted in in a comment some time earlier that are roughly -6.3 and 56.8 with small differences in values, hence the small areas aren't surprising.

ADDENDUM 2

As one last illustration on the area, the worst possible way to compute anything almost always, but we'll trot it out for further confirmation--

>> det([x(1,:)',y(1,:)',ones(3,1)])/2
ans =
 3.4450e-04
>>

[Meghan]

That's not working for me. My node matrix is 40x3 but my coordinates are 40x1, but there are 6 of them. I even tried using 40x3 x-coordinates but I'm still getting an 'Index exceeds matrix dimensions' error.

I get what the code is trying to do, I'm just not sure I can apply it to what I want to do. Or else I'm just not doing something which is blatantly obvious, which is more likely!

[dpb]

All you should have to do is swap the node IDs in the 40x3 array; everything else should be written in terms of lookup of those ID numbers.

[Meghan]

Yea my nodal point ID matrix is still there, eg (Tri.triangles):

13972  13589  13971
13971  14310  14311
13972  13971  14311
14312  13972  14311
13972  14312  13973
13973  14312  14313

But my coordinates are split into x and y matrices:

x=
-6.35311889648438  -6.37179565429688  -6.34118652343750
-6.34118652343750  -6.31210327148438  -6.32128906250000
-6.35311889648438  -6.34118652343750  -6.32128906250000
-6.32095336914063  -6.35311889648438  -6.32128906250000
-6.35311889648438  -6.32095336914063  -6.34356689453125
-6.34356689453125  -6.32095336914063  -6.31997680664063

y= 
56.8188362121582  56.7914733886719  56.7994270324707
56.7994270324707  56.7953987121582  56.8142700195313
56.8188362121582  56.7994270324707  56.8142700195313
56.8274803161621  56.8188362121582  56.8142700195313
56.8188362121582  56.8274803161621  56.8419303894043
56.8419303894043  56.8274803161621  56.8515205383301

Code:

    ix=T_area(Tri.triangles,Tri.triX)<0;
    Tri.triangles(ix,[2:3])=fliplr(Tri.triangles(ix,[2:3]));

Error: Index exceeds matrix dimensions.

Error in Extraction_MRv14_WORKS (line 262)
ix=T_area(Tri.triangles,Tri.triX)<0;

[dpb]

That's a lookup error in the T_area code that isn't shown, it looks like. It doesn't seem there can be an index error in the actual code line itself unless somehow the length of the two arrays isn't commensurate so ix exceeds the size somehow--but that doesn't seem possible on the surface, either.

Can you make a short demo complete test case that can be run by cut 'n paste? I'm thinking the lookup logic is the culprit but I can't see/test that here...

ADDENDUM

OK, I think I see the picture...not sure where you made the error but I was presuming the coordinates were stored only by node ID but I gather from looking at the data you've duplicated the values for every element and are still referring to them by position, not by node.

I'd probably try to revert to not duplicating data unnecessarily but perhaps it is useful in the end for the full storage, I don't know.

BUT, there should be no issues in swapping x and y array columns by index as well to coincide with the node.

The trivial way to write the area (albeit not most rapid, probably) is

polyarea(x.',y.')

This also may produce only positive values, I'm not sure.

I hadn't looked at that part before, indeed the results are different than you're getting. I checked the result several ways and indeed your formula is in error. I've got a meeting; gotta' run, so not sure where it is...

[David Goodmanson]

Along the way, the triangle area formula needs to be corrected. You have four terms and unfortunately the middle two cancel each other, leaving two. There are in fact six terms, which you can get for example in the 'using coordinates' part of the Wikipedia 'triangle' article and dropping the abs value signs there.

[Meghan]

Thanks David, I had just made a slight error in the formula. Fixed it now.

dpb, I've written the code using randomly generated numbers.

%%Determing boundary nodes
Tri.triangles=randn(40,3);

Tri.X=randn(40,3); Tri.Y=randn(40,3);

x1=Tri.X(:,1); y1=Tri.Y(:,1);
x2=Tri.X(:,2); y2=Tri.Y(:,2);
x3=Tri.X(:,3); y3=Tri.Y(:,3);

% direction of triangle circulation, +ive is counterclockwise
for i=1:length(Tri.triangles);
Tri.area(i)=0.5*((x1(i)*y2(i))-(x2(i)*y1(i))+(x2(i)*y3(i))-(x3(i)*y2(i)));
T_area=Tri.area';
t_neg(i)=T_area(i)<0; T_neg=t_neg';
end
T_neg=double(T_neg); 

ix=T_area(Tri.triangles,Tri.triX)<0;
Tri.triangles(ix,[2:3])=fliplr(Tri.triangles(ix,[2:3]));

I am still getting an error, just a different one this time!

Subscript indices must either be real positive integers or logicals.

Error in test (line 18)
ix=T_area(Tri.triangles,Tri.triX)<0;

[dpb]

That's cuz as I said originally "given a function" to calculate area. You've supplied the coordinates as indices into an array above, not as an argument to a function.

See the amended Answer above for the actual paradigm; "you don't need no steenkin' loops!" :) to do the computations or the fixup.

[Meghan]

Thanks so much for taking the time to look at it. I just used polyarea to give me what I want, I must have transcribed the equation wrong.

Hopefully now it will work with the rest of my code! Thanks :)

[dpb]

No problem, glad could help.

Except note that polyarea doesn't return the sign so if the point is the orientation, it doesn't help. I wrote the correct solution above.

[Meghan]

Yea I changed it to be the equation rather than the function, thanks.

I wonder if you can help me with something else actually. Its not really related to the above.

Basically I have a 47x5 matrix with node ID, x, y, z, and a code. I've used unique on it but it shows that I have some node IDs that are the same except for the code number. Do you know how I could remove the duplicate rows, all while keeping the highest code number?

My code and the first few lines of the matrix.

nodes=unique(sortrows([Tri.boundNodes;Tri.interNodes]),'rows');

13589  -6.3717957  56.791473  82.498566  2
13971  -6.3411865  56.799427  48.073643  0
13971  -6.3411865  56.799427  48.073643  2
13972  -6.3531189  56.818836  33.481136  0
13972  -6.3531189  56.818836  33.481136  2
13973  -6.3435669  56.841930  45.153000  0
13973  -6.3435669  56.841930  45.153000  2
14310  -6.3121033  56.795399  46.835625  0
14310  -6.3121033  56.795399  46.835625  2
14311  -6.3212891  56.814270  42.117233  0

[dpb]

Return a list of just the unique nodes (use nodes(:,1) for input) and then either just loop over that list returning the collection for which u(i)==nodes(:,1) and find the max of the code number in each set, returning the optional index variable which will be the index within the set of the desired row. The index array from the unique call will let you relate that row back to the row in the original array.

Or, you can get more sophisticated and take the above outputs from the second unique call and bin them with histc from which you can find the subset of duplicates. accumarray can then be used with an anonymous function including max for each of these groups to do the function of the above "brute force" loop construct. In reality, the loop probably will run faster than the accumarray solution and is simpler to code but by some measures the latter is "more elegant".

[Meghan]

Thanks :) I was thinking along the same lines, just wanted to see if you had a simple solution I might have missed. I tend to overthink a lot of things!

I very much appreciate all the help :)