ExampleRegressionPlane.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 ExampleRegressionPlane.cpp
    @relates RANSAC, SVD
    @ingroup g_examples    
    @brief  this is a class using RANSAC to compute a regression plane with SVD
    @author MIP
*/



#include <iostream>
#include <iomanip>
 
// use BIAS Base Math library and Algos:
#include <Base/Math/Math.hh>

#include <MathAlgo/SVD.hh>


using namespace std;
using namespace BIAS;



void set_example_designmatrix( BIAS::Matrix<double> & A) {
  // fill the 3d  points rowwise into A using the scheme
  const int nr_of_points = 7;
  
  A.newsize( nr_of_points, 4 );
  
  //        -1            X_i             Y_i             Z_i 
  A[0][0] = -1; A[0][1] = 0;    A[0][2] = 0;    A[0][3] = 1; // 1. example points
  A[1][0] = -1; A[1][1] = 1;    A[1][2] = 0;    A[1][3] = 1; // 2.
  A[2][0] = -1; A[2][1] = 0;    A[2][2] = 1;    A[2][3] = 1; // 3.
  A[3][0] = -1; A[3][1] = 1;    A[3][2] = 1;    A[3][3] = 1; // 4.
  A[4][0] = -1; A[4][1] = 0.5;  A[4][2] = 0.5;  A[4][3] = 1; // 5.
  
  A[5][0] = -1; A[5][1] = 0.51;  A[5][2] = 0.59;  A[5][3] = 0.9; // 6.
  A[6][0] = -1; A[6][1] = 3.49;  A[6][2] = 0.55;  A[6][3] = 1.1; // 7.
}


void calc_Regressionsebene( const BIAS::Matrix<double> & A ) {
  cout << "calc_Regressionsebene " << endl;
  cout << "Designmatrix (Matrix<double>) A: " << A << endl;
  
  // get a solution (linear fit) of this equation system by its SVD
  SVD svd( A );
  
  // cout<<"U = " << svd.GetU() <<endl;
  // cout<<"S = " << svd.GetS() <<endl;
  // cout<<"V = " << svd.GetV() <<endl;
  // cout<<"Nullspace dim = " << svd.NullspaceDim() <<endl;;
  
  BIAS::Vector<double> solution(4);
  solution = svd.GetNullvector(); // get the last Nullvector which is the solution 
  cout <<"solution: " << solution << endl;
  
  Vector3<double> n_vec;
  double d;
  // decompose the solution into HNF of a plane
  d = solution[0];
  n_vec[0] = solution[1];
  n_vec[1] = solution[2];
  n_vec[2] = solution[3]; 
  // now (d, n_vec ) is a plane in normalform
  
  // convert to HNF representation: 
  d = d / n_vec.NormL2(); // divide by actual length for HNF
  n_vec.Normalize();// normalize n_vec to length 1 (for HNF)
  
  cout << " plane in HNF:  d = "<<d<<endl;
  cout << " n_vce'0 : " << n_vec << endl;  
}  // ------------------------------------------------------------



// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
int main() { 
    BIAS::Matrix<double> M;
    set_example_designmatrix( M ); // set example points into this nPoints x 4 designmatrix
    calc_Regressionsebene( M ); // calculate the best fit for the regression plane
    return 0;
}