I provide here one possible answer to your question. I had to make a few assumptions, but based on your comments, I believe they are acceptable.
I will describe briefly the assumptions, the algorithm and its properties, and will leave the path for the full mathematical derivation for assertions.
For easy of writing, given $m \in \, \mathbb{N}$, consider the notation $[m] = \{1,\ldots,m\}$.
The assumptions were:
We are sampling integer rectangles contained on the rectangle $[W]\times[H]$. This seems okay since you said you are working with images;
The dimensions $W$ and $H$ are powers of $2$. This seems less feasible, but I believe you could scale your images so that this assumptions is satisfied. On the full derivation, it becomes clear why this assumptions is required.
To get an algorithm to sample rectangles, we actually just need an algorithm to sample intervals from an integer interval $[2^{n}]$. Fix $\alpha \, \in \, (0,1)$, I will now describe a "pseudocode" for $H_\alpha(n)$, a recursive method that receives an integer $n$ and returns $[l, r]$, the left and right index of integer interval from $[2^{n}]$:
If $n = 0$, return $[1, 1]$. If not, go to step (2);
With probability $\alpha$, return $[1, 2^n]$. Else, go to step (3);
Call $H_\alpha(n-1)$ recursively and store its output in $[l_0, r_0]$. With probability $1/2$, return $[l_0, r_0]$. Else, return $[l_0 + 2^{n-1}, r_0 + 2^{n-1}]$.
Note: In the algorithm above, assume independence between all random variables used to make decisions where to halt or not.
Let $W = 2^n$ and $H = 2^m$. To sample a rectangle contained in $[W]\times[H]$, let $(x_1, x_2) = H_\alpha(n)$ and $(y_1, y_2) = H_\alpha(m)$, and assume that $(x_1, x_2)$ and $(y_1, y_2)$ are independent. Your rectangle is $R = [x_1, x_2]\times[y_1, y_2]$.
Here is a python implementation for you:
from random import random
def sample_interval(k):
alpha = 0.4 # Can be any number in (0,1)
if k == 0:
return (0, 0)
# with prob alpha return the whole interval
if random() < alpha:
return (0, 2**k - 1)
l0, r0 = sample_interval(k-1)
# half the times, do not add 2**(k-1)
if random() < 1.0/2:
return (l0, r0)
# half the times, do add 2**(k-1)
else:
return (l0 + 2**(k-1), r0 + 2**(k-1))
# remember that W = 2**n and H = 2**m
def sample_rectangle(n, m):
x1, x2 = sample_interval(n)
y1, y2 = sample_interval(m)
upper = (x1, y1)
lower = (x2, y2)
return (upper, lower)
Now you are surely asking: why sample it like this? The reason is that this algorithm has the properties you want:
(P1) Every integer point (pixel) on $[W]\times[H]$ has the same probability of appearing in the sampled rectangle. This probability is constant for all points, but depends on $n$ and $m$. Asymptotically, if $W$ and $H$ go to infinite, this probability converges to $(\frac{2\alpha}{1+\alpha})^2$;
(P2) The limit average rectangle area divided by the area of the big rectangle is scale free. That is, let $A_{nm}$ be the random variable associated with area a rectangle sampled from $[2^n]\times[2^m]$ according to this algorithm, then
\begin{equation*} a_s = \lim_{n\rightarrow\infty}\lim_{m\rightarrow\infty} \frac{E[A_{nm}]}{W H} = \left(\frac{2\alpha}{1+\alpha}\right)^2 \quad. \end{equation*}
Since these are the properties you asked, this sampling method solves your problem. Note that you can change $\alpha$ freely. Increasing $\alpha$ makes you sample the whole rectangle more. Decreasing it makes you sample smaller rectangles more. If $\alpha$ is very close to $0$, you will probably end up sampling lone pixels.
Mathematical Proofs The proofs are based on mathematical induction since the algorithm is inherently recursive. We also focus on the proof on the interval case. For the rectangles, it follows trivially from the interval case and independence.
Lets prove all points have equal probability of appearing in the random interval. For $n = 0$, it is trivial. Assume that all points in $[2^{n-1}]$ have constant probability $p_{n-1}$ to appear in a random interval (this obviously depends on $n$, but it can not depend on the point).
Now we prove for $k = n$. Consider a point in $[2^k]$ that is smaller or equal to $2^{k-1}$. There are two ways that a random interval can contain this point. First, in the case we sample the whole interval $[2^k]$, which has probability $\alpha$. Second, if, after not sampling in the first step, the interval from recursive call contains the point AND we do not add $2^{k-1}$. This second event has probability $(1-\alpha)\times p_{n-1}\times\frac{1}{2}$ (just use independence and the induction hypothesis). Therefore
$$ p_n = \alpha + \frac{(1-\alpha)}{2}p_{n-1} \quad.$$
You can prove that the probability is the same if you consider a point above $2^{k-1}$. This finishes the proof.
For the average size of the interval, writing $E_n$ as the expected interval size, the same idea gives you the recursion
$$ E_n = \alpha 2^n + (1-\alpha)E_{n-1} \quad.$$
Now divide both sides by $2^n$, we have
$$ \frac{E_n}{2^n} = \alpha + \frac{1-\alpha}{2}\frac{E_{n-1}}{2^{n-1}} \quad.$$
Taking the limit and calling it $e_s$, we have
$$ e_s = \alpha + \frac{1-\alpha}{2}e_s \quad,$$
and a bit of algebra gives you $e_s = \frac{2\alpha}{1+\alpha}$. Since the average area of the rectangle is the product of the average size (again, independence), it follows that
$$ a_s = e_s^2 \quad,$$
giving you your scale free algorithm.