biasshowepiwx.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
*/

#ifdef WIN32
// disable warning noise
#  pragma warning(disable: 4251) // STL MAP DLL
#  pragma warning(disable: 4005) // VS8 wx
#  pragma warning(disable: 4996) // wx localtime

#  include <conio.h> // _getch
#endif //WIN32

#include <Gui/ConsoleRedirectIO.hh>
#include <Gui/ImageCanvas.hh>
#include <Gui/PlainImageCanvas.hh>
#include <Gui/StringConv.hh>
#include <Base/Common/W32Compat.hh>
#include <Base/Image/ImageIO.hh>
#include <Base/ImageUtils/ImageDraw.hh>
#include <Base/Image/ImageConvert.hh>
#include <Base/ImageUtils/Bresenham.hh>
#include <Image/Camera.hh>
#include <Geometry/PMatrix.hh>
#include <Geometry/FMatrix.hh>
#include <Geometry/EpipolarLine.hh>
#include <Geometry/ProjectionParametersSpherical.hh>

#include <wx/wx.h>
#include <wx/splitter.h>

#include <iostream>

using namespace std;
using namespace BIAS;


/**
   @file
   @ingroup g_tools
   @brief show epipolar lines for image pair (WX Version), see biasshowepiwx.cpp
   @author koeser
*/

/*
===============================================================================
Constants
===============================================================================
*/

const int HALF_CUTOUT_SIZE = 30;
const int CUTOUT_SIZE =  2 * HALF_CUTOUT_SIZE;
const int SCALE = 4;

/*
===============================================================================
ComputeCutoutImage
===============================================================================
*/

void ComputeCutoutImage(Image<unsigned char>& src, Image<unsigned char>& dst,
                        int srcX, int srcY, int srcWidth, int srcHeight, int scale)
{
  int dstWidth  = srcWidth  * scale;
  int dstHeight = srcHeight * scale;

  if ((int)dst.GetWidth() != dstWidth || (int)dst.GetHeight() != dstHeight)
    dst.Release();

  if (dst.IsEmpty())
    dst.Init(dstWidth, dstHeight, 3);

  unsigned char** srcIda = src.GetImageDataArray();
  unsigned char** dstIda = dst.GetImageDataArray();

  for (int x = srcX; x < srcX + srcWidth; x++)
  {
    for (int y = srcY; y < srcY + srcHeight; y++)
    {
      unsigned char pixel[3];
      if (x < 0 || y < 0 || x >= (int)src.GetWidth() || y >= (int)src.GetHeight())
      {
        pixel[0] = 0;
        pixel[1] = 0;
        pixel[2] = 0;
      }
      else
      {
        if (src.GetColorModel() == ImageBase::CM_Grey)
        {
          pixel[0] = srcIda[y][x];
          pixel[1] = srcIda[y][x];
          pixel[2] = srcIda[y][x];
        }
        else
        {
          // assuming CM_RGB
          pixel[0] = srcIda[y][3*x+0];
          pixel[1] = srcIda[y][3*x+1];
          pixel[2] = srcIda[y][3*x+2];
        }
      }

      int startX = (x - srcX) * scale;
      int startY = (y - srcY) * scale;

      for (int i = startX; i < startX + scale; i++)
      {
        for (int j = startY; j < startY + scale; j++)
        {
          if (i >= 0 && j >= 0 && i < (int)dst.GetWidth() && j < (int)dst.GetHeight())
          {
            dstIda[j][3*i+0] = pixel[0];
            dstIda[j][3*i+1] = pixel[1];
            dstIda[j][3*i+2] = pixel[2];
          }
        }
      }
    }
  }
}

/** \cond HIDDEN_SYMBOLS
 * \class CutoutFrame - Frame showing a cutout of a larger image
*/
class CutoutFrame : public wxFrame
{
public:
  CutoutFrame(wxWindow* parent, const wxString& title);

  void ShowImage(Image<unsigned char>& image, const string& name);
  void SetTransform(int offsetX, int offsetY, int scale);

private:
  int               _OffsetX, _OffsetY, _Scale;
  PlainImageCanvas* _ImageCanvas;

  void OnMotion(wxMouseEvent& event);
  void OnLeaveWindow(wxMouseEvent& event);
  void OnClose(wxCloseEvent& event);
  void OnKeyDown(wxKeyEvent& keyEvent);
  void OnCharEvt(wxKeyEvent& keyEvent);


  DECLARE_EVENT_TABLE()
};
/** \endcond */

BEGIN_EVENT_TABLE(CutoutFrame, wxFrame)
EVT_MOTION(CutoutFrame::OnMotion)
EVT_LEAVE_WINDOW(CutoutFrame::OnLeaveWindow)
EVT_CLOSE(CutoutFrame::OnClose)
EVT_KEY_DOWN(CutoutFrame::OnKeyDown)
EVT_CHAR(CutoutFrame::OnCharEvt)
END_EVENT_TABLE()

void CutoutFrame::OnKeyDown(wxKeyEvent& event)
{
  if (GetParent() != NULL)
  {
    GetParent()->GetEventHandler()->ProcessEvent(event);
  }
}


void CutoutFrame::OnCharEvt(wxKeyEvent& event)
{
  if (GetParent() != NULL)
  {
    GetParent()->GetEventHandler()->ProcessEvent(event);
  }
}


CutoutFrame::CutoutFrame(wxWindow* parent, const wxString& title)
  : wxFrame(parent, -1, title,
            wxDefaultPosition, wxDefaultSize,
            wxDEFAULT_FRAME_STYLE | wxFRAME_FLOAT_ON_PARENT)
{
  _OffsetX = 0;
  _OffsetY = 0;
  _Scale   = 1;

  _ImageCanvas = new PlainImageCanvas(this);

  CreateStatusBar(1);

  wxBoxSizer* sizer = new wxBoxSizer(wxHORIZONTAL);
  sizer->Add(_ImageCanvas);
  SetSizer(sizer);
}

void CutoutFrame::ShowImage(Image<unsigned char>& image, const string& name)
{
  _ImageCanvas->Show(image, name);
  _ImageCanvas->SetScrollbars(0, 0, 0, 0);
  GetSizer()->SetItemMinSize(_ImageCanvas, image.GetWidth(), image.GetHeight());
  Fit();
}

void CutoutFrame::SetTransform(int offsetX, int offsetY, int scale)
{
  _OffsetX = offsetX;
  _OffsetY = offsetY;
  _Scale   = scale;
}

void CutoutFrame::OnMotion(wxMouseEvent& event)
{
  wxClientDC dc(this);
  PrepareDC(dc);
  wxPoint pos = event.GetPosition();
  int x = dc.DeviceToLogicalX(pos.x) / _Scale + _OffsetX;
  int y = dc.DeviceToLogicalY(pos.y) / _Scale + _OffsetY;
  stringstream text;
  text << "(x, y) = (" << x << ", " << y << ")";
  GetStatusBar()->SetStatusText(AsciiToWx(text.str()));
}

void CutoutFrame::OnLeaveWindow(wxMouseEvent&)
{
  GetStatusBar()->SetStatusText(wxT(""));
}

void CutoutFrame::OnClose(wxCloseEvent& event)
{
  if (event.CanVeto())
  {
    Hide();
    event.Veto();
  }
  else
  {
    Destroy();
  }
}

/*
===============================================================================
MainFrame
===============================================================================
*/

 enum{
   ID_THICKER_LINE= wxID_HIGHEST+100,
   ID_THINNER_LINE,
   ID_SAVE_IMAGES
 };
 /** \cond HIDDEN_SYMBOLS */
class MainFrame : public wxFrame
{
public:
  MainFrame();

  void SetImages(Image<unsigned char>& image1, const Projection& proj1,
                 const string& name1,
                 Image<unsigned char>& image2, const Projection& proj2,
                 const string& name2, bool homographymode = false);
  double GenerateVirtualPerspectiveCamera(ProjectionParametersBase *paramsSpherical1,
                                 ProjectionParametersPerspective *fakeParamsPersp1,
                                 const double& maximagewidth);
private:
  /// if in homography only mode:
  bool onlyHomography_;
  /// this cotains the last clicked point in image1 backprojected onto
  /// plane at infinity
  HomgPoint3D pointatinfinity_;
  PlainImageCanvas*     _ImageCanvas1;
  PlainImageCanvas*     _ImageCanvas2;
  Image<unsigned char>  _Image1, _Image2;
  Image<unsigned char>  PaintedImage1_, PaintedImage2_;
  Projection            _Projection1, _Projection2;
  std::string           _ImageName1, _ImageName2;
  CutoutFrame*          _CutoutFrame1;
  CutoutFrame*          _CutoutFrame2;
  wxPoint               _SelectedPoint1; // selected point in left image
  wxPoint               _SelectedPoint2; // selected point in right image

  void OnLeftDown(wxMouseEvent& event);
  void OnKeyDown(wxKeyEvent& keyEvent);
  void OnCharEvt(wxKeyEvent& keyEvent);
  void ThickerLine(wxCommandEvent& event);
  void ThinnerLine(wxCommandEvent& event);
  void SaveImages(wxCommandEvent& event);

  wxMenuBar* menuBar_;
  wxMenu* file_;
//   wxMenuItem* save_;

  wxMenu* drawing_;
//   wxMenuItem* thicker_;
//   wxMenuItem* thinner_;


  int epiLineWidth_;
  DECLARE_EVENT_TABLE()
};
/** \endcond */
BEGIN_EVENT_TABLE(MainFrame, wxFrame)
EVT_LEFT_DOWN(MainFrame::OnLeftDown)
EVT_MENU(ID_THICKER_LINE, MainFrame::ThickerLine)
EVT_MENU(ID_THINNER_LINE, MainFrame::ThinnerLine)
EVT_MENU(ID_SAVE_IMAGES, MainFrame::SaveImages)
END_EVENT_TABLE()

void MainFrame::SaveImages(wxCommandEvent& event){
  ImageIO::Save("paintedImage1.mip",PaintedImage1_);
  ImageIO::Save("paintedImage2.mip",PaintedImage2_);
}

void MainFrame::ThickerLine(wxCommandEvent& event) 
{ 
  epiLineWidth_ ++; 
  cout<<" "<<epiLineWidth_<<endl;
    int x1, y1, x2, y2;
    _ImageCanvas1->GetViewStart(&x1, &y1);
    _ImageCanvas1->GetViewStart(&x2, &y2);
    _ImageCanvas1->Show(_Image1, _ImageName1);
    _ImageCanvas2->Show(_Image2, _ImageName2);
    _ImageCanvas1->Scroll(x1, y1);
    _ImageCanvas1->Scroll(x2, y2);

}

void MainFrame::ThinnerLine(wxCommandEvent& event) 
{ 
  epiLineWidth_ --; 
  if(epiLineWidth_<1) epiLineWidth_=1;
  cout<<" "<<epiLineWidth_<<endl;
     int x1, y1, x2, y2;
    _ImageCanvas1->GetViewStart(&x1, &y1);
    _ImageCanvas1->GetViewStart(&x2, &y2);
    _ImageCanvas1->Show(_Image1, _ImageName1);
    _ImageCanvas2->Show(_Image2, _ImageName2);
    _ImageCanvas1->Scroll(x1, y1);
    _ImageCanvas1->Scroll(x2, y2);

}




MainFrame::MainFrame()
  : wxFrame(NULL, -1, wxT("BiasShowEpiWX")),
    pointatinfinity_ (0.0, 0.0, 0.0, 0.0), epiLineWidth_(1)
{
  _CutoutFrame1 = new CutoutFrame(this, wxT("Cutout 1"));
  _CutoutFrame2 = new CutoutFrame(this, wxT("Cutout 2"));
  onlyHomography_ = false;
  CreateStatusBar();

  _ImageCanvas1 = new PlainImageCanvas(this, GetStatusBar(), 0);
  _ImageCanvas2 = new PlainImageCanvas(this, GetStatusBar(), 0);

   menuBar_ = new wxMenuBar;
   file_ = new wxMenu;
   file_->Append(ID_SAVE_IMAGES, wxT("&Save"));

   drawing_ = new wxMenu;
   drawing_->Append(ID_THICKER_LINE, wxT("&Increase Line Width"));
   drawing_->Append(ID_THINNER_LINE, wxT("&Decrease Line Width"));


   menuBar_->Append(file_, wxT("&File"));
   menuBar_->Append(drawing_, wxT("&Display"));
   SetMenuBar(menuBar_);


//   wxAcceleratorEntry entries[2];
//   entries[0].Set(wxACCEL_NORMAL, WXK_SPACE, ID_FAKE_UP);
//   //entries[1].Set(wxACCEL_NORMAL, WXK_DOWN, ID_FAKE_DOWN);
//   wxAcceleratorTable accel(2, entries);
//   SetAcceleratorTable(accel);
}



void DrawEpipole(Image<unsigned char>& image1,
                 const Projection& proj1, const Projection& proj2) {
  HomgPoint2D orig_epipole;
  const unsigned char Value[3]={255,0,0};
  bool epipoleIsInImage =
    proj1.DoesPointProjectIntoImage(HomgPoint3D(proj2.GetC()),
                                           orig_epipole);

  if (epipoleIsInImage) {
    ImageDraw<unsigned char>::CircleCenter(image1,
                                           int(rint(orig_epipole[0])),
                                           int(rint(orig_epipole[1])),
                                           5, Value);
    ImageDraw<unsigned char>::CircleCenter(image1,
                                           int(rint(orig_epipole[0])),
                                           int(rint(orig_epipole[1])),
                                           6, Value);
  }

}

void ShowSphericalParameters(Image<unsigned char>& image1,
                             const Projection& proj1) {
  const ProjectionParametersSpherical *pPPS =
    dynamic_cast<const ProjectionParametersSpherical* >(proj1.GetParameters());

  if (pPPS!=NULL) {
    const unsigned char Blue[3]={0,0,255};
    const unsigned char Green[3]={0,255,0};
    int radius = (int(rint(pPPS->GetRadius())));
    HomgPoint2D pp(0,0,1);
    pPPS->GetPrincipal(pp[0], pp[1]);
    ImageDraw<unsigned char>::CircleCenter(image1, int(rint(pp[0])),
                                           int(rint(pp[1])),
                                           int(rint(radius)),
                                           Green);
    ImageDraw<unsigned char>::CircleCenter(image1, int(rint(pp[0])),
                                           int(rint(pp[1])),
                                           int(rint(radius+1)),
                                           Green);
    for (unsigned int i = 0; i<4; i++) {
      Vector3<double> ray(tan(30.0*double(i)*M_PI / 180.0),0,1);
      HomgPoint2D point = pPPS->ProjectLocal(ray);
      if (point.NormL2()<1e-10) break;
      ImageDraw<unsigned char>::CircleCenter(image1, int(rint(pp[0])),
                                             int(rint(pp[1])),
                                             int(rint((point-pp).NormL2())),
                                             Blue);
    }
  }
}

void MainFrame::SetImages(Image<unsigned char>& image1,
                          const Projection& proj1, const string& name1,
                          Image<unsigned char>& image2,
                          const Projection& proj2, const string& name2,
                          bool homographymode)
{
  PaintedImage1_ = _Image1 = image1;
  PaintedImage2_ = _Image2 = image2;
  _Projection1 = proj1;
  _Projection2 = proj2;

  ProjectionParametersPerspective *ppp1 =
    dynamic_cast<ProjectionParametersPerspective *>(_Projection1.GetParameters());
  ProjectionParametersPerspective *ppp2 =
    dynamic_cast<ProjectionParametersPerspective *>(_Projection2.GetParameters());
  if (ppp1!=NULL && ppp2!=NULL) {
    PMatrix P1 = ppp1->GetP();
    PMatrix P2 = ppp2->GetP();
    FMatrix F(P1, P2);
    cout<<"F is "<<F<<endl;
  }







  _ImageName1 = name1;
  _ImageName2 = name2;

  if (((_Projection1.GetC()-_Projection2.GetC()).NormL2()<1e-8)
      || homographymode) {
    BIASWARN(endl<<endl<<endl<<
             "Homography mode or very small baseline, using only H_infinity!"
             <<endl<<endl<<endl);
    onlyHomography_ = true;
  }

  DrawEpipole(_Image1, _Projection1, _Projection2);
  DrawEpipole(_Image2, _Projection2, _Projection1);

  ShowSphericalParameters(_Image1, _Projection1);
  ShowSphericalParameters(_Image2, _Projection2);

  _ImageCanvas1->Show(_Image1, _ImageName1);
  _ImageCanvas2->Show(_Image2, _ImageName2);
  _ImageCanvas1->SetClearDrawingContext(false);
  _ImageCanvas2->SetClearDrawingContext(false);

  wxBoxSizer* sizer = new wxBoxSizer(wxHORIZONTAL);
  sizer->Add(_ImageCanvas1, 1, wxEXPAND);
  sizer->AddSpacer(3);
  sizer->Add(_ImageCanvas2, 1, wxEXPAND);

  sizer->SetItemMinSize(_ImageCanvas1, 100, 100);
  sizer->SetItemMinSize(_ImageCanvas2, 100, 100);

  SetSizer(sizer);
  sizer->SetSizeHints(this);

  SetSize(900, 600); // default main frame size
}

double MainFrame::
GenerateVirtualPerspectiveCamera(ProjectionParametersBase
                                 *paramsSpherical1,
                                 ProjectionParametersPerspective
                                 *fakeParamsPersp1,
                                 const double& maximagewidth) {

  double imageSize = maximagewidth + 1.0;
  // start near 180 degree fov
  double maxangle = M_PI/2.0-0.01;
  KMatrix K(MatrixZero);
  // decrease fov angle until perspective image is smaller than maximagewidth
  while (imageSize>maximagewidth) {
    K = paramsSpherical1->GetFakeKMatrix(imageSize, 1, maxangle);
    maxangle -= 0.01;
    BIASASSERT(maxangle>0.0);
  }
  fakeParamsPersp1->SetC(paramsSpherical1->GetC());
  fakeParamsPersp1->SetQ(paramsSpherical1->GetQ());
  fakeParamsPersp1->SetK(K);
  fakeParamsPersp1->SetImageSize((int)imageSize, (int)imageSize);

  return imageSize;
}

void MainFrame::OnLeftDown(wxMouseEvent& event)
{
  // epipolar lines are nicest handled in perspective images:
  ProjectionParametersPerspective* paramsPersp1 =
    dynamic_cast<ProjectionParametersPerspective*>
    (_Projection1.GetParameters());
  int image1Width  = _Image1.GetWidth();
  int image1Height = _Image1.GetHeight();

  // check whether we must emulate perspective projection:
  ProjectionParametersPerspective fakeParamsPersp1;
  if (paramsPersp1==NULL) {
    // make sure the emulated image coordinates dont go above 100000 otherwise
    // we might get memory problems
    image1Width = image1Height = (int)rint(
                GenerateVirtualPerspectiveCamera(_Projection1.GetParameters(),
                                                 &fakeParamsPersp1,
                                                 image1Width*100.0));
    paramsPersp1 = &fakeParamsPersp1;
    cout<<"Simulating epipolar geometry with perspective image of size "
        <<image1Width<<endl;
  }

  // do the same for the other image
  ProjectionParametersPerspective* paramsPersp2 =
    dynamic_cast<ProjectionParametersPerspective*>(_Projection2.GetParameters());

  int image2Width = _Image2.GetWidth();
  int image2Height = _Image2.GetHeight();

  // check whether we must emulate perspective projection for cam2:
  ProjectionParametersPerspective fakeParamsPersp2;
  if (paramsPersp2==NULL) {

    image2Width = image2Height =
      (int)rint(GenerateVirtualPerspectiveCamera(_Projection2.GetParameters(),
                                                 &fakeParamsPersp2,
                                                 image2Width * 100.0));
    paramsPersp2 = &fakeParamsPersp2;
  }

  PMatrix P1 = paramsPersp1->GetP();
  PMatrix P2 = paramsPersp2->GetP();
  if (event.GetEventObject() == _ImageCanvas1)  {
    wxClientDC dc(this);
    _ImageCanvas1->PrepareDC(dc);
    //wxPoint pos = event.GetPosition();
    _SelectedPoint1.x = dc.DeviceToLogicalX(event.GetX());
    _SelectedPoint1.y = dc.DeviceToLogicalY(event.GetY());


    HomgPoint2D point(_SelectedPoint1.x, _SelectedPoint1.y, 1.0);

    // undistort the clicked point
    if (paramsPersp1== NULL)  {
      Vector3<double> ray, pos;
      _Projection1.GetParameters()->UnProjectLocal(point, pos, ray);
      if (ray.NormL2()<1e-12) {
        BIASERR("clicked on point out of valid area !");
      } else {
        point = paramsPersp1->ProjectLocal(ray);
      }
    }
    else // no fake KMatrix, just perspective distortion
      {
        paramsPersp1->Undistort(point);
      }
    point.Homogenize();

    // images including selected point and epipolar lines
    PaintedImage1_ = _Image1;
    PaintedImage2_ = _Image2;

    // --- LEFT IMAGE ---


    // draw the selected point in PaintedImage1
    ImageDraw<unsigned char>::RectangleCenter(PaintedImage1_,
                                              _SelectedPoint1.x,
                                              _SelectedPoint1.y,
                                              CUTOUT_SIZE);

    Vector3<double> p, dir;
    _Projection1.UnProjectToRay(HomgPoint2D(_SelectedPoint1.x, _SelectedPoint1.y, 1.0), p, dir);
    pointatinfinity_ =
      HomgPoint3D(dir);
    pointatinfinity_[3] = 0.0;
    cout<<"Point at infinity is "<<pointatinfinity_<<endl;
    cout<<"line width "<<epiLineWidth_<<endl;
    HomgPoint2D pointimage2(0,0,0);
    if (_Projection2.DoesPointProjectIntoImage(pointatinfinity_,
                                               pointimage2)) {
      pointimage2.Homogenize();
      // draw the transferred point in PaintedImage2
      ImageDraw<unsigned char>::RectangleCenter(PaintedImage2_,
                                                int(rint(pointimage2[0])),
                                                int(rint(pointimage2[1])),
                                                5);
    } else {
      BIASERR("Point at infinity does not project into image 2:"<<pointimage2);
    }
    if (!onlyHomography_){
      // compute epipoles in big image
      FMatrix F12(P1, P2);
      HomgPoint2D epipole1, epipole2;
      F12.GetEpipolesHomogenized(epipole1, epipole2);

      // compute the projection of the epipolar plane for that point
      EpipolarLine eline1 = epipole1.CrossProduct(point);
      if (eline1.NormL2()>1e-10) {

        HomgPoint2D clickedpoint(_SelectedPoint1.x, _SelectedPoint1.y);
        Vector3<double> pointOnRay, dirRay;
        _Projection1.UnProjectToRay(clickedpoint, pointOnRay, dirRay);

        EpipolarLine::ProjectEpipolarPlane(_Projection1,
                                           PaintedImage1_,
                                           _Projection2.GetC(),
                                           dirRay,
                                           (unsigned char*)NULL,
                                           epiLineWidth_);
        _Projection1.UnProjectToRay(clickedpoint, pointOnRay, dirRay);
        EpipolarLine::ProjectEpipolarPlane(_Projection2,
                                           PaintedImage2_,
                                           _Projection1.GetC(),
                                           dirRay,
                                           (unsigned char*)NULL,
                                           epiLineWidth_);


      }  else {
        BIASERR("You clicked ("<<_SelectedPoint1.x<<";"<<_SelectedPoint1.y
                <<" which is undefined or pretty close to the epipole "
                <<" , in idealized perspective image "<<point<<" and epipole="
                <<epipole1
                <<" resulting in eline "<<eline1);
      }
    }
    // show images (and restore scrollbar positions after that)
    int x1, y1, x2, y2;
    _ImageCanvas1->GetViewStart(&x1, &y1);
    _ImageCanvas1->GetViewStart(&x2, &y2);
    _ImageCanvas1->Show(PaintedImage1_, _ImageName1);
    _ImageCanvas2->Show(PaintedImage2_, _ImageName2);
    _ImageCanvas1->Scroll(x1, y1);
    _ImageCanvas1->Scroll(x2, y2);

    // --- CUTOUT FOR LEFT IMAGE ---
    Image<unsigned char> CutoutImage;
    ComputeCutoutImage(_Image1, CutoutImage,
                       _SelectedPoint1.x - HALF_CUTOUT_SIZE,
                       _SelectedPoint1.y - HALF_CUTOUT_SIZE,
                       CUTOUT_SIZE, CUTOUT_SIZE, SCALE);
    ImageDraw<unsigned char>::
      RectangleCenter(CutoutImage,
               2 + (_SelectedPoint1.x - (_SelectedPoint1.x - HALF_CUTOUT_SIZE)) * SCALE,
               2 + (_SelectedPoint1.y - (_SelectedPoint1.y - HALF_CUTOUT_SIZE)) * SCALE, SCALE);
    ImageDraw<unsigned char>::RectangleCenter(CutoutImage,
               2 + (_SelectedPoint1.x - (_SelectedPoint1.x - HALF_CUTOUT_SIZE)) * SCALE,
               2 + (_SelectedPoint1.y - (_SelectedPoint1.y - HALF_CUTOUT_SIZE)) * SCALE, SCALE * 4);

    _CutoutFrame1->ShowImage(CutoutImage, _ImageName1);
    _CutoutFrame1->SetTransform(_SelectedPoint1.x - HALF_CUTOUT_SIZE,
                                _SelectedPoint1.y - HALF_CUTOUT_SIZE, SCALE);
    _CutoutFrame1->Show(TRUE);
    _CutoutFrame1->Fit();


  }
  // store selected point when clicked in image 2
  if (event.GetEventObject() == _ImageCanvas2)
    {
      wxClientDC dc(this);
      _ImageCanvas2->PrepareDC(dc);
      //wxPoint pos = event.GetPosition();
      _SelectedPoint2.x = dc.DeviceToLogicalX(event.GetX());
      _SelectedPoint2.y = dc.DeviceToLogicalY(event.GetY());
    }

  // update cutout of second image on click or when epipolar line was changed
  if ((event.GetEventObject() == _ImageCanvas1 && _CutoutFrame2->IsShown()) ||
      event.GetEventObject() == _ImageCanvas2)
    {
      Image<unsigned char> CutoutImage;
      ComputeCutoutImage(_Image2, CutoutImage,
                         _SelectedPoint2.x - HALF_CUTOUT_SIZE,
                         _SelectedPoint2.y - HALF_CUTOUT_SIZE,
                         CUTOUT_SIZE, CUTOUT_SIZE, SCALE);

      HomgPoint2D point(_SelectedPoint1.x, _SelectedPoint1.y, 1.0);

      HomgPoint2D epipole1, epipole2;
      Matrix3x3<double> CutoutTransform, CutoutTransformInv;
      CutoutTransform.SetIdentity();
      CutoutTransform[0][0] = CutoutTransform[1][1] = 1.0 / (double)SCALE;
      CutoutTransform[0][2] = (double)(_SelectedPoint2.x - HALF_CUTOUT_SIZE);
      CutoutTransform[1][2] = (double)(_SelectedPoint2.y - HALF_CUTOUT_SIZE);

      HomgPoint2D endofepipolarline(0,0,0);
      if (_Projection2.DoesPointProjectIntoImage(pointatinfinity_,
                                                 endofepipolarline)) {

      #ifdef BIAS_DEBUG
        BIASASSERT(CutoutTransform.GetInverse(CutoutTransformInv)==0);
      #else
        CutoutTransform.GetInverse(CutoutTransformInv);
      #endif
        endofepipolarline = CutoutTransformInv * endofepipolarline;
        endofepipolarline.Homogenize();
        int windowsize = 2;
        if (CutoutImage.IsPositionInImage(int(rint(endofepipolarline[0]))
                                          - windowsize,
                                          int(rint(endofepipolarline[1]))
                                          - windowsize) &&
            CutoutImage.IsPositionInImage(int(rint(endofepipolarline[0]))
                                          + windowsize,
                                          int(rint(endofepipolarline[1]))
                                          + windowsize)) {

          // draw the transferred point in PaintedImage2
          ImageDraw<unsigned char>::
            RectangleCenter(CutoutImage,
                            int(rint(endofepipolarline[0])),
                            int(rint(endofepipolarline[1])),
                            windowsize);
        }
      } else {
        BIASERR("Point at infinity does not project into image 2:"
                <<endofepipolarline);
      }



      if (!onlyHomography_) {
        FMatrix FCutout(P1, P2);
        FCutout = CutoutTransform.Transpose() * FCutout;
        FCutout.GetEpipolesHomogenized(epipole1, epipole2);
        HomgPoint2D::EPosition epos2 =
          epipole2.DeterminePosition(CUTOUT_SIZE*SCALE, CUTOUT_SIZE*SCALE);

        // draw epipolar line
        EpipolarLine eline;
        eline = FCutout * point;
        eline.Recalc(CUTOUT_SIZE*SCALE, CUTOUT_SIZE*SCALE, epipole2, epos2);


        ProjectionParametersBase *pB = _Projection2.GetParameters()->Clone();
        double cx,cy;
        Vector3<double> pointOnRay, dirRay;
        pB->GetPrincipal(cx,cy);
        cx -= _SelectedPoint2.x - HALF_CUTOUT_SIZE;
        cy -= _SelectedPoint2.y - HALF_CUTOUT_SIZE;
        pB->SetPrincipal(cx,cy);
        pB->SetImageSize(CUTOUT_SIZE,CUTOUT_SIZE);
        pB->Rescale(1.0/double(SCALE));
        Projection CutoutProjection(*pB);
        HomgPoint2D clickedpoint(_SelectedPoint1.x, _SelectedPoint1.y);
        _Projection1.UnProjectToRay(clickedpoint, pointOnRay, dirRay);
        EpipolarLine::ProjectEpipolarPlane(CutoutProjection,
                                           CutoutImage,
                                           _Projection1.GetC(),
                                           dirRay,
                                           (unsigned char*)NULL, 
                                           epiLineWidth_
                                           ); 
        delete pB;




      }

      _CutoutFrame2->ShowImage(CutoutImage, _ImageName2);
      _CutoutFrame2->SetTransform(_SelectedPoint2.x - HALF_CUTOUT_SIZE,
                                  _SelectedPoint2.y - HALF_CUTOUT_SIZE, SCALE);
      _CutoutFrame2->Show(TRUE);
      _CutoutFrame2->Fit();
    }
}


void Usage() {
  cout<<endl<<"biasshowepiwx is part of the BIAS library."<<endl;
  cout<<"See www.mip.informatik.uni-kiel.de for new versions."<<endl<<endl;
  cout<<"Usage: "<<endl;
  cout<<"         biasshowepiwx img1 img2 [projection1 projection2]"<<endl;
  cout<<endl<<endl;
  cout<<"This tool loads two images with known camera parameters (intrinsics"
      <<" and extrinsics) and shows them in a window."
      <<" Click a point in the left image and you will see the ray in space"
      <<" projected into the other image (epipolar line). If you have a"
      <<" fisheye image or an image with radial distortion, the line will"
      <<" be curved. The camera parameters can be in BIAS::Camera format"
      <<" within the image meta data (pgm, ppm, mip) or can be given as "
      <<" separate files on the command line. File formats supported are "
      <<" those of BIAS::Projection::Load(), namely"<<endl
      <<"  - projection file in xml written by BIAS::Projection::XMLWrite"
      <<endl
      <<"  - camera parameter file compatible with BIAS::CameraParam class"
      <<endl
      <<"  - .cam file of BBC's FreeD tracking system"<<endl
      <<"  - projection matrix in .mat file compatible with BIAS::PMatrix"
      <<endl<<endl;
  cout<<"The small square shows the vanishing point of the epipolar line "
      <<" where the feature would lie with pure camera rotation, while the "
      <<"circle indicates the other end of the epipolar line (the epipole),"
      <<" which is the projection of the other camera center into the image."
      <<" Both ends and even the line as a whole may or may not be in the"
      <<" image, depending on your camera parameters."<<endl
      <<" Left clicking opens an enlarged window."<<endl;

}



/*
===============================================================================
IntroductionFrame
===============================================================================
*/
enum {
  ID_SELECT_LEFT_IMAGE = wxID_HIGHEST+188,
  ID_SELECT_RIGHT_IMAGE,
  ID_SELECT_LEFT_PROJ,
  ID_SELECT_RIGHT_PROJ,
  ID_TEXT_LEFT_IMAGE,
  ID_TEXT_RIGHT_IMAGE,
  ID_TEXT_LEFT_PROJ,
  ID_TEXT_RIGHT_PROJ,
  ID_BUTTON_OK,
  ID_BUTTON_CANCEL,
  ID_STATIC
};
/** \cond HIDDEN_SYMBOLS */
class IntroductionFrame : public wxDialog
{
public:
  IntroductionFrame(wxFrame* parent,
                    const wxString& title,
                    const wxPoint& pos,
                    const wxSize& size);

  void OnButtonLeftImage( wxCommandEvent& event );
  void OnButtonRightImage( wxCommandEvent& event );
  void OnButtonLeftProj( wxCommandEvent& event );
  void OnButtonRightProj( wxCommandEvent& event );
  void OnButtonOK( wxCommandEvent& event );
  void OnButtonCancel( wxCommandEvent& event );
  void SetStrings(wxString* image1,wxString* image2,
    wxString* proj1,wxString* proj2){
    image1_ = image1;
    image2_ = image2;
    proj1_ = proj1;
    proj2_ = proj2;
  };

protected:
  wxTextCtrl *itemTextCtrlLeftP;
  wxTextCtrl *itemTextCtrlRightP;
  wxTextCtrl *itemTextCtrlRight;
  wxTextCtrl *itemTextCtrlLeft;
  wxString *image1_,*image2_,*proj1_,*proj2_;

  DECLARE_EVENT_TABLE()
  DECLARE_CLASS( IntroductionFrame )
};
/** \endcond */

IMPLEMENT_CLASS( IntroductionFrame, wxFrame )

BEGIN_EVENT_TABLE(IntroductionFrame, wxDialog)
  EVT_BUTTON( ID_SELECT_LEFT_IMAGE, IntroductionFrame::OnButtonLeftImage )
  EVT_BUTTON( ID_SELECT_RIGHT_IMAGE,IntroductionFrame::OnButtonRightImage )
  EVT_BUTTON( ID_SELECT_LEFT_PROJ,  IntroductionFrame::OnButtonLeftProj )
  EVT_BUTTON( ID_SELECT_RIGHT_PROJ, IntroductionFrame::OnButtonRightProj )
  EVT_BUTTON( ID_BUTTON_OK,         IntroductionFrame::OnButtonOK)
  EVT_BUTTON( ID_BUTTON_CANCEL,     IntroductionFrame::OnButtonCancel)
END_EVENT_TABLE()

IntroductionFrame::
IntroductionFrame(wxFrame* parent,const wxString& title,const wxPoint& pos,const wxSize& size):
wxDialog(parent, -1, title, pos, size, wxSTAY_ON_TOP|wxCAPTION|wxRESIZE_BORDER)
{
 // the sizer
  wxFlexGridSizer* itemFlexGridSizer = new wxFlexGridSizer(3, 0, 0);
  this->SetSizer(itemFlexGridSizer);
  wxSize textSize(370,30);
  ///////////////////////////////////////////////////////////////////////
  wxStaticText* itemStaticText0 =
    new wxStaticText( this, ID_STATIC, wxT("First Image"),
                      wxDefaultPosition, wxDefaultSize, 0 );
  wxButton* itemButtonLeft=new wxButton( this, ID_SELECT_LEFT_IMAGE,
                                  wxT("Select"), wxDefaultPosition,
                                  wxDefaultSize, 0 );
  itemTextCtrlLeft = new wxTextCtrl( this, ID_TEXT_LEFT_IMAGE,
                                         wxT(""), wxDefaultPosition,
                                         textSize, 0 );
  itemFlexGridSizer->Add(itemStaticText0, 0,
                         wxALIGN_LEFT|wxALIGN_CENTER_VERTICAL|wxALL,
                         5);
  itemFlexGridSizer->Add(itemTextCtrlLeft, 0,
                         wxALIGN_LEFT|wxALIGN_CENTER_VERTICAL|wxALL,
                         5);
  itemFlexGridSizer->Add(itemButtonLeft, 0,
            wxALIGN_CENTER_HORIZONTAL|wxALIGN_CENTER_VERTICAL|wxALL, 5);
  ///////////////////////////////////////////////////////////////////////
  wxStaticText* itemStaticText01 =
    new wxStaticText( this, ID_STATIC, wxT("Projection"),
                      wxDefaultPosition, wxDefaultSize, 0 );
  wxButton* itemButtonLeftP=new wxButton( this, ID_SELECT_LEFT_PROJ,
                                  wxT("Select"), wxDefaultPosition,
                                  wxDefaultSize, 0 );
  itemTextCtrlLeftP = new wxTextCtrl( this, ID_TEXT_LEFT_PROJ,
                                         wxT(""), wxDefaultPosition,
                                         textSize, 0 );
  itemFlexGridSizer->Add(itemStaticText01, 0,
                         wxALIGN_LEFT|wxALIGN_CENTER_VERTICAL|wxALL,
                         5);
  itemFlexGridSizer->Add(itemTextCtrlLeftP, 0,
                         wxALIGN_LEFT|wxALIGN_CENTER_VERTICAL|wxALL,
                         5);
  itemFlexGridSizer->Add(itemButtonLeftP, 0,
            wxALIGN_CENTER_HORIZONTAL|wxALIGN_CENTER_VERTICAL|wxALL, 5);

  ///////////////////////////////////////////////////////////////////////
  wxStaticText* itemStaticText1 =
    new wxStaticText( this, ID_STATIC, wxT("Second Image"),
                      wxDefaultPosition, wxDefaultSize, 0 );
  wxButton* itemButtonRight=new wxButton( this, ID_SELECT_RIGHT_IMAGE,
                                  wxT("Select"), wxDefaultPosition,
                                  wxDefaultSize, 0 );
  itemTextCtrlRight= new wxTextCtrl( this, ID_TEXT_RIGHT_IMAGE,
                                         wxT(""), wxDefaultPosition,
                                         textSize, 0 );
  itemFlexGridSizer->Add(itemStaticText1, 0,
                         wxALIGN_LEFT|wxALIGN_CENTER_VERTICAL|wxALL,
                         5);
  itemFlexGridSizer->Add(itemTextCtrlRight, 0,
                         wxALIGN_LEFT|wxALIGN_CENTER_VERTICAL|wxALL,
                         5);
  itemFlexGridSizer->Add(itemButtonRight, 0,
            wxALIGN_CENTER_HORIZONTAL|wxALIGN_CENTER_VERTICAL|wxALL, 5);
  ///////////////////////////////////////////////////////////////////////
  wxStaticText* itemStaticText02 =
    new wxStaticText( this, ID_STATIC, wxT("Projection"),
                      wxDefaultPosition, wxDefaultSize, 0 );
  wxButton* itemButtonRightP=new wxButton( this, ID_SELECT_RIGHT_PROJ,
                                  wxT("Select"), wxDefaultPosition,
                                  wxDefaultSize, 0 );
  itemTextCtrlRightP = new wxTextCtrl( this, ID_TEXT_RIGHT_PROJ,
                                         wxT(""), wxDefaultPosition,
                                         textSize, 0 );
  itemFlexGridSizer->Add(itemStaticText02, 0,
                         wxALIGN_LEFT|wxALIGN_CENTER_VERTICAL|wxALL,
                         5);
  itemFlexGridSizer->Add(itemTextCtrlRightP, 0,
                         wxALIGN_LEFT|wxALIGN_CENTER_VERTICAL|wxALL,
                         5);
  itemFlexGridSizer->Add(itemButtonRightP, 0,
            wxALIGN_CENTER_HORIZONTAL|wxALIGN_CENTER_VERTICAL|wxALL, 5);
  ////////////////////////////////////////////////////////////////////////
  wxButton* itemButtonOK = new wxButton( this, ID_BUTTON_OK, _("OK"),
                       wxDefaultPosition, wxDefaultSize, 0 );
  wxButton* itemButtonCancel = new wxButton( this, ID_BUTTON_CANCEL, _("Cancel"),
                       wxDefaultPosition, wxDefaultSize, 0 );

  itemFlexGridSizer->Add(itemButtonOK, 0,
                         wxALIGN_LEFT|wxALIGN_CENTER_VERTICAL|wxALL,
                         5);
  itemFlexGridSizer->Add(itemButtonCancel, 0,
                         wxALIGN_LEFT|wxALIGN_CENTER_VERTICAL|wxALL,
                         5);

  ////////////////////////////////////////////////////////////////////
  Centre();
}

void IntroductionFrame::
OnButtonLeftImage( wxCommandEvent& event ){
  wxString fileName;
  wxFileDialog dialog(this,wxT("Select Image"),wxT(""),wxT(""),
                      wxT("Images (*.mip;*.pgm;*.ppm;*.jpg;*.png)|*.mip;*.pgm;*.ppm;*.jpg;*.png|All files (*.*)|*.*"),
                      wxFD_OPEN);
  int ret = dialog.ShowModal();
  if(ret == wxID_OK){
    fileName = dialog.GetPath();
    itemTextCtrlLeft->SetValue(fileName);
  }
  else
    return;
 }
void IntroductionFrame::
OnButtonRightImage( wxCommandEvent& event ){
  wxString fileName;
  wxFileDialog dialog(this,wxT("Select Image"),wxT(""),wxT(""),
                      wxT("Images (*.mip;*.pgm;*.ppm;*.jpg;*.png)|*.mip;*.pgm;*.ppm;*.jpg;*.png|All files (*.*)|*.*"),
                      wxFD_OPEN);
  int ret = dialog.ShowModal();
  if(ret == wxID_OK){
    fileName = dialog.GetPath();
    itemTextCtrlRight->SetValue(fileName);
  }
  else
    return;
}
void IntroductionFrame::
OnButtonLeftProj( wxCommandEvent& event ){
  wxString fileName;
  wxFileDialog dialog(this,wxT("Select Projection"),wxT(""),wxT(""),
                      wxT("Projections (*.proj;*.xml)|*.proj;*.xml|All files (*.*)|*.*"),
                      wxFD_OPEN);
  int ret = dialog.ShowModal();
  if(ret == wxID_OK){
    fileName = dialog.GetPath();
    itemTextCtrlLeftP->SetValue(fileName);
  }
  else
    return;
}
void IntroductionFrame::
OnButtonRightProj( wxCommandEvent& event ){
  wxString fileName;
  wxFileDialog dialog(this,wxT("Select Projection"),wxT(""),wxT(""),
                      wxT("Projections (*.proj;*.xml)|*.proj;*.xml|All files (*.*)|*.*"),
                      wxFD_OPEN);
  int ret = dialog.ShowModal();
  if(ret == wxID_OK){
    fileName = dialog.GetPath();
    itemTextCtrlRightP->SetValue(fileName);
  }
  else
    return;
}
void IntroductionFrame::
OnButtonOK( wxCommandEvent& event ){
  *image1_= itemTextCtrlLeft->GetValue();
  *image2_= itemTextCtrlRight->GetValue();
  *proj1_= itemTextCtrlLeftP->GetValue();
  *proj2_ = itemTextCtrlRightP->GetValue();

  this->Close();
}
void IntroductionFrame::
OnButtonCancel( wxCommandEvent& event ){
   exit(0);
}

/**
 * \cond HIDDEN_SYMBOLS
 * \class BiasShowEpiApp */
class BiasShowEpiApp : public wxApp
{
public:
  virtual bool OnInit();
};
/** \endcond */

bool BiasShowEpiApp::OnInit()
{
  // utility: create a Redirection console for cout, cerr etc.
  ConsoleRedirectIO();

  // this parameter should be parsed from the command line
  bool homographymode = false;
  MainFrame* frame = new MainFrame();

  wxString image1,image2, proj1, proj2;
  if (argc == 2 && (strcmp(WxToAscii(argv[1]),"-h")==0 ||
      strcmp(WxToAscii(argv[1]),"--help")==0 ||
      strcmp(WxToAscii(argv[1]),"/?")==0)) {
    Usage();
#ifdef WIN32 // avoid closing the window
    cout<<"press any key"; _getch();
#endif
    return false;
  }

  string fileName[2];
  string projName[2];
  // parse command line
  if (argc < 3)
  {
    Usage();
    IntroductionFrame* introframe;
    introframe = new IntroductionFrame(NULL,wxT("Select Images and Projections"),wxDefaultPosition,wxSize(600,250));
    introframe->SetStrings(&image1,&image2,&proj1, &proj2);
    introframe->ShowModal();
    fileName[0] = WxToAscii(image1);
    fileName[1] = WxToAscii(image2);
    projName[0] = WxToAscii(proj1);
    projName[1] = WxToAscii(proj2);
  }
  else{
    fileName[0] = WxToAscii(argv[1]);
    fileName[1] = WxToAscii(argv[2]);
     // use loaded projections, if given
    if (argc >= 5)
    {
      for (int i = 0; i < 2; i++)
      {
        projName[i] = WxToAscii(argv[i+3]);
      }
    }
  }
  // load images
  Camera<unsigned char> image[2];
  for (int i = 0; i < 2; i++)
  {
    ImageBase mytmpimg;
    if (ImageIO::Load(fileName[i],  mytmpimg) != 0)
    {
      cout << "Failed to load image: " << fileName[i] << endl;
      #ifdef WIN32 // avoid closing the window
      cout<<"press any key"; _getch();
      #endif
      return false;
    }
    cout << "loaded image, storage type is " << mytmpimg.GetStorageType() << endl;
    if(mytmpimg.GetStorageType()==ImageBase::ST_float){
      Image<float> mytmpFL(mytmpimg);
      if(mytmpFL.GetMaxPixelValue() <= 1.0){
        mytmpFL.ScaleShift(255, 1.0);
      }
      ImageConvert::ConvertST(mytmpFL, image[i], ImageBase::ST_unsignedchar);
    } else {
      ImageConvert::ConvertST(mytmpimg, image[i], ImageBase::ST_unsignedchar);
    }
    image[i].SetMetaData(*mytmpimg.GetMetaData());
    image[i].ParseMetaData();
    if(projName[i] != ""){
      Projection projection;
      if (projection.Load(projName[i]) != 0)
      {
        cout << "Failed to load projection:" << projName[i] << endl;
        #ifdef WIN32 // avoid closing the window
        cout<<"press any key"; _getch();
        #endif
        return false;
      }
      image[i].SetProj(projection);
    }
  }

  // take projections from images or create them from P matrices
  Projection proj[2];
  for (int i = 0; i < 2; i++)
  {
    if (image[i].IsProjValid())
    {
      proj[i] = image[i].GetProj();
    }
    else if (image[i].IsPValid())
    {
      ProjectionParametersPerspective params;
      params.SetP(image[i].GetP());
      params.SetImageSize(image[i].GetWidth(), image[i].GetHeight());
      proj[i] = Projection(params);
    }
    else
    {
      cout << "No projection and no P matrix for image: " << fileName[i] << endl;
      #ifdef WIN32 // avoid closing the window
      cout<<"press any key"; _getch();
      #endif
      return false;
    }


    ProjectionParametersPerspective* pPPP =
      dynamic_cast<ProjectionParametersPerspective* >(proj[i].GetParameters());
    if (pPPP!=NULL) {
      // have perspective cam, avoid ugly distortion double overflow!
      pPPP->UpdateMinZLocal();
    }
  }

  // convert images to RGB
  Image<unsigned char> rgbImage[2];
  if (ImageConvert::ToRGB(image[0], rgbImage[0]) != 0 ||
      ImageConvert::ToRGB(image[1], rgbImage[1]) != 0)
  {
    cout << "Failed to convert images to RGB" << endl;
    #ifdef WIN32 // avoid closing the window
    cout<<"press any key"; _getch();
    #endif
    return false;
  }


  // show frame
  frame->Show(TRUE);
  frame->SetImages(rgbImage[0], proj[0], fileName[0],
                   rgbImage[1], proj[1], fileName[1], homographymode);
  SetTopWindow(frame);

  return TRUE;
}

IMPLEMENT_APP(BiasShowEpiApp)