We consider a study on regression function estimation over a bounded domain of arbitrary shapes based on triangulation and penalization techniques. A total variation type penalty is imposed to encourage fusion of adjacent triangles, which leads to a partition of the domain consisting of disjointed polygons. The proposed method provides a piecewise linear, and continuous estimator over a data adaptive polygonal partition of the domain. We adopt a coordinate decent algorithm to handle the non-separable structure of the penalty and investigate its convergence property. Regarding the asymptotic results, we establish an oracle type inequality and convergence rate of the proposed estimator. A numerical study is carried out to illustrate the performance of this method. An R software package polygram is available.