CylindricalRectification.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 <Image/CylindricalRectification.hh>
#include <Geometry/ProjectionParametersSpherical.hh>
#include <Geometry/ProjectionParametersPerspective.hh>
#include <Geometry/ProjectionParametersCylindric.hh>
#include <Base/Geometry/HomgPoint3D.hh>
#include <Base/Common/CompareFloatingPoint.hh>
#include <Base/ImageUtils/BresenhamCircle.hh>

using namespace BIAS;
using namespace std;

#define SUPERPARENT RectificationBase<InputStorageType, OutputStorageType>
#define PARENT RectificationViaProjectionMappingBase<InputStorageType, OutputStorageType>


template <class InputStorageType, class OutputStorageType>
BIAS::CylindricalRectification<InputStorageType, OutputStorageType>::
CylindricalRectification(double maxAngle) :
PARENT::RectificationViaProjectionMappingBase("CylindricalRectification"),
maxHalfAngleInX_(maxAngle)
{ }

template <class InputStorageType, class OutputStorageType>
BIAS::CylindricalRectification<InputStorageType, OutputStorageType>::
~CylindricalRectification()
{ }

template <class InputStorageType, class OutputStorageType>
bool BIAS::CylindricalRectification<InputStorageType, OutputStorageType>::
IsInputCameraValid(const BIAS::Image<InputStorageType>& img,
                   const BIAS::ProjectionParametersBase* proj)
{
  const ProjectionParametersPerspective* ppp =
    dynamic_cast<const ProjectionParametersPerspective*>(proj);
  const ProjectionParametersSpherical* pps =
    dynamic_cast<const ProjectionParametersSpherical*>(proj);
  const ProjectionParametersCylindric* ppc =
    dynamic_cast<const ProjectionParametersCylindric*>(proj);

  if (ppp == NULL && pps == NULL && ppc == NULL) {
    return false;
  }

  return true;
}


template <class InputStorageType, class OutputStorageType>
int BIAS::CylindricalRectification<InputStorageType, OutputStorageType>::
DetermineRectificationParameters_()
{

  Pose rectPoseA, rectPoseB;//, relativeTransformA, relativeTransformB;

  if (SUPERPARENT::CalculateRectifiedBases(SUPERPARENT::ppBA_->GetPose(),
                                           SUPERPARENT::ppBB_->GetPose(),
                                           rectPoseA, rectPoseB) != 0) {
    BIASERR("Failed calculating rectified base");
    return -1;
  }
  //  relativeTransformA.BecomeRelativeTransform(cameraA->GetPose(), rectPoseA);
  //  relativeTransformB.BecomeRelativeTransform(cameraB->GetPose(), rectPoseB);

/**
  Determine the camera type and the needed resolutions!
  If perspective determine the min max angles of both cameras in phi and
  theta by projection onto the cylinder. Choose minimal phi range and maximal
  x range. Prohibit theta larger then maxTheta.
*/

  double minX, maxX;
  double minPhi, maxPhi;
  double xMinA, xMinB;
  double xMaxA, xMaxB;
  double phiMinA, phiMinB;
  double phiMaxA, phiMaxB;
  //samples per unit length
  double relativeXResA, relativeXResB;
  //samples per radians
  double relativePhiResA, relativePhiResB;

  if (DetermineCylindricCameraBoundries_(SUPERPARENT::ppBA_,
                                         rectPoseA,
                                         xMinA, xMaxA, phiMinA, phiMaxA,
                                         relativeXResA, relativePhiResA) != 0) {
    BIASERR("Error calculating rectification boundaries.\n");
    return -1;
  }

  if (DetermineCylindricCameraBoundries_(SUPERPARENT::ppBB_,
                                         rectPoseB,
                                         xMinB, xMaxB, phiMinB, phiMaxB,
                                         relativeXResB, relativePhiResB) != 0) {
    BIASERR("Error calculating rectification boundaries.\n");
    return -1;
  }

  //check for overlap
  //  exploit that ordering between minA and maxA respectively minB and maxB
  //  is guaranteed
  if(xMinB > xMaxA || xMinA > xMaxB) {
    BIASERR("it seems that cameras do not have a viewing overlap within epipolar planes!");
    return -1;
  }
  if(phiMinB > phiMaxA || phiMinA > phiMaxB) {
    BIASERR("it seems that cameras do not have shared epipolar lines!");
    return -1;
  }

  //final values
  //the angle for phi and theta has to be used
  if (xMinA > xMinB) minX = xMinB;
  else minX = xMinA;
  if (xMaxA > xMaxB) maxX = xMaxA;
  else maxX = xMaxB;

  // do not use extremal angles here because these angles
  //represent the common epipolar line range
  if(phiMinA > phiMinB) minPhi = phiMinA;
  else minPhi = phiMinB;
  if(phiMaxA < phiMaxB) maxPhi = phiMaxB;
  else maxPhi = phiMaxA;

  double relativeXRes, relativePhiRes;
  if(relativeXResA > relativeXResB) relativeXRes = relativeXResA;
  else relativeXRes = relativeXResB;
  if(relativePhiResA > relativePhiResB) relativePhiRes = relativePhiResA;
  else relativePhiRes = relativePhiResB;

  unsigned int xRes = (unsigned int) ceil(relativeXRes*(maxX - minX));
  unsigned int phiRes = (unsigned int) ceil(relativePhiRes*(maxPhi - minPhi));

  delete PARENT::rectPPA_;
  delete PARENT::rectPPB_;
  PARENT::rectPPA_ = new ProjectionParametersCylindric(minX, maxX,
                                                       minPhi, maxPhi,
                                                       xRes, phiRes);
  PARENT::rectPPB_ = new ProjectionParametersCylindric(minX, maxX,
                                                       minPhi, maxPhi,
                                                       xRes, phiRes);
  PARENT::rectPPA_->SetPose(rectPoseA);
  PARENT::rectPPB_->SetPose(rectPoseB);
/*
#ifdef BIAS_DEBUG
  std::cout<<"orig A:\n";
  std::cout<<*SUPERPARENT::ppBA_;
  std::cout<<"rect A:\n";
  PARENT::rectPPA_->Print();
  std::cout<<"orig B:\n";
  std::cout<<*SUPERPARENT::ppBB_;
  std::cout<<"rect B:\n";
  PARENT::rectPPB_->Print();
#endif
*/

  return 0;
}


/** Helper for DetermineCylindricCameraBoundries_.
 *  Finds the corresponding x and phi coordinates for optical ray and
 *  returns results by rescaling the passed point and updating passed phi.
 *  If point lies on cylinder axis -1 is returned.
 *
 */
template <class InputStorageType, class OutputStorageType>
int BIAS::CylindricalRectification<InputStorageType, OutputStorageType>::
ProjectOntoUnitCylinder_(BIAS::HomgPoint3D& point, double& phi)
{
  double radius= sqrt(point[1]*point[1] + point[2]*point[2]);
  if(Equal(radius, 0.0)) {
    phi = 0.0;
    return -1;
  } else {
    double scale = 1.0 / sqrt(point[1]*point[1] + point[2]*point[2]);
    point *= scale;
    phi = atan2(point[1], point[2]);
  }
  return 0;
}

/** Helper for DetermineCylindricCameraBoundries_.
 *  Finds the corresponding x and phi coordinates for optical ray
 *  and updates min/max values accordingly.
 */
template <class InputStorageType, class OutputStorageType>
void BIAS::CylindricalRectification<InputStorageType, OutputStorageType>::
SetMinMaxFrom_(HomgPoint3D& rectPoint,
                   double& minX, double& maxX,
                   double& minPhi, double& maxPhi)
{
  double phi;
  if(ProjectOntoUnitCylinder_(rectPoint, phi)==0) {
    if(rectPoint[0]>maxX) maxX = rectPoint[0];
    if(rectPoint[0]<minX) minX = rectPoint[0];
    if(phi>maxPhi) maxPhi = phi;
    if(phi<minPhi) minPhi = phi;
  }
}

template <class InputStorageType, class OutputStorageType>
int BIAS::CylindricalRectification<InputStorageType, OutputStorageType>::
DetermineCylindricCameraBoundriesPerspective_(
                                 const ProjectionParametersBase* camera,
                                 const Pose& rectPose,
                                 double& minX, double& maxX,
                                 double& minPhi, double& maxPhi,
                                 double& relativeXRes,
                                 double& relativePhiRes) const
{
  //determine bounding volume on cylinder
  Pose relativeTransform;
  //transforms points from camera coords (local) to cylinder coords (global)
  relativeTransform.BecomeRelativeTransform(camera->GetPose(), rectPose);
  unsigned int width, height;
  camera->GetImageSize(width, height);

  maxX = DBL_MIN;
  minX = DBL_MAX;
  maxPhi = -M_PI;
  minPhi = M_PI;

  Vector3<double> testedPoint;
  HomgPoint3D rectPoint;
  //  double phi;
  //check upper and lower image border
  Vector3<double> p;
  for(unsigned int x=0; x<width; x++) {
    camera->UnProjectLocal(HomgPoint2D(x, 0, 1.0), p, testedPoint);
    if(testedPoint.NormL2() < 1e-10) BIASABORT;
    rectPoint = relativeTransform.LocalToGlobal(HomgPoint3D(testedPoint));
    SetMinMaxFrom_(rectPoint, minX, maxX, minPhi, maxPhi);

    camera->UnProjectLocal(HomgPoint2D(x, height-1, 1.0), p, testedPoint);
    if(testedPoint.NormL2() < 1e-10) BIASABORT;
    rectPoint = relativeTransform.LocalToGlobal(HomgPoint3D(testedPoint));
    SetMinMaxFrom_(rectPoint, minX, maxX, minPhi, maxPhi);
  }

  for(unsigned int y=0; y<height; y++) {
    camera->UnProjectLocal(HomgPoint2D(0, y, 1.0), p, testedPoint);
    if(testedPoint.NormL2() < 1e-10) BIASABORT;
    rectPoint = relativeTransform.LocalToGlobal(HomgPoint3D(testedPoint));
    SetMinMaxFrom_(rectPoint, minX, maxX, minPhi, maxPhi);

    camera->UnProjectLocal(HomgPoint2D(width-1, y, 1.0), p, testedPoint);
    if(testedPoint.NormL2() < 1e-10) BIASABORT;
    rectPoint = relativeTransform.LocalToGlobal(HomgPoint3D(testedPoint));
    SetMinMaxFrom_(rectPoint, minX, maxX, minPhi, maxPhi);
  }

  double clippingMaxX = tan(maxHalfAngleInX_);
// #ifdef BIAS_DEBUG
//   cerr<<"detemined image spaceing:\n";
//   cerr<<"clipMinX/clipMaxX = "<<-clippingMaxX<<"/"<<clippingMaxX<<endl;
//   cerr<<"minX/maxX = "<<minX<<"/"<<maxX<<endl;
//   cerr<<"minPhi/maxPhi = "<<minPhi<<"/"<<maxPhi<<endl;
// #endif


  if(minX < -clippingMaxX) minX = -clippingMaxX;
  if(maxX > clippingMaxX) maxX = clippingMaxX;

  /* replaced by code below:
  //epipole is in image all angles are present
  HomgPoint3D epiDir = relativeTransform.GlobalToLocal(HomgPoint3D(1,0,0));
  HomgPoint2D temp;
  if(camera->DoesPointProjectIntoImage(epiDir, temp)) {
    maxPhi = M_PI;
    minPhi = -M_PI;
  }
  */

  /* find out if the camera is receiving optical rays
     along direction of cylinder axis: */
  HomgPoint3D epiDirHomg = relativeTransform.GlobalToLocal(HomgPoint3D(1,0,0));
  Vector3<double> epiDir;
  epiDirHomg.GetEuclidean(epiDir);
//   cerr<<"epiDir = "<<epiDir<<endl;
  HomgPoint3D epiDirNegHomg =
    relativeTransform.GlobalToLocal(HomgPoint3D(-1,0,0));
  Vector3<double> epiDirNeg;
  epiDirHomg.GetEuclidean(epiDirNeg);
//   cerr<<"epiDirNeg = "<<epiDirNeg<<endl;
  HomgPoint2D temp;



  if(camera->DoesPointProjectIntoImageLocal(epiDir, temp)) {
    maxX = clippingMaxX;//positive x is fully visible!
//      cerr<<"maximizing X...\n";
  }
  if(camera->DoesPointProjectIntoImageLocal(epiDirNeg, temp)) {
    minX = -clippingMaxX;//negative x is fully visible!
//     cerr<<"minimizing X...\n";
  }

  if(maxPhi<=minPhi)
    return -1;

  double deltPhi = maxPhi - minPhi;
  //heuristic number of epilines:
  relativePhiRes = (width+height)/deltPhi;
  //heuristic!
  double deltaX = maxX - minX;
  relativeXRes = sqrt( (double)(width*width+height*height ))/deltaX;
  return 0;
}

template <class InputStorageType, class OutputStorageType>
int BIAS::CylindricalRectification<InputStorageType, OutputStorageType>::
DetermineCylindricCameraBoundriesSpherical_(
                                 const ProjectionParametersBase* camera,
                                 const Pose& rectPose,
                                 double& minX, double& maxX,
                                 double& minPhi, double& maxPhi,
                                 double& relativeXRes,
                                 double& relativePhiRes) const
{
  const ProjectionParametersSpherical* pps =
    dynamic_cast<const ProjectionParametersSpherical*>(camera);

  //  BIASWARN("Spherical camera rectification is not optimized yet.");
  Pose relativeTransform;
  //transforms points from camera coords (local) to cylinder coords (global)
  relativeTransform.BecomeRelativeTransform(camera->GetPose(), rectPose);

  unsigned int width, height;
  camera->GetImageSize(width, height);

  maxX = DBL_MIN;
  minX = DBL_MAX;
  maxPhi = -M_PI;
  minPhi = M_PI;

  //run over max image circle and calculate intersections with unit cylinder
  Vector3<double> testedPoint, p;
  HomgPoint3D rectPoint;
  //  double phi;
  int radius = (int)rint(pps->GetRadius());
  double ppx, ppy;
  pps->GetPrincipal(ppx, ppy);
  int pp[2] = {(int)rint(ppx), (int)rint(ppy)};
  BresenhamCircle circle(pp, radius);
  int ciclePoint[] = {0,0};
  while(circle.GetNext(ciclePoint)) {
    camera->UnProjectLocal(HomgPoint2D(ciclePoint[0], ciclePoint[1], 1.0), p, testedPoint);
    rectPoint = relativeTransform.LocalToGlobal(HomgPoint3D(testedPoint));
    SetMinMaxFrom_(rectPoint, minX, maxX, minPhi, maxPhi);
  }


  double clippingMaxX = tan(maxHalfAngleInX_);
  if(minX < -clippingMaxX) minX = -clippingMaxX;
  if(maxX > clippingMaxX) maxX = clippingMaxX;

  /* find out if the camera is receiving optical rays
     along direction of cylinder axis: */
  HomgPoint3D epiDirHomg = relativeTransform.GlobalToLocal(HomgPoint3D(1,0,0));
  Vector3<double> epiDir;
  epiDirHomg.GetEuclidean(epiDir);
  //   cerr<<"epiDir = "<<epiDir<<endl;

  HomgPoint3D epiDirNegHomg =
    relativeTransform.GlobalToLocal(HomgPoint3D(-1,0,0));
  Vector3<double> epiDirNeg;
  epiDirHomg.GetEuclidean(epiDirNeg);
  //   cerr<<"epiDirNeg = "<<epiDirNeg<<endl;

  HomgPoint2D temp;
  //bool epipoleIsInImage = false;
  if(camera->DoesPointProjectIntoImageLocal(epiDir, temp)) {
    maxX = clippingMaxX;//positive x is fully visible!
    //epipoleIsInImage = true;
    //     cerr<<"maximizing X...\n";
  }
  if(camera->DoesPointProjectIntoImageLocal(epiDirNeg, temp)) {
    minX = -clippingMaxX;//negative x is fully visible!
    //epipoleIsInImage = true;
    //     cerr<<"minimizing X...\n";
  }

  //   cerr<<"minX = "<<minX<<endl;
  //   cerr<<"maxX = "<<maxX<<endl;
  if(maxPhi<=minPhi)
    return -1;

  /*heuristic:
    if worst case is taken the circumference in pixels will give the count
    of distinguishable epipolar curves. If image resolution is not to bad
    those curves will be most likely be only distinguishable by the one pixel
    on the circumference.
    The loss of information by halfing the number of used
    epipolar curves would hence be acceptable.
    Note: the code below represents the worst case assumption.
  */
  double numEpiLines = 2.0*M_PI*pps->GetRadius();
  //  if(epipoleIsInImage) numEpiLines*=2.0;
  double deltaPhi = maxPhi - minPhi;
  relativePhiRes = numEpiLines/deltaPhi;

  /*heuristic:
    The worst case is that an epipolar curve lies on the circumference
    so that number of pixel on half of the circumference are the number
    of samples along an epipolar plane. This is only possible if the
    viewing angle is larger then 2*90Deg. Other wise the arclength of every
    great circle passing through the fov is the same.
 */
  //  double numSamplesOnEpiLine = M_PI*pps->GetRadius()*2.0;

  //we assume that the circumference is the circle with the most samples
  //and no great circle exists that gets more.
  double anglePerSample = 1.0/pps->GetRadius();
  //the minimale stepsize along cylinder not leaving out pixels
  //max length over anglePerSample
  double minXStepSize = tan(anglePerSample);
  relativeXRes = 1.0/minXStepSize;

  //focallengthTheta is the number of pixels per radians in worst case
  //  double focallengthTheta = numSamplesOnEpiLine/pps->GetMaxCamAngle();
  // pps->GetRadius()/pps->GetMaxCamAngle();
  //   cerr<<"focallengthTheta = "<<focallengthTheta<<endl;
  //heuristic!
  //find angle range over maxX and minX
//   double maxTheta = atan2(maxX,1);
//   //   cerr<<"maxX theta = "<<maxTheta<<endl;
//   double minTheta = atan2(minX,1);
//   //   cerr<<"minX theta = "<<minTheta<<endl;
//   double thetaRange = maxTheta-minTheta;


//   /* If sampling is applied along the radius how many pixel a spherical image
//      would have in the viewing angle of the cylinder camera: */
//   double pixelsInThetaRange = thetaRange * focallengthTheta;
  /* Since the cylinder is sampled equidistantly not equiangularily
   * a scale to the number of samples neccessary to reach the same sampling
   */

//   double deltaX = maxX - minX;
//   relativeXRes = pixelsInThetaRange/deltaX;
//   cerr<<"relativeXRes = "<<relativeXRes<<endl;

  return 0;
}

template <class InputStorageType, class OutputStorageType>
int BIAS::CylindricalRectification<InputStorageType, OutputStorageType>::
DetermineCylindricCameraBoundriesCylindric_(
                                 const ProjectionParametersBase* camera,
                                 const Pose& rectPose,
                                 double& minX, double& maxX,
                                 double& minPhi, double& maxPhi,
                                 double& relativeXRes,
                                 double& relativePhiRes) const
{
  //determine bounding volume on cylinder
  Pose relativeTransform;
  //transforms points from camera coords (local) to cylinder coords (global)
  relativeTransform.BecomeRelativeTransform(camera->GetPose(), rectPose);
  unsigned int width, height;
  camera->GetImageSize(width, height);

  maxX = DBL_MIN;
  minX = DBL_MAX;
  maxPhi = -M_PI;
  minPhi = M_PI;

  Vector3<double> testedPoint;
  HomgPoint3D rectPoint;
  Vector3<double> p;
  //    double phi;
  //check upper and lower image border
  for(unsigned int x=0; x<width; x++) {
    camera->UnProjectLocal(HomgPoint2D(x, 0, 1.0), p, testedPoint);
    if(testedPoint.NormL2() < 1e-10) BIASABORT;
    rectPoint = relativeTransform.LocalToGlobal(HomgPoint3D(testedPoint));
    SetMinMaxFrom_(rectPoint, minX, maxX, minPhi, maxPhi);

    camera->UnProjectLocal(HomgPoint2D(x, height-1, 1.0), p, testedPoint);
    if(testedPoint.NormL2() < 1e-10) BIASABORT;
    rectPoint = relativeTransform.LocalToGlobal(HomgPoint3D(testedPoint));
    SetMinMaxFrom_(rectPoint, minX, maxX, minPhi, maxPhi);
  }

  for(unsigned int y=0; y<height; y++) {
    camera->UnProjectLocal(HomgPoint2D(0, y, 1.0), p, testedPoint);
    if(testedPoint.NormL2() < 1e-10) BIASABORT;
    rectPoint = relativeTransform.LocalToGlobal(HomgPoint3D(testedPoint));
    SetMinMaxFrom_(rectPoint, minX, maxX, minPhi, maxPhi);

    camera->UnProjectLocal(HomgPoint2D(width-1, y, 1.0), p, testedPoint);
    if(testedPoint.NormL2() < 1e-10) BIASABORT;
    rectPoint = relativeTransform.LocalToGlobal(HomgPoint3D(testedPoint));
    SetMinMaxFrom_(rectPoint, minX, maxX, minPhi, maxPhi);
  }

  double clippingMaxX = tan(maxHalfAngleInX_);
// #ifdef BIAS_DEBUG
//   cerr<<"detemined image spaceing:\n";
//   cerr<<"clipMinX/clipMaxX = "<<-clippingMaxX<<"/"<<clippingMaxX<<endl;
//   cerr<<"minX/maxX = "<<minX<<"/"<<maxX<<endl;
//   cerr<<"minPhi/maxPhi = "<<minPhi<<"/"<<maxPhi<<endl;
// #endif

  if(minX < -clippingMaxX) minX = -clippingMaxX;
  if(maxX > clippingMaxX) maxX = clippingMaxX;

  /* find out if the camera is receiving optical rays
     along direction of cylinder axis: */
  HomgPoint3D epiDirHomg = relativeTransform.GlobalToLocal(HomgPoint3D(1,0,0));
  Vector3<double> epiDir;
  epiDirHomg.GetEuclidean(epiDir);
  //   cerr<<"epiDir = "<<epiDir<<endl;
  HomgPoint3D epiDirNegHomg =
    relativeTransform.GlobalToLocal(HomgPoint3D(-1,0,0));
  Vector3<double> epiDirNeg;
  epiDirHomg.GetEuclidean(epiDirNeg);
  //   cerr<<"epiDirNeg = "<<epiDirNeg<<endl;
  HomgPoint2D temp;

  if(camera->DoesPointProjectIntoImageLocal(epiDir, temp)) {
    maxX = clippingMaxX;//positive x is fully visible!
    maxPhi = M_PI;
    minPhi = -M_PI;
//     cerr<<"maximizing X...\n";
  }
  if(camera->DoesPointProjectIntoImageLocal(epiDirNeg, temp)) {
    minX = -clippingMaxX;//negative x is fully visible!
    maxPhi = M_PI;
    minPhi = -M_PI;
//     cerr<<"minimizing X...\n";
  }

  if(maxPhi<=minPhi)
    return -1;

  //very! heuristic!
  double deltPhi = maxPhi - minPhi;
  //heuristic number of epilines:
  relativePhiRes = double(height)/deltPhi;

  double deltaX = maxX - minX;
  relativeXRes = double(width)/deltaX;
  return 0;
}

template <class InputStorageType, class OutputStorageType>
int BIAS::CylindricalRectification<InputStorageType, OutputStorageType>::
DetermineCylindricCameraBoundries_(const ProjectionParametersBase* camera,
                                   const Pose& rectPose,
                                   double& minX, double& maxX,
                                   double& minPhi, double& maxPhi,
                                   double& relativeXRes,
                                   double& relativePhiRes) const
{
  const ProjectionParametersPerspective* ppp =
    dynamic_cast<const ProjectionParametersPerspective*>(camera);
  if (ppp == NULL) {
    const ProjectionParametersSpherical* pps =
      dynamic_cast<const ProjectionParametersSpherical*>(camera);
    if (pps == NULL) {
      const ProjectionParametersCylindric* ppc =
        dynamic_cast<const ProjectionParametersCylindric*>(camera);
      if (ppc == NULL) {
        BIASERR("unexpected ProjectionParameter type\n");
        return -1;
      } else {
        return DetermineCylindricCameraBoundriesCylindric_(camera,
                                                           rectPose,
                                                           minX, maxX,
                                                           minPhi, maxPhi,
                                                           relativeXRes,
                                                           relativePhiRes);
      }
    } else {
      return DetermineCylindricCameraBoundriesSpherical_(camera,
                                                         rectPose,
                                                         minX, maxX,
                                                         minPhi, maxPhi,
                                                         relativeXRes,
                                                         relativePhiRes);
    }
  } else {
    return DetermineCylindricCameraBoundriesPerspective_(camera,
                                                         rectPose,
                                                         minX, maxX,
                                                         minPhi, maxPhi,
                                                         relativeXRes,
                                                         relativePhiRes);
  }

  return 0;
}

namespace BIAS{
// TODO FIXME win32 template instantiation problem (JW)
template class BIASImage_EXPORT CylindricalRectification<unsigned char, unsigned char>;
template class BIASImage_EXPORT CylindricalRectification<float, float>;
template class BIASImage_EXPORT CylindricalRectification<unsigned char, float>;
template class BIASImage_EXPORT CylindricalRectification<float, unsigned char>;

#if defined(BUILD_IMAGE_CHAR)
// TODO FIXME win32 template instantiation problem (JW)
template class BIASImage_EXPORT CylindricalRectification<unsigned char, char>;
template class BIASImage_EXPORT CylindricalRectification<char, char>;
#endif

#if defined(BUILD_IMAGE_USHORT)
template class BIASImage_EXPORT CylindricalRectification<unsigned short,unsigned short>;
#endif

#if defined(BUILD_IMAGE_SHORT)
template class BIASImage_EXPORT CylindricalRectification<short, short>;
#endif

#if defined(BUILD_IMAGE_SHORT)&&defined(BUILD_IMAGE_USHORT)
template class BIASImage_EXPORT CylindricalRectification<unsigned short, short>;
#endif

#if defined(BUILD_IMAGE_INT)
template class BIASImage_EXPORT CylindricalRectification<int,int>;
#endif

#if defined(BUILD_IMAGE_USHORT)
template class BIASImage_EXPORT CylindricalRectification<unsigned short,float>;
#endif

#if defined(BUILD_IMAGE_USHORT) && defined(BUILD_IMAGE_INT)
template class BIASImage_EXPORT CylindricalRectification<unsigned short,int>;
#endif

#if defined(BUILD_IMAGE_DOUBLE)
template class BIASImage_EXPORT CylindricalRectification<double,double>;
#endif
}