Functional magnetic resonance imaging (fMRI) functional connectivity between brain regions is often computed using parcellations defined by functional or structural atlases. Typically, some kind of voxel averaging is performed to obtain a single temporal correlation estimate per region pair. However, several estimators can be defined for this task, with various assumptions and degrees of robustness to local noise, global noise, and region size. In this paper, we systematically present and study the properties of 9 different functional connectivity estimators taking into account the spatial structure of fMRI data, based on a simple fMRI data spatial model. These include 3 existing estimators and 6 novel estimators. We demonstrate the empirical properties of the estimators using synthetic, animal, and human data, in terms of graph structure, repeatability and reproducibility, discriminability, dependence on region size, as well as local and global noise robustness. We prove analytically the link between regional intra-correlation and inter-region correlation, and show that the choice of estimator has a strong influence on inter-correlation values. Some estimators, including the commonly used correlation of averages (ca), are positively biased, and have more dependence to region size and intra-correlation than robust alternatives, resulting in spatially-dependent bias. We define the new local correlation of averages estimator with better theoretical guarantees, lower bias, significantly lower dependence on region size (Spearman correlation 0.40 vs 0.55, paired t-test T=27.2, p=1.1e ^-47 ), at negligible cost to discriminative power, compared to the ca estimator. The difference in connectivity pattern between the estimators is not distributed uniformly throughout the brain, but rather shows a clear ventral-dorsal gradient, suggesting that region size and intra-correlation plays an important role in shaping functional networks defined using the ca estimator, and leading to non-trivial differences in their connectivity structure. We provide an open source R package and equivalent Python implementation to facilitate the use of the new estimators, together with preprocessed rat time-series.