Because there have been multiple similar questions about this recently, i create a new post with an example solution.

The code should work robustly without flipping cases for complex curves, it limits twisting and connects a closed curve seamlessly.

It does not deal with problems like self intersections on sharp corners.

I assume one could smooth twisting even more by calculating the curl of the tangent analytically, but it seems good and i'll edit the knot orientations and handles manually, so no need for me to go into this.

		if (visBezier)
		{
			struct BezierSpline4
			{
				static void Evaluate (vec &pos, vec &tangent, const vec &cp1, const vec &cp2, const vec &cp3, const vec &cp4, const float t)
				{
					float u = 1.0f - t;		
	
					float vvv = u*u*u;
					float vvu = u*u*t*3;
					float vuu = u*t*t*3;
					float uuu = t*t*t;
	
					pos = vvv * cp1 + vvu * cp2 + vuu * cp3 + uuu * cp4;

					vvv = -3*u*u;
					vvu = -6*u*t + 3*u*u; 
					vuu =  6*u*t - 3*t*t; 
					uuu =  3*t*t;

					tangent = vvv * cp1 + vvu * cp2 + vuu * cp3 + uuu * cp4;
				}
			};
			
			// track positions
			
			std::vector<vec> loop;
			loop.push_back (vec(0,0,0));
			loop.push_back (vec(0,1,0));
			loop.push_back (vec(1,1,0));
			loop.push_back (vec(1,1,1));
			loop.push_back (vec(1,0,1));
			loop.push_back (vec(0,0,1));
			loop.push_back (vec(0.5f,0.5f,0.5f));
			loop.push_back (vec(1.5f,0.5f,0.5f));
			int count = (int)loop.size();
//static vec finalP (1.5f,0.5f,0.5f);
//ImGui::DragFloat3 ("point", &finalP[0], 0.1);
//loop[count-1] = finalP;

			// create control points for smooth bezier curve  

			struct Knot
			{
				vec pos;
				vec handleL;
				vec handleR;
				vec tangent;
				quat orn;
			};
			std::vector<Knot> knots;
			knots.resize(count);
			
			for (int p=0; p<count; p++)
			{
				int c = (p+1)%count;
				int n = (p+2)%count;

				vec tang = loop[n] - loop[p];
				tang /= 6;

				knots[c].pos = loop[c];
				knots[c].handleL = loop[c] - tang;
				knots[c].handleR = loop[c] + tang;
				knots[c].tangent = tang.Unit();
			}

			// construct initial orientation for the first knot from a given up-vector

			{
				vec up(0,1,0);
				vec upPerpendicularIfFailure(1,0,0);

				vec normal = knots[0].tangent.Cross(up);
				float sql = normal.Dot(normal);
				if (sql < FP_EPSILON)
				{
					normal = knots[0].tangent.Cross(upPerpendicularIfFailure);
					sql = normal.Dot(normal);
				}
				normal /= sqrt(sql);

				matrix temp;
				temp[0] = knots[0].tangent;
				temp[1] = normal;
				temp[2] = temp[0].Cross(temp[1]);

				temp.Rotation()->ToQuat (knots[0].orn);
			}




			// transport initial orientation along the curve

			float angleDefect = 0; // difference between first and last orientation, which we want to be zero for a loop
			for (int p=0; p<count; p++)
			{
				int c = (p+1)%count;

				// rotate previous orientation from previous tangent to current tangent
				quat rot;
				rot.FromVecToVec (knots[p].tangent, knots[c].tangent);
				//quat cOrn = knots[p].orn * rot; // i assume this multiplication order would be more common
				quat cOrn = rot * knots[p].orn; // ...than mine

				if (c==0) // measure angle defect around x axis (tangent)
				{
					matrix m0; m0.Rotation()->FromQuat (cOrn);
					matrix m1; m1.Rotation()->FromQuat (knots[c].orn);
					angleDefect = atan2 (m0[1].Dot(m1[1]), m0[2].Dot(m1[1]));

					//RenderQuat (knots[c].pos, knots[c].orn);
					//RenderQuat (knots[c].pos, cOrn);
					//quat rot;
					//rot.FromAxisAndAngle (vec(1,0,0), -angleDefect);
					//quat proof = cOrn * rot;
					//RenderQuat (knots[c].pos, proof);
				}
				else
				{
					knots[c].orn = cOrn;
				}
			}

			// distribute the angle defect correction over the whole spline - otherwise the last segment could appear twisted

			for (int i=0; i<count; i++)
			{
				float t = float(i) / float(count);
				quat rot;
				rot.FromAxisAndAngle (vec(1,0,0), -angleDefect * t);
				knots[i].orn = knots[i].orn * rot; // reversed order, as we rotate locally around the tangent
			}




			// define cross section shape

			struct Vertex
			{
				vec pos;
				vec norm;
				vec2 uv;
			};
			std::vector<Vertex> crossSection;

			float w = 0.08f;
			float h = 0.04f;
			crossSection.push_back( Vertex {vec(0,-w,-h), vec(0,-0.7f,-0.7f), vec2(0,0)} );
			crossSection.push_back( Vertex {vec(0, w,-h), vec(0, 0.7f,-0.7f), vec2(0.4f,0)} );
			crossSection.push_back( Vertex {vec(0, w, h), vec(0, 0.7f, 0.7f), vec2(0.6f,0)} );
			crossSection.push_back( Vertex {vec(0,-w, h), vec(0,-0.7f, 0.7f), vec2(0.8f,0)} );
			crossSection.push_back( Vertex {vec(0,-w,-h), vec(0,-0.7f,-0.7f), vec2(1,0)} );
			int csCount = (int)crossSection.size();


			// transport shape along spline

			std::vector<Vertex> ring0 = crossSection, ring1 = crossSection;

			for (int l=0; l<=count; l++)
			{
				int c = l % count;
				int n = (c+1) % count;

				int tess = 12;
				for (int i=0; i<=tess; i++)
				{
					float t = float(i) / float(tess);

					vec pos, tangent;
					BezierSpline4::Evaluate (pos, tangent, knots[c].pos, knots[c].handleR, knots[n].handleL, knots[n].pos, t);

					quat orn; orn.Slerp (knots[c].orn, knots[n].orn, t);

					if (1) // optional: snap orientation x axis to spline tangent
					{
						vec x = orn.Rotate(vec(1,0,0));
						quat rot;
						rot.FromVecToVecNonUniformLength (x, tangent);
						orn = rot * orn; // again, flip multiplication order if necessary
					}

					// transform vertices

					matrix xf; 
					xf.Rotation()->FromQuat(orn);
					xf[3] = pos;

					for (int j=0; j<csCount; j++)
					{
						Vertex &vert = ring0[j];
						vert.pos = xf.Transform (crossSection[j].pos);
						vert.norm = xf.Rotate (crossSection[j].norm);
						vert.uv = vec2 (crossSection[j].uv[0], t);
						 
						// some uv hackery for my visuals...
						vert.uv = vec2 (1,1) + vert.uv * 4;
						vec const lightDir = vec(1,-2,-1).Unit();
						float d = lightDir.Dot( vert.norm );
						if (d>=0) d*=d; else d = -sqrt(-d);
						vert.uv[2] = pow(d * 0.5f + 0.5f, 2.2f);
					}

					// render triangles

					if (i) for (int j=0; j<csCount; j++)
					{
						int k = (j+1)%csCount;
						Vertex &v0 = ring0[j];
						Vertex &v1 = ring0[k];
						Vertex &v2 = ring1[k];
						Vertex &v3 = ring1[j];
						
						simpleVis.RenderTriangle (
							(float*)&v0.pos, (float*)&v0.uv, 
							(float*)&v2.pos, (float*)&v2.uv, 
							(float*)&v1.pos, (float*)&v1.uv);
						simpleVis.RenderTriangle (
							(float*)&v0.pos, (float*)&v0.uv, 
							(float*)&v3.pos, (float*)&v3.uv, 
							(float*)&v2.pos, (float*)&v2.uv);
					}

					ring1 = ring0;
				}

			}
		}
                                                    
                                                    
                                                    
// some quaternion code:
                                                    
	__forceinline void FromVecToVec (const sVec3 &dir0, const sVec3 &dir1)
	{
		sScalar dot = dir0.Dot(dir1);
		(*this)[3] = 1.0f + dot;
		if ((*this)[3] >= FP_EPSILON)
		{ 
			*((sVec3*)this) = dir0.Cross(dir1);
			Normalize();
		}
		else
		{
			if (dot<0) *this = sQuat (1.0f, 0, 0, 0); // 180 degrees apart -> 180 degrees rotation around x
			else *this = sQuat (0, 0, 0, 1.0f); // zero vector input -> zero rotation output
		}
	}
	__forceinline void FromVecToVecNonUniformLength (const sVec3 &dir0, const sVec3 &dir1)
	{
		sScalar dot = dir0.Dot(dir1);
		(*this)[3] = sqrt(dir0.SqL() * dir1.SqL()) + dot;
		if ((*this)[3] >= FP_EPSILON)
		{ 
			*((sVec3*)this) = dir0.Cross(dir1);
			Normalize();
		}
		else
		{
			if (dot<0) *this = sQuat (1.0f, 0, 0, 0); // 180 degrees apart -> 180 degrees rotation around x
			else *this = sQuat (0, 0, 0, 1.0f); // zero vector input -> zero rotation output
		}
	}