 Welcome back, CAD from Scratch, episode 17. The topic today is 2D landing triangulations with edge constraints. Before we get into that though, I want to quickly address some of the new viewers. Couple of my other videos went viral for the past few days and my subcount basically quadrupled. So basically this series is my pet project where we're implementing CAD routines from Scratch in C. CAD stands for computer aided design. It's basically a tool that engineers and designers will use to make geometries for analysis, for simulation, or for manufacturing, like maybe 3D printing, or order cutting or whatever. And in these videos I'm basically implementing a CAD program piece by piece, talking about the process. The only caveat being that everything in this series is done entirely from scratch. So there's no dependencies, there's no libraries, stuff like that. So if you run X11 on Linux, take the code off GitHub and compile it. It'll work no problem right out of the box. And all the code is explicitly public domain, so that's no problem. And right now we're in the beginning phases still. We can import models, export models, do some simple operations on models. We can query things, we can change the view in different ways. We can pan the model around. So, simple stuff is working. We have a lot of work to do though. I spent the past few videos working on like a few functions behind the scenes and we'll do that more in the next few videos. But in a few weeks we should have most of the core functionality set up for this software. And so if you're like a Python bro and all you do is you spend six hours a day pipping installing libraries and stuff, you're not gonna enjoy this series I don't think. But if you like knowing how stuff works then maybe give it a chance if it interests you. And one last thing, I always put a timestamp in the description so you can just skip to see the final result. It may help you get a better idea of what's going on if you're a little bit confused in the beginning. Without any further ado, let's get into the actual content. So what is a constrained triangulation? Well, basically let's say you had a point set, these white dots here. In the previous video we talked about how to generate a delineate triangulation of that point set. That would be this yellow grid of triangles here. And one good thing about this triangulation is that all the elements follow this special condition I'll talk about that condition later but basically it makes all the elements have a very nice aspect ratio and be almost equi-angular, be almost equilateral triangles. And so that is very, very nice. The thing is though, there are some cases where you want to impose a constraint on the triangulation that is you want this red line to actually be an edge of the triangulation. Now sometimes that edge will already be part of the triangulation. It will already be a member of the triangles here. When it's not though, you have to somehow cut through triangles and replace them in a way that this red line is now a suitable member of the triangulation. And that's what we're gonna talk about in this video. And to do that we're gonna be kind of adapting some of the work done by Sloan in this paper titled A Fast Algorithm for Generating Constrained Delineate Triangulations which is a very good paper. Take a look at it, it's available online. You can access it, no problem. And I put it upside down here in solidarity with the plight of the Australian people because that's where this guy is from. So, power to people man, good luck. So why would you want to implement this? Well, I can think of a couple reasons why you might want to do this. Some real world examples and some mathematical ones that we're gonna implement in this video but the first one would be contouring through a defined feature. Let's say you went on a hike and you carried your phone around and you recorded the GPS coordinates and the altitude of a couple parts of your hike like on a mountain or on a valley or in a flower field or whatever. And you had these data points here, these white dots. And you went home and you made a contour map of your journey using this method for triangulation. And that's very, very good for that purpose. I mean, that's one of the reasons why people use this kind of triangulation is because it makes a very nice contour map where you don't have very slim triangles. So it's very good for that purpose. However, let's say there was a river here along this red line and you know for a fact that river is at sea level. Well, some of these triangles edges may not be passing through zero in your contour map unless you were to impose that this red line be a triangle edge at sea level and that would fulfill that purpose, right? So that's one reason why you would wanna, you know use this method of triangulation. Number two, what about path planning? So let's say you were, I don't know, the Mars rover and you were roving around Mars taking soil samples at these white points. And you also had to very quickly measure some kind of wind data along this red line in the same journey, right? And so normally you'd make a triangulation to figure out your path planning route. However, you also want this red line to be part of your route because you have to take measurements along that line. And so you can use this methodology to force the red line to be part of your path, right? And the very last thing is what we're gonna do in this video and I'm just slicing a mesh. So when you talk about, you know, constructive solid geometries, you basically have bodies overlapping other bodies and locations where they intersect, those will eventually be constraints on our triangulations that we will have to implement in the next couple of videos. So I'll talk about it later, but know that this is on for a reason. So what is the algorithm to do so? So the first step of the algorithm is just to triangulate normally. To do what we did in episode 16 and just form this yellow grid of triangles from these white point clouds, right? Once you've done that, we can loop over every constraint edge. So every one of those red edges. If the edge is not already part of the triangulation, you know, like for example, this one is not part of the triangulation. If not, record all the edges that it intersects through. So for example, this red edge slices through one, two, three other edges of the triangulation. So you have to record those in this sort of list here. Then for every one of the intersections, you do a couple of things. First, if the triangles that share that edge, just remember every edge inside of this triangulation shares two triangles. So if those triangles make a convex quadrilateral, we will swap their diagonal. I'll talk about more of what that means in a minute, but just know that for now we're gonna swap the diagonals of convex quadrilaterals. Otherwise we skip and go back to step two. Now, if that new diagonal that we just swapped still intersects the constraint, we're gonna add that to the same list of edge intersections. Otherwise, if it doesn't cross the red line, for example, this orange line here is a new edge, right? If it does not cross the red line, we're gonna add that to this of new edges. And for all those new edges, we're gonna loop over them until nothing changes in the triangulation. And we're gonna check if the edge is not a constraint and it still violates the delineary condition, swap the diagonal. That's very complicated words. I don't exactly get what that means, but we'll go through it piece by piece and talk about how it operates. So what is this delineary condition on the right here? Well, we talked about this in a previous video, but basically, let's see how triangles with the points V1, V2, and V3, if you draw a circle through the points of that triangle, and you know the circumcircle, if that circle contains point P, for example, then that triangle is violated by point P in this condition, right? As well, triangle V1, V2, and P, the circle through those points, right? This circle here also includes V3. And so in the same way, V3 violates V1, P, V2. So both of these are bad failing this condition. So the solution to this is actually to swap the diagonal V1, V2 with diagonal P to V3. So you can see here how that would look. Basically, if you draw this as the new diagonal from those two triangles, the triangle through V1, V3, P does not have a circle that includes point V2, or vice versa, P, V3, V2, that circle does not include V1. So this is the good swap. This is the bad swap. So in the previous video, we talked about how we can evaluate this condition, and there's some math to do so. You can take a look at this. It's pretty simple. Basically, you evaluate pieces of it as you go and check if things are true or not as you go, and then if they're true or not, you can kind of stop the computation early. And that's kind of the benefit of using this methodology for evaluating the condition. But take a look at that and see if it makes sense. Now, in terms of swapping the diagonals, what does that actually entail? So let's see if I had two triangles, triangle L and triangle R shown here. If these two triangles require a swap, which these ones will require a swap, you have to do a few things. Obviously, you have to change which vertices are corresponding to which triangles, right? So on the left here, triangle L is P, V2, V1, but on the right here, it's P, V2, V3. You can see that the third vertex there changed identity. But also you have to keep track of the adjacencies. So keeping track of the adjacencies as you go through this methodology allows you to save a lot of time, because if you can know what triangles touch other triangles, you don't have to constantly search the entire face for the triangle that you need, you can just check, oh, what's the triangle next door? And that makes this algorithm much more efficient when you implement it. And so you have to keep track of adjacencies, you have to update them as you go and swap diagonals, you have to again swap out what triangles are touching which other triangles. That's very important to keep that in mind. If you implement this and you have any problems, the problems are gonna be most likely because you messed this part up. Always make sure you keep track of these triangles and their adjacencies very, very accurately. That's where your buckles are gonna be when you implement this. Lastly here, to find the intersections between the constraint edge and your triangulation, it's actually pretty easy. So let's go up to the top again. You can see here, here was our constraint edge and we've cut through one, two, three edges of the existing triangulation. And so the objective now is to figure out which edges we're cutting through. So the first thing you do is you start at a triangle which contains a V1. V1 is one of the vertices of your constraint edge. Now find the triangle that contains V1 because there might be a few. There could be a triangle here, there could be a triangle here, could be a triangle here. Find the triangle that has V1 that opens towards V2. I'll talk about that in a second but basically what we're gonna do is we're gonna use this helper function here which determines what direction is a point from a line segment. Anyway, back to this. Find the triangle that opens from V1 to V2. That's this triangle here. Then go around the edges of the V1 triangle. Remember, we count triangles in counterclockwise order. So the edges will go like this. So this is edge one, this is edge two, this is edge three. Basically, go around and find which edge for which V2 is to the right. So we're not to the right of edge one. We are to the right of edge two. We're not to the right of edge three. Right, as you're going around this triangle. Again, that uses this helper function here, helper function one. What direction is a point from a line segment? So we'll use that a couple times in this algorithm. Then go to this triangle. And this triangle again is stored in the adjacency array of this triangle. So it's very easy to determine what triangle is in this direction just by keeping track of what triangles are next to other triangles in general. And you're gonna repeat that over and over and over and over again until you get to a triangle that includes V2. And once you're there, you're done. We found all the intersections between constraint edge and the given triangulation. So that's the process for that. Now, let's talk about the helper function that we're gonna have to implement to make this work. So as I mentioned before, the first one determines what direction is a point from a line segment. So in this case, let's say you had line segment AB and the point of interest point P. It's very easy to hit the cross product of AB and AP if that's below zero, you're to the right. If that's above zero, you're to the left. I'm very simple. Function number two, do two line segments intersect. We implemented this, I think two videos ago for the triangles, you know, stuff. Basically though, if you could satisfy this condition, basically that requires that your line segments intersect. And the way this works basically is you're checking if the, if sort of this cross product and this cross product are in different directions. At the same time, you're checking if this cross product and this cross product are different directions. If they are, you know, check marks, check marks, then the two line segments do intersect. If they're in the same direction, you know, then there's no intersection that's possible. It's a very simple algorithm to check though. And lastly, help the function number three is a given quadrilateral convex. So this is also simple to check. Basically you go around the quadrilateral. So in terms of the side lengths here in one direction and you take the cross product of side I and side I plus one. If those two are all, sorry, if the cross product of side one and side two is the same as side two and side three and side three and side four and side four and side one. If they're all the same direction, then you are convex. If one of them is the other direction, then you're not convex. And you know, in that algorithm, you can just kind of iterate through until one of them is not the same sign anymore then you just break out of the algorithm and say, nope, you're not convex. With that, we'll go into the actual code here. And all the code is actually, first let me show you the main function. What? CD into a file. Okay, smart. Basically, I had a couple test cases here. I have eight of them. And the way this works is I'm importing to this, I'm basically passing as a parameter a couple of things, a bunch of points, a point cloud, as well as a constraint edge and how many points are there in the point cloud, right? And so for this case one, there's only one constraint edge. It passes through vertex four and six. Test case two has, you know, same point cloud, but different constraint edge, et cetera, et cetera. I have some of these that have multiple edges for constraints just to check and see how they work. And I'll test these cases out at the very end of the video and you can see how that works. So now in terms of the actual implementation of these algorithms here, the two test cases, sorry, the three helper functions are shown here. So the first one is point direction from line segment. This just implements exactly that cross product that I talked about before. If a point is, you know, to the right of AB, this function returns one. This function here returns one. So less than zero returns one. This function here is another bull function, line segments cross 2D. This implements this helper function here. If either of these two are less than zero, we, or sorry, if they're both less than zero, we implement, we say, you know, return one. If not, we return zero. And the last function to help us out is this quadrilateral convex 2D. This returns one if the quadrilateral of P1, P2, P3, P4 is convex. Again, it just completes the cross product here of the sides of the quadrilateral and returns one if it's convex and zero if it's not. Down to the actual meat of the algorithm here. We had this function last time called linear triangulate. I've changed one thing. That is, we're now passing in a list of constraint edges and this structure is defined in our jam.h header file here. Very simple, just a list of arrays. So it's a list of arrays. So you have a number of values in each member of that list and the actual array itself. So very simple. And so the way this works is the first part of this function is from episode 16. So everything from here on down is episode 16 stuff. Very, very simple stuff. Once we get out of episode 16, you know, that's the episode that basically creates this triangulation here. We have a couple of things to do. So the first thing here is we're looping through all of the constraint edges. So that's the first step of the algorithm here. You know, for each constraint edge. So we're looping through this outer while it was looping through all those constraint edges. And we're checking to find the triangles that each vertex is contained in. Now this is useful because we want to determine sort of what triangle contains V1. Oops. Well, triangle contains V1 and V2 because that will kind of determine whether or not we're, you know, at the beginning or end of our intersections. So that's the point of this. So we're doing, we're looping through that. Now we're looking through all of those triangles that contain vertex one, looking for vertex two. Or until we find an edge that crosses vertex one and vertex two. So again, that's the part of this. So we're looking for the first line that crosses this V1 and V2 constraint edge. That would be line two here. That's the point of this. And so we're using again this point direction from line segment to determine if vertex two is still left to the right of both edges on vertex one on this given triangle that we're looking through. Now, if we end up finding the edge, so if miraculously this red line is in that edge, well, we're done, right? So if that happens, we do nothing. If the edge is not found in the vertex array, we will start looking through intersections in our triangle trajectory. So that would be basically, we're starting here and we're passing through every triangle, you know, one by one till we get to vertex two. And again, we're checking if vertex two is to the right of any of the edges in that triangle. So again, we're looping through, like I said here, edge one V2 is not to the right of edge one. It is to the right of edge two. It's not to the right of edge three. So edge two would be the triangle that you'd kind of go along that direction for. That's that here. Now comes the question of looping through edges that need to be removed and adding edges onto the triangulation. So we are in this part of the algorithm here. Yeah, this part of the algorithm here, step two. So basically, we're trying to find triangles that contain a given edge. That's the first step here. And we're trying to evaluate if the quadrilateral is convex. That's quadrilateral convex 2D. We're passing in the points on the quadrilateral that contains two triangles. If they're convex, we swap the diagonal that's shown here. And again, we're keeping track of everything. So not just the vertices, but also the actual adjacencies of the triangles. And that's all of all this that's doing here, all these lines here is keeping track of the triangles and their adjacencies. Now, if that new diagonal still intersects with the constrained edge, we're adding that to our list of intersections. That's what this does here. And if it doesn't, we add that to the new edge list. So that's step two B right here. Now, it also can happen that the quadrilateral is not convex. If that's the case, we go back to this. If it's not convex, we're going back to step two for the next iteration of this for loop. That's all fine. Now, down here, this is step three, I believe. So yes, for each new edge, until nothing happens, nothing changes in the triangulation, we repeat this. So I have this variable here swaps. And while swaps is greater than zero, we will continue this iteration and this process of triangulating and fixing the triangulation. So swaps equals zero. At the end of this, we're going to increment swaps whenever a swap occurs. So the way this works is basically, you loop over all the new edges in the triangulation and if that new edge is not the constrained edge, we identify the vertices. First off, the vertices that are not, let me show an example. So for example, P and V3 are vertices that do not share the boundary of the two triangles. So these would be like the lone vertices of this quadrilateral. So we evaluate those and we can use those to kind of enumerate the different edges and points on the triangles. That's very useful. That's all that this does down here. And finally, we're evaluating, as it says in this step here, if it's not a constraint and it violates the linear condition, we swap that diagonal. So from V1 to V2 to P to V3. So we're evaluating the condition using the same math from the last video for the linear condition, the math being this math. So if we satisfy the condition, we're fine if it doesn't satisfy, that is if triangle zero has a vertex inside triangle one circumcircle or vice versa, we have to swap out the diagonal and that's what all this does down here. Again, it's not very difficult, just take a look and we increment swaps. And at that point, after you go back and loop through your new edge loop, you're basically done. So it's a very simple algorithm to kind of conceptually understand. I will say it's very challenging to implement and it will take you a couple of days to make sure there's no bugs in the way this works. So I'm done with that. Compile. And it should run. And now we can pass an input test case number to pick a test case. I'll just show you how this looks. So passed in test case one or test case four or test case eight, different kinds of results occur. Basically this first block of values here, this is the post-triangulation pre-constraint for text list. So basically that just shows you this sort of result. That's from video episode 16. The second set of indices here, these are the post-constraint vertices. So this will give us this triangulation. And at the bottom here is the actual point cloud itself shows you which index, you know, zero, four or seven, what that corresponds to it for each triangle. So the points in space. And so for this, I have a MATLAB or an octave test file here with all of the point cloud shown here and all the vertex lists here. In MATLAB, you have to have, for a polygon, you have to have the first, second, third, and then back to the first index. When you wanna draw a complete loop of something. So that's why there's four indices here and not three. But I'm also very sure I have a vertex list for test cases, you know, this is the original solution. This is test case one, test case two, three, four, five, six, seven, and eight. And then I'm just plotting the results. So I'll show you how this looks. We'll run it through octave test M. So here is the point cloud. That's what I'm passing into the algorithm every single time. It's the same point cloud for every test case. Hit enter. This is the triangulation. Once we've finished with sort of episode 16's content, this is how the triangles look. It's a very nice triangulation. So this is before any constraints are applied. Now, this is the first test case. So in this test case, I passed in a constrained edge that already existed as part of the triangulation. Because of that, there's no difference between test case one and the pre-constrained triangulation. Cause this is the edge I think that I implemented here as a constrained edge. It already existed. So there's no work to be done. Test case number two, I added a constrained edge that crossed one existing edge. So you can see that this edge here was the constrained edge and it basically replaced this edge. It crossed through and it replaced that edge. So very simple there. It seems to work with test case number two. Test case number three is the same thing, but a different constrained edge. So here this blue line is the constrained edge that I implemented. And you can see it basically replaced these two triangles with these two triangles. So very simple. Test case four, a more extreme example of a constrained edge. So here the edge in maroon here from this diagonal to this diagonal location, that's the constrained edge. And you can see how we've formed these triangles around that constrained edge. So that seems to work. Test case number five, another extreme example perhaps even more extreme. This maroon line here is the constrained edge. And again, we've retrangulated around that. So everything seems to harmonize triangle. Test case number six, here I'm implementing two constrained edges. So this blue one and this maroon one, they are both not in the original triangulation. So we're both implementing those two edges and they crossed through here and here. That's no problem. It seems to work with two different edges. Test case number seven, this is one of the harder ones. I'm implementing this blue line here with a constrained edge as well as the screen line. So they're kind of touching at one location and they're passing through different amounts of the body. And again, you can see it forms very nice triangulation around those two constraints. And lastly, test case number eight, passes I think from this point here to this point here to this point here. And so this is a combination of existing edges and intersecting edges. So again, it seems to work for this case as well. So for all those test cases, we seem to have a functioning triangulation with a constraint. And so I believe we've implemented everything that we need to do to have a constrained or any triangulation in 2D. Thanks for watching. Have a good day.