BackwardMapping.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 "BackwardMapping.hh"
#include <Geometry/RMatrix.hh>
#include <Base/Image/ImageConvert.hh>
#include <Base/Geometry/HomgPoint2D.hh>
#include <Base/Common/BIASpragma.hh>
#include <Image/AffineMapping.hh>
#include <Base/Math/Random.hh>
#include <Base/Image/ImageIO.hh>
#include <Base/Common/FileHandling.hh>



// backward compatibility to gcc < 4.1, all other must have omp.h !
 #ifdef BIAS_HAVE_OPENMP
 #  if  !defined(_OPENMP) && defined(WIN32)
 #    error Please check your OpenMP flags and defines - seem inconsistent.
 #  endif
 #  include <omp.h>
 #endif

#define LOG2 0.30103

using namespace std;

namespace BIAS {

#define MIN_SIGMA_ANISOTROPIC 0.3
#define STEP_THRESH (1.0-MIN_SIGMA_ANISOTROPIC)

template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
GetSourceCoordinates_(const HomgPoint2D& sink,
                      HomgPoint2D& source) const
{
  BIASERR("you have to implement derived class function !");
  // for instance, if you simply want to magnify your source image by 2 do a:
  source[0] = sink[0] * 0.5;
  source[1] = sink[1] * 0.5;
  BIASABORT;
  return -1;
}

template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
GetJacobian_(const HomgPoint2D& sink, Matrix2x2<double>& Jacobian) const
{
  // numerical approximatioin of jacobian
  // compute mapped points for 4-neighborhood
  // this is actually a very good way for the jacobian approximation
  // since we use is exactly for the step size which we assumed for
  // computation

  HomgPoint2D source_left, source_right, source_up, source_down;
  HomgPoint2D locsink(sink);
  locsink[0] -= 0.5;
  if (GetSourceCoordinates(locsink, source_left)<0) return -1;
  locsink[0] += 1.0;
  if (GetSourceCoordinates(locsink, source_right)<0) return -1;
  locsink[0] -= 0.5;
  locsink[1] -= 0.5;
  if (GetSourceCoordinates(locsink, source_up)<0) return -1;
  locsink[1] += 1.0;
  if (GetSourceCoordinates(locsink, source_down)<0) return -1;

  // difference vectors are rows of jacobian
  Jacobian[0][0] = source_right[0] - source_left[0];
  Jacobian[1][0] = source_right[1] - source_left[1];
  Jacobian[0][1] = source_down[0]  - source_up[0];
  Jacobian[1][1] = source_down[1]  - source_up[1];
  return 0;
}

template <class InputStorageType, class OutputStorageType>
void BackwardMapping<InputStorageType, OutputStorageType>::
UpdatePyramidSize(const ROI& ROIsink, int srcwidth, int srcheight) {

  int ulx,uly,lrx,lry;
  ROIsink.GetCorners(ulx,uly,lrx,lry);
  lrx -= 1;
  lry -= 1;
  int midx = (ulx+lrx)/2;
  int midy = (ulx+lrx)/2;
  double width = lrx-ulx;
  double height = lry-uly;

  // compute maximum size of pyramid
  HomgPoint2D testpoint,tmp;
  double maxscale = 0.0, scale = 0.0;

  // map some points and find maximum local scaling
  // ATTENTION: all points that are mapped to (x,y,0) in the source image
  // cause a warping of infinity. It might be incorrect to assume that only
  // the corners define the scale range !
  // Therefore 20 points in the image are mapped and tested
  Random R;
  for (unsigned int c=0; c<20; c++) {
    switch (c) {
      // image corners
    case  0:testpoint = HomgPoint2D(ulx, uly); break;
    case  1:testpoint = HomgPoint2D(ulx, lry); break;
    case  2:testpoint = HomgPoint2D(lrx, lry); break;
    case  3:testpoint = HomgPoint2D(lrx, uly); break;
      // half edges
    case  4:testpoint = HomgPoint2D(midx, uly); break;
    case  5:testpoint = HomgPoint2D(midx, lry); break;
    case  6:testpoint = HomgPoint2D(ulx, midy); break;
    case  7:testpoint = HomgPoint2D(lrx, midy); break;
      // 4 at negative diagonal
    case  8:
    case  9:
    case 10:
    case 11: {
      double fraction(double(c-7)/6.0);
      testpoint = HomgPoint2D(ulx + fraction * width,
                              uly + fraction * height);
      break;
    }
      // 4 at positive diagonal
    case 12:
    case 13:
    case 14:
    case 15: {
      double fraction(double(c-11)/6.0);
      testpoint = HomgPoint2D(ulx + fraction * width,
                              uly + (1.0-fraction) * height);
      break;
    }
    default:testpoint =
              HomgPoint2D(R.GetUniformDistributed(ulx, lrx),
                          R.GetUniformDistributed(uly, lry));

    }
    if (GetSourceCoordinates(testpoint, tmp)!=0) continue;
    // compute local sampling frequency for antialiasing
    Matrix2x2<double> Jacobian;
    GetJacobian(testpoint, Jacobian);

    // this is actually dx = NormL2( Jac * (1;0) )
    //              and dy = NormL2( Jac * (0;1) )
    // It tells how far we move in src if we move by one in sink.
    // This is important for complying with the sampling theorem.
    double
      dx = sqrt(Jacobian[0][0]*Jacobian[0][0]+
                Jacobian[1][0]*Jacobian[1][0]),
      dy = sqrt(Jacobian[0][1]*Jacobian[0][1]+
                Jacobian[1][1]*Jacobian[1][1]);
    if (dx<JACOBI_THRESH) dx = JACOBI_THRESH;
    if (dy<JACOBI_THRESH) dy = JACOBI_THRESH;

    // use larger step to select pyramid level
    scale = (dx > dy)? log(dx)*ONE_OVER_LOG2 : log(dy)*ONE_OVER_LOG2;

    if (scale>maxscale) maxscale = scale;
  }

  // obtain maximal pyramid size
  double maxpossiblepyrsize = log(min(srcwidth, srcheight)/16.0) / LOG2;
  if (maxpossiblepyrsize<1.0) maxpossiblepyrsize = 1.0;

  // + 0.5 (round) + 1.0 (offset of original resolution) + 0.5 layer safety
  SetPyramidSize(min(int(maxscale+2.0), int(maxpossiblepyrsize)));

  // restore automatic mode when this function is invoked
  autoPyramidSize_ = true;

}

template <class InputStorageType, class OutputStorageType>
void BackwardMapping<InputStorageType, OutputStorageType>::
SetPyramidSize(const int newsize)
{
  int oldsize = pyramid_.size();
  if (oldsize>0 && oldsize<newsize) {
    // too small, make extra layers
    pyramid_.CreateAdditionalLayer(newsize-oldsize, 4);
  } else {
    // not initialized yet
    if (oldsize==0) pyramid_.Init(newsize);
    else {
      Image<InputStorageType> backupIm(*pyramid_[0]);
      pyramid_.Init(backupIm, newsize);
      // why do we waste time in making a smaller pyramid ????
      BIASWARN("Potential performance problem, pyramid re-initialized."
               "old size="<<oldsize<<" newsize="<<newsize);
    }
  }
  // user is responsible now
  pyramidSize_ = newsize;
  autoPyramidSize_ = false;
}


template <class InputStorageType, class OutputStorageType>
void BackwardMapping<InputStorageType, OutputStorageType>::
BuildPyramid_(const Image<InputStorageType>& source, bool forceInit,
              int numlayers)
{
  // double check that pyramidSize_ is not too big for the given image
  int thesize = int(log(min(source.GetWidth(), source.GetHeight())/4.0)/LOG2);
  if (thesize<1) thesize = 1;
  if ((int)pyramidSize_>thesize) {
    BIASWARNONCE("Backwardmapping: required/specified pyramid size is "
                 <<pyramidSize_<<", but image size "
                 <<min(source.GetWidth(), source.GetHeight())
                 <<" allows only pyramid size of " <<thesize<<". ");

    pyramidSize_ = thesize;
  }


  // if forceinit we use source to build the pyramid, otherwise
  // we check if there is already a pyramid
  if (forceInit) {
    if (!pyramid_.IsEmpty()) pyramid_.Clear();
  }

  if (pyramid_.IsEmpty()) {
   if (numlayers>0) // this is for init with one layer first...
     pyramid_.InitFromImageBase(source, numlayers);
   else
     pyramid_.InitFromImageBase(source, pyramidSize_);
  }

#ifdef BIAS_DEBUG
  //static int mycounter = 0;
  //pyramid_.WriteImages("BWMpyramid"+toString(mycounter++));
#endif
}

// #####################  dispatcher functions #########################
template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
Map(const Image<InputStorageType>& source,
    Image<OutputStorageType>& sink,
    InterpolationMethod interpolationQuality, bool newSink,
    double SuperSampling)
{
#ifdef BIAS_DEBUG
  // sanity checks
  if (!newSink){
    if (sink.IsEmpty()){
      BIASERR("sink is empty !");
      BIASABORT;
    }
    if (sink.GetChannelCount()!=source.GetChannelCount() ||
        sink.GetColorModel()!=source.GetColorModel()) {
      BIASERR("different image types !");
      BIASABORT;
    }
    if (sink.IsPlanar() && sink.GetChannelCount()>1) {
      BIASERR("not implemented for planar images !");
      BIASABORT;
    }
  }
  if (source.IsPlanar() && source.GetChannelCount()>1){
    BIASERR("not implemented for planar images !");
    BIASABORT;
  }
#endif
  width_  = source.GetWidth();
  height_ = source.GetHeight();
  if (newSink) {
    CalcCoordOffset_(sink, source);
    ChangeImgSize_(sink, source);
  } else {
    offsetX_=0; offsetY_=0;
  }

  // record interpolation mode since the simple map interface doesnt
  interpolationType_ = interpolationQuality;

#ifdef BIAS_DEBUG
  switch (interpolationQuality) {
  case MapNearestNeighbor:
  case MapBilinear:
  case MapBicubic:
  case MapTrilinear:
  case MapAnisotropic:break;
  default: BIASERR("unknown interpolation method."); BIASABORT;
  }
#endif

  Image<OutputStorageType>* pSink = &sink;
  // if supersampling is active generate temporary large image
  if (SuperSampling>1.0) {
    int newwidth = int(rint(double(sink.GetWidth()+2)*SuperSampling)),
      newheight = int(rint(double(sink.GetHeight()+2)*SuperSampling));
    pSink = new Image<OutputStorageType>(newwidth, newheight,
                                         sink.GetChannelCount());
    pSink->SetZero();
    superSampling_ = SuperSampling;
  }
  incomplete_ = false;
  int result = 0;
  if (interpolationQuality>=MapTrilinear) {
    // here we work with a pyramid
    // to avoid computational overhead, we save the pyramid and thus the
    // source image. We dont need to calculate it every time

    // check if pyramid is large enough for current mapper
    //  if (autoPyramidSize_) UpdatePyramidSize(*pSink);
    // init pyramid of size 1 (for auto mode), later, after ROI is determined, actually
    // required pyramid size is determined
    if (autoPyramidSize_) BuildPyramid_(source, true, 1);
    else BuildPyramid_(source, true);

    ROI_ = *source.GetROI();
    result = MapTri_(*pSink);
  } else {
    ROI_ = *source.GetROI();
    result = MapBi_(source, *pSink, interpolationQuality);
  }

  // now map from large image to small image
  if (SuperSampling>1.0) {
    //cout<<"writing supersampled image to ./supersampled.mip ..."<<endl;
    //ImageIO::Save("supersampled.mip", *pSink);
    AffineMapping<OutputStorageType, OutputStorageType> Resizer;
    Resizer.SetAffineTransformation(SuperSampling, 0, 0, SuperSampling,
                                    0.5*SuperSampling, 0.5*SuperSampling);
    int res2 = Resizer.Map(*pSink, sink, MapTrilinear, false, 1.0);
    if (result==0) result = res2;
    delete pSink;
    superSampling_ = 1.0;
  }

  // check if alle pixels could be filled
  if (incomplete_) {
    incomplete_ = false;
    if (result==0) result = 1;
  }
  return result;
}


template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
Map(Image<OutputStorageType>& sink) {
  if (interpolationType_ < MapTrilinear) {
    BIASERR("Call complete map function, no image caching necessary for"
      " simple mapper, dont have your input image any more !");
    BIASABORT;

  }
  incomplete_ = false;
  int result = MapTri_(sink);
  if (incomplete_) {
    incomplete_ = false;
    if (result==0) result = 1;
  }
  return result;
}

// #################### bi-Interpolation ################################

template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
MapBi_(const Image<InputStorageType>& source,
       Image<OutputStorageType>& sink,
       InterpolationMethod interpolationQuality)
{
  HomgPoint2D p_sink(0,0), p_source(0,0);
  p_sink[2] = 1;

  OutputStorageType** Data = sink.GetImageDataArray();
  double tmp = 0;
  const unsigned int ccs = sink.GetChannelCount();
  // calculate interesting region in sink
  int tlx,tly, brx1,bry1;
  GetBoundingBox(source.GetWidth(), source.GetHeight(),
                 sink.GetWidth(), sink.GetHeight(),
                 tlx, tly, brx1, bry1);
  tlx-=offsetX_;  brx1-=offsetX_; tly-=offsetY_; bry1-=offsetY_;
  // clip interesting region with sink roi
  ClipBoundingBoxToROICorners_(sink, tlx, tly, brx1, bry1);

  const int brx = brx1, bry = bry1;

  if (alphaImg_.GetSizeByte()>0) {
    if ( (bry-tly)>(int)alphaImg_.GetHeight() ||
         (brx-tlx)>(int)alphaImg_.GetWidth() )
      BIASERR("alphaImg too small:"<<alphaImg_.GetWidth()<<"x"<<
              alphaImg_.GetHeight()<<"  minimum size:"<<
              (brx-tlx)<<"x"<<(bry-tly));
  }
  const bool alphaMap = (alphaImg_.GetSizeByte()>0);
  //cout<<"offsetX_:"<<offsetX_<<" offsetY_:"<<offsetY_<<endl;
  //cout<<"tl:"<<tlx<<","<<tly<<" br:"<<brx<<","<<bry<<endl;
  OutputStorageType value;
  // iterate over sink and calculate the corresponding pix in source
  // check interpolation modes in speed order
  if (interpolationQuality==MapNearestNeighbor) { // nearest neighbor
    const InputStorageType** ppImageData
      = (const InputStorageType**)source.GetImageDataArray();
    for (int y=tly; y<bry; y++) {
      for (int x=tlx; x<brx; x++) {
        // image point in homogeneous coordinates
        p_sink[0] = x + offsetX_;
        p_sink[1] = y + offsetY_;
        if (GetSourceCoordinates(p_sink, p_source)!=0) {
          incomplete_ = true;
          continue;
        }
        // skip if pixel is not in ROI
        if (!ROI_.IsInROI(int(rint(p_source[0])),
                          int(rint(p_source[1])))) {
          incomplete_ = true;
          continue;
        }
        // run over color channels
        for (unsigned int j=0; j<ccs; j++) {
          //value= Data[y][x*ccs+j];
          ClipValue_(ppImageData[(int)rint(p_source[1])]
                     [(int)rint(p_source[0])*ccs+j], value);
          if (alphaMap) {
            Data[y][x*ccs+j]=
              OutputStorageType(alphaImg_.GetImageDataArray()[y-tly][x-tlx] *
                                float(Data[y][x*ccs+j]) +
                                (1.0-alphaImg_.GetImageDataArray()[y-tly][x-tlx]) *
                                float(value)
                                );
          } else
            Data[y][x*ccs+j]=value;
        }
      }
    }
  } else if (interpolationQuality==MapBilinear) { // bilinear
    for (int y=tly; y<bry; y++) {
      for (int x=tlx; x<brx; x++) {
        // image point in homogeneous coordinates
        p_sink[0] = x + offsetX_;
        p_sink[1] = y + offsetY_;
        if (GetSourceCoordinates(p_sink, p_source)!=0) {
          incomplete_ = true;
          continue;
        }
        // skip if pixel and 2x2 neighbors are not in ROI
        //if (!ROI_.IsInROI(p_source[0], p_source[1])) {
        if (!ROI_.CheckBilinearInterpolation(p_source[0], p_source[1])) {
          incomplete_ = true;
          continue;
        }
        for (unsigned int j=0; j<ccs; j++) {
          // value= Data[y][x*ccs+j];
          ClipValue_(source.BilinearInterpolationRGBInterleaved(p_source[0],
                                                                p_source[1],j),
                     value);
          if (alphaMap) {
            Data[y][x*ccs+j]=
              OutputStorageType(alphaImg_.GetImageDataArray()[y-tly][x-tlx] *
                                float(Data[y][x*ccs+j]) +
                                (1.0-alphaImg_.GetImageDataArray()[y-tly][x-tlx]) *
                                float(value)
                                );
          } else
            Data[y][x*ccs+j]=value;
        }
      }
    }
  } else if (interpolationQuality==MapBicubic) { // bicubic
    for (int y=tly; y<bry; y++) {
      for (int x=tlx; x<brx; x++) {
        // image point in homogeneous coordinates
        p_sink[0] = x + offsetX_;
        p_sink[1] = y + offsetY_;;
        if (GetSourceCoordinates(p_sink, p_source)!=0)  {
          incomplete_ = true;
          continue;
        }
        // skip if pixel and 3x3 neighbors are not in ROI
        //if (!ROI_.IsInROI(p_source[0], p_source[1])) {
        if (!ROI_.CheckBicubicInterpolation(p_source[0], p_source[1])) {
          incomplete_ = true;
          continue;
        }
        // run over color channels
        for (unsigned int j=0; j<ccs; j++) {
          if (GetImageValue_(p_source[0], p_source[1], source, j, tmp)==0) {
            // value= Data[y][x*ccs+j];
            ClipValue_(tmp, value);
            if (alphaMap) {
              Data[y][x*ccs+j]=
                OutputStorageType(alphaImg_.GetImageDataArray()[y-tly][x-tlx] *
                                  float(Data[y][x*ccs+j]) +
                                  (1.0-alphaImg_.GetImageDataArray()[y-tly][x-tlx]) *
                                  float(value)
                                  );
            } else
              Data[y][x*ccs+j]=value;
          }
        }
      }
    }
  }

  return 0;
}

// ########################## tri-Interpolatioon #########################

template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
MapTri_(Image<OutputStorageType>& sink)
{
  // sanity check (who knows ;-) )
  if (pyramid_.IsEmpty()) {
    BIASERR("Called function on uninitialized image !");
    sink.FillImageWithConstValue( (OutputStorageType)(255) );
    return -1;
  }

  // check if we can use fast path or must go the hard way:
  if (alphaImg_.GetSizeByte()==0 && offsetX_==0 && offsetY_==0 &&
      sink.GetChannelCount()==1 && pConcatenation_==NULL &&
      interpolationType_==MapTrilinear && superSampling_==1.0) {
    // fast path / default case: grey without offset and special effects
    return MapTrilinearGreySimple_(sink);
  }

  HomgPoint2D p_sink(0,0), p_source(0,0);

  OutputStorageType** Data = sink.GetImageDataArray();
  double tmp = 0, localscale = 0;
  const unsigned int ccs = sink.GetChannelCount();

  // calculate interesting region in sink by forward mapping of corners
  int tlx,tly, brx1, bry1;
  GetBoundingBox(pyramid_[0]->GetWidth(), pyramid_[0]->GetHeight(),
                 sink.GetWidth(), sink.GetHeight(),
                 tlx,tly,brx1,bry1);
  tlx-=offsetX_;  brx1-=offsetX_; tly-=offsetY_; bry1-=offsetY_;
  // clip interesint region with sink roi
  ClipBoundingBoxToROICorners_(sink, tlx, tly, brx1, bry1);
  const  int brx = brx1, bry = bry1;

  // check if pyramid is still large enough, mapping might have changed,
  // e.g. other homography on same image !
  if (autoPyramidSize_) {
    ROI roisink;
    roisink.Resize(sink.GetWidth(), sink.GetHeight());
    roisink.SetCorners(tlx,tly,brx,bry);
    UpdatePyramidSize(roisink,
                      pyramid_[0]->GetWidth(), pyramid_[0]->GetHeight());
  }

  if (alphaImg_.GetSizeByte()>0) {
    if ( (bry-tly)>(int)alphaImg_.GetHeight() ||
         (brx-tlx)>(int)alphaImg_.GetWidth() )
      BIASERR("alphaImg too small:"<<alphaImg_.GetWidth()<<"x"<<
              alphaImg_.GetHeight()<<"  minimum size:"<<
              (brx-tlx)<<"x"<<(bry-tly));
  }
  const bool alphaMap = (alphaImg_.GetSizeByte()>0);

  Matrix2x2<double> Jacobian;
  OutputStorageType value;

  //cout<<"backwardmapping is running from  "<<tlx<<";"<<tly<<" to "<<brx<<";"<<bry<<endl;

# ifdef BIAS_HAVE_OPENMP
  /*
       OpenMP parallelization:

       - All variables which are first written, the read in each loop
       iteration            should be declared private.
       - Also all loop indices must be declared private except the index of
       the parallelized loop (in our case the outermost), which is private by
       default.
       - All variables declared in the parallel region are private by default.
       - p_sink and p_source are declared firstprivate, so that their
       initialized values are visible in the parallel region (tmp and value
       are overwritten before read).
       - tlx and tly are only read and therefore safe to share.
       - Data is only addressed through the (private) loop indices and is safe
       to share.
       - default(none) disables the default scoping behavior, so that each
       variable must be scoped manually.

       It is important not to use static variables in the parallel region,
       they are not thread-safe and cause race conditions.

       For an OpenMP tutorial visit:
       https://computing.llnl.gov/tutorials/openMP/

  */
#pragma omp parallel for if(omp_get_num_procs() > 1)  \
  shared(tlx,tly,Data) \
    private(tmp,localscale,Jacobian,value) \
    firstprivate(p_sink,p_source)
# endif

    // iterate over sink and calculate the corresponding pix in source
    for ( int y=tly; y<bry; y++) {
      for ( int x=tlx; x<brx; x++) {
        // image point in homogenous coords
        p_sink[0] = x + offsetX_;
        p_sink[1] = y + offsetY_;
        // this is so slow we can make the decision in the innerst loop:
        if (interpolationType_ == MapTrilinear) {
          if (GetSourceCoordinates(p_sink, p_source)!=0)  {
            incomplete_ = true;
            continue;
          }
          // skip if pixel not in roi
          if (!ROI_.IsInROI(p_source[0], p_source[1]))  {
            incomplete_ = true;
            continue;
          }
          // compute layer for trilinear filtering, 0 means best resolution
          ComputeLocalPyramidLayer_(p_sink, localscale);

          // compute color channels
          // this can be speeded up, all the weights are the same for all
          // channels, dont recompute them every time !
          // However for now this code does it ..
          for (unsigned int j=0; j<sink.GetChannelCount(); j++) {
            if (GetTrilinearImageValue_(p_source[0], p_source[1],
                                        localscale, j, tmp) == 0) {
              // value= Data[y][x*ccs+j];
              ClipValue_(tmp, value);
              if (alphaMap) {
                Data[y][x*ccs+j]=
                  OutputStorageType(alphaImg_.GetImageDataArray()[y-tly][x-tlx] *
                                    float(Data[y][x*ccs+j]) +
                                    (1.0-alphaImg_.GetImageDataArray()[y-tly][x-tlx]) *
                                    float(value)
                                    );
              } else {
                Data[y][x*ccs+j]=value;
              }
            } else  {
              incomplete_ = true;
            }
          }
        } else {
          // anisotropic filtering !
          if (GetSourceCoordinates(p_sink, p_source)!=0)  {
            incomplete_ = true;
            continue;
          }

          // skip if pixel not in roi
          if (!ROI_.IsInROI(p_source[0], p_source[1]))  {
            incomplete_ = true;
            continue;
          }

          // compute local sampling frequency matrix
          GetJacobian(p_sink, Jacobian);

          // compute color channels
          // this can be speeded up, all the weights are the same for all
          // channels, dont recompute them every time !
          // However for now this code does it ..
          for (unsigned int j=0; j<sink.GetChannelCount(); j++) {
            if (GetAnisotropicImageValue_(p_source[0], p_source[1],
                                          Jacobian, j, tmp) == 0) {
              //value= Data[y][x*ccs+j];
              ClipValue_(tmp, value);
              if (alphaMap) {

                // should use clipvalue here

                Data[y][x*ccs+j]=
                  OutputStorageType(alphaImg_.GetImageDataArray()[y-tly][x-tlx]*
                                    float(Data[y][x*ccs+j]) +
                                    (1.0 -
                                     alphaImg_.GetImageDataArray()[y-tly][x-tlx])
                                    * float(value)
                                    );
              } else
                Data[y][x*ccs+j]=value;
            } else  {
              incomplete_ = true;
            }
          }
        }
      }
    }
#ifdef BIAS_DEBUG
  if (!rangeChecked_) {
    rangeChecked_ = true;
    BIASWARN("You may have got artefacts in your image because overflow "
             <<"checking of your data types combination "
             <<typeid(InputStorageType).name()<<" -> "
             <<typeid(InputStorageType).name()<<" is not yet implemented."
             <<" Please do it in CheckValueRange_()."<<endl
             <<" Example: If your input image is type int and your output "
             <<" image is type uchar, unpredictable things happen when your"
             <<" input image contains/creates numbers>255 or numbers<0. You"
             <<" may want to round or clip, which is done only heuristically"
             <<" at the moment for your data type combination.");
  }
#endif
  if (aliasing_) {
    aliasing_ = false;
    BIASWARN("You probably got aliasing, because your pyramid was too small."
             "If you did not set the pyramid size manually, the specialized"
             " implementation of this->UpdatePyramidSize() is wrong and you"
             " have such a large local minification that you need a"
             " smaller pyramid top than currently automatically generated "
             " (or the whole image is mapped to only a few pixels)!"
             " Workaround for too little pyramids: You could"
             " try somthing like manual this->SetPyramidSize(log2(imagesize))"
             <<endl<<" Most likely however, your transformation"
             " (this->GetSourceCoordinates, e.g. homography)"
             " is NOT CONTINUOUS (e.g. wrap-around for homography)"
             ", which is a crucial assumption in antialiased BackwardMapping");
  }

  return 0;
}


template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
MapTrilinearGreySimple_(Image<OutputStorageType>& sink) {
  // sanity check (who knows ;-) )
  if (pyramid_.IsEmpty()) {
    BIASERR("Called function on uninitialized image !");
    sink.FillImageWithConstValue( (OutputStorageType)(255) );
    return -1;
  }
  // check if pyramid is still large enough, mapping might have changed,
  // e.g. other homography on same image !
  if (autoPyramidSize_)
    UpdatePyramidSize(*sink.GetROI(),
                      pyramid_[0]->GetWidth(), pyramid_[0]->GetHeight());
  HomgPoint2D p_sink(0,0), p_source(0,0);

  OutputStorageType** Data = sink.GetImageDataArray();
  double localscale = 0;

  // calculate interesting region in sink by forward mapping of corners
  int tlx,tly, brx1, bry1;
  GetBoundingBox(pyramid_[0]->GetWidth(), pyramid_[0]->GetHeight(),
                 sink.GetWidth(), sink.GetHeight(),
                 tlx,tly,brx1,bry1);

  // clip interesint region with sink roi
  ClipBoundingBoxToROICorners_(sink, tlx, tly, brx1, bry1);
  const  int brx = brx1, bry = bry1;

  BIASASSERT(alphaImg_.GetSizeByte()==0);
  BIASASSERT(offsetX_==0);
  BIASASSERT(offsetY_==0);
  BIASASSERT(sink.GetChannelCount()==1);
  BIASASSERT(pConcatenation_==NULL);
  BIASASSERT(superSampling_==1.0);
  Matrix2x2<double> Jacobian;

  // iterate over sink and calculate the corresponding pix in source
  for ( int y=tly; y<bry; y++) {
    OutputStorageType *pLine = Data[y];
    for ( int x=tlx; x<brx; x++) {
      // image point in homogenous coords
      p_sink[0] = x;
      p_sink[1] = y;

      if (GetSourceCoordinates_(p_sink, p_source)!=0) {
        incomplete_ = true;
        continue;
      }
      // skip if pixel not in roi
      if (!ROI_.IsInROI(p_source[0], p_source[1])) {
        incomplete_ = true;
        continue;
      }
      // compute layer for trilinear filtering, 0 means best resolution

      GetJacobian_(p_sink, Jacobian);
      double
        dx = sqrt(Jacobian[0][0]*Jacobian[0][0]+
                  Jacobian[1][0]*Jacobian[1][0]),
        dy = sqrt(Jacobian[0][1]*Jacobian[0][1]+
                  Jacobian[1][1]*Jacobian[1][1]);
      if (dx<JACOBI_THRESH) dx = JACOBI_THRESH;
      if (dy<JACOBI_THRESH) dy = JACOBI_THRESH;
      // use larger step to select pyramid level
      localscale = (dx > dy)? log(dx)*ONE_OVER_LOG2 : log(dy)*ONE_OVER_LOG2;
      // clip scale to zero since we have no better resolution image than 0
      if (localscale<0.0) localscale = 0.0;

      // compute a weight for the smaller and one for the larger pyramid level
      // then compute the colors at these levels and average them according
      // to their weights

      double value1 =  0, value2 = 0;
      double scale_error = localscale - int(localscale);
      int scale1 = int(localscale);
      if (scale1>=(int)pyramid_.size()-1) {
        // pyramid is too small, create new levels on the fly
        if (pyramid_.CreateAdditionalLayer(scale1 - pyramid_.size() + 2)!=0) {
          // we might get aliasing, pyramid is too small !!!
          scale1 = pyramid_.size()-1;
          scale_error = 0.0;
          aliasing_ = true;
        }
      } else if (scale1<0) {
        // use max resolution image, dont need to upsample in pyramid
        scale1 = 0;
        scale_error = 0.0;
      }

      // compute coordinates in smaller image
      double x1 =  p_source[0], y1 =  p_source[1];
      const double offset = pyramid_.GetPositionOffset();
      for (int s=0; s<scale1; s++) {
        // offset necessary because of offset in 2^n downsampling
        x1 = (x1 - offset)  / pyramid_.GetRescaleFactor();
        y1 = (y1 - offset)  / pyramid_.GetRescaleFactor();
      }
      // if scale 1 fits perfectly (at lowest highest pyramid layer)
      // we dont need scale2, set it to same value to avoid access violations
      const int scale2 = (scale_error==0.0)? scale1 : scale1 + 1;
      const double x2  = (scale_error==0.0)? x1     : (x1 - offset) / pyramid_.GetRescaleFactor();
      const double y2  = (scale_error==0.0)? y1     : (y1 - offset) / pyramid_.GetRescaleFactor();

      double weight1 = 1.0 - scale_error, weight2 = scale_error;
      BIAS::Image<InputStorageType> &im =  *pyramid_[scale1];


      // compute value only in image (no pyramid), first try bicubic
      if ((scale1==0) &&
          (x1 <  (double)im.GetROI()->GetCornerLowerRightX() - 2.0) &&
          (y1 <  (double)im.GetROI()->GetCornerLowerRightY() - 2.0) &&
          (x1 >= (double)im.GetROI()->GetCornerUpperLeftX()  + 1.0) &&
          (y1 >= (double)im.GetROI()->GetCornerUpperLeftY()  + 1.0) ) {
        value1 = im.BicubicInterpolation(x1, y1, 0);
      } else if ((x1 < (double)im.GetROI()->GetCornerLowerRightX()-1.0) &&
                 (y1 < (double)im.GetROI()->GetCornerLowerRightY()-1.0) &&
                 (x1 > (double)im.GetROI()->GetCornerUpperLeftX()) &&
                 (y1 > (double)im.GetROI()->GetCornerUpperLeftY()) ) {
        // we are near the border, bilinear
        value1 = im.FastBilinearInterpolationGrey(x1, y1);
      } else {
        // we are AT the border, we could do nearest neighbor here or
        // something else, but for now, skip this pixel
        weight1=0;
      }

      // compute second scale only if necessary
      if (weight2>0.0) {
        const ROI& theROI = *pyramid_[scale2]->GetROI();

        if (x2 < (double)theROI.GetCornerLowerRightX()-1.0 &&
            y2 < (double)theROI.GetCornerLowerRightY()-1.0 &&
            x2 > (double)theROI.GetCornerUpperLeftX() &&
            y2 > (double)theROI.GetCornerUpperLeftY()) {
          // we are near the border, bilinear
          value2 = pyramid_[scale2]->FastBilinearInterpolationGrey(x2, y2);
        } else {
          // we are AT the border, we could do nearest neighbor here or
          // something else, but for now, skip this pixel
          weight2=0;
        }
      }
      // check if we computed any value at all
      const double weightsum = weight1 + weight2;
      if ((weightsum)>DBL_EPSILON) {
        ClipValue_((weight1 * value1 + weight2 * value2)/weightsum, pLine[x]);
      } else {
        incomplete_ = true;
      }
    }
  }
#ifdef BIAS_DEBUG
  if (!rangeChecked_) {
    rangeChecked_ = true;
    BIASERR("You may have got artefacts in your image because overflow "
            <<"checking of your data types combination "
            <<typeid(InputStorageType).name()<<" -> "
            <<typeid(InputStorageType).name()<<" is not yet implemented."
            <<" Please do it in CheckValueRange_().");
  }
#endif

  if (aliasing_) {
    aliasing_ = false;
    BIASERR("You probably got aliasing, because your pyramid was too small."
            "If you did not set the pyramid size manually, the specialized "
            " implementation of this->UpdatePyramidSize() is wrong!");
  }
  return 0;
}

template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
GetTrilinearImageValue_(const double& xsource,
                        const double& ysource,
                        const double& scale,
                        unsigned int channel,
                        double& thereturnvalue)
{
  // compute a weight for the smaller and a weight for the larger pyramid level
  // then compute the colors at these levels and average them according
  // to their weights
  double value1 =  0, value2 = 0;
  double scale_error = scale - int(scale);
  int scale1 = int(scale);

  if (scale1>=(int)pyramid_.size()-1) {
    // pyramid is too small, create new levels on the fly
    if (pyramid_.CreateAdditionalLayer(scale1 - pyramid_.size() + 2)!=0) {
      // we might get aliasing, pyramid is too small !!!
      scale1 = pyramid_.size()-1;
      scale_error = 0.0;
      aliasing_ = true;
    }
  } else if (scale1<0) {
    // use max resolution image, dont need to upsample in pyramid
    scale1 = 0;
    scale_error = 0.0;
  }

  // compute coordinates in smaller image
  double x1 = xsource, y1 = ysource;
  const double offset = pyramid_.GetPositionOffset();
  for (int s=0; s<scale1; s++) {
    // offset necessary because of offset in 2^n downsampling
    x1 = (x1 - offset) / pyramid_.GetRescaleFactor();
    y1 = (y1 - offset) / pyramid_.GetRescaleFactor();
  }

  // if scale 1 fits perfectly (at lowest highest pyramid layer)
  // we dont need scale2, set it to the same value to avoid access violations
  const int scale2 = (scale_error==0.0)? scale1 : scale1 + 1;
  const double x2  = (scale_error==0.0)? x1     : (x1 - offset) / pyramid_.GetRescaleFactor();
  const double y2  = (scale_error==0.0)? y1     : (y1 - offset) / pyramid_.GetRescaleFactor();

  double weight1 = 1.0 - scale_error, weight2 = scale_error;

  if (GetImageValue_(x1, y1, *pyramid_[scale1], channel, value1)!=0) weight1=0;

  // compute second scale only if necessary
  if (weight2>0.0) {
    if (GetImageValue_(x2, y2, *pyramid_[scale2], channel, value2)!=0) {
      weight2=0;
    }
  }
  // check if we computed any value at all
  const double weightsum = weight1 + weight2;
  if ((weightsum)<=DBL_EPSILON) return -1;

  // this is the tri in trilinear:
  thereturnvalue = (weight1 * value1 + weight2 * value2) / weightsum;

  return 0;
}


template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
GetAnisotropicImageValue_(const double& xsource,
                          const double& ysource,
                          Matrix2x2<double>& Jacobian,
                          unsigned int channel,
                          double& thereturnvalue)
{
  // compute eigenvalues of jacobian, which determine the largest step sizes
  Vector2<double> ev1, ev2;

  // make postive definite jTj for gauss distribution
  //Jacobian = Jacobian.Transpose() * Jacobian;
  Matrix2x2<double> Smoother;
  Jacobian.GetSystemMatrix(Smoother);

  // compute eigenvalues = squared step sizes in source
  double d1 = 0.0, d2 = 0.0;
  Smoother.EigenvalueDecomposition(d1, ev1, d2, ev2);

  // remove square, require min eigenvalue 1, because we have no finer image
  // information
  // (should rarely happen for meaningful mappings)
  if (d1<STEP_THRESH) d1 = sqrt(STEP_THRESH); else d1 = sqrt(d1);
  if (d2<STEP_THRESH) d2 = sqrt(STEP_THRESH); else d2 = sqrt(d2);

  // make sure d1 is not smaller than d2
  if (d1<d2) {
    std::swap(d1, d2);
  }

  // compute pyramid levels
  const double largescale = log(d1)*ONE_OVER_LOG2;

  if (largescale<=0.1) {
    // bicubic is sufficient since we are only enlarging
    return GetImageValue_(xsource, ysource,
                          *pyramid_[0], channel, thereturnvalue);
  }


  return pyramid_.GetAnisotropicImageValue(xsource, ysource, Smoother,
                                           thereturnvalue, channel);
}


// ################## data type clipping: double -> char, ... ############

template <class InputStorageType, class OutputStorageType>
void BackwardMapping<InputStorageType, OutputStorageType>::
ClipBoundingBoxToROICorners_(const Image<OutputStorageType>& sink,
                             int& tlx,
                             int& tly,
                             int& brx,
                             int& bry) const
{
  unsigned int minx, miny, maxx, maxy;
  sink.GetROICorners(minx, miny, maxx, maxy);
  if (tlx<(int)minx) tlx = minx;
  if (tly<(int)miny) tly = miny;
  if (brx>(int)maxx) brx = maxx;
  if (bry>(int)maxy) bry = maxy;
}

template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
CalcCoordOffset_(const  Image<OutputStorageType> &sink,
                 const Image<InputStorageType> &source) {
  //sink bounding box
  int xl=0;  int xr=0; int  yl=0; int yr=0;
  GetBoundingBox(sink.GetWidth(), sink.GetHeight(),
                 source.GetWidth(), source.GetHeight(), xl, yl, xr, yr);

  // sink corner coordinates
  double xls=0; //double xrs=sink.GetWidth();
  double yls=0; //double yrs=sink.GetHeight();


  // result bounding box
  double xlr=0;  double  ylr=0;
  //double xrr=0;  double yrr=0;
  if (xl<xls) xlr=xl; else xlr=xls;
  if (yl<yls) ylr=yl; else ylr=yls;
  //if (xr>xrs) xrr=xr; else xrr=xrs;
  //if (yr>yrs) yrr=yr; else yrr=yrs;

  offsetX_= (int)ceil(xlr);
  offsetY_= (int)ceil(ylr);

  return 0;
}

template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
ChangeImgSize_ (Image<OutputStorageType> &sink,
                const Image<InputStorageType> &source) {
  BIAS::Image<OutputStorageType> temp;
  int xl=0;  int xr=0; int  yl=0; int yr=0;

  GetBoundingBox(source.GetWidth(), source.GetHeight(),
                 sink.GetWidth(), sink.GetHeight(), xl, yl, xr, yr);

  // sink corner coordinates
  double xls=0;  double xrs=sink.GetWidth();
  double  yls=0; double yrs=sink.GetHeight();

  // result bounding box
  double xlr=0;  double xrr=0; double  ylr=0; double yrr=0;
  if (xl<xls) xlr=xl; else xlr=xls;
  if (yl<yls) ylr=yl; else ylr=yls;
  if (xr>xrs) xrr=xr; else xrr=xrs;
  if (yr>yrs) yrr=yr; else yrr=yrs;

  int x_abs = (int)ceil(abs(xlr-xrr));
  int y_abs = (int)ceil(abs(ylr-yrr));


  const int numchannels = sink.GetChannelCount();
  temp.Init(x_abs,y_abs,numchannels);
  temp.SetZero();

  int w = sink.GetWidth();
  int h = sink.GetHeight();


  OutputStorageType** data = sink.GetImageDataArray();
  OutputStorageType** datac = temp.GetImageDataArray();
  int xt = (int)abs(offsetX_);
  int yt = (int)abs(offsetY_);
  cout << "translation at tnc " << xt << " " << yt << "\n";

  for (int hi =0; hi<h; hi++) {
    for (int wi=0; wi<w; wi++) {
      for (int c=0; c<numchannels; c++) {
        datac[hi+yt][(wi+xt)*numchannels+c] = data[hi][wi*numchannels+c];
        //cout << hi << " " << wi << "\n";
        // cout <<  hi + 3*c +yt << " " << wi +xt << "\n";
      }
    }
  }
  sink=temp;

  return 0;
}

//Computes Displacement Map, saving the unmapped texture-coordinates of sink
//pixels into an image as described in the header-file
template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
GetDisplacementMap(Image<float>& dismap,
                   int width,
                   int height) {
  Image<InputStorageType> src (width, height);
  return GetDisplacementMap(dismap, src);
}

//Computes Displacement Map, saving the unmapped texture-coordinates of sink
//pixels into an image as described in the header-file
template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
GetDisplacementMap(Image<float>& dismap,
                   Image<InputStorageType> &source) {

  if(dismap.GetChannelCount() != 3){
    cout << "Displacement Map must have 3 Channels but has only "
         << dismap.GetChannelCount() << endl;
    return -1;
  }

  if(!dismap.IsInterleaved()){
    cout << "Displacement map must be INTERLEAVED" << endl;
    return -1;
  }

  float* Data = dismap.GetImageData();
  //TEXTURE Coordinates of unmapped pixels
  double gl_x, gl_y;

  HomgPoint2D sink_p, src_p;
  int count = 0;

  sink_p[2] = 1;
  for(unsigned int j = 0; j < dismap.GetHeight(); j++){
      sink_p[1] = j + offsetY_;

    for(unsigned int i = 0; i < dismap.GetWidth(); i++){
      sink_p[0] = i + offsetX_;

      if(GetSourceCoordinates(sink_p, src_p) != -1 &&
         src_p[0] < source.GetWidth() &&
         src_p[1] < source.GetHeight() &&
         src_p[0] >= 0 &&
         src_p[1] >= 0){
      //  src_p.Homogenize();
        //converting to TEXTURE-Coordinates
        source.BIASToTextureCoordinates(src_p[0], src_p[1],
                                          gl_x,     gl_y);


        Data[count  ] = (float)gl_x;
        Data[count+1] = (float)gl_y;
        Data[count+2] = 1.0;
      }
      else{
        Data[count  ] = -1.0;
        Data[count+1] = -1.0;
        Data[count+2] = -1.0;
      }
      count+=3;
    }
  }

 // dismap.Flip();

  return 0;
}

template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
GetLookupTableSize(unsigned int& width, unsigned int& channels,
                   InterpolationMethod method)
{
  switch (method)
  {
  case MapNearestNeighbor:
    BIASERR("not yet implemented!");
    return -1;
    break;
  case MapBilinear:
    {
      //if lookups were not computed
      if(pLUT_Bil_ == NULL){
        BIASERR("LOOKUP-TABLE FOR BILINEAR IS NOT READY!"
          "Call PrepareLookupTableMapping");
        return -1;
      }

      //prepare image size
      width = pLUT_Bil_size+1;
      channels = BWM_LUT_Entry_Bilinear::GetSerializedSize();
    }
    break;
  case MapTrilinear:
    {
      //if lookups were not computed
      if(pLUT_Tri_ == NULL){
        BIASERR("LOOKUP-TABLE FOR BILINEAR IS NOT READY!"
          "Call PrepareLookupTableMapping");
        return -1;
      }
      width = pLUT_Tri_size+1;
      channels = BWM_LUT_Entry_Trilinear::GetSerializedSize();
    }
    break;
  default:
    BIASERR("not yet implemented!");
    return -1;
    break;
  }

  return 0;
}

template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
GetLookupTable(const string& filename, InterpolationMethod method)
{
  switch (method)
  {
  case MapNearestNeighbor:
    BIASERR("not yet implemented!");
    return -1;
    break;
  case MapBilinear:
    {
      //if lookups were not computed
      if(pLUT_Bil_ == NULL){
        BIASERR("LOOKUP-TABLE FOR BILINEAR IS NOT READY!"
          "Call PrepareLookupTableMapping");
        return -1;
      }

      // use a binary stream for lookup table
      fstream stream;
      stream.open(filename.c_str(), fstream::out |
                                    fstream::binary);

      unsigned int sizeInFloats = pLUT_Bil_size * BWM_LUT_Entry_Bilinear::GetSerializedSize() + 1;

      //allocate buffer for transfer (mem -> file)
      float* buffer = new float[sizeInFloats];

      //write LUT to buffer...
      if (GetLookupTable(buffer, method) != 0) {
        delete[] buffer;
        return -1;
      }

      //first sizeof(unsigned int) bytes in file for LUT-size in floats
      for (unsigned int i=0; i<sizeof(unsigned int); i++)
        stream.put( ((char*)(&sizeInFloats))[i] );

      // write to file and close
      for (unsigned int i = 0; i < sizeInFloats; i++)
          stream.write((const char*)&(buffer[i]), sizeof(float));
      stream.close();

      delete[] buffer;
    }
    break;
  case MapTrilinear:
    {
      //if lookups were not computed
      if(pLUT_Tri_ == NULL){
        BIASERR("LOOKUP-TABLE FOR TRILINEAR IS NOT READY!"
          "Call PrepareLookupTableMapping");
        return -1;
      }

      // use a binary stream for lookup table
      fstream stream;
      stream.open(filename.c_str(), fstream::out |
                                    fstream::binary);

      unsigned int sizeInFloats = pLUT_Tri_size * BWM_LUT_Entry_Trilinear::GetSerializedSize() + 1;

      //allocate buffer for transfer (mem -> file)
      float* buffer = new float[sizeInFloats];

      //write LUT to buffer...
      if (GetLookupTable(buffer, method) != 0) {
        delete[] buffer;
        return -1;
      }

      //first sizeof(unsigned int) bytes in file for LUT-size in floats
      for (unsigned int i=0; i<sizeof(unsigned int); i++)
        stream.put( ((char*)(&sizeInFloats))[i] );

      //...write LUT to file and close
      for (unsigned int i = 0; i < sizeInFloats; i++)
        stream.write((const char*)&(buffer[i]), sizeof(float));
      stream.close();
      delete[] buffer;
    }
    break;
  default:
    BIASERR("not yet implemented!");
    return -1;
    break;
  }

  return 0;
}


template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
GetLookupTable(float* lutImageRow, InterpolationMethod method)
{
  float* pixel = lutImageRow;
  switch (method)
  {
  case MapNearestNeighbor:
    BIASERR("not yet implemented!");
    return -1;
    break;
  case MapBilinear:
    {
      //if lookups were not computed
      if(pLUT_Bil_ == NULL){
        BIASERR("LOOKUP-TABLE FOR BILINEAR IS NOT READY!"
          "Call PrepareLookupTableMapping");
        return -1;
      }

      *pixel = static_cast<float>(method);
      pixel++;

      for(int i = 0; i < pLUT_Bil_size; i++){
        pixel = pLUT_Bil_[i].SerializedOut(pixel);
      }
    }
    break;
  case MapTrilinear:
    {
      //if lookups were not computed
      if(pLUT_Tri_ == NULL){
        BIASERR("LOOKUP-TABLE FOR TRILINEAR IS NOT READY!"
          "Call PrepareLookupTableMapping");
        return -1;
      }
      *pixel = static_cast<float>(method);
      pixel++;

      for(int i = 0; i < pLUT_Tri_size; i++){
        pixel = pLUT_Tri_[i].SerializedOut(pixel);
      }
    }
    break;
  default:
    BIASERR("not yet implemented!");
    return -1;
    break;
  }

  return 0;
}

template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
SetLookupTable(const string& filename, InterpolationMethod& method)
{
  // try to open file
  fstream stream;
  stream.open(filename.c_str(), fstream::in |
                                fstream::binary);
  if (!stream.good())
  {
    BIASERR("Can not read file!");
    return -1;
  }

  //read size of LUT
  unsigned int sizeInFloats;
  for (unsigned int i=0; i<(sizeof(unsigned int)/sizeof(char)); i++) {
    stream.get( ((char*)(&sizeInFloats))[i] );
    if (!stream) {
      BIASERR("Error reading header from file: " << filename);
      return -1;
    }
  }

  //create buffer for LUT
  char* buffer = new char[sizeInFloats*(sizeof(float)/sizeof(char))];

  float *ptrToBuffer = (float*)buffer;

  //fill buffer from file
  unsigned int count = 0;
  while(stream.good()) {
    stream.get(*buffer++);
    count++;
  }
  //check for correct size of filled buffer
  if (count-1 != sizeInFloats*(sizeof(float)/sizeof(char))) {
    BIASERR("Error reading LUT, wrong file size!");
    delete[] buffer;
    return -1;
  }

  //finally create LUT from filled buffer
  int res = SetLookupTable(ptrToBuffer, sizeInFloats, method);
  stream.close();

  return res;
}

template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
SetLookupTable(const float* lutImageRow, unsigned int width,
               InterpolationMethod& method)
{
  const float* pixel = lutImageRow;
  int methodTmp = *pixel;
  method = (InterpolationMethod)methodTmp;
  //const InterpolationMethod method = ;
  pixel++;

  switch (method)
  {
  case MapNearestNeighbor:
    BIASERR("not yet implemented!");
    return -1;
    break;
  case MapBilinear:
    {
      delete[] pLUT_Bil_;
      pLUT_Bil_size = (width-1) / BWM_LUT_Entry_Bilinear::GetSerializedSize(); //subtract one for interpolation identifier
      pLUT_Bil_ = new BWM_LUT_Entry_Bilinear[pLUT_Bil_size];

      for(int i = 0; i < pLUT_Bil_size; i++){
        pixel = pLUT_Bil_[i].SerializedIn(pixel);
      }
    }
    break;
  case MapTrilinear:
    delete[] pLUT_Tri_;
    pLUT_Tri_size = (width-1) / BWM_LUT_Entry_Trilinear::GetSerializedSize(); //subtract one for interpolation identifier
    pLUT_Tri_ = new BWM_LUT_Entry_Trilinear[pLUT_Tri_size];

    for(int i = 0; i < pLUT_Tri_size; i++){
      pixel = pLUT_Tri_[i].SerializedIn(pixel);
    }
    break;
  default:
    BIASERR("not yet implemented!");
    return -1;
    break;
  }

  return 0;
}


template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
PrepareLookupTableMapping(const Image<InputStorageType>& src,
                          Image<OutputStorageType>& sink,
                          InterpolationMethod method, bool newSink) {

  // generate multichannel channel image newsink with x,y and sink.size
  // if trilinear: pyramid (layer in 3rd channel),normalize (use BIASToTexture)
  // call appropriate this->Map()
  // generate pLUT_ with size of sink.pixelcount with correct interpolationtype
  // fill pLut from newsink
  BIASASSERT(sink.GetPixelCount()!=0);

  incomplete_ = false;

  ROI_ = *src.GetROI();

  //USED FOR NN MAPPING
  BWM_LUT_Entry_Nearest*   current_NN = NULL;

  //USED FOR BILINEAR MAPPING
  BWM_LUT_Entry_Bilinear*  current_Bi = NULL;

  //USED FOR TRILINEAR MAPPING
  BWM_LUT_Entry_Trilinear* current_Tri = NULL;

  if (pLUT_Nea_ != NULL) {
    delete[] pLUT_Nea_;
    pLUT_Nea_ = NULL;
  }
  if (pLUT_Bil_ != NULL) {
    delete[] pLUT_Bil_;
    pLUT_Bil_ = NULL;
  }
  if (pLUT_Tri_ != NULL) {
    delete[] pLUT_Tri_;
    pLUT_Tri_ = NULL;
  }

  switch (method)
    {
    case MapNearestNeighbor:
      pLUT_Nea_ = new BWM_LUT_Entry_Nearest[sink.GetPixelCount()];
      current_NN = &pLUT_Nea_[0];
      pLUT_Nea_size = sink.GetPixelCount();
      break;
    case MapBilinear:
      pLUT_Bil_ = new BWM_LUT_Entry_Bilinear[sink.GetPixelCount()];
      current_Bi = &pLUT_Bil_[0];
      pLUT_Bil_size = sink.GetPixelCount();
      break;
    case MapBicubic:
      break;
    case MapTrilinear:
      // check if pyramid is still large enough, mapping might have changed,
      // e.g. other homography on same image! P.S. Check if pyramid exists too
      //pyramid_.SetPyramidSize(1);
      pyramidSize_ = 1;
      if (pyramid_.IsEmpty()) {
        BuildPyramid_(src);
      }
      if (autoPyramidSize_) UpdatePyramidSize(*sink.GetROI(),
                                              src.GetWidth(), src.GetHeight());
      pLUT_Tri_ = new BWM_LUT_Entry_Trilinear[sink.GetPixelCount()];
      current_Tri = &pLUT_Tri_[0];
      pLUT_Tri_size = sink.GetPixelCount();
      break;
    case MapAnisotropic:
      //pLUT_Ani_ = new BWM_LUT_Entry_Anisotropic[sink.GetPixelCount()];
      break;
    case InterpolationMethod_UNKNOWN:
      break;
    }
  HomgPoint2D sink_p, src_p;

  for(unsigned int j = 0; j < sink.GetHeight(); j++){
    for(unsigned int i = 0; i < sink.GetWidth(); i++){
      sink_p[0] = i + offsetX_;
      sink_p[1] = j + offsetY_;
      sink_p[2] = 1;

      //if pixel can be unmapped and its unmapped coordinates are within the
      //source image, then:

      //WARNING: it checks if the pixels coordinates are within
      //         [0,width-1]x[0,height-1] range, which can be a bit confusing
      //         while the unmapped coordinates, like (-0.3; -0.4) theoretical-
      //         -ly can be (and are) used by Nearest Neighbour Mapping, for
      //         example (and it WON't be used by lookup-Nearest Neighbour
      //         mapping in our case, so the images made by Direct-Mapping and
      //         lookuptable-mapping can be slightly different)
      int off=0;
      if(method==MapBilinear)off=1;//avoid reading outside the src image during filtering
                                   //TODO check if this is required for other filter methods too!
      if(GetSourceCoordinates(sink_p, src_p) != -1 &&
         src_p[0]+off < src.GetWidth() &&
         src_p[1]+off < src.GetHeight() &&
         src_p[0] >= 0 &&
         src_p[1] >= 0){
        double xx = src_p[0];
        double yy = src_p[1];

        double xx_dx = xx - floor(xx);
        double yy_dy = yy - floor(yy);

        double local = 0;
        int    local_low = 0;
        double local_dif = 0;

        double offset = pyramid_.GetPositionOffset();

        int src_channels = src.GetChannelCount();
        int src_width    = src.GetWidth();

        int pyr_channels_low, pyr_channels_high;
        int pyr_wid_low, pyr_wid_high;

        switch (method)
          {
          case MapNearestNeighbor:
            //mapped-check
            current_NN->mapped = true;
            current_NN->index =
              src_channels*
              ((int)rint(yy)*src_width+(int)rint(xx));
            current_NN++;
            break;
          case MapBilinear:
            //mapped-check
            current_Bi->mapped = true;

            //pixel-indexes
            current_Bi->index_0_0 =
              src_channels*
              ((int)floor(yy)*src_width+(int)floor(xx));            
            current_Bi->index_0_1 =
              src_channels*
              ((int)floor(yy+1)*src_width+(int)floor(xx)); 
            current_Bi->index_1_0 =
              src_channels*
              ((int)floor(yy)*src_width+(int)floor(xx+1));
            current_Bi->index_1_1 =
              src_channels*
              ((int)floor(yy+1)*src_width+(int)floor(xx+1));

            BIASASSERT(current_Bi->index_0_0<=(int)(src.GetPixelCount()*src_channels));
            BIASASSERT(current_Bi->index_0_1<=(int)(src.GetPixelCount()*src_channels));
            BIASASSERT(current_Bi->index_1_0<=(int)(src.GetPixelCount()*src_channels));
            BIASASSERT(current_Bi->index_1_1<=(int)(src.GetPixelCount()*src_channels));

            //pixel-weights
            current_Bi->XY_weight =    xx_dx *   yy_dy;
            current_Bi->Xy_weight =    xx_dx *(1-yy_dy);
            current_Bi->xY_weight = (1-xx_dx)*   yy_dy;
            current_Bi->xy_weight = (1-xx_dx)*(1-yy_dy);

            current_Bi++;
            break;
          case MapBicubic:
            break;
          case MapTrilinear:
            // WARNING: Direct mapping without LUT creates additional pyramid levels
            // whenever needed to avoid aliasing. LUT mapping clips these extra values
            // resulting in small differences between images mapped with these methods.

            //mapped-check
            current_Tri->mapped = true;

            //compute layer for trilinear filtering, 0 means best resolution
            ComputeLocalPyramidLayer_(sink_p, local);
            if (local > (int)pyramid_.size()-1)
              local = (int)pyramid_.size()-1;
            if (local < 0)
              local = 0;

            local_low = (int)floor(local);
            local_dif = local - (double)local_low;

            pyr_channels_low  = pyramid_[local_low]->GetChannelCount();
            pyr_wid_low       = pyramid_[local_low]->GetWidth();

            // check if local_low+1 is a level of the pyramid else clip
            if (local_low + 1 < (int)pyramid_.size())
            {
              pyr_channels_high = pyramid_[local_low+1]->GetChannelCount();
              pyr_wid_high      = pyramid_[local_low+1]->GetWidth();
              current_Tri->pyr_level_high = local_low+1;
            }
            else
            {
              pyr_channels_high = pyramid_[local_low]->GetChannelCount();
              pyr_wid_high      = pyramid_[local_low]->GetWidth();
              current_Tri->pyr_level_high = local_low;
            }

            //pyramid-level information
            current_Tri->pyr_level_low  = local_low;

            current_Tri->pyr_diff       =     local_dif;
            current_Tri->pyr_diff_inv   = 1.0-local_dif;

            //STAGE 0

            //xx and yy for the lower pyramid-level
            for(int loc = 0; loc < local_low; loc++){
              xx = (xx - offset)*0.5;
              yy = (yy - offset)*0.5;
            }

            // check indixes
            double test;
            if (GetImageValue_(xx, yy, *(pyramid_[local_low]), 0, test))
            {
              current_Tri->mapped = false;
              current_Tri++;
              break;
            }

            //Indexes for the lower level

            current_Tri->stage_0.index_0_0 =
              pyr_channels_low*
              ((int)floor(yy)*pyr_wid_low+(int)floor(xx));
            current_Tri->stage_0.index_0_1 =
              pyr_channels_low*
              ((int)floor(yy+1)*pyr_wid_low+(int)floor(xx));
            current_Tri->stage_0.index_1_0 =
              pyr_channels_low*
              ((int)floor(yy)*pyr_wid_low+(int)floor(xx+1));
            current_Tri->stage_0.index_1_1 =
              pyr_channels_low*
              ((int)floor(yy+1)*pyr_wid_low+(int)floor(xx+1));

            xx_dx = xx - floor(xx);
            yy_dy = yy - floor(yy);

            //Weights for the lower level
            current_Tri->stage_0.XY_weight =    xx_dx *   yy_dy;
            current_Tri->stage_0.Xy_weight =    xx_dx *(1-yy_dy);
            current_Tri->stage_0.xY_weight = (1-xx_dx)*   yy_dy;
            current_Tri->stage_0.xy_weight = (1-xx_dx)*(1-yy_dy);

            //STAGE 1
            if (local != (int)pyramid_.size()-1){
              xx = (xx - offset)*0.5;
              yy = (yy - offset)*0.5;
            }

            // check indixes
            if (current_Tri->pyr_level_high!=local_low) {
              if (GetImageValue_(xx, yy, *(pyramid_[current_Tri->pyr_level_high]), 0, test) != 0) {
                current_Tri->mapped = false;
                current_Tri++;
                break;
              }
            }
            //Indexes for the higher level
            current_Tri->stage_1.index_0_0 =
              pyr_channels_high*
              ((int)floor(yy)*pyr_wid_high+(int)floor(xx));
            current_Tri->stage_1.index_0_1 =
              pyr_channels_high*
              ((int)floor(yy+1)*pyr_wid_high+(int)floor(xx));
            current_Tri->stage_1.index_1_0 =
              pyr_channels_high*
              ((int)floor(yy)*pyr_wid_high+(int)floor(xx+1));
            current_Tri->stage_1.index_1_1 =
              pyr_channels_high*
              ((int)floor(yy+1)*pyr_wid_high+(int)floor(xx+1));

            xx_dx = xx - floor(xx);
            yy_dy = yy - floor(yy);

            //Weights for the higher level
            current_Tri->stage_1.XY_weight =    xx_dx *   yy_dy;
            current_Tri->stage_1.Xy_weight =    xx_dx *(1-yy_dy);
            current_Tri->stage_1.xY_weight = (1-xx_dx)*   yy_dy;
            current_Tri->stage_1.xy_weight = (1-xx_dx)*(1-yy_dy);

            current_Tri++;
            break;
          case MapAnisotropic:
            break;
          case InterpolationMethod_UNKNOWN:
            break;
          }
      } else {
        incomplete_ = true;
        switch (method)
          {
          case MapNearestNeighbor:
            //mapped-check
            current_NN->mapped = false;
            current_NN++;
            break;
          case MapBilinear:
            //mapped-check
            current_Bi->mapped = false;
            current_Bi++;
            break;
          case MapBicubic:
            break;
          case MapTrilinear:
            //mapped-check
            current_Tri->mapped = false;
            current_Tri++;
            break;
          case MapAnisotropic:
            break;
          case InterpolationMethod_UNKNOWN:
            break;
          }
      }
    }
  }
  if (incomplete_) {
    incomplete_ = false;
    return +1;
  }
  return 0;
}

template <class InputStorageType, class OutputStorageType>
int BackwardMapping<InputStorageType, OutputStorageType>::
MapWithLookupTable(const Image<InputStorageType>& src,
                         Image<OutputStorageType>& sink,
                         InterpolationMethod method) {
  //ROI_ = *src.GetROI();
#ifdef BIAS_DEBUG
  if (src.IsPlanar() || sink.IsPlanar()) {
    BIASERR("not implemented for planar images! (not yet ?)");
    BIASABORT;
  }
#endif
  OutputStorageType* Data = sink.GetImageData();

  int wid = sink.GetWidth();
  int hei = sink.GetHeight();
  unsigned int channels = src.GetChannelCount();
  BIASASSERT(channels == sink.GetChannelCount());

  int srcwid = src.GetWidth();
  int srchei = src.GetHeight();
  int biggy2 = srcwid*srchei;
  int biggy = wid*hei;
  int biggy3 = biggy2 * channels;
  //USED FOR NN MAPPING
  BWM_LUT_Entry_Nearest*   current_NN;

  //USED FOR BILINEAR MAPPING
  BWM_LUT_Entry_Bilinear*  current_Bi;

  switch (method)
    {
    case MapNearestNeighbor:
      //if lookups were not computed
      if(pLUT_Nea_ == NULL){
        BIASWARN("LOOKUP-TABLE FOR NN IS NOT READY! generating!!!");
        PrepareLookupTableMapping(src, sink, method);
      }

      //Lookup-Table must be same size as sink
      BIASASSERT(biggy == pLUT_Nea_size);

      current_NN = &pLUT_Nea_[0];

      for(int i = 0; i < biggy; i++){
        if(pLUT_Nea_[i].mapped){
          if(current_NN->index < biggy3) {
        if (channels == 1)
          *Data++ =(OutputStorageType)src.GetImageData()[current_NN->index];
        else {
          for (unsigned int c=0;c<channels;c++) {
        Data[c] =(OutputStorageType)
          src.GetImageData()[(current_NN->index)+c];
          } // for channels
          Data+=channels; 
        } // else if (channels==1)
      } // if (current_NN->index < biggy2)
          else {
            //don't do this if preinitialization of
            //of image background shall be effective
            //(e.g. color unmapped areas in rectified image pair differently)
            //            ClipValue_(0, *Data++);
            Data+=channels;
          }
        }
        else{
            //don't do this if preinitialization of
            //of image background shall be effective
            //(e.g. color unmapped areas in rectified image pair differently)
            //            ClipValue_(0, *Data++);
          Data+=channels;
       }
        current_NN++;
      }
      break;
    case MapBilinear:
      //if lookups were not computed
      if(pLUT_Bil_ == NULL){
        BIASWARN("LOOKUP-TABLE FOR BILINEAR IS NOT READY! generating!!!");
        PrepareLookupTableMapping(src, sink, method);
      }

      //Lookup-Table must be same size as sink
      BIASASSERT(biggy == pLUT_Bil_size);

      current_Bi = &pLUT_Bil_[0];

      for(int i = 0; i < biggy; i++){  
        for(unsigned int c=0; c<channels; c++) {
          if(current_Bi->mapped){    
            ClipValue_(
              (current_Bi->xy_weight*
              src.GetImageData()[current_Bi->index_0_0+c] +
              current_Bi->xY_weight*
              src.GetImageData()[current_Bi->index_0_1+c] +
              current_Bi->Xy_weight*
              src.GetImageData()[current_Bi->index_1_0+c] +
              current_Bi->XY_weight*
               src.GetImageData()[current_Bi->index_1_1+c]),               
              *Data++ ) ;
          }
          else{
            //don't do this if preinitialization of
            //of image background shall be effective
            //(e.g. color unmapped areas in rectified image pair differently)
            //        ClipValue_(0, *Data++);
            Data++;
          }
        }//end of channels
        current_Bi++;
      }
      break;
    case MapBicubic:
      BIASERR("LOOKUP-TABLES FOR BICUBIC ARE NOT IMPLEMENTED!");
      BIASABORT;

      //Lookup-Table must be same size as sink
      //BIASASSERT(biggy == sizeof(pLUT_Bic_));

      break;
    case MapTrilinear: {
      //if lookups were not computed
      if(pLUT_Tri_ == NULL){
        BIASWARN("LOOKUP-TABLE FOR TRILINEAR IS NOT READY! generating!!!");
        PrepareLookupTableMapping(src, sink, method);
      }

      //Lookup-Table must be same size as sink
      BIASASSERT(biggy == pLUT_Tri_size);

      //build pyramid for mapping
      BuildPyramid_(src);
      if (autoPyramidSize_)
        UpdatePyramidSize(*sink.GetROI(), src.GetWidth(), src.GetHeight());

//       cout << sizeof(BWM_LUT_Entry_Bilinear) << endl;
//       cout << sizeof(BWM_LUT_Entry_Trilinear) << endl;

      BWM_LUT_Entry_Trilinear* current_Tri = &pLUT_Tri_[0];
      for(int i = 0; i < biggy; i++){
        for(unsigned int c=0; c<channels; c++) {
          if(current_Tri->mapped) {
            // average bilinear values from two closest pyramid levels
            const BWM_LUT_Entry_Bilinear& pstruct_0 = current_Tri->stage_0;
            const BWM_LUT_Entry_Bilinear& pstruct_1 = current_Tri->stage_1;

            const InputStorageType* pPyrLow  =
              pyramid_[current_Tri->pyr_level_low]->GetImageData();
            const InputStorageType* pPyrHigh =
              pyramid_[current_Tri->pyr_level_high]->GetImageData();
            // shouldnt we use clipvalue here to avoid range errors ?
                        ClipValue_(
              (current_Tri->pyr_diff_inv *
              (pstruct_0.xy_weight * pPyrLow[pstruct_0.index_0_0+c] +
              pstruct_0.xY_weight * pPyrLow[pstruct_0.index_0_1+c] +
              pstruct_0.Xy_weight * pPyrLow[pstruct_0.index_1_0+c] +
              pstruct_0.XY_weight * pPyrLow[pstruct_0.index_1_1+c]) +
              current_Tri->pyr_diff *
              (pstruct_1.xy_weight * pPyrHigh[pstruct_1.index_0_0+c] +
              pstruct_1.xY_weight * pPyrHigh[pstruct_1.index_0_1+c] +
              pstruct_1.Xy_weight * pPyrHigh[pstruct_1.index_1_0+c] +
               pstruct_1.XY_weight * pPyrHigh[pstruct_1.index_1_1+c])),
              *Data++);

          } else {
            // write zero, we might as well only increase pointer here
            // and leave the init up to the caller
            //*Data++ = OutputStorageType(0);
            Data++;
          }
    }
        current_Tri++;
      }
      break;
    }
    case MapAnisotropic:
      BIASERR("LOOKUP-TABLES FOR ANISOTROPIC ARE NOT IMPLEMENTED!");
      BIASABORT;

      //Lookup-Table must be same size as sink
      //BIASASSERT(biggy == sizeof(pLUT_Ani_));

      break;
    default:
      BIASERR("FALSE INTERPOLATION METHOD!");
      return -1;
    }
  return 0;
  //  BIASASSERT(sink.samesize(pLUT_));
  // switch interpolationmethod_
  // cast plut to pluntinterpolationmethod
  // run over all pixels in sink,
}






template <class InputStorageType, class OutputStorageType>
void BackwardMapping<InputStorageType, OutputStorageType>::
GenerateTestImage(Image<InputStorageType>& im, bool highFrequencyCross,
                  InputStorageType dark,
                  InputStorageType bright,
                  const Matrix3x3<double>& Hinv) {
  InputStorageType dbright = bright - dark;
  if (im.IsEmpty()) {
    im.Init(640, 480);
  }
  BIASASSERT(im.GetWidth()>10);
  BIASASSERT(im.GetHeight()>10);
  BIASASSERT(im.GetChannelCount()==1);

  HomgPoint2D point(0.0, 0.0, 1.0);
  Random R;

#define OFFSET 0.5
  for(unsigned int y=0; y<im.GetHeight(); y++) {
    for(unsigned int x=0; x<im.GetWidth(); x++) {
      double angle = 0, angleaverage = 0;
      // avoid aliasing in original image by presmoothing:
      for (unsigned int offset = 0; offset<5; offset++) {
        point[1] = double(y)-(im.GetHeight()/2.0);
        point[0] = double(x)-(im.GetWidth()/2.0);
        switch (offset) {
        case 0: point[0] -= OFFSET; point[1] -= OFFSET; break;
        case 1: point[0] += OFFSET; point[1] -= OFFSET; break;
        case 2: point[0] += OFFSET; point[1] += OFFSET; break;
        case 3: point[0] -= OFFSET; point[1] += OFFSET; break;
        }
        point = Hinv * point;
        point.Homogenize();
        angle = atan2(point[1]-OFFSET, point[0]-OFFSET);
        if (angle<0.0) angle += M_PI*2.0;
        //angle = int(rint(angle / (M_PI*2.0) * 32.0))%2;
        angle = (angle / (M_PI*2.0) * 32.0);
        angleaverage += fabs((angle / 2.0) - rint((angle / 2.0)));

      }
      angleaverage *= 0.2;

      im.GetImageDataArray()[y][x] = dark +
        (InputStorageType)(angleaverage * 2.0*double(dbright));
    }
  }
  if (highFrequencyCross) {
    int crosspos = 2; // 2 mean half image size, 3 means third
    int offset = 12; // dont go closer than this to center
    for(unsigned int y=0; y<im.GetHeight()/crosspos-offset; y++) {
      im.GetImageDataArray()[y][im.GetWidth()/crosspos] =
        (InputStorageType)(dark+bright*((y+1)%2));
      im.GetImageDataArray()[y][im.GetWidth()/crosspos+1] =
        (InputStorageType)(dark+dbright*(y%2));
    }
    for(unsigned int x=0; x<im.GetWidth()/crosspos-offset; x++) {
      im.GetImageDataArray()[im.GetHeight()/crosspos-2][x] =
        im.GetImageDataArray()[im.GetHeight()/crosspos-1][x] =
        im.GetImageDataArray()[im.GetHeight()/crosspos][x] =
        im.GetImageDataArray()[im.GetHeight()/crosspos+1][x] =
        (InputStorageType)(dark+dbright*((x/4)%2));
    }
    for(unsigned int y=im.GetHeight()/crosspos+offset; y<im.GetHeight()-3;
        y++)     {
      im.GetImageDataArray()[y][im.GetWidth()/crosspos-1] =
        im.GetImageDataArray()[y][im.GetWidth()/crosspos-2] =
        im.GetImageDataArray()[y][im.GetWidth()/crosspos] =
        im.GetImageDataArray()[y][im.GetWidth()/crosspos+1] =
        im.GetImageDataArray()[y][im.GetWidth()/crosspos+2] =
        (InputStorageType)(dark+dbright*((y/5)%2));
    }
    for(unsigned int x=im.GetWidth()/crosspos+offset; x<im.GetWidth()-1; x++) {
      im.GetImageDataArray()[im.GetHeight()/crosspos-1][x] =
        im.GetImageDataArray()[im.GetHeight()/crosspos][x] =
        im.GetImageDataArray()[im.GetHeight()/crosspos+1][x] =
        (InputStorageType)(dark+dbright*((x/3)%2));
    }

  }
}

// ######################### template instantiation  ########################


template class BackwardMapping<unsigned char, unsigned char>;
template class BackwardMapping<float, float>;
template class BackwardMapping<unsigned char, float>;
template class BackwardMapping<float, unsigned char>;

#if defined(BUILD_IMAGE_CHAR)
template class BackwardMapping<unsigned char, char>;
template class BackwardMapping<char, char>;
#endif

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

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

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

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

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

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

#if defined(BUILD_IMAGE_DOUBLE)
template class BackwardMapping<double,double>;
#endif

} // namespace BIAS