A polynomial parametric curve spanning known tangent end-points Let $x(t)$ and $y(t)$ be unknown polynomials (of maximum order 3) defining a parametric curve
$$ \mathbf{r}(t)=\begin{bmatrix}x(t)\\y(t)\end{bmatrix} $$
that fits known tangential end-points:
$$ \mathbf{r}(0)=\begin{bmatrix}x_0\\y_0\end{bmatrix} $$
$$ \mathbf{r}(1)=\begin{bmatrix}x_1\\y_1\end{bmatrix} $$
$$ \mathbf{r}'(0)=k_0\begin{bmatrix}x'_0\\y'_0\end{bmatrix} $$
$$ \mathbf{r}'(1)=k_1\begin{bmatrix}x'_1\\y'_1\end{bmatrix} $$
where $x_0$, $y_0$, $x_1$, $y_1$, $x'_0$, $y'_0$, $x'_1$, $y'_1$ are known constants specifying the end-point locations and directions, and $k_0$, $k_1$ are scaling factors we don't care about.
To find $\mathbf{r}(t)$ it is possible to set $k_1=k_2=1$, and then split the problem into two independent cubic Hermite interpolations in x and in y. The resulting $\mathbf{r}(t)$ technically satisfies all constraints above, but when I tested this by interpolating between tangential points drawn from the unit-circle, I got the following result.

Is there an approach that would give a more reasonable result than the above?
My problem specifically concerns data points with pre-determined directions at each data point.
 A: It look like 
the loops you get is due too big tangent vectors,
something like this:

As an heuristic 
(which might not always work well), 
try to scale given tangents 
at both ends of every segment
to fit the distance between given endpoints of the segment,
and you'll get something like this:

More examples with original curve in red:
6 segments:

7 segments: 

The example curves were obtained by
the following procedure:
given $n+1$ points 
$p_i=(x_i,y_i)$ and vectors 
$p'_i=(x'_i,y'_i)$, $i=0\,\ldots\,,n$,
construct $n$ cubic Bezier segments
\begin{align} 
s_i&\left(p_i, 
\ p_i+\tfrac13\left(\frac{\|p_{i+1}-p_i\|}{\|p'_i\|} \right)\,p'_i,
\ p_{i+1}-\tfrac13\left(\frac{\|p_{i+1}-p_i\|}{\|p'_{i+1}\|} \right)\,p'_{i+1},
\ p_{i+1}\right)
,\quad i=0,\,\ldots\,,n-1
.
\end{align}  
A: Try clamped spline interpolation. You'll get a piecewise cubic curve interpolating the given points for which you can prescribe the tangent directions at the endpoints.
A: Here is a suggested solution.  You do not need two independent functions of $t$ because you are just in a two-dimensional space.  Just focus on fitting a third order polynomial to the function $y(t)$.  Then, the vector-valued solution will automatically be $r(t) = t\hat{i} + y(t)\hat{j}$.  Does this make sense?
A: Sorry I have been swamped in a sea of work for a bit, but now as promised my contribution.

You can achieve this by linear fusion of two models,

*

*One set of discretized "snake" vectors, like proposed in question here.
(This approach also partially answers that question.)


*Polynomial regression on both coordinates separately.
You then assign each vector in the snake sum / integral : $$\cases{P_x(\Delta_t k) = \displaystyle\sum_{i=0}^k {\bf v[k]}_x\\P_y(\Delta_t k) = \displaystyle\sum_{i=0}^k {\bf v[k]}_y}$$
Then you choose the nodes in your snake to carry the end-point constraints. Especially constraints about derivatives (tangent vectors) and position are practical to carry such constraints.
Then you build a big linear least squares system out of it and solve it. This will be super general, you can even split and let different polynomials take responsibility of fitting individual areas of the snake and to fuse their end-points (kind of similar to what splines will do), you can also regularize the snake to minimize total curve length, it's first, second partial derivatives and so on.

Here is an example with a black contour and we choose to try and fix the yellow points. Then we fit a snake to the yellow points, and simultaneously force it to be pushed towards a polynomial and simultaneously regularize it's curve length, first derivatives (a bit). As we can see order 3 (red) is not big enough polynomial to fit all points.
