A variety of formation scenarios have been proposed to explain the diversity of properties observed in bulges. Studying their intrinsic shape can help to constrain the dominant mechanisms at the epochs of their assembly. The structural parameters of a magnitude-limited sample of 148 unbarred S0-Sb galaxies were derived in order to study the correlations between bulges and disks, as well as the probability distribution function of the intrinsic equatorial ellipticity of bulges. We present a new fitting algorithm (GASP2D) to perform two-dimensional photometric decomposition of the galaxy surface-brightness distribution. This was assumed to be the sum of the contribution of a bulge and disk component characterized by elliptical and concentric isophotes with constant (but possibly different) ellipticity and position angles. Bulge and disk parameters of the sample galaxies were derived from the J-band images, which were available in the Two Micron All Sky Survey. The probability distribution function of the equatorial ellipticity of the bulges was derived from the distribution of the observed ellipticities of bulges and misalignments between bulges and disks. (1 data file).