 Welcome back cat from scratch episode 16 topic today is Delaney triangulation We're gonna stick with 2d Plainer triangulation that's unconstrained in this video. We will expand that in the next video What is triangulation? It's going from a point cloud to a set of triangles like this and Delaney triangulation Keeps these triangles at a very good aspect ratio So you don't get these thin slivers that you would in other triangulation methods And in this video, we're gonna talk about We're gonna implement pretty much the work done by Sloan here He has a bunch of papers on this topic. You can look them all up. They're all open access And one good thing about Sloan is that he also includes the 4 train routines in the appendix So you can just check that out if you ever get stuck Question why bother with this triangulation method? Well, it's very good to mesh faces with this. So if you're doing finite elements or finite volumes This is a great way to make a very well, you know, shaped elements on those faces It's good for rendering if you recall OpenGL and other rendering engines render triangles. So you need if you need triangles you can do so with this algorithm It's good for contouring and topography. Let's say you were on a hike with your phone and you were recording sort of GPS coordinates and your altitude and you got home You could very easily use this algorithm to create sort of a terrain map of where you were path planning a similar idea and also a bunch of weird math stuff that I have no idea about Now one question is what about video seven in the series we talked about triangulation Well in that video we talked about ear clipping. That's a very efficient way to create a triangulation from a boundary of a polygon It's a very fast one problem though is that it's somewhat challenging to handle interior nodes and holes Which we have to be able to handle in a future video And you can do it with that algorithm But it becomes very complex and unwieldy and it's not worth the effort if you can just do this For a similar amount of effort It's just a better way to handle the problem and then more importantly the triangle quality when you're cutting off these convex faces convex Triangles using ear clipping you end up with a bunch of weird like slivery triangles like this That are not very good. So this algorithm gives you much better triangles and What we're going to talk about today focuses a lot on this condition Basically, you have to evaluate whether or not a given point is inside the trying the circle that Sort of circumscribes Every other triangle so Imagine you had triangle v1 v2 v3 and you drew a circle through those vertices That circle would contain point P And so this would fail that condition. You don't want any Points to follow within the circle that goes through Any other triangle vertices in your point set and also by the way v3 Falls in the triangle for you know these three points to so it feels on both both fronts However, if you swap the diagonal here from v1 v2 vertically to one From P to v3 Horizontally and you draw the circles around those new triangles You see that no point falls inside another triangle's circle So this passes a condition this fails and the way the whole algorithm works is basically you're gonna add points to the Domain and constantly keep checking if this is satisfied everywhere where things changed basically So algorithm so this is one of those things where the algorithm is super simple when you when you hear about it And then when you implement it, it's really hard a lot of small things that you don't think about come up And then after you're done, it's easy again because oh, yeah, that wasn't so bad I mean those things were so small no big deal So, you know, I've been through the whole process. I know that it's like the up and down So the top level algorithm that I'm just gonna talk about top level is basically you create a big triangle That contains your entire point set when I say big I mean humongous like a hundred times larger than the dimension of your point set Then you add in your points one by one You know to that point set and add triangles as you go and you answer a couple questions for each of those points You ask what triangle are we in? Obviously the first point is in your big triangle But after you add in 10 20 30 points and create triangles, you don't know what triangle you're in You can be in triangle 16 or in triangle 45. You have no idea So you figure that out once you figure that out you break that triangle into three That are attached to your point that you just added So imagine you found that point D that you're adding is inside triangle ABC What you do is you break triangle ABC into these triangles DAB DCA and DBC And one important thing that you know, you're probably gonna gloss over because I'm just saying it But it's very important and if you make a mistake like I did it will take you hours to debug You want to make sure that for every triangle that you're adding that D is the first Vertex in your structure because that will help you stay organized about Adjacencies to these triangles as you propagate through the model So keep in mind that you want to keep D at the front of this list Pro tip And then step C here You're gonna swap edges to make sure that Glany condition is not violated anywhere and remember you have to propagate out So for example, if you just if you're adding in D to triangle a b and c you're also checking you know Triangles adjacent to this as well and You know if you know this violates your condition Between this triangle you may have to swap here or swap here Right with this diagonal, right? You're not you have to be able to figure out which diagonal it does not does not fail that condition between you know between This diagonal and this diagonal and this diagonal and this diagonal and this or whatever this diagonal and this diagonal You have to be sure which diagonal is is better for your condition Once you finish that you're done you basically finished the entire algorithm However, you still have the giant triangle that surrounds your point set So you just get rid of that and then you're done So in more detail There are two optional steps that Sloan recommends that are Very good perform performance perspectives the first one improves your precision more than anything basically you It's optional, but you basically scale your point cloud from whatever it is to between zero zero and one one And so it's not just making sure that you know you have y min y max x min x max You know equal to zero and one that's not correct You want to also be able to preserve your relative lens along both axes because this condition uses a circle to determine, you know If you pass or fail not an ellipse So you have to be able to preserve your x and y dimensions, you know the scaling of them. So Don't map between zero and one map Just between the maximum of the two dimensions make that zero to one and let the other one just fall out Whatever it happens to be, you know, maybe it's zero to a point eight or something I don't know and map that and just you know, but zero zero to one one This will help us for the later step and also helps us get better precision Later on Step B another optional step is to sort your points by their proximity This will severely help you speed up your algorithm if you have many many points many many triangles Doesn't matter so much if you have like, you know small amounts of triangles, but you know it it still will it It's it will be noticeable at that level as well Basically what you do that's on our commands is you cut up your domain into this grid shape And you start numbering your points as opposed to having this point be zero this won't be one this might be two This won't be three instead. You basically number your points in a sequence following a continuous You know trace through that grid and it could also be you know a spiral or something You don't have to you don't have to keep, you know a weird grid shape thing like that You can do this But you want to be able to number things like zero one two three four five six That will make it so when we're looking for triangles later on that We can find them very closely To where we just did an operation Previously so that's how this works. No, it's an optional step just like this one Step C the first required step is make your big triangle. So if you map between zero and one look like I recommended This triangle should go between I don't know these coordinates. That's a good idea. That's what Sloan recommends And this will be your first triangle for your algorithm Step D is to loop through all your points that you want to add So be it's 10 points or a thousand points loop through them all and find the triangle that contains them We did this many times in previous videos I don't want to spend too much time talking about it now But you can compute the normals of your triangle of interest and then compute sort of the dot products of You know all of these vectors with those normals and if they're all the same sign basically you're inside of the triangle if not Then you're not inside the triangle one important thing though that you have to know this I didn't know this took me many hours to debug The algorithm they're talking about today it considers a point in the triangle if it's anywhere inside obviously But also anywhere on the boundary so S1 S2 S3 either all with the same sign or To with the same sign and one with zero or within a tolerance So you could you can tolerate one of these to be zero and the other two to have the same sign as well Okay Step D2. This is part of that step It's the optional triangle search step So, you know if we're just looking through triangles in a loop like zero one two three four five And you have a thousand triangles That might take forever because if you're on the thousandth triangle And you have to loop through all those ones and constantly evaluate this that's extremely expensive of an operation so one clever way that they mentioned that we can do this is to basically compute the Dot product of basically say you're looking for point P and you guess triangle, you know Like this triangle right here for example, what you do is you take the dot product of all the edges Outward normals, you know and the vector from this triangle to point P if those are you know positive You go in that direction So you basically go from this triangle to this triangle and you can basically keep repeating that process until you get to where you are So you'll always be moving towards the correct solution Not just guessing randomly, you know, you might just guess when it ran random triangles over and over until you get there That's a very slow method What's loaner what's loan recommends is obviously much much more efficient Step e this one is pretty much the you know the heavy lifting of the actual algorithm is checking that condition You can read the paper. He has a lot of you know equations in there He basically boils it down to evaluating the condition in these you know for four and a half steps This handles a bunch of the numerical problems that you may have when you have very very wide or shallow angles for your Triangles you may have the triangle like this This this sort of sorts those out very well. So we're going to implement these equations here. So x1 x3 x2 those basically refer to The coordinates of points of Adjacent triangles I'll talk about that in a second But just keep in mind that there are Equations to solve this and you can look at Sloan's paper to see how they work. So that's how you check the condition Basically the condition will tell you whether or not you have to swap Diagnose or not because if you recall the whole point of the algorithm is to pick between this diagonal here and This diagonal here So basically swap means you're currently bad try to be good and no swap means you're currently good Don't worry about it. So that's what this swap indicates Now where the condition fails You have to swap Diagnose and this is where a lot of people recommend having these fancy data structures about adjacencies and stuff Personally, I think you can just do it with with simple arrays Have one array that carries which vertices are part of which triangle and then have one array that Indicates what triangles are adjacent to your current to each triangle that you're talking about. So That's that's what you have to basically have for structures to make this work So basically if you want to add point P and you find out that For the triangle that you've determined that you're inside the adjacent triangle has these points v1 v2 and v3 You can you basically have to check If you have a swap that indicates that you have to basically move from this red diagonal to the screen diagonal and You have to update obviously all these triangles. So triangle L is a triangle that you just created Triangle R is one that was opposite it that was already in the model and Then triangles A, B and C are triangles that are adjacent to these that that are going to change adjacencies if You change from red to green So you have to be able to not just keep track of the vertex ordering for these, you know Triangles L and triangles R that you're dealing with directly But you also have to keep track of the adjacencies of triangles A, B and C that are just you know going along for the ride and Then here's a very important thing step G on Triangles A, B and C you very well by doing what you did, you know here. You may have broken delaying this of this Triangle A, B and C for example, so you have to basically keep swapping diagonals and repeat steps F and G until your stack of triangles that you're adding is Empty and Sloan recommends actually using a you know a stack for this and we will do that in our implementation but There may be other ways to handle this as well Then step H as we talked about before is to just straight up delete your big triangle And then step I is to map back from your 0, 0 to 1, 1 To your original domain. So that's the that's the full algorithm Described a very very top-level Perspective so now I'm going to talk about the actual code and so let's just Check the header file first. So only two functions required for this I think Sloan puts like a whole bunch of fortune functions, but you can do them all in two functions here one of them is Well, what only one that you see here in the header file this Triangulate function basically pass it instead of points And this is your point set basically and this is your number of points in the set So we're talking in two dimensions here. So only only two floats in in that dimension now I'll open up the the main Main dot see just so you can see the test cases I'm going to run I have two test cases here. We're passing in a set of nine points like this in in in 2d We're passing them into the function like this and test case number two passing in seven points like this The only difference is test case one is a very simple one just to make sure everything works Test case two is a pretty much It's a more challenging because it has a bunch of collinear points Including points that are collinear with the super triangle the very big triangle So this if this works everything will work So we'll chest both both those out So let's open up the actual Implementation here here. You can see the actual two functions that that we are talking about the first function is This Boolean function, which is called planar point within triangle. This basically evaluates Or let's see it all the way to this condition here if a point P is inside, you know a 1 b 1 c 1 Very simple. We had this many times before however, one key difference is that The condition for a point being in the triangle or not it can all be negative, right? So we have a counterclockwise winding for each triangle or you may have to them be negative and One of those values be below a tolerance. So that basically means that You can have a point on the edge of the triangle not just inside Both of them are valid conditions for a point to be considered inside the triangle that you're talking about That's this function now into the actual function for the triangulation Here it is. So the first thing here is We have to find the min and max boundaries at the pull off the point cloud so that correlates with Finding these values x min x max y min y max very simple You're just looping through then we we remap everything and we preserve that aspect ratio to between zero zero and one one So we just map this find largest dimension in with this way and scale by that to get everything between those two those two points That's basically doing this operation here Next up here is step B just to use that bin sort that Sloan recommends So we implement that right in this way basically we calculate p and q p and q are the IDs of the row in the column that we're in in this Domain here then we use insertion sort to figure out Basically which points belong in which sequence in terms of the the bin order like this So very simple. You can read the code to see what's going on if you'd like Next step here Step C is to make the big triangle. So simply you put right here. We're basically we're we have this let me talk about the actual Arrays we have here. So the first story I'm talking about right now is this points array. So points is basically a Array of all the data points that we have it's actually passed in To this to this function. So it has the x and y coordinates of everything in the in the structure And so we basically you fill the first three vertices with negative 100 negative 100 etc and Increase the number of points in our structure here by three So at this point if we passed in nine points now, we'd have 12 points, right? The first nine being the points of the actual cloud the last three being points of our super triangle Now we have to make these data structures here So I have one structure for vertices and one structure for the triangle adjacencies basically the vertices array just encodes which vertices are part of which triangle and the last Second right here triangles basically encodes which triangles are adjacent to the current one Counting in the same direction as you would for the vertices So the first triangle we have here is triangle, you know, I guess Ten eleven twelve or whatever it would be based off however many points you're passing in This is the super triangle and you see negative one negative one negative one negative one basically means that There is no Adjacent triangle so this basically means that you're on the boundary of space Obviously once you're sort of in the point set you won't have negative ones You'll have actual numbers here, but when you have no other triangles out there. The answer is just negative one so next step here is to Loop through points and you know add them one by one and check which triangle that they're in so This year this loop it into II iterates through every point that's not in the super triangle and then it This loop here this while it basically looks for What triangle are we in and the first guess for a triangle that we're in here is the last triangle that we created? This is also a technique that someone recommends for efficiency is to always check if the previous triangle Contains the point that you want to add now because because we added the points in in this way That is actually pretty likely that the triangle that you just added Includes the next point So we start with that off at the previous triangle And then we have this function here. I won't get into the function yet I just want to talk about this actual while loop. So this while loop Is is what we are using to Continue until we found the triangle that we're actually inside so to find that we're trying they're actually inside We have to go another Condition here. Here's where that other function comes in that we just created planar point within triangle Here we're passing in the point of interest that we're trying to add and the vertices of the triangle that we're guessing that that the point is inside so That basically runs, you know this code from before that I just talked about No problem there At that point it does a couple things so the first thing it does is it's Let's see. Where is it? here it um It deletes Triangle ABC and it adds in these triangles DAB DCA and DBC. That's the first thing it does. That's right here All of this and it creates the adjacencies as they should be and then also It updates adjacencies of the triangles that are nearby. So Recall there are triangles that are Adjacent to this you have to also be able to update the adjacencies for because you're adding new triangles here So this is a new triangle. This is new triangle. This is a new triangle Well, I'm necessarily because you may be keeping one of them the same as the old one just changing the nodes on that triangle But anyway, long story short, you just have to check the adjacencies of all these triangles to make sure they still make sense That's something that you have to do on your own. I can't talk about that. I'll be here for 10 hours You have to just go through and make sure deliberately and you know Very closely pay attention to make sure that every triangle added in has the exactly correct adjacencies and come up with The algorithms you have to have to make that work properly So that's all fine at that point This stack comes into play. So if I scroll on to the stack part basically This this sort of step f and step g are the stack relevant steps So basically the way it works is if you have a triangle with an adjacency That you're affecting By your swap from this vertical diagonal to this horizontal this horde is on the diagonal You have to be able to check if that adjacency is messed up in any way In terms of its passing of the condition of This condition here So you basically add triangles to the stack that have Potential adjacency problem for this condition and you loop through and check all of those And then they probably get out quite a bit You may propagate from this point to this point to this triangle to this triangle as you go Further and further away from you know that point you have to constantly keep checking until no more points have a Have an issue So that's what this loop through stack does I won't get into it except for the part that's right here about talking about the circum-circle thing Here is where you're again checking the condition It's called step e Here's where you're checking if that actually works. So in reality step e f and g are all sort of steps that have to occur on Triangles that are on the stack So you have to do that and then Lastly once you're out of that Entire loop you basically have to Renumber your triangles in the proper order and you have to get rid of Anything to do with the big triangle vertices because again, you don't care about these points here All you care about is the points inside. So Get rid of all these points and also these points will have triangles off of them. Remember, you know a bunch of triangles will be touching these points. So you have to be able to Get rid of all those triangles and then remember the triangle is based off having removed all these triangles So you do that and then at the very end you undo your mapping and finally You plot everything out. So that's how the code works So compile and Run. Here's the results. Don't worry about these actual values. Just know that it works As implemented and I will show you a visual of how this works using octaves. So The first set of data points here is shown this is test case number one You can see these data points all around here and nothing special looking data points They're between negative six and six in X and between negative four and seven and Y or something and there's nothing that's collinear about them They're just a generic set of points when you plug in these triangle vertices and this point cloud values into, you know Functioning out of that plots everything and you you run it you'll get this triangulation And if you look any triangle that's drawn here and if you overdraw circle through the vertices It would contain no other Triangles points. So this passes that Delaney condition the next test was this test case to here Now this one it had a couple points of interest So if you see, you know, it has she has less points only has seven actual vertices versus nine in the previous test However, this vertex here this vertex here this vertex here are all collinear not only that But when you add in the super triangle that super triangle point is further along this line way way up So there's actually four collinear points and that's like a really hard Condition to solve with this algorithm. So basically if this works here, it'll work with any sort of test case If we add the triangulation in it looks like this and again, if you draw a circle through any and all of these Triangle vertices it will contain no other triangles vertex and that's the condition for the algorithm that we just used So at the end of the day long story short half hour video We came up with a way to implement the Delaney triangulation algorithm In a very sort of efficient way. We didn't have to use a bunch of weird data structures We did it with just simple arrays using malloc and realloc Repeatedly and the algorithm is not that complex. It's but it's a there's a lot of small little details here and there That if you don't get perfectly right It won't work perfectly right I took many many hours over the past couple days debugging these small things like a single triangle was off And I didn't know why anyway Just if you want to implement this just take your time and go through step by step by step making sure every single triangle that you're adding is Exactly what it's supposed to be And then you have no problem. So anyway, thanks for watching. I hope you have a great day