The following code draw curve from predefined controlled points, it is correct but the curve smoothness is currently G1, that is it has segments at every step other than the intermediate joint where they touch the controlling point segments. Can you improve it smoothness? CurvePoints Curve::evalBezier(const std::vector< Vector3f >& P, unsigned steps){

    // Check

    if( P.size() < 4 || P.size() % 3 != 1 )
    {
        std::cerr << "evalBezier must be called with 3n+1 control points." << std::endl;
        exit( 0 );
    }

	// 2.3.1 - Create the functionality to process Bezier control points


    // You should implement this function so that it returns a Curve
    // (e.g., a vector< CurvePoint >).  The variable "steps" tells you
    // the number of points to generate on each piece of the spline.
    // At least, that's how the sample solution is implemented and how
    // the SWP files are written.  But you are free to interpret this
    // variable however you want, so long as you can control the
    // "resolution" of the discretized spline curve with it.

    // Make sure that this function computes all the appropriate
    // Vector3fs for each CurvePoint: V,T,N,B.
    // [NBT] should be unit and orthogonal.

    // Also note that you may assume that all Bezier curves that you
    // receive have G1 continuity.  Otherwise, the TNB will not be
    // be defined at points where this does not hold.

    std::cerr << "\t>>> evalBezier has been called with the following input:" << std::endl;

    std::cerr << "\t>>> Control points (type vector< Vector3f >): "<< std::endl;
    for( unsigned i = 0; i < P.size(); ++i )
    {
        std::cerr << "\t>>> " << P[i] << std::endl;
    }

    std::cerr << "\t>>> Steps (type steps): " << steps << std::endl;
    //std::cerr << "\t>>> Returning empty curve." << std::endl;

    // Right now this will just return this empty curve.
    CurvePoints res;
    Matrix4f M_B ( 1, -3,  3, -1,
                   0,  3, -6,  3,
                   0,  0,  3, -3,
                   0,  0,  0,  1);
                   
    Vector3f B_prev(0, 0, 1);
                   
    for (unsigned i = 0; i < P.size() - 3; i += 3) {
        Matrix4f G(
            P[i][0], P[i + 1][0], P[i + 2][0], P[i + 3][0],
            P[i][1], P[i + 1][1], P[i + 2][1], P[i + 3][1],
            P[i][2], P[i + 1][2], P[i + 2][2], P[i + 3][2],
            1      , 1          , 1          , 1
        );
        
        
        for (unsigned j = 0; j <= steps; j++) {
            if ((j == 0 || j == steps) && (i != 0) && (i + 3 != P.size())) continue;
            float t = float(j)/ steps; // t in [0, 1]
            Vector4f M(1, t, t*t, t*t*t); // Monomial basis
            // position
            Vector3f V = (G * M_B * M).xyz();
            
            // tangent
            Vector4f dM(0, 1, 2*t, 3*t*t);
            Vector3f T = (G * M_B * dM).xyz();
            T.normalize();
            
            // Normal
            Vector3f N = Vector3f::cross(B_prev, T);
            N.normalize();
            
            // Update binormal
            Vector3f B = Vector3f::cross(T, N);
            B.normalize();
            res.push_back(CurvePoint{V, T, N, B});
            
            B_prev = B;
        }
    }
      
    return res;

}