ProjectionParametersBase.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 <iostream>
#include <sstream>
#include "Geometry/ProjectionParametersBase.hh"
#include <Geometry/SphericalCoordinates.hh>

using namespace std;
using namespace BIAS;
using namespace std;

#define PROJECTION_STREAM_PRECISION 12

ProjectionParametersBase::
ProjectionParametersBase(const unsigned int width,
                         const unsigned int height)
{
  width_=width;
  height_=height;
  principalX_=0;
  principalY_=0;
  // set default principal point as image center
  if (width > 0)  principalX_ = 0.5*double(width-1);
  if (height > 0) principalY_ = 0.5*double(height-1);
  aspectratio_=1;
  // init with some value:
  ustTransformIntoImage_ = false;
  ustIgnoreDistortion_ = false;
  ustNormalize_ = false;
  QValid_ = CValid_ = false;
  identifier_ = "";
  videoSourceType_ = "";
  // use sigma distance for unscented transform sampling
  // it is not clear whether this is correct, because alpha is not documented
  // in UnscentedTransform and also not in the reference paper !
  SetAlpha(1.0);
  
  //set values to pose as start
  //  SetC(Vector3<double>(0,0,0)); 
  //  SetQ(Quaternion<double>(0,0,0,1));

}

void ProjectionParametersBase::
GetFirstBorderPixel(PixelIterator& it) 
{
  it.x = 0;
  it.y = 0;
}

bool ProjectionParametersBase::
GetNextBorderPixel(PixelIterator& it) 
{
  if (it.x==0) {
    if (it.y==0) {
      // upper left image pixel, move right
      it.x++;
      return true;
    }
    // on way back up at left image border
    it.y--;
    // return false if finished
    return (it.y>0.0);
  } else if (it.x<width_-1) {
    if (it.y==0) {
      // moving left on upper image boundary
      it.x++;
      return true;
    }
    // moving right on lower image boundary
    it.x--;
    return true;
  } else {
    if (it.y<height_-1) {
      // moving down on right image side
      it.y++;return true;
    }
    // lower right corner
    it.x--;
    return true;
  }
}

void ProjectionParametersBase::
GetFirstEdgePosition(PixelIterator& it) 
{
  it.x = -0.5;
  it.y = -0.5;
}

bool ProjectionParametersBase::
GetNextEdgePosition(PixelIterator& it) 
{
  if (it.x==-0.5) {
    if (it.y==-0.5) {
      // upper left image pixel, move right
      it.x++;
      return true;
    }
    // on way back up at left image border
    it.y--;
    // return false if finished
    return (it.y>-0.5);
  } else if (it.x<width_-0.5) {
    if (it.y==-0.5) {
      // moving left on upper image boundary
      it.x++;
      return true;
    }
    // moving right on lower image boundary
    it.x--;
    return true;
  } else {
    if (it.y<height_-0.5) {
      // moving down on right image side
      it.y++;return true;
    }
    // lower right corner
    it.x--;
    return true;
  }
}


int ProjectionParametersBase::
GetSphericalViewingRange(const CoordinateTransform3D& sphericalReferenceFrame,
                         double& minPhi, double& maxPhi, double& centerPhi,
                         double& minTheta, double& maxTheta )
{

  CoordinateTransform3D localSphericalReferenceFrame;
  localSphericalReferenceFrame.BecomeRelativeTransform(sphericalReferenceFrame,
                                                       Pose_);

  if(localSphericalReferenceFrame.GetC().NormL2() > 1e-12) {
    BIASERR("Spherical frame position and camera position are different!"
            << endl << localSphericalReferenceFrame.GetC().NormL2());
    return -1;
  }

  SphericalCoordinates euc2Sphere;
  euc2Sphere.SetAffineBase(localSphericalReferenceFrame);

  minPhi = M_PI;
  maxPhi = -3.1;
  minTheta = M_PI;
  maxTheta = 0;

  double phi, theta, rho;
  HomgPoint2D imgPoint;
  imgPoint[2] = 1.0;
  HomgPoint3D ray;
  PixelIterator it;
  //  GetFirstBorderPixel(it);
  GetFirstEdgePosition(it);
  Vector3<double> p, d;
  do {
    imgPoint[0] = it.x;
    imgPoint[1] = it.y;

    UnProjectLocal(imgPoint, p, d);
    ray.Set(d);
    euc2Sphere.GetSphericalCoordinatesFullPhi(ray, rho, phi, theta);
    if(phi<minPhi) minPhi = phi;
    if(phi>maxPhi) maxPhi = phi;
    if(theta<minTheta) minTheta = theta;
    if(theta>maxTheta) maxTheta = theta;
  }
  //  while(GetNextBorderPixel(it));
  while(GetNextEdgePosition(it));
  //calculate centerPhi:
  imgPoint[0] = (width_-1)/2;
  imgPoint[1] = (height_-1)/2;
  UnProjectLocal(imgPoint, p, d);
  ray.Set(d);
  euc2Sphere.GetSphericalCoordinatesFullPhi(ray, rho, centerPhi, theta);


  //determine whether theta=0 or theta = M_PI do
  //project into the images in order to find minimal/maximal
  // theta angles:
  HomgPoint3D minRay, maxRay;
  euc2Sphere.GetCartesianRayFromFullPhi(0, 0, minRay);
  euc2Sphere.GetCartesianRayFromFullPhi(0, M_PI, maxRay);
  HomgPoint2D tmpImgPoint;
  bool minimumIsIn = DoesPointProjectIntoImageLocal(minRay.GetEuclidean(),
                                                    tmpImgPoint);
  bool maximumIsIn = DoesPointProjectIntoImageLocal(maxRay.GetEuclidean(),
                                                    tmpImgPoint);
  if(maximumIsIn) {
    maxTheta = M_PI;
  }
  if(minimumIsIn) {
    minTheta = 0;
  }

  return 0;
}

bool ProjectionParametersBase::
DoesPointProjectIntoImage(const HomgPoint3D& X,
                          HomgPoint2D& x,
                          bool IgnoreDistortion) const
{
  HomgPoint3D tempX = GetPose().GlobalToLocal(X);
  BIAS::Vector3<double> localX;
  tempX.GetEuclidean(localX);
  return DoesPointProjectIntoImageLocal(localX, x, IgnoreDistortion);
}

const bool ProjectionParametersBase::
DoIntrinsicsDiffer(const ProjectionParametersBase* p) const
{
/* 
  cout << "principalX_ : " << principalX_ << " <-> " << p->principalX_ << endl;
  cout << "principalY_ : " << principalY_ << " <-> " << p->principalY_ << endl;
  cout << "width_ : " << width_ << " <-> " << p->width_ << endl;
  cout << "height_ : " << height_ << " <-> " << p->height_ << endl;
  cout << "aspectratio_ : " << aspectratio_ << " <-> " << p->aspectratio_ << endl;
*/
  if(principalX_ != p->principalX_)
    return true;
  if(principalY_ != p->principalY_)
    return true;
  if(width_ != p->width_)
    return true;
  if(height_ != p->height_)
    return true;
  if(aspectratio_ != p->aspectratio_)
    return true;
  return false;
}

const bool ProjectionParametersBase::
DoExtrinsicsDiffer(const ProjectionParametersBase* p) const
{
/*  
  cout << "CValid_ : " << CValid_ << " <-> " << p->CValid_ << endl;
  cout << "QValid_ : " << QValid_ << " <-> " << p->QValid_ << endl;
  cout << "C : " << GetC() << " <-> " << p->GetC()
       << " (diff = " << (GetC()-p->GetC()).NormL2() << ")" << endl;
  cout << "R : " << GetR() << " <-> " << p->GetR()
       << " (diff = " << (GetR()-p->GetR()).NormL2() << ")" << endl;
  cout << "Q : " << GetQ() << " <-> " << p->GetQ()
       << " (diff = " << (GetQ()-p->GetQ()).NormL2() << ")" << endl;
*/  
  if ((p->CValid_ != CValid_) || (p->QValid_ != QValid_))
    return true;
  if (QValid_ && !BIAS::Equal((GetR()-p->GetR()).NormL2(), 0.0, 1e-9))
    //           GetR() != p->GetR())
    return true;
  if (CValid_ && !BIAS::Equal((GetC()-p->GetC()).NormL2(), 0.0, 1e-9))
    //           GetC() != p->GetC())
    return true;
  return false;
}


RMatrix ProjectionParametersBase::
GetR() const
{
  //if(!QValid_) BIASWARN("Projection is invalid!");
  RMatrix R;
  // unconst, because RMatrix::SetFromQ() in not const enough
  Quaternion<double> qtmp = GetQ();
  R.SetFromQuaternion(qtmp);
  return R;
}


#ifdef BIAS_HAVE_XML2

int ProjectionParametersBase::
XMLGetClassName(std::string& TopLevelTag,
                double& Version) const {
  TopLevelTag = "ProjectionParametersBase";
  Version = 0.2;
  return 0;
}

int ProjectionParametersBase::
XMLOut(const xmlNodePtr Node, XMLIO& XMLObject) const {

  xmlNodePtr childNode;
  childNode = XMLObject.addChildNode(Node, "Identifier");
  string identutf8;
  XMLIO::IsoLatin1ToUtf8( GetIdentifier(),identutf8); 
  XMLObject.addAttribute(childNode, "val", identutf8);

  childNode = XMLObject.addChildNode(Node, "VideoSourceType");
  XMLObject.addAttribute(childNode, "val", GetVideoSourceType());

  childNode = XMLObject.addChildNode(Node, "ImageSize");
  XMLObject.addAttribute(childNode, "width", (int)width_);
  XMLObject.addAttribute(childNode, "height", (int)height_);

  childNode = XMLObject.addChildNode(Node, "PrincipalPoint");
  XMLObject.addAttribute(childNode, "x", principalX_);
  XMLObject.addAttribute(childNode, "y", principalY_);

  childNode = XMLObject.addChildNode(Node, "AspectRatio");
  XMLObject.addAttribute(childNode, "val", aspectratio_);

  // R
  //   Quaternion<double> Q(0.0, 0.0, 0.0, 1.0);
  //   if (PoseValid()){
  //       Q = GetQ();
  //   }
  if (QValid_) {
    Quaternion<double> Q = GetQ();
    stringstream streamR;
    childNode = XMLObject.addChildNode(Node,"Rotation");
    streamR << std::setprecision(PROJECTION_STREAM_PRECISION);
    for (unsigned int i = 0; i < 4; i++) {
      streamR << (double) Q.GetData()[i] <<" ";
    }
    XMLObject.addContent(childNode, streamR.str());
  }

  // C
  //   Vector3<double> C(0.0, 0.0, 0.0);
  //   if (PoseValid()){
  //     C = GetC();
  //   }
  if (CValid_){
    Vector3<double> C = GetC();
    childNode = XMLObject.addChildNode(Node,"Center");
    XMLObject.addAttribute(childNode, "x", C[0]);
    XMLObject.addAttribute(childNode, "y", C[1]);
    XMLObject.addAttribute(childNode, "z", C[2]);
  }
  return 0;
}

int ProjectionParametersBase::
XMLIn(const xmlNodePtr Node, XMLIO& XMLObject) {

  xmlNodePtr childNode;
  if ((childNode = XMLObject.getChild(Node, "Identifier"))!=NULL) {
    string identutf8;
    identutf8 = XMLObject.getAttributeValueString(childNode, "val");
    XMLIO::Utf8ToIsoLatin1(identutf8,identifier_);
  }
 
  if ((childNode = XMLObject.getChild(Node, "VideoSourceType"))!=NULL) {
    videoSourceType_ = XMLObject.getAttributeValueString(childNode, "val");
  }
 
  if ((childNode = XMLObject.getChild(Node, "ImageSize"))!=NULL) {
    width_ = XMLObject.getAttributeValueInt(childNode, "width");
    height_ = XMLObject.getAttributeValueInt(childNode, "height");
  }
  if ((childNode = XMLObject.getChild(Node, "PrincipalPoint"))!=NULL) {
    principalX_ = XMLObject.getAttributeValueDouble(childNode, "x");
    principalY_ = XMLObject.getAttributeValueDouble(childNode, "y");
  }
  if ((childNode = XMLObject.getChild(Node, "CellSize"))!=NULL) {
    double cellSizeX, cellSizeY;
    cellSizeX = XMLObject.getAttributeValueDouble(childNode, "x");
    cellSizeY = XMLObject.getAttributeValueDouble(childNode, "y");
    aspectratio_ = cellSizeX / cellSizeY;
  }
  if ((childNode = XMLObject.getChild(Node, "AspectRatio"))!=NULL) {
    aspectratio_ = XMLObject.getAttributeValueDouble(childNode, "val");
  }

  Quaternion<double> Q;
  // R
  if ((childNode = XMLObject.getChild(Node, "Rotation"))!=NULL) {
    // get data fields
    stringstream ssRotation(XMLObject.getNodeContentString(childNode));
    unsigned int CurrentData = 0;
    while (ssRotation.good() && CurrentData<4) {
      ssRotation >> (Q.GetData()[CurrentData++]);
    }
    if (CurrentData<4) {
      // reached end of stream before expected end of data
      BIASERR("error in xml structure, rotation matrix incomplete.");
      return -1;
    }    
    SetQ(Q);
  }
  else
    QValid_ = false;

  // C
  Vector3<double> C;
  if ((childNode = XMLObject.getChild(Node, "Center"))!=NULL) {
    C[0] = XMLObject.getAttributeValueDouble(childNode, "x");
    C[1] = XMLObject.getAttributeValueDouble(childNode, "y");
    C[2] = XMLObject.getAttributeValueDouble(childNode, "z");
    SetC(C);
  }
  else
    CValid_ = false;
  return 0;
}
#endif

double ProjectionParametersBase::
ViewDifference(const ProjectionParametersBase*) const {
  // there should be an implmentation in each derived class
  BIASERR("ViewDifference not implemented in derived class");
  BIASABORT;
  // visual studio needs this
  return 0.0;
}

Vector3<double> ProjectionParametersBase::
UnProjectToPointLocal(const HomgPoint2D& pos,
                        const double& depth, bool IgnoreDistortion) const {
    Vector3<double> x, p;
    UnProjectLocal(pos, p, x, IgnoreDistortion);
    x *= depth;
#ifdef BIAS_DEBUG
  if (fabs(x.NormL2()-depth)>1e-10) {
    BIASERR("Unprojection failed:" << x << "  Position:"<< pos
             <<"  Norm:" << setw(10) << double(x.NormL2()));
    // BIASABORT;
  }
#endif
  return x;
}

HomgPoint3D ProjectionParametersBase::
UnProjectToPoint(const HomgPoint2D& pos,
                 double depth,
                 bool IgnoreDistortion) const
{
  // cout<<"[ProjectionParametersBase] UnprojectToPoint pos:"<<pos<<endl;
  Vector3<double> point;
  GetR().Mult(UnProjectToPointLocal(pos, depth, IgnoreDistortion), point);
  point.AddIP(GetC());
  return HomgPoint3D(point);
}

HomgPoint3D ProjectionParametersBase::
UnProjectToPoint(const HomgPoint2D& pos,
                 const double depth,
                 const ProjectionParametersBase &proj,
                 bool IgnoreDistortion) const
{
  Vector3<double> point;
  GetR().Mult(proj.UnProjectToPointLocal(pos, depth, IgnoreDistortion), point);
  point.AddIP(GetC());
  return HomgPoint3D(point);
}

int ProjectionParametersBase::
Transform_(const Vector<double>& src, Vector<double>& dst) const {
  BIASASSERT(src.size()==3);BIASASSERT(dst.size()==3);
  HomgPoint2D xsrc(src[0], src[1], src[2]), xdst, p;
  if (ustTransformIntoImage_) {
    xdst = ProjectLocal(xsrc, ustIgnoreDistortion_);
  } else {
    UnProjectLocal(xsrc, p, xdst, ustIgnoreDistortion_);
  }
  if (xdst.NormL2()<1e-10) {
    // point does not transform, failed !
    dst[0] = dst[1] = dst[2] = 0; // FIXME TODO JW
    return -1;
  }
  // we do a homogenize here, for the spherical case the covariance
  // may be interesting on the normalized point.
  // However, for < 160 degrees, homogenize should be no problem also for
  // spherical cameras.
  if (ustNormalize_) {
    xdst.Normalize();
  } else {
    xdst.Homogenize();
  }
  dst[0] = xdst[0];
  dst[1] = xdst[1];
  dst[2] = xdst[2];
  return 0;
}


Matrix3x3<double> ProjectionParametersBase::
UnProjectCovLocal(const HomgPoint2D& pos, const Matrix3x3<double>& cov2D,
               bool IgnoreDistortion, bool Normalize)
{
  // set transformation parameters using class members
  ustIgnoreDistortion_ = IgnoreDistortion;
  ustTransformIntoImage_ = false;
  ustNormalize_ = Normalize;
  // convert into general data structs
  Matrix<double> resCov(3,3,MatrixZero);
  Vector<double> srcPoint(3);
  Vector<double> resPoint(3);
  srcPoint[0] = pos[0];
  srcPoint[1] = pos[1];
  srcPoint[2] = pos[2];
  // use unscented transform to approximate covariance
  Transform(srcPoint, (Matrix<double>)cov2D, resPoint, resCov);
  //cout << "Original cov = " << cov2D << " around " << pos << endl
  //     << "Transformed to " << resCov << " around " << resPoint << endl;
  return (Matrix3x3<double>)resCov;
}


ProjectionParametersBase::~ProjectionParametersBase(){
#ifdef BIAS_DEBUG
  // mark this object as invalid
  width_ = 999999;
  height_ = 999999;
  principalX_ = 999999;
  principalY_ = 999999;
  aspectratio_ = 999999;
#endif
}


KMatrix ProjectionParametersBase::
GetFakeKMatrix(int resolution, const double& maxangle)
  const {
  double dum;
  return GetFakeKMatrix(dum, resolution);
}

KMatrix ProjectionParametersBase::
GetFakeKMatrix(double& imgsize, int resolution, const double& ang) const {
  imgsize = 0;


  Vector3<double> ray1, ray2;
  HomgPoint2D imgCoord1, imgCoord2;
  switch (resolution) {
  case 1:
    // measure 1 deg angle in camera center
    ray1.Set(sin(ang-M_PI/180.0),0.0,cos(ang-M_PI/180.0));
    ray2.Set(sin(ang),0.0,cos(ang));
    imgCoord1 = ProjectLocal(ray1);
    imgCoord2 = ProjectLocal(ray2);
    imgsize=(imgCoord2[0]-imgCoord1[0])*
      tan(ang)/tan(M_PI/180.0);
    break;
  case 2:
    // measure 1 deg angle near max angle
    ray1.Set(sin(ang-M_PI/180.0),0.0,cos(ang-M_PI/180.0));
    ray2.Set(sin(ang),0.0,cos(ang));
    imgCoord1 = ProjectLocal(ray1);
    imgCoord2 = ProjectLocal(ray2);
    imgsize=(imgCoord2[0]-imgCoord1[0])*
      tan(ang)/(tan(ang)-tan(ang-M_PI/180.0));
    break;
  default:
    ray1.Set(sin(ang),0.0,cos(ang));
    imgCoord1 = ProjectLocal(ray1);
    imgsize = imgCoord1[0]-principalX_;
  }
  KMatrix K;
  K.SetIdentity();
  K[0][0] = imgsize/tan(ang);
  K[1][1] = imgsize/tan(ang);
  K[0][2] = imgsize;
  K[1][2] = imgsize;
  imgsize *= 2;
  return K;
}


void ProjectionParametersBase::
Rescale(double ratio, const double offset)
{
  double widthTemp = (double)width_/ratio;
  double heightTemp = (double)height_/ratio;

  /*
  cout << "PPB::Rescale( ratio=" << setprecision(10) << ratio << ", offset=" << offset << "): "
       << "width_ x height_ = (" << width_ << " x " << height_ << "), "
       << "widthTemp x heightTemp = (" << widthTemp << " x " << heightTemp << ")\n";
  */
  width_  = (int)rint(widthTemp); //was floor... which had the undesired sideeffect
                                  //that 639.999999999 became 639 pixels instead
                                  //of 640
  height_ =  (int)rint(heightTemp);
  //cout << "new width_ x height_ = (" << width_ << " x " << height_ << ")\n";

  if( (!BIAS::Equal(widthTemp, (double)width_)) ||
      (!BIAS::Equal(heightTemp, (double)height_))) {
    //         BIASERR("image size has integer type."<<width_<<"x"<<height_
    //                 <<" Passed rescale factor "<< ratio <<" does not conform!\n");
  }
  principalX_ = (principalX_ - offset)/ratio;
  principalY_ = (principalY_ - offset)/ratio;

}



void ProjectionParametersBase::
Rescale(unsigned int width, unsigned int height)
{

  unsigned int oldWidth = width_;
  unsigned int oldHeight = height_;

  width_ = width;
  height_= height;

  principalX_ = principalX_ / ((double)oldWidth/(double)width_);
  principalY_ = principalY_ / ((double)oldHeight/(double)height_);
}


#define JACOBIAN_EPSILON 0.5

int ProjectionParametersBase::
GetUnProjectionJacobian(const HomgPoint2D& x,
                        Matrix2x2<double>& Jac,
                        const bool homogenized) const {
  BIASASSERT(x.IsHomogenized());
  BIASASSERT(homogenized);// normalized not yet implemented
  // numeric approximation in base class, make better in derived classes

  // compute rays for left/right/top/bottom pixels
  Vector3<double> L = UnProjectToPointLocal(HomgPoint2D(x[0]-JACOBIAN_EPSILON,
                                                        x[1]), 1.0);
  if (L.NormL2()<0.5) return -1;
  Vector3<double> R = UnProjectToPointLocal(HomgPoint2D(x[0]+JACOBIAN_EPSILON,
                                                        x[1]), 1.0);

  if (R.NormL2()<0.5) return -1;


  Vector3<double> T = UnProjectToPointLocal(HomgPoint2D(x[0],
                                                        x[1]-JACOBIAN_EPSILON),
                                            1.0);
  if (T.NormL2()<0.5) return -1;
  Vector3<double> B = UnProjectToPointLocal(HomgPoint2D(x[0],
                                                        x[1]+JACOBIAN_EPSILON),
                                            1.0);

  if (B.NormL2()<0.5) return -1;

  //cout<<"Left   is "<<L<<endl;
  //cout<<"Righ   is "<<R<<endl;
  //cout<<"Top    is "<<T<<endl;
  //cout<<"Bottom is "<<B<<endl;

  // check if caller wants tangent plane instead of w=1 plane
  if (!homogenized) {
    // parallel projection into tangent plane, "rotate camera"
    Vector3<double> C = UnProjectToPointLocal(x, 1.0);
    if (C.NormL2()<0.5) return -1;
    //cout<<"feature ray C is "<<C<<endl;
    // here we are faced with the problem of selecting a good x- and y-axis
    // on the tangent plane
    RMatrix Rot(MatrixIdentity);
    // use shortest way on the sphere: rotate around orthogonal axis
    Vector3<double> axis(C.CrossProduct(Vector3<double>(0,0,1)));
    if (axis.NormL2()>1e-7) {
      axis.Normalize();
      Rot.Set(axis, acos(C[2]));
      L = Rot * L;
      R = Rot * R;
      T = Rot * T;
      B = Rot * B;
    }
    //cout<<"new L is "<<L<<endl;
    //cout<<"new R is "<<R<<endl;
    //cout<<"new T is "<<T<<endl;
    //cout<<"new B is "<<B<<endl;
    //cout<<"Rot is "<<Rot<<endl;
    //cout<<"Rot*C="<<Rot*C<<endl;
    BIASASSERT((Rot*C-Vector3<double>(0,0,1)).NormL2()<1e-7);
  } else {
    // homogenize
    L/=L[2];
    R/=R[2];
    T/=T[2];
    B/=B[2];
  }

  // R-L: if we make one step in x direction of the image,
  // what step will we make in the space plane? f(R)-f(L)
  // set into first row of Jac
  Jac[0][0] = (R[0]-L[0]); //   /2.0*JACOBIAN_EPSILON==1
  Jac[1][0] = (R[1]-L[1]); //   /2.0*JACOBIAN_EPSILON==1
  // same with vertical step
  Jac[0][1] = (B[0]-T[0]); //   /2.0*JACOBIAN_EPSILON==1
  Jac[1][1] = (B[1]-T[1]); //   /2.0*JACOBIAN_EPSILON==1

  return 0;
}

bool ProjectionParametersBase::
IsLeftOf(const ProjectionParametersBase& ppb) const
{
  const Vector3<double> baseline = ppb.GetC() - GetC();
  Matrix3x3<double> rMatrixInv;
  GetR().GetInverse(rMatrixInv);
  const Vector3<double> transformedBaseline = rMatrixInv * baseline;

  return transformedBaseline[0] >= 0.0;
}

void ProjectionParametersBase::
SetVideoSourceType(const std::string &name)
{
  videoSourceType_ = name;
}