$
\def \J {\Delta}
\def \half {\frac{1}{2}}
\def \kwart {\frac{1}{4}}
\def \hieruit {\quad \Longrightarrow \quad}
\def \slechts {\quad \Longleftrightarrow \quad}
\def \norm {\frac{1}{\sigma \sqrt{2\pi}} \; }
\def \EN {\quad \mbox{and} \quad}
\def \OF {\quad \mbox{or} \quad}
\def \wit {\quad \mbox{;} \quad}
\def \spekhaken {\iint}
\newcommand{\dq}[2]{\displaystyle \frac{\partial #1}{\partial #2}}
\newcommand{\oq}[2]{\partial #1 / \partial #2}
\newcommand{\qq}[3]{\frac{\partial^2 #1}{{\partial #2}{\partial #3}}}
\def \erf {\operatorname{Erf}}
$
Sharpened Points and Contours
A point in the plane is an idealized pixel, without size. A pixel is a fuzzyfied
point, smeared out over a small domain in the plane. Cast in more mathematical
terms: a point at $(p,q)$ is a delta function at that place. The fuzzyfication
of that point, the pixel, can be modelled by a 2-D bell shaped function, which
is centered at $(p,q)$:
\begin{eqnarray*}
\mbox{point} &=& \delta(x-p,y-q) \quad = \quad \delta(x-p)\delta(y-q) \\
\mbox{fuzzy} &=& \frac{1}{2\pi\sigma^2}\:
e^{-\half\left[(x-p)^2+(y-q)^2\right]/\sigma^2} \\
&=& \norm e^{-\half(x-p)^2/\sigma^2}
\norm e^{-\half(y-q)^2/\sigma^2}
\end{eqnarray*}
The point can be fuzzyfied to a pixel and the pixel can be sharpened to a point.
The latter is accomplished by a limit where the spread called $\sigma$
becomes zero:
$$
\lim_{\sigma\rightarrow 0}
\frac{1}{2\pi\sigma^2}\:e^{-\half\left[(x-p)^2+(y-q)^2\right]/\sigma^2} =
$$ $$
\lim_{\sigma\rightarrow 0}
\norm e^{-\half(x-p)^2/\sigma^2}
\lim_{\sigma\rightarrow 0}
\norm e^{-\half(y-q)^2/\sigma^2}
$$ $$
= \delta(x-p)\delta(y-q) = \delta(x-p,y-q)
$$
The above becomes somewhat less trivial as soon as (the fuzzyfication of) a
point is integrated over a certain area, which is enclosed by a contour curve:
$$
\spekhaken \delta(x-p,y-q) \: dx\,dy
$$
The outcome of this integral depends on whether the point $(p,q)$
is inside or outside the contour curve which encloses the area.
In the former case it is obviously one, in the latter case it is zero:
$$
\spekhaken \delta(x-p,y-q) \: dx\,dy
= \left\{ \begin{array}{ll} 1 \\ 0 \end{array} \right.
\mbox{if} \; (p,q) \begin{array}{cc} \mbox{inside } \\
\mbox{outside} \end{array}
$$
If the delta function is "cut into pieces", which happens if $(p,q)$ is exactly
at the boundary curve, then the outcome of the integral is a undefined
value, though somewhere between $0$ and $1$ .
But let's discard this special case for the moment being. And consider instead
the fuzzyfication of the above:
$$
\spekhaken
\frac{1}{2\pi\sigma^2}\:e^{-\half\left[(x-p)^2+(y-q)^2\right]/\sigma^2}
\: dx\,dy
$$ $$
= \spekhaken
\norm e^{-\half(x-p)^2/\sigma^2}
\norm e^{-\half(y-q)^2/\sigma^2} \: dx\,dy
$$
According to Green's Theorem, any such area integral may be converted into a
contour integral. As follows:
$$
\spekhaken \left( \dq{Q}{x} - \dq{P}{y} \right) \, dx \, dy =
\oint \left( P\,dx + Q\,dy \right)
$$
Where the line integral must be evaluated along the contour enclosing the area
of interest. Converting a line integral into an area integral is easy. But the
reverse, converting an area integral into a line integral, requires some more
ingenuity. In our case, though, there exists more than one possibility to do it.
We opt for the following. Let the function $P(x,y)$ be zero everywhere and let
the function $Q(x,y)$ be defined by:
$$
Q(x,y) = \erf\left(\frac{x-p}{\sigma}\right) \;
\norm e^{-\half(y-q)^2/\sigma^2} \hieruit
$$ $$
\dq{Q}{x} = \norm e^{-\half(x-p)^2/\sigma^2}
\norm e^{-\half(y-q)^2/\sigma^2}
$$
Which is exactly what we want. Now the result of Green's Theorem becomes:
$$
\spekhaken \dq{Q}{x} \,dx\,dy = \oint Q\,dy = \oint
\erf\left(\frac{x-p}{\sigma}\right) \norm e^{-\half(y-q)^2/\sigma^2} \,dy
$$
In the limiting case, when $\sigma \rightarrow 0$ , the exponential function
becomes a delta function and the $\erf$ function becomes an integrated delta
function, namely the Heaviside function. Thus we find:
$$
\lim_{\sigma\rightarrow 0} \oint
\erf\left(\frac{x-p}{\sigma}\right) \norm e^{-\half(y-q)^2/\sigma^2} \,dy
= \oint H(x-p) \, \delta(y-q) \, dy
$$
But there are a myriad other ways to look at the integral:
$$
\spekhaken
\frac{1}{2\pi\sigma^2}\:e^{-\half\left[(x-p)^2+(y-q)^2\right]/\sigma^2}
\: dx\,dy
$$
For example, take the origin of your coordinates at $(p,q) = (0,0)$ and replace
the Cartesian $(x,y)$ system by polar coordinates $(r,\phi)$ :
$$
\spekhaken
\frac{1}{2\pi\sigma^2}\: e^{-\half\left[x^2+y^2\right]/\sigma^2}
\: dx\,dy
= \spekhaken
\frac{1}{2\pi\sigma^2}\: e^{-\half r^2/\sigma^2}
\: r\,dr\,d\phi
$$
The latter integral can be evaluated as follows:
$$
\spekhaken
\frac{1}{2\pi\sigma^2}\: e^{-\half r^2/\sigma^2}
\: r\,dr\,d\phi =
\frac{1}{2\pi} \int_{\phi} \left[ - \int_r
e^{-\half (r/\sigma)^2} d\left\{-\half(r/\sigma)^2\right\} \right] d\phi
$$
If we take the limit for $\sigma \rightarrow 0$ of the integral between square
brackets, then it is reduced to unity:
$$
\lim_{\sigma \rightarrow 0} - \int_0^{-\half R^2/\sigma^2} e^t \, dt =
\lim_{\sigma \rightarrow 0} \left[ 1 - e^{-\half (R/\sigma)^2} \right] = 1
$$
Here $R$ denotes the distance of the boundary curve to the origin. If it may be
assumed that the origin is distinct from the boundary, then for $\sigma
\rightarrow 0$ the quotient $R/\sigma$ approaches infinity and the exponential
function becomes zero indeed. So we are left with:
$$
\spekhaken
\frac{1}{2\pi\sigma^2}\: e^{-\half r^2/\sigma^2}
\: r\,dr\,d\phi = \frac{1}{2\pi} \oint d\phi
$$
Where the integral at the right is recognized as the winding number of
the curve enclosing the area of interest. This leads to the following
concise conclusion: \large
$$ \overline{\underline{\left|
\oint H(x) \, \delta(y) \, dy = \frac{1}{2\pi} \oint d\phi
\right| }} $$ \normalsize
Where the point of interest, formerly called $(p,q)$ , is identified with the
origin $(0,0)$ on a permanent basis. I feel that this conclusion is valid for
an arbitrary bunch of mutually and eventually self-intersecting closed
curves. With the sole condition that the origin is not coincident with any place
on these.
The rest of the story is technology, not a theory. I want to spend
a few words on it, though. The integral with $H(x)$ and $\delta(y)$ in it
is actually equivalent with the following well-known statement. Draw a line
from the point of interest $P$, in the positive x-direction, towards infinity.
Then count the number of intersections of that line with the curves that form
the boundary of the area of interest $A$.
If the number of intersections is odd then $P$ is inside $A$.
If the number of intersections is even then $P$ is outside $A$.
However, the technical implementation of this is far more tricky that one might
expect. Special cases like being smaller $(<)$, equal $(=)$ or greater $(>)$
than zero should be distinguished carefully, eventually resulting in
$3^2 \times 3^2 = 81$
different combinations (!). It should be mentioned, as well, that the Heaviside
function $H(x)$ algorithmically means that the infamous clipping problem
(see: Computer Graphics) has an instance at $x=0$. Last but not least, if you
really want to improve efficiency, then a sorting and searching problem
with respect to the $y$-coordinates will become part of the project.