We discuss an algorithm (implemented in Matlab) that computes numerically total-degree bivariate orthogonal polynomials (OPs) given an algebraic cubature formula with positive weights, and constructs the orthogonal projection (hyperinterpolation) of a function sampled at the cubature nodes. The method is applicable to nonstandard regions where OPs are not known analytically, for example convex and concave polygons, or circular sections such as sectors, lenses and lunes
