HomgPlane3D.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*/
#include "HomgPlane3D.hh"

#include <Base/Common/CompareFloatingPoint.hh>

using namespace BIAS;
using namespace std;

void HomgPlane3D::
Set(const HomgPoint3D& PointOnPlane, const Vector3<double>& NormalVector)
{
  HomgPoint3D poph(PointOnPlane); // de-const
  // normalise for numerical stability
  // zero point is not a valid point in projective space
  BIASASSERT(!Equal(poph.NormL2(), 0.0)); 
  poph /= poph.NormL2();

  // normalise for numerical stability, only the direction of the 
  // normal vector is important, not its length
  Vector3<double> normal(NormalVector);
  BIASASSERT(!Equal(normal.NormL2(), 0.0)); 
  normal /= normal.NormL2();

  // the plane with normal parallel to NormalVector is given by 
  //   [ NormalVector' p4 ]'
  // since the plane is an entity in projective space, the same plane is
  // described by 
  //   lambda * [ NormalVector' p4 ]' 
  // for each lambda != 0
  // The PointOnPlane only resides on the plane when the scalar product
  // between point and plane vanishes
  //   plane' * PointOnPlane = 0
  // or alternatively
  //   pop' * lambda * NormalVector + w * lambda * p4 = 0 (i)
  // with PointOnPlane = [ pop' w ]'
  Vector3<double> pop(poph[0], poph[1], poph[2]);
  double sp = normal.ScalarProduct(pop);
  if ( fabs(poph[3]) > fabs(sp) ){ 
    // setting lambda = 1, equation (i) can easily be resolved
    //    p4 = - pop' * NormalVector / w
    data_[0] = normal[0]; 
    data_[1] = normal[1]; 
    data_[2] = normal[2]; 
    data_[3] = -sp / poph[3];
  } else {
    // setting lambda * p4 = 1
    // results in 
    //    lambda = -w /  pop' * NormalVector
    double lambda = -poph[3] / sp;
    data_[0] = lambda * normal[0]; 
    data_[1] = lambda * normal[1]; 
    data_[2] = lambda * normal[2]; 
    data_[3] = 1.0;
  }
}