ImageBlender.cpp

#include <Image/ImageBlender.hh>
#include <math.h>
#include <MathAlgo/SVD.hh>
#include <Base/Common/FileHandling.hh>
#include <Base/Image/ImageIO.hh>
#include <Utils/ThreeDOut.hh>
#include <Base/Image/ImageConvert.hh>

using namespace std;
using namespace BIAS;


// --- constructor ---
// -------------------
ImageBlender::ImageBlender() {
  cylinderHeight_ = 1.0;
  SetOuputImageSize((unsigned int)(1200));

  writeVrml_ = false;
  drawImageBorders_ = false;
  horizonAlignment_ = HORIZON_ALIGNMENT_X;
}


// -------------------------------------
// --- adds a camera to the database ---
// -------------------------------------
int ImageBlender::AddCamera(const BIAS::Camera<unsigned char>& camera,
                            unsigned int weightType) {

  if (!camera.IsProjValid()) {
    BIASWARN("Projection is invalid! I'll ignore this image...");
    return -1;
  }

  BIAS::UUID uuid = camera.GetUID();
  if (!uuid.IsValid()) {
    BIASWARN("UUID is invalid! I'll ignore this image...");
    return -1;
  }


  // convert to rgb
  Camera<unsigned char> temp(camera);
  ImageConvert::ToRGB(camera, temp);

  // compute alpha channel weight
  Camera<float> temp2;
  temp2.Init(temp.GetWidth(), temp.GetHeight(), 3);

  // TODO: use ImageConvert::ConvertST()
  for (unsigned int y = 0; y < temp.GetHeight(); y++) {
    for (unsigned int x = 0; x < temp.GetWidth(); x++) {
      unsigned char *pD = &temp.GetImageDataArray()[y][3 * x];
      float *pD2 = &temp2.GetImageDataArray()[y][3 * x];

      pD2[0] = pD[0];
      pD2[1] = pD[1];
      pD2[2] = pD[2];
    }
  }

  ComputeAlphaChannelWeight(temp2, weightType);

  // copy projection data
  temp2.SetProj(camera.GetProj());

  // save image to database
  inputImages_[uuid] = temp2;
  imageIDs_.push_back(uuid);

  return 0;
}


// -------------------------------
// --- blends all added images ---
// -------------------------------
bool ImageBlender::BlendImages(BIAS::Camera<unsigned char> &destination,
                               const ProjectionParametersBase& ppc,
                               const Image<float>* depthmap,
                               double gaussSigma) {

  if (imageIDs_.size() < 1) {
    BIASWARN("I need at least one image to blend!");
    return false;
  }

  Camera<float> lowPassCam;
  Camera<float> highPassCam;

  Projection TargetProjection(ppc);

  ProjectionMapping<float, float> pm;
  // set sink cam and depth map pointer (or null to disable 3D)
  pm.SetSinkCam(TargetProjection, depthmap);
  InterpolationMethod mappingQuality = MapTrilinear;
  double superSampling = 1.0; // set to 2.0 for double width supersampling

  // prepare visualization cyl cam in 3D
  ThreeDOut scene3D;

  // init destination image:
  // - set size
  // - convert to RGBA
  // - set alpha values to transparent
  cout << "size of final image is "
       << cylindricImageWidth_ << "x" << cylindricImageHeight_ << endl;
  destination.Release();
  destination.Init(cylindricImageWidth_, cylindricImageHeight_, 4);

  for (unsigned int y = 0; y < destination.GetHeight(); y++) {
    for (unsigned int x = 0; x < destination.GetWidth(); x++) {
      unsigned char *pDA = &destination.GetImageDataArray()[y][4 * x];
      pDA[0] = 0;
      pDA[1] = 0;
      pDA[2] = 0;
      pDA[3] = ALPHA_TRANSPARENT;
    }
  }


  // ------------ FILTERING

  BIASCDOUT((D_IMAGEBLENDER_MINIMAL | D_IMAGEBLENDER_FILTERING),
            "filtering images..." << endl << flush);

  gaussFilter_.SetSigma(gaussSigma);

  // low-pass and high-pass filtering for all images
  for (unsigned int i = 0; i < imageIDs_.size(); i++) {
    BIASCDOUT((D_IMAGEBLENDER_MINIMAL | D_IMAGEBLENDER_FILTERING),
            "processing image " << i << endl << flush);

    Projection proj(inputImages_[imageIDs_[i]].GetProj());
    pm.SetSourceCam(proj);

    // map orig images to cyl
    Camera<float> camCyl;
    camCyl.Init(cylindricImageWidth_, cylindricImageHeight_, 4);
    for (unsigned int y = 0; y < camCyl.GetHeight(); y++) {
      for (unsigned int x = 0; x < camCyl.GetWidth(); x++) {
        float *pD = &camCyl.GetImageDataArray()[y][4 * x];
        pD[3] = ALPHA_TRANSPARENT;
      }
    }
    pm.Map(inputImages_[imageIDs_[i]],
           camCyl,
           mappingQuality, false,
           superSampling);

    // --- perform low-pass filtering ---
    BIASCDOUT((D_IMAGEBLENDER_MINIMAL | D_IMAGEBLENDER_FILTERING),
            "low-pass filtering with sigma=" << gaussSigma << "..."
            << endl << flush);

    // low-pass filtering
    gaussFilter_.Filter(camCyl, lowPassCam);

    // weight channel is now blurred, too -> we need to fix this
    for (unsigned int y = 0; y < camCyl.GetHeight(); y++) {
      for (unsigned int x = 0; x < camCyl.GetWidth(); x++) {
        float *pC = &camCyl.GetImageDataArray()[y][4 * x];
        float *pL = &lowPassCam.GetImageDataArray()[y][4 * x];
        pL[3] = pC[3];
      }
    }

    // store low-pass filtered image to disk
    BIASCDOUT(D_IMAGEBLENDER_FILTERING,
            "writing temp low-pass image to disk... " << endl << flush);
    //ImageIO::Save(string("tmp_pano_blender_low-pass_") + FileHandling::toString(i)
    //                    + ".mip",
    //                    lowPassCam);
    ImageIO::Save(string("tmp_pano_blender_low-pass_") 
                        + FileHandling::toString(i)
                        + ".mip",
                  lowPassCam);

    if (DebugLevelIsSet(D_IMAGEBLENDER_FILTERING)) {
#ifdef BLENDER_FIX_ME
      Camera<unsigned int> temp;
      temp.Init(lowPassCam.GetWidth(), lowPassCam.GetHeight(), 4);
      ImageConvert::ConvertST(dynamic_cast<const BIAS::ImageBase>(lowPassCam),
                              dynamic_cast<const BIAS::ImageBase>(temp),
                              ImageBase::CM_RGBA);

      ostringstream s;
      s << "100_low-pass_image_" << i;
      cout << "writing " << s.str() << " to disk..." << endl;

      s << ".png";
      ImageIO::Save(s.str(), temp);
#endif
    }

    // --- perform high-pass filtering ---
    BIASCDOUT((D_IMAGEBLENDER_MINIMAL | D_IMAGEBLENDER_FILTERING),
            "high-pass filtering... "<< endl << flush);

    // high-pass filtering
    highPassCam.Release();
    highPassCam.Init(lowPassCam.GetWidth(), lowPassCam.GetHeight(), 4);

    // subtract low-pass filtered image from normal image to
    // get high frequencies
    for (unsigned int y = 0; y < highPassCam.GetHeight(); y++) {
      for (unsigned int x = 0; x < highPassCam.GetWidth(); x++) {
        float *pDC = &camCyl.GetImageDataArray()[y][4 * x];
        float *pDL = &lowPassCam.GetImageDataArray()[y][4 * x];
        float *pDH = &highPassCam.GetImageDataArray()[y][4 * x];

        pDH[0] = pDC[0] - pDL[0];
        pDH[1] = pDC[1] - pDL[1];
        pDH[2] = pDC[2] - pDL[2];
        pDH[3] = pDC[3]; // take alpha value from weighted input image
      }
    }

    // store high-pass filtered image to disk
    BIASCDOUT(D_IMAGEBLENDER_FILTERING,
            "writing temp high-pass image to disk... "<< endl << flush);
    //ImageIO::Save(string("tmp_pano_blender_high-pass_") + FileHandling::toString(i)
    //                    + ".mip",
    //                    highPassCam);
    ImageIO::Save(string("tmp_pano_blender_high-pass_") 
                        + FileHandling::toString(i)
                        + ".mip",
                  highPassCam);

    if (DebugLevelIsSet(D_IMAGEBLENDER_FILTERING)) {
#ifdef BLENDER_FIX_ME
      Camera<unsigned int> temp;
      temp.Init(highPassCam.GetWidth(), highPassCam.GetHeight(), 4);
      ImageConvert::ConvertST(highPassCam, temp, ImageBase::CM_RGBA);

      ostringstream s;
      s << "100_high-pass_image_" << i;
      cout << "writing " << s.str() << " to disk..." << endl;

      s << ".png";
      ImageIO::Save(s.str(), temp);
#endif
    }
  }

  // input images can be erased to save memory (only if drawImageBorders_ = false)
  if (!drawImageBorders_) {
    inputImages_.clear();
  }

  // ----------- BLENDING
  BIASCDOUT((D_IMAGEBLENDER_MINIMAL | D_IMAGEBLENDER_BLENDING),
            endl << "blending images");

  float r, g, b, w, rLow, gLow, bLow, wLow;

  Camera<float> destLow;
  Camera<float> destHigh;
  destLow.Init(destination.GetWidth(), destination.GetHeight(),
               destination.GetChannelCount());
  destHigh.Init(destination.GetWidth(), destination.GetHeight(),
                destination.GetChannelCount());
  destination.Clear(0);
  destLow.Clear(0.0);
  destHigh.Clear(0.0);
  unsigned char *pDest = destination.GetImageData();
  float *pDLow = destLow.GetImageData();
  float *pDHigh = destHigh.GetImageData();

  // run over all images
  for (unsigned int imageCount = 0; imageCount < imageIDs_.size(); imageCount++) {

    // load current low/high pair
    ImageIO::Load(string("tmp_pano_blender_low-pass_") + FileHandling::toString(imageCount)
                   + ".mip",
                  lowPassCam);
    ImageIO::Load(string("tmp_pano_blender_high-pass_") + FileHandling::toString(imageCount)
                   + ".mip",
                  highPassCam);

    // run over all pixels
    for (unsigned int i = 0;
         i < destination.GetWidth() * destination.GetHeight() * 4;
         i += 4) {

      // produce some output every some pixels
      if ((i % 30000) == 0) {
        BIASCDOUT((D_IMAGEBLENDER_MINIMAL | D_IMAGEBLENDER_BLENDING), ".");
      }

      // draw low frequency with linear weight
      float *pL = lowPassCam.GetImageData();

      if (pL[i + 3] > 0) {
        r = pDLow[i + 0];
        g = pDLow[i + 1];
        b = pDLow[i + 2];
        w = pDLow[i + 3];

        rLow = float(pL[i + 0]);
        gLow = float(pL[i + 1]);
        bLow = float(pL[i + 2]);
        wLow = float(pL[i + 3]);

        r = (w * r + wLow * rLow) / (w + wLow);
        g = (w * g + wLow * gLow) / (w + wLow);
        b = (w * b + wLow * bLow) / (w + wLow);
        w = w + wLow;

        pDLow[i + 0] = r;
        pDLow[i + 1] = g;
        pDLow[i + 2] = b;
        pDLow[i + 3] = w;
      }

      // add high frequency from image with highest weight
      float *pH = highPassCam.GetImageData();

      if (pH[i + 3] > pDHigh[i + 3]) {
        pDHigh[i + 0] = pH[i + 0];
        pDHigh[i + 1] = pH[i + 1];
        pDHigh[i + 2] = pH[i + 2];
        pDHigh[i + 3] = pH[i + 3];
      }
    }
  }

  // store result in destination image
  for (unsigned int i = 0;
         i < destination.GetWidth() * destination.GetHeight() * 4;
         i += 4) {

    if (pDLow[i + 3] > 0.0) {
      // avoid overflow and underflow
      float sumR = pDLow[i + 0] + pDHigh[i + 0];
      float sumG = pDLow[i + 1] + pDHigh[i + 1];
      float sumB = pDLow[i + 2] + pDHigh[i + 2];
      if (sumR > 255.0) sumR = 255.0;
      if (sumG > 255.0) sumG = 255.0;
      if (sumB > 255.0) sumB = 255.0;
      if (sumR < 0.0) sumR = 0.0;
      if (sumG < 0.0) sumG = 0.0;
      if (sumB < 0.0) sumB = 0.0;

      // store result in image
      pDest[i + 0] = (unsigned char)(sumR);
      pDest[i + 1] = (unsigned char)(sumG);
      pDest[i + 2] = (unsigned char)(sumB);
      pDest[i + 3] = ALPHA_OPAQUE;
    }
  }

  BIASCDOUT((D_IMAGEBLENDER_MINIMAL | D_IMAGEBLENDER_BLENDING),
            endl << flush);

  // draw border of images into cylinder
  if (drawImageBorders_) {
    for (unsigned int i = 0; i < imageIDs_.size(); i++) {
      cout << "drawing border of image " << i << endl;

      PixelIterator it;

      Projection projection;
      projection = inputImages_[imageIDs_[i]].GetProj();

      ProjectionParametersPerspective* ppp;
      ppp = (ProjectionParametersPerspective*)(projection.GetParameters());

      // run over all border pixels of image
      ppp->GetFirstBorderPixel(it);
      do {
        HomgPoint2D p2DHom(it.x, it.y);

        // compute ray and project it to cylinder
        Vector3<double> ray, p;
        ppp->UnProjectToRay(p2DHom, p, ray);
        HomgPoint3D p3DHom(ray); // interface conversion

        // if the cylinder geometry is correct, all borders should be in the
        // cylinder image, but you never know:
        if (!ppc.DoesPointProjectIntoImage(p3DHom, p2DHom)) continue;

        // paint border pixel in red
        destination.SetPixel(255,
                             (unsigned int)(p2DHom[0]),
                             (unsigned int)(p2DHom[1]),
                             0);
        destination.SetPixel(0,
                             (unsigned int)(p2DHom[0]),
                             (unsigned int)(p2DHom[1]),
                             1);
        destination.SetPixel(0,
                             (unsigned int)(p2DHom[0]),
                             (unsigned int)(p2DHom[1]),
                             2);

        // ensure opaque alpha channel
        destination.SetPixel(ALPHA_OPAQUE,
                             (unsigned int)(p2DHom[0]),
                             (unsigned int)(p2DHom[1]),
                             3);
      } while (ppp->GetNextBorderPixel(it));
    }
  }

  // write cyl projection to output image
  Projection proj(ppc);
  destination.SetProj(proj);
  destination.UpdateMetaData();

  // visualize cyl cam in 3D
  // TODO: dafuer sorgen, dass TriangleMesh auch RGBA nimmt
  if (writeVrml_) {
    Image<unsigned char> rgbim;
    ImageConvert::Convert(destination, rgbim, ImageBase::CM_RGB);

    TriangleMesh triangleMesh;
    triangleMesh.GenerateTexturedCamera(&ppc,
                                       rgbim,
                                       1.0,
                                       1.0,
                                       1.0);
    scene3D.AddTriangleMesh(triangleMesh);

    scene3D.VRMLOut("cylindricPanorama.wrl");
  }

  return true;
}


void ImageBlender::
ComputeCylCamGeometry(ProjectionParametersCylindric &ppc) {

  if (imageIDs_.size() < 1) {
    BIASWARN("I need at least one image to compute destination geometry!");
    return;
  }

  // init projection parameters
  Vector3<double> C(inputImages_[imageIDs_[0]].GetProj().GetC());
  RMatrix RTarget(inputImages_[imageIDs_[0]].GetProj().GetR());

  // compute orientation R:
  // get axis of each image that is assumed to lie in a plane
  // and save them in a matrix
  Matrix<double> M(2 * imageIDs_.size(), 3, 0.0);



  unsigned int testAlignment = horizonAlignment_;
  if (horizonAlignment_== HORIZON_ALIGNMENT_AUTO)
    testAlignment = 0;

  double bestResiduum = HUGE_VAL;
  int selectedmode = -1;
  Vector3<double> besth(0,1,0), h(0,1,0);
  do {
    // regard alignment of horizon
    if (testAlignment == HORIZON_ALIGNMENT_X) {
      for (unsigned int i = 0; i < imageIDs_.size(); i++) {
        RMatrix R = inputImages_[imageIDs_[i]].GetProj().GetR();
        for (unsigned int c=0;c<3;c++) M[i][c] = R[c][0];
      }
    } else if (testAlignment == HORIZON_ALIGNMENT_Y) {
      for (unsigned int i = 0; i < imageIDs_.size(); i++) {
        RMatrix R = inputImages_[imageIDs_[i]].GetProj().GetR();
        for (unsigned int c=0;c<3;c++) M[i][c] = R[c][1];
      }
    } else if (testAlignment == HORIZON_ALIGNMENT_UNKNOWN) {
      // HORIZON_ALIGNMENT_UNKNOWN
      for (unsigned int i = 0; i < imageIDs_.size(); i++) {
        RMatrix R = inputImages_[imageIDs_[i]].GetProj().GetR();
        for (unsigned int c=0;c<3;c++) M[i][c] = R[c][2];
      }
    }

    // now we can solve the LES M*h=0, where h will become
    // the y axis of the cylindrical camera
    SVD svd(M, 0.1, false);
    // the smallest eigenvalue may become quite large, so we grab
    // the solution vector "by hand"
    Matrix<double> vt;
    vt = svd.GetVT();

    // now get the first, the second and the third best null vectors
    Vector3<double> secondh(0,0,1), thirdh(0,0,1);
    for (unsigned int i = 0; i < 3; i++) {
      h[i] = vt[vt.num_rows() - 1][i];
      secondh[i] =vt[vt.num_rows() - 2][i];
      thirdh[i] =vt[0][i];
    }
    // normalize them for fair comparison
    h.Normalize();
    secondh.Normalize();
    thirdh.Normalize();
    // M*h is a vector containing the scalar products (cosines) of all camera
    // vectors with the solution. Its norm is an "average" cosine and should be
    // zero since we search a perpendicular vector cos(90)=0
    double residuum = (M*h).NormL2()/double(h.Size());
    double secondresiduum = (M*secondh).NormL2()/double(h.Size());
    double thirdresiduum = (M*thirdh).NormL2()/double(h.Size());
    // avoid pathological case and division by zero
    if (residuum<1e-10) residuum = 1e-10;
    if (secondresiduum<1e-10) secondresiduum=1e-10;
    if (thirdresiduum<1e-10) thirdresiduum=1e-10;
    //cout<<"r1="<<residuum<<" r2="<<secondresiduum<<" r3="<<thirdresiduum<<endl;
    //cout<<"S is "<<svd.GetS()<<endl;
    // we dont want to get the rotation axis, therefore if the null space is
    // two-dimensional, reject this solution! Two dimensional means here that
    // the residuum of the second solution is closer to the best than to the
    // worst solution on an order of magnitude scale
    bool nullspacedim2 = thirdresiduum / secondresiduum > 10.0;
    // (log(thirdresiduum / residuum) - 2 * secondresiduum / residuum)>0;
    //if (nullspacedim2) cout<<"2-dim. null space!"<<endl;
    //else  cout<<"1-dim. null space!"<<endl;
    if (residuum<bestResiduum && !nullspacedim2) {
      // found better axis!
      besth = h;
      bestResiduum = residuum;
      selectedmode = testAlignment;
    }

    if (horizonAlignment_== HORIZON_ALIGNMENT_AUTO) testAlignment++;
  } while (testAlignment!=horizonAlignment_);
  h = besth;
  if (horizonAlignment_ == HORIZON_ALIGNMENT_AUTO) switch (selectedmode) {
    case HORIZON_ALIGNMENT_X:cout<<"Applying horizontal alignment"<<endl;break;
    case HORIZON_ALIGNMENT_Y:cout<<"Applying vertical alignment"<<endl;break;
    default:cout<<"Applying optical axis alignment"<<endl;break;
  }

  horizonAlignment_ = selectedmode;











  // project the choosen axes of all cams onto the plane defined
  // by its normal vector h using parallel projection:
  vector< Vector3<double> > projectedVectors;

  for (unsigned int imCount = 0; imCount < imageIDs_.size(); imCount++) {
    Vector3<double> orig, proj, temp;
    RMatrix Rot = inputImages_[imageIDs_[imCount]].GetProj().GetR();
    //orig[0] = Rot[0][2];
    //orig[1] = Rot[1][2];
    //orig[2] = Rot[2][2];

    // regard alignment of horizon
    if (horizonAlignment_ == HORIZON_ALIGNMENT_X) {
      orig[0] = Rot[0][0];
      orig[1] = Rot[1][0];
      orig[2] = Rot[2][0];
    }
    else if (horizonAlignment_ == HORIZON_ALIGNMENT_Y) {
      orig[0] = Rot[0][1];
      orig[1] = Rot[1][1];
      orig[2] = Rot[2][1];
    }
    else { // HORIZON_ALIGNMENT_UNKNOWN and illegal values
      orig[0] = Rot[0][2];
      orig[1] = Rot[1][2];
      orig[2] = Rot[2][2];
    }

    h.Mult(orig.ScalarProduct(h), temp);
    orig.Sub(temp, proj);

    projectedVectors.push_back(proj);
  }

  // choose vectors a and v that span the plane
  // a is the first camera's projected optical axis
  BIASASSERT(!projectedVectors.empty());
  Vector3<double> a(projectedVectors[0]);
  a.Normalize();

  // v must be orthogonal to h and a
  Vector3<double> v;
  a.CrossProduct(h, v);
  v.Normalize();

  // now calc rotation so that pano is centered in the cyl
  // for this we need the angles between successive optical axis to find the
  // biggest "gap"

  // first transform the projected vectors into coordinate system of the plane.
  // {a, v} is a base of the plane space, so multliply each vec with a
  // 2x3-matrix that contains a and v als row vectors
  vector< Vector2<double> > vectorsInPlaneSpace(projectedVectors.size());

  Matrix<double> planeTransform(2, 3, 0);
  planeTransform[0][0] = a[0];
  planeTransform[0][1] = a[1];
  planeTransform[0][2] = a[2];
  planeTransform[1][0] = v[0];
  planeTransform[1][1] = v[1];
  planeTransform[1][2] = v[2];

  for (unsigned int i = 0; i < projectedVectors.size(); i++) {
    // interface conversion
    Vector<double> v1(projectedVectors[i]);
    Vector<double> v2;

    // transform each vec to plane's 2d coord system
    planeTransform.Mult(v1, v2);
    vectorsInPlaneSpace[i] = v2;

    // normalize vec
    double norm = vectorsInPlaneSpace[i].NormL2();
    vectorsInPlaneSpace[i] /= norm;

    cout << "vec " << i << " " << vectorsInPlaneSpace[i]
         << " norm " << vectorsInPlaneSpace[i].NormL2()<< endl;
  }

  // we choose vectorsInPlaneSpace[0] as reference.
  // note: vectorsInPlaneSpace[0] points in direction of x axis

  // calc angles between ref vec and all other vecs
  vector<double> angles(vectorsInPlaneSpace.size());

  for (unsigned int i = 0; i < vectorsInPlaneSpace.size(); i++) {
    angles[i] = CalcAngleToXAxis(vectorsInPlaneSpace[i], true);
    cout << i << ". angle relative to reference vector " << angles[i] << endl;
  }

  // determine the order of the vecs and calc angles of successive vecs
  vector<bool> visited(vectorsInPlaneSpace.size());
  vector<double> successiveAngles;
  vector<int> indices;
  double sumAngles = 0.0;
  // init vector
  for (unsigned int i = 0; i < visited.size(); i++) {
    visited[i] = false;
  }

  // mark ref axis as visited, since angle between ref vec and ref vec is 0
  visited[0] = true;

  for (unsigned int i = 1; i < vectorsInPlaneSpace.size(); i++) {
    double smallestAngle = 361;
    int index = -1;

    // look for next axis
    for (unsigned int j = 0; j < vectorsInPlaneSpace.size(); j++) {
      if (!visited[j] && (angles[j] < smallestAngle)) {
        smallestAngle = angles[j];
        index = j;
      }
    }
    BIASASSERT(index != -1);

    visited[index] = true;
    double angle;
    if (successiveAngles.empty()) {
      angle = smallestAngle;
    }
    else {
      angle = smallestAngle - successiveAngles[successiveAngles.size() - 1];
    }
    successiveAngles.push_back(angle);
    indices.push_back(index);

    sumAngles += angle;
  }

  double sumAnglesRemainder = sumAngles - int(sumAngles);
  sumAngles = int(sumAngles) % 360;
  sumAngles += sumAnglesRemainder;

  cout << "sum of all angles is " << sumAngles << endl;

  BIASASSERT(sumAngles <= 360);

  // add angle between last axis und our reference axis, i.e. close the circle
  successiveAngles.push_back(360 - sumAngles);
  indices.push_back(0);

  for (unsigned int i = 0; i < indices.size(); i++) {
    cout << "index " << indices[i]
         << " angle " << successiveAngles[i] << endl;
  }

  // search for vecs with greates "gap" between them
  double biggestAngle = 0.0;
  int indexL = -1, indexR = -1;

  for (unsigned int i = 0; i < indices.size(); i++) {
    if (successiveAngles[i] > biggestAngle) {
      biggestAngle = successiveAngles[i];
      if (i == 0) {
        indexL = indices[indices.size() - 1];
      }
      else {
        indexL = indices[i - 1];
      }
      indexR = indices[i];
    }
  }

  cout << "biggest gap is " << biggestAngle << " degrees between axis "
       << indexL << " and " << indexR << endl;

  // calc vector that lies in the middle of this two vecs
  Vector2<double> midVec;
  if (biggestAngle != 180.0) {
    midVec = (vectorsInPlaneSpace[indexL] + vectorsInPlaneSpace[indexR]) / 2.0;
  }
  else { // biggest gap is exactly 180 degrees
    // take vectorsInPlaneSpace[indexL] and rotate 90 degrees to get midVec
    Matrix2x2<double> rotTmp;

    rotTmp[0][0] =  0;
    rotTmp[0][1] = 1;
    rotTmp[1][0] = -1;
    rotTmp[1][1] =  0;

    cout << "found 180 deg gap - using rot " << rotTmp << endl;

    cout << "with vec " << vectorsInPlaneSpace[indexL] << endl;

    rotTmp.Mult(vectorsInPlaneSpace[indexL], midVec);

    cout << "mid vec is " << midVec << endl;
  }
  // normalize vec
  midVec /= midVec.NormL2();
  cout << "optical axis after norm " << midVec << endl;

  // transform midVec back to 3d space and use it as cylinder's optical axis
  Matrix<double> planeTransformInv(3, 2, 0);
  SVD svdInvert;
  planeTransformInv = svdInvert.Invert(planeTransform);

  // interface conversion
  Vector<double> v1(midVec);
  Vector<double> v2;

  planeTransformInv.Mult(v1, v2);
  a = v2;

  cout << "final optical axis is " << a << endl;

  // v must be orthogonal to h and a
  a.CrossProduct(h, v);
  v.Normalize();

  // modify RTarget
  RTarget[0][0] = h[0];
  RTarget[1][0] = h[1];
  RTarget[2][0] = h[2];
  RTarget[0][1] = v[0];
  RTarget[1][1] = v[1];
  RTarget[2][1] = v[2];
  RTarget[0][2] = a[0];
  RTarget[1][2] = a[1];
  RTarget[2][2] = a[2];

  ppc.SetC(C);
  ppc.SetR(RTarget);

  CheckFov(ppc);

  // CheckFov alters ppc so we need to set this again
  // TODO: dies hier mit nach CheckFov()
  ppc.SetC(C);
  ppc.SetR(RTarget);
}


inline double ImageBlender::CalcAngleToXAxis(const Vector2<double> &v,
                                             bool wantDegrees) {
  Vector2<double> xAxis(1.0, 0.0);
  double ret = 0.0;
  double temp = 0.0;

  // the angle of two vecs x, y is computed as follows:
  // angle(x, y) = arccos ( (x * y) / (|x| * |y|) )
  temp = xAxis.ScalarProduct(v);
  temp /= xAxis.NormL2() * v.NormL2();
  ret = acos(temp);

  if (wantDegrees) { // want degrees, not rad
    ret = ret * 180 / M_PI;
  }
  // if y component < 0 then this vec is "to the left", i.e. ccw
  // since this equation only yields values in [0, ... , 180[, we need to
  // consider the algebraic sign of the y component of the other vec,
  // which determines the relative direction.
  if (v[1] < 0) {
    if (wantDegrees) {
      ret = 360.0 - ret;
    }
    else {
      ret = 2 * M_PI - ret;
    }
  }

  return ret;
}


inline double ImageBlender::CalcAngleToYAxis(const Vector2<double> &v,
                                             bool wantDegrees) {
  Vector2<double> yAxis(0.0, 1.0);
  double ret = 0.0;
  double temp = 0.0;

  // the angle of two vecs x, y is computed as follows:
  // angle(x, y) = arccos ( (x * y) / (|x| * |y|) )
  temp = yAxis.ScalarProduct(v);
  temp /= yAxis.NormL2() * v.NormL2();
  ret = acos(temp);

  if (wantDegrees) { // want degrees, not rad
    ret = ret * 180 / M_PI;
  }
  // if x-component <= 0 then this vec is "to the left", i.e. ccw
  // since this equation only yields values in [0, ... , 180[, we need to
  // consider the algebraic sign of the x component of the other vec,
  // which determines the relative direction.
  if (v[0] <= 0) {
    if (wantDegrees) {
      ret = 360.0 - ret;
    }
    else {
      ret = 2 * M_PI - ret;
    }
  }

  return ret;
}


inline void ImageBlender::CheckFov(ProjectionParametersCylindric &ppc) {

  cout << "checking fov of cameras..." << endl;

  double minX=1e100, minY=1e100,
    maxX=-1e100, maxY=-1e100; // the borders of panorama in cylinder coords
  Vector3<double> p;

  for (unsigned int i = 0; i < imageIDs_.size(); i++) {
    double minZ = 1.0;
    PixelIterator it;

    Projection proj;
    proj = inputImages_[imageIDs_[i]].GetProj();

    ProjectionParametersPerspective* ppp;

    ppp = (ProjectionParametersPerspective*)(proj.GetParameters());

    // run over all border pixels of image
    ppp->GetFirstBorderPixel(it);
    do {

      HomgPoint2D x(it.x, it.y);

      // --- check FOV of input cam

      // compute ray in camera coordinate system
      Vector3<double> ray;
      ppp->UnProjectLocal(x, p, ray);

#ifdef BIAS_DEBUG
      double norm = 0;
      norm = ray.NormL2();
      if (!Equal(norm, 1.0)) {
        BIASERR("ray not unit length: " << ray);
        BIASABORT;
      }
#endif
      if (ray[2] < minZ) {
        minZ = ray[2];
      }

      // --- check FOV of cyl cam

      // compute ray and project it to cylinder
      ppp->UnProjectToRay(x, p, ray);
      HomgPoint3D p3DHom(ray); // interface conversion

      // if the cylinder geometry is correct, all borders should be in the
      // cylinder image, but you never know:
      if (!ppc.DoesPointProjectIntoImage(p3DHom, x)) continue;

      if (x[0] < minX) minX = x[0];
      if (x[1] < minY) minY = x[1];
      if (x[0] > maxX) maxX = x[0];
      if (x[1] > maxY) maxY = x[1];

    } while (ppp->GetNextBorderPixel(it));

    //double angle = asin(minZ) * 360.0 / M_PI;
#ifdef BIAS_DEBUG
    double angle = 0.0;
    angle = asin(minZ);

    BIASCDOUT(D_IMAGEBLENDER_GEOMETRY,
              "Max field of view of cam "
              << i << " is " << angle << endl << flush);
#endif
    // add some safety margin
    minZ *= 0.8;

    // prevent strange cases
    if (minZ < 0.01) {
      minZ = 0.01;
    }

    // clip everything with smaller z as behind camera
    ppp->SetMinZLocal(minZ);
    inputImages_[imageIDs_[i]].SetProj(Projection(*ppp));
  }

  // calc FOV of cyl
  double halfWidth = cylindricImageWidth_ / 2.0;
  double halfHeight = cylindricImageHeight_ / 2.0;

  double xMin = -1.0 * ((halfWidth - minX) / halfWidth);
  double xMax = (maxX - halfWidth) / halfWidth;

  double phiMin = -M_PI * ((halfHeight - minY) / halfHeight);
  double phiMax = M_PI * ((maxY - halfHeight) / halfHeight);

  double factor = cylindricImageHeight_ / (maxY - minY);
  cylindricImageWidth_ = (unsigned int)((maxX - minX) * factor);

  ProjectionParametersCylindric ppc2(xMin, xMax,
                                     phiMin, phiMax,
                                     cylindricImageWidth_,
                                     cylindricImageHeight_);

  ppc = ppc2;
}


inline double CalculateWeightCircular(int dx,
                                             int dy,
                                             double distCenter) {
  double ret = (1.0 - sqrt(double(dx * dx + dy * dy)) / distCenter) * 255.0;

  if (ret < 0.0)
    return 0.0;
  else
    return ret;
}


inline double CalculateWeightRectangular(unsigned int w,
                                                unsigned int h,
                                                unsigned int x,
                                                unsigned int y,
                                                unsigned int mx,
                                                unsigned int my) {
  // if pixel is near border assign weight = 0
  //if (x == 0 || y == 0 || x == (w - 1) || y == (h - 1)) return 0;
  //if (x == 1 || y == 1 || x == (w - 2) || y == (h - 2)) return 0;

  double m, mp;
  bool aboveAC, aboveBD;
  double alphaValue = 0.0;

  // determine in wich of the four triangles 1, 2, 3, 4 the pixel lies:
  // A-----------D
  // | \   1   / |
  // |   \   /   |
  // | 4   X   2 |
  // |   /   \   |
  // | /   3   \ |
  // B-----------C

  m = double(h) / double(w);

  // determine if the pixel is above or below the line AC
  mp = double(y) / double(x);
  if (mp > m) {
    aboveAC = false;
  }
  else {
    aboveAC = true;
  }

  // determine if the pixel is above or below the line BD
  mp = double(y) / double(mx + (mx - x));
  if (mp > m) {
    aboveBD = false;
  }
  else {
    aboveBD = true;
  }

  // compute weight
  if (aboveAC && aboveBD) { // triangle 1
    alphaValue = double(y) / double(my);
  }
  else if (aboveAC && !aboveBD) { // triangle 2
    alphaValue = double(w - x) / double(mx);
  }
  else if (!aboveAC && !aboveBD) { // triangle 3
    alphaValue = double(h - y) / double(my);
  }
  else if (!aboveAC && aboveBD) { // triangle 4
    alphaValue = double(x) / double(mx);
  }

  return alphaValue * 255.0;
}


// -----------------------------------------
// --- computes the weight of each pixel ---
// -----------------------------------------
void ImageBlender::ComputeAlphaChannelWeight(BIAS::Image<float> &image,
                                             unsigned int weightType) {
  // TODO: calculate weight only for one quarter of the image, rest is symmetric

  ConvertImageToRGBA(image);

  const unsigned int w = image.GetWidth();
  const unsigned int h = image.GetHeight();

  const unsigned int mx = w / 2;
  const unsigned int my = h / 2;

  // calc length of the line from (0, 0) to (mx, my)
  const double distCenter = sqrt(double(mx * mx + my * my));


  // if WEIGHT_TYPE_MIX is used, mixRatio defines the ratio between
  // rectangular and circular weight
  // e.g. if mixRation = 0.3 then use 30 % rect weight and 70 % circ weight
  // if mixRatio is less than 1.0, gain ensures that the highest weights in the
  // final image (which are at the center of the image) are 1.0
  const double mixRatio = 0.5; // must be in [0, 1]
  const double gain = 1.0 / (mixRatio * (1.0 - mixRatio));

  for (unsigned int x = 0; x < w; x++) {
    for (unsigned int y = 0; y < h; y++) {

      unsigned int dx = mx - x;
      unsigned int dy = my - y;

      double alphaValue = 0.0;

      double rectVal, circVal;

      switch (weightType) {
        case WEIGHT_TYPE_NONE:
        alphaValue = ALPHA_OPAQUE;
        break;

        case WEIGHT_TYPE_CIRCULAR:
        alphaValue = CalculateWeightCircular(dx, dy, distCenter);
        break;

        case WEIGHT_TYPE_CIRCULAR_FIT:
        if (mx < my) {
          alphaValue = CalculateWeightCircular(dx, dy, mx);
        }
        else {
          alphaValue = CalculateWeightCircular(dx, dy, my);
        }
        break;

        case WEIGHT_TYPE_RECTANGULAR:
        alphaValue = CalculateWeightRectangular(w, h, x, y, mx, my);
        break;

        case WEIGHT_TYPE_MIX:
        rectVal = CalculateWeightRectangular(w, h, x, y, mx, my) / 255.0;
        circVal = CalculateWeightCircular(dx, dy, distCenter) / 255.0;
        alphaValue = (mixRatio * rectVal)
                     * ((1.0 - mixRatio) * circVal)
                     * gain
                     * 255.0;
        break;

        case WEIGHT_TYPE_RECT_NONLINEAR:
        rectVal = CalculateWeightRectangular(w, h, x, y, mx, my) / 255.0;
        alphaValue = rectVal * rectVal * 255.0;
        break;

        default:
        BIASERR("Unknown weight type!");
      }

      // avoid clipping
      if (alphaValue < 0.0)   alphaValue =   0.0;
      if (alphaValue > 255.0) alphaValue = 255.0;

      image.GetImageDataArray()[y][4 * x + 3] = alphaValue;
    }
  }
}


// ------------------------------
// --- converts image to RGBA ---
// ------------------------------
void ImageBlender::ConvertImageToRGBA(BIAS::Image<float> &image) {

  // convert image to RGB
  if (image.GetColorModel() != ImageBase::CM_RGB) {
    Image<float> tempImage;
    ImageConvert::Convert(image, tempImage, ImageBase::CM_RGB);
    image = tempImage;
  }

  // add alphachannel
  Image<float> alphaImage(image.GetWidth(),
                          image.GetHeight(), 4);

  alphaImage.SetUID(image.GetUID());

  alphaImage.SetColorModel(ImageBase::CM_RGBA);
  for (unsigned int y = 0; y < alphaImage.GetHeight(); y++) {
    for (unsigned int x = 0; x < alphaImage.GetWidth(); x++) {
      float *pD = &image.GetImageDataArray()[y][3 * x];
      float *pDA = &alphaImage.GetImageDataArray()[y][4 * x];
      *pDA++ = *pD++;
      *pDA++ = *pD++;
      *pDA++ = *pD++;
      *pDA++ = ALPHA_OPAQUE;
    }
  }

  image = alphaImage;
}