In this paper we propose an algorithm that computes a low cardinality PItype (Positive weights and Interior nodes) algebraic cubature rule of degree n, with at most (n + 1)(n + 2)/2 nodes, over curvilinear polygons defined by piecewise rational functions. Typical examples are domains whose boundary is defined piecewise by NURBS curves or by composite Bezier curves. The key tools are a relevant but overlooked theorem of 1976 by Wilhelmsen on Tchakaloff sets, a specific in-domain algorithm for such curvilinear polygons and the sparse nonnegative solution of underdetermined moment matching systems by the Lawson-Hanson NonNegative Least Squares solver. Many numerical tests are performed to show the flexibility of this approach and the implemented MATLAB toolbox is freely available to the users, in particular for possible applications within FEM/VEM with NURBS-shaped curvilinear elements. Keywords low-cardinality algebraic cubature • PI (Positive Interior) rules • curvilinear polygons, NURBS • crossing number • winding number • Tchakaloff sets • underdetermined moment matching systems • NonNegative Least Squares • NEFEM (NURBS-Enhanced FEM) • VEM (Virtual Element Method).