Misunderstanding of Monte Carlo Pi Estimation I am fairly sure that I understand the how Monte Carlo integration works but I am not understanding the formulation of how it is used to estimate Pi.  I am going by the procedure outlined in the 5th slide of this presentation http://homepages.inf.ed.ac.uk/imurray2/teaching/09mlss/slides.pdf
I understand the preliminary steps.  Pi is equal to 4 times the area of a quarter of the unit circle.  And the area of the top-right quarter of the unit circle centered at (0,0) is equivalent to the integral of the curve that is the top-right quarter of the unit circle in $0<x<1$ and $0<y<1$.  
What I don't understand is how this integral is
$\iint I((x^2+y^2)<1)P(x,y)dxdy$
where $P(x,y)$ is uniformly distributed in the unit square around the quarter circle (i.e. it is always equal to 1 if $0<x<1$ and $0<y<1$ and 0 otherwise).  So this would mean that
 $I((x^2+y^2)<1)P(x,y)$
is the function that is the top-right quadrant of the unit circle at $0<x<1$ and $0<y<1$ but I do not understand how this is true since the indicator function can only be 1 or 0.  I understand that it is probably written in this way to make Monte Carlo sampling easy (i.e. it is an expectation so just sample from $P(x,y)$ and get the average of the samples applied to $I((x^2+y^2)<1)$) but it just does not make intuitive sense to me why that integral represents the area under that curve.
Could someone provide an intuitive explanation of this.  Maybe show how that integral was derived in a step-by-step way?
EDIT:
I was able to gain a better understanding by relating the expectation to an area.  I will explain it here in case it helps anyone.  First start with relating Pi to the area of the top-right quadrant of the unit circle
$\pi=4\times A_{tr}$
Then we place the top-right quadrant into the unit square.  And under a uniform distribution over the unit square, the area of the circle quadrant is proportional to the probability of obtaining a sample from it.  It follows that the following equality holds
$P(x^2+y^2<1)=\frac{A_{tr}}{A_{square}}$
and $A_{square}=1$ so
$P(x^2+y^2<1)=A_{tr}$
And substituting into the original equation
$\pi=4\times P(x^2+y^2<1)$
and it is also true that $P(x^2+y^2<1)=E[I(x^2+y^2<1)]$ which is equal to the original double integral.
So I understood it by relating the area to a probability then relating that probability to an expectation that is equivalent to the integral.  Let me know if I have made any mistakes.
 A: The area of a  circle circle of radius $l$ is equal to $\pi l^2$. It means that a quarter of circle has area $l^2\pi/4$. This means that the square with side the radius of the circle as $area=l^2$. 
This means that the ratio between the area of a quarter of circle and the area of the square is $\pi/4$.  
A point  $(x,y) $ is in the square if $ 0<x<1, 0<y<1$. 
and it is in the quarter of circle if  $ 0<x<1, 0<y<1 ,x^2+y^2<1$.  
Your integral is so $∬I((x^2+y^2)<1)P(x,y)= ∬I((x^2+y^2)<1) I(0<x<1)I(0<y<1)$ That is exactly the area described by a quarter of circle

A: The simplest intuitive explanation relies on understanding that $E(I(A)) = P(A)$.  Thus, $\int \int I(x^2+y^2 < 1)dxdy = P(x^2 + y^2 < 1)$.  Once you realize the double integal is simply a probability, it should make intuitive sense that you could sample $x$ and $y$ from the unit square and compute the proportion of draws for which $x^2 + y^2 <1$.   
Perhaps the other piece of intuition missing from your understanding is the connection between area and probability.  Since the area of the entire unit square is 1 and points $(x,y)$ are uniformly distributed within the square, the area of any region $A$ within the unit square would correspond to the probability that a randomly chosen point would be within $A$.
A: I landed on this surfing CV, and I see that the code of the Monte Carlo is in Octave. I happen to have a simulation in R that makes the idea of deriving the number $\pi$ as a bivariate uniform distribution in the $[0,1]$ plane under the constraints of the integrals in the OP very intuitive:
Given that the quarter of a circle is enclosed in a 1-unit square, the area is $\pi/4$. So generating uniformly distributed points in the square $(x,y)$ will end up carpeting the entire square, and calculating the fraction fulfilling $1 < \sqrt{(x^2+y^2)}$ will be tantamount to integrating $∬\textbf{1}((x^2+y^2)<1) \,\textbf{1}(0<x<1)\,\textbf{1}(0<y<1)$ since we are just selecting the fraction of dots within the circle in relation to the unit square:
x <- runif(1e4); y <- runif(1e4)
radius <- sqrt(x^2 + y^2)
# Selecting those values within the circle is obtained with radius[radius < 1]:
(pi = length(radius[radius < 1]) / length(radius)) * 4     =    3.1272

We can plot the values falling within the radius among 10,000 draws:

And we can, naturally, get closer and closer approximation by selecting more points. With 1 million points we get:
(pi = length(radius[radius < 1]) / length(radius)) * 4 [1] 3.141644
a very approximate result. Here is the plot:

