We present a 3D method to numerically design a realizable metasurface, which transforms a given incident field into a radiated field that satisfies mask-type (inequality) constraints. The method is based on an integral equation formulation, with local impedance boundary condition (IBC) approximation. The procedure yields the spatial distribution of the impedance, yet the process involves the synthesis of the equivalent current only. This current is constrained to correspond to a realizable surface impedance, i.e., passive, lossless, and with reactance values bounded by practical realizability limits. The current-based design avoids any solution of the forward problem, and the impedance is obtained from the synthesized current only at the end of the process. The procedure is gradient-based, with the gradient expressed in closed form. This allows handling large metasurfaces, with full spatial variability of the impedance in two dimensions. The method requires no a priori information, and all relevant operations in the iterative process can be evaluated with O ( N log N ) complexity. Application examples concentrate on the case of on-surface excitation and far-field pattern specifications; they show designs of circular and rectangular metasurface antennas of 20 wavelengths in size, with pencil- and shaped-beam patterns, and for both circular and linear polarization.