Skip to content
GitLab
Commit bd4d773f authored by nanored's avatar nanored
Browse files

location error

parent 82c43f1e
Hide whitespace changes
Inline Side-by-side
TPs/harris.cpp
View file @ bd4d773f
......@@ -9,19 +9,31 @@
using namespace std;
using namespace Imagine;
const static double focal = 2630;
/************************/
/****** PARAMETERS ******/
/************************/
const static double focal = 2600;
// Standart deviation of the gaussian used for the gradient
const static double sigmaD = 1.0;
const static int kD = 2;
const static int kD = 2; // (2*kD+1) is the size of the convolution filter
// Standart deviation of the gaussian used for the correlation
const static double sigmaI = 4.0;
const static int kI = 5;
const static int kI = 5;// (2*kI+1) is the size of the convolution filter
// Constant in the Harris response
const static double kappa = 0.045;
// Threshold before applying NMS and ANMS
const static double threshold = 0.005;
// Radius for NMS
const static double rNMS = 5;
// The constant c in NMS and ANMS
const static double c = 0.9;
// Number of points to plot
const static int N = 200;
/************************/
// Return a 1D gaussian kernel of standard deviation sigma of size 2k+1
vector<double> gaussianKernel(double sigma, int k) {
......@@ -84,6 +96,10 @@ Image<double> applyFilterY(const Image<double> &I, const vector<double> &K) {
return IK;
}
// Compute the Harris responses. R is an image of size (I.w - 2*offset, I.h - 2*offset) storing the Harris
// responses for each pixel. A pixel (i, j) in R corresponds to (i+offset, j+offset) in I.
// Returns a sorted list of pair (minus the response, pixel index). Then the first element corresponds to
// the highest response. The pixel index for the pixel (i, j) is given by i + j*w.
vector<pair<double, int>> Harris(const Image<double> &I, double threshold, Image<double> &R, int &offset) {
int w = I.width(), h = I.height();
......@@ -139,6 +155,7 @@ vector<pair<double, int>> Harris(const Image<double> &I, double threshold, Image
return responses;
}
// Apply the homography H to the image I and return the result. (x0, y0) are the true coordinates of the pixel (0,0)
Image<double> applyHomography(const Image<double> &I, const FMatrix<double,3,3> &H, double &x0, double &y0) {
double w = I.width(), h = I.height();
FMatrix<double,3,3> invH = inverse(H);
......@@ -175,6 +192,9 @@ Image<double> applyHomography(const Image<double> &I, const FMatrix<double,3,3>
return I2;
}
// Remove the points near of the border in the Harris responses on the image transform by the homography H.
// Important because outside the image there are black pixel, then this add a gradient that doesn't exist.
// w and h are the size of the orignal image and w2 is the width of the transformed image.
void cleanResponses(vector<pair<double, int>> &points, int w2, double x0, double y0,
int w, int h, const FMatrix<double,3,3> &H)
{
......@@ -202,7 +222,8 @@ void cleanResponses(vector<pair<double, int>> &points, int w2, double x0, double
points.resize(j);
}
FMatrix<double,3,3> getRotattion(double theta) {
// Return the homography corresponding to a rotation of an angle theta of the image
FMatrix<double,3,3> getRotation(double theta) {
FMatrix<double,3,3> H;
theta *= -1;
H(0, 0) = cos(theta), H(0, 1) = -sin(theta), H(0, 2) = 0;
......@@ -211,17 +232,19 @@ FMatrix<double,3,3> getRotattion(double theta) {
return H;
}
// Return the homography corresponding to a scale of the image
FMatrix<double,3,3> getScale(double scale) {
FVector<double,3> d;
d[0] = scale, d[1] =scale, d[2] = 1;
return Diagonal(d);
}
// Return a random homography defined by a rotation of the camera
FMatrix<double,3,3> getRandomHomography(double w, double h) {
FMatrix<double,3,3> Rx, Ry, Rz, K, H;
Rx = getRotattion(1.0 * (doubleRandom() - 0.5));
Ry = getRotattion(1.0 * (doubleRandom() - 0.5));
Rz = getRotattion(0.3 * (doubleRandom() - 0.5));
Rx = getRotation(1.0 * (doubleRandom() - 0.5));
Ry = getRotation(1.0 * (doubleRandom() - 0.5));
Rz = getRotation(0.3 * (doubleRandom() - 0.5));
swap(Ry(0, 1), Ry(0, 2));
swap(Ry(1, 0), Ry(2, 0));
swap(Ry(1, 1), Ry(2, 2));
......@@ -235,6 +258,7 @@ FMatrix<double,3,3> getRandomHomography(double w, double h) {
return H;
}
// Apply gaussian nois of standard deviation sigma on I
Image<double> applyNoise(const Image<double> &I, double sigma) {
int w = I.width(), h = I.height();
Image<double> I2(w, h);
......@@ -245,6 +269,7 @@ Image<double> applyNoise(const Image<double> &I, double sigma) {
return I2;
}
// Convert a double image in byte image
Image<byte> convertIm(const Image<double> &I) {
int w = I.width(), h = I.height();
Image<byte> Ib(w, h);
......@@ -254,10 +279,13 @@ Image<byte> convertIm(const Image<double> &I) {
return Ib;
}
// Save a double Image
bool save(const Image<double> &I, string name) {
return save(convertIm(I), name);
}
// Apply threshold after computing responses. Usefull because we set the threshold to 0 when computing
// Harris response to plot the number of points in function of the threshold
void applyThreshold(vector<pair<double, int>> &responses) {
double t = responses[0].first * threshold;
int s = 1;
......@@ -318,6 +346,7 @@ vector<pair<double, int>> ANMS(const vector<pair<double, int>> &points, int w) {
return res;
}
// Plot the N best points
void plotPoints(const vector<pair<double, int>> &points, int N, int w, double fact=1.0, int ox=0, int oy=0) {
for(int i = 0; i < N; i++) {
int p = points[i].second;
......@@ -325,10 +354,11 @@ void plotPoints(const vector<pair<double, int>> &points, int N, int w, double fa
}
}
bool save_responses(const vector<pair<double, int>> &res, string title, string name) {
// Save the responses in res in a file 'name' to plot the number of detection in funcion of the threshold
bool save_responses(const vector<pair<double, int>> &res, string title) {
int n = res.size();
ofstream f;
f.open(name);
f.open(title + ".txt");
if(!f.is_open()) return false;
f << title << "\n" << n << "\n";
for(int i = n-1; i >= 0; i--)
......@@ -337,26 +367,101 @@ bool save_responses(const vector<pair<double, int>> &res, string title, string n
return true;
}
void compute_NMS_and_display(Window W, int i, Image<byte> Ib, vector<pair<double, int>> &responses,
int w, int h, const Image<double> &R, int offset,
int N, double fact=1.0, int ox=0, int oy=0)
bool location_error(const vector<pair<double, int>> &res, vector<pair<double, int>> &res2,
double x0, double y0, FMatrix<double,3,3> H, int w, int w2, string title) {
// We sort points in function of the abcisse
vector<pair<double, double>> points;
for(auto &a : res) points.push_back({a.second % w, a.second / w});
sort(points.begin(), points.end());
int n = points.size();
FVector<double,3> v;
v[2] = 1;
vector<double> ts;
FMatrix<double,3,3> invH = inverse(H);
for(auto p : res2) {
v[0] = p.second % w2 + x0, v[1] = p.second / w2 + y0;
v = invH * v;
v /= v[2];
double x = v[0], y = v[1];
double dist = 1e8;
// dichotomie to find the nearest point in abscisse
int a = 0, b = n-1;
while(a < b) {
int c = (a + b) / 2;
if(x < points[c].first) b = c-1;
else a = c+1;
}
// Explore both sides of c
bool stop;
do {
stop = true;
if(a < n) {
double dx = points[a].first - x;
dx *= dx;
if(dx < dist) {
stop = false;
double dy = points[a].second - y;
dist = min(dist, dx + dy*dy);
a ++;
}
}
if(b >= 0) {
double dx = points[b].first - x;
dx *= dx;
if(dx < dist) {
stop = false;
double dy = points[b].second - y;
dist = min(dist, dx + dy*dy);
b --;
}
}
} while(!stop);
ts.push_back(dist);
}
sort(ts.begin(), ts.end());
// Then save
int n2 = res2.size();
ofstream f;
f.open(title + "_loc.txt");
if(!f.is_open()) return false;
f << title << "\n" << n2 << "\n";
for(double t : ts) f << t << "\n";
f.close();
return true;
}
// Function to compute NMS, ANMS and display results
void compute_NMS_and_display(Window W, int i, Image<byte> Ib, vector<pair<double, int>> &responses, string name,
vector<pair<double, int>> &res0, vector<pair<double, int>> &nms0, vector<pair<double, int>> &anms0,
int w, int h, const Image<double> &R, int offset, int N,
int w0=0, double x0=0, double y0=0, FMatrix<double,3,3> H=FMatrix<double,3,3>::Identity(),
double fact=1.0, int ox=0, int oy=0)
{
setActiveWindow(W, i);
display(Ib, ox, oy, false, fact);
plotPoints(responses, N, w, fact, ox, oy);
save_responses(responses, name);
applyThreshold(responses);
if(res0.empty()) res0 = responses;
else location_error(res0, responses, x0, y0, H, w0, w, name);
// best NMS responses
setActiveWindow(W, i+5);
display(Ib, ox, oy, false, fact);
vector<pair<double, int>> nms = NMS(responses, w, h, R, offset);
plotPoints(nms, N, w, fact, ox, oy);
if(nms0.empty()) nms0 = nms;
else location_error(nms0, nms, x0, y0, H, w0, w, name+"_nms");
// best ANMS responses
setActiveWindow(W, i+10);
display(Ib, ox, oy, false, fact);
vector<pair<double, int>> anms = ANMS(responses, w);
plotPoints(anms, N, w, fact, ox, oy);
if(anms0.empty()) anms0 = anms;
else location_error(anms0, anms, x0, y0, H, w0, w, name+"_anms");
}
int main() {
......@@ -384,25 +489,25 @@ int main() {
// Some variables
Image<double> R; // Image of responses
int offset; // Offset between original image and R
vector<pair<double, int>> responses; // Sorted responses
FMatrix<double,3,3> H;
Image<double> I2;
Image<byte> I2b;
double x0, y0;
int w2, h2;
double fact;
int ox, oy;
vector<pair<double, int>> responses, responses0, nms0, anms0; // Sorted responses
FMatrix<double,3,3> H; // Homography
Image<double> I2; // Transformed image
Image<byte> I2b; // Byte image to display result
double x0, y0; // true coordinates of pixel (0, 0) in I2
int w2, h2; // Size of I2
double fact; // factor of display for I2
int ox, oy; // display offset for I2
// best Harris responses
cout << "Original" << endl;
responses = Harris(I, 0, R, offset);
save_responses(responses, "Original", "original.txt");
compute_NMS_and_display(W, 0, I0, responses, w, h, R, offset, N);
compute_NMS_and_display(W, 0, I0, responses, "original",
responses0, nms0, anms0, w, h, R, offset, N);
// best Harris responses on rotated image
cout << "Rotation" << endl;
double angle = 0.4;
H = getRotattion(angle);
H = getRotation(angle);
I2 = applyHomography(I, H, x0, y0);
I2b = convertIm(I2);
w2 = I2.width(), h2 = I2.height();
......@@ -410,17 +515,18 @@ int main() {
ox = 0.5 * (w - w2*fact), oy = 0.5 * (h - h2*fact);
responses = Harris(I2, 0, R, offset);
cleanResponses(responses, w2, x0, y0, w, h, H);
save_responses(responses, "Rotation", "rotation.txt");
compute_NMS_and_display(W, 1, I2b, responses, w2, h2, R, offset, N, fact, ox, oy);
compute_NMS_and_display(W, 1, I2b, responses, "rotation", responses0, nms0, anms0,
w2, h2, R, offset, N, w, x0, y0, H, fact, ox, oy);
// best Harris responses on noised image
cout << "Noise" << endl;
H = FMatrix<double,3,3>::Identity();
double sigma = 20;
I2 = applyNoise(I, sigma);
I2b = convertIm(I2);
responses = Harris(I2, 0, R, offset);
save_responses(responses, "Noise", "noise.txt");
compute_NMS_and_display(W, 2, I2b, responses, w, h, R, offset, N);
compute_NMS_and_display(W, 2, I2b, responses, "noise", responses0, nms0, anms0,
w, h, R, offset, N, w);
// best Harris responses on scaled image
cout << "Scale" << endl;
......@@ -432,8 +538,8 @@ int main() {
fact = min(1.0, min(double(w) / double(w2), double(h) / double(h2)));
ox = 0.5 * (w - w2*fact), oy = 0.5 * (h - h2*fact);
responses = Harris(I2, 0, R, offset);
save_responses(responses, "Scale", "scale.txt");
compute_NMS_and_display(W, 3, I2b, responses, w2, h2, R, offset, N, fact, ox, oy);
compute_NMS_and_display(W, 3, I2b, responses, "scale", responses0, nms0, anms0,
w2, h2, R, offset, N, w, x0, y0, H, fact, ox, oy);
// best Harris responses on viewpoint changed image
cout << "Viewpoint" << endl;
......@@ -445,8 +551,8 @@ int main() {
ox = 0.5 * (w - w2*fact), oy = 0.5 * (h - h2*fact);
responses = Harris(I2, 0, R, offset);
cleanResponses(responses, w2, x0, y0, w, h, H);
save_responses(responses, "Viewpoint change", "viewpoint.txt");
compute_NMS_and_display(W, 4, I2b, responses, w2, h2, R, offset, N, fact, ox, oy);
compute_NMS_and_display(W, 4, I2b, responses, "viewpoint", responses0, nms0, anms0,
w2, h2, R, offset, N, w, x0, y0, H, fact, ox, oy);
// Stop the programm when clicked
endGraphics();
......
Supports Markdown
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