A statistical mechanical integral equation theory is developed to describe the average structure of a polar liquid around a complex molecular solute of irregular shape. The integral equation is formulated in three-dimensional Cartesian coordinates from the hypernetted chain (HNC) equation for a solute at infinite dilution. The direct correlation function of the pure solvent used in the theory is taken from the analytical solution of the mean spherical approximation (MSA) equation for a liquid constituted of nonpolarizable hard spheres with an embedded dipole at their center. It is demonstrated explicitly that, in the limit where the size of the solvent particles becomes very small, the present theory reduces to the well-known equations for macroscopic electrostatics in which the solvent is represented in terms of a dielectric continuum. A linearized version of the integral equation corresponds to a three-dimensional extension of the familiar MSA equation. This 3D-MSA integral equation is illustrated with numerical applications to the case of an ion, a water molecule, and N-methylacetamide. The numerical solution, obtained on a discrete three-dimensional cubic grid in Cartesian coordinates, yields the average solvent density and polarization density at all the points x,y,z around the solute. All spatial convolutions appearing in the theory are calculated using three-dimensional numerical fast Fourier transforms (FFT).