The cyclide spline surface is a G1 smooth piecewise surface composed of Dupin cyclide patches, thus inheriting several favorable geometric properties of the Dupin cyclide, such as the closeness under offset operation. Due to the lack of shape flexibility of Dupin cyclides, it has been an outstanding problem to use a G1 smooth piecewise Dupin cyclide surface to model a free-form shape. We solve this problem by proposing a novel method for approximating a given free-form smooth surface with a cyclide spline surface. The key to our solution is increasing the fitting flexibility of a cyclide spline by relaxing and optimizing both the vertices and surface frames of the cyclide spline simultaneously in a global manner. Furthermore, we apply patch subdivision for local surface refinement and use spherical patches to fill in umbilical regions. The effectiveness of our method is demonstrated on a variety of surface shapes.