Thank you guys very much for all the suggestions, really appreciated. So I proceeded following the method that Zakwayda described, and it partially works, meaning that on a set of about 1000 splines, roughly 5% of those seem to still come out with wrong vectors, thus messing up the tube generation (the tube gets twisted as in my OP). I suppose I made some mistake in my implementation.
I made a ring silhouette of n points, located at origin; then I compute front, up and side vectors of some points along the spline, and later on I rotate and move it to the desired points along the spline.
Here is the pseudo-C++ code I have:
vec3 up(0,1,0);
vec3 currentPoint; // the current point being accounted
vec3 nextPoint = currentPoint + aTinyIncrement; // calculated interpolating a small increment on the curve
vec3 frontVector = (nextPoint - currentPoint).Normalize();
vec3 splineLeft = crossProduct(frontVector, up);
vec3 splineUp = crossProduct(splineLeft, frontVector);
So at this point, I perform the rotation and positioning of the silhouette ring:
float angle = acos(dotProduct(frontVector, up));
matrix myMatrix = matrixRotate(newMatrix, angle, splineLeft);
for (int t=0; t<nPoints; t++)
newPoint = currentPoint + myMatrix * vec4(cylinderShape.at(t)) ;
I use the glm library for all the calculations for matrices and vectors (so you know I don't use some custom-made and possibly flawed routines); I suppose the method I came up so far is still missing something, because as I said before it works pretty well but not 100% fail-safe. So on about 1000 splines (well, tubes), I see some are still twisted at some points.
Cheers