With respect to the double integral
$4 \iint I(x^2 + y^2 < 1) P(x, y) dx~dy$,
we can easily see that this evaluates to $\pi$ if we actually perform the calculation. First notice that the function P is 1 only on the inside of the unit square in the upper right quadrant. Consequently, we can rewrite the integral as $4 \int_0^1 \int_0^1 I(x^2 + y^2 < 1) dx~dy$.Now we notice that the indicator function is 1 when ever the coordinate (x, y) is inside the unit circle. This means that multiplying the two values will yield 1 when ever (x, y) is in the upper right quadrant of the unit circle, and will be 0 everywhere else. We can therefore transform our double integral into
$4 \int_0^1 \int_0^{\sqrt{1 - y^2}} 1~dx~dy = 4 \int_0^1 \sqrt{1-y^2} dy = 4 \frac{\pi}{4} = \pi$.
If we think carefully about what the indicator and P are doing, intuitively it should be obvious what is going on. These two functions will be 1 in the red area of your diagram. If we think about what the area of a quarter of a unit circle is, it should be obvious that it is $\pi/4$ since $\pi~r^2 = \pi~1^2 = \pi$. This integral is just a round about way of calculating the area of a quarter of the unit circle.
With respect to the octave code, if we go through it piece by piece everything should be clear.
First we have a=rand(S, 2). All this does is generate S samples from the unit square in the upper right quadrant. In other words, it generates samples from the square shown in his diagram. Note that P will evaluate to 1 for all the S samples in a.
Next we have this mysterious expression a.*a. All this does is multiply each entry in the matrix by itself. The sum with the 2 says to add each row in a, so in effect, sum(a.*a, 2) says compute $x^2 + y^2$ for each of my S samples. The < operator will return a 1 if true or a 0 if false, so it basically serves like our indicator function in the integral. In other words, the result of sum(a.*a, 2) < 1 will be a bunch of 1's and 0's, where a 1 indicates that a sample was inside the upper right quadrant of the unit circle and a 0 indicates that the sample was outside the unit circle (but still inside the unit square of course).
The key observation is that the ratio of 1's to 0's in this sampling scheme will be identical to the integral $\int_0^1 \int_0^1 I(x^2 + y^2 < 1) dx~dy$. Thus all we have to do is take the mean of the 1's and 0's, which is exactly the ratio we want. Finally we multiply by 4, and we have $\pi$. Thus the code 4*mean(sum(a.*a, 2)<1) really does approximate the integral we care about, which as we just showed, is exactly the value of $\pi$.