A new method for extracting neutron densities from intermediate energy elastic proton-nucleus scattering observables uses a global Dirac phenomenological (DP) approach based on the Relativistic Impulse Approximation (RIA). Data sets for 40 Ca, 48 Ca and 208 Pb in the energy range from 500 MeV to 1040 MeV are considered. The global fits are successful in reproducing the data and in predicting data sets not included in the analysis. Using this global approach, energy independent neutron densities are obtained. The vector point proton density distribution, ρ p v , is determined from the empirical charge density after unfolding the proton form factor. The other densities, ρ n v , ρ p s , ρ n s , are parametrized. This work provides energy independent values for the RMS neutron radius, R n and the neutron skin thickness, S n , in contrast to the energy dependent values obtained by previous studies. In addition, the results presented in paper show that the expected rms neutron radius and skin thickness for 40 Ca is accurately reproduced. The values of R n and S n obtained from the global fits that we consider to be the most reliable are given as follows: for 40 Ca, 3.314 > R n > 3.310 fm and −0.063 > S n > −0.067 fm; for 48 Ca, 3.459 > R n > 3.413 fm and 0.102 > S n > 0.056 fm; and for 208 Pb 5.550 > R n > 5.522 fm and 0.111 > S n > 0.083 fm. These values are in reasonable agreement with nonrelativistic Skyrme Hartree-Fock models and with relativistic Hartree-Bogoliubov models with density-dependent meson-nucleon couplings. The results from the global fits for 48 Ca and 208 Pb are generally not in agreement with the usual relativistic mean-field models.