ExampleConic.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 ExampleConic.cpp
   @relates Conic2D
   @brief Example for conic usage
   @ingroup g_examples
   @author MIP
*/

#include <Base/Common/W32Compat.hh>
#include <Base/Geometry/HomgPoint2D.hh>
#include <Geometry/Conic2D.hh>
#include <Base/Image/Image.hh>
#include <Base/Image/ImageIO.hh>
#include <Geometry/Quadric3D.hh>
#include <Geometry/PMatrix.hh>

using namespace BIAS;
using namespace std;

#define RADIUS 50.0

int main() 
{

  cout << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ExampleConic $$$$$$$$$$$$$$$$$$$$$$"
       << endl;
  PMatrix P_current;
  P_current.SetIdentity();

  Conic2D Conic, Circle;
  HomgPoint2D C(320,240), C2(550,240);
  Conic.SetEllipse(C, 35.0 *M_PI/180.0, 180.0, 70);
  
  // ****** intersection of conic and circle *******
  for (unsigned int i=0; i<15; i++) {
    C2[0] -= 25.0;
    Circle.SetEllipse(C2, 0.0, RADIUS, RADIUS);
    cout << "Circle center is now at "<<C2<<endl;
    if (Conic.IntersectsCircle(C2, RADIUS)) cout << "Intersection :"<<i<<endl;
    else cout << "No Intersection !"<<endl;
    Image<unsigned char> TestImage(640,480);
    TestImage.Clear(200);
    Circle.Draw(TestImage);
    Conic.Draw(TestImage);
    char filename[1000];
    sprintf(filename,"testimage%03d",i);
    ImageIO::Save(filename,TestImage,ImageIO::FF_pgm);
  }
  
  
  // ********** conic from quadric ****************
  HomgPoint3D QC(0,0,20); // quadric center

  Vector3<double> S, G_0(1,0,0), H_0(0,1,0), I_0(0,0,1);
  BIAS::HomgPoint3D G(QC), H(QC), I(QC);

 
  CovMatrix3x3 Cov;
  Cov.SetIdentity();
  Cov[0][0] = 100;
  Cov[1][1] = 64;

  Matrix3x3<double> VT, tmp;
#ifdef BIASASSERT_ISACTIVE
  int resPCA =
#endif
  Cov.GetPCA(S, VT);
  BIASASSERT(resPCA == 0);
  cout <<"Singular values of Cov are "<<S<<endl;
  
  tmp = VT;
  tmp.GetInverse(VT);

  // get main axes in wcs
  G_0 = VT * G_0; G_0.Normalize();
  H_0 = VT * H_0; H_0.Normalize();
  I_0 = VT * I_0; I_0.Normalize();

    // VT = VT.Transpose();
  // BIAS::Vector3<double> D0,D1,D2;
  // D0 =  VT.GetRow(0); D0.Normalize();
  // D1 =  VT.GetRow(1); D1.Normalize();
  // D2 =  VT.GetRow(2); D2.Normalize();

  for (unsigned int l=0; l<3; l++) {
    G[l] +=  sqrt(S[0]) * G_0[l] * GAUSS3D_CONFIDENCE_99_PERCENT;
    H[l] +=  sqrt(S[1]) * H_0[l] * GAUSS3D_CONFIDENCE_99_PERCENT;
    I[l] +=  sqrt(S[2]) * I_0[l] * GAUSS3D_CONFIDENCE_99_PERCENT;
  }
      
  // G, H and I are now extremal points of the ellipsoid, i.e. they are
  // each of them is the farest point from the center in the main dirs


  // check if ellipsoid center projects into conic, MUST always be true
  Quadric3D TestQ(QC, Cov, GAUSS3D_CONFIDENCE_99_PERCENT );
  cout << "Quadric is "<< TestQ << endl;
      
  Conic2D TestC; TestC.SetQuadricProjection(TestQ, P_current);
  cout << "Conic is "<<(TestC.IsEllipse()?"an ":"no ")
       << "ellipse: "<< TestC << endl;

  BIAS::HomgPoint3D *pPoint = &QC;
  BIAS::HomgPoint2D Projection;
  for (int u=0; u<4; u++) {
    switch (u) {
    case 0:pPoint = &QC;break;
    case 1:pPoint = &G; break;
    case 2:pPoint = &H; break;
    case 3:pPoint = &I; break;
    default:BIASERR("Debug loop out of range !");
    }
    Projection  = P_current * (*pPoint);
    Projection.Homogenize();

    double Location = TestQ.LocatePoint(*pPoint);
    cout <<"Computing point on quadric ! Location is "<<Location
         <<", pointindex[u] is "<<u<<endl;;
     
    Location = TestC.LocatePoint(Projection);
    if (Location>DBL_EPSILON) {
      BIASERR("Covariance ellipsoid projection: Point["<<u<< "] not in"
              <<" projection area !"<<endl<<(*pPoint)<<" -> "
              << Projection<<endl<<" Location is "<<Location
              <<", Conic is "<<TestC
              <<", Covariance matrix is "<< Cov 
              <<" sigmas are "<<S[0]<<" "<<S[1]<<" "<<S[2]<<endl
              <<"Ellipsoid center is at "<< QC);
      if (TestC.IsEllipse()) {
        BIAS::HomgPoint2D Center;
        double dAngle, radius_a, radius_b;
        if (TestC.GetEllipseParameters(Center,dAngle,radius_a,radius_b)) {
          BIASERR("Ellipse at "<<Center<<" with "<<radius_a<<"/"<<radius_b
                  <<" rotated by "<<dAngle/M_PI*180.0
                  <<" degrees counterclockwise.");
        } else {
          BIASERR("Error retrieving explicit ellipse parameters.");
        }
      } else {
        BIASERR("Projection is NO ellipse !");
      }
    } 
  }


  return 0;
}