Future sixth-generation (6G) networks are envisioned to provide both sensing and communications functionalities by using densely deployed base stations (BSs) with massive antennas operating in millimeter wave (mmWave) and terahertz (THz). Due to the large number of antennas and the high frequency band, the sensing and communications will operate within the near-field region, thus making the conventional designs based on the far-field channel models inapplicable. This paper studies a near-field multiple-input-multiple-output (MIMO) radar sensing system, in which the transceivers with massive antennas aim to localize multiple near-field targets in the three-dimensional (3D) space. In particular, we adopt a general wavefront propagation model by considering the exact spherical wavefront with both channel phase and amplitude variations over different antennas. Besides, we consider the general transmit signal waveforms and also consider the unknown cluttered environments. Under this setup, the unknown parameters to estimate include the 3D coordinates and the complex reflection coefficients of the multiple targets, as well as the noise and interference covariance matrix. Accordingly, we derive the Cramér-Rao bound (CRB) for estimating the target coordinates. Next, to facilitate practical localization, we propose an efficient estimator based on the 3D approximate cyclic optimization (3D-ACO), which is obtained following the maximum likelihood (ML) criterion. Finally, numerical results show that considering the exact antenna-varying channel amplitudes achieves more accurate CRB as compared to prior works based on constant channel amplitudes across antennas, especially when the targets are close to the transceivers. It is also shown that the proposed estimator achieves localization performance close to the derived CRB, thus validating its superior performance.