Commit 521d1dc4 authored by nanored's avatar nanored

harris implemented

parent b7a7af8a
#include <Imagine/Images.h>
#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;
using namespace Imagine;
const static double sigmaD = 1.0;
const static int kD = 2;
const static double sigmaI = 4.0;
const static int kI = 5;
const static double kappa = 0.05;
// Return a 1D gaussian kernel of standard deviation sigma of size 2k+1
vector<double> gaussianKernel(double sigma, int k) {
vector<double> K;
double sig2 = sigma * sigma;
double sum = 0;
for(int x = -k; x <= k; x++) {
double v = exp(-0.5 * x*x / sig2);
K.push_back(v);
sum += v;
}
for(int i = 0; i < 2*k+1; i++)
K[i] /= sum;
return K;
}
// Compute the convolution of two 1D kernels K1 and K2
vector<double> convKernels(const vector<double> &K1, const vector<double> &K2) {
int n1 = K1.size(), n2 = K2.size();
int k1 = (n1-1) / 2, k2 = (n2-1) / 2;
int n = abs(n1 - n2) + 1;
int k = (n-1) / 2;
int jmin = k1+k2-k, jmax = k1+k2+k;
vector<double> K(n, 0);
for(int i = 0; i < n1; i++)
for(int j = max(0, jmin-i); j <= min(n2-1, jmax-i); j++)
K[i+j-jmin] += K1[i] * K2[j];
return K;
}
// Apply a filter of kernel K on the rows of an image I
Image<double> applyFilterX(const Image<double> &I, const vector<double> &K) {
int k = (K.size()-1) / 2;
int w = I.width()-2*k, h = I.height();
Image<double> IK(w, h);
#pragma omp paralle for
for(int x = 0; x < w; x++) {
for(int y = 0; y < h; y++) {
IK(x, y) = 0;
for(int i = 0; i < 2*k+1; i++)
IK(x, y) += I(x+i, y) * K[2*k-i];
}
}
return IK;
}
// Apply a filter of kernel K on the colons of an image I
Image<double> applyFilterY(const Image<double> &I, const vector<double> &K) {
int k = (K.size()-1) / 2;
int w = I.width(), h = I.height()-2*k;
Image<double> IK(w, h);
#pragma omp paralle for
for(int x = 0; x < w; x++) {
for(int y = 0; y < h; y++) {
IK(x, y) = 0;
for(int i = 0; i < 2*k+1; i++)
IK(x, y) += I(x, y+i) * K[2*k-i];
}
}
return IK;
}
int main() {
Image<byte> I0;
if(! load(I0, srcPath("im1.jpg"))) {
cerr << "Unable to load the image !" << endl;
return 1;
}
int w = I0.width(), h = I0.height();
Image<double> I(w, h);
for(int i = 0; i < w; i++)
for(int j = 0; j < h; j++)
I(i, j) = I0(i, j);
// Compute Image derivatives
vector<double> deriv = {-1.0, 0, 1.0};
vector<double> G1 = gaussianKernel(sigmaD, kD+1);
vector<double> GD = convKernels(deriv, G1);
int wD = w-2*kD, hD = h-2*kD;
Image<double> Ix = applyFilterX(I, GD);
Image<double> Iy = applyFilterY(I, GD);
// Compute product images
Image<double> Ixx(wD, hD), Iyy(wD, hD), Ixy(wD, hD);
for(int i = 0; i < wD; i++)
for(int j = 0; j < hD; j++) {
Ixx(i, j) = Ix(i, j+kD) * Ix(i, j+kD);
Iyy(i, j) = Iy(i+kD, j) * Iy(i+kD, j);
Ixy(i, j) = Ix(i, j+kD) * Iy(i+kD, j);
}
vector<double> GI = gaussianKernel(sigmaI, kI);
Ixx = applyFilterX(applyFilterY(Ixx, GI), GI);
Iyy = applyFilterX(applyFilterY(Iyy, GI), GI);
Ixy = applyFilterX(applyFilterY(Ixy, GI), GI);
int w2 = Ixx.width(), h2 = Ixx.height();
int offset = (w - w2) / 2;
// Corner Response
Image<double> R(w2, h2);
vector<pair<double, int>> responses;
for(int i = 0; i < w2; i++)
for(int j = 0; j < h2; j++) {
double det = Ixx(i, j) * Iyy(i, j) - Ixy(i, j) * Ixy(i, j);
double trace = Ixx(i, j) + Iyy(i, j);
R(i, j) = det - kappa * trace;
if(i > 1 && j > 1) {
double r = R(i-1, j-1);
bool max = true;
for(int k = i-2; k <= i; k++)
for(int l = j-2; l <= j; l++)
if(R(k, l) > r) max = false;
if(max) responses.push_back({-r, i-1+offset + (j-1+offset)*w});
}
}
cout << "Harris found " << responses.size() << " corner points" << endl;
sort(responses.begin(), responses.end());
// Display
openWindow(w, h);
display(I0);
for(int i = 0; i < 300; i++) {
int p = responses[i].second;
drawCircle(p % w, p / w, 4, RED, 3);
}
click();
endGraphics();
return 0;
}
\ No newline at end of file
......@@ -84,7 +84,7 @@ $$ M_w = \sum_u \sum_v w(u, v) \pmat{I_x^2 & I_x I_y \\ I_x I_y & I_y^2}(u, v) $
Cette matrice étant symétrique, on peut l'orthodiagonalisé pour obtenir les valeurs propres $\lambda_1$ et $\lambda_2$. Finalement notre nouvelle caractérisation pour que $(0, 0)$ soit un coin est que ces deux valeurs propres soient grandes. Pour quantifier le fait qu'elles soient toutes deux grandes sans que l'une soit vraiment plus grande que l'autre, Harris et Stephens ont introduit la valeur :
$$ R_w = \lambda_1 \lambda_2 - \kappa (\lambda_1 + \lambda_2)^2 = \det(M_w) - \kappa \, \trace(M_w)^2 $$
$\kappa < 1/4$ et généralement on choisi $\kappa \in [0.04, 0.06]$. Ici on a fait le travail pour le point (0, 0) mais de manière générale pour le point $(i, j)$ il suffit de prendre :
$$ M_w(i, j) = \left[ w \otimes \pmat{I_x^2 & I_x I_y \\ I_x I_y & I_y^2} \right] (i, j) $$
$$ M_w(i, j) = \left[ w * \pmat{I_x^2 & I_x I_y \\ I_x I_y & I_y^2} \right] (i, j) $$
On prend souvent une gaussienne pour $w$. Enfin une fois qu'on a le score $R_w(i, j)$ pour tous les pixels $(i, j)$, on ne garde plus que ceux qui ont un score supérieur à un certain seuil et qui sont des maximas locaux dans leur voisinage de 8 pixels.
\subs{Maximalité}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment