There's a relation between two-dimensional Brownian motion and conformal maps, see e.g. Thurston's answer to this question. Given two non-empty simply-connected domains $U$ and $V$ in the complex plane which are not equal to the complex plane, the Riemann mapping theorem guarantees the existence of a conformal bijection $U\to V$. Is there a way to use Brownian motion to construct such a map?

Let $U$ be a Jordan domain and $b_0,b_1,b_2$ be points appearing in order on $\partial U$. Let $f : U \to \mathbb{H}$ be the Riemann map normalized so that $f(b_0)=0$, $f(b_1)=1$ and $f(b_2)=\infty$. Chris Bishop taught me a poor man's way to construct $f$ using Brownian motion.

Let $A_k$ be the arc of $\partial U$ between $b_k$ and $b_{k+1}$, where the index is taken modulo 3. Similarly, let $B_0$ denote the interval $(0,1)$, $B_1$ the interval $(1,+\infty)$ and $B_2$ the interval $(-\infty,0)$, so that $f$ maps $A_k$ to $B_k$ for each $k$. For $z\in U$, let $p_k$ denote the probability that Brownian motion started at $z$ first hits $\partial U$ along $A_k$. The probability $p_k=\omega(z,U)(A_k)$ is known as the harmonic measure of $A_k$ with respect to $(z,U)$ and is invariant under conformal maps. One can show there is a unique point $w \in \mathbb{H}$ such that $\omega(w,\mathbb{H})(B_k)=p_k$ for each $k$. By conformal invariance, we have $w=f(z)$.

Here is how to find $w$ given $(p_0,p_1,p_2)$. If $C_0,C_1,C_2$ are consecutive arcs on the unit circle such that $C_k$ measures $2\pi p_k$ radians, then the harmonic measure $\omega(0,\mathbb{D})(C_k)$ is obviously $p_k$. Now, let $c_k$ be the point at the junction between $C_{k-1}$ and $C_k$ and let $M$ be the Moebius transformation sending $c_0,c_1,c_2$ to $0,1,\infty$ respectively. You can write $M$ explicitly as $$ M(\zeta) = \frac{(c_1-c_2)(\zeta-c_0)}{(c_1-c_0)(\zeta-c_2)}. $$ It is easy to see that $M$ maps the unit disk to the upper half-plane. By conformal invariance, we have $$\omega(M(0),\mathbb{H})(B_k)=\omega(0,\mathbb{D})(C_k)=p_k$$ for each $k$ and hence $M(0)=w=f(z)$.

I implemented this in Matlab some years ago. What you do is for each $z$ where you want to evaluate the Riemann map, throw a bunch of Brownian motions to estimate the probabilities $p_0$,$p_1$,$p_2$. It was not very efficient, but was quite easy to program.

I believe that a way of doing this is discussed in Chapter 5 of Richard Bass' "Probabilistic Techniques in Analysis." If I remember correctly, the main place Brownian motion is used is to construct the Green function for the domain.

This chapter of this book discusses how the proofs of many theorems in complex analysis can be proven with Brownian motion... Picards theorem and Carleson's corona theorem in particular. Most of the proofs stem from Thurston's observation that Brownian motion is conformally invariant, but it is a good read for those interested.