As you say, four points (without their associated "$t$" values) are not enough to define a planar cubic Bezier curve.
In fact, 6 points are needed. But, on the other hand, the 6 points have to satisfy certain conditions, or else the required curve does not exist. In your case, you know that the curve exists, because it was used to generate the points in the first place.
The full story is given in this paper: J. Kozak, M. Krajnc, Geometric interpolation by planar cubic polynomial curves, Comput. Aided Geom. Des., 24 (2007), pp. 67-78. You can get a copy of the paper at the authors website.
Another approach (instead of the clever techniques suggested in the paper) is to use brute-force numerical methods. Start with 6 points $X_1, \ldots, X_6$ and 6 parameter values $t_1, \ldots, t_6$. We can assume that $t_1 = 0$ and $t_6 = 1$, but $t_2, t_3, t_4, t_5$ are unknown.We want to find control points $P_1, P_2, P_3, P_4$ for our curve. Obviously we can set $P_1 = X_1$ and $P_4 = X_6$, but $P_2$ and $P_3$ are unknown. We can (numerically) find values of $t_2, t_3, t_4, t_5$ and $P_2$ and $P_3$ that minimize the quantity
$$ \sum_{i=2}^{i=5}\Vert C(t_i) - X_i \Vert^2 $$
where $C(t)$ is the Bezier curve with control points $P_1, P_2, P_3, P_4$. This is a fairly nasty non-linear problem, but a good numerical algorithm will usually converge to the desired solution, given a decent starting point.
Another numerical approach: with the same notation we used above, solve the equations:
$$ C(t_i) = X_i \quad (i = 2,3,4,5)$$
Since the $X_i$ are 2D points, there are 8 equations here. We have 8 unknowns -- $t_2, t_3, t_4, t_5$ and the $x$ and $y$ coordinates of $P_2$ and $P_3$. The equations are non-linear, but a good numerical root-finding package should be able to handle them.
A much better approach is to get the Bezier curves from the CAD package. If the design guys won't give you the curves, you should charge them more money for your services, since they are making your life so much more difficult. Their choice -- curves or cash :-)
What you need is something called "Boehm's algorithm" (after its originator, Wolfgang Boehm). It has a simple geometric interpretation, and drawing a few pictures should make it clear. There is a pretty good explanation (with pictures) in this document.
The algorithm is based on a process called "knot insertion". You keep inserting knots into the b-spline curve until each knot has multiplicity 3. Then, the b-spline control points of this refined curve give you the Bezier control points of its segments.
So, if you're writing code to do this, one approach is to write a knot insertion function first, and then call it repeatedly.
There is knot insertion code here.
Best Answer
Let $F(s,t)=(s/2+t)^2-t$. The equation $F(s,t)=0$ describes the locus of the points of the plane lying on the given quadratic Bézier curve (which is nothing but a parabola).
If $A$ and $B$ are two points of the plane such that $F(A)>0$ and $F(B)<0$ (I hope the notation is clear), then any continuous curve joining $A$ with $B$ must intersect the parabola by Bolzano's theorem, because $F$ is a continuous function. Hence $A$ and $B$ lie on opposite sides of the parabola: equation $F(s,t)<0$ describes the locus of points which are on the same side of the parabola as $B$, while $F(s,t)>0$ describes the locus of points which are on the same side of the parabola as $A$.
Observe now that the midpoint $M$ of $P_1P_2$ lies on the region enclosed by the curve and the triangle, and its barycentric coordinates are $s=0$, $t=1/2$, so that $F(M)=-1/4<0$. By the above reasoning, all the other points in that region are on the same side of the parabola as $M$ and have then $F(s,t)<0$.
EDIT.
For a quadratic Bézier curve joining $P_1$ with $P_2$ and having $P_0$ as control point, we want to prove that $(s/2+t)^2=t$, where $s$ and $t$ are the barycentric coordinates (relative to $P_0$ and $P_1$)) of any point $P$ on the curve.
I'll consider the simple case where $P_1=(-1,1)$, $P_2=(1,1)$ and $P_0=(-1,0)$, so that the Bézier curve is the parabola of equation $y=x^2$ (see figure below). The general case can be obtained from that by an affine transformation, which preserves both the barycentric coordinates of point $P$ and its lying on the Bézier curve.
Let then $P=(x,x^2)$ be any point on the parabola. It's easy to find that $$ PA=1-x^2,\quad MC=MB={1\over2}MP_0={1+x^2\over2},\quad PB={1+x^2\over2}-x={(1-x)^2\over2}. $$ Barycentric coordinates $s$ and $t$ are respectively the areas of triangles $PP_1P_2$ and $PP_0P_2$, divided by the area of $P_0P_1P_2$ (which is $2$). Hence: $$ s={1\over2}PA={1-x^2\over2},\quad t={1\over2}PB={(1-x)^2\over4}. $$ It is now easy to check that: $$ \left({s\over2}+t\right)^2=\left({1-x^2\over4}+{1-2x+x^2\over4}\right)^2 =\left({1-x\over2}\right)^2=t. $$