One of the main focus of scattered data interpolation is fitting a smooth surface to a set of non-uniformly distributed data points which extends to all positions in a prescribed domain. In this paper, given a set of scattered data V ={(x i ,y i ) , i=1,...,n} R 2 over a polygonal domain and a corresponding set of real numbers n i i z 1 , we wish to construct a surface S which has continuous varying tangent plane everywhere (G 1 ) such that S(x i ,y i ) = z i . Specifically, the polynomial being considered belong to G 1 quartic Bézier functions over a triangulated domain. In order to construct the surface, we need to construct the triangular mesh spanning over the unorganized set of points, V which will then have to be covered with Bézier patches with coefficients satisfying the G 1 continuity between patches and the minimized sum of squares of principal curvatures. Examples are also presented to show the effectiveness of our proposed method.