   Closest point on polyhedron defined by intersection of halfspaces to a point
  Hi, first post here :)Mathematically formulated, my question is: Given any point p in 3D and a set of normal vectors n_1,...,n_k, each defining a halfspace H_i := {x in R^3 | dot(n_i, x) <= 0}, what is the point q in the intersection volume V = intersection(H_1,...,H_k) that is closest to p?Specifically I'm interested in how to compute q _fast_, since I believe I already know how.If V contains p (which can be decided in O(k)), of course q = p, so let's assume that p is not in V. Note that all planes defining V contain the origin, so V basically looks like an infinite cone, with infinite edges starting from the origin, and flat surfaces between those edges. The crux is that we don't know these edges in advance, which would allow for an O(k) implementation by doint a simple feature test against the faces and edges of the cone.So one way to solve the problem is to build the mesh first, which I believe can be done in O(k*log(k)), and then feature test p against the mesh to get q. I have an algorithm in mind that performs at O(k*log(k)) by exploiting that all edges start at the origin, so I can incrementally intersect the planes which gives me an intermediate set of edges each time, but the edges can be ordered eg. clockwise and stored in an array. Then computing the next plane intsersction can be done in O(log(k)) by binary searching in the edges array. So doing that for k planes gives O(k*log(k)).But this is just my first guess, so I hope this can be done even faster (Any ideas?). In my research on this topic I got the notion that the "halfspace intersection mesh building" is the dual problem of computing the convex hull for a set of points, but I couldn't find an in-depth explanation for that yet. So any pointers on that topic are appreciated, since the dual formulation may give me a better understanding of the problem.Some random remarks: I need to solve the problem to implement a character controller to move a character (sphere, capsule, ...) through a polygon soup. I know there are a lot posts about this topic on the forum, but I couldn't find one that handles multiple contact points simultaneously. In this context, the point p is the desired motion vector where I want to move the character, and the normals n_1,...,n_k are the contact normals pointing away from the character. To prevent interpenetration I need to compute q, which will be the actual motion vector. So no interiour point methods please, since I need real time computation speed ;)   
  I am not quite sure if I got you right.If I think about your description then one thingflashs up my mind. You try to take several normalsand ONE point to create n plane.By definition your point p will always be part of thehalfspace dot(n_i, x) <= 0.So your point p will always be contained in V.I think you want to describe another problem.Please describe your problem more detailed.   ;@Spynacker:Quote:Original post by SpynackerYou try to take several normalsand ONE point to create n plane.[...]I think you want to describe another problem.Please describe your problem more detailed.His problem is well defined.He has several planes, each with a different normal, but all passing through the origin.  Or, more accurately, he has a number of half-planes (a half-plane is the set of all points on a particular side of a plane).The point he refers to is not on any of the planes; rather, he wants to find the closest point in the intersection of these half-spaces to the given point.More precisely, given vectors n1,n2,...,nN in R^3 and p in R^3, he wants to find  argmin  ||p - q||^2    q  s.t.    n1^T q > 0    n2^T q > 0   .   .   .    nN^T q > 0  .Here I am using standard mathematical notation for expressing constrained minimization problems; if you're unfamiliar with it, I explain it here.This minimization problem is a quadratic programming (QP) problem, and mejiwa's question is how to solve it efficiently.  The biggest piece of structure he has which may enable him to do better than using general-purpose QP algorithms is that he's working in three dimensions; this enables him to do things like sort edges in clockwise vs. counterclockwise order.@mejiwa:We're getting a little out of my area of expertise -- and your algorithm seems quite fast so I doubt I can improve on it -- but I'll try to be of some small help.First, I notice that the problem you've formulated is a QP, whereas in most of the literature I've seen on collision handling what people are actually concerned with are LMCPs, which they solve using Lemke's algorithm.  These IIRC can be solved as QPs, but have additional structure (hence the specialized algorithm).  Are you sure that you don't want to be doing what they're doing?Second: To formulate the dual problem (if by "dual" they mean "dual" in the optimization sense),1 - write the Lagrangian for the optimization problem above (y'know, Lagrange multipliers).2 - Solve for q in dL/dq=0 to express it entirely in terms of normals and Lagrange multipliers and plug this back in to eliminate 'q.'3 - Now you have an expression entirely in terms of normals and Lagrange multipliers.  The dual problem is to maximize this expression s.t. all Lagrange multipliers are positive.I've got to go; if that's not a good enough explanation of the dual problem I can use some actual math symbols later. :-)...later:Ok, here's how it's done.The Lagrangian is (rescaling the cost function by 1/2),L = (1/2)||p - q||^2 + y1 n1^T q + y2 n2^T q + ... + yN nN^T qwhere y1,...,yN are the Lagrange multipliers.  Then,dL/dq = p^T + q^T + y1 n1^T + ... + yN nN^T = 0soq = -p - y1 n1 - ... - yN nN .Substituting this back into the Lagrangian,L = (1/2)||2 p + y1 n1 + ... + yn nN|| +    y1 n1^T[-p - y1 n1 - ... - yN nN] + ...    yN nN^T[-p - y1 n1 - ... - yN nN]= (1/2)||2 p + y1 n1 + ... + yn nN||  - (y1 n1 + ... + yN nN)^T p  - y1 y1 n1 - ... - y1 yN nN  ...  - yN y1 n1 - ... - yN yN nNwhich you now want to maximize w.r.t. y1,...,yN.  I'll leave it to you to simplify this. :-)I don't immediately see the connection to convex hulls, so maybe this is not what is meant.  I'll think about this and post back...[Edited by - Emergent on April 12, 2010 9:33:25 AM];  Oh sorry I was a bit in hurry and messed up x and p.So just ignore my post as it is related to some other problem. ^^   ;@Emergent: Sorry for the late reply, I was in university all the whole day.First, thanks for transforming it to the dual problem in terms of quadratic programming. I don't understand yet what exactly the dual problem is, so I will have to read in some books to fully understand _why_ you performed these transformation steps.What I had in mind when speaking about the dual problem was actually something different (but the dual problem you provided could be as valuable as well). As I said I have not found anything that really defines this duality, just a few hints like this description of the qhalf algorithm of the qhull program, which says "Qhull computes a halfspace intersection by the geometric duality between points and halfspaces", and the book Computational geometry: algorithms and applications also talks about point-line duality, just in 2D (see page 253).Quote:Here I am using standard mathematical notation for expressing constrained minimization problems; if you're unfamiliar with it, I explain it here.I understand them, but I still have to learn a lot about optimization. This topic just catched my interest :)Quote:This minimization problem is a quadratic programming (QP) problem, and mejiwa's question is how to solve it efficiently. The biggest piece of structure he has which may enable him to do better than using general-purpose QP algorithms is that he's working in three dimensions; this enables him to do things like sort edges in clockwise vs. counterclockwise order.Like you said this is actually a QP problem, and yeah, working in 3D  seems to be important, and I think it's also important that we're trying to minimize over a cone (in 2D it's easy to find an O(k) algorithm, just project p on the plane where dot(p,n_i) is maximal to obtain a point q0, then set q:=q0 if the distance |p-0| > |p-q0|, otherwise set q:=0). I tried to apply the 2D algorithm to 3D, but when in that case q0 is not on the hull of the volume (which might easily happen), I'm stuck. I only know that in that case q has to be on an edge of the hull, and not on a surface (at least I couldn't construct a counterexample in my head yet). But there are many edges to test, so I'm kinda lost there.Quote:We're getting a little out of my area of expertise -- and your algorithm seems quite fast so I doubt I can improve on it -- but I'll try to be of some small help.You definitely are :)Quote:First, I notice that the problem you've formulated is a QP, whereas in most of the literature I've seen on collision handling what people are actually concerned with are LMCPs, which they solve using Lemke's algorithm. These IIRC can be solved as QPs, but have additional structure (hence the specialized algorithm). Are you sure that you don't want to be doing what they're doing?I think LMCP is used for more involved things like friction and joints. In my case an LCP formulation should suffice, if I'm not mistaken. I know Lemke is fast in most cases, but it's actually an O(3) algorithm, so theoretically my algorithm is already faster. But I will try to convert it to an LCP problem, maybe I can spot some invariant there that wasn't obvious before :)This morning I also stumbled over this page, they call the problem "Projection on Polyhedral Cone". I didn't have the time to read more than the introduction yet, which says that it's still an open problem to find a formula that solves my problem (formula != algorithm (!)), I will update this post when I find anything about an efficient algorithm on that page.   ;  I found this paper by John E. Lloyd which solves the problem for contact forces when velocity and friction are involved by transforming it into an LCP and using Lemke's Algorithm to solve it. It claims to have an expected performance of O(n + m*n) where m is the number of bodies and n the number of contacts, so in my case (one character object at zero speed and a static world geometry, no friction) this has an expected complexity of O(n) (I couldn't find the worst case complexity yet), which is the best you can hope for.It's a bit lengthy (30 pages), so I will need some time to understand all the details and find the bits of information that I need. If in the meantime anyone thinks to have found another clever algorithm other than that in the paper, please let me know [smile]   