Minimal distance to a cube in 2D and 3D from a point lying outside This is kind of a geometrical question. For my program I want to compute the minimal distance $r$ from a given point to the cube. Here is a drawing which shows what I mean:

I have two vectors $\vec{p}$ and $\vec{q}$ which indicate the position of my two points. Point $p$ can be anywhere outside the cube. Point $q$ is exactly in the middle of the cube. The distance from point $q$ to the cubes surface is always $d$. I can easily compute $R$ which is the distance from point $q$ to point $p$. But what I need is the minimal distance $r$ from point $p$ to the cube. I am sure that I have to distinguish several cases depending on where point $p$ is located. I think there are three cases: 
1) The minimal distance $r$ from point $p$ is to the edge of the cube (as drawn in the picture)
2) The minimal distance $r$ from point $p$ is to the corner of the cube
3) The minimal distance $r$ from point $p$ is to the surface of the cube
After hours of trying to find a nice solution I hope someone can give me a hint.
 A: For an axis aligned cube, there are nice tricks.  Consider, for example, the axis aligned cube with corners $(\pm 1,\pm 1,\pm 1)$ (after a scaling and a shift, everything reduces to this case).
For a point $(x,y,z)$:


*

*If $|x|\leq 1$, $|y|\leq 1$ and $|z|>1$, then the distance is $|z|-1$ (the face is closest).

*If $|x|\leq 1$, $|y|>1$, and $|z|>1$, then the distance is $\sqrt{(|y|-1)^2+(|z|-1)^2}$ (the edge is closest).

*If $|x|>1$, $|y|>1$, and $|z|>1$, then the distance is $\sqrt{(|x|-1)^2+(|y|-1)^2+(|z|-1)^2}$ (the vertex is closest).
All other cases are similar.  To visualize what is going on, draw a square, but extend the edges into (infinite) lines.  This breaks up the space outside the box into $8$ regions, all points in each region are closest to the same edge or point.  Now, doing the same thing in three dimensions results in $26$ regions and you need to figure out which region you're in.
Pseudocode:
 if |x|<=1
 {
      if |y|<=1
           d=|z|-1
      else
      {
           if |z|<=1
                 d=|y|-1
           else
                 d=sqrt((|y|-1)^2+(|z|-1)^2)
      }
 }
 else
 {
      if |y|<=1
      {
           if |z|<=1
                d=|x|-1
           else
                d=sqrt((|x|-1)^2+(|z|-1)^2)
      }
      else
           if |z|<=1
                d=sqrt((|x|-1)^2+(|y|-1)^2)
           else
                d=sqrt((|x|-1)^2+(|y|-1)^2+(|z|-1)^2)
      }
 }

Also as @YvesDaoust mentions in the comments, this can be rewritten as 
$$
\sqrt{\max\{0,|x|-1\}^2+\max\{0,|y|-1\}^2+\max\{0,|z|-1\}^2}
$$
although, after unpacking the maximums, you get, essentially, the series of inequalities above.
A: Formulation as QP:
$$
\min \{ x^\top Q x + c^\top x \mid Ax \le b \}
$$
where
$$
d(x,p)^2 = (x-p)^2 = x^2 - 2 x^\top p + p^2
$$
so
$$
Q = I \quad c = -2p
$$
and we minimize $d(x,p)^2 - p^2$.
Further we have six half-spaces
$$
n_i^\top x \le d_i
$$
where the $n_i$ are unit normal vectors of the planes through the facets and the $d_i$ is the (signed distance) of the plane to the origin. We have $n_i = - n_j$ for opposite facets.
