I am trying to understand the following code from ImageJ: http://pastebin.com/tXfhNxqf
The problem:
When computing the gaussian kernel we use the gaussian function
$$ f(x) = e^{-\dfrac{x^2}{2 \sigma^2}} $$
for a given $\sigma$, and it is common to take a radius $R = R(\sigma)$ and compute only the components until that $R$ (in order to make computations finite). In the code $R$ is computed in the line 27.
If we take as the new function $\hat{f}(x) = f(x), \, |x| < R$ and $\hat{f}(x) = 0$ elsewhere, we have a step at $x = R$ which leads to problems. In order to fix this, a common aproach is to replace the near-edge values by a 2nd order polynomial $p(x)$ which is $0$ at $x = R$ and satisfy $p(r) = f(r), \, p'(r) = f'(r)$ for some $0 < r < R$.
(This is also explained in the source code, as a comment).
How to find this polynomial?
My solution:
Let $p(x) = a x^2 + bx + c $ , then we have to solve the system:
$$ \begin{cases} p(r) = f(r) \\ p'(r) = f'(r) \\ p(R) = 0 \\ p'(R) = 0 \end{cases}$$
for the unknowns $a,b,c,r.$
What I think the code does:
Let $p(x) = a^2 (x-R)^2$, this already satisfy $p(R)=0, \,p'(R) = 0$. Then we have $p(r) = f(r)$ which implies $a^2 (x-R)^2 = f(r)$ solving for $a$, $a(r) = \frac{\sqrt{f(r)}}{R-r}$ (line 38, where kernel[0][r] is $f(r)$ and kRadius is $R$).
Then it takes r such that a(r) is minimum.
How is this a solution to the original problem?
The correct formula for $a(r)$ is $$ a(r)=\frac{f(r)^{\tfrac{1}{2}}}{R-r} $$ (note that $R > r$). Now $a(r)$ minimal implies $a'(r)=0$. Computing the derivative results in $$ a'(r) = a(r)\left(\frac{f'(r)}{2f(r)}+\frac{1}{R-r}\right).$$ So $a(r)$ minimal implies that $$ f'(r)=\frac{-2f(r)}{R-r}. $$ The derivative of $p$ is $p'(x) = -2a(r)^2(R-x)$ and setting $x= r$ results in $$ p'(r)=\frac{-2f(r)}{R-r}=f'(r) $$ so the derivatives of $p$ and $f$ coincide at $r$ as required.