ExampleSVD3x3.cpp

/* This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

/** @example ExampleSVD3x3.cpp
    @relates SVD3x3, SVD
    @ingroup g_examples    
    @brief   Example for singular value decomposition of a 3x3 matrix
    @author  woelk 07/2004 
*/

#include <MathAlgo/SVD.hh>
#include <MathAlgo/SVD3x3.hh>
#include <Base/Math/Random.hh>
#include <Base/Math/Matrix3x3.hh>
#include <Base/Math/Matrix.hh>

using namespace BIAS;
using namespace std;

#define F_SVD_THRESH 0.1

int main()
{
  // Define 3x3 matrix (fundamental matrix)
  Matrix<double> F(3, 3, MatrixZero);
  F[0][0] =0.00075284066373225704;    
  F[0][1] =-0.099906426438442303;      
  F[0][2] =0.16267766990357244;
  F[1][0] =0.10333163045064053;
  F[1][1] =0.0014727467665194025;
  F[1][2] =0.6803859609102173;
  F[2][0] =-0.16234020366691171;      
  F[2][1] =-0.68087420277635458;     
  F[2][2] =0.0021896762792443826;
  F.PrintPretty(cout, "F");

  // Compute with SVD3x3
  SVD3x3 svdF3x3(F, F_SVD_THRESH);
  Matrix<double> S(3, 3, MatrixZero);
  for (int i = 0; i < 3; i++)
    S[i][i] = svdF3x3.GetS()[i];

  cout << endl << "-- SVD3x3 --" << endl;
  svdF3x3.GetU().PrintPretty(cout, "U");
  S.PrintPretty(cout, "S");
  svdF3x3.GetVT().PrintPretty(cout, "VT");

  Matrix<double> USVT = svdF3x3.GetU() * S * svdF3x3.GetVT();
  USVT.PrintPretty(cout, "U*S*VT");

  cout << endl << "|F - U*S*VT| = " << Matrix<double>(F-USVT).NormL2() << endl
       << endl << "NullspaceDim() = " << svdF3x3.NullspaceDim() << endl;
  
  // Compare results with SVD
  SVD svdF(F, F_SVD_THRESH);
  for (int i = 0; i < 3; i++)
    S[i][i] = svdF.GetS()[i];

  cout << endl << "-- SVD --" << endl;
  svdF.GetU().PrintPretty(cout, "U");
  S.PrintPretty(cout, "S");
  svdF.GetVT().PrintPretty(cout, "VT");

  USVT =  svdF.GetU() * S * svdF.GetVT();
  USVT.PrintPretty(cout, "U*S*VT");

  cout << endl << "|F - U*S*VT| = " << Matrix<double>(F-USVT).NormL2() << endl
       << endl << "NullspaceDim() = " << svdF3x3.NullspaceDim() << endl
       << endl;
  
  return 0;
}