Curve fitting to unorganized data points is a very challenging problem that arises in a wide variety of scientific and engineering applications. Given a set of scattered and noisy data points, the goal is to construct a curve that corresponds to the best estimate of the unknown underlying relationship between two variables. Although many papers have addressed the problem, this remains very challenging. In this paper we propose to solve the curve fitting problem to noisy scattered data using a parallel hierarchical genetic algorithm and B-splines. We use a novel hierarchical structure to represent both the model structure and the model parameters. The best B-spline model is searched using bi-objective fitness function. As a result, our method determines the number and locations of the knots, and the B-spline coefficients simultaneously and automatically. In addition, to accelerate the estimation of B-spline parameters the algorithm is implemented with two levels of parallelism, taking advantages of the new hardware platforms. Finally, to validate our approach, we fitted curves from scattered noisy points and results were compared through numerical simulations with several methods, which are widely used in fitting tasks. Results show a better performance on the reference methods.