We propose a new approach for the joint estimation of aberration parameters and unknown object from diversity images with applications in imaging systems with extended objects as astronomical ground-based observations or solar telescopes. The motivation behind our idea is to decrease the computational complexity of the conventional phase diversity (PD) algorithm and avoid the convergence to local minima due to the use of nonlinear estimation algorithms. Our approach is able to give a good starting point for an iterative algorithm or it can be used as a new wavefront estimation method. When the wavefront aberrations are small, the wavefront can be approximated with a linear term which leads to a quadratic point-spread function (PSF) in the aberration parameters. The presented approach involves recording two or more diversity images and, based on the before mentioned approximation estimates the aberration parameters and the object by solving a system of bilinear equations, which is obtained by subtracting from each diversity image the focal plane image. Moreover, using the quadratic PSFs gives improved performance to the conventional PD algorithm through the fact that the gradients of the PSFs have simple analytical formulas.