Commit 750c887c authored by nanored's avatar nanored
Browse files

TP 3

parent a0d67b86
......@@ -11,7 +11,8 @@
using namespace Imagine;
using namespace std;
static const double BETA = 0.01f; // Probability of failure
static const bool REFINIG = false; // If true the refining is used
static const double BETA = 0.001f; // Probability of failure
struct Match {
double x1, y1, x2, y2;
......@@ -50,7 +51,7 @@ void algoSIFT(Image<Color,2> I1, Image<Color,2> I2,
// Return a good normalization factor for matches
double getNormalizationFactor(vector<Match>& matches) {
// the factor of normalization minimizes
// SUM[on m in matches] log(fact * m.x1)^2 * ... * log(fact * m.y2)^2
// SUM[on m in matches] log(fact * m.x1)^2 + ... + log(fact * m.y2)^2
double sum = 0;
for(const Match &m : matches) {
sum += log(m.x1 + 1);
......@@ -79,10 +80,53 @@ vector<int> getInliers(const vector<Match>& matches, const FMatrix<double,3,3> F
return inliers;
}
// Compute the fundamental matrix using the n first matches and using the normalization matrix N
// (Without RANSAC)
FMatrix<double,3,3> tempF(const vector<Match>& matches, int n, const FMatrix<double,3,3>& N) {
Matrix<double> A(max(n,9), 9); // The matrix A of the linear problem
Matrix<double> U, Vt; // Matrixes used in SVD decomposition of A
Vector<double> S; // Eigenvalues in the SVD decomposition of A
// We fill the last row of A with zeros if n=8
if(n == 8)
for(int i = 0; i < 9; i++) A(8, i) = 0;
FMatrix<double,3,3> F; // Fundamental matrix computed in the iteration
FMatrix<double,3,3> U2, Vt2; // Matrixes used in SVD decomposition of F
FVector<double,3> S2; // Eigenvalues in the SVD decomposition of F
// We compute A
for(int i = 0; i < n; i++) {
// Don't forget the normalization factor
double x1 = matches[i].x1 * N(0, 0);
double y1 = matches[i].y1 * N(1, 1);
double x2 = matches[i].x2 * N(0, 0);
double y2 = matches[i].y2 * N(1, 1);
A(i, 0) = x1*x2; A(i, 1) = x1*y2; A(i, 2) = x1;
A(i, 3) = y1*x2; A(i, 4) = y1*y2; A(i, 5) = y1;
A(i, 6) = x2; A(i, 7) = y2; A(i, 8) = 1;
}
// f is the last row of Vt in the SVD of A
svd(A, U, S, Vt);
Vector<double> V9 = Vt.getRow(8);
for(int i = 0; i < 3; i++)
for(int j = 0; j < 3; j++)
F(i, j) = V9[3*i+j];
// Add the constraint det F = 0
svd(F, U2, S2, Vt2);
S2[2] = 0;
F = U2 * Diagonal(S2) * Vt2;
// Again, don't forget normalization
F = N * F * N;
return F;
}
// RANSAC algorithm to compute F from point matches (8-point algorithm)
// Parameter matches is filtered to keep only inliers as output.
FMatrix<double,3,3> computeF(vector<Match>& matches) {
const double distMax = 1.2f; // Pixel error for inlier/outlier discrimination
const double distMax = 1.15f; // Pixel error for inlier/outlier discrimination
int Niter=100000; // Adjusted dynamically
FMatrix<double,3,3> bestF;
vector<Match> all=matches;
......@@ -96,16 +140,6 @@ FMatrix<double,3,3> computeF(vector<Match>& matches) {
return bestF;
}
FMatrix<double,9,9> A; // The matrix A of the linear problem
FMatrix<double,9,9> U, Vt; // Matrixes used in SVD decomposition of A
FVector<double,9> S; // Eigenvalues in the SVD decomposition of A
// We fill the row line of A with zeros
for(int i = 0; i < 9; i++) A(8, i) = 0;
FMatrix<double,3,3> F; // Fundamental matrix computed in the iteration
FMatrix<double,3,3> U2, Vt2; // Matrixes used in SVD decomposition of F
FVector<double,3> S2; // Eigenvalues in the SVD decomposition of F
// Factor of normalization
double norm_fact = getNormalizationFactor(matches);
FVector<double,3> diagN;
......@@ -120,33 +154,12 @@ FMatrix<double,3,3> computeF(vector<Match>& matches) {
int niter = 0; // Number of iteration realized
while(niter < Niter) {
// We compute A
for(int i = 0; i < 8; i++) {
// Draw a match
// Draw 8 matches
for(int i = 0; i < 8; i++)
swap(matches[i], matches[unifs[i](dre)]);
// Don't forget the normalization factor
double x1 = matches[i].x1 * norm_fact;
double y1 = matches[i].y1 * norm_fact;
double x2 = matches[i].x2 * norm_fact;
double y2 = matches[i].y2 * norm_fact;
A(i, 0) = x1*x2; A(i, 1) = x1*y2; A(i, 2) = x1;
A(i, 3) = y1*x2; A(i, 4) = y1*y2; A(i, 5) = y1;
A(i, 6) = x2; A(i, 7) = y2; A(i, 8) = 1;
}
// f is the last row of Vt in the SVD of A
svd(A, U, S, Vt);
FVector<double, 9> V9 = Vt.getRow(8);
for(int i = 0; i < 3; i++)
for(int j = 0; j < 3; j++)
F(i, j) = V9[3*i+j];
// Add the constraint det F = 0
svd(F, U2, S2, Vt2);
S2[2] = 0;
F = U2 * Diagonal(S2) * Vt2;
// Again, don't forget normalization
F = N * F * N;
// Compute F for these matches
FMatrix<double,3,3> F = tempF(matches, 8, N);
// Compute inliers
vector<int> inliers = getInliers(all, F, distMax);
......@@ -163,12 +176,18 @@ FMatrix<double,3,3> computeF(vector<Match>& matches) {
niter ++;
}
/************ END OF TODO ************/
// Updating matches with inliers only
matches.clear();
for(size_t i=0; i<bestInliers.size(); i++)
matches.push_back(all[bestInliers[i]]);
// Recompute F if the REFINING variable is set to true
if(REFINIG)
bestF = tempF(matches, matches.size(), N);
/************ END OF TODO ************/
return bestF;
}
......@@ -184,18 +203,19 @@ void displayEpipolar(Image<Color> I1, Image<Color> I2,
// --------------- TODO ------------
int w = I1.width(), w2 = I2.width();
FVector<double, 3> point, line;
Color c(rand()%256,rand()%256,rand()%256);
if(x < w) {
point[0] = x; point[1] = y; point[2] = 1;
line = transpose(F) * point;
line /= -line[1];
drawLine(w, line[2], w+w2, line[0]*w2+line[2], RED, 1);
drawLine(w, line[2], w+w2, line[0]*w2+line[2], c, 1);
} else {
point[0] = x-w; point[1] = y; point[2] = 1;
line = F * point;
line /= -line[1];
drawLine(0, line[2], w, line[0]*w+line[2], RED, 1);
drawLine(0, line[2], w, line[0]*w+line[2], c, 1);
}
drawCircle(x, y, 4, BLUE, 3);
drawCircle(x, y, 4, c, 3);
}
}
......
// Imagine++ project
// Project: Seeds
// Author: Pascal Monasse
#include <Imagine/Graphics.h>
#include <Imagine/Images.h>
#include <queue>
#include <iostream>
using namespace Imagine;
using namespace std;
/// Min and max disparities
static const int dMin=-30, dMax=-7;
/// Min NCC for a seed
static const float nccSeed=0.95;
/// Radius of patch for correlation
static const int win=(9-1)/2;
/// To avoid division by 0 for constant patch
static const float EPS=0.01f;
/// A seed
struct Seed {
Seed(int x0, int y0, int d0, float ncc0)
: x(x0), y(y0), d(d0), ncc(ncc0) {}
int x,y, d;
float ncc;
};
/// Order by NCC
bool operator<(const Seed& s1, const Seed& s2) {
return (s1.ncc<s2.ncc);
}
/// 4-neighbors
static const int dx[]={+1, 0, -1, 0};
static const int dy[]={ 0, -1, 0, +1};
/// Display disparity map
static void displayDisp(const Image<int> disp, Window W, int subW) {
Image<Color> im(disp.width(), disp.height());
for(int j=0; j<disp.height(); j++)
for(int i=0; i<disp.width(); i++) {
if(disp(i,j)<dMin || disp(i,j)>dMax)
im(i,j)=Color(0,255,255);
else {
int g = 255*(disp(i,j)-dMin)/(dMax-dMin);
im(i,j)= Color(g,g,g);
}
}
setActiveWindow(W,subW);
display(im);
save(im, "res"+to_string(subW)+".jpg");
showWindow(W,subW);
}
/// Show 3D window
static void show3D(const Image<Color> im, const Image<int> disp) {
#ifdef IMAGINE_OPENGL // Imagine++ must have been built with OpenGL support...
// Intrinsic parameters given by Middlebury website
const float f=3740;
const float d0=-200; // Doll images cropped by this amount
const float zoom=2; // Half-size images, should double measured disparity
const float B=0.160; // Baseline in m
FMatrix<float,3,3> K(0.0f);
K(0,0)=-f/zoom; K(0,2)=disp.width()/2;
K(1,1)= f/zoom; K(1,2)=disp.height()/2;
K(2,2)=1.0f;
K = inverse(K);
K /= K(2,2);
std::vector<FloatPoint3> pts;
std::vector<Color> col;
for(int j=0; j<disp.height(); j++)
for(int i=0; i<disp.width(); i++)
if(dMin<=disp(i,j) && disp(i,j)<=dMax) {
float z = B*f/(zoom*disp(i,j)+d0);
FloatPoint3 pt((float)i,(float)j,1.0f);
pts.push_back(K*pt*z);
col.push_back(im(i,j));
}
Mesh mesh(&pts[0], pts.size(), 0,0,0,0,VERTEX_COLOR);
mesh.setColors(VERTEX, &col[0]);
Window W = openWindow3D(512,512,"3D");
setActiveWindow(W);
showMesh(mesh);
#else
std::cout << "No 3D: Imagine++ not built with OpenGL support" << std::endl;
#endif
}
/// Correlation between patches centered on (i1,j1) and (i2,j2). The values
/// m1 or m2 are subtracted from each pixel value.
static float correl(const Image<byte>& im1, int i1,int j1,float m1,
const Image<byte>& im2, int i2,int j2,float m2) {
float dist=0.0f;
// ------------- TODO -------------
float norm1=EPS, norm2=EPS;
for(int di = -win; di <= win; di ++) {
for(int dj = -win; dj <= win; dj++) {
float dim1 = im1(i1+di, j1+dj) - m1;
float dim2 = im2(i2+di, j2+dj) - m2;
dist += dim1 * dim2;
norm1 += dim1 * dim1;
norm2 += dim2 * dim2;
}
}
dist /= (sqrt(norm1) * sqrt(norm2));
// ---------- END OF TODO ----------
return dist;
}
/// Sum of pixel values in patch centered on (i,j).
static float sum(const Image<byte>& im, int i, int j) {
float s=0.0f;
// ------------- TODO -------------
for(int x = i-win; x <= i+win; x ++)
for(int y = j-win; y <= j+win; y++)
s += im(x, y);
// ----------- END OF TODO --------
return s;
}
/// Centered correlation of patches of size 2*win+1.
static float ccorrel(const Image<byte>& im1,int i1,int j1,
const Image<byte>& im2,int i2,int j2) {
float m1 = sum(im1,i1,j1);
float m2 = sum(im2,i2,j2);
int w = 2*win+1;
return correl(im1,i1,j1,m1/(w*w), im2,i2,j2,m2/(w*w));
}
/// Compute disparity map from im1 to im2, but only at points where NCC is
/// above nccSeed. Set to true the seeds and put them in Q.
static void find_seeds(Image<byte> im1, Image<byte> im2,
float nccSeed,
Image<int>& disp, Image<bool>& seeds,
std::priority_queue<Seed>& Q) {
disp.fill(dMin-1);
seeds.fill(false);
while(! Q.empty())
Q.pop();
const int maxy = std::min(im1.height(),im2.height());
const int refreshStep = (maxy-2*win)*5/100;
for(int y=win; y+win<maxy; y++) {
if((y-win-1)/refreshStep != (y-win)/refreshStep)
std::cout << "Seeds: " << 5*(y-win)/refreshStep <<"%\r"<<std::flush;
for(int x=win; x+win<im1.width(); x++) {
// ------------- TODO -------------
float ncc = -1; // Best correlation found
int dmax = 0; // displacement of the best correlation
for(int d = dMin; d <= dMax; d++) {
int x2 = x+d;
if(x2 < win || x2+win >= im2.width()) continue;
float cor = ccorrel(im1, x, y, im2, x2, y);
if(cor > ncc) {
ncc = cor;
dmax = d;
}
}
if(ncc >= nccSeed) {
Q.emplace(x, y, dmax, ncc);
disp(x, y) = dmax;
seeds(x, y) = true;
}
// ----------- END OF TODO --------
}
}
std::cout << std::endl;
}
/// Propagate seeds
static void propagate(Image<byte> im1, Image<byte> im2,
Image<int>& disp, Image<bool>& seeds,
std::priority_queue<Seed>& Q) {
while(! Q.empty()) {
Seed s=Q.top();
Q.pop();
for(int i=0; i<4; i++) {
float x=s.x+dx[i], y=s.y+dy[i];
if(0 <= x-win && 0 <= y-win &&
x+win < im2.width() && y+win < im2.height() && !seeds(x,y)) {
// ------------- TODO -------------
float ncc = -1; // Best correlation found
int dmax = s.d; // displacement of the best correlation
for(int d = max(dMin, s.d-1); d <= min(dMax, s.d+1); d++) {
int x2 = x+d;
if(x2 < win || x2+win >= im2.width()) continue;
float cor = ccorrel(im1, x, y, im2, x2, y);
if(cor > ncc) {
ncc = cor;
dmax = d;
}
}
Q.emplace(x, y, dmax, ncc);
disp(x, y) = dmax;
seeds(x, y) = true;
// ----------- END OF TODO --------
}
}
}
}
int main()
{
// Load and display images
Image<Color> I1, I2;
if( ! load(I1, srcPath("im1.jpg")) ||
! load(I2, srcPath("im2.jpg")) ) {
cerr<< "Unable to load images" << endl;
return 1;
}
std::string names[5]={"image 1","image 2","dense","seeds","propagation"};
Window W = openComplexWindow(I1.width(), I1.height(), "Seeds propagation",
5, names);
setActiveWindow(W,0);
display(I1,0,0);
setActiveWindow(W,1);
display(I2,0,0);
Image<int> disp(I1.width(), I1.height());
Image<bool> seeds(I1.width(), I1.height());
std::priority_queue<Seed> Q;
// Dense disparity
find_seeds(I1, I2, -1.0f, disp, seeds, Q);
displayDisp(disp,W,2);
// Only seeds
find_seeds(I1, I2, nccSeed, disp, seeds, Q);
displayDisp(disp,W,3);
// Propagation of seeds
propagate(I1, I2, disp, seeds, Q);
displayDisp(disp,W,4);
// Show 3D (use shift click to animate)
show3D(I1,disp);
endGraphics();
return 0;
}
......@@ -140,4 +140,11 @@ Dans notre cas la précision est donné par les distances $d(\bx_i', F^\trans \b
$$ d(\bx_i', F^\trans \bx_i) = \dfrac{\left| \left( F^\trans \bx_i \right)_1 x_i' + \left( F^\trans \bx_i \right)_2 y_i' + \left( F^\trans \bx_i \right)_3 \right|}{\sqrt{\left( F^\trans \bx_i \right)_1^2 + \left( F^\trans \bx_i \right)_2^2}} $$
Supposons que nous avons $m$ paires de points corrects. Alors la probabilité de prendre $k$ paires de points corrects est $(m / n)^k$. On veut que la probabilité de $N_{iter}$ itérations de RANSAC donnent que des sous-ensembles contenant une mauvaise paire de points soit inférieur à $\beta = 1\%$ :
$$ \left( 1 - (m / n)^k \right)^{N_{iter}} \leqslant \beta \qquad \Leftrightarrow \qquad N_{iter} \geqslant \dfrac{\log \beta}{\log \left( 1 - (m / n)^k \right)} $$
$m$ est inconnu mais on connaît une borne inférieur qui est le plus grand nombre de paires de points expliqué par un modèle obtenu jusque là. A chaque fois qu'on trouve un meilleur modèle on peut alors recalculer $m$.
\ No newline at end of file
$m$ est inconnu mais on connaît une borne inférieur qui est le plus grand nombre de paires de points expliqué par un modèle obtenu jusque là. A chaque fois qu'on trouve un meilleur modèle on peut alors recalculer $m$.
\sect{Rectification épipôlaire}
\paragraph{}
Il est commode d'avoir les lignes épipôlaires parallèles et à la même ordonné dans les deux images. La conséquence est que les épipôles sont les points à l'infini horizontalement.
$$ e = e' = \pmat{1 \\ 0 \\ 0} $$
On peut toujours se retrouver dans cette situation en appliquant une rotation virtuelle à la caméra qui résulte en une homographie.
\ No newline at end of file
\chapter{Carte de disparités}
\myminitoc
\sect{Triangulation}
\paragraph{}
Commençons par reprendre les formules binoculaire pour observer un point de l'espace $X$ avec des caméras de matrices de projections $P$ et $P'$. Les coordonnées homogènes de l'images de $X$ par les deux caméras sont notées $x$ et $x'$ :
$$ x = PX \Leftrightarrow [x]_\times PX = 0 \qquad x' = P'X \Leftrightarrow [x']_\times P'X = 0 $$
L'objectif de la triangulation est à partir des deux images $x$ et $x'$, de retrouver le point dans l'espace $X$. Cela est faisable en calculant une SVD pour trouver le noyau d'une matrice :
$$ X \in \ker \pmat{[x]_\times P \\ [x']_\times P'} $$
\paragraph{}
En se plaçant dans le repère de la seconde caméra et en notant $R$ la matrice de rotation du repère de la seconde caméra au repère de la première caméra et $T$ la translation entre les deux caméras. On peut réécrire les équations précédentes avec :
$$ \lambda x = K(RX + T) \qquad \lambda' x' = K' X $$
On pose alors le vecteur $Y = \pmat{X^\trans & 1 & \lambda & \lambda'}^\trans$ qui vérifie le système suivant :
$$ \pmat{KR & KT & -x & 0 \\ K' & 0 & 0 & -x'} Y = 0 $$
On retrouve alors $X$ à l'aide de la SVD de la matrice de droite qui nous donne le noyau.
\paragraph{Récupération de $R$ et $T$}
Supposons que nous connaissons les matrices de calibration $K$ et $K'$ en plus de la matrice essentielle $E$ (obtenable à partir de la matrice fondamentale $F$ grâce à $K$ et $K'$). On essaye alors de retrouver $R$ et $T$. On sait que $E = \left[ T \right]_\times R$.
$$ E^\trans E y = R^\trans [T]_\times^\trans [T]_\times R y = - R^\trans [T]_\times \left([T]_\times (R y) \right) = - R^\trans \left( (T^\trans (R y)) T - (T^\trans T) (Ry) \right) $$
Ce qui donne :
$$ E^\trans E = R^\trans \left( T T^\trans - \| T \|_2^2 I_3 \right) R = - x x^\trans + \| T \|_2^2 I_3 \qquad \text{Avec } x = R^\trans T $$
On observe alors que $E^\trans E x = 0$ et que pour $y$ orthogonal à $x$, $E^\trans E y = \| T \|_2^2 y$. Ainsi $E^\trans E$ est une matrice ortho-diagonalisable de valeur propres $\| T \|_2^2$ et $0$ et de rang 2. Ainsi les valeurs propres de $E$ sont $\sigma_1 = \sigma_2 = \| T \|_2 = \sigma$ et $\sigma_3 = 0$. La SVD de $E$ donne :
$$ E = U \pmat{\sigma \\ & \sigma \\ & & 0 } V^\trans = \sigma U \pmat{0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0} U^\trans U \pmat{0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 1} V^\trans = \left[ T \right]_\times R $$
On identifie alors :
$$ \left[ T \right]_\times = \pm \sigma U \pmat{0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0} U^\trans \qquad R = \pm U \pmat{0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 1} V^\trans $$
Ce qui donne 4 solutions possibles :
$$ \left\{ \begin{array}{lll}
T & = & \pm \sigma U [e_3]_\times U^\trans \\[1mm]
R & = & \pm U R_z \left( \pm \frac{\pi}{2} \right) V^\trans
\end{array} \right. $$
\ No newline at end of file
......@@ -151,5 +151,6 @@
\input{chap1.tex}
\input{chap2.tex}
\input{chap3.tex}
\end{document}
\ No newline at end of file
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