Continuous functions on $\mathbb{R}^2$ with special property The following problem is from Miklos Schweitzer competition (Year 1983, Problem 7):

Prove that if the function $f: \mathbb{R}^{2}\to [0, 1]$ is
  continuous, and its average on every circle of radius one equals the
  function value at the center, then $f$ is constant.

There is some discussion here but no elementary solution (from my point of view). I would like to see a solution that avoids Fourier transforms. But perhaps someone could elaborate on Kent Merryfield's observations and give a complex-analytic solution?
From the number of upvotes this meta question has received, I will assume that its okay to ask such questions here.
 A: A proof from the martingale convergence theorem.  (I got this from Persi Diaconis back when we were in graduate school.)  
Consider i.i.d. random variables $X_n$, distributed according to normalized arc length on the unit circle in $\mathbb R^2$.    These are defined on some sample space $(\Omega,\mathcal F, \mathbb P)$.  Define a stochastic basis $(\mathcal F_n)$, defining $\mathcal F_n \subset \mathcal F$ to be the sigma-field generated by $X_1,\cdots,X_n$.  Define $S_n = X_1+\cdots +X_n$, a random walk in $\mathbb R^2$ suggested by this problem.  
Almost surely, the sequence $(S_n)$ is dense in the plane.  Such random walks in $2$ dimensions are recurrent: another fact needed here.  (The probability of returning to a small neighborhood of the origin after $n$ steps is asymptotically proportional to $1/n$ and the series $\sum 1/n$ diverges, so by Borel-Cantelli, we return infinitely often.)
Now suppose we have our continuous function $f : \mathbb R^2 \to [0,1]$ with the averaging property.  Consider the real stochastic process
$$
Y_n = f(S_n)
$$
The conditional expectation obeys
$$
\mathbb E[Y_{n+1} | \mathcal F_n] = Y_n .
$$
This is exactly the averaging property, and the reason $X_n$ was defined as it was.  Thus: $(Y_n)$ is a martingale.  The values are in $[0,1]$ so it is a bounded martingale.   Therefore by the martingale convergence theorem, $Y_n$ converges a.s.  But, with probability one, $S_n$ is dense in the plane.  So (recall $f$ is continuous) the only way $f(S_n)$ can converge is for $f$ to be constant in the plane.
A: An almost-proof from the central limit theorem. This write up took far more words than I expected. The point is that $f(x_0, y_0)$ is the average of $f$ on a circle centered around $(x_0, y_0)$. The value of $f$ at each point on that circle is, in turn, the average of $f$ at circles centered round the points of the unit circle. This means that $f(x_0,y_0)$ can be written as some sort of two dimensional average of the values of $f$ on a disc of radius $2$ around $(x_0, y_0)$. Repeating this argument, $f$ can be written as a weighted average of the values of $f$ on a disc of radius $n$ and, by the central limit theorem, as $n \to \infty$, the weighting function will start to look like a bell curve of width $\approx \sqrt{n}$. So $f(x_1,y_1) - f(x_2,y_2)$ is the integral of $f$ against a difference of two bell curves, whose centers are displaced from each other by $(x_1-x_2, y_1-y_2)$. As $n \to \infty$, the two bell curves will get very wide and flat, so the difference between their centers becomes insignificant and the integral goes to zero. Okay, now for the proof.

The hypothesis is that
$$f(x,y) = \frac{1}{2 \pi} \int_{\theta} f(x+\cos \theta, y+\sin \theta) d \theta \quad (\ast)$$
Plugging formula $(\ast)$ into itself,
$$f(x,y) = \frac{1}{4 \pi^2} \int_{\theta} \int_{\phi} f(x+\cos \theta+\cos \phi, y+\sin \theta + \sin \phi) d \theta d \phi$$
which, if I got the Jacobian factors right, turns into 
$$f(x,y) = \frac{1}{2 \pi^2} \int_{u,v} f(x+u,y+v) \sqrt{(u^2+v^2)(2-u^2-v^2)} du dv$$
I don't care about the details of this formula, what I care is that it looks like 
$$f(x,y) = \int_{u,v} f(x+u, y+v) \mu_2(u,v) du dv$$
for some function $\mu_2$. 
If I plug $(\ast)$ into itself to get a triple integral, I will get some different formula
$$f(x,y) = \int_{u,v} f(x+u, y+v) \mu_3(u,v) du dv$$
and, in general, we have
$$f(x,y) = \int_{u,v} f(x+u, y+v) \mu_n(u,v) du dv$$
where $\mu_n$ is the $n$-fold convolution of the measure of the unit circle with itself. If you are uncomfortable convolving measures and want to convolve functions, take $n$ even and define $\mu_n = \mu_2 \ast \mu_2 \ast \cdots \ast \mu_2$, where I have convolved $n/2$ times.
By the central limit theorem in $\mathbb{R}^2$, 
$$\mu_n(x,  y) \approx \frac{C}{n} \exp\left(-\frac{x^2+y^2}{n \sigma} \right)$$
Here $C$ and $\sigma$ are positive real numbers I didn't feel like computing. This proof will be nonrigorous to the extent that I don't replace that $\approx$ by some precise limiting notion; I was too lazy to scan through the different notions of convergence for which the central limit theorem has been proved and figure out which one I want.
Now, fix $(x_1, y_1)$ and $(x_2, y_2) \in \mathbb{R}^2$. So
$$f(x_1, y_1) - f(x_2, y_2) =$$
$$ \int_{(u,v) \in \mathbb{R}^2} f(x_1+u, y_1+v) \mu_n(u,v) dx dy - \int_{(u,v) \in \mathbb{R}^2} f(x_2+u, y_2+v) \mu_n(u,v) dx dy=$$
$$\int_{(a,b) \in \mathbb{R}^2} f(a,b) \left( \mu_n(a-x_1, b-y_1) - \mu_n(a-x_2, b-y_2) \right) da db \approx$$
$$ \int_{(a,b) \in \mathbb{R}^2} f(a,b) \frac{C}{n} \left( \exp\left( - \frac{(a-x_1)^2+(b-y_1)^2}{\sigma n} \right) - \exp\left( - \frac{(a-x_2)^2+(b-y_2)^2}{\sigma n} \right) \right) da db.$$
Acting as if the $\approx$ were an $=$, we have
$$|f(x_1,y_1) - f(x_2, y_2)| \leq$$
$$\frac{C}{n} \int_{(a,b) \in \mathbb{R}^2} \left| \exp\left( - \frac{(a-x_1)^2+(b-y_1)^2}{\sigma n} \right) - \exp\left( - \frac{(a-x_2)^2+(b-y_2)^2}{\sigma n} \right) \right| da db.$$
We will show that
$$\lim_{n \to \infty} \frac{C}{n} \int_{(a,b) \in \mathbb{R}^2} \left| \exp\left( - \frac{(a-x_1)^2+(b-y_1)^2}{\sigma n} \right) - \exp\left( - \frac{(a-x_2)^2+(b-y_2)^2}{\sigma n} \right) \right| da db = 0$$
and, thus, once we get the right form of the central limit theorem to cite, $f$ is constant.
So it remains to compute 
$$\lim_{n \to \infty} \frac{C}{n} \int_{(a,b) \in \mathbb{R}^2} \left| \exp\left( - \frac{(a-x_1)^2+(b-y_1)^2}{\sigma n} \right) - \exp\left( - \frac{(a-x_2)^2+(b-y_2)^2}{\sigma n} \right) \right| da db $$
By translational and rotational symmetry, we may assume that $(x_1, y_1)$ and $(x_2, y_2)$ are of the form $(r,0)$ and $(-r,0)$ for some positive $r$.
$$\frac{C}{n} \int_{(a,b) \in \mathbb{R}^2} \left| \exp\left( - \frac{(a-r)^2+b^2}{\sigma n} \right) - \exp\left( - \frac{(a+r)^2+b^2}{\sigma n} \right) \right| da db =$$
$$\frac{C}{n} 2 \int_{b=-\infty}^{\infty} \int_{a=0}^{\infty}  \exp\left( - \frac{(a-r)^2+b^2}{\sigma n} \right) - \exp\left( - \frac{(a+r)^2+b^2}{\sigma n} \right)  da db =$$
$$\frac{2C}{n} \int_{b=-\infty}^{\infty} \int_{a=0}^{2r}  \exp\left( - \frac{(a-r)^2+b^2}{\sigma n} \right) da db = \frac{2C}{n} \int_{c=-r}^r \exp \frac{-c^2}{\sigma n} dc \int_{b=-\infty}^{\infty} \exp\left( \frac{-b^2}{\sigma n} \right) db$$
$$=\frac{2 C \sqrt{\pi \sigma}}{\sqrt{n}} \int_{c=-r}^r \exp \frac{-c^2}{\sigma n} dc = 2 C \sqrt{\pi \sigma} \int_{x=-r/\sqrt{n}}^{r/\sqrt{n}} \exp\left(- \frac{x^2}{\sigma}  \right) dx.$$
(We have made the subsititutions $c=a-r$ and $x=c/\sqrt{n}$.) Clearly, the last integral goes to $0$.
