An edge can be defined as the change of color between two contiguous area of pixels, bigger the change of color, more noticeable is the edge. This definition perfectly match the definition of derivative:
$$f'_{x(x,y)}=\lim_{\Delta_x\to0}\frac{f_{(x+\Delta_x,y)}-f_{(x,y)}}{\Delta_x}$$
$$f'_{y(x,y)}=\lim_{\Delta_y\to0}\frac{f_{(x,y+\Delta_y)}-f_{(x,y)}}{\Delta_y}$$
Edge detection filters are in fact an application of the Derivative. Since we can treat an image as a discrete function, we will use an approximation and instead of $\Delta_x\to0$ and $\Delta_y\to0$ we will use $\Delta_x=1$ and $\Delta_y=1$, so now we have:
$$G_x=f_{(x+1,y)}-f_{(x,y)}$$
$$G_y=f_{(x,y+1)}-f_{(x,y)}$$
$G_x$ and $G_y$ are the partial derivatives of $f_{(x,y)}$ and the gradient and direction can be calculated as:
$$G=\sqrt{G_x^2+G_y^2}$$
$$dir=\arctan{\frac{G_y}{G_x}}$$
Sometimes $G=|G_y|+|G_x|$ is used as an approximation for the gradient, for fast calculation. With all the theory presented we can start hacking some code.
#include <iostream>
#include <cmath>
#include <QCoreApplication>
#include <QImage>
enum GradienType {
GradienTypeX,
GradienTypeY,
GradienTypeMod,
GradienTypeSum
};
GradienType gradientType = GradienTypeMod;
inline int gradient(int diffX, int diffY, GradienType gradientType)
{
switch (gradientType) {
case GradienTypeX:
return qAbs(diffX);
case GradienTypeY:
return qAbs(diffY);
case GradienTypeMod:
return sqrt(diffX * diffX + diffY * diffY);
case GradienTypeSum:
return qAbs(diffX) + qAbs(diffY);
default:
break;
}
return 0;
}
int main(int argc, char *argv[])
{
QCoreApplication a(argc, argv);
Q_UNUSED(a)
QImage inImage("lena.png");
inImage = inImage.convertToFormat(QImage::Format_RGB32);
QImage outImage(inImage.size(), inImage.format());
for (int y = 0; y < inImage.height(); y++) {
const QRgb *iLine = (const QRgb *) inImage.constScanLine(y);
const QRgb *iLinePrev = (const QRgb *) inImage.constScanLine(y - 1);
QRgb *oLine = (QRgb *) outImage.scanLine(y);
for (int x = 0; x < inImage.width(); x++) {
int diffX = qGray(iLine[x]);
int diffY = diffX;
if (x > 0)
diffX -= qGray(iLine[x - 1]);
if (y > 0)
diffY -= qGray(iLinePrev[x]);
int grad = gradient(diffX, diffY, gradientType);
quint8 c = qBound(0, grad, 255);
oLine[x] = qRgba(c, c, c, qAlpha(iLine[x]));
}
}
outImage.save("edge.png");
return EXIT_SUCCESS;
}
Also, I have included the option to see the gradient just in one of the directions. And here we can compare the results:
This filter is a good approach to edge detection, but suffer a serious problem, its very sensible to noise. This filter will mark garbage as if it were a real edge. But we are lucky because we had seen a trick to remove noise in an image. We can just apply a gaussian filter to the image and then the edge detection algorithm, but this will require us to apply two operations to the image, do you remember my post about image convolution? do you remember this kernel?:
int kw = 3;
int kh = 3;
qreal kernel[] = {-1, -1, -1,
-1, 8, -1,
-1, -1, -1};
This kernel produce the same result as we are finding. But, why this kernel works? why choose that combination of numbers? can we use a bigger kernel?, continue reading for the explanation.
One of the properties of the convolution is that the derivative of a convolution is equal to convolution between of one of the functions and the derivative of the other, this is:
$$\delta(f*g)=\delta(f)*g=f*\delta(g)$$
so, we can apply the derivative to the gaussian filter, and then apply the resulting filter to the image. Let say we have:
$$f(U)=k\exp^U$$
where:
$$k=\frac{1}{\sigma\sqrt{2\pi}}$$
$$U=-\frac{x^2+y^2}{2\sigma^2}$$
The first and second derivatives are:
$$ f'_{(U)}=kU'\exp^U \\ f''_{(U)}=k(U''+U'^2)\exp^U $$
The partial derivatives of U are:
$$ \begin{matrix} U'_x=-\frac{x}{\sigma^2} & U'_y=-\frac{y}{\sigma^2} \\ U''_x=-\frac{1}{\sigma^2} & U''_y=-\frac{1}{\sigma^2} \end{matrix} $$
so, replacing U and it's derivatives we have:
$$ \begin{matrix} \frac{\partial f_{(x,y)}}{\partial x}=-\frac{x}{\sigma^3\sqrt{2\pi}}\exp^{-\frac{x^2+y^2}{2\sigma^2}} & \frac{\partial f_{(x,y)}}{\partial y}=-\frac{y}{\sigma^3\sqrt{2\pi}}\exp^{-\frac{x^2+y^2}{2\sigma^2}} \\ \frac{\partial^2 f_{(x,y)}}{\partial x^2}=\frac{x^2-\sigma^2}{\sigma^5\sqrt{2\pi}}\exp^{-\frac{x^2+y^2}{2\sigma^2}} & \frac{\partial^2 f_{(x,y)}}{\partial y^2}=\frac{y^2-\sigma^2}{\sigma^5\sqrt{2\pi}}\exp^{-\frac{x^2+y^2}{2\sigma^2}} \end{matrix} $$
The first derivative is commonly known as Differential of Gaussian, while the second derivative is commonly known as Laplacian of Gaussian.
In the following graph you can see the gaussian filter, and it's first and second derivatives in $\mathbb{R}^2$:
inline QVector<qreal> edgeKernel(int radius,
qreal sigma,
qreal scaleXY,
qreal scaleW,
bool axysY=false,
bool round=false,
int *kl=NULL)
{
int kw = 2 * radius + 1;
QVector<qreal> kernel(kw * kw);
for (int j = 0; j < kw; j++)
for (int i = 0; i < kw; i++) {
qreal x = scaleXY * (i - radius);
qreal y = scaleXY * (j - radius);
qreal k = - scaleW * (axysY? y: x) / (pow(sigma, 3) * sqrt(2 * M_PI));
qreal sigma2 = -2 * sigma * sigma;
qreal weight = k * exp((x * x + y * y) / sigma2);
kernel[i + j * kw] = round? qRound(weight): weight;
}
*kl = kw;
return kernel;
}
Using this code and formulas, we can practically create any kernel in the form:$$ \begin{vmatrix} -a & 0 & a \\ -b & 0 & b \\ -a & 0 & a \end{vmatrix} $$
We can tweak the parameters to obtain the most common edge detection filters as for example the Sobel, Prewitt or Scharr operators, or we can even extend the kernel for bigger radius. Here are some examples:
| Operator | KernelX | KernelY | Parameters |
|---|---|---|---|
| Sobel | $$\begin{vmatrix}-1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1\end{vmatrix}$$ | $$\begin{vmatrix}1 & 2 & 1 \\ 0 & 0 & 0 \\ -1 & -2 & -1\end{vmatrix}$$ | radius = 1 sigma = 1 scaleXY = 1 scaleW = 10 round = true |
| Prewitt | $$\begin{vmatrix}-1 & 0 & 1 \\ -1 & 0 & 1 \\ -1 & 0 & 1\end{vmatrix}$$ | $$\begin{vmatrix}1 & 1 & 1 \\ 0 & 0 & 0 \\ -1 & -1 & -1\end{vmatrix}$$ | radius = 1 sigma = 2 scaleXY = 1 scaleW = 25 round = true |
| Scharr | $$\begin{vmatrix}-3 & 0 & 3 \\ -10 & 0 & 10 \\ -3 & 0 & 3\end{vmatrix}$$ | $$\begin{vmatrix}3 & 10 & 3 \\ 0 & 0 & 0 \\ -3 & -10 & -3\end{vmatrix}$$ | radius = 1 sigma = 0.6 scaleXY = 1 scaleW = 22 round = true |
You can try tweaking the parameters of the filter, for example, a bigger radius gives as result a stronger edge, but also remarks spurious pixels, bigger sigma values reduces spurious pixel but gives thiner edges, you can disable round for an accurate but slower result.
Here is an example for radius = 2, sigma = 10, scaleXY = 0.5, scaleW = 200 and round = false:
The source code for the Laplacian of Gaussian kernel is as follow:
inline QVector<qreal> edgeKernel(int radius,
qreal sigma,
qreal scaleW,
qreal offsetW,
bool round=false,
int *kl=NULL)
{
int kw = 2 * radius + 1;
QVector<qreal> kernel(kw * kw);
for (int j = 0; j < kw; j++)
for (int i = 0; i < kw; i++) {
qreal x = i - radius;
qreal y = j - radius;
qreal k = scaleW * (x * x + y * y - sigma * sigma) / (pow(sigma, 5) * sqrt(2 * M_PI));
qreal sigma2 = -2 * sigma * sigma;
qreal weight = k * exp((x * x + y * y) / sigma2) + offsetW;
kernel[i + j * kw] = round? qRound(weight): weight;
}
*kl = kw;
return kernel;
}
Some LoG kernel examples:
| Kernel | Parameters |
|---|---|
| $$\begin{vmatrix}0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0\end{vmatrix}$$ | radius = 1 sigma = 0.46 scaleW = 1 offsetXY = 0 round = true |
| $$\begin{vmatrix}1 & 1 & 1 \\ 1 & -8 & 1 \\ 1 & 1 & 1\end{vmatrix}$$ | radius = 1 sigma = 0.21 scaleW = 0.2 offsetW = 1 round = true |
and this is the result for the second kernel:
Visit my github to get the full source code.
No comments:
Post a Comment