Distance between 2 points on a Bézier surface 16 control points as input for bicubic Bézier surface and then any 2 points on that surface. These are the data I have. 
I need to compute distance between those 2 points on the Bézier surface. Please keep in mind that the distance between those two points must not be euclidean distance. It should be the distance along the Bézier surface.
Someone please provide me a general method to do this work.
 A: Look up the term "geodesic", on this site and elsewhere. You will find various numerical techniques. They are typically techniques that work with any parametric surface, not just with Bezier surfaces. As far as I know, there are no properties of bicubic Bezier surfaces that simplify the problem, so you lose nothing by using a general algorithm.
A good place to start might be the answer to this question, or this web page, or this paper.
I added the "geodesic" tag to your question.
A: Let $\{P_{mn} \in \mathbb{R}^3 \, | \, m,m = 0,1,2,3\}$ bet the $16$ points in three space, where 
$$P_{mn} = \left(\matrix{P_{mn}^1 \\ P_{mn}^2\\ P_{mn}^3}\right) \, \in \mathbb{R}^3 \,\,\, \text{ where } \,\,\, m,n = 0,1,2,3$$
$$F \, : \, [0,1]\times[0,1] \, \to \,  \mathbb{R}^3$$
$$u = \left(\matrix{u_1\\ u_2}\right) \mapsto  \left(\matrix{F^1\\ F^2\\ F^3}\right)$$
$$F(u_1,u_2) = \sum_{n=0}^3\sum_{m=0}^3 \,\, \binom{3}{m}\binom{3}{n} \, u_1^m\,(1-u_1)^{3-m} \, u_2^n \, (1-u_2)^{3-n} \,\, P_{mn}$$ so we can think of $F$ as a vector of three bicubic functions
$$F(u) = F(u_1,u_2) = \left(\matrix{F^1(u_1,u_2) \\ F^2(u_1,u_2) \\F^3(u_1,u_2)} \right) \, \in \mathbb{R}^3 \,\,\, \text{ for } \,\,\, u = (u_1,u_2) \in [0,1]\times[0,1]$$
where
$$F^k(u) = F^k(u_1,u_2) = \sum_{n=0}^3\sum_{m=0}^3 \,\, \binom{3}{m}\binom{3}{n} \, u_1^m\,(1-u_1)^{3-m} \, u_2^n \, (1-u_2)^{3-n} \,\, P^k_{mn}$$
for $k = 1,2,3$


*

*Compute the two vector function $$\partial\frac{\partial F}{\partial u_1}(u)= \left(\matrix{\frac{\partial F^1}{\partial u_1}(u) \\ \frac{\partial F^2}{\partial u}(u) \\  \frac{\partial F^3}{\partial u}(u) }\right) \,\,\,\text{ and }\,\,\,\, \frac{\partial F}{\partial u_2}(u)= \left(\matrix{\frac{\partial F^1}{\partial u_2}(u) \\ \frac{\partial F^2}{\partial u_2}(u) \\  \frac{\partial F^3}{\partial u_2}(u) }\right)$$ where $\frac{\partial F}{\partial u_1}(u)$ and $\frac{\partial F}{\partial u_2}(u)$ are two transverse tangent vectors to the surface $F$ at the point $F(u)$. 

*Calculate the dot products 
$$g_{11}(u) = \left(\frac{\partial F}{\partial u_1}(u) \cdot \frac{\partial F}{\partial u_1}(u)\right) =  \left(\frac{\partial F}{\partial u_1}(u)\right)^T \, \frac{\partial F}{\partial u_1}(u)$$
$$g_{12}(u) = g_{21}(u) = \left(\frac{\partial F}{\partial u_1}(u) \cdot \frac{\partial F}{\partial u_2}(u)\right) = \left(\frac{\partial F}{\partial u_1}(u)\right)^T \, \frac{\partial F}{\partial u_2}(u)$$
$$g_{22}(u) = \left(\frac{\partial F}{\partial u_2}(u) \cdot \frac{\partial F}{\partial u_2}(u)\right) = \left(\frac{\partial F}{\partial u_2}(u)\right)^T \, \frac{\partial F}{\partial u_2}(u) $$
This provides you with the Riemanina metric 
$$g(u) =  \left(\, \matrix{g_{11}(u) & g_{12}(u) \\ g_{21}(u) & g_{22}(u)}\,\right) $$ on the coordinate patch $u = (u_1,u_2) \in [0,1]\times[0,1]$

*Next I am going to reference myself and I am going to just copy and paste an answer I gave to a question here: Computing geodesics (or shortest paths).
I hope it helps you. If you questions, feel free to ask. Just have in mind that there the notation $S$ is the same as the notation $F$ in this current post.
