Heat Methods in Geometry Processing

thank you all for being here so I'm going to talk about a algorithms called heat methods which are for practical geometry processing and actually these are going to be related to a lot of things we've already heard about in the workshop PDE based signal processing computing geodesic distance even statistics on manifolds and lots of other things the central idea is to use short time diffusion as a building block for geometric computation and the reason for this is that diffusion processes capture a lot of information a lot of geometric information about the manifold and work I'll talk about is represented by these two papers one with Clarisse title and maxford etske who's here at the workshop and another one by my PhD student nicolas sharp and an undergrad at CMU use of sullivan who's now at Caltech what we want to do is build up all sorts of basic geometry processing operations computing distance doing parallel transport and so on based on this principle of heat diffusion and actually although some of these seem like really basic things you'd imagine there are already algorithms to do that something like compute a karcher mean of a distribution okay that's something somebody surely has figured out actually we'll see that in many of these cases these are really the first practical algorithms for doing this on general surfaces and so I'll talk about that since publishing the original heat method paper there have been a lot of real-world applications of this technique for things ranging from an atomic data visualization to various kinds of computational fabrication doing optimal transport on curved domains there you soon picked at Pixar and animation software and you can also find it in the computational geometry algorithms library Siegel so why is this principle of short time diffusion so useful for computation I think there are really kind of two big reasons one is generality so to formulate these methods we really need to know very little about your manifold your discrete manifold just a inner product or mass matrix and some kind of first order differential operator a gradient operator or covariant derivative and so what that means is we can really apply this any kind of geometric data structure polygon meshes point clouds voxels ations we're really not just talking about one specific kind the other reason is efficiency so the computation is ultimately going to boil down to just solving some standard linear partial differential equations heat equation Poisson equation so things that are extremely well studied in scientific computing and for this reason any algorithm that's formulated in terms of these these heat methods can immediately sort of piggyback on fast implementations parallel solvers and so forth right if you build your algorithm this way you're not responsible for having to optimize all that code I think another reason that often gets overlooked is how hard is it to implement so especially for people who are already working in this setting where they're solving linear PDE s it's kind of a no-brainer to implement these operations this way you already have your solver you already have your Laplace operator or whatever why not compute things this way and I think that's why so many people have used it for for real applications ok so I'm gonna start with the original heat method which really nicely explains the basic principles and the motivating problem here is computing geodesic distance so imagine I have some curved manifold and I want to know all the distances to one given point and just to be really clear here what I'm talking about is not the straight line Euclidean distance through space which would be very easy to compute but rather the length of the shortest path along the domain or the geodesic distance in fact for everything I'm going to say I never even really need to know where this manifold isn't sitting in space all these algorithms are purely intrinsic ok so how can we get our hands on the geodesic distance well one thing is we could solve the geodesic equation which says a curve gamma to be geodesic well it's tangent vector shouldn't turn as we walk along the curve and this is a non linear ordinary differential equation in classical computational geometry the approach to solving this is to really explicitly track these geodesics as they expand out from the source and this can give you really really accurate solutions on polyhedral meshes but also comes at a cost very expensive computation and takes a whole lot of work to make this more efficient another approach is to instead solve the hamilton-jacobi equation associated with the geodesic equations so the so-called eye Conal equation which says that if he is the distance function describing the distance to the source point or source set then the gradient has unit magnitude so intuitively this is just saying the change in distance the rate of change in distance is 1 meter per meter and so if I have a source set let's say a doesn't have to be a point I could have some other set I just set feet to 0 along this set and I compute the viscosity solution to the I Colonel equation everywhere else that gives me this this distance function this equation is nonlinear and it's hyperbolic and algorithms for solving this equation generally work by propagating outward from the source so the kind of canonical example for curved surfaces is the fast marching method of Kimmelman Safian right and so all of these algorithms in one way or another kind of have the spirit of Dijkstra's algorithm you start at some source point to compute distances you march outward in some deterministic order and so they all share a common set of challenges so for instance if I now want to compute the distance to another source point I kind of have to start all over from scratch there's not really a whole lot of computation I can reuse from that first computation another challenge is because there's kind of a serial dependence on these these updates it's really hard to parallel eyes these methods if you have a multi-core system you know GPU you have a cluster you have to think really hard about how to make these things parallel and finally a lot of these traditional algorithms really start with a triangle mesh as or some kind of simplicial mesh as the starting hypothesis and formulate the the algorithm from that perspective which means it's not immediately obvious how to always generalize these things to point clouds or polygon meshes or other geometric data structures ok so because of these challenges let's go and take a look at a different equation the diffusion equation which simply says if I have some distribution of heat bumps in this distribution or getting it smoothed out over time and unlike the eye kernel equation which is nonlinear and hyperbolic this is a linear parabolic equation computationally this means you can do something like well just repeated local averaging with your neighbors to try to get a solution and in general there are lots of different numerical strategies of course for solving this kind of equation going back at least two piron's an important special case is when my initial heat distribution is just a spike so Dirac Delta centered at some point X and as this distribution spreads out on a Romanian manifold it looks really something like a Gaussian and this Gaussian like function will call this the heat kernel so the heat kernel just says how much heat is diffused from the source point X to any other point Y after time T okay okay so this is all very nice some nice linear equation but what does this have to do with distance why is this helpful well there's this very beautiful formula I srinivasan farad on which says actually you know I can really think okay if this thing really looks like a Gaussian then what I can do is I can kind of invert that Gaussian and as time goes to zero as I let the heat flow less and less and less and less actually this will really give me exactly the geodesic distance all right I'll really approach there's some straightforward intuition for this which is okay if I think about diffusion as a bunch of particles taking random walks well if I find a particle very far from the source after only a very short amount of time the only thing it could have done is traveled on an extremely straight path right or a geodesic path so that kind of gives some intuition for why this should be true so from here you think okay well now computing geodesic distance doesn't sound so hard I should just be able to compute a heat kernel solve some linear equation plug it into this formula and I'm done all right unfortunately it's not quite so easy because of approximation error I'm always solving this for some finite time T and so if I really just go ahead and plug in Vera Don's formula I get something that doesn't quite look like distance so on this spiral I'm thinking about the distance going from the center and then intrinsically outward along the spiral it really should have evenly spaced level sets right not these non uniformly spaced level sets we see above and the reason is well I have some error in my approximation of the heat kernel if I now try to just pretend it's a Gaussian and invert it well if it doesn't look like a Gaussian I don't get something that looks like distance okay so we need to think a little bit deeper about this problem rather than just blindly applying this right on and so what we can do is actually go back to the I Connell equation this this provides a really useful insight which is oh well the gradient should have unit norm at every point okay so as long as I can construct some function that's monotonic in the geodesic distance then I can take the gradient of that function and just normalize it to get something that looks like the gradient of distance okay doesn't have to be the heat kernel but just something something that's monotonic in distance and from there if I want to recover the distance itself I'm just looking for the scalar potential whose gradient agrees with this normalized vector field okay and that amounts to solving a Poisson equation so all together this gives us a recipe for computing distance that's exactly the heat method so what do we do again we solve the heat equation for short time T we compute its normalized gradient X and then we solve a Poisson equation for the distance V and the beautiful thing about this picture is we've taken this problem that is traditionally nonlinear whether you think about the geodesic equation or the icono equation and we split it up into two easy linear pieces by just doing a trivial nonlinear operation in the middle right we've sort of hidden all the non-linearity in this middle step and one thing that was really surprising to us when we tried this okay it sounds good in principle when you go and try it in practice it works surprisingly well just a simple scheme so here we're seeing the distance to a set of points on these models on the left or as I said it doesn't have to be distance or points you could also say to a curve so maybe distance to the boundary curve of course to generate results like this I have to discretize right I have to turn this into a real computational algorithm and well really the algorithm we've described is an algorithm in the smooth setting it's involving you know PDEs and this kind of thing vector fields and so there's many many different ways we could try to discretize it in terms of time discretization we do something very simple we just say let's take one step of backward Euler and what that does is it transforms our parabolic equation into an even easier linear elliptic equation right identity minus T times the laplacian applied to these so we equals the initial financial values maybe something like this spike of heat now one thing you have to ask is what do i set t - I'm taking a finite step of backward Euler how big should I make that step very Don's formula would seem to imply that if I make that time step smaller and smaller and smaller to get more and more accurate of course it's hard to believe that that's that's really going to be true and of course and what happens is there's discretization error so as we make T smaller if it gets too small we start to recover something that looks like the graph distance rather than the geodesic distance you can really prove that this happens so based on just some basic scaling arguments about the heat equation we can say really we should be setting the time step to the square of the length scale so something like the average edge length in the mesh and in terms of the constant of proportionality just empirically you go out you try this on a huge number of examples and you see I yet you know really that constant of proportionality that should just be one to get good accuracy okay so that's it for time discretization in terms of spatial discretization as I mention in the beginning we don't need to know much about the domain we just need some kind of discreet gradient operator and an inner product or mass matrix because if we have this information well then the divergence is the adjoint of the gradients the laplacian is the divergence of the gradient one thing you could do from here is say okay I know how to compute all that on let's say piecewise linear finite elements on triangle meshes that's one common example I'll talk about a lot but you don't have to do it that way all right you can easily get this information the gradient operator the inner product on a point cloud on some general polygonal mesh on a voxel ization all sorts of other things and I think that's that's again why this is a really nice way of approaching this this distance problem because if you go back to these traditional computational geometry approaches you know it's really not a not always clear how to apply them to something like a point cloud if you do follow this piecewise linear path okay you can ask about accuracy right of course just like any finite element method the quality of your triangulation will impact the accuracy of the solution but when you're solving elliptic equations life is somehow a little bit easier than when you're solving hyperbolic equations so when you solve a hyperbolic equation what kind of triangulation do you need well let's think about something like this very obtuse triangulation and a source that's gonna leave from a it's gonna reach B first and then C and D but because of the way things are connected up C and D you're going to get connected in the wrong order so I violate causality and it's so to really be kind of safe when I'm doing hyperbolic problems I'd really like to have a non obtuse triangulation unfortunately not obtuse and acute triangulations are very hard to come by even in the plane took a long time for people to come up with good algorithms for surface meshes you know cute regulations I'm not really sure what you can even do and in higher dimensions there's situations where you can't even construct an acute triangulation so this is a very strong condition for elliptic problems the thing you care about is a maximum principle so if you have a harmonic function meaning something that's in the kernel of your discrete Laplace operator what you what you would like to be true what you should know it know should be true from the smooth setting is the value at every interior node should be a convex combination of the neighboring values and this can be ensured by just making sure that your mesh is de l'année it satisfies the de l'année property this is very very easy to obtain constructing de l'année meshes you can do it very easily and quickly in any dimension although I should say it's not even necessarily necessary so it's a it's a sufficient condition if you have a de l'année mesh then you're guaranteed to get a maximum principle you can often get this maximum principle even if your mesh is quite bad so here's some examples we're running this heat method on meshes that are very far from de l'année and still recover quite a good approximation of distance now if you want a little bit more peace of mind you could say okay how do I get a de l'année surface mesh for instance this might sound like a pain but actually there's a very beautiful perspective on this this problem which is that you can always achieve this still any criterion for any surface mesh without changing the geometry no change to the geometry at all this may be sounds a little crazy until you see this picture so the idea is you use something called an intrinsic triangulation so normally when you work with triangulations you think okay if I'm going to do something like let's say I want to flip one of those edges well I get a picture like this I can connect the edge along a different direction between two different vertices and I've really distorted the geometry right the surface is no longer where it used to be well something I can do instead is an intrinsic edge flip I can just imagine that the geometry stays where it is and now I have an edge that goes across the original geometry intrinsically this is still straight intrinsically I still just have three Euclidean lengths which just describe a triangle I can do all my build my operators in the same way right and in fact this lets you get really nice results even on very very poor quality meshes so here's a mesh that was designed not for finite element simulation but for 3d printing so just a completely horrific mesh really really nasty angles and if I try to run the heat method directly on this mesh guess what it doesn't work out so well on the other hand if I run a simple edge flipping algorithm which of these intrinsic edge flips then I can get much higher quality results without changing the geometry without changing the vertex count and without doing much computation about the same amount of computation I need to do for matrix pre factorization so that's really nice if I want to get even better then I can do an intrinsic version of de l'année refinement so not only flipping edges but also inserting vertices and now I get very very close to the reference solution without increasing the number of degrees of freedom very much so long story short there's a lot of good ways to get good meshes when you need the de l'année criterion right it's pretty easy to come by in terms of order of convergence you can check that empirically the heat method seems to at least in the strangulated setting exhibit linear convergence so same as fast marching on triangle mesh is an important thing to keep in mind is that actually you can't possibly do better than second order accuracy and the reason is you've already introduced discretization error just by triangulating a smooth surface and so if I take a smooth sphere and I triangulate it and then I compute the exact distance along the polyhedral surface this isn't only a second order approximation of the true distance on the sphere of the arc length okay so if you're doing linear you're doing pretty well although you'd like to do a little better one thing that we can do that's that's quite nice is okay we don't have to use triangle meshes right I said you can use any discretization you like so how about running the heat method on higher-order surface representations so people have worked this out or c1 finite elements or situ subdivision elements and again because I can just build a Laplace operator or gradient operator on these things and apply the same method you can get better accuracy so closer to second order rather than first order different way you can improve accuracy this method is to say well let's go back again to this icon elision and I have kind of a variational principle but if I just take the l2 residual of the I Conal equation then it has an Euler Lagrange equation that looks like this Laplace of the distance is the divergence of the normalized distance gradient well that sounds very familiar that sounds like the final step of the heat equation are sorry the heat method and so you one way to view the heat method is that it's actually the first step in in fixed point iteration for this distance but using an extremely intelligent initial guess right that's that's key so if you want to increase the accuracy a little bit you can do a couple more of these fixed point iterations using the same pre factored matrices no change to the matrix and you see that the accuracy gets a little bit better but actually the the initial gets already gets you most of the way there in terms of performance one really nice thing about working with linear equations is you can pre factor the Associated matrices so if I want to compute the distance to many many different points on the same domain I can pre factor the matrices for my heat equation and my Poisson equation just once and then iteratively or just repeatedly apply back substitution back substitution is some tiny tiny fraction of the cost of solving the original system so I get about an order of magnitude speed-up by doing this just to give a sense of performance compared to kind of traditional methods if you're using these computational geometry methods that are very exact but but slow this could take a matter of hours to compute this all pairs distances using something like fast marching which has the serial ordering takes maybe around 55 minutes on this example and heat method takes around five baths it in five minutes so this this really does give you pretty big amortized gains on the other hand factorization can be a problem if you're working with very very big problems all right you start to run out of memory very quickly so in this case you can switch over this is not our work but work by others who say let's switch over and replace this factorization with an 80mm scheme why not and you can show that this exhibits really good scaling roughly linear scaling all the way up to 200 million triangles just on our ordinary CPU they're very GPU friendly algorithm you can imagine extending this to other parallel architectures okay and another nice thing about the method is pretty much any kind of platform you're on you can find some free implementation of the method including one that is just fun is you can do this in the web browser so here's you can go to our page with lots of nice little geometry algorithms on it including the heat method and you can see here's our favorite shape the bunny so fairly really good number of triangles right in the web browser you can compute distance all right so that's it for the that's it for the the ordinary heat method and so far we've been able to compute geodesic distance so what can we do about computing other kinds of quantities yes yeah no I won't get that good yeah so so the way that your source set shows up is just on the right hand side so you basically you have some you imagine you your curve is some kind of house Dorf measure and you discretize this onto whatever your finite element space is that works pretty well you can get what's up right true yeah so so you can get kind of more accuracy by actually changing the mesh but comes at a cost so it's really just a speed accuracy trade-off yeah right in the limit they will be the minimum of the function that you compute numerically they may not all be 0 so that's what I mean by you kind of have this speed performance trade-off or speed accuracy trade-off we can talk we can talk later ok so so given given the set up a natural question to ask is what what might happen if I now diffuse vector data rather than scalar data and the short answer just kind of the the preview of what happens is well we obtain parallel transport rather than geodesic distance more precisely what happens is if I start with a vector at a single point and I do this kind of short time diffusion of vector vector data apply similar kind of method what I get is a vector field where at each points I have the vector obtained by parallel transport along the shortest geodesic to that point ok so this field is the field obtained by parallel transport of that initial vector to every other point of the surface along shortest geodesics initially this sounds like kind of a strange quantity why do I want to compute this never heard of this before why is it why is it useful well we'll see that it turns out actually to be remarkably useful for all sorts of standard geometric queries it's kind of a basic geometric kernel that lets you do all sorts of other things ok so what is vector diffusion okay here's a picture of scalar diffusion alright again I'm starting with this spike of heat and letting it diffuse out of our time vector diffusion looks something like this I start out with a vector concentrated at a single point and I diffuse it over the manifold now in the Euclidean setting I could just do this by diffusing the individual scalar components and I would get a vector field that is constant modulated by a Gaussian right looks something like this on a curved manifold I have to think a little harder because these tangent vectors are going from one tangent space to another so exactly how do they how do they move as I go from one point to the other so this information is essentially encoded by what's called the connection laplacian there are many different ways to understand the connection laplacian it's not so different from the ordinary laplacian okay it's a second-order elliptic operator it's acting now on tangent vector field rather than scalar functions you can view it as well as we just discussed the infinitesimal generator of vector diffusion but that's a little circular you can you can think of it as the composition of the covariant derivative with the adjoint of the covariant derivative you can write it as the trace of the second covariant derivative or the way i like to think about it is is a Hessian of the vector dershlit energy so I have an energy that just measures basically how parallel is this vector field and take the Hessian of that energy in physics this is also known as the magnetic Schrodinger operator so it describes the behavior of a particle in a magnetic field and just as a little side note this is not the same as the Hodge laplacian although there is a close relationship to the Hodge laplacian via what's called the veidt's and BOC identity so it's related by the curvature of the manifold this is the two dimensional case okay for computation what we really want understand is what's the what is the heat kernel of this operator look like what does that tell us what geometric information does it in code so the scalar kernel we can write as this asymptotic expansion we can write it as a polynomial in the diffusion time T and these functions fee are called the heat kernel coefficients since we really only care about small T we really only care about the lowest order coefficient which in this case is just the identity essentially just one so it doesn't doesn't really matter to us very much and what does matter is this leading coefficient which is the Gaussian of the geodesic distance this is what we took advantage of in very Dawn's formula to get the geodesic distance and as T goes to zero we recover this kind of Gaussian function the vector u curl looks just the same except that these heat kernel coefficients of course like this difference and again the most important one is the first one which turns out to be the parallel transport map so it's the map that takes a vector at a point X to a point Y by parallel transport along shortest geodesics really beautiful okay so now if we let time go zero will go towards zero everything else falls away and we really just have these the identity map the parallel transport map and this Gaussian of distance so if we want to recover parallel transport a natural seeming thing to do would be just take the quotient if I take the quotient I eliminate this Gaussian term I just get the parallel transport map actually I should say this is a little bit of a simplified picture there's another there's another term in here when x and y are very far apart that shows up that also cancels out so we kind of kill everything we don't care about and just get the parallel transport map okay so this gives us kind of a first attempt at an algorithm just like we we tried in the in the scalar case to just apply very Don's formula now let's just try taking this quotient so candidate algorithm is computer short time scale or you can all compute a short time vector heat kernel take their quotient the problem is again you have approximation error so here's one example where you can see what goes wrong if I have my sources not just a single vector but I have two vectors what should happen is I should get a vector field where the magnitude of the vector field at each point is the magnitude of the closer of the two vectors so it should look something like the thing on top a piecewise constant function what happens instead is I get some smearing out of magnitude blurring together of these vectors because this isn't quite quite exact the colors are just the magnitudes so I don't know this this color is the length of this vector this color is the length of that vector okay so we have to think a little a little bit harder about this so let's say for the moment that all I really care about and forget about the the vector part of this I just want to do a closest point interpolation of magnitude so I have two points they each have a magnitude associated with them and at every other point I want to function that's equal to the magnitude of the closer of the two points in the Euclidean setting what I could do is I could think of these as two deltas with different magnitudes I could convey over the Gaussian to get a function U and then I could normalize by a sum of Gaussian centered at those same points but not scaled by the magnitudes just the gaussians yeah and so what will happen is as T goes to Z I get exactly this kind of closest point approximation or extrapolation that I wanted alright so if I'm if I'm very close to you know this point all of the stuff that's going on at the other point falls away and I'm really just dividing a Gaussian scaled by the magnitude by a Gaussian and I get the magnitude right and I do that as T goes to zero I get something that starts to look more and more like kind of a Voronoi diagram on a curved surface I can replace convolution with a Gaussian by heat diffusion short time heat diffusion and I get something that looks like this so if I have points on my manifold then I get something again looks kind of like a geodesic Voronoi diagram or if I could have a curve with a varying value along it and I'll get the closest the value of the closest point along that curve okay and actually we already have all the ingredients now to formulate this vector heat method so what am I going to do I'm gonna take my source vectors I'll diffuse them a little bit with this vector diffusion and now I have at least I have some information about what's the right direction for this vector field I then scale or I then diffuse the magnitudes of the vectors so forget about the direction I just do a scalar heat equation diffuse the magnitudes and I diffuse the indicator function where these points are without the magnitudes just as we saw in the previous slide if I now take the ratio of these two scalar functions I get this magnitude of the closest vector and to get my final vector field I say okay I take the diffuse vector and I give it the magnitude that I computed in these in these intermediate steps and it just scaled them by the magnitude of the closest vector again something very beautiful happens we took this whole process of finding shortest geodesics and doing parallel transport along these geodesics which is very nonlinear thing and we've written it down as just linear PDE s and some very simple non-linearity just normalization taking a quotient at each point here's what looks like on curved surfaces so I get a vectors of different magnitudes the vector field I get is again the vector field of the closest vector parallel transported along shortest geodesics I can do that for points I can do that for curves yes there's still uh there's still a cut locus of the I think that you kind of imagine that the Mobius strip gets cut along the cut locus and now it's an orientable surface with boundary I think that should work okay right so again we can we can compute this quantity without explicitly finding geodesics without actually explicitly doing parallel transport why is this useful computationally well here's a little experiment so I took the same problem they said well what what would be the cost if I simply I didn't I didn't even go through the whole process of doing the trip the parallel transport I just went ahead and I traced out all the shortest paths that I would need to do parallel transport along if I were going to do it in this very explicit way that's already 15 times slower than just solving the the linear equations so forgetting any of cost of actually computing the geodesic distance reporting the parallel transport even just tracing these pads is already way slower than just solving the Pease well you need the parallel transport you need to know how the vector gets transport along the shortest path it's not well we'll pull a trick like this okay I think I see what you're saying about you you will pull the trick like this in a minute and maybe this will get your question okay so as far as discritization everything's just like the ordinary heat method except we now need a dis critias ation of this connection laplacian which we didn't have before the surprisingly little work on discretizing the connection laplacian few references here but there is a pretty simple approach that will work on any kind of discrete manifold at least give you some connectional question so suppose you already know how to define a discreet laplacian on your domain in your point cloud at your polygon mesh or whatever all I'm going to do is measure the smallest rotation between neighboring tangent spaces and use those smallest rotations to adjust the off-diagonal entries of my scalar laplacian so for instance for a triangulated surface I can encode these rotations as unit complex numbers and so the what would normally be so the local stiffness matrix would normally be a cement real symmetric matrix I'm going to take these little rotations and I'm gonna put them on the off diagonals or IJ an R IJ bar for the IJ entries now I get a complex permission matrix one way to motivate this is to say well this is really just the Hessian of a very simple discrete dear so energy that measures how smooth the field is as I take a vector from I to J by rotating it how much does that differ from the vector that I find at J the Hessian of that energy is one of these connection options so again we this nice property that if you have some nice already discrete laplace operator well you can go ahead and implement the vector heat method on that domain so here we're doing this on point clouds on polygonal meshes on box locations whatever we can also easily generalize this idea to other vector bundles to other kinds of vector valued data so symmetric direction fields cross fields differential forms and so forth there's not there's really not much more work to do here but I'll defer to the paper for the details the accuracy and convergence behavior again very similar to the the scalar case so even though we're doing this transport by kind of diffusion you'd kind of think oh we're sort of smearing things out you really do get some very accurate notion of parallel transport only only different from the exact parallel transport by it some small fraction of an angle and convergence again at least in the piecewise linear setting is again looks linear empirically we can go back to this question of de l'année triangulations so we can say again okay if I have a very poor quality triangulation but I still want to run this algorithm what do I do well we can extend this classic de l'année flipping algorithm just a very simple algorithm for getting a double triangulation to also keep track of tangent data so now rather than our skeletal laplacian having a maximum principal our connection with caution has a maximum like principle so if you have a bad match something that can happen let's say have a vector field that's in the kernel of this connection laplacian but have a bad triangulation it doesn't look very smooth the vector can get flipped around and so now I get instead on this intrinsic triangulation I get something that looks like this all the vectors are pointing in the right direction to be a little more precise what happens is at every at every point if I parallel transport my neighbors to where I'm sitting I will be in the convex cone of my neighboring vectors I was a lot more to say here which which we talked a lot about in this this paper will present later this this summer okay in terms of applications you know what do you do with this what is it useful to be able to efficiently compute parallel transport on curved surfaces well one basic thing is you're doing some kind of level set simulation and fluid simulation or whatever and you want to extrapolate the velocity at the interface to the rest of the domain here it's really nice because you get a global extrapolation at very little cost your back substitution rather than just giving you the velocity in a narrow band gives it to you over the whole surface and one thing that's really interesting is because of this closest point property because the vector at every point is the parallel transport to the closest point on the interface you actually nicely preserve signed distance functions without doing any kind of read distancing so just naturally the field you're advocating and will try to preserve this signed distance property so after 200 iterations the signed distance function is is still looking pretty much also as usual oops what I do here works directly on curved domains point clouds all these contexts where traditional algorithms might not be applicable another thing we can do is compute a what's called logarithmic map so if the exponential map is the thing where if I'm sitting at a point and I have a tangent vector it's going to tell me oh if you walk along that direction here's the point you're going to end up at the log map is the inverse of that if I sitting at a point let's say P and I have a point Q that I'm interested in what direction should I walk and how far should I walk to get to Q ok this turns out to be a useful query for geometry processing algorithms we'll see that in a second and it's very easy to compute via the vector heat method so what we do is we start with just any unit vector at X we parallel transport it to the whole surface and now I have a vector field that kind of represents the horizontal direction in a local coordinate system next I take a radial vector field centered around my point I also extend that to the whole surface and now I have something that looks like you can imagine this is like the the gradient of the distance function I compute the angle between these two alright the angle between the horizontal direction in the radial direction gives me theta or I guess Phi in my polar coordinates I can recover the radial component the distance by just solving the same Poisson equation that I had in the ordinary heat method and altogether that gives me R cos Phi R sine Phi okay so again solving a lot of linear equations doing a little bit of non linearity at each point at the end here's what this looks like so if I have a point P on this surface then I'm plotting the level sets of the radius function and the angle P over the whole surface what you notice is really important is I don't just get an accurate log map right in the vicinity of that point I really get an accurate logarithm back map over the whole surface all the way up to the cut locus and this is gonna let us do some interesting geometry processing now just as a point of comparison there are previous methods that compute approximations of the log map but they do this using this kind of local Dijkstra like traversals he start at the center you start marching out and one of the one of the issues with hyperbolic processes hadn't really discussed before is well if you make some error along the way that error is just going to continue to propagate as you go further and further along the surface so even though these things are fairly accurate right at the beginning you know they all look similar to the reference solution by the time you get out to the arms they can be very poor approximation so whereas we get something that's very very close to the exact solution okay why is this useful let's take a look at some applications now of the log map so here's a good motivating question what is the center of a set of points on a curved manifold you have a set of points you want to find the center well if you look at this picture you think oh it's easy I just take the extrinsic average and I project it back onto the surface well if you have more complicated shapes this doesn't work out so well right if I take the extrinsic average and projected there was actually a lot of nice discussion about that on Tuesday I think okay so instead what we can do is define the center as a point that minimizes the sum of geodesic distances right I can define an energy that says sum up the distance between my candidate center point and all the points I want to find the center of raised to some power when this power is two we get that Karcher oarfish a mean when power is 1 we get the what's called the geometric median you think these kind of l/2 vs. l1 distances the latter one will be more robust to outliers the naive algorithm for this ok there is a simple algorithm just evaluate this sum for every point on the meant on the manifold and take the minimizer yeah you can do that it's pretty expensive so we can do something a lot better or a lot more efficient by using the log map and the idea is just a simple iterative algorithm not ours which is I pick some random initial starting point compute the log of all the points that I want to find the center of I still have vectors pointing towards the the points that I want to take the center of and I compute the mean of that vector if that vector is not zero that means there's some direction I can travel to get closer to the center so I move along the exponential map in that direction and repeat until I'm done this can be easily modified I won't write it but it can be easily modified for the case P equals 1 and I think for all P greater than 1 it's not too bad what's the cost well we just have to compute this this line in blue we just had to compute this is the bottleneck of the computation compute the log map once per iteration and that's nice that's all we have to do because it means it's no more expensive to compute the mean of a very large set of points than a very small set of points or really if I have a distribution it's really easy to compute the mean of that so yes sorry yeah mm-hmm yeah totally yep absolutely so but but that's just true of the definition of karcher means in general karcher means not unique so this will find one of them yeah you can yep right that's right yeah so we use a line search yeah yeah and that really helps yeah on the surface yeah okay so simple algorithm but it's actually quite effective so really tons to converge in about ten iterations even on some very complicated shapes so here we start you know here's here's the set of points these orange points I want to compute the mean of I started some point way out here and very quickly it zips over to the center the steps you can barely see but there are these little dots so not many steps here's an example of doing this on a distribution rather than a set of points so I'm plotting now both the cart romina and the geometric median and here you just see the robustness of this geometric media and if I add some little splotch of probability over here it doesn't shift did you metric median very very much and I kind of mentioned at the beginning that these are sort of the you know in some sense the first real algorithm to do this you might think yeah there's there's you know surely was it a way to do it before but if you really go and try to use these local approximations of the log map to implement the same algorithm to be they behaved very very so they might just wander around randomly and get stuck or maybe they wander around randomly and eventually they just get lucky and get near the cart Ramin and okay they get sucked in but but we really are able to do this with a lot of robustness and very quickly there are other ways of computing karcher like centers also but they generally just don't get the answer right okay from there okay now we can build here's a stack of applications we can build on top of the karcher the karcher means to do interesting stuff so one is to compute centroidal Voronoi diagrams okay so an ordinary Voronoi diagram is you have a set of sites on the surface and you partition it into the sets that are closest to those points centroidal just means that you are moving these sites so that they eventually end up at the center of mass of their self and you want all the sites to be at the center of mass of their cell this is actually surprisingly tricky on curved manifolds because it really just the Voronoi diagrams just don't look like they do in Euclidian space so for one thing cells don't even have to be simply connected right you can of cells with interesting topology here's like a cell that wraps around a handle they don't have to contain their centroid they can be non convex and points that are very close in our end might be very far intrinsically on the surface so if you look at algorithms for computing these centroidal or no diagrams what they normally do is they say oh well locally you know the Euclidean distance and the geodesic distance are pretty close so let's just use that but you really can get into trouble if you do that so in some sense you know this is really the first wave of computing these centroidal Voronoi diagrams for very large cells on curved surfaces okay we can also use the karcher mean for or the I guess in this case the geometric median for doing landmarks on surfaces you want to match different surfaces so what do you do well one thing to kind of think about is that shape matching algorithms often require not only salient points on the surface but also a correspondence between these points it's not just enough to just know you know where the fingertips are you have to know which one is the thumb which one is the finger and so forth the index finger and so forth so otherwise you might have to solve like a pretty hard matching problem so what are some canonical points that depend only on the geometry right we have no no auxilary data one is you take the geometric median of the entire surface so you have a uniform probability distribution on the surface what's it centered well on this on this model it's something somewhere around the bellybutton actually there's some subtlety here it's not unique again so it might be one on the back as well but you find these canonical centers and then you say from there you do furthest point sampling okay I take those that initial set of centers what's the furthest point from that one what's the furthest point from those with the furthest point from those and you get a canonical ordering of these landmarks that can be used for matching so that's kind of a nice nice idea okay so we could go on and on I mean if they think hopefully I've convinced you at this point that this quantity is useful for geometry processing just a few kind of high-level remarks just to say again you know what is the whole takeaway here there's a lot of ideas the key idea is just to say reducing geometry processing to solving linear PDE s specifically linear PDE s really opens the door to a lot of new approaches to fast geometric algorithms and that's because you can immediately benefit from anybody else's developments in numerical linear algebra in parallel solvers and all this kind of thing you can also take advantage of people who come up with better discrete differential operators right this will immediately improve the accuracy of all these algorithms that are built on top of this machinery again just want to stress the heat method is not one particular discrete algorithm right we're not just talking about you know triangle meshes that's not the heat method the heat method is kind of a smooth algorithm or general principle for generating lots of different possible discrete algorithms that can be solved in many different ways it's already fairly mature there's lots of different descritization x' ways to improve accuracy fill up fairly fast solvers widely available code but of course there's a lot of questions also that still could be answered so for one thing given the close connection with the eye kernel equation you might ask are the ways of kind of hybrid hybridizing fast marching type solvers and heat method type solvers right in terms of getting good accuracy good performance all this kind of thing can you deal with extreme and I saw aunt Inez Trippi we thought a nice talk about that on Tuesday you know I'm not not sure how well this this holds up but you think about that there's all sorts of other elliptic operators you could throw in there and just ask okay instead of the connection laplacian I throw in Dirac operator or something what geometric information can I now extract from my manifold questions of formal convergence everything I said was numerical we have no results about formal convergence under refinement how do you do this in very high dimensions right I can't if I can't really directly solve these linear systems but can i still take advantage of this principle to get information about my space doing higher-order statistics that compute it means you know you could talk about covariance all this kind of thing and I hate to say it but machine learning also right not only can you use can use machine learning to improve accuracy as I think has been done for fast marketing but the other way around can you use the same kind of idea to do machine learning on curved surfaces because if you think about what we did it's you know we keep saying we do something like a convolution right heat diffusion and then we have a little non-linearity at each point well maybe this reminds some of you of deep neural networks right okay you don't yeah you don't have this translational invariance of the of the convolution but I don't I don't think that's a that's really a precursor for doing CNN's I mean it's it's something that helps with efficiency but it's not fun of it all okay so I'm gonna end but before we get too much into machine learning thank you thank you very much [Applause]

Channel: Keenan Crane

Views: 5,923

Keywords: science, technology, geometry, heat diffusion, shortest path, geodesic, geodesic distance, parallel transport, Karcher mean, Fréchet mean, exponential map, log map, SIGGRAPH, computer graphics, algorithms, Carnegie Mellon, computer science, diffusion, Voronoi, geometric median, shape analysis, machine learning, convolution, medial axis, parallel computing

Id: 4IZ-ykGnIRc

Length: 49min 29sec (2969 seconds)

Published: Wed May 29 2019