ImageBlenderIncremental.cpp

#include <Image/ImageBlenderIncremental.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>

#ifndef WIN32
#include <sys/sysinfo.h> // determine available memeory
#endif

using namespace std;
using namespace BIAS;


// --- constructor ---
// -------------------
ImageBlenderIncremental::ImageBlenderIncremental(bool incremental) {
  cylinderHeight_ = 1.0;
  SetOuputImageSize((unsigned int)(1200));
  blendIncremental_ = incremental;
  writeVrml_ = false;
  drawImageBorders_ = false;
  horizonAlignment_ = HORIZON_ALIGNMENT_X;
  ppOutput_ = NULL;
  gaussFilter_.SetSigma(1.5);
  plane_.SetZero();
}


void ImageBlenderIncremental::
SetOutputParameters(const ProjectionParametersBase& imagegeometry,
                    RGBAuc bgcolor) {
  if (ppOutput_!=NULL) {
    delete ppOutput_;
  }
  ppOutput_ = imagegeometry.Clone();
  ppOutput_->GetImageSize(outputImageWidth_, outputImageHeight_);

#ifndef WIN32
  struct sysinfo mysysinfo;
  sysinfo(&mysysinfo);

  cout<<"output image size requested: "<<outputImageWidth_
      <<"X"<< outputImageHeight_<<endl;

  cout<<"This will require more than "
      << outputImageWidth_ * outputImageHeight_ * 64 / 1000000
      <<"MB of RAM, you have "<<mysysinfo.freeram / 1000000
      <<"MB available..."<<endl;

  if (outputImageWidth_ * outputImageHeight_ * 70 > mysysinfo.freeram) {
    // later, we should tile the output image in this situation
    BIASERR(endl<<endl<<endl<<"MEMORY SHORTAGE !!! ABORTING..."<<endl<<endl<<endl);
    BIASABORT;
  }
#endif

  BIASASSERT(outputImageHeight_>0);
  BIASASSERT(outputImageWidth_>0);
  // set sink cam and depth map pointer (or null to disable 3D)
  Projection P(*ppOutput_);
  pm_.SetSinkCam(P, Depth_.IsEmpty()?NULL:&Depth_);
  if (plane_.NormL2()>0.1) {
      cout<<"mapping across plane "<<plane_<<endl;
      pm_.SetPlane(plane_);
  }
  if (!latestMosaicHigh_.IsEmpty()) latestMosaicHigh_.Release();
  if (!latestMosaicLow_.IsEmpty()) latestMosaicLow_.Release();
  latestMosaicHigh_.Init(outputImageWidth_, outputImageHeight_, 4u, ImageBase::ST_float);
  latestMosaicHigh_.SetColorModel(ImageBase::CM_RGBA);
  latestMosaicLow_.Init(outputImageWidth_, outputImageHeight_, 4u, ImageBase::ST_float);
  latestMosaicLow_.SetColorModel(ImageBase::CM_RGBA);

  // normally, one wants to use zero weight here for bg color
  bgcolor[3] = 0;
  FillBGImage_(latestMosaicHigh_,bgcolor);
  FillBGImage_(latestMosaicLow_,bgcolor);




}

void ImageBlenderIncremental::
FillBGImage_(Image<float>& BGImage, const RGBAuc& bgcolor) {
  const unsigned int w = BGImage.GetWidth();
  const unsigned int h = BGImage.GetHeight();

  const float R=bgcolor[0];
  const float G=bgcolor[1];
  const float B=bgcolor[2];
  const float A=bgcolor[3];

  // initialize with background color
  float *pImage = BGImage.GetImageData();
  for (unsigned int x = 0; x < w; x++) {
    for (unsigned int y = 0; y < h; y++) {
      *pImage++ = R;
      *pImage++ = G;
      *pImage++ = B;
      *pImage++ = A;
    }
  }
}




void ImageBlenderIncremental::
SetUserWeightImage(const Image<float>& weights) {
  weightimage_ = weights;
}



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

  if (!camera.IsProjValid()) {
    BIASWARN("Projection 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());

  if (!blendIncremental_) {
    // batch mode. wait until all data is there
    // no work now, just store in vector ...
    BIAS::UUID uuid = camera.GetUID();
    BIASASSERT(uuid.IsValid());

    // save image to database
    inputImages_[uuid] = temp2;
    imageIDs_.push_back(uuid);
  } else {
    // ok, the real work has to be done now
    pm_.SetSourceCam(camera.GetProj());
    if (plane_.NormL2()>0.1) {
        //cout<<"mapping across plane "<<plane_<<endl;
        pm_.SetPlane(plane_);
    } else {
      //cout<<"not using a plane "<<plane_<<endl;
    }

    // map orig images to cyl
    Camera<float> camCyl;
    camCyl.Init(outputImageWidth_, outputImageHeight_, 4);

    FillBGImage_(camCyl, RGBAuc(0,0,0,0));
    //cout<<"output image size used in PM: "<<outputImageWidth_
      //  <<"X"<< outputImageHeight_<<endl;

    pm_.Map(temp2, camCyl);// mappingQuality, false, superSampling);

//#define BIAS_DEBUG_IBI
#ifdef BIAS_DEBUG_IBI
    ImageIO::Save("IBI_lastin.mip", temp2);
    ImageIO::Save("IBI_lastout.mip", camCyl);
#endif

    Camera<float> lowPassCam, highPassCam;
    GetLowPassAndHighPassImage_(camCyl, lowPassCam, highPassCam);
#ifdef BIAS_DEBUG_IBI
    ImageIO::Save("IBI_lastlowpass.mip", lowPassCam);
    ImageIO::Save("IBI_lasthighpass.mip", highPassCam);
#endif
    LowPassFusiondAndHighPassMax_(latestMosaicLow_, latestMosaicHigh_,
                                  lowPassCam, highPassCam);
#ifdef BIAS_DEBUG_IBI
    ImageIO::Save("IBI_lasthighpass_fused.mip", latestMosaicHigh_);
    ImageIO::Save("IBI_lastlowpass_fused.mip", latestMosaicLow_);
#endif

  }

  return 0;
}

void ImageBlenderIncremental::
GetMosaicRGBA(Image<unsigned char>& mosaic){
  // check if image already has correct dimensions
  if (mosaic.GetWidth()!=latestMosaicHigh_.GetWidth() ||
      mosaic.GetHeight()!=latestMosaicHigh_.GetHeight() ||
      mosaic.GetChannelCount()!=4)
    mosaic.Release();

  // clear output image
  if (mosaic.IsEmpty()) {
    mosaic.Init(latestMosaicHigh_.GetWidth(), latestMosaicHigh_.GetHeight(), 4);
  }

  // fuse low pass and high pass parts so far
  AddLowPassAndHighPassImage_(latestMosaicLow_, latestMosaicHigh_,
                              mosaic);
}

void ImageBlenderIncremental::
GetMosaicRGB(Image<unsigned char>& mosaic){
  // check if image already has correct dimensions
  if (mosaic.GetWidth()!=latestMosaicHigh_.GetWidth() ||
      mosaic.GetHeight()!=latestMosaicHigh_.GetHeight() ||
      mosaic.GetChannelCount()!=3)
    mosaic.Release();

  // clear output image
  if (mosaic.IsEmpty()) {
    mosaic.Init(latestMosaicHigh_.GetWidth(), latestMosaicHigh_.GetHeight(), 3);
  }

  // fuse low pass and high pass parts so far
  AddLowPassAndHighPassImage_(latestMosaicLow_, latestMosaicHigh_,
                              mosaic);
}

void ImageBlenderIncremental::
GetMosaicRGBA(Image<float>& mosaic){
  // check if image already has correct dimensions
  if (mosaic.GetWidth()!=latestMosaicHigh_.GetWidth() ||
      mosaic.GetHeight()!=latestMosaicHigh_.GetHeight() ||
      mosaic.GetChannelCount()!=4)
    mosaic.Release();

  // clear output image
  if (mosaic.IsEmpty()) {
    mosaic.Init(latestMosaicHigh_.GetWidth(), latestMosaicHigh_.GetHeight(), 4);
  }

  // fuse low pass and high pass parts so far
  AddLowPassAndHighPassImage_(latestMosaicLow_, latestMosaicHigh_,
                              mosaic);
}

void ImageBlenderIncremental::
GetMosaicRGB(Image<float>& mosaic){
  // check if image already has correct dimensions
  if (mosaic.GetWidth()!=latestMosaicHigh_.GetWidth() ||
      mosaic.GetHeight()!=latestMosaicHigh_.GetHeight() ||
      mosaic.GetChannelCount()!=3)
    mosaic.Release();

  // clear output image
  if (mosaic.IsEmpty()) {
    mosaic.Init(latestMosaicHigh_.GetWidth(), latestMosaicHigh_.GetHeight(), 3);
  }

  // fuse low pass and high pass parts so far
  AddLowPassAndHighPassImage_(latestMosaicLow_, latestMosaicHigh_,
                              mosaic);
}



void ImageBlenderIncremental::
CopyAlphaChannel_(const Image<float>& camCyl,
                      Image<float>& lowPassCam) {
  const float *pC = camCyl.GetImageData()+3;
  float *pL = lowPassCam.GetImageData()+3;
  for (unsigned int y = 0; y < camCyl.GetHeight(); y++) {
    for (unsigned int x = 0; x < camCyl.GetWidth(); x++) {
      *pL = *pC;
      pL += 4;
      pC += 4;
    }
  }
}

void ImageBlenderIncremental::
GetLowPassAndHighPassImage_(const Image<float>& camCyl,
                           Image<float>& lowPassCam,
                           Image<float>& highPassCam) {

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

  // weight channel is now blurred, too, undo:
  CopyAlphaChannel_(camCyl,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
  const float *pDC = camCyl.GetImageData();
  float *pDL = lowPassCam.GetImageData();
  float *pDH = highPassCam.GetImageData();
  for (unsigned int y = 0; y < highPassCam.GetHeight(); y++) {
    for (unsigned int x = 0; x < highPassCam.GetWidth(); x++) {
      *pDH++ = *pDC++ - *pDL++;
      *pDH++ = *pDC++ - *pDL++;
      *pDH++ = *pDC++ - *pDL++;
      *pDH++ = *pDC++; // take alpha value from weighted input image
      pDL++;
    }
  }
}


void ImageBlenderIncremental::
LowPassFusiondAndHighPassMax_(Image<float>& destLow, Image<float>& destHigh,
                              Image<float>& lowPassCam, Image<float>& highPassCam) {
  BIASASSERT(destLow.GetWidth() == destHigh.GetWidth());
  BIASASSERT(destLow.GetWidth() == lowPassCam.GetWidth());
  BIASASSERT(destLow.GetWidth() == highPassCam.GetWidth());
  BIASASSERT(destLow.GetChannelCount() == destHigh.GetChannelCount());
  BIASASSERT(destLow.GetChannelCount() == lowPassCam.GetChannelCount());
  BIASASSERT(destLow.GetChannelCount() == highPassCam.GetChannelCount());

  float w, wLow, weightsum;
  float *pDLow = destLow.GetImageData();
  float *pDHigh = destHigh.GetImageData();
  // draw low frequency with linear weight
  float *pL = lowPassCam.GetImageData();
  // add high frequency from image with highest weight
  float *pH = highPassCam.GetImageData();

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

    // produce some output every some pixels
    if ((i % 3000) == 0) {
      BIASCDOUT((D_IMAGEBLENDER_MINIMAL | D_IMAGEBLENDER_BLENDING), ".");
    }
    if (pL[3] > 0.0f) {
      w = pDLow[3];
      wLow = pL[3];
      weightsum = w + wLow;

      *pDLow = (w * *pDLow + wLow * *pL++) / weightsum;
      pDLow++;
      *pDLow = (w * *pDLow + wLow * *pL++) / weightsum;
      pDLow++;
      *pDLow = (w * *pDLow + wLow * *pL++) / weightsum;
      pDLow++;
      *pDLow++ = weightsum;
      pL++;

      // update only if pixel-alpha is better:
      if (pH[3] > pDHigh[3]) {
        *pDHigh++ = *pH++;
        *pDHigh++ = *pH++;
        *pDHigh++ = *pH++;
        *pDHigh++ = *pH++;
      } else {
        pDHigh += 4;
        pH += 4;
      }
    } else {
      pDLow += 4;
      pL += 4;
      pDHigh += 4;
      pH += 4;
    }
  }
}

void ImageBlenderIncremental::
AddLowPassAndHighPassImage_(Image<float>& destLow, Image<float>& destHigh,
                            Image<unsigned char>& destination) {
  int channelcount = destination.GetChannelCount();
  //cout<<"have "<<channelcount<<" channels"<<endl;
  unsigned char *pDest = destination.GetImageData();
  float *pDLow = destLow.GetImageData();
  float *pDHigh = destHigh.GetImageData();
  if (channelcount==4) {
    // 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;
      // }
    }
#ifdef BIAS_DEBUG_IBI
    ImageIO::Save("IBI_last4channellow.mip", destLow);
    ImageIO::Save("IBI_last4channelhigh.mip", destHigh);
    ImageIO::Save("IBI_last4channelsum.mip", destination);
#endif
  } else {
    BIASASSERT(channelcount==3);
    // store result in destination image
    for (unsigned int i = 0;
         i < destination.GetWidth() * destination.GetHeight();
         i ++) {
   //   if (pDLow[3] > 0.0) {
        // avoid overflow and underflow
        float sumR = *pDLow++ + *pDHigh++;
        float sumG = *pDLow++ + *pDHigh++;
        float sumB = *pDLow++ + *pDHigh++;
        pDLow++; pDHigh++; // alpha channel

        if (sumR > 255.0) sumR = 255.0;
        if (sumR < 0.0) sumR = 0.0;

        if (sumG > 255.0) sumG = 255.0;
        if (sumG < 0.0) sumG = 0.0;

        if (sumB > 255.0) sumB = 255.0;
        if (sumB < 0.0) sumB = 0.0;

        // store result in image
        *pDest++ = (unsigned char)(rint(sumR));
        *pDest++ = (unsigned char)(rint(sumG));
        *pDest++ = (unsigned char)(rint(sumB));

     // } else {
     //   pDest  +=3;
     //   pDLow  +=4;
     //   pDHigh +=4;
     //  }
    }
#ifdef BIAS_DEBUG_IBI
    ImageIO::Save("IBI_last3channellow.mip", destLow);
    ImageIO::Save("IBI_last3channelhigh.mip", destHigh);
    ImageIO::Save("IBI_last3channelsum.mip", destination);
#endif
  }
}

void ImageBlenderIncremental::
AddLowPassAndHighPassImage_(Image<float>& destLow, Image<float>& destHigh,
                            Image<float>& destination) {
  int channelcount = destination.GetChannelCount();
  //cout<<"have "<<channelcount<<" channels"<<endl;
  float *pDest = destination.GetImageData();
  float *pDLow = destLow.GetImageData();
  float *pDHigh = destHigh.GetImageData();
  if (channelcount==4) {
    // store result in destination image
    for (unsigned int i = 0;
         i < destination.GetWidth() * destination.GetHeight();
         i ++) {
      *pDest++ = *pDLow++ + *pDHigh++;
      *pDest++ = *pDLow++ + *pDHigh++;
      *pDest++ = *pDLow++ + *pDHigh++;
      *pDest++ = ALPHA_OPAQUE;
      pDLow++; pDHigh++;

    }
#ifdef BIAS_DEBUG_IBI
    ImageIO::Save("IBI_last4channellow.mip", destLow);
    ImageIO::Save("IBI_last4channelhigh.mip", destHigh);
    ImageIO::Save("IBI_last4channelsum.mip", destination);
#endif
  } else {
    BIASASSERT(channelcount==3);
    // store result in destination image
    for (unsigned int i = 0;
         i < destination.GetWidth() * destination.GetHeight();
         i ++) {
      *pDest++ = *pDLow++ + *pDHigh++;
      *pDest++ = *pDLow++ + *pDHigh++;
      *pDest++ = *pDLow++ + *pDHigh++;
      pDLow++; pDHigh++; // alpha channel
    }
#ifdef BIAS_DEBUG_IBI
    ImageIO::Save("IBI_last3channellow.mip", destLow);
    ImageIO::Save("IBI_last3channelhigh.mip", destHigh);
    ImageIO::Save("IBI_last3channelsum.mip", destination);
#endif
  }
}



// ----------------------------------------------
// --- blends all added images in batch mode ---
// ----------------------------------------------
bool ImageBlenderIncremental::
BlendImages(BIAS::Camera<unsigned char> &destination,
            const ProjectionParametersBase& ppc,
            const Image<float>* depthmap,
            double gaussSigma) {
  if (this->blendIncremental_) {
    BIASERR("Call GetMosaic() to get the result image.");
    BIASABORT;
  }
  if (imageIDs_.size() < 1) {
    BIASWARN("I need at least one image to blend!");
    return false;
  }

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

  Projection TargetProjection(ppc);
  pm_.SetSinkCam(TargetProjection, Depth_.IsEmpty()? NULL : &Depth_);



  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 "
       << outputImageWidth_ << "x" << outputImageHeight_ << endl;
  destination.Release();
  destination.Init(outputImageWidth_, outputImageHeight_, 4);

  destination.SetZero();


  // ------------ 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(outputImageWidth_, outputImageHeight_, 4);
    FillBGImage_(camCyl, RGBAuc(0,0,0,0));


    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);
    GetLowPassAndHighPassImage_(camCyl, lowPassCam, highPassCam);

    // 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);

    // 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");



  Camera<float> destLow;
  Camera<float> destHigh;
  destLow.Init(destination.GetWidth(), destination.GetHeight(),
               destination.GetChannelCount());
  destHigh.Init(destination.GetWidth(), destination.GetHeight(),
                destination.GetChannelCount());

  destLow.Clear(0.0);
  destHigh.Clear(0.0);


  // 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);

    LowPassFusiondAndHighPassMax_(destLow, destHigh, lowPassCam, highPassCam);

  }


  destination.Clear(0);
  AddLowPassAndHighPassImage_(destLow, destHigh, destination);


  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 ImageBlenderIncremental::
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 ImageBlenderIncremental::
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 ImageBlenderIncremental::
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 ImageBlenderIncremental::
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 = outputImageWidth_ / 2.0;
  double halfHeight = outputImageHeight_ / 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 = outputImageHeight_ / (maxY - minY);
  outputImageWidth_ = (unsigned int)((maxX - minX) * factor);

  ProjectionParametersCylindric ppc2(xMin, xMax,
                                     phiMin, phiMax,
                                     outputImageWidth_,
                                     outputImageHeight_);

  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 ImageBlenderIncremental::
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 (weightType==WEIGHT_TYPE_USERIMAGE) {
      // in case the user specified a weight image, copy the weights from this
      // image into the alpha channel of the source
      // the weight image has to be of exactly the same size as the curent
      // input image
      float *pWeight = weightimage_.GetImageData();
      float *pImage = image.GetImageData()+3;
    BIASASSERT(image.GetWidth()==weightimage_.GetWidth());
    BIASASSERT(image.GetHeight()==weightimage_.GetHeight());
    for (unsigned int x = 0; x < w; x++) {
      for (unsigned int y = 0; y < h; y++) {
        *pImage = *pWeight++;
        pImage += 4; // jump to next alpha
      }
    }
  } else {

  // 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 ImageBlenderIncremental::
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;
}