Does Pearson correlation require removal of bivariate or univariate outliers? Does the Pearson's correlation estimator require no bivariate outliers, or no outliers in each of two individual vectors of data? The answer will impact on how I winsorize outliers before calculating my Pearson estimate. 
 A: You will have to remove both. But that will not be enough. You will also have to remove those observations that are outlying on any projections of your data along any direction in $\mathbb{R}^2$. This is because of so-called multivariate outliers.
It is possible for outliers to depart significantly from the pattern of the majority of the data without necessarily standing out on any of your variables taken individually (As an example, consider the cluster of red dots in the plot attached to this answer).
Because of this, multi (in this case bi)-variate outliers can in general not be reliably detected using a coordinate-wise approaches. The only reliable way of finding them is in fact to use a multivariate trimming/winsorizing approach (i.e. one that considers univariate projections of your data along all directions). Many such method exists and you will find good implementations in most modern statistical packages (R, MATLAB, STATA, SAS,....).
It is important to recognize that multivariate outliers are equally effective at wrecking the Pearson correlation as their better known coordinate-wise cousins. An example below illustrates this.
Geometrically, coordinate-wise trimming/winsorizing amounts to drawing a rectangle around the majority of your data and considering any observations outside that rectangle as outlying. In contrast, multivariate trimming amounts to drawing an ellipse around the majority of your data and considering any point outside of that ellipse as an outlier. 
The former approach can only detect outliers if they are outlying on at least one of the coordinates (or in other words along a direction parallel to an axis of a scatter-plot of your data). The second approach, in contrast, does not suffer from this limitation. In other words, multivariate trimming approaches can detect outliers regardless of the multivariate direction in which they are outlying (and this include directions parallel to the an axis of a scatter-plot of your data).
Consider this example (the code to reproduce it is below this post):

The Pearson correlation computed on the full data (that is for the black+red points considered together) is 0.5. In this case, none of the observations will be down-weighted by a coordinate-wise winsorizing approach since for all observations $1\leq i\leq n$, the distance 
$$\max\left(\frac{|x_{i1}-\text{median}(x_1)|}{\text{mad}(x_1)},\frac{|x_{i2}-\text{median}(x_2)|}{\text{mad}(x_2)}\right).$$ 
is smaller than 3. Therefore, in this case, the Pearson correlation computed on the winsorized observations would be identical to the Pearson correlation computed on the original data.
Now, 0.5 is very far from the correlation of the good part of the data (e.g. considering only the black dots) which in this case is 84%. 
Contrast this with the results obtained by carrying a bivariate trimming and estimating the correlation of the remaining (untrimmed) observations. In this  case, the multivariate trimming was done  using the FastMCD(1) algorithm, perhaps the most popular multivariate trimming algorithm. FastMCD correctly identifies the red dots as being too far from the ellipse enclosing the majority of the data and flags them as outliers. Then, the correlation estimated on the remaining observations is now 85%, which is close enough to the correct result.

(1) Rousseeuw P. J. and Van Driessen K. (1999). A Fast Algorithm for the
  Minimum Covariance Determinant Estimator. Technometrics, 41, 212--223.

library(MASS)
library(rrcov)
n<-100
p<-2
set.seed(123)
A<-matrix(rnorm((p+1)*p),p+1,p)
A<-eigen(var(A))$vector
B<-A%*%diag(c(16,1))%*%t(A)
C<-t(A)%*%diag(c(4,1))%*%A
x<-mvrnorm(n,rep(0,p),B)
y<-mvrnorm(n,rep(0,p),C)
d<-which.max(mahalanobis(y,rep(0,p),B))
y<-mvrnorm(floor(n/5),y[d,],diag(2)/100)
z<-rbind(x,y)

plot(z,asp=1,type="n")
points(x,col="black",pch=16)
points(y,col="red",pch=16)

cor(z)

          [,1]      [,2]
[1,] 1.0000000 0.5018708
[2,] 0.5018708 1.0000000

cov2cor(CovMcd(z)@cov)
         [,1]      [,2]
[1,] 1.0000000 0.8592597
[2,] 0.8592597 1.0000000

cor(x)
          [,1]      [,2]
[1,] 1.0000000 0.8485822
[2,] 0.8485822 1.0000000

d1<-which((abs(z[,1]-median(z[,1]))/mad(z[,1])<3) & (abs(z[,2]-median(z[,2]))/mad(z[,2])<3))
length(d1)
[1] 120

